EgammaMCTools/src/ElectronMCTruthFinder.cc

0001 #include "RecoEgamma/EgammaMCTools/interface/ElectronMCTruthFinder.h"
0002
0003 #include <algorithm>
0004
0005 ElectronMCTruthFinder::ElectronMCTruthFinder() {}
0006
0007 ElectronMCTruthFinder::~ElectronMCTruthFinder() {}
0008
0009 std::vector<ElectronMCTruth> ElectronMCTruthFinder::find(const std::vector<SimTrack>& theSimTracks,
0010                                                          const std::vector<SimVertex>& theSimVertices) {
0011   //std::cout << "  ElectronMCTruthFinder::find " << std::endl;
0012
0013   std::vector<ElectronMCTruth> result;
0014
0015   // Local variables
0016   //const int SINGLE=1;
0017   //const int DOUBLE=2;
0018   //const int PYTHIA=3;
0019   //const int ELECTRON_FLAV=1;
0020   //const int PIZERO_FLAV=2;
0021   //const int PHOTON_FLAV=3;
0022
0023   //int ievtype=0;
0024   //int ievflav=0;
0025
0026   std::vector<SimTrack> electronTracks;
0027   SimVertex primVtx;
0028
0029   fill(theSimTracks, theSimVertices);
0030
0031   int iPV = -1;
0032   //int partType1=0;
0033   //int partType2=0;
0034   std::vector<SimTrack>::const_iterator iFirstSimTk = theSimTracks.begin();
0035   if (!(*iFirstSimTk).noVertex()) {
0036     iPV = (*iFirstSimTk).vertIndex();
0037
0038     int vtxId = (*iFirstSimTk).vertIndex();
0039     primVtx = theSimVertices[vtxId];
0040
0041     //partType1 = (*iFirstSimTk).type();
0042
0043     //std::cout <<  " Primary vertex id " << iPV << " first track type " << (*iFirstSimTk).type() << std::endl;
0044   } else {
0045     //std::cout << " First track has no vertex " << std::endl;
0046   }
0047
0048   // CLHEP::HepLorentzVector primVtxPos= primVtx.position();
0049   math::XYZTLorentzVectorD primVtxPos(
0050       primVtx.position().x(), primVtx.position().y(), primVtx.position().z(), primVtx.position().e());
0051
0052   // Look at a second track
0053   iFirstSimTk++;
0054   //if ( iFirstSimTk!=  theSimTracks.end() ) {
0055   //
0056   //  if (  (*iFirstSimTk).vertIndex() == iPV) {
0057   //    partType2 = (*iFirstSimTk).type();
0058   //
0059   //    //std::cout <<  " second track type " << (*iFirstSimTk).type() << " vertex " <<  (*iFirstSimTk).vertIndex() << std::endl;
0060   //
0061   //  } else {
0062   //    //std::cout << " Only one kine track from Primary Vertex " << std::endl;
0063   //  }
0064   //}
0065
0066   //std::cout << " Loop over all particles " << std::endl;
0067
0068   for (std::vector<SimTrack>::const_iterator iSimTk = theSimTracks.begin(); iSimTk != theSimTracks.end(); ++iSimTk) {
0069     if ((*iSimTk).noVertex())
0070       continue;
0071
0072     //int vertexId = (*iSimTk).vertIndex();
0073     //SimVertex vertex = theSimVertices[vertexId];
0074
0075     //std::cout << " Particle type " <<  (*iSimTk).type() << " Sim Track ID " << (*iSimTk).trackId() << " momentum " << (*iSimTk).momentum() <<  " vertex position " << vertex.position() << " vertex ID " << vertexId  << std::endl;
0076     if ((*iSimTk).vertIndex() == iPV) {
0077       if (std::abs((*iSimTk).type()) == 11) {
0078         //std::cout << " Found a primary electron with ID  " << (*iSimTk).trackId() << " momentum " << (*iSimTk).momentum() <<  std::endl;
0079         electronTracks.push_back(*iSimTk);
0080       }
0081     }
0082   }
0083
0084   /// Now store the electron truth
0085   std::vector<CLHEP::Hep3Vector> bremPos;
0086   std::vector<CLHEP::HepLorentzVector> pBrem;
0087   std::vector<float> xBrem;
0088
0089   for (std::vector<SimTrack>::iterator iEleTk = electronTracks.begin(); iEleTk != electronTracks.end(); ++iEleTk) {
0090     //std::cout << " Looping on the primary electron pt  " << std::sqrt((*iEleTk).momentum().perp2()) << " electron track ID " << (*iEleTk).trackId() << std::endl;
0091
0092     SimTrack trLast = (*iEleTk);
0093     unsigned int eleId = (*iEleTk).trackId();
0094     float remainingEnergy = trLast.momentum().e();
0095     //    CLHEP::HepLorentzVector motherMomentum = (*iEleTk).momentum();
0096     //    CLHEP::HepLorentzVector primEleMom = (*iEleTk).momentum();
0097     math::XYZTLorentzVectorD motherMomentum(
0098         (*iEleTk).momentum().x(), (*iEleTk).momentum().y(), (*iEleTk).momentum().z(), (*iEleTk).momentum().e());
0099     math::XYZTLorentzVectorD primEleMom(motherMomentum);
0100     int eleVtxIndex = (*iEleTk).vertIndex();
0101
0102     bremPos.clear();
0103     pBrem.clear();
0104     xBrem.clear();
0105
0106     for (std::vector<SimTrack>::const_iterator iSimTk = theSimTracks.begin(); iSimTk != theSimTracks.end(); ++iSimTk) {
0107       if ((*iSimTk).noVertex())
0108         continue;
0109       if ((*iSimTk).vertIndex() == iPV)
0110         continue;
0111
0112       //std::cout << " (*iEleTk)->trackId() " << (*iEleTk).trackId() << " (*iEleTk)->vertIndex() "<< (*iEleTk).vertIndex()  << " (*iSimTk).vertIndex() "  <<  (*iSimTk).vertIndex() << " (*iSimTk).type() " <<   (*iSimTk).type() << " (*iSimTk).trackId() " << (*iSimTk).trackId() << std::endl;
0113
0114       int vertexId1 = (*iSimTk).vertIndex();
0115       const SimVertex& vertex1 = theSimVertices[vertexId1];
0116       int vertexId2 = trLast.vertIndex();
0117       //SimVertex vertex2 = theSimVertices[vertexId2];
0118
0119       int motherId = -1;
0120
0121       if ((vertexId1 == vertexId2) && ((*iSimTk).type() == (*iEleTk).type()) && trLast.type() == 22) {
0122         //std::cout << " Here a e/gamma brem vertex " << std::endl;
0123
0124         //std::cout << " Secondary from electron:  particle1  type " << (*iSimTk).type() << " trackId " << (*iSimTk).trackId() << " vertex ID " << vertexId1 << " vertex position " << std::sqrt(vertex1.position().perp2()) << " parent index "<< vertex1.parentIndex() << std::endl;
0125
0126         //std::cout << " Secondary from electron:  particle2  type " << trLast.type() << " trackId " <<  trLast.trackId()<< " vertex ID " << vertexId2 << " vertex position " << std::sqrt(vertex2.position().perp2()) << " parent index " << vertex2.parentIndex() << std::endl;
0127
0128         //std::cout << " Electron pt " << std::sqrt((*iSimTk).momentum().perp2()) << " photon pt " <<  std::sqrt(trLast.momentum().perp2()) << "Mother electron pt " <<  sqrt(motherMomentum.perp2()) << std::endl;
0129         //std::cout << " eleId " << eleId << std::endl;
0130         float eLoss = remainingEnergy - ((*iSimTk).momentum() + trLast.momentum()).e();
0131         //std::cout << " eLoss " << eLoss << std::endl;
0132
0133         if (vertex1.parentIndex()) {
0134           unsigned motherGeantId = vertex1.parentIndex();
0135           std::map<unsigned, unsigned>::iterator association = geantToIndex_.find(motherGeantId);
0136           if (association != geantToIndex_.end())
0137             motherId = association->second;
0138
0139           //int motherType = motherId == -1 ? 0 : theSimTracks[motherId].type();
0140           //std::cout << " Parent to this vertex   motherId " << motherId << " mother type " <<  motherType << " Sim track ID " <<  theSimTracks[motherId].trackId() << std::endl;
0141           if (theSimTracks[motherId].trackId() == eleId) {
0142             //std::cout << "  ***** Found the Initial Mother Electron ****   theSimTracks[motherId].trackId() " <<  theSimTracks[motherId].trackId() << " eleId " <<  eleId << std::endl;
0143             eleId = (*iSimTk).trackId();
0144             remainingEnergy = (*iSimTk).momentum().e();
0145             motherMomentum = (*iSimTk).momentum();
0146
0147             pBrem.push_back(CLHEP::HepLorentzVector(
0148                 trLast.momentum().px(), trLast.momentum().py(), trLast.momentum().pz(), trLast.momentum().e()));
0149             bremPos.push_back(CLHEP::HepLorentzVector(
0150                 vertex1.position().x(), vertex1.position().y(), vertex1.position().z(), vertex1.position().t()));
0151             xBrem.push_back(eLoss);
0152           }
0153
0154         } else {
0155           //std::cout << " This vertex has no parent tracks " <<  std::endl;
0156         }
0157       }
0158       trLast = (*iSimTk);
0159
0160     }  // End loop over all SimTracks
0161     //std::cout << " Going to build the ElectronMCTruth: pBrem size " << pBrem.size() << std::endl;
0162     /// here fill the electron
0163     CLHEP::HepLorentzVector tmpEleMom(primEleMom.px(), primEleMom.py(), primEleMom.pz(), primEleMom.e());
0164     CLHEP::HepLorentzVector tmpVtxPos(primVtxPos.x(), primVtxPos.y(), primVtxPos.z(), primVtxPos.t());
0165     result.push_back(ElectronMCTruth(tmpEleMom, eleVtxIndex, bremPos, pBrem, xBrem, tmpVtxPos, (*iEleTk)));
0166
0167   }  // End loop over primary electrons
0168
0169   return result;
0170 }
0171
0172 void ElectronMCTruthFinder::fill(const std::vector<SimTrack>& simTracks, const std::vector<SimVertex>& simVertices) {
0173   //std::cout << "  ElectronMCTruthFinder::fill " << std::endl;
0174
0175   unsigned nVtx = simVertices.size();
0176   unsigned nTks = simTracks.size();
0177
0178   // Empty event, do nothin'
0179   if (nVtx == 0)
0180     return;
0181
0182   // create a map associating geant particle id and position in the
0183   // event SimTrack vector
0184   for (unsigned it = 0; it < nTks; ++it) {
0185     geantToIndex_[simTracks[it].trackId()] = it;
0186     //std::cout << " ElectronMCTruthFinder::fill it " << it << " simTracks[it].trackId() " <<  simTracks[it].trackId() << std::endl;
0187   }
0188 }