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.empty() ) {
  
  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.empty() ) {
       // 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;
 
  }  




}