Reference Manual

Go to the documentation of this file.
00001 #include "RecoEgamma/EgammaMCTools/interface/ElectronMCTruthFinder.h"
00002 
00003 
00004 #include <algorithm>                                                          
00005 
00006 ElectronMCTruthFinder::ElectronMCTruthFinder() {
00007   
00008   
00009 }
00010 
00011 ElectronMCTruthFinder::~ElectronMCTruthFinder() 
00012 {
00013   
00014 }
00015 
00016 
00017 std::vector<ElectronMCTruth> ElectronMCTruthFinder::find(std::vector<SimTrack> theSimTracks, std::vector<SimVertex> theSimVertices ) {
00018   //std::cout << "  ElectronMCTruthFinder::find " << std::endl;
00019   
00020   
00021   
00022   
00023   std::vector<ElectronMCTruth> result;
00024   
00025   // Local variables  
00026   const int SINGLE=1;
00027   const int DOUBLE=2;
00028   const int PYTHIA=3;
00029   const int ELECTRON_FLAV=1;
00030   const int PIZERO_FLAV=2;
00031   const int PHOTON_FLAV=3;
00032   
00033   int ievtype=0;
00034   int ievflav=0;
00035   
00036   
00037   std::vector<SimTrack> electronTracks;
00038   SimVertex primVtx;   
00039   
00040   fill(theSimTracks,  theSimVertices);
00041   
00042   int iPV=-1;   
00043   int partType1=0;
00044   int partType2=0;
00045   std::vector<SimTrack>::iterator iFirstSimTk = theSimTracks.begin();
00046   if (  !(*iFirstSimTk).noVertex() ) {
00047     iPV =  (*iFirstSimTk).vertIndex();
00048     
00049     int vtxId =   (*iFirstSimTk).vertIndex();
00050     primVtx = theSimVertices[vtxId];
00051     
00052     partType1 = (*iFirstSimTk).type();
00053     
00054     
00055     //std::cout <<  " Primary vertex id " << iPV << " first track type " << (*iFirstSimTk).type() << std::endl;  
00056   } else {
00057     //std::cout << " First track has no vertex " << std::endl;
00058   }
00059   
00060   // CLHEP::HepLorentzVector primVtxPos= primVtx.position(); 
00061   math::XYZTLorentzVectorD primVtxPos(primVtx.position().x(),
00062                                       primVtx.position().y(),
00063                                       primVtx.position().z(),
00064                                       primVtx.position().e());
00065   
00066   // Look at a second track
00067   iFirstSimTk++;
00068   if ( iFirstSimTk!=  theSimTracks.end() ) {
00069     
00070     if (  (*iFirstSimTk).vertIndex() == iPV) {
00071       partType2 = (*iFirstSimTk).type();  
00072   
00073       //std::cout <<  " second track type " << (*iFirstSimTk).type() << " vertex " <<  (*iFirstSimTk).vertIndex() << std::endl;  
00074       
00075     } else {
00076       //std::cout << " Only one kine track from Primary Vertex " << std::endl;
00077     }
00078   }
00079   
00080   //std::cout << " Loop over all particles " << std::endl;
00081   
00082   int npv=0;
00083   for (std::vector<SimTrack>::iterator iSimTk = theSimTracks.begin(); iSimTk != theSimTracks.end(); ++iSimTk){
00084     if (  (*iSimTk).noVertex() ) continue;
00085     
00086     int vertexId = (*iSimTk).vertIndex();
00087     SimVertex vertex = theSimVertices[vertexId];
00088     
00089     //std::cout << " Particle type " <<  (*iSimTk).type() << " Sim Track ID " << (*iSimTk).trackId() << " momentum " << (*iSimTk).momentum() <<  " vertex position " << vertex.position() << " vertex ID " << vertexId  << std::endl;  
00090     if ( (*iSimTk).vertIndex() == iPV ) {
00091       npv++;
00092       if ( std::abs((*iSimTk).type() ) == 11) {
00093         //std::cout << " Found a primary electron with ID  " << (*iSimTk).trackId() << " momentum " << (*iSimTk).momentum() <<  std::endl;
00094         electronTracks.push_back( *iSimTk );
00095       }
00096     }
00097   }
00098   //std::cout << " There are " << npv << " particles originating in the PV " << std::endl;
00099     
00100   
00101   if(npv > 4) {
00102     ievtype = PYTHIA;
00103   } else if(npv == 1) {
00104     if( std::abs(partType1) == 11 ) {
00105       ievtype = SINGLE; ievflav = ELECTRON_FLAV;
00106     } else if(partType1 == 111) {
00107       ievtype = SINGLE; ievflav = PIZERO_FLAV;
00108     } else if(partType1 == 22) {
00109       ievtype = SINGLE; ievflav = PHOTON_FLAV;
00110     }
00111   } else if(npv == 2) {
00112     if (  std::abs(partType1) == 11 && std::abs(partType2) == 11 ) {
00113       ievtype = DOUBLE; ievflav = ELECTRON_FLAV;
00114     } else if(partType1 == 111 && partType2 == 111)   {
00115       ievtype = DOUBLE; ievflav = PIZERO_FLAV;
00116     } else if(partType1 == 22 && partType2 == 22)   {
00117       ievtype = DOUBLE; ievflav = PHOTON_FLAV;
00118     }
00119   }
00120 
00121 
00122 
00124   std::vector<CLHEP::Hep3Vector> bremPos;  
00125   std::vector<CLHEP::HepLorentzVector> pBrem;
00126   std::vector<float> xBrem;
00127  
00128   
00129   for (std::vector<SimTrack>::iterator iEleTk = electronTracks.begin(); iEleTk != electronTracks.end(); ++iEleTk){
00130     //std::cout << " Looping on the primary electron pt  " << std::sqrt((*iEleTk).momentum().perp2()) << " electron track ID " << (*iEleTk).trackId() << std::endl;
00131     
00132     
00133     
00134     SimTrack trLast =(*iEleTk); 
00135     unsigned int eleId = (*iEleTk).trackId();
00136     float remainingEnergy =trLast.momentum().e();
00137 //    CLHEP::HepLorentzVector motherMomentum = (*iEleTk).momentum();
00138 //    CLHEP::HepLorentzVector primEleMom = (*iEleTk).momentum();
00139     math::XYZTLorentzVectorD motherMomentum((*iEleTk).momentum().x(),
00140                                             (*iEleTk).momentum().y(),
00141                                             (*iEleTk).momentum().z(),
00142                                             (*iEleTk).momentum().e());
00143     math::XYZTLorentzVectorD primEleMom(motherMomentum);
00144     int eleVtxIndex= (*iEleTk).vertIndex();
00145     
00146     bremPos.clear();
00147     pBrem.clear();
00148     xBrem.clear();     
00149 
00150    
00151     for (std::vector<SimTrack>::iterator iSimTk = theSimTracks.begin(); iSimTk != theSimTracks.end(); ++iSimTk){
00152         
00153       if (  (*iSimTk).noVertex() )                    continue;
00154       if ( (*iSimTk).vertIndex() == iPV )             continue;
00155 
00156       //std::cout << " (*iEleTk)->trackId() " << (*iEleTk).trackId() << " (*iEleTk)->vertIndex() "<< (*iEleTk).vertIndex()  << " (*iSimTk).vertIndex() "  <<  (*iSimTk).vertIndex() << " (*iSimTk).type() " <<   (*iSimTk).type() << " (*iSimTk).trackId() " << (*iSimTk).trackId() << std::endl;
00157       
00158       int vertexId1 = (*iSimTk).vertIndex();
00159       SimVertex vertex1 = theSimVertices[vertexId1];
00160       int vertexId2 = trLast.vertIndex();
00161       SimVertex vertex2 = theSimVertices[vertexId2];
00162       
00163       
00164       int motherId=-1;
00165       
00166       if(  (  vertexId1 ==  vertexId2 ) && ( (*iSimTk).type() == (*iEleTk).type() ) && trLast.type() == 22   ) {
00167         //std::cout << " Here a e/gamma brem vertex " << std::endl;
00168         
00169         //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;
00170         
00171         //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;
00172         
00173         //std::cout << " Electron pt " << std::sqrt((*iSimTk).momentum().perp2()) << " photon pt " <<  std::sqrt(trLast.momentum().perp2()) << "Mother electron pt " <<  sqrt(motherMomentum.perp2()) << std::endl;
00174         //std::cout << " eleId " << eleId << std::endl;
00175         float eLoss = remainingEnergy - ( (*iSimTk).momentum() + trLast.momentum()).e();
00176         //std::cout << " eLoss " << eLoss << std::endl;              
00177 
00178 
00179         if ( vertex1.parentIndex()  ) {
00180           
00181           unsigned  motherGeantId = vertex1.parentIndex(); 
00182           std::map<unsigned, unsigned >::iterator association = geantToIndex_.find( motherGeantId );
00183           if(association != geantToIndex_.end() )
00184             motherId = association->second;
00185           
00186           //int motherType = motherId == -1 ? 0 : theSimTracks[motherId].type();
00187           //std::cout << " Parent to this vertex   motherId " << motherId << " mother type " <<  motherType << " Sim track ID " <<  theSimTracks[motherId].trackId() << std::endl; 
00188           if ( theSimTracks[motherId].trackId() == eleId ) {
00189             
00190             //std::cout << "  ***** Found the Initial Mother Electron ****   theSimTracks[motherId].trackId() " <<  theSimTracks[motherId].trackId() << " eleId " <<  eleId << std::endl;
00191             eleId= (*iSimTk).trackId();
00192             remainingEnergy =   (*iSimTk).momentum().e();
00193             motherMomentum = (*iSimTk).momentum();
00194             
00195             
00196             pBrem.push_back( CLHEP::HepLorentzVector(trLast.momentum().px(),trLast.momentum().py(),
00197                                               trLast.momentum().pz(),trLast.momentum().e()) );
00198             bremPos.push_back( CLHEP::HepLorentzVector(vertex1.position().x(),vertex1.position().y(),
00199                                                 vertex1.position().z(),vertex1.position().t()) );
00200             xBrem.push_back(eLoss);
00201             
00202           }
00203           
00204           
00205           
00206           
00207         } else {
00208           //std::cout << " This vertex has no parent tracks " <<  std::endl;
00209         }
00210         
00211       }
00212       trLast=(*iSimTk);
00213       
00214     } // End loop over all SimTracks 
00215     //std::cout << " Going to build the ElectronMCTruth: pBrem size " << pBrem.size() << std::endl;
00217     CLHEP::HepLorentzVector tmpEleMom(primEleMom.px(),primEleMom.py(),
00218                                primEleMom.pz(),primEleMom.e() ) ;
00219     CLHEP::HepLorentzVector tmpVtxPos(primVtxPos.x(),primVtxPos.y(),primVtxPos.z(),primVtxPos.t());
00220     result.push_back ( ElectronMCTruth( tmpEleMom, eleVtxIndex,  bremPos, pBrem, xBrem,  tmpVtxPos,(*iEleTk)  )  ) ;
00221 
00222     
00223   } // End loop over primary electrons 
00224   
00225     
00226    
00227    return result;
00228 }
00229 
00230 
00231 
00232 void ElectronMCTruthFinder::fill(std::vector<SimTrack>& simTracks, 
00233                                  std::vector<SimVertex>& simVertices ) {
00234   //std::cout << "  ElectronMCTruthFinder::fill " << std::endl;
00235 
00236   unsigned nVtx = simVertices.size();
00237   unsigned nTks = simTracks.size();
00238 
00239   // Empty event, do nothin'
00240   if ( nVtx == 0 ) return;
00241 
00242   // create a map associating geant particle id and position in the 
00243   // event SimTrack vector
00244   for( unsigned it=0; it<nTks; ++it ) {
00245     geantToIndex_[ simTracks[it].trackId() ] = it;
00246     //std::cout << " ElectronMCTruthFinder::fill it " << it << " simTracks[it].trackId() " <<  simTracks[it].trackId() << std::endl;
00247  
00248   }  
00249 
00250 
00251 }
CMSSW_4_4_3_patch1/src/RecoEgamma/EgammaMCTools/src/ElectronMCTruthFinder.cc
