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