Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:24:57 

0001 #include "RecoEgamma/EgammaMCTools/interface/PhotonMCTruthFinder.h"
0002 #include "RecoEgamma/EgammaMCTools/interface/PhotonMCTruth.h"
0003 #include "RecoEgamma/EgammaMCTools/interface/ElectronMCTruth.h"
0004 //
0005 //
0006 #include "SimDataFormats/CrossingFrame/interface/CrossingFrame.h"
0007 #include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
0008 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
0009 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h"
0010 
0011 #include <algorithm>
0012 
0013 PhotonMCTruthFinder::PhotonMCTruthFinder() {
0014   //std::cout << " PhotonMCTruthFinder CTOR " << std::endl;
0015 }
0016 
0017 std::vector<PhotonMCTruth> PhotonMCTruthFinder::find(const std::vector<SimTrack>& theSimTracks,
0018                                                      const std::vector<SimVertex>& theSimVertices) {
0019   //  std::cout << "  PhotonMCTruthFinder::find " << std::endl;
0020 
0021   std::vector<PhotonMCTruth> result;
0022 
0023   //const float pi = 3.141592653592;
0024   //const float twopi=2*pi;
0025 
0026   // Local variables
0027   const int SINGLE = 1;
0028   const int DOUBLE = 2;
0029   const int PYTHIA = 3;
0030   const int ELECTRON_FLAV = 1;
0031   const int PIZERO_FLAV = 2;
0032   const int PHOTON_FLAV = 3;
0033 
0034   int ievtype = 0;
0035   int ievflav = 0;
0036 
0037   std::vector<SimTrack*> photonTracks;
0038 
0039   std::vector<SimTrack> trkFromConversion;
0040   std::vector<ElectronMCTruth> electronsFromConversions;
0041   SimVertex primVtx;
0042 
0043   fill(theSimTracks, theSimVertices);
0044 
0045   //  std::cout << " After fill " << theSimTracks.size() << " " << theSimVertices.size() << std::endl;
0046   if (!theSimTracks.empty()) {
0047     int iPV = -1;
0048     int partType1 = 0;
0049     int partType2 = 0;
0050     std::vector<SimTrack>::const_iterator iFirstSimTk = theSimTracks.begin();
0051     if (!(*iFirstSimTk).noVertex()) {
0052       iPV = (*iFirstSimTk).vertIndex();
0053 
0054       int vtxId = (*iFirstSimTk).vertIndex();
0055       primVtx = theSimVertices[vtxId];
0056 
0057       partType1 = (*iFirstSimTk).type();
0058       //    std::cout <<  " Primary vertex id " << iPV << " first track type " << (*iFirstSimTk).type() << std::endl;
0059     } else {
0060       //std::cout << " First track has no vertex " << std::endl;
0061     }
0062 
0063     math::XYZTLorentzVectorD primVtxPos(
0064         primVtx.position().x(), primVtx.position().y(), primVtx.position().z(), primVtx.position().e());
0065 
0066     // Look at a second track
0067     iFirstSimTk++;
0068     if (iFirstSimTk != theSimTracks.end()) {
0069       if ((*iFirstSimTk).vertIndex() == iPV) {
0070         partType2 = (*iFirstSimTk).type();
0071         //      std::cout <<  " second track type " << (*iFirstSimTk).type() << " vertex " <<  (*iFirstSimTk).vertIndex() << std::endl;
0072 
0073       } else {
0074         // std::cout << " Only one kine track from Primary Vertex " << std::endl;
0075       }
0076     }
0077 
0078     //std::cout << " Loop over all particles " << std::endl;
0079 
0080     int npv = 0;
0081     //int iPho=0;
0082     //int iPizero=0;
0083     //   theSimTracks.reset();
0084     for (std::vector<SimTrack>::const_iterator iSimTk = theSimTracks.begin(); iSimTk != theSimTracks.end(); ++iSimTk) {
0085       if ((*iSimTk).noVertex())
0086         continue;
0087 
0088       //int vertexId = (*iSimTk).vertIndex();
0089       //SimVertex vertex = theSimVertices[vertexId];
0090 
0091       //    std::cout << " Particle type " <<  (*iSimTk).type() << " Sim Track ID " << (*iSimTk).trackId() << " momentum " << (*iSimTk).momentum() <<  " vertex position " << vertex.position() << " vertex index " << (*iSimTk).vertIndex() << std::endl;
0092       if ((*iSimTk).vertIndex() == iPV) {
0093         npv++;
0094         if ((*iSimTk).type() == 22) {
0095           //    std::cout << " Found a primary photon with ID  " << (*iSimTk).trackId() << " momentum " << (*iSimTk).momentum() <<  std::endl;
0096 
0097           photonTracks.push_back(&(const_cast<SimTrack&>(*iSimTk)));
0098         }
0099       }
0100     }
0101 
0102     //  std::cout << " There are " << npv << " particles originating in the PV " << std::endl;
0103 
0104     if (npv >= 3) {
0105       ievtype = PYTHIA;
0106     } else if (npv == 1) {
0107       if (std::abs(partType1) == 11) {
0108         ievtype = SINGLE;
0109         ievflav = ELECTRON_FLAV;
0110       } else if (partType1 == 111) {
0111         ievtype = SINGLE;
0112         ievflav = PIZERO_FLAV;
0113       } else if (partType1 == 22) {
0114         ievtype = SINGLE;
0115         ievflav = PHOTON_FLAV;
0116       }
0117     } else if (npv == 2) {
0118       if (std::abs(partType1) == 11 && std::abs(partType2) == 11) {
0119         ievtype = DOUBLE;
0120         ievflav = ELECTRON_FLAV;
0121       } else if (partType1 == 111 && partType2 == 111) {
0122         ievtype = DOUBLE;
0123         ievflav = PIZERO_FLAV;
0124       } else if (partType1 == 22 && partType2 == 22) {
0125         ievtype = DOUBLE;
0126         ievflav = PHOTON_FLAV;
0127       }
0128     }
0129 
0130     //////  Look into converted photons
0131 
0132     int isAconversion = 0;
0133     int phoMotherType = -1;
0134     int phoMotherVtxIndex = -1;
0135     int phoMotherId = -1;
0136     if (ievflav == PHOTON_FLAV || ievflav == PIZERO_FLAV || ievtype == PYTHIA) {
0137       //     std::cout << " It's a primary PHOTON or PIZERO or PYTHIA event with " << photonTracks.size() << " photons ievtype " << ievtype << " ievflav " << ievflav<<  std::endl;
0138 
0139       //    for (std::vector<SimTrack*>::iterator iPhoTk = photonTracks.begin(); iPhoTk != photonTracks.end(); ++iPhoTk){
0140       //      std::cout << " All gamma found from PV " << (*iPhoTk)->momentum() << " photon track ID " << (*iPhoTk)->trackId() << " vertex index " << (*iPhoTk)->vertIndex() << std::endl;
0141       //  }
0142 
0143       for (std::vector<SimTrack>::const_iterator iPhoTk = theSimTracks.begin(); iPhoTk != theSimTracks.end();
0144            ++iPhoTk) {
0145         trkFromConversion.clear();
0146         electronsFromConversions.clear();
0147 
0148         if ((*iPhoTk).type() != 22)
0149           continue;
0150         int photonVertexIndex = (*iPhoTk).vertIndex();
0151         int phoTrkId = (*iPhoTk).trackId();
0152         //std::cout << " Looping on gamma looking for conversions " << (*iPhoTk).momentum() << " photon track ID " << (*iPhoTk).trackId() << std::endl;
0153 
0154         // check who is his mother
0155         const SimVertex& vertex = theSimVertices[photonVertexIndex];
0156         phoMotherId = -1;
0157         if (vertex.parentIndex() != -1) {
0158           unsigned motherGeantId = vertex.parentIndex();
0159           std::map<unsigned, unsigned>::iterator association = geantToIndex_.find(motherGeantId);
0160           if (association != geantToIndex_.end())
0161             phoMotherId = association->second;
0162           phoMotherType = phoMotherId == -1 ? 0 : theSimTracks[phoMotherId].type();
0163 
0164           if (phoMotherType == 111 || phoMotherType == 221 || phoMotherType == 331) {
0165             //std::cout << " Parent to this vertex   motherId " << phoMotherId << " mother type " <<  phoMotherType << " Sim track ID " <<  theSimTracks[phoMotherId].trackId() << std::endl;
0166             //std::cout << " Son of a pizero or eta " << phoMotherType << std::endl;
0167           }
0168         }
0169 
0170         for (std::vector<SimTrack>::const_iterator iEleTk = theSimTracks.begin(); iEleTk != theSimTracks.end();
0171              ++iEleTk) {
0172           if ((*iEleTk).noVertex())
0173             continue;
0174           if ((*iEleTk).vertIndex() == iPV)
0175             continue;
0176           if (std::abs((*iEleTk).type()) != 11)
0177             continue;
0178 
0179           int vertexId = (*iEleTk).vertIndex();
0180           const SimVertex& vertex = theSimVertices[vertexId];
0181           int motherId = -1;
0182 
0183           //std::cout << " Secondary from photons particle type " << (*iEleTk).type() << " trackId " <<  (*iEleTk).trackId() << " vertex ID " << vertexId << std::endl;
0184           if (vertex.parentIndex() != -1) {
0185             unsigned motherGeantId = vertex.parentIndex();
0186             std::map<unsigned, unsigned>::iterator association = geantToIndex_.find(motherGeantId);
0187             if (association != geantToIndex_.end())
0188               motherId = association->second;
0189 
0190             //int motherType = motherId == -1 ? 0 : theSimTracks[motherId].type();
0191 
0192             //std::cout << " Parent to this vertex   motherId " << motherId << " mother type " <<  motherType << " Sim track ID " <<  theSimTracks[motherId].trackId() << std::endl;
0193 
0194             std::vector<CLHEP::Hep3Vector> bremPos;
0195             std::vector<CLHEP::HepLorentzVector> pBrem;
0196             std::vector<float> xBrem;
0197 
0198             if (theSimTracks[motherId].trackId() == (*iPhoTk).trackId()) {
0199               //std::cout << " Found the Mother Photon " << std::endl;
0200               /// find truth about this electron and store it since it's from a converted photon
0201 
0202               trkFromConversion.push_back((*iEleTk));
0203               SimTrack trLast = (*iEleTk);
0204               float remainingEnergy = trLast.momentum().e();
0205               //HepLorentzVector primEleMom=(*iEleTk).momentum();
0206               //HepLorentzVector motherMomentum=(*iEleTk).momentum();
0207               math::XYZTLorentzVectorD primEleMom((*iEleTk).momentum().x(),
0208                                                   (*iEleTk).momentum().y(),
0209                                                   (*iEleTk).momentum().z(),
0210                                                   (*iEleTk).momentum().e());
0211               math::XYZTLorentzVectorD motherMomentum(primEleMom);
0212               unsigned int eleId = (*iEleTk).trackId();
0213               int eleVtxIndex = (*iEleTk).vertIndex();
0214 
0215               bremPos.clear();
0216               pBrem.clear();
0217               xBrem.clear();
0218 
0219               for (std::vector<SimTrack>::const_iterator iSimTk = theSimTracks.begin(); iSimTk != theSimTracks.end();
0220                    ++iSimTk) {
0221                 if ((*iSimTk).noVertex())
0222                   continue;
0223                 if ((*iSimTk).vertIndex() == iPV)
0224                   continue;
0225 
0226                 //std::cout << " (*iEleTk)->trackId() " << (*iEleTk).trackId() << " (*iEleTk)->vertIndex() "<< (*iEleTk).vertIndex()  << " (*iSimTk).vertIndex() "  <<  (*iSimTk).vertIndex() << " (*iSimTk).type() " <<   (*iSimTk).type() << " (*iSimTk).trackId() " << (*iSimTk).trackId() << std::endl;
0227 
0228                 int vertexId1 = (*iSimTk).vertIndex();
0229                 const SimVertex& vertex1 = theSimVertices[vertexId1];
0230                 int vertexId2 = trLast.vertIndex();
0231                 //SimVertex vertex2 = theSimVertices[vertexId2];
0232 
0233                 int motherId = -1;
0234 
0235                 if ((vertexId1 == vertexId2) && ((*iSimTk).type() == (*iEleTk).type()) && trLast.type() == 22) {
0236                   //std::cout << " Here a e/gamma brem vertex " << std::endl;
0237 
0238                   //std::cout << " Secondary from electron:  particle1  type " << (*iSimTk).type() << " trackId " <<  (*iSimTk).trackId() << " vertex ID " << vertexId1 << " vertex position " << sqrt(vertex1.position().perp2()) << " parent index "<< vertex1.parentIndex() << std::endl;
0239 
0240                   //std::cout << " Secondary from electron:  particle2  type " << trLast.type() << " trackId " <<  trLast.trackId() << " vertex ID " << vertexId2 << " vertex position " << sqrt(vertex2.position().perp2()) << " parent index " << vertex2.parentIndex() << std::endl;
0241 
0242                   //std::cout << " Electron pt " << sqrt((*iSimTk).momentum().perp2()) << " photon pt " << sqrt(trLast.momentum().perp2()) << " Mother electron pt " <<  sqrt(motherMomentum.perp2()) << std::endl;
0243                   //std::cout << " eleId " << eleId << std::endl;
0244                   float eLoss = remainingEnergy - ((*iSimTk).momentum() + trLast.momentum()).e();
0245                   //std::cout << " eLoss " << eLoss << std::endl;
0246 
0247                   if (vertex1.parentIndex() != -1) {
0248                     unsigned motherGeantId = vertex1.parentIndex();
0249                     std::map<unsigned, unsigned>::iterator association = geantToIndex_.find(motherGeantId);
0250                     if (association != geantToIndex_.end())
0251                       motherId = association->second;
0252 
0253                     //int motherType = motherId == -1 ? 0 : theSimTracks[motherId].type();
0254                     //std::cout << " Parent to this vertex   motherId " << motherId << " mother type " <<  motherType << " Sim track ID " <<  theSimTracks[motherId].trackId() << std::endl;
0255                     if (theSimTracks[motherId].trackId() == eleId) {
0256                       //std::cout << "  ***** Found the Initial Mother Electron ****   theSimTracks[motherId].trackId() " <<  theSimTracks[motherId].trackId() << " eleId " <<  eleId << std::endl;
0257                       eleId = (*iSimTk).trackId();
0258                       remainingEnergy = (*iSimTk).momentum().e();
0259                       motherMomentum = (*iSimTk).momentum();
0260 
0261                       pBrem.push_back(CLHEP::HepLorentzVector(trLast.momentum().px(),
0262                                                               trLast.momentum().py(),
0263                                                               trLast.momentum().pz(),
0264                                                               trLast.momentum().e()));
0265                       bremPos.push_back(CLHEP::HepLorentzVector(vertex1.position().x(),
0266                                                                 vertex1.position().y(),
0267                                                                 vertex1.position().z(),
0268                                                                 vertex1.position().t()));
0269                       xBrem.push_back(eLoss);
0270                     }
0271 
0272                   } else {
0273                     //std::cout << " This vertex has no parent tracks " <<  std::endl;
0274                   }
0275                 }
0276                 trLast = (*iSimTk);
0277 
0278               }  // End loop over all SimTracks
0279               //std::cout << " Going to build the ElectronMCTruth for this electron from converted photon: pBrem size " << pBrem.size() << std::endl;
0280               /// here fill the electron
0281 
0282               CLHEP::HepLorentzVector tmpEleMom(primEleMom.px(), primEleMom.py(), primEleMom.pz(), primEleMom.e());
0283               CLHEP::HepLorentzVector tmpVtxPos(primVtxPos.x(), primVtxPos.y(), primVtxPos.z(), primVtxPos.t());
0284               electronsFromConversions.push_back(ElectronMCTruth(
0285                   tmpEleMom, eleVtxIndex, bremPos, pBrem, xBrem, tmpVtxPos, const_cast<SimTrack&>(*iEleTk)));
0286             }  //// Electron from conversion found
0287 
0288           } else {
0289             //std::cout << " This vertex has no parent tracks " <<  std::endl;
0290           }
0291 
0292         }  // End of loop over the SimTracks
0293 
0294         //std::cout << " DEBUG trkFromConversion.size() " << trkFromConversion.size() << " electronsFromConversions.size() " << electronsFromConversions.size() << std::endl;
0295 
0296         math::XYZTLorentzVectorD motherVtxPosition(0., 0., 0., 0.);
0297         CLHEP::HepLorentzVector phoMotherMom(0., 0., 0., 0.);
0298         CLHEP::HepLorentzVector phoMotherVtx(0., 0., 0., 0.);
0299 
0300         if (phoMotherId >= 0) {
0301           phoMotherVtxIndex = theSimTracks[phoMotherId].vertIndex();
0302           const SimVertex& motherVtx = theSimVertices[phoMotherVtxIndex];
0303           motherVtxPosition = math::XYZTLorentzVectorD(
0304               motherVtx.position().x(), motherVtx.position().y(), motherVtx.position().z(), motherVtx.position().e());
0305 
0306           phoMotherMom.setPx(theSimTracks[phoMotherId].momentum().x());
0307           phoMotherMom.setPy(theSimTracks[phoMotherId].momentum().y());
0308           phoMotherMom.setPz(theSimTracks[phoMotherId].momentum().z());
0309           phoMotherMom.setE(theSimTracks[phoMotherId].momentum().t());
0310           // std::cout << " PhotonMCTruthFinder mother " << phoMotherId << " type " << phoMotherType << " Momentum" <<  phoMotherMom.et() << std::endl;
0311 
0312           phoMotherVtx.setX(motherVtxPosition.x());
0313           phoMotherVtx.setY(motherVtxPosition.y());
0314           phoMotherVtx.setZ(motherVtxPosition.z());
0315           phoMotherVtx.setT(motherVtxPosition.t());
0316         }
0317 
0318         if (!electronsFromConversions.empty()) {
0319           // if ( trkFromConversion.size() > 0 ) {
0320           isAconversion = 1;
0321           //std::cout  << " CONVERTED photon " <<   "\n";
0322 
0323           //int convVtxId =  trkFromConversion[0].vertIndex();
0324           int convVtxId = electronsFromConversions[0].vertexInd();
0325           const SimVertex& convVtx = theSimVertices[convVtxId];
0326           // CLHEP::HepLorentzVector vtxPosition = convVtx.position();
0327           math::XYZTLorentzVectorD vtxPosition(
0328               convVtx.position().x(), convVtx.position().y(), convVtx.position().z(), convVtx.position().e());
0329 
0330           //result.push_back( PhotonMCTruth(isAconversion, (*iPhoTk).momentum(), photonVertexIndex, phoTrkId, vtxPosition,   primVtx.position(), trkFromConversion ));
0331           CLHEP::HepLorentzVector tmpPhoMom((*iPhoTk).momentum().px(),
0332                                             (*iPhoTk).momentum().py(),
0333                                             (*iPhoTk).momentum().pz(),
0334                                             (*iPhoTk).momentum().e());
0335           CLHEP::HepLorentzVector tmpVertex(vtxPosition.x(), vtxPosition.y(), vtxPosition.z(), vtxPosition.t());
0336           CLHEP::HepLorentzVector tmpPrimVtx(primVtxPos.x(), primVtxPos.y(), primVtxPos.z(), primVtxPos.t());
0337 
0338           result.push_back(PhotonMCTruth(isAconversion,
0339                                          tmpPhoMom,
0340                                          photonVertexIndex,
0341                                          phoTrkId,
0342                                          phoMotherType,
0343                                          phoMotherMom,
0344                                          phoMotherVtx,
0345                                          tmpVertex,
0346                                          tmpPrimVtx,
0347                                          electronsFromConversions));
0348 
0349         } else {
0350           isAconversion = 0;
0351           //std::cout  << " UNCONVERTED photon " <<   "\n";
0352           CLHEP::HepLorentzVector vtxPosition(0., 0., 0., 0.);
0353           CLHEP::HepLorentzVector tmpPhoMom((*iPhoTk).momentum().px(),
0354                                             (*iPhoTk).momentum().py(),
0355                                             (*iPhoTk).momentum().pz(),
0356                                             (*iPhoTk).momentum().e());
0357           CLHEP::HepLorentzVector tmpPrimVtx(primVtxPos.x(), primVtxPos.y(), primVtxPos.z(), primVtxPos.t());
0358           //     result.push_back( PhotonMCTruth(isAconversion, (*iPhoTk).momentum(),  photonVertexIndex, phoTrkId, vtxPosition,   primVtx.position(), trkFromConversion ));
0359           result.push_back(PhotonMCTruth(isAconversion,
0360                                          tmpPhoMom,
0361                                          photonVertexIndex,
0362                                          phoTrkId,
0363                                          phoMotherType,
0364                                          phoMotherMom,
0365                                          phoMotherVtx,
0366                                          vtxPosition,
0367                                          tmpPrimVtx,
0368                                          electronsFromConversions));
0369         }
0370 
0371       }  // End loop over the primary photons
0372 
0373     }  // Event with one or two photons
0374 
0375     //std::cout << "  PhotonMCTruthFinder photon size " << result.size() << std::endl;
0376   }
0377 
0378   return result;
0379 }
0380 
0381 void PhotonMCTruthFinder::fill(const std::vector<SimTrack>& simTracks, const std::vector<SimVertex>& simVertices) {
0382   //  std::cout << "  PhotonMCTruthFinder::fill " << std::endl;
0383 
0384   // Watch out there ! A SimVertex is in mm (stupid),
0385 
0386   unsigned nVtx = simVertices.size();
0387   unsigned nTks = simTracks.size();
0388 
0389   //  std::cout << "  PhotonMCTruthFinder::fill " << nVtx << " " << nTks << std::endl;
0390 
0391   // Empty event, do nothin'
0392   if (nVtx == 0)
0393     return;
0394 
0395   // create a map associating geant particle id and position in the
0396   // event SimTrack vector
0397   for (unsigned it = 0; it < nTks; ++it) {
0398     geantToIndex_[simTracks[it].trackId()] = it;
0399     //    std::cout << " PhotonMCTruthFinder::fill it " << it << " simTracks[it].trackId() " <<  simTracks[it].trackId() << std::endl;
0400   }
0401 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!