STARLIGHT/starlight/src/starlight.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 //
  20 ///////////////////////////////////////////////////////////////////////////
  21 //
  22 // File and Version Information:
  23 // $Rev:: 164                         $: revision of last commit
  24 // $Author:: odjuvsla                 $: author of last commit
  25 // $Date:: 2013-10-06 16:18:08 +0200 #$: date of last commit
  26 //
  27 // Description:
  28 //
  29 //
  30 //
  31 ///////////////////////////////////////////////////////////////////////////
  32
  33
  34 #include <iostream>
  35 #include <fstream>
  36 #include <cstdlib>
  37
  38 #include "starlightconfig.h"
  39
  40 #ifdef ENABLE_PYTHIA
  41 #include "PythiaStarlight.h"
  42 #endif
  43
  44 #ifdef ENABLE_DPMJET
  45 #include "starlightdpmjet.h"
  46 #endif
  47
  48 #ifdef ENABLE_PYTHIA6
  49 #include "starlightpythia.h"
  50 #endif
  51
  52 #include "reportingUtils.h"
  53 #include "inputParameters.h"
  54 #include "eventchannel.h"
  55 #include "gammagammaleptonpair.h"
  56 #include "gammagammasingle.h"
  57 #include "gammaavm.h"
  58 #include "psifamily.h"
  59 #include "twophotonluminosity.h"
  60 #include "gammaaluminosity.h"
  61 #include "incoherentPhotonNucleusLuminosity.h"
  62 #include "upcevent.h"
  63 #include "eventfilewriter.h"
  64 #include "starlight.h"
  65
  66
  67 using namespace std;
  68 using namespace starlightConstants;
  69
  70
  71 starlight::starlight() :
  72                 _beam0                 (0),
  73                 _beam1                 (0),
  74                 _beamSystem            (0),
  75                 _eventChannel          (0),
  76                 _nmbEventsPerFile      (100),
  77                 _nmbEventsToGenerate   (10),
  78                 _configFileName        ("slight.in"),
  79                 _eventDataFileName     ("slight.out"),
  80                 _lumLookUpTableFileName("slight.txt"),
  81                 _isInitialised         (false)
  82 { }
  83
  84
  85 starlight::~starlight()
  86 { }
  87
  88
  89 bool
  90 starlight::init()
  91 {
  92         if(Starlight_VERSION_MAJOR == 9999)
  93         {
  94                 cout << "##################################" << endl
  95                 << " Initialising Starlight version: trunk..." << endl
  96                 << "##################################" << endl;
  97         }
  98         else
  99         {
 100                 cout << "##################################" << endl
 101                 << " Initialising Starlight version: " << Starlight_VERSION_MAJOR << "."
 102                 << Starlight_VERSION_MINOR << "." << Starlight_VERSION_MINOR_MINOR << "..." << endl
 103                 << "##################################" << endl;
 104         }
 105
 106         _nmbEventsPerFile    = inputParametersInstance.nmbEvents();  // for now we write only one file...
 107
 108         _beamSystem = new beamBeamSystem;
 109
 110 //      cout << "Created beam system with beam lorentz gamma: " << _beamSystem->beamLorentzGamma() << endl;
 111
 112         // streamsize precision(15);
 113         cout.setf(ios_base::fixed,ios_base::floatfield);
 114         cout.precision(15);
 115         const bool lumTableIsValid = luminosityTableIsValid();
 116         bool createChannel = true;
 117         switch (inputParametersInstance.interactionType())      {
 118         case PHOTONPHOTON:
 119                 if (!lumTableIsValid) {
 120                         printInfo << "creating luminosity table for photon-photon channel" << endl;
 121                         twoPhotonLuminosity(_beamSystem->beam1(), _beamSystem->beam2());
 122                 }
 123                 break;
 124         case PHOTONPOMERONNARROW:  // narrow and wide resonances use
 125         case PHOTONPOMERONWIDE:    // the same luminosity function
 126                 if (!lumTableIsValid) {
 127                         printInfo << "creating luminosity table for coherent photon-Pomeron channel" <<endl;
 128                         photonNucleusLuminosity lum(*_beamSystem);
 129                 }
 130                 break;
 131         case PHOTONPOMERONINCOHERENT:  // narrow and wide resonances use
 132                 if (!lumTableIsValid) {
 133                         printInfo << "creating luminosity table for incoherent photon-Pomeron channel" << endl;
 134                         incoherentPhotonNucleusLuminosity lum(*_beamSystem);
 135                 }
 136                 break;
 137 #ifdef ENABLE_DPMJET
 138         case PHOTONUCLEARSINGLE:
 139                 createChannel = false;
 140                 _eventChannel = new starlightDpmJet(*_beamSystem);
 141                 std::cout << "CREATING PHOTONUCLEAR/DPMJET SINGLE" << std::endl;
 142                 dynamic_cast<starlightDpmJet*>(_eventChannel)->setSingleMode();
 143                 dynamic_cast<starlightDpmJet*>(_eventChannel)->setMinGammaEnergy(inputParametersInstance.minGammaEnergy());
 144                 dynamic_cast<starlightDpmJet*>(_eventChannel)->setMaxGammaEnergy(inputParametersInstance.maxGammaEnergy());
 145                 dynamic_cast<starlightDpmJet*>(_eventChannel)->init();
 146                 break;
 147         case PHOTONUCLEARDOUBLE:
 148                 createChannel = false;
 149                 _eventChannel = new starlightDpmJet(*_beamSystem);
 150                 std::cout << "CREATING PHOTONUCLEAR/DPMJET DOUBLE" << std::endl;
 151                 dynamic_cast<starlightDpmJet*>(_eventChannel)->setDoubleMode();
 152                 dynamic_cast<starlightDpmJet*>(_eventChannel)->setMinGammaEnergy(inputParametersInstance.minGammaEnergy());
 153                 dynamic_cast<starlightDpmJet*>(_eventChannel)->setMaxGammaEnergy(inputParametersInstance.maxGammaEnergy());
 154                 dynamic_cast<starlightDpmJet*>(_eventChannel)->init();
 155                 break;
 156         case PHOTONUCLEARSINGLEPA:
 157                 createChannel = false;
 158                 _eventChannel = new starlightDpmJet(*_beamSystem);
 159                 std::cout << "CREATING PHOTONUCLEAR/DPMJET SINGLE" << std::endl;
 160                 dynamic_cast<starlightDpmJet*>(_eventChannel)->setSingleMode();
 161                 dynamic_cast<starlightDpmJet*>(_eventChannel)->setProtonMode();
 162                 dynamic_cast<starlightDpmJet*>(_eventChannel)->setMinGammaEnergy(inputParametersInstance.minGammaEnergy());
 163                 dynamic_cast<starlightDpmJet*>(_eventChannel)->setMaxGammaEnergy(inputParametersInstance.maxGammaEnergy());
 164                 dynamic_cast<starlightDpmJet*>(_eventChannel)->init();
 165                 break;
 166 #endif
 167 #ifdef ENABLE_PYTHIA6
 168         case PHOTONUCLEARSINGLEPAPY:
 169                 createChannel = false;
 170                 _eventChannel = new starlightPythia(*_beamSystem);
 171                 std::cout << "CREATING PHOTONUCLEAR/PYTHIA SINGLE" << std::endl;
 172                 dynamic_cast<starlightPythia*>(_eventChannel)->setSingleMode();
 173                 dynamic_cast<starlightPythia*>(_eventChannel)->setMinGammaEnergy(inputParametersInstance.minGammaEnergy());
 174                 dynamic_cast<starlightPythia*>(_eventChannel)->setMaxGammaEnergy(inputParametersInstance.maxGammaEnergy());
 175                 dynamic_cast<starlightPythia*>(_eventChannel)->init(inputParametersInstance.pythiaParams(), inputParametersInstance.pythiaFullEventRecord());
 176                 break;
 177 #endif
 178         default:
 179                 {
 180                         printWarn << "unknown interaction type '" << inputParametersInstance.interactionType() << "'."
 181                                   << " cannot initialize starlight." << endl;
 182                         return false;
 183                 }
 184         }
 185
 186         if(createChannel)
 187         {
 188           if (!createEventChannel())
 189                   return false;
 190         }
 191
 192         _isInitialised = true;
 193         return true;
 194 }
 195
 196
 197 upcEvent
 198 starlight::produceEvent()
 199 {
 200         if (!_isInitialised) {
 201                 printErr << "trying to generate event but Starlight is not initialised. aborting." << endl;
 202                 exit(-1);
 203         }
 204         return _eventChannel->produceEvent();
 205 }
 206
 207
 208 bool
 209 starlight::luminosityTableIsValid() const
 210 {
 211         printInfo << "using random seed = " << inputParametersInstance.randomSeed() << endl;
 212
 213         ifstream lumLookUpTableFile(_lumLookUpTableFileName.c_str());
 214         lumLookUpTableFile.precision(15);
 215         if ((!lumLookUpTableFile) || (!lumLookUpTableFile.good())) {
 216                 // printWarn << "cannot open file '" << _lumLookUpTableFileName << "'" << endl;
 217                 return false;
 218         }
 219
 220         unsigned int beam1Z, beam1A, beam2Z, beam2A;
 221         double       beamLorentzGamma = 0;
 222         double       maxW = 0, minW = 0;
 223         unsigned int nmbWBins;
 224         double       maxRapidity = 0;
 225         unsigned int nmbRapidityBins;
 226         int          productionMode, beamBreakupMode;
 227         bool         interferenceEnabled = false;
 228         double       interferenceStrength = 0;
 229         bool         coherentProduction = false;
 230         double       incoherentFactor = 0, deuteronSlopePar = 0, maxPtInterference = 0;
 231         int          nmbPtBinsInterference;
 232         std::string  validationKey;
 233         if (!(lumLookUpTableFile
 234               >> validationKey
 235               >> beam1Z >> beam1A
 236               >> beam2Z >> beam2A
 237               >> beamLorentzGamma
 238               >> maxW >> minW >> nmbWBins
 239               >> maxRapidity >> nmbRapidityBins
 240               >> productionMode
 241               >> beamBreakupMode
 242               >> interferenceEnabled >> interferenceStrength
 243               >> coherentProduction >> incoherentFactor
 244               >> deuteronSlopePar
 245               >> maxPtInterference
 246               >> nmbPtBinsInterference))
 247                 // cannot read parameters from lookup table file
 248                 return false;
 249
 250         std::string validationKeyEnd;
 251         while(!lumLookUpTableFile.eof())
 252         {
 253           lumLookUpTableFile >> validationKeyEnd;
 254         }
 255         lumLookUpTableFile.close();
 256         return (validationKey == inputParametersInstance.parameterValueKey() && validationKeyEnd == validationKey);
 257         return true;
 258 }
 259
 260
 261 bool
 262 starlight::createEventChannel()
 263 {
 264         switch (inputParametersInstance.prodParticleType()) {
 265         case ELECTRON:
 266         case MUON:
 267         case TAUON:
 268                 {
 269                         _eventChannel = new Gammagammaleptonpair(*_beamSystem);
 270                         if (_eventChannel)
 271                                 return true;
 272                         else {
 273                                 printWarn << "cannot construct Gammagammaleptonpair event channel." << endl;
 274                                 return false;
 275                         }
 276                 }
 277         case A2:        // jetset
 278         case ETA:       // jetset
 279         case ETAPRIME:  // jetset
 280         case ETAC:      // jetset
 281         case F0:        // jetset
 282                 {
 283 #ifdef ENABLE_PYTHIA
 284                         // PythiaOutput = true;
 285                         cout<<"Pythia is enabled!"<<endl;
 286 //                      return true;
 287 #else
 288                         printWarn << "Starlight is not compiled against Pythia8; "
 289                                   << "jetset event channel cannot be used." << endl;
 290                         return false;
 291 #endif
 292                 }
 293         case F2:
 294         case F2PRIME:
 295         case ZOVERZ03:
 296                 {
 297                   //  #ifdef ENABLE_PYTHIA
 298                         cout<<" This is f2, f2prim, zoverz03 "<<endl;
 299                         _eventChannel= new Gammagammasingle(*_beamSystem);
 300                         if (_eventChannel)
 301                                 return true;
 302                         else {
 303                                 printWarn << "cannot construct Gammagammasingle event channel." << endl;
 304                                 return false;
 305                         }
 306                         // #endif
 307                         //                      printWarn << "Starlight is not compiled against Pythia8; "
 308                         //          << "Gammagammasingle event channel cannot be used." << endl;
 309                         // return false;
 310                 }
 311         case RHO:
 312         case RHOZEUS:
 313         case FOURPRONG:
 314         case OMEGA:
 315         case PHI:
 316         case JPSI:
 317         case JPSI2S:
 318         case JPSI2S_ee:
 319         case JPSI2S_mumu:
 320         case JPSI_ee:
 321         case JPSI_mumu:
 322         case UPSILON:
 323         case UPSILON_ee:
 324         case UPSILON_mumu:
 325         case UPSILON2S:
 326         case UPSILON2S_ee:
 327         case UPSILON2S_mumu:
 328         case UPSILON3S:
 329         case UPSILON3S_ee:
 330         case UPSILON3S_mumu:
 331                 {
 332                         if (inputParametersInstance.interactionType() == PHOTONPOMERONNARROW) {
 333                                 _eventChannel = new Gammaanarrowvm(*_beamSystem);
 334                                 if (_eventChannel)
 335                                         return true;
 336                                 else {
 337                                         printWarn << "cannot construct Gammaanarrowvm event channel." << endl;
 338                                         return false;
 339                                 }
 340                         }
 341
 342                         if (inputParametersInstance.interactionType() == PHOTONPOMERONWIDE) {
 343                                 _eventChannel = new Gammaawidevm(*_beamSystem);
 344                                 if (_eventChannel)
 345                                         return true;
 346                                 else {
 347                                         printWarn << "cannot construct Gammaawidevm event channel." << endl;
 348                                         return false;
 349                                 }
 350                         }
 351
 352                         if (inputParametersInstance.interactionType() == PHOTONPOMERONINCOHERENT) {
 353                                 _eventChannel = new Gammaaincoherentvm(*_beamSystem);
 354                                 if (_eventChannel)
 355                                         return true;
 356                                 else {
 357                                         printWarn << "cannot construct Gammaanarrowvm event channel." << endl;
 358                                         return false;
 359                                 }
 360                         }
 361
 362                         printWarn << "interaction type '" << inputParametersInstance.interactionType() << "' "
 363                                   << "cannot be used with particle type '" << inputParametersInstance.prodParticleType() << "'. "
 364                                   << "cannot create event channel." << endl;
 365                         return false;
 366                 }
 367         default:
 368                 {
 369                         printWarn << "unknown event channel '" << inputParametersInstance.prodParticleType() << "'."
 370                                   << " cannot create event channel." << endl;
 371                         return false;
 372                 }
 373         }
 374 }