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