Photos.cxx
1 #include <stdarg.h>
2 #include <iostream>
3 #include <vector>
4 
5 #include "PhotosRandom.h"
6 #include "PhotosEvent.h"
7 #include "Photos.h"
8 #include "Log.h"
9 using std::vector;
10 using std::cout;
11 using std::endl;
12 using std::ios_base;
13 
14 namespace Photospp
15 {
16 
17 Photos Photos::_instance;
18 
19 vector<vector<int>* > *Photos::supBremList = 0;
20 vector<vector<int>* > *Photos::forceBremList = 0;
21 vector<pair<int,double>* > *Photos::forceMassList = 0;
22 vector<int> *Photos::ignoreStatusCodeList = 0;
23 bool Photos::isSuppressed=false;
24 bool Photos::massFrom4Vector=true;
31 bool Photos::IfPair=false;
32 bool Photos::IfPhot=true;
33 int Photos::EventNo=0;
34 double (*Photos::randomDouble)() = PhotosRandom::randomReal;
35 
36 Photos::MomentumUnits Photos::momentumUnit = Photos::DEFAULT_MOMENTUM;
37 
38 Photos::Photos()
39 {
40  setAlphaQED (0.00729735039);
41  setInfraredCutOff (0.01);
42  setInterference (true);
43  setDoubleBrem (true);
44  setQuatroBrem (false);
46  setCorrectionWtForW (true);
47 
48  // setExponentiation(true) moved to initialize() due to communication
49  // problems with Fortran under MacOS.
50  phokey.iexp = 1;
51 }
52 
54 {
55 // Should return if already initialized?
56 
57 /*******************************************************************************
58  All the following parameter setters can be called after PHOINI.
59  Initialization of kinematic correction against rounding errors.
60  The set values will be used later if called with zero.
61  Default parameter is 1 (no correction) optionally 2, 3, 4
62  In case of exponentiation new version 5 is needed in most cases.
63  Definition given here will be thus overwritten in such a case
64  below in routine PHOCIN
65 *******************************************************************************/
66  if(!phokey.iexp) initializeKinematicCorrections(1);
67  else setExponentiation(true);
68 
69 // Initialize status counter for warning messages
70  for(int i=0;i<10;i++) phosta.status[i]=0;
71 // elementary security level, should remain 1 but we may want to have a method to change.
72  phosta.ifstop=1;
73 
74  pholun.phlun=6; // Logical output unit for printing error messages
75 
76 // Further initialization done automatically
77 // see places with - VARIANT A - VARIANT B - all over to switch between options
78 
79 #ifndef VARIANTB
80 //----------- SLOWER VARIANT A, but stable ------------------
81 //--- it is limiting choice for small XPHCUT in fixed orer
82 //--- modes of operation
83 
84 // Best choice is if FINT=2**N where N+1 is maximal number
85 // of charged daughters
86 // see report on overweihted events
87  if(phokey.interf) maxWtInterference(2.0);
88  else maxWtInterference(1.0);
89 #else
90 
91 //----------- FASTER VARIANT B ------------------
92 //-- it is good for tests of fixed order and small XPHCUT
93 //-- but is less promising for more complex cases of interference
94 //-- sometimes fails because of that
95 
96  if(phokey.interf) maxWtInterference(1.8);
97  else maxWtInterference(0.0);
98 #endif
99 //------------END VARIANTS A B -----------------------
100 
101 //------------------------------------------------------------------------------
102 // Print PHOTOS header
103 //------------------------------------------------------------------------------
104  int coutPrec = cout.precision(6);
105  ios_base::fmtflags flags = cout.setf(ios_base::scientific, ios_base::floatfield);
106  cout<<endl;
107  cout<<"********************************************************************************"<<endl<<endl;
108 
109  cout<<" ========================="<<endl;
110  cout<<" PHOTOS, Version: "<<VER_MAJOR<<"."<<VER_MINOR<<endl;
111  cout<<" Released at: "<<DAT_DAY<<"/"<<DAT_MONTH<<"/"<<DAT_YEAR<<endl;
112  cout<<" ========================="<<endl<<endl;
113 
114  cout<<" Photos QED corrections in Particle Decays"<<endl<<endl;
115 
116  cout<<" Monte Carlo Program - by E. Barberio, B. van Eijk and Z. Was"<<endl;
117  cout<<" From version 2.09 - by P. Golonka and Z. Was"<<endl;
118  cout<<" From version 3.00 - by N. Davidson, T. Przedzinski and Z. Was"<<endl;
119 
120  cout<<"********************************************************************************"<<endl<<endl;
121 
122  cout<<" Internal (default) input parameters: "<<endl<<endl;
123  cout<<" INTERF= "<<phokey.interf<<" ISEC= " <<phokey.isec <<" ITRE= "<<phokey.itre
124  <<" IEXP= " <<phokey.iexp <<" IFTOP= "<<phokey.iftop<<" IFW= " <<phokey.ifw <<endl;
125  cout<<" ALPHA_QED= "<<phocop.alpha<<" XPHCUT= "<<phocop.xphcut<<endl<<endl;
126 
127  if(phokey.interf) cout<<" Option with interference is active"<<endl;
128  if(phokey.isec) cout<<" Option with double photons is active"<<endl;
129  if(phokey.itre) cout<<" Option with triple/quatric photons is active"<<endl;
130  if(phokey.iexp) cout<<" Option with exponentiation is active EPSEXP="<<phokey.expeps<<endl;
131  if(phokey.iftop) cout<<" Emision in t tbar production is active"<<endl;
132  if(phokey.ifw) cout<<" Correction wt in decay of W is active"<<endl;
133  if(meCorrectionWtForZ) cout<<" ME in decay of Z is active"<<endl;
134  if(meCorrectionWtForW) cout<<" ME in decay of W is active"<<endl;
135  if(IfPair) cout<<" emission of pairs is active"<<endl;
136  if(!IfPhot) cout<<" emission of photons is inactive"<<endl;
137 
138  cout<<endl<<" WARNING: /HEPEVT/ is not anymore used."<<endl<<endl;
139 /*
140  cout<<endl<<" WARNING (1): /HEPEVT/ is not anymore the standard common block"<<endl<<endl;
141 
142  cout<<" PHOTOS expects /HEPEVT/ to have REAL*8 variables. To change to"<<endl;
143  cout<<" REAL*4 modify its declaration in subr. PHOTOS_GET PHOTOS_SET:"<<endl;
144  cout<<" REAL*8 d_h_phep, d_h_vhep"<<endl;
145  cout<<" WARNING (2): check dims. of /hepevt/ /phoqed/ /ph_hepevt/."<<endl;
146  cout<<" HERE: d_h_nmxhep=10000 and NMXHEP=10000"<<endl<<endl;
147 */
148  cout<<"********************************************************************************"<<endl;
149  // Revert output stream flags and precision
150  cout.precision(coutPrec);
151  cout.flags (flags);
152 
153  // Add reggeon, pomeron and its diffractive states to the list
154  // of decays where bremsstrahlung is suppressed
157  Photos::suppressBremForDecay(0,9902110);
158  Photos::suppressBremForDecay(0,9902210);
159  Photos::suppressBremForDecay(0,9900110);
160  Photos::suppressBremForDecay(0,9900210);
161  Photos::suppressBremForDecay(0,9900220);
162  Photos::suppressBremForDecay(0,9900330);
163  Photos::suppressBremForDecay(0,9900440);
164 
165  // Set suppression of all pi0 decays and K_L -> gamma e+ e- ...
166  // Previously initialization in Fortran IPHEKL(i) routine and used in PHOCHK
167  // i=1 was emission allowed, i=2 was blocked 0 was when the option was used.
168  // now in IPHEKL_setPi0KLnoEmmision we have only 1 to allow emissions
169  // and 2 to block.
170  // Can be overriden by using 'Photos::IPHEKL_setPi0KLnoEmmision(0)'
171  // method several times use Photos::forceBremForDecay() and can be
172  // over-ruled in part.
173 
174  Photos::IPHEKL_setPi0KLnoEmission(2); // Blocks emission in pi0 and in Kl to gamma e+ e- if parameter is !1 (enables otherwise)
175  Photos::IPHQRK_setQarknoEmission (1,0);// Blocks emission from quarks if buf=1 (default); enables if buf=2
176  // option 2 is not realistic and for tests only
177 
178 // Initialize Marsaglia and Zaman random number generator
179  PhotosRandom::initialize();
180 }
182 {
183 // prints infomation like Photos::initialize; may be called after reinitializations.
184 
185 /*******************************************************************************
186  Once parameter setters are called after PHOINI one may want to know and write
187  into output info including all reinitializations.
188 *******************************************************************************/
189 
190 
191 //------------------------------------------------------------------------------
192 // Print PHOTOS header again
193 //------------------------------------------------------------------------------
194  int coutPrec = cout.precision(6);
195  ios_base::fmtflags flags = cout.setf(ios_base::scientific, ios_base::floatfield);
196  cout<<endl;
197  cout<<"********************************************************************************"<<endl<<endl;
198  cout<<" ========================================="<<endl;
199  cout<<" PHOTOS, information routine"<<endl;
200  cout<<" Input parameters after reinitialization: "<<endl<<endl;
201  cout<<" ========================================="<<endl<<endl;
202  cout<<"********************************************************************************"<<endl<<endl;
203  cout<<" INTERF= "<<phokey.interf<<" ISEC= " <<phokey.isec <<" ITRE= "<<phokey.itre
204  <<" IEXP= " <<phokey.iexp <<" IFTOP= "<<phokey.iftop<<" IFW= " <<phokey.ifw <<endl;
205  cout<<" ALPHA_QED= "<<phocop.alpha<<" XPHCUT= "<<phocop.xphcut<<endl<<endl;
206 
207  if(phokey.interf) cout<<" Option with interference is active"<<endl;
208  if(phokey.isec) cout<<" Option with double photons is active"<<endl;
209  if(phokey.itre) cout<<" Option with triple/quatric photons is active"<<endl;
210  if(phokey.iexp) cout<<" Option with exponentiation is active EPSEXP="<<phokey.expeps<<endl;
211  if(phokey.iftop) cout<<" Emision in t tbar production is active"<<endl;
212  if(phokey.ifw) cout<<" Correction wt in decay of W is active"<<endl;
213  if(meCorrectionWtForZ) cout<<" ME in decay of Z is active"<<endl;
214  if(meCorrectionWtForW) cout<<" ME in decay of W is active"<<endl;
215  if(meCorrectionWtForScalar) cout<<" ME in decay of Scalar is active"<<endl;
216  if(IfPair) cout<<" emission of pairs is active"<<endl;
217  if(!IfPhot) cout<<" emission of photons is inactive"<<endl;
218 
219  cout<<endl<<" WARNING: /HEPEVT/ is not anymore used."<<endl<<endl;
220  // Revert output stream flags and precision
221  cout.precision(coutPrec);
222  cout.flags (flags);
223 }
224 
226 {
227  PhotosBranch b(p);
228  if(!b.getSuppressionStatus()) b.process();
229 }
230 
232 {
233  vector<PhotosParticle *> particles = p->getDecayTree();
234  vector<PhotosBranch *> branches = PhotosBranch::createBranches(particles);
235  for(int i=0;i<(int)branches.size();i++) branches.at(i)->process();
236 }
237 
238 void Photos::suppressBremForDecay(int count, int motherID, ... )
239 {
240  va_list arg;
241  va_start(arg, motherID);
242  vector<int> *v = new vector<int>();
243  v->push_back(motherID);
244  for(int i = 0;i<count;i++)
245  {
246  v->push_back(va_arg(arg,int));
247  }
248  va_end(arg);
249  v->push_back(0);
250  if(!supBremList) supBremList = new vector< vector<int>* >();
251  supBremList->push_back(v);
252 }
253 
254 void Photos::suppressBremForBranch(int count, int motherID, ... )
255 {
256  va_list arg;
257  va_start(arg, motherID);
258  vector<int> *v = new vector<int>();
259  v->push_back(motherID);
260  for(int i = 0;i<count;i++)
261  {
262  v->push_back(va_arg(arg,int));
263  }
264  va_end(arg);
265  v->push_back(1);
266  if(!supBremList) supBremList = new vector< vector<int>* >();
267  supBremList->push_back(v);
268 }
269 
270 void Photos::forceBremForDecay(int count, int motherID, ... )
271 {
272  va_list arg;
273  va_start(arg, motherID);
274  vector<int> *v = new vector<int>();
275  v->push_back(motherID);
276  for(int i = 0;i<count;i++)
277  {
278  v->push_back(va_arg(arg,int));
279  }
280  va_end(arg);
281  v->push_back(0);
282  if(!forceBremList) forceBremList = new vector< vector<int>* >();
283  forceBremList->push_back(v);
284 }
285 
286 void Photos::forceBremForBranch(int count, int motherID, ... )
287 {
288  va_list arg;
289  va_start(arg, motherID);
290  vector<int> *v = new vector<int>();
291  v->push_back(motherID);
292  for(int i = 0;i<count;i++)
293  {
294  v->push_back(va_arg(arg,int));
295  }
296  va_end(arg);
297  v->push_back(1);
298  if(!forceBremList) forceBremList = new vector< vector<int>* >();
299  forceBremList->push_back(v);
300 }
301 
302  // Previously this functionality was encoded in FORTRAN routine
303  // PHOCHK which was having some other functionality as well
305 {
306  if(m==1)
307  {
308  cout << "MODOP=1 -- enables emission in pi0 to gamma e+e- : TEST " << endl ;
309  cout << "MODOP=1 -- enables emission in Kl to gamma e+e- : TEST " << endl ;
311  Photos::forceBremForDecay(3, 130,22,11,-11);
312  Photos::forceBremForDecay(3,-130,22,11,-11);
313  }
314  else if(m!=1)
315  {
316  cout << "MODOP=2 -- blocks emission in Kl to gamma e+e-: DEFAULT" << endl ;
317  cout << "MODOP=2 -- blocks emission in pi0 to gamma e+e-: DEFAULT" << endl ;
319  Photos::suppressBremForDecay(3, 130,22,11,-11);
320  Photos::suppressBremForDecay(3,-130,22,11,-11);
321  }
322 }
323 
324  // Previously this functionality was encoded in FORTRAN routine
325  // PHOCHK which was having some other functionality as well
326 bool Photos::IPHQRK_setQarknoEmission(int MODCOR, int PDGID)
327 {
328  static int IPHQRK_MODOP=-1;
329  if(IPHQRK_MODOP==-1 && MODCOR==0){
330  cout << "stop from IPHQRK_setQarknoEmission lack of initialization" << endl ;
331  exit(-1);
332  }
333  else if (MODCOR != 0){
334  IPHQRK_MODOP = MODCOR;
335  if(MODCOR ==1) cout << " IPHQRK_setQarknoEmission MODOP=1 -- blocks emission from light quarks: DEFAULT" << endl ;
336  if(MODCOR !=1) cout << " IPHQRK_setQarknoEmission MODOP=2 -- emission from light quarks allowed: TEST " << endl ;
337  }
338  if(IPHQRK_MODOP!=1) return true;
339 
340  return PDGID>9;
341 }
342 
343 void Photos::createHistoryEntries(bool flag, int status)
344 {
345  if(status<3)
346  {
347  Log::Warning()<<"Photos::createHistoryEntries: status must be >=3"<<endl;
348  return;
349  }
350 
351  isCreateHistoryEntries = flag;
352  historyEntriesStatus = status;
353  ignoreParticlesOfStatus(status);
354 }
355 
357 {
358  if(status<3)
359  {
360  Log::Warning()<<"Photos::ignoreParticlesOfStatus: status must be >=3"<<endl;
361  return;
362  }
363 
364  if(!ignoreStatusCodeList) ignoreStatusCodeList = new vector<int>();
365 
366  // Do not add duplicate entries to the list
367  for(unsigned int i=0;i<ignoreStatusCodeList->size();i++)
368  if( status==ignoreStatusCodeList->at(i) ) return;
369 
370  ignoreStatusCodeList->push_back(status);
371 }
372 
374 {
375  if(!ignoreStatusCodeList) return;
376 
377  for(unsigned int i=0;i<ignoreStatusCodeList->size();i++)
378  {
379  if( status==ignoreStatusCodeList->at(i) )
380  {
381  ignoreStatusCodeList->erase(ignoreStatusCodeList->begin()+i);
382  return;
383  }
384  }
385 }
386 
388 {
389  if(!ignoreStatusCodeList) return false;
390 
391  for(unsigned int i=0;i<ignoreStatusCodeList->size();i++)
392  if( status==ignoreStatusCodeList->at(i) ) return true;
393 
394  return false;
395 }
396 
397 void Photos::setRandomGenerator( double (*gen)() )
398 {
399  if(gen==NULL) randomDouble = PhotosRandom::randomReal;
400  else randomDouble = gen;
401 }
402 
404 {
405  phokey.iexp = (int) expo;
406  if(expo)
407  {
408  setDoubleBrem(false);
409  setQuatroBrem(false);
410  setInfraredCutOff(0.0000001);
412  phokey.expeps=0.0001;
413  }
414 }
415 
416  void Photos::setPairEmission(bool corr)
417 {
418  IfPair=corr;
419 }
420 
422 {
423  IfPhot=corr;
424 }
425 
427 {
428  meCorrectionWtForW=corr;
429 }
430 
432 {
433  meCorrectionWtForZ=corr;
434 }
436 {
438 }
439 
440 void Photos::setStopAtCriticalError(bool stop)
441 {
442  phosta.ifstop=(int)stop;
443  if(!stop)
444  {
445  Log::Info()<<"PHOTOS production mode. Elementary test of data flow from event record disabled. "<<endl
446  <<"Prior checks of the complete configuration "<<endl
447  <<"(for the particular set of input parameters) must have been done! "<<endl;
448  }
449 }
450 
451 
453 {
454  if(!forceMassList) forceMassList = new vector<pair<int,double>* >();
455  forceMassList->push_back( new pair<int,double>(pdgid, -1.0) );
456 }
457 
458 void Photos::forceMass(int pdgid, double mass)
459 {
460  if(mass<0.0)
461  {
462  Log::Warning()<<"Photos::forceMass: Mass must be > 0.0"<<endl;
463  return;
464  }
465 
466  if(!forceMassList) forceMassList = new vector<pair<int,double>* >();
467  forceMassList->push_back( new pair<int,double>(pdgid, mass) );
468 }
469 
470 } // namespace Photospp
static void setInfraredCutOff(double cut_off)
Definition: Photos.h:100
static void createHistoryEntries(bool flag, int status)
Definition: Photos.cxx:343
static void suppressBremForDecay(int count, int motherID,...)
Definition: Photos.cxx:238
static bool meCorrectionWtForW
Definition: Photos.h:212
static void setPairEmission(bool ifpair)
Definition: Photos.cxx:416
std::vector< PhotosParticle * > getDecayTree()
static bool IfPair
Definition: Photos.h:218
static void ignoreParticlesOfStatus(int status)
Definition: Photos.cxx:356
static bool isCreateHistoryEntries
Definition: Photos.h:215
static void processBranch(PhotosParticle *p)
Definition: Photos.cxx:231
static void setAlphaQED(double alpha)
Definition: Photos.h:103
static bool meCorrectionWtForZ
Definition: Photos.h:209
static vector< vector< int > * > * supBremList
Definition: Photos.h:191
static void setQuatroBrem(bool quatroBrem)
Definition: Photos.h:112
static void deIgnoreParticlesOfStatus(int status)
Definition: Photos.cxx:373
static void forceBremForDecay(int count, int motherID,...)
Definition: Photos.cxx:270
static void forceMassFromEventRecord(int pdgid)
Definition: Photos.cxx:452
static int EventNo
Definition: Photos.h:182
static bool massFrom4Vector
Definition: Photos.h:188
static bool IfPhot
Definition: Photos.h:221
static void IPHEKL_setPi0KLnoEmission(int m)
Definition: Photos.cxx:304
static void setMeCorrectionWtForZ(bool corr)
Definition: Photos.cxx:431
static vector< int > * ignoreStatusCodeList
Definition: Photos.h:200
static double momentum_conservation_threshold
Definition: Photos.h:203
static void iniInfo()
Definition: Photos.cxx:181
static void setPhotonEmission(bool ifphot)
Definition: Photos.cxx:421
static int historyEntriesStatus
Definition: Photos.h:225
Controls the configuration and initialization of Photos.
static bool isSuppressed
Definition: Photos.h:185
static bool isStatusCodeIgnored(int status)
Definition: Photos.cxx:387
static void setRandomGenerator(double(*gen)())
Definition: Photos.cxx:397
static void initialize()
Definition: Photos.cxx:53
static void setMeCorrectionWtForScalar(bool corr)
Definition: Photos.cxx:435
static void processParticle(PhotosParticle *p)
Definition: Photos.cxx:225
static void setInterference(bool interference)
Definition: Photos.h:106
static void initializeKinematicCorrections(int flag)
Definition: Photos.h:149
static void setMeCorrectionWtForW(bool corr)
Definition: Photos.cxx:426
static bool meCorrectionWtForScalar
Definition: Photos.h:206
static void maxWtInterference(double interference)
Definition: Photos.h:97
static void setExponentiation(bool expo)
Definition: Photos.cxx:403
static double(* randomDouble)()
Definition: Photos.h:228
static vector< pair< int, double > * > * forceMassList
Definition: Photos.h:197
static void suppressBremForBranch(int count, int motherID,...)
Definition: Photos.cxx:254
static void setTopProcessRadiation(bool top)
Definition: Photos.h:136
static vector< PhotosBranch * > createBranches(vector< PhotosParticle * > particles)
static void setDoubleBrem(bool doub)
Definition: Photos.h:109
static void forceMass(int pdgid, double mass)
Definition: Photos.cxx:458
static void forceBremForBranch(int count, int motherID,...)
Definition: Photos.cxx:286
static vector< vector< int > * > * forceBremList
Definition: Photos.h:194