CMS 3D CMS Logo

List of all members | Public Member Functions | Private Member Functions | Private Attributes
PhotonMCTruthFinder Class Reference

#include <PhotonMCTruthFinder.h>

Public Member Functions

void clear ()
 
std::vector< PhotonMCTruthfind (const std::vector< SimTrack > &simTracks, const std::vector< SimVertex > &simVertices)
 
 PhotonMCTruthFinder ()
 
virtual ~PhotonMCTruthFinder ()
 

Private Member Functions

void fill (const std::vector< SimTrack > &theSimTracks, const std::vector< SimVertex > &theSimVertices)
 

Private Attributes

std::map< unsigned, unsigned > geantToIndex_
 

Detailed Description

Author
N. Marinelli Notre Dame

Definition at line 20 of file PhotonMCTruthFinder.h.

Constructor & Destructor Documentation

◆ PhotonMCTruthFinder()

PhotonMCTruthFinder::PhotonMCTruthFinder ( )

Definition at line 13 of file PhotonMCTruthFinder.cc.

13  {
14  //std::cout << " PhotonMCTruthFinder CTOR " << std::endl;
15 }

◆ ~PhotonMCTruthFinder()

virtual PhotonMCTruthFinder::~PhotonMCTruthFinder ( )
inlinevirtual

Definition at line 23 of file PhotonMCTruthFinder.h.

23 {}

Member Function Documentation

◆ clear()

void PhotonMCTruthFinder::clear ( void  )
inline

Definition at line 27 of file PhotonMCTruthFinder.h.

References geantToIndex_.

Referenced by TkConvValidator::analyze().

27 { geantToIndex_.clear(); }
std::map< unsigned, unsigned > geantToIndex_

◆ fill()

void PhotonMCTruthFinder::fill ( const std::vector< SimTrack > &  theSimTracks,
const std::vector< SimVertex > &  theSimVertices 
)
private

Definition at line 381 of file PhotonMCTruthFinder.cc.

References geantToIndex_, L1TkPrimaryVertex_cfi::nVtx, TrackCandidateProducer_cfi::simTracks, and HGCalValidator_cfi::simVertices.

Referenced by find().

381  {
382  // std::cout << " PhotonMCTruthFinder::fill " << std::endl;
383 
384  // Watch out there ! A SimVertex is in mm (stupid),
385 
386  unsigned nVtx = simVertices.size();
387  unsigned nTks = simTracks.size();
388 
389  // std::cout << " PhotonMCTruthFinder::fill " << nVtx << " " << nTks << std::endl;
390 
391  // Empty event, do nothin'
392  if (nVtx == 0)
393  return;
394 
395  // create a map associating geant particle id and position in the
396  // event SimTrack vector
397  for (unsigned it = 0; it < nTks; ++it) {
398  geantToIndex_[simTracks[it].trackId()] = it;
399  // std::cout << " PhotonMCTruthFinder::fill it " << it << " simTracks[it].trackId() " << simTracks[it].trackId() << std::endl;
400  }
401 }
std::map< unsigned, unsigned > geantToIndex_

◆ find()

std::vector< PhotonMCTruth > PhotonMCTruthFinder::find ( const std::vector< SimTrack > &  simTracks,
const std::vector< SimVertex > &  simVertices 
)

find truth about this electron and store it since it's from a converted photon

here fill the electron

Definition at line 17 of file PhotonMCTruthFinder.cc.

References funct::abs(), MillePedeFileConverter_cfg::e, simKBmtfDigis_cfi::eLoss, fill(), geantToIndex_, caHitNtupletGeneratorKernels::if(), CoreSimTrack::momentum(), objects.autophobj::motherId, SimVertex::parentIndex(), CoreSimVertex::position(), mps_fire::result, submitPVValidationJobs::t, CoreSimTrack::type(), bphysicsOniaDQM_cfi::vertex, SimTrack::vertIndex(), x, y, and z.

Referenced by TkConvValidator::analyze(), and PizeroMCTruthFinder::find().

18  {
19  // std::cout << " PhotonMCTruthFinder::find " << std::endl;
20 
21  std::vector<PhotonMCTruth> result;
22 
23  //const float pi = 3.141592653592;
24  //const float twopi=2*pi;
25 
26  // Local variables
27  const int SINGLE = 1;
28  const int DOUBLE = 2;
29  const int PYTHIA = 3;
30  const int ELECTRON_FLAV = 1;
31  const int PIZERO_FLAV = 2;
32  const int PHOTON_FLAV = 3;
33 
34  int ievtype = 0;
35  int ievflav = 0;
36 
37  std::vector<SimTrack*> photonTracks;
38 
39  std::vector<SimTrack> trkFromConversion;
40  std::vector<ElectronMCTruth> electronsFromConversions;
41  SimVertex primVtx;
42 
43  fill(theSimTracks, theSimVertices);
44 
45  // std::cout << " After fill " << theSimTracks.size() << " " << theSimVertices.size() << std::endl;
46  if (!theSimTracks.empty()) {
47  int iPV = -1;
48  int partType1 = 0;
49  int partType2 = 0;
50  std::vector<SimTrack>::const_iterator iFirstSimTk = theSimTracks.begin();
51  if (!(*iFirstSimTk).noVertex()) {
52  iPV = (*iFirstSimTk).vertIndex();
53 
54  int vtxId = (*iFirstSimTk).vertIndex();
55  primVtx = theSimVertices[vtxId];
56 
57  partType1 = (*iFirstSimTk).type();
58  // std::cout << " Primary vertex id " << iPV << " first track type " << (*iFirstSimTk).type() << std::endl;
59  } else {
60  //std::cout << " First track has no vertex " << std::endl;
61  }
62 
63  math::XYZTLorentzVectorD primVtxPos(
64  primVtx.position().x(), primVtx.position().y(), primVtx.position().z(), primVtx.position().e());
65 
66  // Look at a second track
67  iFirstSimTk++;
68  if (iFirstSimTk != theSimTracks.end()) {
69  if ((*iFirstSimTk).vertIndex() == iPV) {
70  partType2 = (*iFirstSimTk).type();
71  // std::cout << " second track type " << (*iFirstSimTk).type() << " vertex " << (*iFirstSimTk).vertIndex() << std::endl;
72 
73  } else {
74  // std::cout << " Only one kine track from Primary Vertex " << std::endl;
75  }
76  }
77 
78  //std::cout << " Loop over all particles " << std::endl;
79 
80  int npv = 0;
81  //int iPho=0;
82  //int iPizero=0;
83  // theSimTracks.reset();
84  for (std::vector<SimTrack>::const_iterator iSimTk = theSimTracks.begin(); iSimTk != theSimTracks.end(); ++iSimTk) {
85  if ((*iSimTk).noVertex())
86  continue;
87 
88  //int vertexId = (*iSimTk).vertIndex();
89  //SimVertex vertex = theSimVertices[vertexId];
90 
91  // std::cout << " Particle type " << (*iSimTk).type() << " Sim Track ID " << (*iSimTk).trackId() << " momentum " << (*iSimTk).momentum() << " vertex position " << vertex.position() << " vertex index " << (*iSimTk).vertIndex() << std::endl;
92  if ((*iSimTk).vertIndex() == iPV) {
93  npv++;
94  if ((*iSimTk).type() == 22) {
95  // std::cout << " Found a primary photon with ID " << (*iSimTk).trackId() << " momentum " << (*iSimTk).momentum() << std::endl;
96 
97  photonTracks.push_back(&(const_cast<SimTrack&>(*iSimTk)));
98  }
99  }
100  }
101 
102  // std::cout << " There are " << npv << " particles originating in the PV " << std::endl;
103 
104  if (npv >= 3) {
105  ievtype = PYTHIA;
106  } else if (npv == 1) {
107  if (std::abs(partType1) == 11) {
108  ievtype = SINGLE;
109  ievflav = ELECTRON_FLAV;
110  } else if (partType1 == 111) {
111  ievtype = SINGLE;
112  ievflav = PIZERO_FLAV;
113  } else if (partType1 == 22) {
114  ievtype = SINGLE;
115  ievflav = PHOTON_FLAV;
116  }
117  } else if (npv == 2) {
118  if (std::abs(partType1) == 11 && std::abs(partType2) == 11) {
119  ievtype = DOUBLE;
120  ievflav = ELECTRON_FLAV;
121  } else if (partType1 == 111 && partType2 == 111) {
122  ievtype = DOUBLE;
123  ievflav = PIZERO_FLAV;
124  } else if (partType1 == 22 && partType2 == 22) {
125  ievtype = DOUBLE;
126  ievflav = PHOTON_FLAV;
127  }
128  }
129 
131 
132  int isAconversion = 0;
133  int phoMotherType = -1;
134  int phoMotherVtxIndex = -1;
135  int phoMotherId = -1;
136  if (ievflav == PHOTON_FLAV || ievflav == PIZERO_FLAV || ievtype == PYTHIA) {
137  // std::cout << " It's a primary PHOTON or PIZERO or PYTHIA event with " << photonTracks.size() << " photons ievtype " << ievtype << " ievflav " << ievflav<< std::endl;
138 
139  // for (std::vector<SimTrack*>::iterator iPhoTk = photonTracks.begin(); iPhoTk != photonTracks.end(); ++iPhoTk){
140  // std::cout << " All gamma found from PV " << (*iPhoTk)->momentum() << " photon track ID " << (*iPhoTk)->trackId() << " vertex index " << (*iPhoTk)->vertIndex() << std::endl;
141  // }
142 
143  for (std::vector<SimTrack>::const_iterator iPhoTk = theSimTracks.begin(); iPhoTk != theSimTracks.end();
144  ++iPhoTk) {
145  trkFromConversion.clear();
146  electronsFromConversions.clear();
147 
148  if ((*iPhoTk).type() != 22)
149  continue;
150  int photonVertexIndex = (*iPhoTk).vertIndex();
151  int phoTrkId = (*iPhoTk).trackId();
152  //std::cout << " Looping on gamma looking for conversions " << (*iPhoTk).momentum() << " photon track ID " << (*iPhoTk).trackId() << std::endl;
153 
154  // check who is his mother
155  const SimVertex& vertex = theSimVertices[photonVertexIndex];
156  phoMotherId = -1;
157  if (vertex.parentIndex() != -1) {
158  unsigned motherGeantId = vertex.parentIndex();
159  std::map<unsigned, unsigned>::iterator association = geantToIndex_.find(motherGeantId);
160  if (association != geantToIndex_.end())
161  phoMotherId = association->second;
162  phoMotherType = phoMotherId == -1 ? 0 : theSimTracks[phoMotherId].type();
163 
164  if (phoMotherType == 111 || phoMotherType == 221 || phoMotherType == 331) {
165  //std::cout << " Parent to this vertex motherId " << phoMotherId << " mother type " << phoMotherType << " Sim track ID " << theSimTracks[phoMotherId].trackId() << std::endl;
166  //std::cout << " Son of a pizero or eta " << phoMotherType << std::endl;
167  }
168  }
169 
170  for (std::vector<SimTrack>::const_iterator iEleTk = theSimTracks.begin(); iEleTk != theSimTracks.end();
171  ++iEleTk) {
172  if ((*iEleTk).noVertex())
173  continue;
174  if ((*iEleTk).vertIndex() == iPV)
175  continue;
176  if (std::abs((*iEleTk).type()) != 11)
177  continue;
178 
179  int vertexId = (*iEleTk).vertIndex();
180  const SimVertex& vertex = theSimVertices[vertexId];
181  int motherId = -1;
182 
183  //std::cout << " Secondary from photons particle type " << (*iEleTk).type() << " trackId " << (*iEleTk).trackId() << " vertex ID " << vertexId << std::endl;
184  if (vertex.parentIndex() != -1) {
185  unsigned motherGeantId = vertex.parentIndex();
186  std::map<unsigned, unsigned>::iterator association = geantToIndex_.find(motherGeantId);
187  if (association != geantToIndex_.end())
188  motherId = association->second;
189 
190  //int motherType = motherId == -1 ? 0 : theSimTracks[motherId].type();
191 
192  //std::cout << " Parent to this vertex motherId " << motherId << " mother type " << motherType << " Sim track ID " << theSimTracks[motherId].trackId() << std::endl;
193 
194  std::vector<CLHEP::Hep3Vector> bremPos;
195  std::vector<CLHEP::HepLorentzVector> pBrem;
196  std::vector<float> xBrem;
197 
198  if (theSimTracks[motherId].trackId() == (*iPhoTk).trackId()) {
199  //std::cout << " Found the Mother Photon " << std::endl;
201 
202  trkFromConversion.push_back((*iEleTk));
203  SimTrack trLast = (*iEleTk);
204  float remainingEnergy = trLast.momentum().e();
205  //HepLorentzVector primEleMom=(*iEleTk).momentum();
206  //HepLorentzVector motherMomentum=(*iEleTk).momentum();
207  math::XYZTLorentzVectorD primEleMom((*iEleTk).momentum().x(),
208  (*iEleTk).momentum().y(),
209  (*iEleTk).momentum().z(),
210  (*iEleTk).momentum().e());
211  math::XYZTLorentzVectorD motherMomentum(primEleMom);
212  unsigned int eleId = (*iEleTk).trackId();
213  int eleVtxIndex = (*iEleTk).vertIndex();
214 
215  bremPos.clear();
216  pBrem.clear();
217  xBrem.clear();
218 
219  for (std::vector<SimTrack>::const_iterator iSimTk = theSimTracks.begin(); iSimTk != theSimTracks.end();
220  ++iSimTk) {
221  if ((*iSimTk).noVertex())
222  continue;
223  if ((*iSimTk).vertIndex() == iPV)
224  continue;
225 
226  //std::cout << " (*iEleTk)->trackId() " << (*iEleTk).trackId() << " (*iEleTk)->vertIndex() "<< (*iEleTk).vertIndex() << " (*iSimTk).vertIndex() " << (*iSimTk).vertIndex() << " (*iSimTk).type() " << (*iSimTk).type() << " (*iSimTk).trackId() " << (*iSimTk).trackId() << std::endl;
227 
228  int vertexId1 = (*iSimTk).vertIndex();
229  const SimVertex& vertex1 = theSimVertices[vertexId1];
230  int vertexId2 = trLast.vertIndex();
231  //SimVertex vertex2 = theSimVertices[vertexId2];
232 
233  int motherId = -1;
234 
235  if ((vertexId1 == vertexId2) && ((*iSimTk).type() == (*iEleTk).type()) && trLast.type() == 22) {
236  //std::cout << " Here a e/gamma brem vertex " << std::endl;
237 
238  //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;
239 
240  //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;
241 
242  //std::cout << " Electron pt " << sqrt((*iSimTk).momentum().perp2()) << " photon pt " << sqrt(trLast.momentum().perp2()) << " Mother electron pt " << sqrt(motherMomentum.perp2()) << std::endl;
243  //std::cout << " eleId " << eleId << std::endl;
244  float eLoss = remainingEnergy - ((*iSimTk).momentum() + trLast.momentum()).e();
245  //std::cout << " eLoss " << eLoss << std::endl;
246 
247  if (vertex1.parentIndex() != -1) {
248  unsigned motherGeantId = vertex1.parentIndex();
249  std::map<unsigned, unsigned>::iterator association = geantToIndex_.find(motherGeantId);
250  if (association != geantToIndex_.end())
251  motherId = association->second;
252 
253  //int motherType = motherId == -1 ? 0 : theSimTracks[motherId].type();
254  //std::cout << " Parent to this vertex motherId " << motherId << " mother type " << motherType << " Sim track ID " << theSimTracks[motherId].trackId() << std::endl;
255  if (theSimTracks[motherId].trackId() == eleId) {
256  //std::cout << " ***** Found the Initial Mother Electron **** theSimTracks[motherId].trackId() " << theSimTracks[motherId].trackId() << " eleId " << eleId << std::endl;
257  eleId = (*iSimTk).trackId();
258  remainingEnergy = (*iSimTk).momentum().e();
259  motherMomentum = (*iSimTk).momentum();
260 
261  pBrem.push_back(CLHEP::HepLorentzVector(trLast.momentum().px(),
262  trLast.momentum().py(),
263  trLast.momentum().pz(),
264  trLast.momentum().e()));
265  bremPos.push_back(CLHEP::HepLorentzVector(vertex1.position().x(),
266  vertex1.position().y(),
267  vertex1.position().z(),
268  vertex1.position().t()));
269  xBrem.push_back(eLoss);
270  }
271 
272  } else {
273  //std::cout << " This vertex has no parent tracks " << std::endl;
274  }
275  }
276  trLast = (*iSimTk);
277 
278  } // End loop over all SimTracks
279  //std::cout << " Going to build the ElectronMCTruth for this electron from converted photon: pBrem size " << pBrem.size() << std::endl;
281 
282  CLHEP::HepLorentzVector tmpEleMom(primEleMom.px(), primEleMom.py(), primEleMom.pz(), primEleMom.e());
283  CLHEP::HepLorentzVector tmpVtxPos(primVtxPos.x(), primVtxPos.y(), primVtxPos.z(), primVtxPos.t());
284  electronsFromConversions.push_back(ElectronMCTruth(
285  tmpEleMom, eleVtxIndex, bremPos, pBrem, xBrem, tmpVtxPos, const_cast<SimTrack&>(*iEleTk)));
286  }
287 
288  } else {
289  //std::cout << " This vertex has no parent tracks " << std::endl;
290  }
291 
292  } // End of loop over the SimTracks
293 
294  //std::cout << " DEBUG trkFromConversion.size() " << trkFromConversion.size() << " electronsFromConversions.size() " << electronsFromConversions.size() << std::endl;
295 
296  math::XYZTLorentzVectorD motherVtxPosition(0., 0., 0., 0.);
297  CLHEP::HepLorentzVector phoMotherMom(0., 0., 0., 0.);
298  CLHEP::HepLorentzVector phoMotherVtx(0., 0., 0., 0.);
299 
300  if (phoMotherId >= 0) {
301  phoMotherVtxIndex = theSimTracks[phoMotherId].vertIndex();
302  const SimVertex& motherVtx = theSimVertices[phoMotherVtxIndex];
303  motherVtxPosition = math::XYZTLorentzVectorD(
304  motherVtx.position().x(), motherVtx.position().y(), motherVtx.position().z(), motherVtx.position().e());
305 
306  phoMotherMom.setPx(theSimTracks[phoMotherId].momentum().x());
307  phoMotherMom.setPy(theSimTracks[phoMotherId].momentum().y());
308  phoMotherMom.setPz(theSimTracks[phoMotherId].momentum().z());
309  phoMotherMom.setE(theSimTracks[phoMotherId].momentum().t());
310  // std::cout << " PhotonMCTruthFinder mother " << phoMotherId << " type " << phoMotherType << " Momentum" << phoMotherMom.et() << std::endl;
311 
312  phoMotherVtx.setX(motherVtxPosition.x());
313  phoMotherVtx.setY(motherVtxPosition.y());
314  phoMotherVtx.setZ(motherVtxPosition.z());
315  phoMotherVtx.setT(motherVtxPosition.t());
316  }
317 
318  if (!electronsFromConversions.empty()) {
319  // if ( trkFromConversion.size() > 0 ) {
320  isAconversion = 1;
321  //std::cout << " CONVERTED photon " << "\n";
322 
323  //int convVtxId = trkFromConversion[0].vertIndex();
324  int convVtxId = electronsFromConversions[0].vertexInd();
325  const SimVertex& convVtx = theSimVertices[convVtxId];
326  // CLHEP::HepLorentzVector vtxPosition = convVtx.position();
327  math::XYZTLorentzVectorD vtxPosition(
328  convVtx.position().x(), convVtx.position().y(), convVtx.position().z(), convVtx.position().e());
329 
330  //result.push_back( PhotonMCTruth(isAconversion, (*iPhoTk).momentum(), photonVertexIndex, phoTrkId, vtxPosition, primVtx.position(), trkFromConversion ));
331  CLHEP::HepLorentzVector tmpPhoMom((*iPhoTk).momentum().px(),
332  (*iPhoTk).momentum().py(),
333  (*iPhoTk).momentum().pz(),
334  (*iPhoTk).momentum().e());
335  CLHEP::HepLorentzVector tmpVertex(vtxPosition.x(), vtxPosition.y(), vtxPosition.z(), vtxPosition.t());
336  CLHEP::HepLorentzVector tmpPrimVtx(primVtxPos.x(), primVtxPos.y(), primVtxPos.z(), primVtxPos.t());
337 
338  result.push_back(PhotonMCTruth(isAconversion,
339  tmpPhoMom,
340  photonVertexIndex,
341  phoTrkId,
342  phoMotherType,
343  phoMotherMom,
344  phoMotherVtx,
345  tmpVertex,
346  tmpPrimVtx,
347  electronsFromConversions));
348 
349  } else {
350  isAconversion = 0;
351  //std::cout << " UNCONVERTED photon " << "\n";
352  CLHEP::HepLorentzVector vtxPosition(0., 0., 0., 0.);
353  CLHEP::HepLorentzVector tmpPhoMom((*iPhoTk).momentum().px(),
354  (*iPhoTk).momentum().py(),
355  (*iPhoTk).momentum().pz(),
356  (*iPhoTk).momentum().e());
357  CLHEP::HepLorentzVector tmpPrimVtx(primVtxPos.x(), primVtxPos.y(), primVtxPos.z(), primVtxPos.t());
358  // result.push_back( PhotonMCTruth(isAconversion, (*iPhoTk).momentum(), photonVertexIndex, phoTrkId, vtxPosition, primVtx.position(), trkFromConversion ));
359  result.push_back(PhotonMCTruth(isAconversion,
360  tmpPhoMom,
361  photonVertexIndex,
362  phoTrkId,
363  phoMotherType,
364  phoMotherMom,
365  phoMotherVtx,
366  vtxPosition,
367  tmpPrimVtx,
368  electronsFromConversions));
369  }
370 
371  } // End loop over the primary photons
372 
373  } // Event with one or two photons
374 
375  //std::cout << " PhotonMCTruthFinder photon size " << result.size() << std::endl;
376  }
377 
378  return result;
379 }
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > XYZTLorentzVectorD
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:14
const math::XYZTLorentzVectorD & position() const
Definition: CoreSimVertex.h:21
const math::XYZTLorentzVectorD & momentum() const
Definition: CoreSimTrack.h:19
int type() const
particle type (HEP PDT convension)
Definition: CoreSimTrack.h:22
void fill(const std::vector< SimTrack > &theSimTracks, const std::vector< SimVertex > &theSimVertices)
int vertIndex() const
index of the vertex in the Event container (-1 if no vertex)
Definition: SimTrack.h:33
std::tuple< layerClusterToCaloParticle, caloParticleToLayerCluster > association
std::map< unsigned, unsigned > geantToIndex_
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
int parentIndex() const
Definition: SimVertex.h:29

Member Data Documentation

◆ geantToIndex_

std::map<unsigned, unsigned> PhotonMCTruthFinder::geantToIndex_
private

Definition at line 31 of file PhotonMCTruthFinder.h.

Referenced by clear(), fill(), and find().