Reference Manual

Go to the documentation of this file.
00001 #include "RecoEgamma/EgammaMCTools/interface/PhotonMCTruthFinder.h"
00002 #include "RecoEgamma/EgammaMCTools/interface/PhotonMCTruth.h"
00003 #include "RecoEgamma/EgammaMCTools/interface/ElectronMCTruth.h"
00004 //
00005 //
00006 #include "SimDataFormats/CrossingFrame/interface/CrossingFrame.h"
00007 #include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
00008 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
00009 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h"
00010 
00011 #include <algorithm>
00012 
00013 
00014 PhotonMCTruthFinder::PhotonMCTruthFinder( ) {
00015  //std::cout << " PhotonMCTruthFinder CTOR " << std::endl;
00016 
00017  
00018 }
00019 
00020 std::vector<PhotonMCTruth> PhotonMCTruthFinder::find(std::vector<SimTrack> theSimTracks, std::vector<SimVertex> theSimVertices ) {
00021   //  std::cout << "  PhotonMCTruthFinder::find " << std::endl;
00022 
00023   std::vector<PhotonMCTruth> result;
00024 
00025   //const float pi = 3.141592653592;
00026   //const float twopi=2*pi;
00027 
00028 // Local variables  
00029   const int SINGLE=1;
00030   const int DOUBLE=2;
00031   const int PYTHIA=3;
00032   const int ELECTRON_FLAV=1;
00033   const int PIZERO_FLAV=2;
00034   const int PHOTON_FLAV=3;
00035   
00036   int ievtype=0;
00037   int ievflav=0;
00038   
00039   
00040   std::vector<SimTrack*> photonTracks;
00041   
00042   std::vector<SimTrack> trkFromConversion;
00043   std::vector<ElectronMCTruth> electronsFromConversions;
00044   SimVertex primVtx;   
00045 
00046   
00047   fill(theSimTracks,  theSimVertices);
00048 
00049   //  std::cout << " After fill " << theSimTracks.size() << " " << theSimVertices.size() << std::endl;
00050   if (    theSimTracks.size() != 0 ) {
00051   
00052   int iPV=-1;   
00053   int partType1=0;
00054   int partType2=0;
00055   std::vector<SimTrack>::iterator iFirstSimTk = theSimTracks.begin();
00056   if (  !(*iFirstSimTk).noVertex() ) {
00057     iPV =  (*iFirstSimTk).vertIndex();
00058     
00059     int vtxId =   (*iFirstSimTk).vertIndex();
00060     primVtx = theSimVertices[vtxId];
00061     
00062     
00063     
00064     partType1 = (*iFirstSimTk).type();
00065     //    std::cout <<  " Primary vertex id " << iPV << " first track type " << (*iFirstSimTk).type() << std::endl;  
00066   } else {
00067     //std::cout << " First track has no vertex " << std::endl;
00068     
00069   }
00070   
00071 
00072 
00073   math::XYZTLorentzVectorD primVtxPos(primVtx.position().x(),
00074                                        primVtx.position().y(),
00075                                        primVtx.position().z(),
00076                                        primVtx.position().e());           
00077  
00078 
00079   // Look at a second track
00080   iFirstSimTk++;
00081   if ( iFirstSimTk!=  theSimTracks.end() ) {
00082     
00083     if (  (*iFirstSimTk).vertIndex() == iPV) {
00084       partType2 = (*iFirstSimTk).type();  
00085       //      std::cout <<  " second track type " << (*iFirstSimTk).type() << " vertex " <<  (*iFirstSimTk).vertIndex() << std::endl;  
00086       
00087     } else {
00088       // std::cout << " Only one kine track from Primary Vertex " << std::endl;
00089     }
00090   }
00091   
00092   //std::cout << " Loop over all particles " << std::endl;
00093   
00094   int npv=0;
00095   //int iPho=0;
00096   //int iPizero=0;
00097   //   theSimTracks.reset();
00098   for (std::vector<SimTrack>::iterator iSimTk = theSimTracks.begin(); iSimTk != theSimTracks.end(); ++iSimTk){
00099     if (  (*iSimTk).noVertex() ) continue;
00100 
00101 
00102     int vertexId = (*iSimTk).vertIndex();
00103     SimVertex vertex = theSimVertices[vertexId];
00104  
00105     //    std::cout << " Particle type " <<  (*iSimTk).type() << " Sim Track ID " << (*iSimTk).trackId() << " momentum " << (*iSimTk).momentum() <<  " vertex position " << vertex.position() << " vertex index " << (*iSimTk).vertIndex() << std::endl;  
00106     if ( (*iSimTk).vertIndex() == iPV ) {
00107       npv++;
00108       if ( (*iSimTk).type() == 22) {
00109         //      std::cout << " Found a primary photon with ID  " << (*iSimTk).trackId() << " momentum " << (*iSimTk).momentum() <<  std::endl; 
00110 
00111         
00112         photonTracks.push_back( &(*iSimTk) );
00113 
00114       } 
00115       
00116     }
00117     
00118   } 
00119   
00120   //  std::cout << " There are " << npv << " particles originating in the PV " << std::endl;
00121   
00122    if(npv >= 3) {
00123      ievtype = PYTHIA;
00124    } else if(npv == 1) {
00125      if( std::abs(partType1) == 11 ) {
00126        ievtype = SINGLE; ievflav = ELECTRON_FLAV;
00127      } else if(partType1 == 111) {
00128        ievtype = SINGLE; ievflav = PIZERO_FLAV;
00129      } else if(partType1 == 22) {
00130        ievtype = SINGLE; ievflav = PHOTON_FLAV;
00131      }
00132    } else if(npv == 2) {
00133      if (  std::abs(partType1) == 11 && std::abs(partType2) == 11 ) {
00134        ievtype = DOUBLE; ievflav = ELECTRON_FLAV;
00135      } else if(partType1 == 111 && partType2 == 111)   {
00136        ievtype = DOUBLE; ievflav = PIZERO_FLAV;
00137      } else if(partType1 == 22 && partType2 == 22)   {
00138        ievtype = DOUBLE; ievflav = PHOTON_FLAV;
00139      }
00140    }      
00141    
00143 
00144    int isAconversion=0;   
00145    int phoMotherType=-1;
00146    int phoMotherVtxIndex=-1;
00147    int phoMotherId=-1;
00148    if( ievflav == PHOTON_FLAV || ievflav== PIZERO_FLAV || ievtype == PYTHIA ) {
00149 
00150      //     std::cout << " It's a primary PHOTON or PIZERO or PYTHIA event with " << photonTracks.size() << " photons ievtype " << ievtype << " ievflav " << ievflav<<  std::endl;
00151 
00152      //    for (std::vector<SimTrack*>::iterator iPhoTk = photonTracks.begin(); iPhoTk != photonTracks.end(); ++iPhoTk){
00153      //      std::cout << " All gamma found from PV " << (*iPhoTk)->momentum() << " photon track ID " << (*iPhoTk)->trackId() << " vertex index " << (*iPhoTk)->vertIndex() << std::endl;  
00154      //  }        
00155 
00156      for (std::vector<SimTrack>::iterator iPhoTk = theSimTracks.begin(); iPhoTk != theSimTracks.end(); ++iPhoTk){
00157 
00158        trkFromConversion.clear();           
00159        electronsFromConversions.clear();
00160 
00161        if ( (*iPhoTk).type() != 22 ) continue;
00162        int photonVertexIndex= (*iPhoTk).vertIndex();
00163        int phoTrkId= (*iPhoTk).trackId();
00164        //std::cout << " Looping on gamma looking for conversions " << (*iPhoTk).momentum() << " photon track ID " << (*iPhoTk).trackId() << std::endl;
00165      
00166        // check who is his mother
00167        SimVertex vertex = theSimVertices[ photonVertexIndex];
00168        phoMotherId=-1;
00169        if ( vertex.parentIndex() != -1 ) {
00170          
00171          unsigned  motherGeantId = vertex.parentIndex(); 
00172          std::map<unsigned, unsigned >::iterator association = geantToIndex_.find( motherGeantId );
00173          if(association != geantToIndex_.end() )
00174            phoMotherId = association->second;
00175          phoMotherType = phoMotherId == -1 ? 0 : theSimTracks[phoMotherId].type();
00176 
00177          
00178 
00179          if ( phoMotherType == 111 || phoMotherType == 221 || phoMotherType == 331 ) {
00180          
00181            //std::cout << " Parent to this vertex   motherId " << phoMotherId << " mother type " <<  phoMotherType << " Sim track ID " <<  theSimTracks[phoMotherId].trackId() << std::endl;
00182            //std::cout << " Son of a pizero or eta " << phoMotherType << std::endl;
00183          }
00184          
00185        }
00186 
00187 
00188 
00189        for (std::vector<SimTrack>::iterator iEleTk = theSimTracks.begin(); iEleTk != theSimTracks.end(); ++iEleTk){
00190          if (  (*iEleTk).noVertex() )                    continue;
00191          if ( (*iEleTk).vertIndex() == iPV )             continue; 
00192          if ( std::abs((*iEleTk).type()) != 11  )             continue;
00193 
00194          int vertexId = (*iEleTk).vertIndex();
00195          SimVertex vertex = theSimVertices[vertexId];
00196          int motherId=-1;
00197         
00198 
00199          //std::cout << " Secondary from photons particle type " << (*iEleTk).type() << " trackId " <<  (*iEleTk).trackId() << " vertex ID " << vertexId << std::endl;
00200          if ( vertex.parentIndex() != -1 ) {
00201 
00202            unsigned  motherGeantId = vertex.parentIndex(); 
00203            std::map<unsigned, unsigned >::iterator association = geantToIndex_.find( motherGeantId );
00204            if(association != geantToIndex_.end() )
00205              motherId = association->second;
00206            
00207            
00208            //int motherType = motherId == -1 ? 0 : theSimTracks[motherId].type();
00209            
00210            //std::cout << " Parent to this vertex   motherId " << motherId << " mother type " <<  motherType << " Sim track ID " <<  theSimTracks[motherId].trackId() << std::endl;
00211 
00212            std::vector<CLHEP::Hep3Vector> bremPos;  
00213            std::vector<CLHEP::HepLorentzVector> pBrem;
00214            std::vector<float> xBrem;
00215            
00216            if ( theSimTracks[motherId].trackId() == (*iPhoTk).trackId() ) {
00217              //std::cout << " Found the Mother Photon " << std::endl;
00219 
00220 
00221 
00222              trkFromConversion.push_back( (*iEleTk ) );
00223              SimTrack trLast =(*iEleTk); 
00224              float remainingEnergy =trLast.momentum().e();
00225              //HepLorentzVector primEleMom=(*iEleTk).momentum();  
00226              //HepLorentzVector motherMomentum=(*iEleTk).momentum();  
00227              math::XYZTLorentzVectorD primEleMom((*iEleTk).momentum().x(),
00228                                              (*iEleTk).momentum().y(),
00229                                              (*iEleTk).momentum().z(),
00230                                              (*iEleTk).momentum().e());  
00231              math::XYZTLorentzVectorD motherMomentum(primEleMom);  
00232              unsigned int eleId = (*iEleTk).trackId();     
00233              int eleVtxIndex= (*iEleTk).vertIndex();
00234            
00235     
00236              bremPos.clear();
00237              pBrem.clear();
00238              xBrem.clear();   
00239                      
00240              for (std::vector<SimTrack>::iterator iSimTk = theSimTracks.begin(); iSimTk != theSimTracks.end(); ++iSimTk){
00241                
00242                if (  (*iSimTk).noVertex() )                    continue;
00243                if ( (*iSimTk).vertIndex() == iPV )             continue;
00244                
00245                //std::cout << " (*iEleTk)->trackId() " << (*iEleTk).trackId() << " (*iEleTk)->vertIndex() "<< (*iEleTk).vertIndex()  << " (*iSimTk).vertIndex() "  <<  (*iSimTk).vertIndex() << " (*iSimTk).type() " <<   (*iSimTk).type() << " (*iSimTk).trackId() " << (*iSimTk).trackId() << std::endl;
00246                
00247                int vertexId1 = (*iSimTk).vertIndex();
00248                SimVertex vertex1 = theSimVertices[vertexId1];
00249                int vertexId2 = trLast.vertIndex();
00250                SimVertex vertex2 = theSimVertices[vertexId2];
00251                
00252                
00253                int motherId=-1;
00254                
00255                if(  (  vertexId1 ==  vertexId2 ) && ( (*iSimTk).type() == (*iEleTk).type() ) && trLast.type() == 22   ) {
00256                  //std::cout << " Here a e/gamma brem vertex " << std::endl;
00257                  
00258                  //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;
00259                  
00260                  //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;
00261                  
00262                  //std::cout << " Electron pt " << sqrt((*iSimTk).momentum().perp2()) << " photon pt " << sqrt(trLast.momentum().perp2()) << " Mother electron pt " <<  sqrt(motherMomentum.perp2()) << std::endl;
00263                  //std::cout << " eleId " << eleId << std::endl;
00264                  float eLoss = remainingEnergy - ( (*iSimTk).momentum() + trLast.momentum()).e();
00265                  //std::cout << " eLoss " << eLoss << std::endl;              
00266                  
00267                  
00268                  if ( vertex1.parentIndex() != -1  ) {
00269                    
00270                    unsigned  motherGeantId = vertex1.parentIndex(); 
00271                    std::map<unsigned, unsigned >::iterator association = geantToIndex_.find( motherGeantId );
00272                    if(association != geantToIndex_.end() )
00273                      motherId = association->second;
00274                    
00275                    //int motherType = motherId == -1 ? 0 : theSimTracks[motherId].type();
00276                    //std::cout << " Parent to this vertex   motherId " << motherId << " mother type " <<  motherType << " Sim track ID " <<  theSimTracks[motherId].trackId() << std::endl; 
00277                    if (theSimTracks[motherId].trackId() == eleId ) {
00278                      
00279                      //std::cout << "  ***** Found the Initial Mother Electron ****   theSimTracks[motherId].trackId() " <<  theSimTracks[motherId].trackId() << " eleId " <<  eleId << std::endl;
00280                      eleId= (*iSimTk).trackId();
00281                      remainingEnergy =   (*iSimTk).momentum().e();
00282                      motherMomentum = (*iSimTk).momentum();
00283                      
00284                      
00285                      pBrem.push_back( CLHEP::HepLorentzVector(trLast.momentum().px(),
00286                                                        trLast.momentum().py(),
00287                                                        trLast.momentum().pz(),
00288                                                        trLast.momentum().e()) );
00289                      bremPos.push_back( CLHEP::HepLorentzVector(vertex1.position().x(),
00290                                                          vertex1.position().y(),
00291                                                          vertex1.position().z(),
00292                                                          vertex1.position().t()) );
00293                      xBrem.push_back(eLoss);
00294                      
00295                    }
00296                    
00297                    
00298                    
00299                    
00300                  } else {
00301                    //std::cout << " This vertex has no parent tracks " <<  std::endl;
00302                  }
00303                  
00304                }
00305                trLast=(*iSimTk);
00306                
00307              } // End loop over all SimTracks 
00308              //std::cout << " Going to build the ElectronMCTruth for this electron from converted photon: pBrem size " << pBrem.size() << std::endl;
00310 
00311              CLHEP::HepLorentzVector tmpEleMom(primEleMom.px(),primEleMom.py(),
00312                                         primEleMom.pz(),primEleMom.e() );
00313              CLHEP::HepLorentzVector tmpVtxPos(primVtxPos.x(),primVtxPos.y(),
00314                                         primVtxPos.z(),primVtxPos.t() );
00315              electronsFromConversions.push_back ( 
00316                 ElectronMCTruth( tmpEleMom, eleVtxIndex,  bremPos, pBrem, 
00317                                  xBrem,  tmpVtxPos , (*iEleTk)  )  ) ;
00318            }   
00319  
00320            
00321 
00322 
00323          } else {
00324            //std::cout << " This vertex has no parent tracks " <<  std::endl;
00325          }
00326          
00327          
00328        } // End of loop over the SimTracks      
00329        
00330        //std::cout << " DEBUG trkFromConversion.size() " << trkFromConversion.size() << " electronsFromConversions.size() " << electronsFromConversions.size() << std::endl;
00331 
00332        math::XYZTLorentzVectorD motherVtxPosition(0.,0.,0.,0.);
00333        CLHEP::HepLorentzVector phoMotherMom(0.,0.,0.,0.);
00334        CLHEP::HepLorentzVector phoMotherVtx(0.,0.,0.,0.); 
00335 
00336        if ( phoMotherId >= 0) {
00337         
00338          phoMotherVtxIndex = theSimTracks[phoMotherId].vertIndex();
00339          SimVertex motherVtx = theSimVertices[ phoMotherVtxIndex];
00340          motherVtxPosition =math::XYZTLorentzVectorD (motherVtx.position().x(),
00341                                                       motherVtx.position().y(),
00342                                                       motherVtx.position().z(),
00343                                                       motherVtx.position().e());
00344          
00345          phoMotherMom.setPx( theSimTracks[phoMotherId].momentum().x());
00346          phoMotherMom.setPy( theSimTracks[phoMotherId].momentum().y());
00347          phoMotherMom.setPz( theSimTracks[phoMotherId].momentum().z() );
00348          phoMotherMom.setE( theSimTracks[phoMotherId].momentum().t());
00349          // std::cout << " PhotonMCTruthFinder mother " << phoMotherId << " type " << phoMotherType << " Momentum" <<  phoMotherMom.et() << std::endl;  
00350  
00351          phoMotherVtx.setX ( motherVtxPosition.x());
00352          phoMotherVtx.setY ( motherVtxPosition.y());
00353          phoMotherVtx.setZ ( motherVtxPosition.z());
00354          phoMotherVtx.setT ( motherVtxPosition.t());
00355        
00356        }
00357 
00358 
00359        if ( electronsFromConversions.size() > 0 ) {
00360        // if ( trkFromConversion.size() > 0 ) {
00361          isAconversion=1;
00362          //std::cout  << " CONVERTED photon " <<   "\n";    
00363 
00364          //int convVtxId =  trkFromConversion[0].vertIndex();
00365          int convVtxId =electronsFromConversions[0].vertexInd();
00366          SimVertex convVtx = theSimVertices[convVtxId];
00367          // CLHEP::HepLorentzVector vtxPosition = convVtx.position();
00368          math::XYZTLorentzVectorD vtxPosition(convVtx.position().x(),
00369                                           convVtx.position().y(),
00370                                           convVtx.position().z(),
00371                                           convVtx.position().e());
00372          
00373            
00374          //result.push_back( PhotonMCTruth(isAconversion, (*iPhoTk).momentum(), photonVertexIndex, phoTrkId, vtxPosition,   primVtx.position(), trkFromConversion ));
00375          CLHEP::HepLorentzVector tmpPhoMom( (*iPhoTk).momentum().px(), (*iPhoTk).momentum().py(),
00376                                      (*iPhoTk).momentum().pz(), (*iPhoTk).momentum().e() ) ;
00377          CLHEP::HepLorentzVector tmpVertex( vtxPosition.x(), vtxPosition.y(), vtxPosition.z(), vtxPosition.t() );
00378          CLHEP::HepLorentzVector tmpPrimVtx( primVtxPos.x(), primVtxPos.y(), primVtxPos.z(), primVtxPos.t() ) ;
00379 
00380 
00381 
00382          result.push_back( PhotonMCTruth(isAconversion, tmpPhoMom, photonVertexIndex, phoTrkId, phoMotherType,phoMotherMom, phoMotherVtx, tmpVertex,  
00383          tmpPrimVtx, electronsFromConversions ));
00384 
00385        } else {
00386          isAconversion=0;
00387          //std::cout  << " UNCONVERTED photon " <<   "\n";    
00388          CLHEP::HepLorentzVector vtxPosition(0.,0.,0.,0.);
00389          CLHEP::HepLorentzVector tmpPhoMom( (*iPhoTk).momentum().px(), (*iPhoTk).momentum().py(),
00390                                      (*iPhoTk).momentum().pz(), (*iPhoTk).momentum().e() ) ;
00391          CLHEP::HepLorentzVector tmpPrimVtx( primVtxPos.x(), primVtxPos.y(), primVtxPos.z(), primVtxPos.t() ) ;
00392          //      result.push_back( PhotonMCTruth(isAconversion, (*iPhoTk).momentum(),  photonVertexIndex, phoTrkId, vtxPosition,   primVtx.position(), trkFromConversion ));
00393          result.push_back( PhotonMCTruth(isAconversion, tmpPhoMom,  photonVertexIndex, phoTrkId, phoMotherType, phoMotherMom, phoMotherVtx, vtxPosition,   
00394          tmpPrimVtx, electronsFromConversions ));
00395         
00396        }
00397        
00398 
00399 
00400      } // End loop over the primary photons
00401      
00402      
00403    }   // Event with one or two photons 
00404 
00405 
00406    //std::cout << "  PhotonMCTruthFinder photon size " << result.size() << std::endl;
00407   }
00408 
00409    return result;
00410 
00411 }
00412 
00413 void PhotonMCTruthFinder::fill(std::vector<SimTrack>& simTracks, std::vector<SimVertex>& simVertices ) {
00414   //  std::cout << "  PhotonMCTruthFinder::fill " << std::endl;
00415 
00416 
00417 
00418  // Watch out there ! A SimVertex is in mm (stupid), 
00419  
00420   unsigned nVtx = simVertices.size();
00421   unsigned nTks = simTracks.size();
00422 
00423   //  std::cout << "  PhotonMCTruthFinder::fill " << nVtx << " " << nTks << std::endl;
00424 
00425   // Empty event, do nothin'
00426   if ( nVtx == 0 ) return;
00427 
00428   // create a map associating geant particle id and position in the 
00429   // event SimTrack vector
00430   for( unsigned it=0; it<nTks; ++it ) {
00431     geantToIndex_[ simTracks[it].trackId() ] = it;
00432     //    std::cout << " PhotonMCTruthFinder::fill it " << it << " simTracks[it].trackId() " <<  simTracks[it].trackId() << std::endl;
00433  
00434   }  
00435 
00436 
00437 
00438 
00439 }
/data/refman/pasoursint/CMSSW_4_1_8_patch13/src/RecoEgamma/EgammaMCTools/src/PhotonMCTruthFinder.cc
