PizeroMCTruthFinder.cc

Go to the documentation of this file.

#include "RecoEgamma/EgammaMCTools/interface/PizeroMCTruthFinder.h"
#include "RecoEgamma/EgammaMCTools/interface/PhotonMCTruthFinder.h"
#include "RecoEgamma/EgammaMCTools/interface/ElectronMCTruthFinder.h"


#include <algorithm>                                                          

PizeroMCTruthFinder::PizeroMCTruthFinder() {

  thePhotonMCTruthFinder_ = new PhotonMCTruthFinder();
  theElectronMCTruthFinder_ = new ElectronMCTruthFinder();


}

PizeroMCTruthFinder::~PizeroMCTruthFinder() 
{

  delete thePhotonMCTruthFinder_;
  delete theElectronMCTruthFinder_;
  std::cout << "~PizeroMCTruthFinder" << std::endl;
}

std::vector<PizeroMCTruth> PizeroMCTruthFinder::find(std::vector<SimTrack> theSimTracks, std::vector<SimVertex> theSimVertices ) {
  std::cout << "  PizeroMCTruthFinder::find " << std::endl;

  
  std::vector<PhotonMCTruth> mcPhotons=thePhotonMCTruthFinder_->find (theSimTracks,  theSimVertices);  

  std::vector<PizeroMCTruth> result;
  std::vector<PhotonMCTruth> photonsFromPizero;
  std::vector<ElectronMCTruth> electronsFromPizero;

  // 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> pizeroTracks;
  SimVertex primVtx;   
  
    
  fill(theSimTracks,  theSimVertices);
  
  int iPV=-1;   
  int partType1=0;
  int partType2=0;
  std::vector<SimTrack>::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;
  }
  
  // CLHEP::HepLorentzVector primVtxPos= primVtx.position(); 
  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;
    }
  }
  
  int npv=0;
  
  for (std::vector<SimTrack>::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() << std::endl;  

    if ( (*iSimTk).vertIndex() == iPV ) {
      npv++;
      if ( std::abs((*iSimTk).type() ) == 111) {

    std::cout << " Found a primary pizero with ID  " << (*iSimTk).trackId() << " momentum " << (*iSimTk).momentum() <<  std::endl;
    
    pizeroTracks.push_back( *iSimTk );

    // CLHEP::HepLorentzVector momentum = (*iSimTk).momentum();
    math::XYZTLorentzVectorD momentum((*iSimTk).momentum().x(),
                                      (*iSimTk).momentum().y(),
                                      (*iSimTk).momentum().z(),
                                      (*iSimTk).momentum().e());


      } 
    }
  }
 
  
  std::cout << " There are " << npv << " particles originating in the PV " << std::endl;
  
  if(npv > 4) {
    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;
    }
  }


  
  for (std::vector<SimTrack>::iterator iPizTk = pizeroTracks.begin(); iPizTk != pizeroTracks.end(); ++iPizTk){
    std::cout << " Looping on the primary pizero pt  " << sqrt((*iPizTk).momentum().perp2()) << " pizero track ID " << (*iPizTk).trackId() << std::endl;
    
    photonsFromPizero.clear();
    std::cout << " mcPhotons.size " << mcPhotons.size() << std::endl;
    for ( std::vector<PhotonMCTruth>::iterator iPho=mcPhotons.begin(); iPho !=mcPhotons.end(); ++iPho ){
      int phoVtxIndex = (*iPho).vertexInd();
      SimVertex phoVtx = theSimVertices[phoVtxIndex];
      unsigned int phoParentInd= phoVtx.parentIndex();
      std::cout << " photon parent vertex index " << phoParentInd << std::endl;
      
      if (phoParentInd == (*iPizTk).trackId())  {
    std::cout << "Matched Photon ID " << (*iPho).trackId() << "  vtx " << phoParentInd << " with pizero " << (*iPizTk).trackId() << std::endl;
    photonsFromPizero.push_back( *iPho);
    
      }
    }
    std::cout << " Photon matching the pizero vertex " << photonsFromPizero.size() <<std::endl;
    
    
    // build pizero MC thruth
    CLHEP::HepLorentzVector tmpMom( (*iPizTk).momentum().px(), (*iPizTk).momentum().py(),
                             (*iPizTk).momentum().pz(), (*iPizTk).momentum().e() ) ;
    CLHEP::HepLorentzVector tmpPos( primVtx.position().x(), primVtx.position().y(),
                             primVtx.position().z(), primVtx.position().t() ) ;
    result.push_back( PizeroMCTruth (  tmpMom, photonsFromPizero, tmpPos ) );
    
    
  }   // end loop over primary pizeros

    
  std::cout << " Pizero size " << result.size() <<  std::endl;
  
  
  return result;
}



void PizeroMCTruthFinder::fill(std::vector<SimTrack>& simTracks, 
  std::vector<SimVertex>& simVertices ) {


unsigned nVtx = simVertices.size();
unsigned nTks = simTracks.size();

// 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 << " PizeroMCTruthFinder::fill it " << it << " simTracks[it].trackId() " <<  simTracks[it].trackId() << std::endl;
 
}  


}