CMS 3D CMS Logo
PhotonMCTruthFinder.cc
Go to the documentation of this file.
1 #include "RecoEgamma/EgammaMCTools/interface/PhotonMCTruthFinder.h"
2 #include "RecoEgamma/EgammaMCTools/interface/PhotonMCTruth.h"
3 #include "RecoEgamma/EgammaMCTools/interface/ElectronMCTruth.h"
4 //
5 //
6 #include "SimDataFormats/CrossingFrame/interface/CrossingFrame.h"
7 #include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
8 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
9 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h"
10 
11 #include <algorithm>
12 
13 PhotonMCTruthFinder::PhotonMCTruthFinder() {
14  //std::cout << " PhotonMCTruthFinder CTOR " << std::endl;
15 }
16 
17 std::vector<PhotonMCTruth> PhotonMCTruthFinder::find(const std::vector<SimTrack>& theSimTracks,
18  const std::vector<SimVertex>& theSimVertices) {
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  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  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  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  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  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 }
380 
381 void PhotonMCTruthFinder::fill(const std::vector<SimTrack>& simTracks, const std::vector<SimVertex>& simVertices) {
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 }
DDAxes::y
CoreSimTrack::momentum
const math::XYZTLorentzVectorD & momentum() const
Definition: CoreSimTrack.h:19
SimVertex
Definition: SimVertex.h:5
SimVertex::parentIndex
int parentIndex() const
Definition: SimVertex.h:29
math::XYZTLorentzVectorD
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > XYZTLorentzVectorD
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:14
CrossingFrame.h
PhotonMCTruthFinder::PhotonMCTruthFinder
PhotonMCTruthFinder()
Definition: PhotonMCTruthFinder.cc:13
TrackCandidateProducer_cfi.simTracks
simTracks
Definition: TrackCandidateProducer_cfi.py:15
DDAxes::x
PhotonMCTruthFinder::find
std::vector< PhotonMCTruth > find(const std::vector< SimTrack > &simTracks, const std::vector< SimVertex > &simVertices)
Definition: PhotonMCTruthFinder.cc:17
simKBmtfDigis_cfi.eLoss
eLoss
Definition: simKBmtfDigis_cfi.py:9
TrackingVertexContainer.h
ElectronMCTruth.h
MixCollection.h
PhotonMCTruth.h
DDAxes::z
PhotonMCTruth
Definition: PhotonMCTruth.h:20
OrderedSet.t
t
Definition: OrderedSet.py:90
PhotonMCTruthFinder.h
bphysicsOniaDQM_cfi.vertex
vertex
Definition: bphysicsOniaDQM_cfi.py:7
CoreSimTrack::type
int type() const
particle type (HEP PDT convension)
Definition: CoreSimTrack.h:22
CoreSimVertex::position
const math::XYZTLorentzVectorD & position() const
Definition: CoreSimVertex.h:21
TrackValidation_cff.association
association
Definition: TrackValidation_cff.py:286
PhotonMCTruthFinder::fill
void fill(const std::vector< SimTrack > &theSimTracks, const std::vector< SimVertex > &theSimVertices)
Definition: PhotonMCTruthFinder.cc:381
objects.autophobj.motherId
motherId
Definition: autophobj.py:237
TrackingParticleFwd.h
SimTrack
Definition: SimTrack.h:6
ElectronMCTruth
Definition: ElectronMCTruth.h:18
SimTrack::vertIndex
int vertIndex() const
index of the vertex in the Event container (-1 if no vertex)
Definition: SimTrack.h:30
mps_fire.result
result
Definition: mps_fire.py:303
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
HGCalValidator_cfi.simVertices
simVertices
Definition: HGCalValidator_cfi.py:43
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
PhotonMCTruthFinder::geantToIndex_
std::map< unsigned, unsigned > geantToIndex_
Definition: PhotonMCTruthFinder.h:31
Generated for CMSSW Reference Manual by doxygen 1.8.16