23 std::vector<PhotonMCTruth>
result;
32 const int ELECTRON_FLAV=1;
33 const int PIZERO_FLAV=2;
34 const int PHOTON_FLAV=3;
40 std::vector<SimTrack*> photonTracks;
42 std::vector<SimTrack> trkFromConversion;
43 std::vector<ElectronMCTruth> electronsFromConversions;
47 fill(theSimTracks, theSimVertices);
50 if ( theSimTracks.size() != 0 ) {
55 std::vector<SimTrack>::const_iterator iFirstSimTk = theSimTracks.begin();
56 if ( !(*iFirstSimTk).noVertex() ) {
57 iPV = (*iFirstSimTk).vertIndex();
59 int vtxId = (*iFirstSimTk).vertIndex();
60 primVtx = theSimVertices[vtxId];
64 partType1 = (*iFirstSimTk).type();
81 if ( iFirstSimTk!= theSimTracks.end() ) {
83 if ( (*iFirstSimTk).vertIndex() == iPV) {
84 partType2 = (*iFirstSimTk).type();
98 for (std::vector<SimTrack>::const_iterator iSimTk = theSimTracks.begin(); iSimTk != theSimTracks.end(); ++iSimTk){
99 if ( (*iSimTk).noVertex() )
continue;
106 if ( (*iSimTk).vertIndex() == iPV ) {
108 if ( (*iSimTk).type() == 22) {
111 photonTracks.push_back(&(const_cast<SimTrack&>(*iSimTk)));
123 }
else if(npv == 1) {
125 ievtype = SINGLE; ievflav = ELECTRON_FLAV;
126 }
else if(partType1 == 111) {
127 ievtype = SINGLE; ievflav = PIZERO_FLAV;
128 }
else if(partType1 == 22) {
129 ievtype = SINGLE; ievflav = PHOTON_FLAV;
131 }
else if(npv == 2) {
133 ievtype = DOUBLE; ievflav = ELECTRON_FLAV;
134 }
else if(partType1 == 111 && partType2 == 111) {
135 ievtype = DOUBLE; ievflav = PIZERO_FLAV;
136 }
else if(partType1 == 22 && partType2 == 22) {
137 ievtype = DOUBLE; ievflav = PHOTON_FLAV;
144 int phoMotherType=-1;
145 int phoMotherVtxIndex=-1;
147 if( ievflav == PHOTON_FLAV || ievflav== PIZERO_FLAV || ievtype == PYTHIA ) {
155 for (std::vector<SimTrack>::const_iterator iPhoTk = theSimTracks.begin(); iPhoTk != theSimTracks.end(); ++iPhoTk){
157 trkFromConversion.clear();
158 electronsFromConversions.clear();
160 if ( (*iPhoTk).type() != 22 )
continue;
161 int photonVertexIndex= (*iPhoTk).vertIndex();
162 int phoTrkId= (*iPhoTk).trackId();
166 SimVertex vertex = theSimVertices[ photonVertexIndex];
173 phoMotherId = association->second;
174 phoMotherType = phoMotherId == -1 ? 0 : theSimTracks[phoMotherId].type();
178 if ( phoMotherType == 111 || phoMotherType == 221 || phoMotherType == 331 ) {
188 for (std::vector<SimTrack>::const_iterator iEleTk = theSimTracks.begin(); iEleTk != theSimTracks.end(); ++iEleTk){
189 if ( (*iEleTk).noVertex() )
continue;
190 if ( (*iEleTk).vertIndex() == iPV )
continue;
191 if (
std::abs((*iEleTk).type()) != 11 )
continue;
193 int vertexId = (*iEleTk).vertIndex();
194 SimVertex vertex = theSimVertices[vertexId];
202 std::map<unsigned, unsigned >::iterator association =
geantToIndex_.find( motherGeantId );
204 motherId = association->second;
211 std::vector<CLHEP::Hep3Vector> bremPos;
212 std::vector<CLHEP::HepLorentzVector> pBrem;
213 std::vector<float> xBrem;
215 if ( theSimTracks[motherId].trackId() == (*iPhoTk).trackId() ) {
221 trkFromConversion.push_back( (*iEleTk ) );
223 float remainingEnergy =trLast.
momentum().e();
227 (*iEleTk).momentum().y(),
228 (*iEleTk).momentum().z(),
229 (*iEleTk).momentum().e());
231 unsigned int eleId = (*iEleTk).trackId();
232 int eleVtxIndex= (*iEleTk).vertIndex();
239 for (std::vector<SimTrack>::const_iterator iSimTk = theSimTracks.begin(); iSimTk != theSimTracks.end(); ++iSimTk){
241 if ( (*iSimTk).noVertex() )
continue;
242 if ( (*iSimTk).vertIndex() == iPV )
continue;
246 int vertexId1 = (*iSimTk).vertIndex();
247 SimVertex vertex1 = theSimVertices[vertexId1];
254 if( ( vertexId1 == vertexId2 ) && ( (*iSimTk).type() == (*iEleTk).type() ) && trLast.
type() == 22 ) {
263 float eLoss = remainingEnergy - ( (*iSimTk).momentum() + trLast.
momentum()).
e();
270 std::map<unsigned, unsigned >::iterator association =
geantToIndex_.find( motherGeantId );
272 motherId = association->second;
276 if (theSimTracks[motherId].trackId() == eleId ) {
279 eleId= (*iSimTk).trackId();
280 remainingEnergy = (*iSimTk).momentum().e();
281 motherMomentum = (*iSimTk).momentum();
284 pBrem.push_back( CLHEP::HepLorentzVector(trLast.
momentum().px(),
288 bremPos.push_back( CLHEP::HepLorentzVector(vertex1.
position().x(),
292 xBrem.push_back(eLoss);
310 CLHEP::HepLorentzVector tmpEleMom(primEleMom.px(),primEleMom.py(),
311 primEleMom.pz(),primEleMom.e() );
312 CLHEP::HepLorentzVector tmpVtxPos(primVtxPos.x(),primVtxPos.y(),
313 primVtxPos.z(),primVtxPos.t() );
314 electronsFromConversions.push_back (
316 xBrem, tmpVtxPos , const_cast<SimTrack&>(*iEleTk) ) ) ;
332 CLHEP::HepLorentzVector phoMotherMom(0.,0.,0.,0.);
333 CLHEP::HepLorentzVector phoMotherVtx(0.,0.,0.,0.);
335 if ( phoMotherId >= 0) {
337 phoMotherVtxIndex = theSimTracks[phoMotherId].vertIndex();
338 SimVertex motherVtx = theSimVertices[ phoMotherVtxIndex];
344 phoMotherMom.setPx( theSimTracks[phoMotherId].momentum().
x());
345 phoMotherMom.setPy( theSimTracks[phoMotherId].momentum().
y());
346 phoMotherMom.setPz( theSimTracks[phoMotherId].momentum().
z() );
347 phoMotherMom.setE( theSimTracks[phoMotherId].momentum().
t());
350 phoMotherVtx.setX ( motherVtxPosition.x());
351 phoMotherVtx.setY ( motherVtxPosition.y());
352 phoMotherVtx.setZ ( motherVtxPosition.z());
353 phoMotherVtx.setT ( motherVtxPosition.t());
358 if ( electronsFromConversions.size() > 0 ) {
364 int convVtxId =electronsFromConversions[0].vertexInd();
365 SimVertex convVtx = theSimVertices[convVtxId];
374 CLHEP::HepLorentzVector tmpPhoMom( (*iPhoTk).momentum().px(), (*iPhoTk).momentum().py(),
375 (*iPhoTk).momentum().pz(), (*iPhoTk).momentum().e() ) ;
376 CLHEP::HepLorentzVector tmpVertex( vtxPosition.x(), vtxPosition.y(), vtxPosition.z(), vtxPosition.t() );
377 CLHEP::HepLorentzVector tmpPrimVtx( primVtxPos.x(), primVtxPos.y(), primVtxPos.z(), primVtxPos.t() ) ;
381 result.push_back(
PhotonMCTruth(isAconversion, tmpPhoMom, photonVertexIndex, phoTrkId, phoMotherType,phoMotherMom, phoMotherVtx, tmpVertex,
382 tmpPrimVtx, electronsFromConversions ));
387 CLHEP::HepLorentzVector vtxPosition(0.,0.,0.,0.);
388 CLHEP::HepLorentzVector tmpPhoMom( (*iPhoTk).momentum().px(), (*iPhoTk).momentum().py(),
389 (*iPhoTk).momentum().pz(), (*iPhoTk).momentum().e() ) ;
390 CLHEP::HepLorentzVector tmpPrimVtx( primVtxPos.x(), primVtxPos.y(), primVtxPos.z(), primVtxPos.t() ) ;
392 result.push_back(
PhotonMCTruth(isAconversion, tmpPhoMom, photonVertexIndex, phoTrkId, phoMotherType, phoMotherMom, phoMotherVtx, vtxPosition,
393 tmpPrimVtx, electronsFromConversions ));
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > XYZTLorentzVectorD
Lorentz vector with cylindrical internal representation using pseudorapidity.
void fill(const std::vector< SimTrack > &theSimTracks, const std::vector< SimVertex > &theSimVertices)
std::map< unsigned, unsigned > geantToIndex_
Abs< T >::type abs(const T &t)
const math::XYZTLorentzVectorD & position() const
int vertIndex() const
index of the vertex in the Event container (-1 if no vertex)
int type() const
particle type (HEP PDT convension)
const math::XYZTLorentzVectorD & momentum() const