]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STARLIGHT/starlight/src/inputParameters.cpp~
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
STARLIGHT code and interface
[u/mrichter/AliRoot.git] / STARLIGHT / starlight / src / inputParameters.cpp~
CommitLineData
da32329d
AM		 1///////////////////////////////////////////////////////////////////////////
2//
3// Copyright 2010
4//
5// This file is part of starlight.
6//
7// starlight is free software: you can redistribute it and/or modify
8// it under the terms of the GNU General Public License as published by
9// the Free Software Foundation, either version 3 of the License, or
10// (at your option) any later version.
11//
12// starlight is distributed in the hope that it will be useful,
13// but WITHOUT ANY WARRANTY; without even the implied warranty of
14// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15// GNU General Public License for more details.
16//
17// You should have received a copy of the GNU General Public License
18// along with starlight. If not, see <http://www.gnu.org/licenses/>.
19//
20///////////////////////////////////////////////////////////////////////////
21//
22// File and Version Information:
23// $Rev:: 166 $: revision of last commit
24// $Author:: odjuvsla $: author of last commit
25// $Date:: 2013-10-06 16:18:12 +0200 #$: date of last commit
26//
27// Description:
28//
29//
30//
31///////////////////////////////////////////////////////////////////////////
32
33
34#include <iostream>
35#include <fstream>
36
37#include "reportingUtils.h"
38#include "starlightconstants.h"
39#include "inputParameters.h"
40#include "inputParser.h"
41#include "starlightconfig.h"
42#include <cmath>
43#include <cstring>
44#include "randomgenerator.h"
45
46using namespace std;
47using namespace starlightConstants;
48
49parameterlist parameterbase::_parameters;
50
51#define REQUIRED true
52#define NOT_REQUIRED false
53
54//______________________________________________________________________________
55inputParameters::inputParameters()
56 : _configFileName ("slight.in"),
57 _beam1Z ("BEAM_1_Z",0),
58 _beam1A ("BEAM_1_A",0),
59 _beam2Z ("BEAM_2_Z",0),
60 _beam2A ("BEAM_2_A",0),
61 _beam1LorentzGamma ("BEAM_1_GAMMA",0),
62 _beam2LorentzGamma ("BEAM_2_GAMMA",0),
63 _maxW ("W_MAX",0),
64 _minW ("W_MIN",0),
65 _nmbWBins ("W_N_BINS",0),
66 _maxRapidity ("RAP_MAX",0),
67 _nmbRapidityBins ("RAP_N_BINS",0),
68 _ptCutEnabled ("CUT_PT",false),
69 _ptCutMin ("PT_MIN",0),
70 _ptCutMax ("PT_MAX",0),
71 _etaCutEnabled ("CUT_ETA",false),
72 _etaCutMin ("ETA_MIN",0),
73 _etaCutMax ("ETA_MAX",0),
74 _productionMode ("PROD_MODE",0),
75 _nmbEventsTot ("N_EVENTS",0),
76 _prodParticleId ("PROD_PID",0),
77 _randomSeed ("RND_SEED",0),
78 _outputFormat ("OUTPUT_FORMAT",0),
79 _beamBreakupMode ("BREAKUP_MODE",0),
80 _interferenceEnabled ("INTERFERENCE",false),
81 _interferenceStrength ("IF_STRENGTH",0),
82 _coherentProduction ("COHERENT",false),
83 _incoherentFactor ("INCO_FACTOR",0),
84 _deuteronSlopePar ("BFORD",0),
85 _maxPtInterference ("INT_PT_MAX",0),
86 _nmbPtBinsInterference ("INT_PT_N_BINS",0),
87 _ptBinWidthInterference("INT_PT_WIDTH",0),
88 _protonEnergy ("PROTON_ENERGY",0, NOT_REQUIRED),
89 _minGammaEnergy ("MIN_GAMMA_ENERGY",0, NOT_REQUIRED),
90 _maxGammaEnergy ("MAX_GAMMA_ENERGY",0, NOT_REQUIRED),
91 _pythiaParams ("PYTHIA_PARAMS","", NOT_REQUIRED),
92 _pythiaFullEventRecord ("PYTHIA_FULL_EVENTRECORD",false, NOT_REQUIRED),
93 _xsecCalcMethod ("XSEC_METHOD",0, NOT_REQUIRED),
94 _nThreads ("N_THREADS",1, NOT_REQUIRED),
95 _nBinsQKniehl ("N_BINS_Q_KNIEHL", 0, NOT_REQUIRED),
96 _nBinsEKniehl ("N_BINS_E_KNIEHL", 0, NOT_REQUIRED),
97 _nBinsBKniehl ("N_BINS_B_KNIEHL", 0, NOT_REQUIRED),
98 _qMaxKniehl ("Q_MAX_KNIEHL", 0, NOT_REQUIRED),
99 _eGammaMinKniehl ("E_GAMMA_MIN_KNIEHL", 0, NOT_REQUIRED),
100 _eGammaMaxKniehl ("E_GAMMA_MAX_KNIEHL", 0, NOT_REQUIRED),
101 _bMinKniehl ("B_MIN_KNIEHL", 0, NOT_REQUIRED),
102 _bMaxKniehl ("B_MAX_KNIEHL", 0, NOT_REQUIRED)
103{
104 // All parameters must be initialised in initialisation list!
105 // If not: error: 'parameter<T, validate>::parameter() [with T = unsigned int, bool validate = true]' is private
106 // or similar
107
108 _ip.addParameter(_beam1Z);
109 _ip.addParameter(_beam2Z);
110 _ip.addParameter(_beam1A);
111 _ip.addParameter(_beam2A);
112
113 _ip.addParameter(_beam1LorentzGamma);
114 _ip.addParameter(_beam2LorentzGamma);
115
116 _ip.addParameter(_maxW);
117 _ip.addParameter(_minW);
118
119 _ip.addParameter(_nmbWBins);
120
121 _ip.addParameter(_maxRapidity);
122 _ip.addParameter(_nmbRapidityBins);
123
124 _ip.addParameter(_ptCutEnabled);
125 _ip.addParameter(_ptCutMin);
126 _ip.addParameter(_ptCutMax);
127
128 _ip.addParameter(_etaCutEnabled);
129 _ip.addParameter(_etaCutMax);
130 _ip.addParameter(_etaCutMin);
131
132 _ip.addParameter(_productionMode);
133 _ip.addParameter(_nmbEventsTot);
134 _ip.addParameter(_prodParticleId);
135 _ip.addParameter(_randomSeed);
136 _ip.addParameter(_outputFormat);
137 _ip.addParameter(_beamBreakupMode);
138 _ip.addParameter(_interferenceEnabled);
139 _ip.addParameter(_interferenceStrength);
140 _ip.addParameter(_coherentProduction);
141 _ip.addParameter(_incoherentFactor);
142 _ip.addParameter(_deuteronSlopePar);
143 _ip.addParameter(_maxPtInterference);
144 _ip.addParameter(_nmbPtBinsInterference);
145 _ip.addParameter(_minGammaEnergy);
146 _ip.addParameter(_maxGammaEnergy);
147 _ip.addParameter(_pythiaParams);
148 _ip.addParameter(_pythiaFullEventRecord);
149 _ip.addParameter(_xsecCalcMethod);
150 _ip.addParameter(_nThreads);
151 _ip.addParameter(_nBinsBKniehl);
152 _ip.addParameter(_nBinsQKniehl);
153 _ip.addParameter(_nBinsEKniehl);
154 _ip.addParameter(_qMaxKniehl);
155 _ip.addParameter(_eGammaMaxKniehl);
156 _ip.addParameter(_eGammaMinKniehl);
157 _ip.addParameter(_bMaxKniehl);
158 _ip.addParameter(_bMinKniehl);
159}
160
161
162//______________________________________________________________________________
163inputParameters::~inputParameters()
164{ }
165
166
167//______________________________________________________________________________
168bool
169inputParameters::configureFromFile(const std::string &configFileName)
170{
171 // open config file
172 _configFileName = configFileName;
173
174 int nParameters = _ip.parseFile(_configFileName);
175
176 if(nParameters == -1)
177 {
178 printWarn << "could not open file '" << _configFileName << "'" << endl;
179 return false;
180 }
181
182 //ip.printParameterInfo(cout);
183
184 if(_ip.validateParameters(cerr))
185 printInfo << "successfully read input parameters from '" << _configFileName << "'" << endl;
186 else {
187 printWarn << "problems reading input parameters from '" << _configFileName << "'" << endl
188 << *this;
189 return false;
190 }
191 return true;
192}
193 bool inputParameters::init()
194 {
195 // Set the seed for the random generator
196 randyInstance.SetSeed(_randomSeed.value());
197
198 // Calculate beam gamma in CMS frame
199 double rap1 = acosh(beam1LorentzGamma());
200 double rap2 = -acosh(beam2LorentzGamma());
201 _beamLorentzGamma = cosh((rap1-rap2)/2);
202
203 std::cout << "Rapidity beam 1: " << rap1 << ", rapidity beam 2: " << rap2 << ", rapidity CMS system: " << (rap1+rap2)/2 << ", beam gamma in CMS: " << _beamLorentzGamma<< std::endl;
204 _ptBinWidthInterference = maxPtInterference() / nmbPtBinsInterference();
205 _protonEnergy = _beamLorentzGamma * protonMass;
206
207 // define interaction type
208 switch (productionMode()) {
209 case 1:
210 _interactionType = PHOTONPHOTON;
211 break;
212 case 2:
213 _interactionType = PHOTONPOMERONNARROW;
214 break;
215 case 3:
216 _interactionType = PHOTONPOMERONWIDE;
217 break;
218 case 4:
219 _interactionType = PHOTONPOMERONINCOHERENT;
220 break;
221 case 5:
222 _interactionType = PHOTONUCLEARSINGLE;
223 break;
224 case 6:
225 _interactionType = PHOTONUCLEARDOUBLE;
226 break;
227 case 7:
228 _interactionType = PHOTONUCLEARSINGLEPA;
229 break;
230 case 8:
231 _interactionType = PHOTONUCLEARSINGLEPAPY;
232 break;
233// case 9:
234// _interactionType = PHOTONPHOTONKNIEHL;
235// break;
236// case 10:
237// _interactionType = PHOTONPHOTONKNIEHLMODIFIED;
238// break;
239 default:
240 printWarn << "unknown production mode '" << _productionMode << "'" << endl;
241 return false;
242 }
243
244 //Trying to define the proper Wmins and Wmaxs. a TEMPORARY fix....Better solution=??
245 double mass = 0;
246 double width = 0;
247 double defaultMinW = 0; // default for _minW, unless it is defined later [GeV/c^2]
248 switch (prodParticleId()) {
249 case 11: // e+e- pair
250 _particleType = ELECTRON;
251 _decayType = LEPTONPAIR;
252 defaultMinW = 0.01; // default is 0.01; up to 0.15 is safe for Summer 2000 triggering for e+e- pairs
253 break;
254 case 13: // mu+mu- pair
255 _particleType = MUON;
256 _decayType = LEPTONPAIR;
257 defaultMinW = 2 * muonMass;
258 break;
259 case 15: // tau+tau- pair
260 _particleType = TAUON;
261 _decayType = LEPTONPAIR;
262 defaultMinW = 2 * tauMass;
263 break;
264// case 24: // W+W- pair
265// _particleType = W;
266// _decayType = WW;
267// defaultMinW = 2 * muonMass;
268// break;
269 case 115: // a_2(1320)
270 _particleType = A2;
271 _decayType = SINGLEMESON;
272 break;
273 case 221: // eta
274 _particleType = ETA;
275 _decayType = SINGLEMESON;
276 break;
277 case 225: // f_2(1270)
278 _particleType = F2;
279 defaultMinW = 2*pionChargedMass;
280 _decayType = SINGLEMESON;
281 break;
282 case 331: // eta'(958)
283 _particleType = ETAPRIME;
284 _decayType = SINGLEMESON;
285 break;
286 case 335: // f_2'(1525)
287 _particleType = F2PRIME;
288 _decayType = SINGLEMESON;
289 break;
290 case 441: // eta_c(1s)
291 _particleType = ETAC;
292 _decayType = SINGLEMESON;
293 defaultMinW = etaCMass - 5 * 0.0267;
294 break;
295 case 9010221: // f_0(980), was orginally called 10221? updated to standard number
296 _particleType = F0;
297 _decayType = SINGLEMESON;
298 defaultMinW = 2*pionNeutralMass;
299 break;
300 case 33: // Z"/Z03
301 _particleType = ZOVERZ03;
302 _decayType = SINGLEMESON;
303 break;
304// case 25: // Higgs
305// _particleType = HIGGS;
306// _decayType = SINGLEMESON;
307// break;
308 case 113: // rho(770)
309 _particleType = RHO;
310 _decayType = WIDEVMDEFAULT;
311 mass = 0.7685;
312 width = 0.1507;
313 defaultMinW = 2 * pionChargedMass;
314 _maxW = mass + 5 * width;
315 break;
316 case 913: // rho(770) with direct pi+pi- decay, interference given by ZEUS data
317 _particleType = RHOZEUS;
318 _decayType = WIDEVMDEFAULT;
319 mass = 0.7685;
320 width = 0.1507;
321 defaultMinW = 2 * pionChargedMass;
322 _maxW = mass + 5 * width; // use the same 1.5GeV max mass as ZEUS
323 break;
324 case 999: // pi+pi-pi+pi- phase space decay
325 _particleType = FOURPRONG;
326 _decayType = WIDEVMDEFAULT;
327 mass = 1.350;
328 width = 0.360;
329 defaultMinW = 4 * pionChargedMass;
330 _maxW = 3;
331 break;
332 case 223: // omega(782)
333 _particleType = OMEGA;
334 _decayType = NARROWVMDEFAULT; // will probably be moved to 3-body decay
335 mass = 0.78194;
336 width = 0.00843;
337 defaultMinW = mass - 5 * width;
338 _maxW = mass + 5 * width;
339 break;
340 case 333: // phi(1020)
341 _particleType = PHI;
342 _decayType = NARROWVMDEFAULT;
343 mass = 1.019413;
344 width = 0.00443;
345 defaultMinW = 2 * kaonChargedMass;
346 _maxW = mass + 5 * width;
347 break;
348 case 443: // J/psi
349 _particleType = JPSI;
350 _decayType = NARROWVMDEFAULT;
351 mass = 3.09692; // JN 3.09688;
352 width = 0.000091; // JN 0.000087;
353 defaultMinW = mass - 5 * width;
354 _maxW = mass + 5 * width;
355 break;
356 case 443011: // J/psi
357 _particleType = JPSI_ee;
358 _decayType = NARROWVMDEFAULT;
359 mass = 3.09692; // JN 3.09688;
360 width = 0.000091; // JN 0.000087;
361 defaultMinW = mass - 5 * width;
362 _maxW = mass + 5 * width;
363 break;
364 case 443013: // J/psi
365 _particleType = JPSI_mumu;
366 _decayType = NARROWVMDEFAULT;
367 mass = 3.09692; // JN 3.09688;
368 width = 0.000091; // JN 0.000087;
369 defaultMinW = mass - 5 * width;
370 _maxW = mass + 5 * width;
371 break;
372 case 444: // J/psi
373 _particleType = JPSI2S;
374 _decayType = NARROWVMDEFAULT;
375 mass = 3.686093;
376 width = 0.000337;
377 defaultMinW = mass - 5 * width;
378 _maxW = mass + 5 * width;
379 break;
380 case 444011: // J/psi
381 _particleType = JPSI2S_ee;
382 _decayType = NARROWVMDEFAULT;
383 mass = 3.686093;
384 width = 0.000337;
385 defaultMinW = mass - 5 * width;
386 _maxW = mass + 5 * width;
387 break;
388 case 444013: // J/psi
389 _particleType = JPSI2S_mumu;
390 _decayType = NARROWVMDEFAULT;
391 mass = 3.686093;
392 width = 0.000337;
393 defaultMinW = mass - 5 * width;
394 _maxW = mass + 5 * width;
395 break;
396 case 553: // Upsilon
397 _particleType = UPSILON;
398 _decayType = NARROWVMDEFAULT;
399 mass = 9.46030;
400 width = 0.00005402;
401 defaultMinW = mass - 5 * width;
402 _maxW = mass + 5 * width;
403 break;
404 case 553011: // Upsilon
405 _particleType = UPSILON_ee;
406 _decayType = NARROWVMDEFAULT;
407 mass = 9.46030;
408 width = 0.00005402;
409 defaultMinW = mass - 5 * width;
410 _maxW = mass + 5 * width;
411 break;
412 case 553013: // Upsilon
413 _particleType = UPSILON_mumu;
414 _decayType = NARROWVMDEFAULT;
415 mass = 9.46030;
416 width = 0.00005402;
417 defaultMinW = mass - 5 * width;
418 _maxW = mass + 5 * width;
419 break;
420 case 554: // Upsilon(2S)
421 _particleType = UPSILON2S;
422 _decayType = NARROWVMDEFAULT;
423 mass = 10.02326;
424 width = 0.00003198;
425 defaultMinW = mass - 5 * width;
426 _maxW = mass + 5 * width;
427 break;
428 case 554011: // Upsilon(2S)
429 _particleType = UPSILON2S_ee;
430 _decayType = NARROWVMDEFAULT;
431 mass = 10.02326;
432 width = 0.00003198;
433 defaultMinW = mass - 5 * width;
434 _maxW = mass + 5 * width;
435 break;
436 case 554013: // Upsilon(2S)
437 _particleType = UPSILON2S_mumu;
438 _decayType = NARROWVMDEFAULT;
439 mass = 10.02326;
440 width = 0.00003198;
441 defaultMinW = mass - 5 * width;
442 _maxW = mass + 5 * width;
443 break;
444 case 555: // Upsilon(3S)
445 mass = 10.3552;
446 width = 0.00002032;
447 defaultMinW = mass - 5 * width;
448 _maxW = mass + 5 * width;
449 _particleType = UPSILON3S;
450 _decayType = NARROWVMDEFAULT;
451 break;
452 case 555011: // Upsilon(3S)
453 _particleType = UPSILON3S_ee;
454 _decayType = NARROWVMDEFAULT;
455 mass = 10.3552;
456 width = 0.00002032;
457 defaultMinW = mass - 5 * width;
458 _maxW = mass + 5 * width;
459 break;
460 case 555013: // Upsilon(3S)
461 _particleType = UPSILON3S_mumu;
462 _decayType = NARROWVMDEFAULT;
463 mass = 10.3552;
464 width = 0.00002032;
465 defaultMinW = mass - 5 * width;
466 _maxW = mass + 5 * width;
467 break;
468 default:
469 printWarn << "unknown particle ID " << _prodParticleId << endl;
470 return false;
471 } // _prodParticleId
472
473 if (_minW.value() == -1)
474 _minW = defaultMinW;
475
476 printInfo << "using the following " << *this;
477
478 return true;
479}
480
481
482//______________________________________________________________________________
483ostream&
484inputParameters::print(ostream& out) const
485{
486 out << "starlight parameters:" << endl
487 << " config file name ...................... '" << _configFileName << "'" << endl
488 << " beam 1 atomic number ................... " << _beam1Z.value() << endl
489 << " beam 1 atomic mass number .............. " << _beam1A.value() << endl
490 << " beam 2 atomic number ................... " << _beam2Z.value() << endl
491 << " beam 2 atomic mass number .............. " << _beam2A.value() << endl
492 << " Lorentz gamma of beams in CM frame ..... " << _beamLorentzGamma << endl
493 << " mass W of produced hadronic system ..... " << _minW.value() << " < W < " << _maxW.value() << " GeV/c^2" << endl
494 << " # of W bins ............................ " << _nmbWBins.value() << endl
495 << " maximum absolute value for rapidity .... " << _maxRapidity.value() << endl
496 << " # of rapidity bins ..................... " << _nmbRapidityBins.value() << endl
497 << " cut in pT............................... " << yesNo(_ptCutEnabled.value()) << endl
498 << " minumum pT.......................... " << _ptCutMin.value() << " GeV/c" << endl
499 << " maximum pT.......................... " << _ptCutMax.value() << " GeV/c" << endl
500 << " cut in eta.............................. " << yesNo(_etaCutEnabled.value()) << endl
501 << " minumum eta......................... " << _etaCutMin.value() << endl
502 << " maximum eta......................... " << _etaCutMax.value() << endl
503 << " meson production mode .................. " << _productionMode.value() << endl
504 << " number of events to generate ........... " << _nmbEventsTot.value() << endl
505 << " PDG ID of produced particle ............ " << _prodParticleId.value() << endl
506 << " seed for random generator .............. " << _randomSeed.value() << endl
507 << " output format .......................... " << _outputFormat.value() << endl
508 << " breakup mode for beam particles ........ " << _beamBreakupMode.value() << endl
509 << " interference enabled ................... " << yesNo(_interferenceEnabled.value()) << endl
510 << " interference strength .................. " << _interferenceStrength.value() << endl
511 << " coherent scattering off nucleus ........ " << yesNo(_coherentProduction.value()) << endl
512 << " scaling factor for incoh. VM prod. ..... " << _incoherentFactor.value() << endl
513 << " deuteron slope parameter ............... " << _deuteronSlopePar.value() << " (GeV/c)^{-2}" << endl
514 << " maximum p_T for interference calc. ..... " << _maxPtInterference.value() << " GeV/c" << endl
515 << " # of p_T bins for interference calc. ... " << _nmbPtBinsInterference.value() << endl;
516 return out;
517}
518
519
520//______________________________________________________________________________
521ostream&
522inputParameters::write(ostream& out) const
523{
524
525 out << "BEAM_1_Z" << beam1Z () <<endl
526 << "BEAM_2_Z" << beam1A () <<endl
527 << "BEAM_1_A" << beam2Z () <<endl
528 << "BEAM_2_A" << beam2A () <<endl
529 << "BEAM_GAMMA" << beamLorentzGamma () <<endl
530 << "W_MAX" << maxW () <<endl
531 << "W_MIN" << minW () <<endl
532 << "W_N_BINS" << nmbWBins () <<endl
533 << "RAP_MAX" << maxRapidity () <<endl
534 << "RAP_N_BINS" << nmbRapidityBins () <<endl
535 << "CUT_PT" << ptCutEnabled () <<endl
536 << "PT_MIN" << ptCutMin () <<endl
537 << "PT_MAX" << ptCutMax () <<endl
538 << "CUT_ETA" << etaCutEnabled () <<endl
539 << "ETA_MIN" << etaCutMin () <<endl
540 << "ETA_MAX" << etaCutMax () <<endl
541 << "PROD_MODE" << productionMode () <<endl
542 << "N_EVENTS" << nmbEvents () <<endl
543 << "PROD_PID" << prodParticleId () <<endl
544 << "RND_SEED" << randomSeed () <<endl
545 << "OUTPUT_FORMAT" << outputFormat () <<endl
546 << "BREAKUP_MODE" << beamBreakupMode () <<endl
547 << "INTERFERENCE" << interferenceEnabled () <<endl
548 << "IF_STRENGTH" << interferenceStrength () <<endl
549 << "COHERENT" << coherentProduction () <<endl
550 << "INCO_FACTOR" << incoherentFactor () <<endl
551 << "BFORD" << deuteronSlopePar () <<endl
552 << "INT_PT_MAX" << maxPtInterference () <<endl
553 << "INT_PT_N_BINS" << nmbPtBinsInterference() <<endl;
554
555 return out;
556}
557
558std::string
559inputParameters::parameterValueKey() const
560{
561
562 // std::stringstream s;
563
564// s <<_beam1A<<_beam1Z<<_beam2A<<_beam1LorentzGamma<<_beam2LorentzGamma<<_maxW<<_minW;
565// <<_nmbWBins<<_maxRapidity<<_nmbRapidityBins<<_
566
567
568 //return s;
569
570 return parameterbase::_parameters.validationKey()
571 ;
572}