EgammaMCTools/src/PizeroMCTruthFinder.cc

0001 #include "RecoEgamma/EgammaMCTools/interface/PizeroMCTruthFinder.h"
0002 #include "RecoEgamma/EgammaMCTools/interface/PhotonMCTruthFinder.h"
0003 #include "RecoEgamma/EgammaMCTools/interface/ElectronMCTruthFinder.h"
0004
0005 #include <algorithm>
0006
0007 PizeroMCTruthFinder::PizeroMCTruthFinder() {
0008   thePhotonMCTruthFinder_ = new PhotonMCTruthFinder();
0009   theElectronMCTruthFinder_ = new ElectronMCTruthFinder();
0010 }
0011
0012 PizeroMCTruthFinder::~PizeroMCTruthFinder() {
0013   delete thePhotonMCTruthFinder_;
0014   delete theElectronMCTruthFinder_;
0015   std::cout << "~PizeroMCTruthFinder" << std::endl;
0016 }
0017
0018 std::vector<PizeroMCTruth> PizeroMCTruthFinder::find(const std::vector<SimTrack>& theSimTracks,
0019                                                      const std::vector<SimVertex>& theSimVertices) {
0020   std::cout << "  PizeroMCTruthFinder::find " << std::endl;
0021
0022   std::vector<PhotonMCTruth> mcPhotons = thePhotonMCTruthFinder_->find(theSimTracks, theSimVertices);
0023
0024   std::vector<PizeroMCTruth> result;
0025   std::vector<PhotonMCTruth> photonsFromPizero;
0026   std::vector<ElectronMCTruth> electronsFromPizero;
0027
0028   // Local variables
0029   //const int SINGLE=1;
0030   //const int DOUBLE=2;
0031   //const int PYTHIA=3;
0032   //const int ELECTRON_FLAV=1;
0033   //const int PIZERO_FLAV=2;
0034   //const int PHOTON_FLAV=3;
0035
0036   //int ievtype=0;
0037   //int ievflav=0;
0038
0039   std::vector<SimTrack> pizeroTracks;
0040   SimVertex primVtx;
0041
0042   fill(theSimTracks, theSimVertices);
0043
0044   int iPV = -1;
0045   //int partType1=0;
0046   //int partType2=0;
0047   std::vector<SimTrack>::const_iterator iFirstSimTk = theSimTracks.begin();
0048   if (!(*iFirstSimTk).noVertex()) {
0049     iPV = (*iFirstSimTk).vertIndex();
0050
0051     int vtxId = (*iFirstSimTk).vertIndex();
0052     primVtx = theSimVertices[vtxId];
0053
0054     //partType1 = (*iFirstSimTk).type();
0055     std::cout << " Primary vertex id " << iPV << " first track type " << (*iFirstSimTk).type() << std::endl;
0056   } else {
0057     std::cout << " First track has no vertex " << std::endl;
0058   }
0059
0060   // CLHEP::HepLorentzVector primVtxPos= primVtx.position();
0061   math::XYZTLorentzVectorD primVtxPos(
0062       primVtx.position().x(), primVtx.position().y(), primVtx.position().z(), primVtx.position().e());
0063
0064   // Look at a second track
0065   iFirstSimTk++;
0066   if (iFirstSimTk != theSimTracks.end()) {
0067     if ((*iFirstSimTk).vertIndex() == iPV) {
0068       //partType2 = (*iFirstSimTk).type();
0069       std::cout << " second track type " << (*iFirstSimTk).type() << " vertex " << (*iFirstSimTk).vertIndex()
0070                 << std::endl;
0071
0072     } else {
0073       std::cout << " Only one kine track from Primary Vertex " << std::endl;
0074     }
0075   }
0076
0077   int npv = 0;
0078
0079   for (std::vector<SimTrack>::const_iterator iSimTk = theSimTracks.begin(); iSimTk != theSimTracks.end(); ++iSimTk) {
0080     if ((*iSimTk).noVertex())
0081       continue;
0082
0083     int vertexId = (*iSimTk).vertIndex();
0084     const SimVertex& vertex = theSimVertices[vertexId];
0085
0086     std::cout << " Particle type " << (*iSimTk).type() << " Sim Track ID " << (*iSimTk).trackId() << " momentum "
0087               << (*iSimTk).momentum() << " vertex position " << vertex.position() << std::endl;
0088
0089     if ((*iSimTk).vertIndex() == iPV) {
0090       npv++;
0091       if (std::abs((*iSimTk).type()) == 111) {
0092         std::cout << " Found a primary pizero with ID  " << (*iSimTk).trackId() << " momentum " << (*iSimTk).momentum()
0093                   << std::endl;
0094
0095         pizeroTracks.push_back(*iSimTk);
0096
0097         // CLHEP::HepLorentzVector momentum = (*iSimTk).momentum();
0098         math::XYZTLorentzVectorD momentum(
0099             (*iSimTk).momentum().x(), (*iSimTk).momentum().y(), (*iSimTk).momentum().z(), (*iSimTk).momentum().e());
0100       }
0101     }
0102   }
0103
0104   std::cout << " There are " << npv << " particles originating in the PV " << std::endl;
0105
0106   //if(npv > 4) {
0107   //  ievtype = PYTHIA;
0108   //} else if(npv == 1) {
0109   //  if( std::abs(partType1) == 11 ) {
0110   //    ievtype = SINGLE; ievflav = ELECTRON_FLAV;
0111   //  } else if(partType1 == 111) {
0112   //    ievtype = SINGLE; ievflav = PIZERO_FLAV;
0113   //  } else if(partType1 == 22) {
0114   //    ievtype = SINGLE; ievflav = PHOTON_FLAV;
0115   //  }
0116   //} else if(npv == 2) {
0117   //  if (  std::abs(partType1) == 11 && std::abs(partType2) == 11 ) {
0118   //    ievtype = DOUBLE; ievflav = ELECTRON_FLAV;
0119   //  } else if(partType1 == 111 && partType2 == 111)   {
0120   //    ievtype = DOUBLE; ievflav = PIZERO_FLAV;
0121   //  } else if(partType1 == 22 && partType2 == 22)   {
0122   //    ievtype = DOUBLE; ievflav = PHOTON_FLAV;
0123   //  }
0124   //}
0125
0126   for (std::vector<SimTrack>::iterator iPizTk = pizeroTracks.begin(); iPizTk != pizeroTracks.end(); ++iPizTk) {
0127     std::cout << " Looping on the primary pizero pt  " << sqrt((*iPizTk).momentum().perp2()) << " pizero track ID "
0128               << (*iPizTk).trackId() << std::endl;
0129
0130     photonsFromPizero.clear();
0131     std::cout << " mcPhotons.size " << mcPhotons.size() << std::endl;
0132     for (std::vector<PhotonMCTruth>::iterator iPho = mcPhotons.begin(); iPho != mcPhotons.end(); ++iPho) {
0133       int phoVtxIndex = (*iPho).vertexInd();
0134       const SimVertex& phoVtx = theSimVertices[phoVtxIndex];
0135       unsigned int phoParentInd = phoVtx.parentIndex();
0136       std::cout << " photon parent vertex index " << phoParentInd << std::endl;
0137
0138       if (phoParentInd == (*iPizTk).trackId()) {
0139         std::cout << "Matched Photon ID " << (*iPho).trackId() << "  vtx " << phoParentInd << " with pizero "
0140                   << (*iPizTk).trackId() << std::endl;
0141         photonsFromPizero.push_back(*iPho);
0142       }
0143     }
0144     std::cout << " Photon matching the pizero vertex " << photonsFromPizero.size() << std::endl;
0145
0146     // build pizero MC thruth
0147     CLHEP::HepLorentzVector tmpMom(
0148         (*iPizTk).momentum().px(), (*iPizTk).momentum().py(), (*iPizTk).momentum().pz(), (*iPizTk).momentum().e());
0149     CLHEP::HepLorentzVector tmpPos(
0150         primVtx.position().x(), primVtx.position().y(), primVtx.position().z(), primVtx.position().t());
0151     result.push_back(PizeroMCTruth(tmpMom, photonsFromPizero, tmpPos));
0152
0153   }  // end loop over primary pizeros
0154
0155   std::cout << " Pizero size " << result.size() << std::endl;
0156
0157   return result;
0158 }
0159
0160 void PizeroMCTruthFinder::fill(const std::vector<SimTrack>& simTracks, const std::vector<SimVertex>& simVertices) {
0161   unsigned nVtx = simVertices.size();
0162   unsigned nTks = simTracks.size();
0163
0164   // Empty event, do nothin'
0165   if (nVtx == 0)
0166     return;
0167
0168   // create a map associating geant particle id and position in the
0169   // event SimTrack vector
0170   for (unsigned it = 0; it < nTks; ++it) {
0171     geantToIndex_[simTracks[it].trackId()] = it;
0172     std::cout << " PizeroMCTruthFinder::fill it " << it << " simTracks[it].trackId() " << simTracks[it].trackId()
0173               << std::endl;
0174   }
0175 }