[u/mrichter/AliRoot.git] / STARLIGHT / starlight / src / inputParameters.cpp

///////////////////////////////////////////////////////////////////////////
//
//    Copyright 2010
//
//    This file is part of starlight.
//
//    starlight is free software: you can redistribute it and/or modify
//    it under the terms of the GNU General Public License as published by
//    the Free Software Foundation, either version 3 of the License, or
//    (at your option) any later version.
//
//    starlight is distributed in the hope that it will be useful,
//    but WITHOUT ANY WARRANTY; without even the implied warranty of
//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
//    GNU General Public License for more details.
//
//    You should have received a copy of the GNU General Public License
//    along with starlight. If not, see <http://www.gnu.org/licenses/>.
//
///////////////////////////////////////////////////////////////////////////
//
// File and Version Information:
// $Rev:: 182                         $: revision of last commit
// $Author:: jnystrand                $: author of last commit
// $Date:: 2014-09-12 00:58:56 +0200 #$: date of last commit
//
// Description:
//
//
//
///////////////////////////////////////////////////////////////////////////


#include <iostream>
#include <fstream>

#include "reportingUtils.h"
#include "starlightconstants.h"
#include "inputParameters.h"
#include "inputParser.h"
#include "starlightconfig.h"
#include <cmath>
#include <cstring>
#include "randomgenerator.h"

using namespace std;
using namespace starlightConstants;

parameterlist parameterbase::_parameters;

#define REQUIRED true
#define NOT_REQUIRED false

//______________________________________________________________________________
inputParameters::inputParameters()
	: _configFileName        ("slight.in"),
 	  _beam1Z                ("BEAM_1_Z",0),
	  _beam1A                ("BEAM_1_A",0),
	  _beam2Z                ("BEAM_2_Z",0),
	  _beam2A                ("BEAM_2_A",0),
	  _beam1LorentzGamma     ("BEAM_1_GAMMA",0),
	  _beam2LorentzGamma     ("BEAM_2_GAMMA",0),
	  _maxW                  ("W_MAX",0),
	  _minW                  ("W_MIN",0),
	  _nmbWBins              ("W_N_BINS",0),
	  _maxRapidity           ("RAP_MAX",0),
	  _nmbRapidityBins       ("RAP_N_BINS",0),
	  _ptCutEnabled          ("CUT_PT",false),
	  _ptCutMin              ("PT_MIN",0),
	  _ptCutMax              ("PT_MAX",0),
	  _etaCutEnabled         ("CUT_ETA",false),
	  _etaCutMin             ("ETA_MIN",0),
	  _etaCutMax             ("ETA_MAX",0),
	  _productionMode        ("PROD_MODE",0),
	  _nmbEventsTot          ("N_EVENTS",0),
	  _prodParticleId        ("PROD_PID",0),
	  _randomSeed            ("RND_SEED",0),
	  _beamBreakupMode       ("BREAKUP_MODE",0),
	  _interferenceEnabled   ("INTERFERENCE",false),
	  _interferenceStrength  ("IF_STRENGTH",0),
	  _maxPtInterference     ("INT_PT_MAX",0),
	  _nmbPtBinsInterference ("INT_PT_N_BINS",0),
	  _ptBinWidthInterference("INT_PT_WIDTH",0),
	  _coherentProduction    ("COHERENT",false),
	  _incoherentFactor      ("INCO_FACTOR",0),
	  _protonEnergy          ("PROTON_ENERGY",0, NOT_REQUIRED),
	  _minGammaEnergy	 ("MIN_GAMMA_ENERGY",0, NOT_REQUIRED),
	  _maxGammaEnergy	 ("MAX_GAMMA_ENERGY",0, NOT_REQUIRED),
	  _pythiaParams          ("PYTHIA_PARAMS","", NOT_REQUIRED),
	  _pythiaFullEventRecord ("PYTHIA_FULL_EVENTRECORD",false, NOT_REQUIRED),
	  _xsecCalcMethod	 ("XSEC_METHOD",0, NOT_REQUIRED)
{
  // All parameters must be initialised in initialisation list! 
  // If not: error: 'parameter<T, validate>::parameter() [with T = unsigned int, bool validate = true]' is private
  // or similar
	
	_ip.addParameter(_beam1Z);
	_ip.addParameter(_beam2Z);
	_ip.addParameter(_beam1A);
	_ip.addParameter(_beam2A);

	_ip.addParameter(_beam1LorentzGamma);
	_ip.addParameter(_beam2LorentzGamma);
	
	_ip.addParameter(_maxW);
	_ip.addParameter(_minW);
	
	_ip.addParameter(_nmbWBins);

	_ip.addParameter(_maxRapidity);
	_ip.addParameter(_nmbRapidityBins);
	
	_ip.addParameter(_ptCutEnabled);
	_ip.addParameter(_ptCutMin);
	_ip.addParameter(_ptCutMax);
	
	_ip.addParameter(_etaCutEnabled);
	_ip.addParameter(_etaCutMax);
	_ip.addParameter(_etaCutMin);
	
	_ip.addParameter(_productionMode);
	_ip.addParameter(_nmbEventsTot);
	_ip.addParameter(_prodParticleId);
	_ip.addParameter(_randomSeed);
	_ip.addParameter(_beamBreakupMode);
	_ip.addParameter(_interferenceEnabled);
	_ip.addParameter(_interferenceStrength);
	_ip.addParameter(_maxPtInterference);
	_ip.addParameter(_nmbPtBinsInterference);
	_ip.addParameter(_coherentProduction);
	_ip.addParameter(_incoherentFactor);
	_ip.addParameter(_minGammaEnergy);
	_ip.addParameter(_maxGammaEnergy);
	_ip.addParameter(_pythiaParams);
	_ip.addParameter(_pythiaFullEventRecord);
	_ip.addParameter(_xsecCalcMethod);
}


//______________________________________________________________________________
inputParameters::~inputParameters()
{ }


//______________________________________________________________________________
bool
inputParameters::configureFromFile(const std::string &configFileName)
{
	// open config file
	_configFileName = configFileName;
	
	int nParameters = _ip.parseFile(_configFileName);
	
	if(nParameters == -1) 
	{
		printWarn << "could not open file '" << _configFileName << "'" << endl;
	  return false;
	}
	
	//ip.printParameterInfo(cout);
	
	if(_ip.validateParameters(cerr))
		printInfo << "successfully read input parameters from '" << _configFileName << "'" << endl;
	else {
 		printWarn << "problems reading input parameters from '" << _configFileName << "'" << endl
 		          << *this;
 		return false;
 	}
 	return true;
}
 bool inputParameters::init()
 {
	// Set the seed for the random generator
	randyInstance.SetSeed(_randomSeed.value());
	
 	// Calculate beam gamma in CMS frame
 	double rap1 = acosh(beam1LorentzGamma());
	double rap2 = -acosh(beam2LorentzGamma());
	_beamLorentzGamma = cosh((rap1-rap2)/2);
	
	std::cout << "Rapidity beam 1: " << rap1 << ", rapidity beam 2: " << rap2 << ", rapidity CMS system: " << (rap1+rap2)/2 << ", beam gamma in CMS: " << _beamLorentzGamma<< std::endl;
	_ptBinWidthInterference = maxPtInterference() / nmbPtBinsInterference();
	_protonEnergy           = _beamLorentzGamma * protonMass;

	// define interaction type
	switch (productionMode()) {
	case 1:
		_interactionType = PHOTONPHOTON;
		break;
	case 2:
		_interactionType = PHOTONPOMERONNARROW;
		break;
	case 3:
		_interactionType = PHOTONPOMERONWIDE;
		break;
        case 4:
                _interactionType = PHOTONPOMERONINCOHERENT;
                break;
	case 5:
		_interactionType = PHOTONUCLEARSINGLE;
		break;
	case 6:
		_interactionType = PHOTONUCLEARDOUBLE;
		break;
	case 7:
		_interactionType = PHOTONUCLEARSINGLEPA;
		break;
	case 8:
		_interactionType = PHOTONUCLEARSINGLEPAPY;
		break;
// 	case 9:
// 		_interactionType = PHOTONPHOTONKNIEHL;
// 		break;
//         case 10:
//                 _interactionType = PHOTONPHOTONKNIEHLMODIFIED;
//                 break;
	default:
		printWarn << "unknown production mode '" << _productionMode << "'" << endl;
		return false;
	}
  
	//Trying to define the proper Wmins and Wmaxs. a TEMPORARY fix....Better solution=??
	double mass        = 0;
	double width       = 0;
	double defaultMinW = 0;  // default for _minW, unless it is defined later [GeV/c^2]
	switch (prodParticleId()) {
	case 11:  // e+e- pair
		_particleType = ELECTRON;
		_decayType    = LEPTONPAIR;
		defaultMinW   = 0.01; // default is 0.01; up to 0.15 is safe for Summer 2000 triggering for e+e- pairs
		break;
	case 13:  // mu+mu- pair
		_particleType = MUON;
		_decayType    = LEPTONPAIR;
		defaultMinW   = 2 * muonMass;
		break;
	case 15:  // tau+tau- pair
		_particleType = TAUON;
		_decayType    = LEPTONPAIR;
		defaultMinW   = 2 * tauMass;
		break;
	case 10015:  // tau+tau- pair
		_particleType = TAUONDECAY;
		_decayType    = LEPTONPAIR;
		defaultMinW   = 2 * tauMass;
		break;
// 	case 24:  // W+W- pair
// 		_particleType = W;
// 		_decayType    = WW;
// 		defaultMinW   = 2 * muonMass;
// 		break;
	case 115:  // a_2(1320)
		_particleType = A2;
		_decayType    = SINGLEMESON;
		break;
	case 221:  // eta
		_particleType = ETA;
		_decayType    = SINGLEMESON;
		break;
	case 225:  // f_2(1270)
		_particleType = F2;
                defaultMinW   = 2*pionChargedMass;
		_decayType    = SINGLEMESON;
		break;
	case 331:  // eta'(958)
		_particleType = ETAPRIME;
		_decayType    = SINGLEMESON;
		break;
	case 335:  // f_2'(1525)
		_particleType = F2PRIME;
		_decayType    = SINGLEMESON;
		break;
	case 441:  // eta_c(1s)
		_particleType = ETAC;
		_decayType    = SINGLEMESON;
                defaultMinW   = etaCMass - 5 * 0.0267;
		break;
	case 9010221:  // f_0(980), was orginally called 10221? updated to standard number
		_particleType = F0;
		_decayType    = SINGLEMESON;
                defaultMinW   = 2*pionNeutralMass;
		break;
	case 33:  // Z"/Z03
		_particleType = ZOVERZ03;
		_decayType    = SINGLEMESON;
		break;
// 	case 25: // Higgs
// 		_particleType = HIGGS;
// 		_decayType    = SINGLEMESON;
// 		break;
	case 113:  // rho(770)
		_particleType = RHO;
		_decayType    = WIDEVMDEFAULT;
		mass          = 0.7685;
		width         = 0.1507;
		defaultMinW   = 2 * pionChargedMass;
		_maxW         = mass + 5 * width;
		break;
	case 913:  // rho(770) with direct pi+pi- decay, interference given by ZEUS data
		_particleType = RHOZEUS;
		_decayType    = WIDEVMDEFAULT;
		mass          = 0.7685;                    
		width         = 0.1507;
		defaultMinW   = 2 * pionChargedMass;
		_maxW         = mass + 5 * width;  // use the same 1.5GeV max mass as ZEUS
		break;
	case 999:  // pi+pi-pi+pi- phase space decay
		_particleType = FOURPRONG;
		_decayType    = WIDEVMDEFAULT;
		mass          = 1.350;
		width         = 0.360;
		defaultMinW   = 4 * pionChargedMass;
		_maxW         = 3;
		break;
	case 223:  // omega(782)
		_particleType = OMEGA;
		_decayType    = NARROWVMDEFAULT;  // will probably be moved to 3-body decay
		mass          = 0.78194;
		width         = 0.00843;
		defaultMinW   = mass - 5 * width;
		_maxW         = mass + 5 * width;
		break;
	case 333:  // phi(1020)
		_particleType = PHI;
		_decayType    = NARROWVMDEFAULT;
		mass          = 1.019413;
		width         = 0.00443;
		defaultMinW   = 2 * kaonChargedMass;
		_maxW         = mass + 5 * width;
		break;
	case 443:  // J/psi
		_particleType = JPSI;
		_decayType    = NARROWVMDEFAULT;
		mass          = 3.09692;   // JN  3.09688;
		width         = 0.000091;  // JN  0.000087;
		defaultMinW   = mass - 5 * width;
		_maxW         = mass + 5 * width;
		break;
	case 443011:  // J/psi
		_particleType = JPSI_ee;
		_decayType    = NARROWVMDEFAULT;
		mass          = 3.09692;   // JN  3.09688;
		width         = 0.000091;  // JN  0.000087;
		defaultMinW   = mass - 5 * width;
		_maxW         = mass + 5 * width;
		break;
	case 443013:  // J/psi
		_particleType = JPSI_mumu;
		_decayType    = NARROWVMDEFAULT;
		mass          = 3.09692;   // JN  3.09688;
		width         = 0.000091;  // JN  0.000087;
		defaultMinW   = mass - 5 * width;
		_maxW         = mass + 5 * width;
		break;
	case 444:  // J/psi
		_particleType = JPSI2S;
		_decayType    = NARROWVMDEFAULT;
		mass          = 3.686093;
		width         = 0.000337;
		defaultMinW   = mass - 5 * width;
		_maxW         = mass + 5 * width;
		break;
	case 444011: // J/psi
		_particleType = JPSI2S_ee;
		_decayType    = NARROWVMDEFAULT;
		mass          = 3.686093;     
		width         = 0.000337;
		defaultMinW   = mass - 5 * width;
		_maxW         = mass + 5 * width;
		break;
	case 444013:  // J/psi
		_particleType = JPSI2S_mumu;
		_decayType    = NARROWVMDEFAULT;
		mass          = 3.686093;
		width         = 0.000337;
		defaultMinW   = mass - 5 * width;
		_maxW         = mass + 5 * width;
		break;
	case 553:  // Upsilon
		_particleType = UPSILON;
		_decayType    = NARROWVMDEFAULT;
		mass          = 9.46030;
		width         = 0.00005402;
		defaultMinW   = mass - 5 * width;
		_maxW         = mass + 5 * width;
		break;
	case 553011:  // Upsilon
		_particleType = UPSILON_ee;
		_decayType    = NARROWVMDEFAULT;
		mass          = 9.46030;
		width         = 0.00005402;
		defaultMinW   = mass - 5 * width;
		_maxW         = mass + 5 * width;
		break;
	case 553013:  // Upsilon
		_particleType = UPSILON_mumu;
		_decayType    = NARROWVMDEFAULT;
		mass          = 9.46030;
		width         = 0.00005402;
		defaultMinW   = mass - 5 * width;
		_maxW         = mass + 5 * width;
		break;
	case 554:  // Upsilon(2S)
		_particleType = UPSILON2S;
		_decayType    = NARROWVMDEFAULT;
		mass          = 10.02326;
		width         = 0.00003198;
		defaultMinW   = mass - 5 * width;
		_maxW         = mass + 5 * width;
		break;
	case 554011:  // Upsilon(2S)
		_particleType = UPSILON2S_ee;
		_decayType    = NARROWVMDEFAULT;
		mass          = 10.02326;
		width         = 0.00003198;
		defaultMinW   = mass - 5 * width;
		_maxW         = mass + 5 * width;
		break;
	case 554013:  // Upsilon(2S)
		_particleType = UPSILON2S_mumu;
		_decayType    = NARROWVMDEFAULT;
		mass          = 10.02326;
		width         = 0.00003198;
		defaultMinW   = mass - 5 * width;
		_maxW         = mass + 5 * width;
		break;
	case 555:  // Upsilon(3S)
		mass          = 10.3552;
		width         = 0.00002032;
		defaultMinW   = mass - 5 * width;
		_maxW         = mass + 5 * width;
		_particleType = UPSILON3S;
		_decayType    = NARROWVMDEFAULT;
		break;
	case 555011:  // Upsilon(3S)
		_particleType = UPSILON3S_ee;
		_decayType    = NARROWVMDEFAULT;
		mass          = 10.3552;
		width         = 0.00002032;
		defaultMinW   = mass - 5 * width;
		_maxW         = mass + 5 * width;
		break;
	case 555013:  // Upsilon(3S)
		_particleType = UPSILON3S_mumu;
		_decayType    = NARROWVMDEFAULT;
		mass          = 10.3552;
		width         = 0.00002032;
		defaultMinW   = mass - 5 * width;
		_maxW         = mass + 5 * width;
		break;
	default:
		printWarn << "unknown particle ID " << _prodParticleId << endl;
		return false;
	}  // _prodParticleId

	if (_minW.value() == -1)
		_minW = defaultMinW;

	printInfo << "using the following " << *this;
	
	return true;
}


//______________________________________________________________________________
ostream&
inputParameters::print(ostream& out) const
{
	out << "starlight parameters:" << endl
	    << "    config file name  ...................... '" << _configFileName << "'" << endl
	    << "    beam 1 atomic number ................... " << _beam1Z.value() << endl
	    << "    beam 1 atomic mass number .............. " << _beam1A.value() << endl
	    << "    beam 2 atomic number ................... " << _beam2Z.value() << endl
	    << "    beam 2 atomic mass number .............. " << _beam2A.value() << endl
	    << "    Lorentz gamma of beams in CM frame ..... " << _beamLorentzGamma << endl
	    << "    mass W of produced hadronic system ..... " << _minW.value() << " < W < " << _maxW.value() << " GeV/c^2" << endl
	    << "    # of W bins ............................ " << _nmbWBins.value() << endl
	    << "    maximum absolute value for rapidity .... " << _maxRapidity.value() << endl
	    << "    # of rapidity bins ..................... " << _nmbRapidityBins.value() << endl
	    << "    cut in pT............................... " << yesNo(_ptCutEnabled.value()) << endl;
    if (_ptCutEnabled.value()) {
	out << "        minumum pT.......................... " << _ptCutMin.value() << " GeV/c" << endl
	    << "        maximum pT.......................... " << _ptCutMax.value() << " GeV/c" << endl;}
	out << "    cut in eta.............................. " << yesNo(_etaCutEnabled.value()) << endl;
    if (_etaCutEnabled.value()) {
	out << "        minumum eta......................... " << _etaCutMin.value() << endl
	    << "        maximum eta......................... " << _etaCutMax.value() << endl;}
        out << "    production mode ........................ " << _productionMode.value() << endl
	    << "    number of events to generate ........... " << _nmbEventsTot.value() << endl
	    << "    PDG ID of produced particle ............ " << _prodParticleId.value() << endl
	    << "    seed for random generator .............. " << _randomSeed.value() << endl
	    << "    breakup mode for beam particles ........ " << _beamBreakupMode.value() << endl
	    << "    interference enabled ................... " << yesNo(_interferenceEnabled.value()) << endl;
    if (_interferenceEnabled.value()) {
	out << "    interference strength .................. " << _interferenceStrength.value() << endl
	    << "    maximum p_T for interference calc. ..... " << _maxPtInterference.value() << " GeV/c" << endl
	    << "    # of p_T bins for interference calc. ... " << _nmbPtBinsInterference.value() << endl;}
    if (_productionMode.value()!=1) {
	out  << "    coherent scattering off nucleus ........ " << yesNo(_coherentProduction.value()) << endl;
    	if (!_coherentProduction.value()) {
	 out << "    scaling factor for incoh. VM prod. ..... " << _incoherentFactor.value() << endl;}
	}
	return out;
}


//______________________________________________________________________________
ostream&
inputParameters::write(ostream& out) const
{
  
	out << "BEAM_1_Z"      << beam1Z               () <<endl
	    << "BEAM_2_Z"      << beam1A               () <<endl
	    << "BEAM_1_A"      << beam2Z               () <<endl
	    << "BEAM_2_A"      << beam2A               () <<endl
	    << "BEAM_GAMMA"    << beamLorentzGamma     () <<endl
	    << "W_MAX"         << maxW                 () <<endl
	    << "W_MIN"         << minW                 () <<endl
	    << "W_N_BINS"      << nmbWBins             () <<endl
	    << "RAP_MAX"       << maxRapidity          () <<endl
	    << "RAP_N_BINS"    << nmbRapidityBins      () <<endl
	    << "CUT_PT"        << ptCutEnabled         () <<endl
	    << "PT_MIN"        << ptCutMin             () <<endl
	    << "PT_MAX"        << ptCutMax             () <<endl
	    << "CUT_ETA"       << etaCutEnabled        () <<endl
	    << "ETA_MIN"       << etaCutMin            () <<endl
	    << "ETA_MAX"       << etaCutMax            () <<endl
	    << "PROD_MODE"     << productionMode       () <<endl
	    << "N_EVENTS"      << nmbEvents            () <<endl
	    << "PROD_PID"      << prodParticleId       () <<endl
	    << "RND_SEED"      << randomSeed           () <<endl
	    << "BREAKUP_MODE"  << beamBreakupMode      () <<endl
	    << "INTERFERENCE"  << interferenceEnabled  () <<endl
	    << "IF_STRENGTH"   << interferenceStrength () <<endl
	    << "INT_PT_MAX"    << maxPtInterference    () <<endl
	    << "INT_PT_N_BINS" << nmbPtBinsInterference() <<endl
	    << "COHERENT"      << coherentProduction   () <<endl
	    << "INCO_FACTOR"   << incoherentFactor     () <<endl;

	return out;
}

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