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(std::vector<SimTrack> theSimTracks, 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>::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>::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 
112  photonTracks.push_back( &(*iSimTk) );
113 
114  }
115 
116  }
117 
118  }
119 
120  // std::cout << " There are " << npv << " particles originating in the PV " << std::endl;
121 
122  if(npv >= 3) {
123  ievtype = PYTHIA;
124  } else if(npv == 1) {
125  if( std::abs(partType1) == 11 ) {
126  ievtype = SINGLE; ievflav = ELECTRON_FLAV;
127  } else if(partType1 == 111) {
128  ievtype = SINGLE; ievflav = PIZERO_FLAV;
129  } else if(partType1 == 22) {
130  ievtype = SINGLE; ievflav = PHOTON_FLAV;
131  }
132  } else if(npv == 2) {
133  if ( std::abs(partType1) == 11 && std::abs(partType2) == 11 ) {
134  ievtype = DOUBLE; ievflav = ELECTRON_FLAV;
135  } else if(partType1 == 111 && partType2 == 111) {
136  ievtype = DOUBLE; ievflav = PIZERO_FLAV;
137  } else if(partType1 == 22 && partType2 == 22) {
138  ievtype = DOUBLE; ievflav = PHOTON_FLAV;
139  }
140  }
141 
143 
144  int isAconversion=0;
145  int phoMotherType=-1;
146  int phoMotherVtxIndex=-1;
147  int phoMotherId=-1;
148  if( ievflav == PHOTON_FLAV || ievflav== PIZERO_FLAV || ievtype == PYTHIA ) {
149 
150  // std::cout << " It's a primary PHOTON or PIZERO or PYTHIA event with " << photonTracks.size() << " photons ievtype " << ievtype << " ievflav " << ievflav<< std::endl;
151 
152  // for (std::vector<SimTrack*>::iterator iPhoTk = photonTracks.begin(); iPhoTk != photonTracks.end(); ++iPhoTk){
153  // std::cout << " All gamma found from PV " << (*iPhoTk)->momentum() << " photon track ID " << (*iPhoTk)->trackId() << " vertex index " << (*iPhoTk)->vertIndex() << std::endl;
154  // }
155 
156  for (std::vector<SimTrack>::iterator iPhoTk = theSimTracks.begin(); iPhoTk != theSimTracks.end(); ++iPhoTk){
157 
158  trkFromConversion.clear();
159  electronsFromConversions.clear();
160 
161  if ( (*iPhoTk).type() != 22 ) continue;
162  int photonVertexIndex= (*iPhoTk).vertIndex();
163  int phoTrkId= (*iPhoTk).trackId();
164  //std::cout << " Looping on gamma looking for conversions " << (*iPhoTk).momentum() << " photon track ID " << (*iPhoTk).trackId() << std::endl;
165 
166  // check who is his mother
167  SimVertex vertex = theSimVertices[ photonVertexIndex];
168  phoMotherId=-1;
169  if ( vertex.parentIndex() != -1 ) {
170 
171  unsigned motherGeantId = vertex.parentIndex();
172  std::map<unsigned, unsigned >::iterator association = geantToIndex_.find( motherGeantId );
173  if(association != geantToIndex_.end() )
174  phoMotherId = association->second;
175  phoMotherType = phoMotherId == -1 ? 0 : theSimTracks[phoMotherId].type();
176 
177 
178 
179  if ( phoMotherType == 111 || phoMotherType == 221 || phoMotherType == 331 ) {
180 
181  //std::cout << " Parent to this vertex motherId " << phoMotherId << " mother type " << phoMotherType << " Sim track ID " << theSimTracks[phoMotherId].trackId() << std::endl;
182  //std::cout << " Son of a pizero or eta " << phoMotherType << std::endl;
183  }
184 
185  }
186 
187 
188 
189  for (std::vector<SimTrack>::iterator iEleTk = theSimTracks.begin(); iEleTk != theSimTracks.end(); ++iEleTk){
190  if ( (*iEleTk).noVertex() ) continue;
191  if ( (*iEleTk).vertIndex() == iPV ) continue;
192  if ( std::abs((*iEleTk).type()) != 11 ) continue;
193 
194  int vertexId = (*iEleTk).vertIndex();
195  SimVertex vertex = theSimVertices[vertexId];
196  int motherId=-1;
197 
198 
199  //std::cout << " Secondary from photons particle type " << (*iEleTk).type() << " trackId " << (*iEleTk).trackId() << " vertex ID " << vertexId << std::endl;
200  if ( vertex.parentIndex() != -1 ) {
201 
202  unsigned motherGeantId = vertex.parentIndex();
203  std::map<unsigned, unsigned >::iterator association = geantToIndex_.find( motherGeantId );
204  if(association != geantToIndex_.end() )
205  motherId = association->second;
206 
207 
208  //int motherType = motherId == -1 ? 0 : theSimTracks[motherId].type();
209 
210  //std::cout << " Parent to this vertex motherId " << motherId << " mother type " << motherType << " Sim track ID " << theSimTracks[motherId].trackId() << std::endl;
211 
212  std::vector<CLHEP::Hep3Vector> bremPos;
213  std::vector<CLHEP::HepLorentzVector> pBrem;
214  std::vector<float> xBrem;
215 
216  if ( theSimTracks[motherId].trackId() == (*iPhoTk).trackId() ) {
217  //std::cout << " Found the Mother Photon " << std::endl;
219 
220 
221 
222  trkFromConversion.push_back( (*iEleTk ) );
223  SimTrack trLast =(*iEleTk);
224  float remainingEnergy =trLast.momentum().e();
225  //HepLorentzVector primEleMom=(*iEleTk).momentum();
226  //HepLorentzVector motherMomentum=(*iEleTk).momentum();
227  math::XYZTLorentzVectorD primEleMom((*iEleTk).momentum().x(),
228  (*iEleTk).momentum().y(),
229  (*iEleTk).momentum().z(),
230  (*iEleTk).momentum().e());
231  math::XYZTLorentzVectorD motherMomentum(primEleMom);
232  unsigned int eleId = (*iEleTk).trackId();
233  int eleVtxIndex= (*iEleTk).vertIndex();
234 
235 
236  bremPos.clear();
237  pBrem.clear();
238  xBrem.clear();
239 
240  for (std::vector<SimTrack>::iterator iSimTk = theSimTracks.begin(); iSimTk != theSimTracks.end(); ++iSimTk){
241 
242  if ( (*iSimTk).noVertex() ) continue;
243  if ( (*iSimTk).vertIndex() == iPV ) continue;
244 
245  //std::cout << " (*iEleTk)->trackId() " << (*iEleTk).trackId() << " (*iEleTk)->vertIndex() "<< (*iEleTk).vertIndex() << " (*iSimTk).vertIndex() " << (*iSimTk).vertIndex() << " (*iSimTk).type() " << (*iSimTk).type() << " (*iSimTk).trackId() " << (*iSimTk).trackId() << std::endl;
246 
247  int vertexId1 = (*iSimTk).vertIndex();
248  SimVertex vertex1 = theSimVertices[vertexId1];
249  int vertexId2 = trLast.vertIndex();
250  //SimVertex vertex2 = theSimVertices[vertexId2];
251 
252 
253  int motherId=-1;
254 
255  if( ( vertexId1 == vertexId2 ) && ( (*iSimTk).type() == (*iEleTk).type() ) && trLast.type() == 22 ) {
256  //std::cout << " Here a e/gamma brem vertex " << std::endl;
257 
258  //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;
259 
260  //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;
261 
262  //std::cout << " Electron pt " << sqrt((*iSimTk).momentum().perp2()) << " photon pt " << sqrt(trLast.momentum().perp2()) << " Mother electron pt " << sqrt(motherMomentum.perp2()) << std::endl;
263  //std::cout << " eleId " << eleId << std::endl;
264  float eLoss = remainingEnergy - ( (*iSimTk).momentum() + trLast.momentum()).e();
265  //std::cout << " eLoss " << eLoss << std::endl;
266 
267 
268  if ( vertex1.parentIndex() != -1 ) {
269 
270  unsigned motherGeantId = vertex1.parentIndex();
271  std::map<unsigned, unsigned >::iterator association = geantToIndex_.find( motherGeantId );
272  if(association != geantToIndex_.end() )
273  motherId = association->second;
274 
275  //int motherType = motherId == -1 ? 0 : theSimTracks[motherId].type();
276  //std::cout << " Parent to this vertex motherId " << motherId << " mother type " << motherType << " Sim track ID " << theSimTracks[motherId].trackId() << std::endl;
277  if (theSimTracks[motherId].trackId() == eleId ) {
278 
279  //std::cout << " ***** Found the Initial Mother Electron **** theSimTracks[motherId].trackId() " << theSimTracks[motherId].trackId() << " eleId " << eleId << std::endl;
280  eleId= (*iSimTk).trackId();
281  remainingEnergy = (*iSimTk).momentum().e();
282  motherMomentum = (*iSimTk).momentum();
283 
284 
285  pBrem.push_back( CLHEP::HepLorentzVector(trLast.momentum().px(),
286  trLast.momentum().py(),
287  trLast.momentum().pz(),
288  trLast.momentum().e()) );
289  bremPos.push_back( CLHEP::HepLorentzVector(vertex1.position().x(),
290  vertex1.position().y(),
291  vertex1.position().z(),
292  vertex1.position().t()) );
293  xBrem.push_back(eLoss);
294 
295  }
296 
297 
298 
299 
300  } else {
301  //std::cout << " This vertex has no parent tracks " << std::endl;
302  }
303 
304  }
305  trLast=(*iSimTk);
306 
307  } // End loop over all SimTracks
308  //std::cout << " Going to build the ElectronMCTruth for this electron from converted photon: pBrem size " << pBrem.size() << std::endl;
310 
311  CLHEP::HepLorentzVector tmpEleMom(primEleMom.px(),primEleMom.py(),
312  primEleMom.pz(),primEleMom.e() );
313  CLHEP::HepLorentzVector tmpVtxPos(primVtxPos.x(),primVtxPos.y(),
314  primVtxPos.z(),primVtxPos.t() );
315  electronsFromConversions.push_back (
316  ElectronMCTruth( tmpEleMom, eleVtxIndex, bremPos, pBrem,
317  xBrem, tmpVtxPos , (*iEleTk) ) ) ;
318  }
319 
320 
321 
322 
323  } else {
324  //std::cout << " This vertex has no parent tracks " << std::endl;
325  }
326 
327 
328  } // End of loop over the SimTracks
329 
330  //std::cout << " DEBUG trkFromConversion.size() " << trkFromConversion.size() << " electronsFromConversions.size() " << electronsFromConversions.size() << std::endl;
331 
332  math::XYZTLorentzVectorD motherVtxPosition(0.,0.,0.,0.);
333  CLHEP::HepLorentzVector phoMotherMom(0.,0.,0.,0.);
334  CLHEP::HepLorentzVector phoMotherVtx(0.,0.,0.,0.);
335 
336  if ( phoMotherId >= 0) {
337 
338  phoMotherVtxIndex = theSimTracks[phoMotherId].vertIndex();
339  SimVertex motherVtx = theSimVertices[ phoMotherVtxIndex];
340  motherVtxPosition =math::XYZTLorentzVectorD (motherVtx.position().x(),
341  motherVtx.position().y(),
342  motherVtx.position().z(),
343  motherVtx.position().e());
344 
345  phoMotherMom.setPx( theSimTracks[phoMotherId].momentum().x());
346  phoMotherMom.setPy( theSimTracks[phoMotherId].momentum().y());
347  phoMotherMom.setPz( theSimTracks[phoMotherId].momentum().z() );
348  phoMotherMom.setE( theSimTracks[phoMotherId].momentum().t());
349  // std::cout << " PhotonMCTruthFinder mother " << phoMotherId << " type " << phoMotherType << " Momentum" << phoMotherMom.et() << std::endl;
350 
351  phoMotherVtx.setX ( motherVtxPosition.x());
352  phoMotherVtx.setY ( motherVtxPosition.y());
353  phoMotherVtx.setZ ( motherVtxPosition.z());
354  phoMotherVtx.setT ( motherVtxPosition.t());
355 
356  }
357 
358 
359  if ( electronsFromConversions.size() > 0 ) {
360  // if ( trkFromConversion.size() > 0 ) {
361  isAconversion=1;
362  //std::cout << " CONVERTED photon " << "\n";
363 
364  //int convVtxId = trkFromConversion[0].vertIndex();
365  int convVtxId =electronsFromConversions[0].vertexInd();
366  SimVertex convVtx = theSimVertices[convVtxId];
367  // CLHEP::HepLorentzVector vtxPosition = convVtx.position();
368  math::XYZTLorentzVectorD vtxPosition(convVtx.position().x(),
369  convVtx.position().y(),
370  convVtx.position().z(),
371  convVtx.position().e());
372 
373 
374  //result.push_back( PhotonMCTruth(isAconversion, (*iPhoTk).momentum(), photonVertexIndex, phoTrkId, vtxPosition, primVtx.position(), trkFromConversion ));
375  CLHEP::HepLorentzVector tmpPhoMom( (*iPhoTk).momentum().px(), (*iPhoTk).momentum().py(),
376  (*iPhoTk).momentum().pz(), (*iPhoTk).momentum().e() ) ;
377  CLHEP::HepLorentzVector tmpVertex( vtxPosition.x(), vtxPosition.y(), vtxPosition.z(), vtxPosition.t() );
378  CLHEP::HepLorentzVector tmpPrimVtx( primVtxPos.x(), primVtxPos.y(), primVtxPos.z(), primVtxPos.t() ) ;
379 
380 
381 
382  result.push_back( PhotonMCTruth(isAconversion, tmpPhoMom, photonVertexIndex, phoTrkId, phoMotherType,phoMotherMom, phoMotherVtx, tmpVertex,
383  tmpPrimVtx, electronsFromConversions ));
384 
385  } else {
386  isAconversion=0;
387  //std::cout << " UNCONVERTED photon " << "\n";
388  CLHEP::HepLorentzVector vtxPosition(0.,0.,0.,0.);
389  CLHEP::HepLorentzVector tmpPhoMom( (*iPhoTk).momentum().px(), (*iPhoTk).momentum().py(),
390  (*iPhoTk).momentum().pz(), (*iPhoTk).momentum().e() ) ;
391  CLHEP::HepLorentzVector tmpPrimVtx( primVtxPos.x(), primVtxPos.y(), primVtxPos.z(), primVtxPos.t() ) ;
392  // result.push_back( PhotonMCTruth(isAconversion, (*iPhoTk).momentum(), photonVertexIndex, phoTrkId, vtxPosition, primVtx.position(), trkFromConversion ));
393  result.push_back( PhotonMCTruth(isAconversion, tmpPhoMom, photonVertexIndex, phoTrkId, phoMotherType, phoMotherMom, phoMotherVtx, vtxPosition,
394  tmpPrimVtx, electronsFromConversions ));
395 
396  }
397 
398 
399 
400  } // End loop over the primary photons
401 
402 
403  } // Event with one or two photons
404 
405 
406  //std::cout << " PhotonMCTruthFinder photon size " << result.size() << std::endl;
407  }
408 
409  return result;
410 
411 }
412 
413 void PhotonMCTruthFinder::fill(std::vector<SimTrack>& simTracks, std::vector<SimVertex>& simVertices ) {
414  // std::cout << " PhotonMCTruthFinder::fill " << std::endl;
415 
416 
417 
418  // Watch out there ! A SimVertex is in mm (stupid),
419 
420  unsigned nVtx = simVertices.size();
421  unsigned nTks = simTracks.size();
422 
423  // std::cout << " PhotonMCTruthFinder::fill " << nVtx << " " << nTks << std::endl;
424 
425  // Empty event, do nothin'
426  if ( nVtx == 0 ) return;
427 
428  // create a map associating geant particle id and position in the
429  // event SimTrack vector
430  for( unsigned it=0; it<nTks; ++it ) {
431  geantToIndex_[ simTracks[it].trackId() ] = it;
432  // std::cout << " PhotonMCTruthFinder::fill it " << it << " simTracks[it].trackId() " << simTracks[it].trackId() << std::endl;
433 
434  }
435 
436 
437 
438 
439 }
math::XYZTLorentzVectorD
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > XYZTLorentzVectorD
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:15
abs
#define abs(x)
Definition: mlp_lapack.h:159
detailsBasic3DVector::z
double double double z
Definition: newBasic3DVector.h:17
PhotonMCTruthFinder::fill
void fill(std::vector< SimTrack > &theSimTracks, std::vector< SimVertex > &theSimVertices)
Definition: PhotonMCTruthFinder.cc:413
SimVertex::parentIndex
int parentIndex() const
Definition: SimVertex.h:33
PhotonMCTruthFinder::geantToIndex_
std::map< unsigned, unsigned > geantToIndex_
Definition: PhotonMCTruthFinder.h:40
query.result
tuple result
Definition: query.py:137
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(std::vector< SimTrack > simTracks, std::vector< SimVertex > simVertices)
Definition: PhotonMCTruthFinder.cc:20
CoreSimTrack::type
int type() const
particle type (HEP PDT convension)
Definition: CoreSimTrack.h:40
CoreSimTrack::momentum
const math::XYZTLorentzVectorD & momentum() const
particle info...
Definition: CoreSimTrack.h:36
if
perl if(1 lt scalar(@::datatypes))
Definition: edlooper.cc:31
x
Definition: DDAxes.h:10
Generated for CMSSW Reference Manual by doxygen 1.8.5