]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STARLIGHT/starlight/src/starlight.cpp
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
STARLIGHT code and interface
[u/mrichter/AliRoot.git] / STARLIGHT / starlight / src / starlight.cpp
CommitLineData
da32329d
AM		 1///////////////////////////////////////////////////////////////////////////
2//
3// Copyright 2010
4//
5// This file is part of starlight.
6//
7// starlight is free software: you can redistribute it and/or modify
8// it under the terms of the GNU General Public License as published by
9// the Free Software Foundation, either version 3 of the License, or
10// (at your option) any later version.
11//
12// starlight is distributed in the hope that it will be useful,
13// but WITHOUT ANY WARRANTY; without even the implied warranty of
14// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15// GNU General Public License for more details.
16//
17// You should have received a copy of the GNU General Public License
18// along with starlight. If not, see <http://www.gnu.org/licenses/>.
19//
20///////////////////////////////////////////////////////////////////////////
21//
22// File and Version Information:
23// $Rev:: 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
67using namespace std;
68using namespace starlightConstants;
69
70
71starlight::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
85starlight::~starlight()
86{ }
87
88
89bool
90starlight::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
197upcEvent
198starlight::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
208bool
209starlight::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
261bool
262starlight::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}