|
|
|
#include <ElectronMCTruthFinder.h>
Public Member Functions | |
| ElectronMCTruthFinder () | |
| std::vector< ElectronMCTruth > | find (std::vector< SimTrack > simTracks, std::vector< SimVertex > simVertices) |
| virtual | ~ElectronMCTruthFinder () |
Private Member Functions | |
| void | fill (std::vector< SimTrack > &theSimTracks, std::vector< SimVertex > &theSimVertices) |
Private Attributes | |
| std::map< unsigned, unsigned > | geantToIndex_ |
| PhotonMCTruthFinder * | thePhotonMCTruthFinder_ |
Definition at line 13 of file ElectronMCTruthFinder.h.
| ElectronMCTruthFinder::ElectronMCTruthFinder | ( | ) |
Definition at line 6 of file ElectronMCTruthFinder.cc.
{
}
| ElectronMCTruthFinder::~ElectronMCTruthFinder | ( | ) | [virtual] |
Definition at line 11 of file ElectronMCTruthFinder.cc.
{
}
| void ElectronMCTruthFinder::fill | ( | std::vector< SimTrack > & | theSimTracks, |
| std::vector< SimVertex > & | theSimVertices | ||
| ) | [private] |
Definition at line 232 of file ElectronMCTruthFinder.cc.
References geantToIndex_.
Referenced by find().
{
//std::cout << " ElectronMCTruthFinder::fill " << std::endl;
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 << " ElectronMCTruthFinder::fill it " << it << " simTracks[it].trackId() " << simTracks[it].trackId() << std::endl;
}
}
| std::vector< ElectronMCTruth > ElectronMCTruthFinder::find | ( | std::vector< SimTrack > | simTracks, |
| std::vector< SimVertex > | simVertices | ||
| ) |
Now store the electron truth
here fill the electron
Definition at line 17 of file ElectronMCTruthFinder.cc.
References abs, fill(), geantToIndex_, if(), CoreSimTrack::momentum(), SimVertex::parentIndex(), CoreSimVertex::position(), query::result, CoreSimTrack::trackId(), CoreSimTrack::type(), and SimTrack::vertIndex().
Referenced by MCElectronAnalyzer::analyze().
{
//std::cout << " ElectronMCTruthFinder::find " << std::endl;
std::vector<ElectronMCTruth> result;
// 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> electronTracks;
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;
}
}
//std::cout << " Loop over all particles " << 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() << " vertex ID " << vertexId << std::endl;
if ( (*iSimTk).vertIndex() == iPV ) {
npv++;
if ( std::abs((*iSimTk).type() ) == 11) {
//std::cout << " Found a primary electron with ID " << (*iSimTk).trackId() << " momentum " << (*iSimTk).momentum() << std::endl;
electronTracks.push_back( *iSimTk );
}
}
}
//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;
}
}
std::vector<CLHEP::Hep3Vector> bremPos;
std::vector<CLHEP::HepLorentzVector> pBrem;
std::vector<float> xBrem;
for (std::vector<SimTrack>::iterator iEleTk = electronTracks.begin(); iEleTk != electronTracks.end(); ++iEleTk){
//std::cout << " Looping on the primary electron pt " << std::sqrt((*iEleTk).momentum().perp2()) << " electron track ID " << (*iEleTk).trackId() << std::endl;
SimTrack trLast =(*iEleTk);
unsigned int eleId = (*iEleTk).trackId();
float remainingEnergy =trLast.momentum().e();
// CLHEP::HepLorentzVector motherMomentum = (*iEleTk).momentum();
// CLHEP::HepLorentzVector primEleMom = (*iEleTk).momentum();
math::XYZTLorentzVectorD motherMomentum((*iEleTk).momentum().x(),
(*iEleTk).momentum().y(),
(*iEleTk).momentum().z(),
(*iEleTk).momentum().e());
math::XYZTLorentzVectorD primEleMom(motherMomentum);
int eleVtxIndex= (*iEleTk).vertIndex();
bremPos.clear();
pBrem.clear();
xBrem.clear();
for (std::vector<SimTrack>::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 " << std::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 " << std::sqrt(vertex2.position().perp2()) << " parent index " << vertex2.parentIndex() << std::endl;
//std::cout << " Electron pt " << std::sqrt((*iSimTk).momentum().perp2()) << " photon pt " << std::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() ) {
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: 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());
result.push_back ( ElectronMCTruth( tmpEleMom, eleVtxIndex, bremPos, pBrem, xBrem, tmpVtxPos,(*iEleTk) ) ) ;
} // End loop over primary electrons
return result;
}
std::map<unsigned, unsigned> ElectronMCTruthFinder::geantToIndex_ [private] |
Definition at line 28 of file ElectronMCTruthFinder.h.
Definition at line 29 of file ElectronMCTruthFinder.h.
1.7.3