estarlight/src/e_starlight.cpp

0001 ///////////////////////////////////////////////////////////////////////////
0002 //
0003 //    Copyright 2017
0004 //
0005 //    This file is part of estarlight.
0006 //
0007 //    starlight is free software: you can redistribute it and/or modify
0008 //    it under the terms of the GNU General Public License as published by
0009 //    the Free Software Foundation, either version 3 of the License, or
0010 //    (at your option) any later version.
0011 //
0012 //    starlight is distributed in the hope that it will be useful,
0013 //    but WITHOUT ANY WARRANTY; without even the implied warranty of
0014 //    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0015 //    GNU General Public License for more details.
0016 //
0017 //    You should have received a copy of the GNU General Public License
0018 //    along with starlight. If not, see <http://www.gnu.org/licenses/>.
0019 //
0020 ///////////////////////////////////////////////////////////////////////////
0021 //
0022 // File and Version Information:
0023 // $Rev:: 263                         $: revision of last commit
0024 // $Author:: mlomnitz                  $: author of last commit
0025 // $Date:: 02/28/2017 #$: date of last commit
0026 //
0027 // Description:
0028 //
0029 // Main interface for eSTARlight code inherited from STARlight
0030 //
0031 //
0032 ///////////////////////////////////////////////////////////////////////////
0033
0034
0035 #include <iostream>
0036 #include <fstream>
0037 #include <cstdlib>
0038
0039 #include "starlightconfig.h"
0040
0041 #ifdef ENABLE_PYTHIA
0042 #include "PythiaStarlight.h"
0043 #endif
0044
0045 #ifdef ENABLE_PYTHIA6
0046 #include "starlightpythia.h"
0047 #endif
0048
0049 #include "reportingUtils.h"
0050 #include "inputParameters.h"
0051 #include "eventchannel.h"
0052 #include "gammagammasingle.h"
0053 #include "gammaavm.h"
0054 #include "twophotonluminosity.h"
0055 #include "gammaaluminosity.h"
0056 #include "gammaeluminosity.h"
0057 #include "incoherentPhotonNucleusLuminosity.h"
0058 #include "e_starlight.h"
0059
0060
0061 using namespace std;
0062 using namespace starlightConstants;
0063
0064
0065 e_starlight::e_starlight() :
0066         _beamSystem            (0),
0067         _eventChannel          (0),
0068         _nmbEventsPerFile      (100),
0069         _isInitialised         (false),
0070         _inputParameters       (0)
0071 { }
0072
0073
0074 e_starlight::~e_starlight()
0075 { }
0076
0077
0078 bool
0079 e_starlight::init()
0080 {
0081   // ???
0082         if(Starlight_VERSION_MAJOR == 9999)
0083     {
0084       cout  << endl << "#########################################" << endl
0085             << " Initialising Starlight version: trunk..." << endl
0086             << "#########################################" << endl << endl;
0087     }
0088     else
0089     {
0090       cout  << endl << "#########################################" << endl
0091             << " Initialising Starlight version: " << Starlight_VERSION_MAJOR << "."
0092             << Starlight_VERSION_MINOR << "." << Starlight_VERSION_MINOR_MINOR << "..." << endl
0093             << "#########################################" << endl << endl;
0094     }
0095
0096     _nmbEventsPerFile    = _inputParameters->nmbEvents();  // for now we write only one file...
0097
0098     _beamSystem = new beamBeamSystem(*_inputParameters);
0099
0100     // This sets the precision used by cout for printing
0101     cout.setf(ios_base::fixed,ios_base::floatfield);
0102     cout.precision(3);
0103
0104         std::string _baseFileName;
0105         _baseFileName = _inputParameters->baseFileName();
0106        _lumLookUpTableFileName = _baseFileName + ".txt";
0107
0108     const bool lumTableIsValid = luminosityTableIsValid();
0109
0110     // Do some sanity checks of the input parameters here.
0111     if( _inputParameters->targetBeamA() == 0 ){
0112       //      printErr << endl << "One of the two beams must be an electron in eSTARlight" << endl;
0113       return false;
0114     }
0115
0116     bool createChannel = true;
0117     switch (_inputParameters->interactionType())    {
0118       // Photon-photon scattering is not implemented in eSTARlight as of yet
0119       /*
0120         case PHOTONPHOTON:
0121         if (!lumTableIsValid) {
0122             printInfo << "creating luminosity table for photon-photon channel" << endl;
0123             twoPhotonLuminosity(*_inputParameters, _beamSystem->beam1(), _beamSystem->beam2());
0124         }
0125         break;
0126       */
0127     case E_PHOTONPOMERONNARROW:  // narrow and wide resonances use
0128     case E_PHOTONPOMERONWIDE:    // the same luminosity function
0129         if (!lumTableIsValid) {
0130             printInfo << "creating luminosity table for coherent photon-Pomeron channel" <<endl;
0131             photonElectronLuminosity lum(*_inputParameters, *_beamSystem);
0132         }
0133         break;
0134         // Incoherent electro-production is still pending
0135         /*        case PHOTONPOMERONINCOHERENT:  // narrow and wide resonances use
0136                 if (!lumTableIsValid) {
0137                         printInfo << "creating luminosity table for incoherent photon-Pomeron channel" << endl;
0138                         incoherentPhotonNucleusLuminosity lum(*_inputParameters, *_beamSystem);
0139             }
0140             break;*/
0141     default:
0142         {
0143             printWarn << "unknown interaction type '" << _inputParameters->interactionType() << "'."
0144                       << " cannot initialize eSTARlight." << endl;
0145             return false;
0146         }
0147     }
0148
0149     if(createChannel)
0150     {
0151       if (!createEventChannel())
0152           return false;
0153     }
0154
0155     _isInitialised = true;
0156     return true;
0157 }
0158
0159
0160 eXEvent
0161 e_starlight::produceEvent()
0162 {
0163     if (!_isInitialised) {
0164       //        printErr << "trying to generate event but Starlight is not initialised. aborting." << endl;
0165         exit(-1);
0166     }
0167     return _eventChannel->e_produceEvent();
0168 }
0169
0170
0171 bool
0172 e_starlight::luminosityTableIsValid() const
0173 {
0174     printInfo << "using random seed = " << _inputParameters->randomSeed() << endl;
0175
0176     ifstream lumLookUpTableFile(_lumLookUpTableFileName.c_str());
0177     lumLookUpTableFile.precision(15);
0178     if ((!lumLookUpTableFile) || (!lumLookUpTableFile.good())) {
0179         return false;
0180     }
0181
0182     unsigned int targetBeamZ, targetBeamA;
0183     double       beamLorentzGamma = 0;
0184     double       maxW = 0, minW = 0;
0185     unsigned int nmbWBins;
0186     double       maxRapidity = 0;
0187     unsigned int nmbRapidityBins;
0188     int          productionMode, beamBreakupMode;
0189     //Remove interference options for eX collisions (no logner symmetric)
0190     //  bool         interferenceEnabled = false;
0191     //  double       interferenceStrength = 0;
0192     bool         coherentProduction = false;
0193     double       incoherentFactor = 0, deuteronSlopePar = 0, maxPtInterference = 0;
0194     int          nmbPtBinsInterference;
0195     std::string  validationKey;
0196     if (!(lumLookUpTableFile
0197           >> validationKey
0198           >> targetBeamZ >> targetBeamA
0199           >> beamLorentzGamma
0200           >> maxW >> minW >> nmbWBins
0201           >> maxRapidity >> nmbRapidityBins
0202           >> productionMode
0203           >> beamBreakupMode
0204           >> maxPtInterference
0205           >> nmbPtBinsInterference
0206           >> coherentProduction >> incoherentFactor
0207           >> deuteronSlopePar))
0208       // cannot read parameters from lookup table file
0209         return false;
0210
0211     std::string validationKeyEnd;
0212     while(!lumLookUpTableFile.eof())
0213     {
0214       lumLookUpTableFile >> validationKeyEnd;
0215     }
0216     lumLookUpTableFile.close();
0217     return (validationKey == _inputParameters->parameterValueKey() && validationKeyEnd == validationKey);
0218     return true;
0219 }
0220
0221
0222 bool
0223 e_starlight::createEventChannel()
0224 {
0225     switch (_inputParameters->prodParticleType()) {
0226     case ELECTRON:
0227     case MUON:
0228     case TAUON:
0229         case TAUONDECAY:
0230       printWarn << "cannot construct Gammagammaleptonpair event channel." << endl;
0231       return false;
0232
0233     case A2:        // jetset/pythia
0234     case ETA:       // jetset/pythia
0235     case ETAPRIME:  // jetset/pythia
0236     case ETAC:      // jetset/pythia
0237     case F0:        // jetset/pythia
0238         {
0239
0240           //#ifdef ENABLE_PYTHIA
0241             // PythiaOutput = true;
0242           //                cout<<"Pythia is enabled!"<<endl;
0243 //          return true;
0244 //#else
0245 //          printWarn << "Starlight is not compiled against Pythia8; "
0246 //                    << "jetset event channel cannot be used." << endl;
0247 //          return false;
0248 //#endif
0249         }
0250     case F2:
0251     case F2PRIME:
0252     case ZOVERZ03:
0253     case AXION: // AXION HACK
0254         {
0255           //  #ifdef ENABLE_PYTHIA
0256                 cout<<" This is f2, f2prim, rho^0 rho^0, or axion "<<endl;
0257             _eventChannel= new Gammagammasingle(*_inputParameters, *_beamSystem);
0258             if (_eventChannel)
0259                 return true;
0260             else {
0261                 printWarn << "cannot construct Gammagammasingle event channel." << endl;
0262                 return false;
0263             }
0264         }
0265     case RHO:
0266     case RHOZEUS:
0267     case FOURPRONG:
0268     case OMEGA:
0269     case OMEGA_pi0gamma:
0270     case PHI:
0271     case JPSI:
0272     case JPSI2S:
0273     case JPSI2S_ee:
0274     case JPSI2S_mumu:
0275     case JPSI_ee:
0276     case JPSI_mumu:
0277     case UPSILON:
0278     case UPSILON_ee:
0279     case UPSILON_mumu:
0280     case UPSILON2S:
0281     case UPSILON2S_ee:
0282     case UPSILON2S_mumu:
0283     case UPSILON3S:
0284     case UPSILON3S_ee:
0285     case UPSILON3S_mumu:
0286         {
0287             if (_inputParameters->interactionType() == E_PHOTONPOMERONNARROW) {
0288                 _eventChannel = new e_Gammaanarrowvm(*_inputParameters, *_beamSystem);
0289                 if (_eventChannel)
0290                     return true;
0291                 else {
0292                     printWarn << "cannot construct Gammaanarrowvm event channel." << endl;
0293                     return false;
0294                 }
0295             }
0296
0297             if (_inputParameters->interactionType() == E_PHOTONPOMERONWIDE) {
0298                 _eventChannel = new e_Gammaawidevm(*_inputParameters, *_beamSystem);
0299                 if (_eventChannel)
0300                     return true;
0301                 else {
0302                     printWarn << "cannot construct Gammaawidevm event channel." << endl;
0303                     return false;
0304                 }
0305             }
0306
0307
0308             printWarn << "interaction type '" << _inputParameters->interactionType() << "' "
0309                       << "cannot be used with particle type '" << _inputParameters->prodParticleType() << "'. "
0310                       << "cannot create event channel." << endl;
0311             return false;
0312         }
0313     default:
0314         {
0315             printWarn << "unknown event channel '" << _inputParameters->prodParticleType() << "'."
0316                       << " cannot create event channel." << endl;
0317             return false;
0318         }
0319     }
0320 }