#include <RecoEgamma/EgammaMCTools/interface/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 | ( | ) |
ElectronMCTruthFinder::~ElectronMCTruthFinder | ( | ) | [virtual] |
void ElectronMCTruthFinder::fill | ( | std::vector< SimTrack > & | theSimTracks, | |
std::vector< SimVertex > & | theSimVertices | |||
) | [private] |
Definition at line 232 of file ElectronMCTruthFinder.cc.
References geantToIndex_, and it.
Referenced by find().
00233 { 00234 //std::cout << " ElectronMCTruthFinder::fill " << std::endl; 00235 00236 unsigned nVtx = simVertices.size(); 00237 unsigned nTks = simTracks.size(); 00238 00239 // Empty event, do nothin' 00240 if ( nVtx == 0 ) return; 00241 00242 // create a map associating geant particle id and position in the 00243 // event SimTrack vector 00244 for( unsigned it=0; it<nTks; ++it ) { 00245 geantToIndex_[ simTracks[it].trackId() ] = it; 00246 //std::cout << " ElectronMCTruthFinder::fill it " << it << " simTracks[it].trackId() " << simTracks[it].trackId() << std::endl; 00247 00248 } 00249 00250 00251 }
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 funct::abs(), e, fill(), geantToIndex_, CoreSimTrack::momentum(), SimVertex::parentIndex(), CoreSimVertex::position(), HLT_VtxMuL3::result, CoreSimTrack::type(), and SimTrack::vertIndex().
Referenced by MCElectronAnalyzer::analyze().
00017 { 00018 //std::cout << " ElectronMCTruthFinder::find " << std::endl; 00019 00020 00021 00022 00023 std::vector<ElectronMCTruth> result; 00024 00025 // Local variables 00026 const int SINGLE=1; 00027 const int DOUBLE=2; 00028 const int PYTHIA=3; 00029 const int ELECTRON_FLAV=1; 00030 const int PIZERO_FLAV=2; 00031 const int PHOTON_FLAV=3; 00032 00033 int ievtype=0; 00034 int ievflav=0; 00035 00036 00037 std::vector<SimTrack> electronTracks; 00038 SimVertex primVtx; 00039 00040 fill(theSimTracks, theSimVertices); 00041 00042 int iPV=-1; 00043 int partType1=0; 00044 int partType2=0; 00045 std::vector<SimTrack>::iterator iFirstSimTk = theSimTracks.begin(); 00046 if ( !(*iFirstSimTk).noVertex() ) { 00047 iPV = (*iFirstSimTk).vertIndex(); 00048 00049 int vtxId = (*iFirstSimTk).vertIndex(); 00050 primVtx = theSimVertices[vtxId]; 00051 00052 partType1 = (*iFirstSimTk).type(); 00053 00054 00055 //std::cout << " Primary vertex id " << iPV << " first track type " << (*iFirstSimTk).type() << std::endl; 00056 } else { 00057 //std::cout << " First track has no vertex " << std::endl; 00058 } 00059 00060 // HepLorentzVector primVtxPos= primVtx.position(); 00061 math::XYZTLorentzVectorD primVtxPos(primVtx.position().x(), 00062 primVtx.position().y(), 00063 primVtx.position().z(), 00064 primVtx.position().e()); 00065 00066 // Look at a second track 00067 iFirstSimTk++; 00068 if ( iFirstSimTk!= theSimTracks.end() ) { 00069 00070 if ( (*iFirstSimTk).vertIndex() == iPV) { 00071 partType2 = (*iFirstSimTk).type(); 00072 00073 //std::cout << " second track type " << (*iFirstSimTk).type() << " vertex " << (*iFirstSimTk).vertIndex() << std::endl; 00074 00075 } else { 00076 //std::cout << " Only one kine track from Primary Vertex " << std::endl; 00077 } 00078 } 00079 00080 //std::cout << " Loop over all particles " << std::endl; 00081 00082 int npv=0; 00083 for (std::vector<SimTrack>::iterator iSimTk = theSimTracks.begin(); iSimTk != theSimTracks.end(); ++iSimTk){ 00084 if ( (*iSimTk).noVertex() ) continue; 00085 00086 int vertexId = (*iSimTk).vertIndex(); 00087 SimVertex vertex = theSimVertices[vertexId]; 00088 00089 //std::cout << " Particle type " << (*iSimTk).type() << " Sim Track ID " << (*iSimTk).trackId() << " momentum " << (*iSimTk).momentum() << " vertex position " << vertex.position() << " vertex ID " << vertexId << std::endl; 00090 if ( (*iSimTk).vertIndex() == iPV ) { 00091 npv++; 00092 if ( fabs((*iSimTk).type() ) == 11) { 00093 //std::cout << " Found a primary electron with ID " << (*iSimTk).trackId() << " momentum " << (*iSimTk).momentum() << std::endl; 00094 electronTracks.push_back( *iSimTk ); 00095 } 00096 } 00097 } 00098 //std::cout << " There are " << npv << " particles originating in the PV " << std::endl; 00099 00100 00101 if(npv > 4) { 00102 ievtype = PYTHIA; 00103 } else if(npv == 1) { 00104 if( abs(partType1) == 11 ) { 00105 ievtype = SINGLE; ievflav = ELECTRON_FLAV; 00106 } else if(partType1 == 111) { 00107 ievtype = SINGLE; ievflav = PIZERO_FLAV; 00108 } else if(partType1 == 22) { 00109 ievtype = SINGLE; ievflav = PHOTON_FLAV; 00110 } 00111 } else if(npv == 2) { 00112 if ( abs(partType1) == 11 && abs(partType2) == 11 ) { 00113 ievtype = DOUBLE; ievflav = ELECTRON_FLAV; 00114 } else if(partType1 == 111 && partType2 == 111) { 00115 ievtype = DOUBLE; ievflav = PIZERO_FLAV; 00116 } else if(partType1 == 22 && partType2 == 22) { 00117 ievtype = DOUBLE; ievflav = PHOTON_FLAV; 00118 } 00119 } 00120 00121 00122 00124 std::vector<Hep3Vector> bremPos; 00125 std::vector<HepLorentzVector> pBrem; 00126 std::vector<float> xBrem; 00127 00128 00129 for (std::vector<SimTrack>::iterator iEleTk = electronTracks.begin(); iEleTk != electronTracks.end(); ++iEleTk){ 00130 //std::cout << " Looping on the primary electron pt " << std::sqrt((*iEleTk).momentum().perp2()) << " electron track ID " << (*iEleTk).trackId() << std::endl; 00131 00132 00133 00134 SimTrack trLast =(*iEleTk); 00135 int eleId = (*iEleTk).trackId(); 00136 float remainingEnergy =trLast.momentum().e(); 00137 // HepLorentzVector motherMomentum = (*iEleTk).momentum(); 00138 // HepLorentzVector primEleMom = (*iEleTk).momentum(); 00139 math::XYZTLorentzVectorD motherMomentum((*iEleTk).momentum().x(), 00140 (*iEleTk).momentum().y(), 00141 (*iEleTk).momentum().z(), 00142 (*iEleTk).momentum().e()); 00143 math::XYZTLorentzVectorD primEleMom(motherMomentum); 00144 int eleVtxIndex= (*iEleTk).vertIndex(); 00145 00146 bremPos.clear(); 00147 pBrem.clear(); 00148 xBrem.clear(); 00149 00150 00151 for (std::vector<SimTrack>::iterator iSimTk = theSimTracks.begin(); iSimTk != theSimTracks.end(); ++iSimTk){ 00152 00153 if ( (*iSimTk).noVertex() ) continue; 00154 if ( (*iSimTk).vertIndex() == iPV ) continue; 00155 00156 //std::cout << " (*iEleTk)->trackId() " << (*iEleTk).trackId() << " (*iEleTk)->vertIndex() "<< (*iEleTk).vertIndex() << " (*iSimTk).vertIndex() " << (*iSimTk).vertIndex() << " (*iSimTk).type() " << (*iSimTk).type() << " (*iSimTk).trackId() " << (*iSimTk).trackId() << std::endl; 00157 00158 int vertexId1 = (*iSimTk).vertIndex(); 00159 SimVertex vertex1 = theSimVertices[vertexId1]; 00160 int vertexId2 = trLast.vertIndex(); 00161 SimVertex vertex2 = theSimVertices[vertexId2]; 00162 00163 00164 int motherId=-1; 00165 00166 if( ( vertexId1 == vertexId2 ) && ( (*iSimTk).type() == (*iEleTk).type() ) && trLast.type() == 22 ) { 00167 //std::cout << " Here a e/gamma brem vertex " << std::endl; 00168 00169 //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; 00170 00171 //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; 00172 00173 //std::cout << " Electron pt " << std::sqrt((*iSimTk).momentum().perp2()) << " photon pt " << std::sqrt(trLast.momentum().perp2()) << "Mother electron pt " << sqrt(motherMomentum.perp2()) << std::endl; 00174 //std::cout << " eleId " << eleId << std::endl; 00175 float eLoss = remainingEnergy - ( (*iSimTk).momentum() + trLast.momentum()).e(); 00176 //std::cout << " eLoss " << eLoss << std::endl; 00177 00178 00179 if ( vertex1.parentIndex() ) { 00180 00181 unsigned motherGeantId = vertex1.parentIndex(); 00182 std::map<unsigned, unsigned >::iterator association = geantToIndex_.find( motherGeantId ); 00183 if(association != geantToIndex_.end() ) 00184 motherId = association->second; 00185 00186 int motherType = motherId == -1 ? 0 : theSimTracks[motherId].type(); 00187 //std::cout << " Parent to this vertex motherId " << motherId << " mother type " << motherType << " Sim track ID " << theSimTracks[motherId].trackId() << std::endl; 00188 if ( theSimTracks[motherId].trackId() == eleId ) { 00189 00190 //std::cout << " ***** Found the Initial Mother Electron **** theSimTracks[motherId].trackId() " << theSimTracks[motherId].trackId() << " eleId " << eleId << std::endl; 00191 eleId= (*iSimTk).trackId(); 00192 remainingEnergy = (*iSimTk).momentum().e(); 00193 motherMomentum = (*iSimTk).momentum(); 00194 00195 00196 pBrem.push_back( HepLorentzVector(trLast.momentum().px(),trLast.momentum().py(), 00197 trLast.momentum().pz(),trLast.momentum().e()) ); 00198 bremPos.push_back( HepLorentzVector(vertex1.position().x(),vertex1.position().y(), 00199 vertex1.position().z(),vertex1.position().t()) ); 00200 xBrem.push_back(eLoss); 00201 00202 } 00203 00204 00205 00206 00207 } else { 00208 //std::cout << " This vertex has no parent tracks " << std::endl; 00209 } 00210 00211 } 00212 trLast=(*iSimTk); 00213 00214 } // End loop over all SimTracks 00215 //std::cout << " Going to build the ElectronMCTruth: pBrem size " << pBrem.size() << std::endl; 00217 HepLorentzVector tmpEleMom(primEleMom.px(),primEleMom.py(), 00218 primEleMom.pz(),primEleMom.e() ) ; 00219 HepLorentzVector tmpVtxPos(primVtxPos.x(),primVtxPos.y(),primVtxPos.z(),primVtxPos.t()); 00220 result.push_back ( ElectronMCTruth( tmpEleMom, eleVtxIndex, bremPos, pBrem, xBrem, tmpVtxPos,(*iEleTk) ) ) ; 00221 00222 00223 } // End loop over primary electrons 00224 00225 00226 00227 return result; 00228 }
std::map<unsigned, unsigned> ElectronMCTruthFinder::geantToIndex_ [private] |
Definition at line 29 of file ElectronMCTruthFinder.h.