/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_7/src/RecoEgamma/EgammaMCTools/src/PhotonMCTruthFinder.cc

Go to the documentation of this file.

#include "RecoEgamma/EgammaMCTools/interface/PhotonMCTruthFinder.h"
#include "RecoEgamma/EgammaMCTools/interface/PhotonMCTruth.h"
#include "RecoEgamma/EgammaMCTools/interface/ElectronMCTruth.h"
//
//
#include "SimDataFormats/CrossingFrame/interface/CrossingFrame.h"
#include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
#include "SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h"

#include <algorithm>


PhotonMCTruthFinder::PhotonMCTruthFinder( ) {
 //std::cout << " PhotonMCTruthFinder CTOR " << std::endl;

 
}

std::vector<PhotonMCTruth> PhotonMCTruthFinder::find(const std::vector<SimTrack>& theSimTracks, const std::vector<SimVertex>& theSimVertices ) {
  //  std::cout << "  PhotonMCTruthFinder::find " << std::endl;

  std::vector<PhotonMCTruth> result;

  //const float pi = 3.141592653592;
  //const float twopi=2*pi;

// Local variables  
  const int SINGLE=1;
  const int DOUBLE=2;
  const int PYTHIA=3;
  const int ELECTRON_FLAV=1;
  const int PIZERO_FLAV=2;
  const int PHOTON_FLAV=3;
  
  int ievtype=0;
  int ievflav=0;
  
  
  std::vector<SimTrack*> photonTracks;
  
  std::vector<SimTrack> trkFromConversion;
  std::vector<ElectronMCTruth> electronsFromConversions;
  SimVertex primVtx;   

  
  fill(theSimTracks,  theSimVertices);

  //  std::cout << " After fill " << theSimTracks.size() << " " << theSimVertices.size() << std::endl;
  if (    theSimTracks.size() != 0 ) {
  
  int iPV=-1;   
  int partType1=0;
  int partType2=0;
  std::vector<SimTrack>::const_iterator iFirstSimTk = theSimTracks.begin();
  if (  !(*iFirstSimTk).noVertex() ) {
    iPV =  (*iFirstSimTk).vertIndex();
    
    int vtxId =   (*iFirstSimTk).vertIndex();
    primVtx = theSimVertices[vtxId];
    
    
    
    partType1 = (*iFirstSimTk).type();
    //    std::cout <<  " Primary vertex id " << iPV << " first track type " << (*iFirstSimTk).type() << std::endl;  
  } else {
    //std::cout << " First track has no vertex " << std::endl;
    
  }
  


  math::XYZTLorentzVectorD primVtxPos(primVtx.position().x(),
                                       primVtx.position().y(),
                                       primVtx.position().z(),
                                       primVtx.position().e());           
 

  // Look at a second track
  iFirstSimTk++;
  if ( iFirstSimTk!=  theSimTracks.end() ) {
    
    if (  (*iFirstSimTk).vertIndex() == iPV) {
      partType2 = (*iFirstSimTk).type();  
      //      std::cout <<  " second track type " << (*iFirstSimTk).type() << " vertex " <<  (*iFirstSimTk).vertIndex() << std::endl;  
      
    } else {
      // std::cout << " Only one kine track from Primary Vertex " << std::endl;
    }
  }
  
  //std::cout << " Loop over all particles " << std::endl;
  
  int npv=0;
  //int iPho=0;
  //int iPizero=0;
  //   theSimTracks.reset();
  for (std::vector<SimTrack>::const_iterator iSimTk = theSimTracks.begin(); iSimTk != theSimTracks.end(); ++iSimTk){
    if (  (*iSimTk).noVertex() ) continue;


    //int vertexId = (*iSimTk).vertIndex();
    //SimVertex vertex = theSimVertices[vertexId];
 
    //    std::cout << " Particle type " <<  (*iSimTk).type() << " Sim Track ID " << (*iSimTk).trackId() << " momentum " << (*iSimTk).momentum() <<  " vertex position " << vertex.position() << " vertex index " << (*iSimTk).vertIndex() << std::endl;  
    if ( (*iSimTk).vertIndex() == iPV ) {
      npv++;
      if ( (*iSimTk).type() == 22) {
        //      std::cout << " Found a primary photon with ID  " << (*iSimTk).trackId() << " momentum " << (*iSimTk).momentum() <<  std::endl; 

        photonTracks.push_back(&(const_cast<SimTrack&>(*iSimTk)));

      } 
      
    }
    
  } 
  
  //  std::cout << " There are " << npv << " particles originating in the PV " << std::endl;
  
   if(npv >= 3) {
     ievtype = PYTHIA;
   } else if(npv == 1) {
     if( std::abs(partType1) == 11 ) {
       ievtype = SINGLE; ievflav = ELECTRON_FLAV;
     } else if(partType1 == 111) {
       ievtype = SINGLE; ievflav = PIZERO_FLAV;
     } else if(partType1 == 22) {
       ievtype = SINGLE; ievflav = PHOTON_FLAV;
     }
   } else if(npv == 2) {
     if (  std::abs(partType1) == 11 && std::abs(partType2) == 11 ) {
       ievtype = DOUBLE; ievflav = ELECTRON_FLAV;
     } else if(partType1 == 111 && partType2 == 111)   {
       ievtype = DOUBLE; ievflav = PIZERO_FLAV;
     } else if(partType1 == 22 && partType2 == 22)   {
       ievtype = DOUBLE; ievflav = PHOTON_FLAV;
     }
   }      
   

   int isAconversion=0;   
   int phoMotherType=-1;
   int phoMotherVtxIndex=-1;
   int phoMotherId=-1;
   if( ievflav == PHOTON_FLAV || ievflav== PIZERO_FLAV || ievtype == PYTHIA ) {

     //     std::cout << " It's a primary PHOTON or PIZERO or PYTHIA event with " << photonTracks.size() << " photons ievtype " << ievtype << " ievflav " << ievflav<<  std::endl;

     //    for (std::vector<SimTrack*>::iterator iPhoTk = photonTracks.begin(); iPhoTk != photonTracks.end(); ++iPhoTk){
     //      std::cout << " All gamma found from PV " << (*iPhoTk)->momentum() << " photon track ID " << (*iPhoTk)->trackId() << " vertex index " << (*iPhoTk)->vertIndex() << std::endl;  
     //  }        

     for (std::vector<SimTrack>::const_iterator iPhoTk = theSimTracks.begin(); iPhoTk != theSimTracks.end(); ++iPhoTk){

       trkFromConversion.clear();           
       electronsFromConversions.clear();

       if ( (*iPhoTk).type() != 22 ) continue;
       int photonVertexIndex= (*iPhoTk).vertIndex();
       int phoTrkId= (*iPhoTk).trackId();
       //std::cout << " Looping on gamma looking for conversions " << (*iPhoTk).momentum() << " photon track ID " << (*iPhoTk).trackId() << std::endl;
     
       // check who is his mother
       SimVertex vertex = theSimVertices[ photonVertexIndex];
       phoMotherId=-1;
       if ( vertex.parentIndex() != -1 ) {
         
         unsigned  motherGeantId = vertex.parentIndex(); 
         std::map<unsigned, unsigned >::iterator association = geantToIndex_.find( motherGeantId );
         if(association != geantToIndex_.end() )
           phoMotherId = association->second;
         phoMotherType = phoMotherId == -1 ? 0 : theSimTracks[phoMotherId].type();

         

         if ( phoMotherType == 111 || phoMotherType == 221 || phoMotherType == 331 ) {
         
           //std::cout << " Parent to this vertex   motherId " << phoMotherId << " mother type " <<  phoMotherType << " Sim track ID " <<  theSimTracks[phoMotherId].trackId() << std::endl;
           //std::cout << " Son of a pizero or eta " << phoMotherType << std::endl;
         }
         
       }



       for (std::vector<SimTrack>::const_iterator iEleTk = theSimTracks.begin(); iEleTk != theSimTracks.end(); ++iEleTk){
         if (  (*iEleTk).noVertex() )                    continue;
         if ( (*iEleTk).vertIndex() == iPV )             continue; 
         if ( std::abs((*iEleTk).type()) != 11  )             continue;

         int vertexId = (*iEleTk).vertIndex();
         SimVertex vertex = theSimVertices[vertexId];
         int motherId=-1;
        

         //std::cout << " Secondary from photons particle type " << (*iEleTk).type() << " trackId " <<  (*iEleTk).trackId() << " vertex ID " << vertexId << std::endl;
         if ( vertex.parentIndex() != -1 ) {

           unsigned  motherGeantId = vertex.parentIndex(); 
           std::map<unsigned, unsigned >::iterator association = geantToIndex_.find( motherGeantId );
           if(association != geantToIndex_.end() )
             motherId = association->second;
           
           
           //int motherType = motherId == -1 ? 0 : theSimTracks[motherId].type();
           
           //std::cout << " Parent to this vertex   motherId " << motherId << " mother type " <<  motherType << " Sim track ID " <<  theSimTracks[motherId].trackId() << std::endl;

           std::vector<CLHEP::Hep3Vector> bremPos;  
           std::vector<CLHEP::HepLorentzVector> pBrem;
           std::vector<float> xBrem;
           
           if ( theSimTracks[motherId].trackId() == (*iPhoTk).trackId() ) {
             //std::cout << " Found the Mother Photon " << std::endl;



             trkFromConversion.push_back( (*iEleTk ) );
             SimTrack trLast =(*iEleTk); 
             float remainingEnergy =trLast.momentum().e();
             //HepLorentzVector primEleMom=(*iEleTk).momentum();  
             //HepLorentzVector motherMomentum=(*iEleTk).momentum();  
             math::XYZTLorentzVectorD primEleMom((*iEleTk).momentum().x(),
                                             (*iEleTk).momentum().y(),
                                             (*iEleTk).momentum().z(),
                                             (*iEleTk).momentum().e());  
             math::XYZTLorentzVectorD motherMomentum(primEleMom);  
             unsigned int eleId = (*iEleTk).trackId();     
             int eleVtxIndex= (*iEleTk).vertIndex();
           
    
             bremPos.clear();
             pBrem.clear();
             xBrem.clear();   
                     
             for (std::vector<SimTrack>::const_iterator iSimTk = theSimTracks.begin(); iSimTk != theSimTracks.end(); ++iSimTk){
               
               if (  (*iSimTk).noVertex() )                    continue;
               if ( (*iSimTk).vertIndex() == iPV )             continue;
               
               //std::cout << " (*iEleTk)->trackId() " << (*iEleTk).trackId() << " (*iEleTk)->vertIndex() "<< (*iEleTk).vertIndex()  << " (*iSimTk).vertIndex() "  <<  (*iSimTk).vertIndex() << " (*iSimTk).type() " <<   (*iSimTk).type() << " (*iSimTk).trackId() " << (*iSimTk).trackId() << std::endl;
               
               int vertexId1 = (*iSimTk).vertIndex();
               SimVertex vertex1 = theSimVertices[vertexId1];
               int vertexId2 = trLast.vertIndex();
               //SimVertex vertex2 = theSimVertices[vertexId2];
               
               
               int motherId=-1;
               
               if(  (  vertexId1 ==  vertexId2 ) && ( (*iSimTk).type() == (*iEleTk).type() ) && trLast.type() == 22   ) {
                 //std::cout << " Here a e/gamma brem vertex " << std::endl;
                 
                 //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;
                 
                 //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;
                 
                 //std::cout << " Electron pt " << sqrt((*iSimTk).momentum().perp2()) << " photon pt " << sqrt(trLast.momentum().perp2()) << " Mother electron pt " <<  sqrt(motherMomentum.perp2()) << std::endl;
                 //std::cout << " eleId " << eleId << std::endl;
                 float eLoss = remainingEnergy - ( (*iSimTk).momentum() + trLast.momentum()).e();
                 //std::cout << " eLoss " << eLoss << std::endl;              
                 
                 
                 if ( vertex1.parentIndex() != -1  ) {
                   
                   unsigned  motherGeantId = vertex1.parentIndex(); 
                   std::map<unsigned, unsigned >::iterator association = geantToIndex_.find( motherGeantId );
                   if(association != geantToIndex_.end() )
                     motherId = association->second;
                   
                   //int motherType = motherId == -1 ? 0 : theSimTracks[motherId].type();
                   //std::cout << " Parent to this vertex   motherId " << motherId << " mother type " <<  motherType << " Sim track ID " <<  theSimTracks[motherId].trackId() << std::endl; 
                   if (theSimTracks[motherId].trackId() == eleId ) {
                     
                     //std::cout << "  ***** Found the Initial Mother Electron ****   theSimTracks[motherId].trackId() " <<  theSimTracks[motherId].trackId() << " eleId " <<  eleId << std::endl;
                     eleId= (*iSimTk).trackId();
                     remainingEnergy =   (*iSimTk).momentum().e();
                     motherMomentum = (*iSimTk).momentum();
                     
                     
                     pBrem.push_back( CLHEP::HepLorentzVector(trLast.momentum().px(),
                                                       trLast.momentum().py(),
                                                       trLast.momentum().pz(),
                                                       trLast.momentum().e()) );
                     bremPos.push_back( CLHEP::HepLorentzVector(vertex1.position().x(),
                                                         vertex1.position().y(),
                                                         vertex1.position().z(),
                                                         vertex1.position().t()) );
                     xBrem.push_back(eLoss);
                     
                   }
                   
                   
                   
                   
                 } else {
                   //std::cout << " This vertex has no parent tracks " <<  std::endl;
                 }
                 
               }
               trLast=(*iSimTk);
               
             } // End loop over all SimTracks 
             //std::cout << " Going to build the ElectronMCTruth for this electron from converted photon: pBrem size " << pBrem.size() << std::endl;

             CLHEP::HepLorentzVector tmpEleMom(primEleMom.px(),primEleMom.py(),
                                        primEleMom.pz(),primEleMom.e() );
             CLHEP::HepLorentzVector tmpVtxPos(primVtxPos.x(),primVtxPos.y(),
                                        primVtxPos.z(),primVtxPos.t() );
             electronsFromConversions.push_back ( 
                ElectronMCTruth( tmpEleMom, eleVtxIndex,  bremPos, pBrem, 
                                 xBrem,  tmpVtxPos , const_cast<SimTrack&>(*iEleTk)  )  ) ;
           }   
 
           


         } else {
           //std::cout << " This vertex has no parent tracks " <<  std::endl;
         }
         
         
       } // End of loop over the SimTracks      
       
       //std::cout << " DEBUG trkFromConversion.size() " << trkFromConversion.size() << " electronsFromConversions.size() " << electronsFromConversions.size() << std::endl;

       math::XYZTLorentzVectorD motherVtxPosition(0.,0.,0.,0.);
       CLHEP::HepLorentzVector phoMotherMom(0.,0.,0.,0.);
       CLHEP::HepLorentzVector phoMotherVtx(0.,0.,0.,0.); 

       if ( phoMotherId >= 0) {
        
         phoMotherVtxIndex = theSimTracks[phoMotherId].vertIndex();
         SimVertex motherVtx = theSimVertices[ phoMotherVtxIndex];
         motherVtxPosition =math::XYZTLorentzVectorD (motherVtx.position().x(),
                                                      motherVtx.position().y(),
                                                      motherVtx.position().z(),
                                                      motherVtx.position().e());
         
         phoMotherMom.setPx( theSimTracks[phoMotherId].momentum().x());
         phoMotherMom.setPy( theSimTracks[phoMotherId].momentum().y());
         phoMotherMom.setPz( theSimTracks[phoMotherId].momentum().z() );
         phoMotherMom.setE( theSimTracks[phoMotherId].momentum().t());
         // std::cout << " PhotonMCTruthFinder mother " << phoMotherId << " type " << phoMotherType << " Momentum" <<  phoMotherMom.et() << std::endl;  
 
         phoMotherVtx.setX ( motherVtxPosition.x());
         phoMotherVtx.setY ( motherVtxPosition.y());
         phoMotherVtx.setZ ( motherVtxPosition.z());
         phoMotherVtx.setT ( motherVtxPosition.t());
       
       }


       if ( electronsFromConversions.size() > 0 ) {
       // if ( trkFromConversion.size() > 0 ) {
         isAconversion=1;
         //std::cout  << " CONVERTED photon " <<   "\n";    

         //int convVtxId =  trkFromConversion[0].vertIndex();
         int convVtxId =electronsFromConversions[0].vertexInd();
         SimVertex convVtx = theSimVertices[convVtxId];
         // CLHEP::HepLorentzVector vtxPosition = convVtx.position();
         math::XYZTLorentzVectorD vtxPosition(convVtx.position().x(),
                                          convVtx.position().y(),
                                          convVtx.position().z(),
                                          convVtx.position().e());
         
           
         //result.push_back( PhotonMCTruth(isAconversion, (*iPhoTk).momentum(), photonVertexIndex, phoTrkId, vtxPosition,   primVtx.position(), trkFromConversion ));
         CLHEP::HepLorentzVector tmpPhoMom( (*iPhoTk).momentum().px(), (*iPhoTk).momentum().py(),
                                     (*iPhoTk).momentum().pz(), (*iPhoTk).momentum().e() ) ;
         CLHEP::HepLorentzVector tmpVertex( vtxPosition.x(), vtxPosition.y(), vtxPosition.z(), vtxPosition.t() );
         CLHEP::HepLorentzVector tmpPrimVtx( primVtxPos.x(), primVtxPos.y(), primVtxPos.z(), primVtxPos.t() ) ;



         result.push_back( PhotonMCTruth(isAconversion, tmpPhoMom, photonVertexIndex, phoTrkId, phoMotherType,phoMotherMom, phoMotherVtx, tmpVertex,  
         tmpPrimVtx, electronsFromConversions ));

       } else {
         isAconversion=0;
         //std::cout  << " UNCONVERTED photon " <<   "\n";    
         CLHEP::HepLorentzVector vtxPosition(0.,0.,0.,0.);
         CLHEP::HepLorentzVector tmpPhoMom( (*iPhoTk).momentum().px(), (*iPhoTk).momentum().py(),
                                     (*iPhoTk).momentum().pz(), (*iPhoTk).momentum().e() ) ;
         CLHEP::HepLorentzVector tmpPrimVtx( primVtxPos.x(), primVtxPos.y(), primVtxPos.z(), primVtxPos.t() ) ;
         //      result.push_back( PhotonMCTruth(isAconversion, (*iPhoTk).momentum(),  photonVertexIndex, phoTrkId, vtxPosition,   primVtx.position(), trkFromConversion ));
         result.push_back( PhotonMCTruth(isAconversion, tmpPhoMom,  photonVertexIndex, phoTrkId, phoMotherType, phoMotherMom, phoMotherVtx, vtxPosition,   
         tmpPrimVtx, electronsFromConversions ));
        
       }
       


     } // End loop over the primary photons
     
     
   }   // Event with one or two photons 


   //std::cout << "  PhotonMCTruthFinder photon size " << result.size() << std::endl;
  }

   return result;

}

void PhotonMCTruthFinder::fill(const std::vector<SimTrack>& simTracks, const std::vector<SimVertex>& simVertices ) {
  //  std::cout << "  PhotonMCTruthFinder::fill " << std::endl;



 // Watch out there ! A SimVertex is in mm (stupid), 
 
  unsigned nVtx = simVertices.size();
  unsigned nTks = simTracks.size();

  //  std::cout << "  PhotonMCTruthFinder::fill " << nVtx << " " << nTks << std::endl;

  // Empty event, do nothin'
  if ( nVtx == 0 ) return;

  // create a map associating geant particle id and position in the 
  // event SimTrack vector
  for( unsigned it=0; it<nTks; ++it ) {
    geantToIndex_[ simTracks[it].trackId() ] = it;
    //    std::cout << " PhotonMCTruthFinder::fill it " << it << " simTracks[it].trackId() " <<  simTracks[it].trackId() << std::endl;
 
  }  




}