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