EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
inputParameters.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file inputParameters.cpp
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 //
21 //
22 // File and Version Information:
23 // $Rev:: 276 $: revision of last commit // $Author:: jnystra#$: author of last commit
24 // $Date:: 2016-09-13 19:54:42 +0100 #$: date of last commit //
25 // Description:
26 //
27 //
28 //
30 
31 
32 #include <iostream>
33 #include <fstream>
34 
35 #include "reportingUtils.h"
36 #include "starlightconstants.h"
37 #include "inputParameters.h"
38 #include "inputParser.h"
39 #include "starlightconfig.h"
40 #include <cmath>
41 #include <cstring>
42 #include "randomgenerator.h"
43 
44 using namespace std;
45 using namespace starlightConstants;
46 
47 parameterlist parameterbase::_parameters;
48 
49 #define REQUIRED true
50 #define NOT_REQUIRED false
51 
52 //______________________________________________________________________________
53 inputParameters::inputParameters()
54  : _baseFileName ("baseFileName","slight"),
55  _targetBeamZ ("TARGET_BEAM_Z",0),
56  _targetBeamA ("TARGET_BEAM_A",0),
57  _electronBeamLorentzGamma ("ELECTRON_BEAM_GAMMA",0),
58  _targetBeamLorentzGamma ("TARGET_BEAM_GAMMA",0),
59  _maxW ("W_MAX",0),
60  _minW ("W_MIN",0),
61  _maxW_GP ("W_GP_MAX", 0,NOT_REQUIRED),
62  _minW_GP ("W_GP_MIN", 0,NOT_REQUIRED),
63  _nmbWBins ("W_N_BINS",0),
64  _maxRapidity ("RAP_MAX",9,NOT_REQUIRED),
65  _nmbRapidityBins ("RAP_N_BINS",100,NOT_REQUIRED),
66  _nmbEnergyBins ("EGA_N_BINS",0),
67  _ptCutEnabled ("CUT_PT",false, NOT_REQUIRED),
68  _ptCutMin ("PT_MIN",0, NOT_REQUIRED),
69  _ptCutMax ("PT_MAX",0, NOT_REQUIRED),
70  _etaCutEnabled ("CUT_ETA",false, NOT_REQUIRED),
71  _etaCutMin ("ETA_MIN",0, NOT_REQUIRED),
72  _etaCutMax ("ETA_MAX",0, NOT_REQUIRED),
73  _productionMode ("PROD_MODE",0),
74  _nmbEventsTot ("N_EVENTS",0),
75  _prodParticleId ("PROD_PID",0),
76  _randomSeed ("RND_SEED",0),
77  _beamBreakupMode ("BREAKUP_MODE",0),
78  _interferenceEnabled ("INTERFERENCE",false,NOT_REQUIRED),
79  _interferenceStrength ("IF_STRENGTH",0,NOT_REQUIRED),
80  _maxPtInterference ("INT_PT_MAX",0,NOT_REQUIRED),
81  _nmbPtBinsInterference ("INT_PT_N_BINS",0),
82  _ptBinWidthInterference("INT_PT_WIDTH",0),
83  _protonEnergy ("PROTON_ENERGY",0, NOT_REQUIRED),
84  _electronEnergy ("ELECTRON_ENERGY",0, NOT_REQUIRED),
85  _minGammaEnergy ("MIN_GAMMA_ENERGY",6.0, NOT_REQUIRED),
86  _maxGammaEnergy ("MAX_GAMMA_ENERGY",600000.0, NOT_REQUIRED),
87  _minGammaQ2 ("MIN_GAMMA_Q2",0,NOT_REQUIRED),
88  _maxGammaQ2 ("MAX_GAMMA_Q2",0,NOT_REQUIRED),
89  _nmbGammaQ2Bins ("INT_GAMMA_Q2_BINS",400,NOT_REQUIRED),
90  _pythiaParams ("PYTHIA_PARAMS","", NOT_REQUIRED),
91  _outputFormat ("OUTPUT_FORMAT",0, NOT_REQUIRED),
92  _backwardsProduction ("BACKWARDS_PRODUCTION",false, NOT_REQUIRED),
93 
94  _xsecCalcMethod ("XSEC_METHOD",0, NOT_REQUIRED),
95  _axionMass ("AXION_MASS",50, NOT_REQUIRED), // AXION HACK
96  _bslopeDefinition ("BSLOPE_DEFINITION",0, NOT_REQUIRED),
97  _bslopeValue ("BSLOPE_VALUE",4.0,NOT_REQUIRED),
98  _impulseVM ("SELECT_IMPULSE_VM",0,NOT_REQUIRED),
99  _quantumGlauber ("QUANTUM_GLAUBER",0,NOT_REQUIRED)
100 {
101  // All parameters must be initialised in initialisation list!
102  // If not: error: 'parameter<T, validate>::parameter() [with T = unsigned int, bool validate = true]' is private
103  // or similar
104 
105  _ip.addParameter(_baseFileName);
106 
107  _ip.addParameter(_targetBeamZ);
108  _ip.addParameter(_targetBeamA);
109 
110  _ip.addParameter(_targetBeamLorentzGamma);
111  _ip.addParameter(_electronBeamLorentzGamma);
112 
113  _ip.addParameter(_maxW);
114  _ip.addParameter(_minW);
115 
116  _ip.addParameter(_maxW_GP);
117  _ip.addParameter(_minW_GP);
118 
119  _ip.addParameter(_nmbWBins);
120 
121  _ip.addParameter(_maxRapidity);
122  _ip.addParameter(_nmbRapidityBins);
123  //
124  _ip.addParameter(_nmbEnergyBins);
125 
126  _ip.addParameter(_ptCutEnabled);
127  _ip.addParameter(_ptCutMin);
128  _ip.addParameter(_ptCutMax);
129 
130  _ip.addParameter(_etaCutEnabled);
131  _ip.addParameter(_etaCutMax);
132  _ip.addParameter(_etaCutMin);
133 
134  _ip.addParameter(_productionMode);
135  _ip.addParameter(_nmbEventsTot);
136  _ip.addParameter(_prodParticleId);
137  _ip.addParameter(_randomSeed);
138  _ip.addParameter(_beamBreakupMode);
139  _ip.addParameter(_interferenceEnabled);
140  _ip.addParameter(_interferenceStrength);
141  _ip.addParameter(_maxPtInterference);
142  _ip.addParameter(_nmbPtBinsInterference);
143  _ip.addParameter(_minGammaEnergy);
144  _ip.addParameter(_maxGammaEnergy);
145  _ip.addParameter(_minGammaQ2);
146  _ip.addParameter(_maxGammaQ2);
147  _ip.addParameter(_nmbGammaQ2Bins);
148  _ip.addParameter(_pythiaParams);
149  _ip.addParameter(_outputFormat);
150  _ip.addParameter(_backwardsProduction);
151  _ip.addParameter(_xsecCalcMethod);
152  _ip.addParameter(_axionMass); // AXION HACK
153  _ip.addParameter(_bslopeDefinition);
154  _ip.addParameter(_bslopeValue);
155  _ip.addParameter(_impulseVM);
156  _ip.addParameter(_quantumGlauber);
157 }
158 
159 
160 //______________________________________________________________________________
161 inputParameters::~inputParameters()
162 { }
163 
164 
165 //______________________________________________________________________________
166 bool
167 inputParameters::configureFromFile(const std::string &_configFileName)
168 {
169 
170  int nParameters = _ip.parseFile(_configFileName);
171 
172  if(nParameters == -1)
173  {
174  printWarn << "could not open file '" << _configFileName << "'" << endl;
175  return false;
176  }
177 
178 
179  if(_ip.validateParameters(cerr))
180  printInfo << "successfully read input parameters from '" << _configFileName << "'" << endl;
181  else {
182  printWarn << "problems reading input parameters from '" << _configFileName << "'" << endl
183  << *this;
184  return false;
185  }
186  return true;
187 }
188  bool inputParameters::init()
189  {
190 
191  // Calculate beam gamma in CMS frame
192  double rap1 = acosh(electronBeamLorentzGamma());
193  double rap2 = -acosh(targetBeamLorentzGamma());
194 
195  double _electronEnergy_lab = (electronBeamLorentzGamma())*(starlightConstants::mel);
196  double _ionEnergy_lab = targetBeamLorentzGamma()*(targetBeamA())*protonMass;
197  double _ionPz_lab = -1.0*sqrt(_ionEnergy_lab*_ionEnergy_lab - (targetBeamA()*targetBeamA())*protonMass*protonMass);
198  double _electronPz_lab = 1.0*sqrt(_electronEnergy_lab*_electronEnergy_lab - starlightConstants::mel*(starlightConstants::mel));
199  double _totalEnergy_lab = 1.0*_electronEnergy_lab + _ionEnergy_lab;
200  double _totalPz_lab = _ionPz_lab + _electronPz_lab;
201  _rap_CM = (1.0/2.0)*log((_totalEnergy_lab + _totalPz_lab)/(_totalEnergy_lab - _totalPz_lab));
202  double _totalEnergy_COM;
203  _totalEnergy_COM = sqrt((_totalEnergy_lab) * (_totalEnergy_lab) - (_totalPz_lab) * (_totalPz_lab)); //This is Lorentz invariant.
204 
205  _beamLorentzGamma = cosh(_rap_CM-rap2);
206  _targetLorentzGamma = cosh(rap1-rap2);
207 
208  if( targetBeamA() == 1) //proton case 0.87 fm = 4.4 GeV^{-1}
209  _targetR = 4.4;
210  else
211  _targetR = 6.1 * pow(targetBeamA(), 1./3.);
212  //_targetR = 4.4/2.;
213 
214  _fixedQ2Range = false;
215  std::cout << "Rapidity electron beam: " << rap1 << ", rapidity target beam: " << rap2 << ", rapidity CMS system: " << _rap_CM << ", beam gamma in CMS: " << _beamLorentzGamma<< std::endl;
216  std::cout << "Rapidity beam 1 in beam 2 frame: " << rap1-rap2 << ", beam 1 gamma in beam 2 frame: " << _targetLorentzGamma<< std::endl;
217  _ptBinWidthInterference = maxPtInterference() / nmbPtBinsInterference();
218  _protonEnergy = _beamLorentzGamma * protonMass;
219  // Electron energy in the target frame of reference
220  _electronEnergy = _targetLorentzGamma * starlightConstants::mel;
221 
222  if( (targetBeamZ()==1) && (targetBeamA()==3) ){
223  printWarn << "tritium is not currently supported" << endl;
224  return false;}
225  // check that rho production uses wide resonance option
226  if(_prodParticleId.value()==113 && _productionMode.value()==2){
227  printWarn << endl<< "For rho meson production, you should choose the wide resonance option (production mode = 3)" << endl;
228  return false;}
229 
230  // define interaction type
231  switch (productionMode()) {
232  case 1:
233  _interactionType = PHOTONPHOTON;
234  break;
235  case 2:
236  _interactionType = PHOTONPOMERONNARROW;
237  break;
238  case 3:
239  _interactionType = PHOTONPOMERONWIDE;
240  break;
241  case 4:
242  _interactionType = PHOTONPOMERONINCOHERENT;
243  break;
244  case 5:
245  _interactionType = PHOTONUCLEARSINGLE;
246  break;
247  case 6:
248  _interactionType = PHOTONUCLEARDOUBLE;
249  break;
250  case 7:
251  _interactionType = PHOTONUCLEARSINGLEPA;
252  break;
253  case 8:
254  _interactionType = PHOTONUCLEARSINGLEPAPY;
255  break;
256 // case 9:
257 // _interactionType = PHOTONPHOTONKNIEHL;
258 // break;
259 // case 10:
260 // _interactionType = PHOTONPHOTONKNIEHLMODIFIED;
261 // break;
262  case 12:
263  _interactionType = E_PHOTONPOMERONNARROW;
264  break;
265  case 13:
266  _interactionType = E_PHOTONPOMERONWIDE;
267  break;
268 
269  default:
270  printWarn << "unknown production mode '" << _productionMode << "'" << endl;
271  return false;
272  }
273  if( (_productionMode.value() ==4) && (_interferenceEnabled.value())) {
274  printWarn << " cannot enable interference for incoherent production " << endl;
275  return false;
276  }
277 
278  double mass = 0;
279  double width = 0;
280  double defaultMinW = 0; // default for _minW, unless it is defined later [GeV/c^2]
281  double defaultMaxW = 0; // default for _maxW, unless it is defined later [GeV/c^2]
282  switch (prodParticleId()) {
283 
284 // case 24: // W+W- pair
285 // _particleType = W;
286 // _decayType = WW;
287 // defaultMinW = 2 * muonMass;
288 // break;
289  case 115: // a_2(1320)
290  _particleType = A2;
291  _decayType = SINGLEMESON;
292  mass = starlightConstants::a2Mass;
293  width = starlightConstants::a2Width;
294  defaultMinW = mass - 5*width; // JES 6.17.2015 to avoid problems with default of 0
295  defaultMaxW = mass + 5*width; // JES 6.17.2015 to avoid problems with no default
296  _inputBranchingRatio = 1.0;
297  break;
298  case 221: // eta
299  _particleType = ETA;
300  _decayType = SINGLEMESON;
301  mass = starlightConstants::etaMass;
302  width = starlightConstants::etaWidth;
303  defaultMinW = mass - 5*width; // JES 6.17.2015 to avoid problems with default of 0
304  defaultMaxW = mass + 5*width; // JES 6.17.2015 to avoid problems with no default
305  _inputBranchingRatio = 1.0;
306  break;
307  case 225: // f_2(1270)
308  _particleType = F2;
309  _decayType = SINGLEMESON;
310  mass = starlightConstants::f2Mass;
311  width = starlightConstants::f2Width;
312  defaultMinW = mass - 5*width; // JES 6.17.2015 to avoid problems with default of 0
313  defaultMaxW = mass + 5*width; // JES 6.17.2015 to avoid problems with no default
314  _inputBranchingRatio = starlightConstants::f2BrPiPi;
315  break;
316  case 331: // eta'(958)
317  _particleType = ETAPRIME;
318  _decayType = SINGLEMESON;
319  mass = starlightConstants::etaPrimeMass;
320  width = starlightConstants::etaPrimeWidth;
321  defaultMinW = mass - 5*width; // JES 6.17.2015 to avoid problems with default of 0
322  defaultMaxW = mass + 5*width; // JES 6.17.2015 to avoid problems with no default
323  _inputBranchingRatio = 1.0;
324  break;
325  case 335: // f_2'(1525)
326  _particleType = F2PRIME;
327  _decayType = SINGLEMESON;
328  mass = starlightConstants::f2PrimeMass;
329  width = starlightConstants::f2PrimeWidth;
330  defaultMinW = mass - 5*width; // JES 6.17.2015 to avoid problems with default of 0
331  defaultMaxW = mass + 5*width; // JES 6.17.2015 to avoid problems with no default
332  _inputBranchingRatio = starlightConstants::f2PrimeBrKK;
333  break;
334  case 441: // eta_c(1s)
335  _particleType = ETAC;
336  _decayType = SINGLEMESON;
337  mass = starlightConstants::etaCMass;
338  width = starlightConstants::etaCWidth;
339  defaultMinW = mass - 5*width; // JES 6.17.2015 to avoid problems with default of 0
340  defaultMaxW = mass + 5*width; // JES 6.17.2015 to avoid problems with no default
341  _inputBranchingRatio = 1.0;
342  break;
343  case 9010221: // f_0(980), was orginally called 10221? updated to standard number
344  _particleType = F0;
345  _decayType = SINGLEMESON;
346  mass = starlightConstants::f0Mass;
347  width = starlightConstants::f0Width;
348  defaultMinW = mass - 5*width; // JES 6.17.2015 to avoid problems with default of 0
349  defaultMaxW = mass + 5*width; // JES 6.17.2015 to avoid problems with no default
350  _inputBranchingRatio = starlightConstants::f0BrPiPi;
351  break;
352  case 33: // Z"/Z0 This is the rho^0 rho^0 final state SRK
353  _particleType = ZOVERZ03;
354  _decayType = SINGLEMESON;
355  defaultMinW = 4*pionChargedMass;
356  defaultMaxW = 1.6; // JES 6.17.2015 to avoid problems with no default
357  _inputBranchingRatio = 1.0;
358  break;
359  case 88: // axion// AXION HACK, till break statement
360  _particleType = AXION;
361  _decayType = SINGLEMESON;
362  mass = _axionMass.value();
363  width = 1/(64*starlightConstants::pi)*mass*mass*mass/(1000*1000);//Fix Lambda=1000 GeV, rescaling is trivial.
364  defaultMinW = mass - 5*width; // JES 6.17.2015 to avoid problems with default of 0
365  defaultMaxW = mass + 5*width; // JES 6.17.2015 to avoid problems with no default
366  break; // AXION HACK, end
367 // case 25: // Higgs
368 // _particleType = HIGGS;
369 // _decayType = SINGLEMESON;
370 // break;
371  case 113: // rho(770)
372  _particleType = RHO;
373  _decayType = WIDEVMDEFAULT;
374  mass = starlightConstants::rho0Mass;
375  width = starlightConstants::rho0Width;
376  defaultMinW = 2 * pionChargedMass;
377  defaultMaxW = mass + 5 * width;
378  _inputBranchingRatio = starlightConstants::rho0BrPiPi;
379  break;
380  case 913: // rho(770) with direct pi+pi- decay, interference given by ZEUS data
381  _particleType = RHOZEUS;
382  _decayType = WIDEVMDEFAULT;
383  mass = starlightConstants::rho0Mass;
384  width = starlightConstants::rho0Width;
385  defaultMinW = 2 * pionChargedMass;
386  defaultMaxW = mass + 5 * width; // use the same 1.5GeV max mass as ZEUS
387  _inputBranchingRatio = starlightConstants::rho0BrPiPi;
388  break;
389  // case 999: // pi+pi-pi+pi- phase space decay
390  // _particleType = FOURPRONG;
391  // _decayType = WIDEVMDEFAULT;
392  // mass = starlightConstants::rho0PrimeMass;
393  // width = starlightConstants::rho0PrimeWidth;
394  // defaultMinW = 4 * pionChargedMass;
395  // defaultMaxW = sqrt(beam1LorentzGamma()*beam2LorentzGamma())*2*(starlightConstants::hbarc)/(1.2*pow(float(beam1A()),1./6.)*pow(float(beam2A()),1./6.)); // JES 6.17.2015 to avoid problems with no default
396  // _inputBranchingRatio = 1.0;
397  // break;
398  case 223: // omega(782)
399  _particleType = OMEGA;
400  _decayType = NARROWVMDEFAULT;
401  mass = starlightConstants::OmegaMass;
402  width = starlightConstants::OmegaWidth;
403  defaultMinW = mass - 5 * width;
404  defaultMaxW = mass + 5 * width;
405  _inputBranchingRatio = starlightConstants::OmegaBrPiPi;
406  break;
407  case 223022: // omega(782) -> pi0 + gamma
408  _particleType = OMEGA_pi0gamma;
409  _decayType = NARROWVMDEFAULT;
410  mass = starlightConstants::OmegaMass;
411  width = starlightConstants::OmegaWidth;
412  defaultMinW = mass - 5 * width;
413  defaultMaxW = mass + 5 * width;
414  _inputBranchingRatio = starlightConstants::OmegaBrPi0Gamma;
415  break;
416  case 333: // phi(1020)
417  _particleType = PHI;
418  _decayType = NARROWVMDEFAULT;
419  mass = starlightConstants::PhiMass;
420  width = starlightConstants::PhiWidth;
421  defaultMinW = 2 * kaonChargedMass;
422  defaultMaxW = mass + 5 * width;
423  _inputBranchingRatio = starlightConstants::PhiBrKK;
424  break;
425  case 443: // J/psi
426  _particleType = JPSI;
427  _decayType = NARROWVMDEFAULT;
428  mass = starlightConstants::JpsiMass;
429  width = starlightConstants::JpsiWidth;
430  defaultMinW = mass - 5 * width;
431  defaultMaxW = mass + 5 * width;
432  _inputBranchingRatio = (starlightConstants::JpsiBree + starlightConstants::JpsiBrmumu)/2.;
433  break;
434  case 443011: // J/psi
435  _particleType = JPSI_ee;
436  _decayType = NARROWVMDEFAULT;
437  mass = starlightConstants::JpsiMass;
438  width = starlightConstants::JpsiWidth;
439  defaultMinW = mass - 5 * width;
440  defaultMaxW = mass + 5 * width;
441  _inputBranchingRatio = starlightConstants::JpsiBree;
442  break;
443  case 443013: // J/psi
444  cout<<"In inputParameters setting J/psi mass!"<<endl;
445  _particleType = JPSI_mumu;
446  _decayType = NARROWVMDEFAULT;
447  mass = starlightConstants::JpsiMass;
448  width = starlightConstants::JpsiWidth;
449  defaultMinW = mass - 5 * width;
450  defaultMaxW = mass + 5 * width;
451  _inputBranchingRatio = starlightConstants::JpsiBrmumu;
452  break;
453  case 444: // psi(2S)
454  _particleType = JPSI2S;
455  _decayType = NARROWVMDEFAULT;
456  mass = starlightConstants::Psi2SMass;
457  width = starlightConstants::Psi2SWidth;
458  defaultMinW = mass - 5 * width;
459  defaultMaxW = mass + 5 * width;
460  _inputBranchingRatio = (starlightConstants::Psi2SBree + starlightConstants::Psi2SBrmumu)/2.;
461  break;
462  case 444011: // psi(2S)
463  _particleType = JPSI2S_ee;
464  _decayType = NARROWVMDEFAULT;
465  mass = starlightConstants::Psi2SMass;
466  width = starlightConstants::Psi2SWidth;
467  defaultMinW = mass - 5 * width;
468  defaultMaxW = mass + 5 * width;
469  _inputBranchingRatio = starlightConstants::Psi2SBree;
470  break;
471  case 444013: // psi(2S)
472  _particleType = JPSI2S_mumu;
473  _decayType = NARROWVMDEFAULT;
474  mass = starlightConstants::Psi2SMass;
475  width = starlightConstants::Psi2SWidth;
476  defaultMinW = mass - 5 * width;
477  defaultMaxW = mass + 5 * width;
478  _inputBranchingRatio = starlightConstants::Psi2SBrmumu;
479  break;
480  case 553: // Upsilon(1S)
481  _particleType = UPSILON;
482  _decayType = NARROWVMDEFAULT;
483  mass = starlightConstants::Upsilon1SMass;
484  width = starlightConstants::Upsilon1SWidth;
485  defaultMinW = mass - 5 * width;
486  defaultMaxW = mass + 5 * width;
487  _inputBranchingRatio = (starlightConstants::Upsilon1SBree + starlightConstants::Upsilon1SBrmumu)/2.;
488  break;
489  case 553011: // Upsilon
490  _particleType = UPSILON_ee;
491  _decayType = NARROWVMDEFAULT;
492  mass = starlightConstants::Upsilon1SMass;
493  width = starlightConstants::Upsilon1SWidth;
494  defaultMinW = mass - 5 * width;
495  defaultMaxW = mass + 5 * width;
496  _inputBranchingRatio = starlightConstants::Upsilon1SBree;
497  break;
498  case 553013: // Upsilon
499  _particleType = UPSILON_mumu;
500  _decayType = NARROWVMDEFAULT;
501  mass = starlightConstants::Upsilon1SMass;
502  width = starlightConstants::Upsilon1SWidth;
503  defaultMinW = mass - 5 * width;
504  defaultMaxW = mass + 5 * width;
505  _inputBranchingRatio = starlightConstants::Upsilon1SBrmumu;
506  break;
507  case 554: // Upsilon(2S)
508  _particleType = UPSILON2S;
509  _decayType = NARROWVMDEFAULT;
510  mass = starlightConstants::Upsilon2SMass;
511  width = starlightConstants::Upsilon2SWidth;
512  defaultMinW = mass - 5 * width;
513  defaultMaxW = mass + 5 * width;
514  _inputBranchingRatio = (starlightConstants::Upsilon2SBree + starlightConstants::Upsilon2SBrmumu)/2.;
515  break;
516  case 554011: // Upsilon(2S)
517  _particleType = UPSILON2S_ee;
518  _decayType = NARROWVMDEFAULT;
519  mass = starlightConstants::Upsilon2SMass;
520  width = starlightConstants::Upsilon2SWidth;
521  defaultMinW = mass - 5 * width;
522  defaultMaxW = mass + 5 * width;
523  _inputBranchingRatio = starlightConstants::Upsilon2SBree;
524  break;
525  case 554013: // Upsilon(2S)
526  _particleType = UPSILON2S_mumu;
527  _decayType = NARROWVMDEFAULT;
528  mass = starlightConstants::Upsilon2SMass;
529  width = starlightConstants::Upsilon2SWidth;
530  defaultMinW = mass - 5 * width;
531  defaultMaxW = mass + 5 * width;
532  _inputBranchingRatio = starlightConstants::Upsilon2SBrmumu;
533  break;
534  case 555: // Upsilon(3S)
535  _particleType = UPSILON3S;
536  _decayType = NARROWVMDEFAULT;
537  mass = starlightConstants::Upsilon3SMass;
538  width = starlightConstants::Upsilon3SWidth;
539  defaultMinW = mass - 5 * width;
540  defaultMaxW = mass + 5 * width;
541  _inputBranchingRatio = (starlightConstants::Upsilon3SBree + starlightConstants::Upsilon3SBrmumu)/2.;
542  break;
543  case 555011: // Upsilon(3S)
544  _particleType = UPSILON3S_ee;
545  _decayType = NARROWVMDEFAULT;
546  mass = starlightConstants::Upsilon3SMass;
547  width = starlightConstants::Upsilon3SWidth;
548  defaultMinW = mass - 5 * width;
549  defaultMaxW = mass + 5 * width;
550  _inputBranchingRatio = starlightConstants::Upsilon3SBree;
551  break;
552  case 555013: // Upsilon(3S)
553  _particleType = UPSILON3S_mumu;
554  _decayType = NARROWVMDEFAULT;
555  mass = starlightConstants::Upsilon3SMass;
556  width = starlightConstants::Upsilon3SWidth;
557  defaultMinW = mass - 5 * width;
558  defaultMaxW = mass + 5 * width;
559  _inputBranchingRatio = starlightConstants::Upsilon3SBrmumu;
560  break;
561  default:
562  printWarn << "unknown particle ID " << _prodParticleId << endl;
563  return false;
564  } // _prodParticleId
565 
566  if (_minW.value() == -1)
567  _minW = defaultMinW;
568  if (_maxW.value() == -1)
569  _maxW = defaultMaxW;
570  if ( _maxW.value() <= _minW.value() ) {
571  printWarn << "maxW must be greater than minW" << endl;
572  return false;
573  }
574 
575  if (_minW_GP.value() == -1)
576  _minW_GP = 0;
577  if (_maxW_GP.value() == -1)
578  _maxW_GP = 1e7;
579  if ( _maxW_GP.value() <= _minW_GP.value() ) {
580  printWarn << "maxW_GA must be greater than minW_GP" << endl;
581  printWarn <<"The value of minW_GP is " << _minW_GP << endl;
582  printWarn <<"The value of maxW_GP is " << _maxW_GP << endl;
583  return false;
584  }
585  if(_minW_GP.value() > _totalEnergy_COM) {
586  /* If the minimum COM energy of gamma and nucleon is set by user
587  to be greater than the total energy supply by electron and target, code will terminate. */
588  printWarn << "ERROR: The current input minimum CoM energy (_minW_GP) is: "<< _minW_GP
589  << " which is larger than the total center-of-mass energy of the electron-ion system: "
590  <<_totalEnergy_COM << endl;
591  cout<<"Exiting now..." << endl;
592  return false;
593  }
594 
595  // Sanity check on Q2 range in case it is specified by user
596  if( _minGammaQ2.value() != 0 || _maxGammaQ2.value() != 0){
597  if( _minGammaQ2.value() <0 || _maxGammaQ2.value() <=_minGammaQ2.value() )
598  printWarn << "Input values for min and max Q2 are inconsistent: "<<_minGammaQ2.value()<<","<<_maxGammaQ2.value()
599  <<". Continuing with default "<<endl;
600  else
601  _fixedQ2Range = true;
602  }
603  //Photon energy limits in C.M.S frame, used for some basic safety chekcs
604  _cmsMinPhotonEnergy = 0.5*(((mass+protonMass)*(mass+protonMass)-
605  protonMass*protonMass)/(_protonEnergy.value()+sqrt(_protonEnergy.value()*_protonEnergy.value()-protonMass*protonMass)));
606  _cmsMaxPhotonEnergy = 0.5*mass*exp(9);
607  // Photon limits in target frame: this is where the photon flux is well described
608  _targetMaxPhotonEnergy = (_targetLorentzGamma - 10. ) *starlightConstants::mel;
609  _targetMinPhotonEnergy = mass;
610  printInfo << "using the following " << *this;
611 
612  return true;
613 }
614 
615 
616 //______________________________________________________________________________
617 ostream&
618 inputParameters::print(ostream& out) const
619 {
620  out << "starlight parameters:" << endl
621  << " base file name ...................... '" << _baseFileName.value() << "'" << endl
622  << " target beam atomic number .............. " << _targetBeamZ.value() << endl
623  << " target beam atomic mass number ......... " << _targetBeamA.value() << endl
624  << " Lorentz gamma of beams in CM frame ..... " << _beamLorentzGamma << endl
625  << " mass W of produced hadronic system ..... " << _minW.value() << " < W < " << _maxW.value() << " GeV/c^2" << endl
626  << " # of W bins ............................ " << _nmbWBins.value() << endl
627  << " maximum absolute value for rapidity .... " << _maxRapidity.value() << endl
628  << " # of rapidity bins ..................... " << _nmbRapidityBins.value() << endl
629  << " # of Egamma bins ....................... " << _nmbEnergyBins.value() << endl
630  << " cut in pT............................... " << yesNo(_ptCutEnabled.value()) << endl;
631  if (_ptCutEnabled.value()) {
632  out << " minumum pT.......................... " << _ptCutMin.value() << " GeV/c" << endl
633  << " maximum pT.......................... " << _ptCutMax.value() << " GeV/c" << endl;}
634  out << " cut in eta.............................. " << yesNo(_etaCutEnabled.value()) << endl;
635  if (_etaCutEnabled.value()) {
636  out << " minumum eta......................... " << _etaCutMin.value() << endl
637  << " maximum eta......................... " << _etaCutMax.value() << endl;}
638  out << " production mode ........................ " << _productionMode.value() << endl
639  << " number of events to generate ........... " << _nmbEventsTot.value() << endl
640  << " PDG ID of produced particle ............ " << _prodParticleId.value() << endl
641  << " seed for random generator .............. " << _randomSeed.value() << endl
642  << " breakup mode for beam particles ........ " << _beamBreakupMode.value() << endl
643  << " interference enabled ................... " << yesNo(_interferenceEnabled.value()) << endl;
644  if (_interferenceEnabled.value()) {
645  out << " interference strength .................. " << _interferenceStrength.value() << endl
646  << " maximum p_T for interference calc. ..... " << _maxPtInterference.value() << " GeV/c" << endl
647  << " # of p_T bins for interference calc. ... " << _nmbPtBinsInterference.value() << endl;}
648  if (_productionMode.value()!=1) {
649  if (_productionMode.value()==4) {
650  out << " coherent scattering off nucleus ........ no" << endl;}
651  else {
652  out << " coherent scattering off nucleus ........ yes" << endl;
653  }
654  }
655  if (_fixedQ2Range == true ){
656  out << " fixed photon Q2 range .................. " << _minGammaQ2.value() << " < Q2 < "
657  <<_maxGammaQ2.value() << " GeV/c^2 "<<endl;
658  }
659  out<< " Q2_BINS" <<nmbGammaQ2Bins () <<endl;
660  out <<" Quantum Glauber parameter............... "<<_quantumGlauber.value()<<endl;
661  out <<" Impulse VM parameter.................... "<<_impulseVM.value()<<endl;
662  return out;
663 }
664 
665 
666 //______________________________________________________________________________
667 ostream&
668 inputParameters::write(ostream& out) const
669 {
670 
671  out << "baseFileName" << baseFileName () <<endl
672  << "TARGET_BEAM_Z" << targetBeamZ () <<endl
673  << "TARGET_BEAM_A" << targetBeamA () <<endl
674  << "BEAM_GAMMA" << beamLorentzGamma () <<endl
675  << "W_MAX" << maxW () <<endl
676  << "W_MIN" << minW () <<endl
677  << "W_GP_MAX" << maxW_GP () <<endl
678  << "W_GP_MIN" << minW_GP () <<endl
679  << "W_N_BINS" << nmbWBins () <<endl
680  << "RAP_MAX" << maxRapidity () <<endl
681  << "RAP_N_BINS" << nmbRapidityBins () <<endl
682  << "CUT_PT" << ptCutEnabled () <<endl
683  << "PT_MIN" << ptCutMin () <<endl
684  << "PT_MAX" << ptCutMax () <<endl
685  << "CUT_ETA" << etaCutEnabled () <<endl
686  << "ETA_MIN" << etaCutMin () <<endl
687  << "ETA_MAX" << etaCutMax () <<endl
688  << "PROD_MODE" << productionMode () <<endl
689  << "N_EVENTS" << nmbEvents () <<endl
690  << "PROD_PID" << prodParticleId () <<endl
691  << "RND_SEED" << randomSeed () <<endl
692  << "BREAKUP_MODE" << beamBreakupMode () <<endl
693  << "INTERFERENCE" << interferenceEnabled () <<endl
694  << "IF_STRENGTH" << interferenceStrength () <<endl
695  << "INT_PT_MAX" << maxPtInterference () <<endl
696  << "INT_PT_N_BINS" << nmbPtBinsInterference() <<endl;
697  out<< "USING FIXED RANGE"<<_fixedQ2Range<<endl;
698  if( _fixedQ2Range == true ){
699  out << "Q2_MIN" << minGammaQ2 () <<endl
700  << "Q2_MAX" << maxGammaQ2 () <<endl;
701  }
702  out<< " Q2_BINS" <<nmbGammaQ2Bins () <<endl;
703  return out;
704 }
705 
706 std::string
707 inputParameters::parameterValueKey() const
708 {
709  return parameterbase::_parameters.validationKey();
710 }