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"
9using std::vector;
10using std::cout;
11using std::endl;
12using std::ios_base;
13
14namespace Photospp
15{
16
17Photos Photos::_instance;
18
19vector<vector<int>* > *Photos::supBremList = 0;
20vector<vector<int>* > *Photos::forceBremList = 0;
21vector<pair<int,double>* > *Photos::forceMassList = 0;
22vector<int> *Photos::ignoreStatusCodeList = 0;
23bool Photos::isSuppressed=false;
31bool Photos::IfPair=false;
32bool Photos::IfPhot=true;
33int Photos::EventNo=0;
34double (*Photos::randomDouble)() = PhotosRandom::randomReal;
35
36Photos::MomentumUnits Photos::momentumUnit = Photos::DEFAULT_MOMENTUM;
37
38Photos::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
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
238void 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
254void 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
270void 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
286void 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
326bool 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
343void 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
352 historyEntriesStatus = 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 {
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
397void 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
417{
418 IfPair=corr;
419}
420
422{
423 IfPhot=corr;
424}
425
427{
429}
430
432{
434}
436{
438}
439
440void 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
458void 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
Controls the configuration and initialization of Photos.
static vector< PhotosBranch * > createBranches(vector< PhotosParticle * > particles)
std::vector< PhotosParticle * > getDecayTree()
static void forceBremForDecay(int count, int motherID,...)
Definition Photos.cxx:270
static void setQuatroBrem(bool quatroBrem)
static void forceMass(int pdgid, double mass)
Definition Photos.cxx:458
static bool isCreateHistoryEntries
static void initialize()
Definition Photos.cxx:53
static bool isStatusCodeIgnored(int status)
Definition Photos.cxx:387
static void suppressBremForBranch(int count, int motherID,...)
Definition Photos.cxx:254
static void setMeCorrectionWtForW(bool corr)
Definition Photos.cxx:426
static vector< pair< int, double > * > * forceMassList
static void setRandomGenerator(double(*gen)())
Definition Photos.cxx:397
static void suppressBremForDecay(int count, int motherID,...)
Definition Photos.cxx:238
static void createHistoryEntries(bool flag, int status)
Definition Photos.cxx:343
static void setExponentiation(bool expo)
Definition Photos.cxx:403
static void setInfraredCutOff(double cut_off)
static void forceBremForBranch(int count, int motherID,...)
Definition Photos.cxx:286
static void setTopProcessRadiation(bool top)
static bool meCorrectionWtForScalar
static void setInterference(bool interference)
static void ignoreParticlesOfStatus(int status)
Definition Photos.cxx:356
static void forceMassFromEventRecord(int pdgid)
Definition Photos.cxx:452
static void maxWtInterference(double interference)
static void setMeCorrectionWtForZ(bool corr)
Definition Photos.cxx:431
static void setPhotonEmission(bool ifphot)
Definition Photos.cxx:421
static double(* randomDouble)()
static vector< vector< int > * > * supBremList
static vector< vector< int > * > * forceBremList
static void setMeCorrectionWtForScalar(bool corr)
Definition Photos.cxx:435
static vector< int > * ignoreStatusCodeList
static void IPHEKL_setPi0KLnoEmission(int m)
Definition Photos.cxx:304
static void processBranch(PhotosParticle *p)
Definition Photos.cxx:231
static void processParticle(PhotosParticle *p)
Definition Photos.cxx:225
static void setPairEmission(bool ifpair)
Definition Photos.cxx:416
static void setDoubleBrem(bool doub)
static double momentum_conservation_threshold
static void initializeKinematicCorrections(int flag)
static void iniInfo()
Definition Photos.cxx:181
static void deIgnoreParticlesOfStatus(int status)
Definition Photos.cxx:373
static void setAlphaQED(double alpha)