CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
PhotosInterface.cc
Go to the documentation of this file.
1 
2 #include <iostream>
3 
4 // #include "GeneratorInterface/Pythia6Interface/interface/Pythia6Service.h"
5 
7 
9 
10 #include "HepMC/GenEvent.h"
11 #include "HepMC/IO_HEPEVT.h"
12 #include "HepMC/HEPEVT_Wrapper.h"
13 
14 using namespace gen;
15 using namespace edm;
16 using namespace std;
17 
18 
19 extern "C"{
20 
21  void phoini_( void );
22  void photos_( int& );
23 
24  double phoran_(int *idummy)
25  {
26  return decayRandomEngine->flat();
27  }
28 /*
29  double phoranc_(int *idummy)
30  {
31  return decayRandomEngine->flat();
32  }
33 */
34 
35  extern struct {
36  // bool qedrad[NMXHEP];
37  bool qedrad[4000]; // hardcoded for now...
38  } phoqed_;
39 
40 }
41 
42 
44  : fOnlyPDG(-1)
45 {
46  fSpecialSettings.push_back("QED-brem-off:all");
48  fIsInitialized = false;
49 }
50 
52  : fOnlyPDG(-1)
53 {
54  fSpecialSettings.push_back("QED-brem-off:all");
55  fIsInitialized = false;
56 }
57 
59 {
60 
61  fOnlyPDG = ipdg;
62 // std::ostringstream command;
63 // command << "QED-brem-off:" << fOnlyPDG ;
64  fSpecialSettings.clear();
65 // fSpecialSettings.push_back( command.str() );
66 
67  return;
68 
69 }
70 
72 {
73 
74  if ( fIsInitialized ) return; // do init only once
75 
76  phoini_();
77 
78  fIsInitialized = true;
79 
80  return;
81 }
82 
83 HepMC::GenEvent* PhotosInterface::apply( HepMC::GenEvent* evt )
84 {
85 
86  if ( !fIsInitialized ) return evt; // conv.read_next_event();
87 
88  // loop over HepMC::GenEvent, find vertices
89 
90  // for ( int ip=0; ip<evt->particles_size(); ip++ )
91  for ( int ip=0; ip<4000; ip++ ) // 4000 is the max size of the array
92  {
93  phoqed_.qedrad[ip]=true;
94  }
95 
96  //
97  // now do actual job
98  //
99 
100  for ( int iv=1; iv<=evt->vertices_size(); iv++ )
101  {
102 
103  bool legalVtx = false;
104 
105  fSecVtxStore.clear();
106 
107  HepMC::GenVertex* vtx = evt->barcode_to_vertex( -iv ) ;
108 
109  if ( vtx->particles_in_size() != 1 ) continue; // more complex than we need
110  if ( vtx->particles_out_size() <= 1 ) continue; // no outcoming particles
111 
112  if ( (*(vtx->particles_in_const_begin()))->pdg_id() == 111 ) continue; // pi0 decay vtx - no point to try
113 
114  if ( fOnlyPDG != 1 && (*(vtx->particles_in_const_begin()))->pdg_id() != fOnlyPDG )
115  {
116  continue;
117  }
118  else
119  {
120  // requested for specific PDG ID only, typically tau (15)
121  //
122  // first check if a brem vertex, where outcoming are the same pdg id and a photon
123  //
124  bool same = false;
125  for ( HepMC::GenVertex::particle_iterator pitr=vtx->particles_begin(HepMC::children);
126  pitr != vtx->particles_end(HepMC::children); ++pitr)
127  {
128  if ( (*pitr)->pdg_id() == fOnlyPDG )
129  {
130  same = true;
131  break;
132  }
133  }
134  if ( same ) continue;
135 
136  // OK, we get here if incoming fOnlyPDG and something else outcoming
137  // call it for the whole branch starting at vtx barcode iv, and go on
138  // NOTE: theoretically, it has a danger of double counting in vertices
139  // down the decay branch originating from fOnlyPDG, but in practice
140  // it's unlikely that down the branchg there'll be more fOnlyPDG's
141 
142  // cross-check printout
143  // vtx->print();
144 
145  // overprotection...
146  //
147  if ( fOnlyPDG == 15 && fAvoidTauLeptonicDecays && isTauLeptonicDecay( vtx ) ) continue;
148 
149  applyToBranch( evt, -iv );
150  continue;
151  }
152 
153  // configured for all types of particles
154  //
155  for ( HepMC::GenVertex::particle_iterator pitr=vtx->particles_begin(HepMC::children);
156  pitr != vtx->particles_end(HepMC::children); ++pitr)
157  {
158 
159  // quark or gluon out of this vertex - no good !
160  if ( abs((*pitr)->pdg_id()) >=1 && abs((*pitr)->pdg_id()) <=8 ) break;
161  if ( abs((*pitr)->pdg_id()) == 21 ) break;
162 
163  if ( (*pitr)->status() == 1 || (*pitr)->end_vertex() )
164  {
165  // OK, legal already !
166  legalVtx = true;
167  break;
168  }
169  }
170 
171  if ( !legalVtx ) continue;
172 
173  applyToVertex( evt, -iv );
174 
175  } // end of master loop
176 
177  // restore event number in HEPEVT (safety measure, somehow needed by Hw6)
178  HepMC::HEPEVT_Wrapper::set_event_number( evt->event_number() );
179 
180  return evt;
181 
182 }
183 
184 void PhotosInterface::applyToVertex( HepMC::GenEvent* evt, int vtxbcode )
185 {
186 
187  HepMC::GenVertex* vtx = evt->barcode_to_vertex( vtxbcode );
188 
189  if ( fAvoidTauLeptonicDecays && isTauLeptonicDecay( vtx ) ) return;
190 
191  // cross-check printout
192  //
193  // vtx->print();
194 
195  // first, flush out HEPEVT & tmp barcode storage
196  //
197  HepMC::HEPEVT_Wrapper::zero_everything();
198  fBarcodes.clear();
199 
200  // add incoming particle
201  //
202  int index = 1;
203  HepMC::HEPEVT_Wrapper::set_id( index, (*(vtx->particles_in_const_begin()))->pdg_id() );
204  HepMC::FourVector vec4;
205  vec4 = (*(vtx->particles_in_const_begin()))->momentum();
206  HepMC::HEPEVT_Wrapper::set_momentum( index, vec4.x(), vec4.y(), vec4.z(), vec4.e() );
207  HepMC::HEPEVT_Wrapper::set_mass( index, (*(vtx->particles_in_const_begin()))->generated_mass() );
208  HepMC::HEPEVT_Wrapper::set_position( index, vtx->position().x(), vtx->position().y(),
209  vtx->position().z(), vtx->position().t() );
210  HepMC::HEPEVT_Wrapper::set_status( index, (*(vtx->particles_in_const_begin()))->status() );
211  HepMC::HEPEVT_Wrapper::set_parents( index, 0, 0 );
212  fBarcodes.push_back( (*(vtx->particles_in_const_begin()))->barcode() );
213 
214  // add outcoming particles (decay products)
215  //
216  int lastDau = 1;
217  for ( HepMC::GenVertex::particle_iterator pitr=vtx->particles_begin(HepMC::children);
218  pitr != vtx->particles_end(HepMC::children); ++pitr)
219  {
220 
221  if ( (*pitr)->status() == 1 || (*pitr)->end_vertex() )
222  {
223  index++;
224  vec4 = (*pitr)->momentum();
225  HepMC::HEPEVT_Wrapper::set_id( index, (*pitr)->pdg_id() );
226  HepMC::HEPEVT_Wrapper::set_momentum( index, vec4.x(), vec4.y(), vec4.z(), vec4.e() );
227  HepMC::HEPEVT_Wrapper::set_mass( index, (*pitr)->generated_mass() );
228  vec4 = (*pitr)->production_vertex()->position();
229  HepMC::HEPEVT_Wrapper::set_position( index, vec4.x(), vec4.y(), vec4.z(), vec4.t() );
230  HepMC::HEPEVT_Wrapper::set_status( index, (*pitr)->status() );
231  HepMC::HEPEVT_Wrapper::set_parents( index, 1, 1 );
232  fBarcodes.push_back( (*pitr)->barcode() );
233  lastDau++;
234  }
235  if ( (*pitr)->end_vertex() )
236  {
237  fSecVtxStore.push_back( (*pitr)->end_vertex()->barcode() );
238  }
239  }
240 
241  // store, further to set NHEP in HEPEVT
242  //
243  int nentries = index;
244 
245  // reset master pointer to mother
246  index = 1;
247  HepMC::HEPEVT_Wrapper::set_children ( index, 2, lastDau ); // FIXME: need to check
248  // if last daughter>=2 !!!
249 
250  // finally, set number of entries (NHEP) in HEPEVT
251  //
252  HepMC::HEPEVT_Wrapper::set_number_entries( nentries );
253 
254  // cross-check printout HEPEVT
255  // HepMC::HEPEVT_Wrapper::print_hepevt();
256 
257  // OK, 1-level vertex is formed - now, call PHOTOS
258  //
259  photos_( index ) ;
260 
261  // another cross-check printout HEPEVT - after photos
262  // HepMC::HEPEVT_Wrapper::print_hepevt();
263 
264 
265  // now check if something has been generated
266  // and make all adjustments to underlying vtx/parts
267  //
268  attachParticles( evt, vtx, nentries );
269 
270  // ugh, done with this vertex !
271 
272  return;
273 
274 }
275 
276 void PhotosInterface::applyToBranch( HepMC::GenEvent* evt, int vtxbcode )
277 {
278 
279 
280  fSecVtxStore.clear();
281 
282  // 1st level vertex
283  //
284  applyToVertex( evt, vtxbcode );
285 
286  // now look down the branch for more vertices, if any
287  //
288  // Note: fSecVtxStore gets filled up in applyToVertex, if necessary
289  //
290  unsigned int vcounter = 0;
291 
292  while ( vcounter < fSecVtxStore.size() )
293  {
294  applyToVertex( evt, fSecVtxStore[vcounter] );
295  vcounter++;
296  }
297 
298  fSecVtxStore.clear();
299 
300  return;
301 
302 }
303 
304 void PhotosInterface::attachParticles( HepMC::GenEvent* evt, HepMC::GenVertex* vtx, int nentries )
305 {
306 
307  if ( HepMC::HEPEVT_Wrapper::number_entries() > nentries )
308  {
309  // yes, need all corrections and adjustments -
310  // figure out how many photons and what particles in
311  // the decay branch have changes;
312  // also, follow up each one and correct accordingly;
313  // at the same time, add photon(s) to the GenVertex
314  //
315 
316  // vtx->print();
317 
318  int largestBarcode = -1;
319  int Nbcodes = fBarcodes.size();
320 
321  for ( int ip=1; ip<Nbcodes; ip++ )
322  {
323 
324  int bcode = fBarcodes[ip];
325  HepMC::GenParticle* prt = evt->barcode_to_particle( bcode );
326  if ( bcode > largestBarcode ) largestBarcode = bcode;
327  double px = HepMC::HEPEVT_Wrapper::px(ip+1);
328  double py = HepMC::HEPEVT_Wrapper::py(ip+1);
329  double pz = HepMC::HEPEVT_Wrapper::pz(ip+1);
330  double e = HepMC::HEPEVT_Wrapper::e(ip+1);
331  double m = HepMC::HEPEVT_Wrapper::m(ip+1);
332 
333  if ( prt->end_vertex() )
334  {
335 
336  HepMC::GenVertex* endVtx = prt->end_vertex();
337 
338  std::vector<int> secVtxStorage;
339  secVtxStorage.clear();
340 
341  secVtxStorage.push_back( endVtx->barcode() );
342 
343  HepMC::FourVector mom4 = prt->momentum();
344 
345  // now rescale all descendants
346  double bet1[3], bet2[3], gam1, gam2, pb;
347  double mass = mom4.m();
348  bet1[0] = -(mom4.px()/mass);
349  bet1[1] = -(mom4.py()/mass);
350  bet1[2] = -(mom4.pz()/mass);
351  bet2[0] = px/m;
352  bet2[1] = py/m;
353  bet2[2] = pz/m;
354  gam1 = mom4.e()/mass;
355  gam2 = e/m;
356 
357  unsigned int vcounter = 0;
358 
359  while ( vcounter < secVtxStorage.size() )
360  {
361 
362  HepMC::GenVertex* theVtx = evt->barcode_to_vertex( secVtxStorage[vcounter] );
363 
364  for ( HepMC::GenVertex::particle_iterator pitr=theVtx->particles_begin(HepMC::children);
365  pitr != theVtx->particles_end(HepMC::children); ++pitr)
366  {
367 
368  if ( (*pitr)->end_vertex() )
369  {
370  secVtxStorage.push_back( (*pitr)->end_vertex()->barcode() );
371  }
372 
373  if ( theVtx->particles_out_size() == 1 && (*pitr)->pdg_id() == prt->pdg_id() )
374  {
375  // carbon copy
376  (*pitr)->set_momentum( HepMC::FourVector(px,py,pz,e) );
377  continue;
378  }
379 
380  HepMC::FourVector dmom4 = (*pitr)->momentum();
381 
382  // Boost vector to parent rest frame...
383  pb = bet1[0]*dmom4.px() + bet1[1]*dmom4.py() + bet1[2]*dmom4.pz();
384  double dpx = dmom4.px() + bet1[0] * (dmom4.e() + pb/(gam1+1.) );
385  double dpy = dmom4.py() + bet1[1] * (dmom4.e() + pb/(gam1+1.) );
386  double dpz = dmom4.pz() + bet1[2] * (dmom4.e() + pb/(gam1+1.) );
387  double de = gam1*dmom4.e() + pb;
388  // ...and boost back to modified parent frame
389  pb = bet2[0]*dpx + bet2[1]*dpy + bet2[2]*dpz;
390  dpx += bet2[0] * ( de + pb/(gam2+1.) );
391  dpy += bet2[1] * ( de + pb/(gam2+1.) );
392  dpz += bet2[2] * ( de + pb/(gam2+1.) );
393  de *= gam2;
394  de += pb;
395 
396  (*pitr)->set_momentum( HepMC::FourVector(dpx,dpy,dpz,de) );
397 
398  }
399  vcounter++;
400  }
401 
402  secVtxStorage.clear();
403  }
404 
405  prt->set_momentum( HepMC::FourVector(px,py,pz,e) );
406 
407  } // ok, all affected particles update, but the photon(s) still not inserted
408 
409  int newlyGen = HepMC::HEPEVT_Wrapper::number_entries() - nentries;
410 
411  if ( largestBarcode < evt->particles_size() )
412  {
413  // need to adjust barcodes down from the affected vertex/particles
414  // such that we "free up" barcodes for newly generated photons
415  // in the middle of the event record
416  //
417  for ( int ipp=evt->particles_size(); ipp>largestBarcode; ipp-- )
418  {
419  (evt->barcode_to_particle(ipp))->suggest_barcode( ipp+newlyGen );
420  }
421  }
422 
423  // now attach new generated photons to the vertex
424  //
425  for ( int ipnw=1; ipnw<=newlyGen; ipnw++ )
426  {
427  int nbcode = largestBarcode+ipnw;
428  int pdg_id = HepMC::HEPEVT_Wrapper::id( nentries+ipnw );
429  int status = HepMC::HEPEVT_Wrapper::status( nentries+ipnw );
430  double px = HepMC::HEPEVT_Wrapper::px( nentries+ipnw );
431  double py = HepMC::HEPEVT_Wrapper::py( nentries+ipnw );
432  double pz = HepMC::HEPEVT_Wrapper::pz( nentries+ipnw );
433  double e = HepMC::HEPEVT_Wrapper::e( nentries+ipnw );
434  double m = HepMC::HEPEVT_Wrapper::m( nentries+ipnw );
435 
436  HepMC::GenParticle* NewPart = new HepMC::GenParticle( HepMC::FourVector(px,py,pz,e),
437  pdg_id, status);
438  NewPart->set_generated_mass( m );
439  NewPart->suggest_barcode( nbcode );
440  vtx->add_particle_out( NewPart ) ;
441  }
442 
443  //vtx->print();
444  //std::cout << " leaving attachParticles() " << std::endl;
445 
446  } // end of global if-statement
447 
448  return;
449 }
450 
451 bool PhotosInterface::isTauLeptonicDecay( HepMC::GenVertex* vtx )
452 {
453 
454  if ( abs((*(vtx->particles_in_const_begin()))->pdg_id()) != 15 ) return false;
455 
456  for ( HepMC::GenVertex::particle_iterator pitr=vtx->particles_begin(HepMC::children);
457  pitr != vtx->particles_end(HepMC::children); ++pitr)
458  {
459  if ( abs((*pitr)->pdg_id()) == 11 || abs((*pitr)->pdg_id()) == 13 )
460  {
461  return true;
462  }
463  }
464 
465  return false;
466 
467 }
468 
469 /* very first version... but Phptos seems to want a SINGLE vertex, nor a branch
470 
471 void PhotosInterface::applyToBranch( HepMC::GenEvent* evt, int vtxbcode )
472 {
473 
474  HepMC::GenVertex* vtx = evt->barcode_to_vertex( vtxbcode );
475 
476  vtx->print();
477 
478  // special case - do nothing
479  // we don't brem off tau since it'll be done by master generator,
480  // and within tau leptonic decays the brem is done by tauola
481  //
482  if ( fAvoidTauLeptonicDecays && isTauLeptonicDecay( vtx ) ) return;
483 
484  std::vector<int> secVtxStore;
485  secVtxStore.clear();
486 
487  // first, flush out HEPEVT & tmp barcode storage
488  //
489  HepMC::HEPEVT_Wrapper::zero_everything();
490  fBarcodes.clear();
491 
492  // form 1st level vertex
493  //
494  // add incoming particle
495  //
496  int index = 1;
497  HepMC::HEPEVT_Wrapper::set_id( index, (*(vtx->particles_in_const_begin()))->pdg_id() );
498  HepMC::FourVector vec4;
499  vec4 = (*(vtx->particles_in_const_begin()))->momentum();
500  HepMC::HEPEVT_Wrapper::set_momentum( index, vec4.x(), vec4.y(), vec4.z(), vec4.e() );
501  HepMC::HEPEVT_Wrapper::set_mass( index, (*(vtx->particles_in_const_begin()))->generated_mass() );
502  HepMC::HEPEVT_Wrapper::set_position( index, vtx->position().x(), vtx->position().y(),
503  vtx->position().z(), vtx->position().t() );
504  HepMC::HEPEVT_Wrapper::set_status( index, (*(vtx->particles_in_const_begin()))->status() );
505  HepMC::HEPEVT_Wrapper::set_parents( index, 0, 0 );
506  fBarcodes.push_back( (*(vtx->particles_in_const_begin()))->barcode() );
507 
508  int lastDau = 1;
509  for ( HepMC::GenVertex::particle_iterator pitr=vtx->particles_begin(HepMC::children);
510  pitr != vtx->particles_end(HepMC::children); ++pitr)
511  {
512 
513  // put particles into HEPEVT - form 1st level vertex
514  //
515  if ( (*pitr)->status() == 1 || (*pitr)->end_vertex() )
516  {
517  index++;
518  vec4 = (*pitr)->momentum();
519  HepMC::HEPEVT_Wrapper::set_id( index, (*pitr)->pdg_id() );
520  HepMC::HEPEVT_Wrapper::set_momentum( index, vec4.x(), vec4.y(), vec4.z(), vec4.e() );
521  HepMC::HEPEVT_Wrapper::set_mass( index, (*pitr)->generated_mass() );
522  vec4 = (*pitr)->production_vertex()->position();
523  HepMC::HEPEVT_Wrapper::set_position( index, vec4.x(), vec4.y(), vec4.z(), vec4.t() );
524  HepMC::HEPEVT_Wrapper::set_status( index, (*pitr)->status() );
525  HepMC::HEPEVT_Wrapper::set_parents( index, 1, 1 );
526  fBarcodes.push_back( (*pitr)->barcode() );
527  lastDau++;
528  }
529 
530  if ( (*pitr)->end_vertex() )
531  {
532  secVtxStore.push_back( (*pitr)->end_vertex()->barcode() );
533  }
534  }
535 
536  if ( lastDau < 2 ) lastDau = 2;
537  HepMC::HEPEVT_Wrapper::set_children ( 1, 2, lastDau );
538 
539  // now look down the branch for more vertices, if any
540  //
541  unsigned int vcounter = 0;
542  int firstDau = lastDau + 1;
543  int index1 = index;
544  while ( vcounter < secVtxStore.size() )
545  {
546  HepMC::GenVertex* theVtx = evt->barcode_to_vertex( secVtxStore[vcounter] );
547  for ( HepMC::GenVertex::particle_iterator pitr1=theVtx->particles_begin(HepMC::children);
548  pitr1 != theVtx->particles_end(HepMC::children); ++pitr1)
549  {
550  if ( (*pitr1)->status() == 1 || (*pitr1)->end_vertex() )
551  {
552  index++;
553  vec4 = (*pitr1)->momentum();
554  HepMC::HEPEVT_Wrapper::set_id( index, (*pitr1)->pdg_id() );
555  HepMC::HEPEVT_Wrapper::set_momentum( index, vec4.x(), vec4.y(), vec4.z(), vec4.e() );
556  HepMC::HEPEVT_Wrapper::set_mass( index, (*pitr1)->generated_mass() );
557  vec4 = (*pitr1)->production_vertex()->position();
558  HepMC::HEPEVT_Wrapper::set_position( index, vec4.x(), vec4.y(), vec4.z(), vec4.t() );
559  HepMC::HEPEVT_Wrapper::set_status( index, (*pitr1)->status() );
560  HepMC::HEPEVT_Wrapper::set_parents( index, index1, index1 );
561  fBarcodes.push_back( (*pitr1)->barcode() );
562  lastDau++;
563  }
564  if ( (*pitr1)->end_vertex() )
565  {
566  secVtxStore.push_back( (*pitr1)->end_vertex()->barcode() );
567  }
568  }
569  index1 += 1;
570  HepMC::HEPEVT_Wrapper::set_children ( index1, firstDau, lastDau );
571  index1 = index;
572  firstDau = lastDau + 1;
573  vcounter++;
574  }
575 
576  // finally, set number of entries (NHEP) in HEPEVT
577  //
578  int nentries = index;
579  HepMC::HEPEVT_Wrapper::set_number_entries( nentries );
580 
581  // test printout
582  HepMC::HEPEVT_Wrapper::print_hepevt();
583 
584  // don't start from the first one since it's going to be
585  // the incoming particles (e.g. tau) which is already
586  // treated by the brem from master generator
587  //
588  index = 2;
589  // oh well, maybe 2 was a bad idea...
590  index = 1;
591  photos_( index );
592 
593  HepMC::HEPEVT_Wrapper::print_hepevt();
594 
595  attachParticles( evt, vtx, nentries );
596 
597  return;
598 
599 }
600 */
601 
void phoini_(void)
void attachParticles(HepMC::GenEvent *, HepMC::GenVertex *, int)
double phoran_(int *idummy)
void photos_(int &)
std::vector< int > fBarcodes
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
std::vector< vec3 > vec4
Definition: HCALResponse.h:18
void applyToBranch(HepMC::GenEvent *, int)
std::vector< std::string > fSpecialSettings
bool qedrad[4000]
struct @360 phoqed_
std::vector< int > fSecVtxStore
void applyToVertex(HepMC::GenEvent *, int)
CLHEP::HepRandomEngine * decayRandomEngine
tuple status
Definition: ntuplemaker.py:245
bool isTauLeptonicDecay(HepMC::GenVertex *)
HepMC::GenEvent * apply(HepMC::GenEvent *)