Generation of dictionaries and rootmaps with Root6
[u/mrichter/AliRoot.git] / TPHIC / AliGenTPHIC.cxx
CommitLineData
0b5dd071 1/**************************************************************************
2 * Copyright(c) 1998-2003, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 * *
15 **************************************************************************/
16
17/* $Id$ */
18
0b5dd071 19// Event generator of two-photon processes
20// in ultra-peripheral ion collisions.
21// 5 two-photon process are implemented, see comments to SetProcess().
22//%
23// The example of the generator initialization for the process
24// gamma gamma -> X in CaCa collisions at 7000 TeV/nucleon is the following:
25//%
26// AliGenTPHIC *gener = new AliGenTPHIC();
27// gener->SetProcess(1);
28// gener->SetBeamEnergy(3500.);
29// gener->SetBeamZ(20);
30// gener->SetBeamA(40);
31// gener->SetMggRange(70.,200.);
32// gener->SetYggRange(-5.,5.);
33// gener->SetLumFunName("lum_ca_70_200.dat");
34// gener->SetLumFunFlag(-1);
35// gener->Init();
36//%
37// The two-photon luminosity function needed for the process cross section
38// calculation is time-consuming, therefore it can be calculated once for a
39// selected two-photon energy and rapidity range and selected ion species,
40// saved into a file and then can be reused in further calculation.
41//%
42// The integral cross section of the process is calculated in each event
43// by the Monte Carlo integration procedure, the differential cross section
44// of each thrown event is assigned as a weight to each track of the event
45// which can be then used during the event analysis by retreiving this weight
46// from any of the event tracks.
47//
48// The manual of the fortran generator is available in
49// $ALICE_ROOT/TPHIC/TPHIC_doc.ps.gz
50// For the two-photon physics in heavy ion collisions see
51// G.Baur et al, Phys.Rep. 364 (2002), 359.
52//%
53// Author: Yuri.Kharlov@cern.ch
54// 15 April 2003
55
56#include <TParticle.h>
57#include <TParticlePDG.h>
58#include <TDatabasePDG.h>
37b09b91 59#include <TClonesArray.h>
60
0b5dd071 61#include "AliPythia.h"
62#include "AliRun.h"
63#include <AliGenTPHIC.h>
64#include <TPHICgen.h>
65#include <TPHICcommon.h>
66
67ClassImp(AliGenTPHIC)
68
69//------------------------------------------------------------
70
f893d4a5 71AliGenTPHIC::AliGenTPHIC() :
72 AliGenMC(),
73 fTPHICgen(0x0),
74 fPythia(0x0),
75 fParticles(0x0),
76 fEvent(-1),
77 fDebug(0),
78 fDebugEventFirst(-1),
79 fDebugEventLast(-1)
0b5dd071 80{
81 // Constructor: create generator instance,
82 // create particle array
83 // Set TPHIC parameters to default values:
84 // eta_b production in Ca-Ca collisions at 7 A*TeV
85
86 SetMC(new TPHICgen());
27d01a0d 87 fTPHICgen = (TPHICgen*) fMCEvGen;
0b5dd071 88 fPythia = AliPythia::Instance();
89 fParticles = new TClonesArray("TParticle",100);
90
91 SetProcess ();
92 SetBeamEnergy();
93 SetBeamZ ();
94 SetBeamA ();
95 SetYggRange ();
96 SetMggRange ();
97 SetNgridY ();
98 SetNgridM ();
99 SetLumFunName();
100 SetLumFunFlag();
101 SetKfOnium ();
102}
103
104//____________________________________________________________
0b5dd071 105AliGenTPHIC::~AliGenTPHIC()
106{
107 // Destroys the object, deletes and disposes all TParticles currently on list.
108 if (fParticles) {
109 fParticles->Delete();
110 delete fParticles;
111 fParticles = 0;
112 }
113}
114
115//____________________________________________________________
116void AliGenTPHIC::Init()
117{
118 // Initialize the generator TPHIC
119
120 fTPHICgen->Initialize();
121 fEvent = 0;
122}
123
124//____________________________________________________________
125void AliGenTPHIC::Generate()
126{
127 // Generate one event of two-photon process.
128 // Gaussian smearing on the vertex is done if selected.
129 // All particles from the TPHIC/PYTHIA event listing
130 // are stored in TreeK, and only final-state particles are tracked.
131 // The event differectial cross section is assigned as a weight
132 // to each track of the event.
133
8feee6f7 134 Float_t polar[3] = {0,0,0};
135 Float_t origin0[3] = {0,0,0};
21391258 136 Float_t time0 = 0.;
8feee6f7 137 Float_t origin[3] = {0,0,0};
0b5dd071 138 Float_t p[3], tof;
139 Double_t weight;
140
141 Int_t ks,kf,iparent,nt, trackIt;
142 Float_t random[6];
143 const Float_t kconv=0.001/2.999792458e8;
144
145 fTPHICgen->GenerateEvent();
146 weight = GetXSectionCurrent();
147 if (gAlice->GetEvNumber()>=fDebugEventFirst &&
148 gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(1);
149 fPythia->ImportParticles(fParticles,"All");
150
151 if (fDebug == 1)
152 Info("Generate()","one event is produced");
153
154 Int_t j;
21391258 155 for (j=0;j<3;j++) origin0[j]=fOrigin[j];
156 time0 = fTimeOrigin;
0b5dd071 157 if(fVertexSmear==kPerEvent) {
158 Rndm(random,6);
159 for (j=0;j<3;j++) {
160 origin0[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
161 TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
162 }
21391258 163 Rndm(random,2);
164 time0 += fOsigma[2]/TMath::Ccgs()*
165 TMath::Cos(2*random[0]*TMath::Pi())*
166 TMath::Sqrt(-2*TMath::Log(random[1]));
0b5dd071 167 }
168
169 Int_t ip;
170 Int_t np = fParticles->GetEntriesFast();
171 TParticle *iparticle;
172// Int_t* pParent = new Int_t[np];
173 for (ip=0; ip<np; ip++) {
174 iparticle = (TParticle *) fParticles->At(ip);
175 ks = iparticle->GetStatusCode();
176// // No initial state partons
177// if (ks==21) continue;
178 p[0] = iparticle->Px();
179 p[1] = iparticle->Py();
180 p[2] = iparticle->Pz();
181 origin[0] = origin0[0]+iparticle->Vx()/10.;
182 origin[1] = origin0[1]+iparticle->Vy()/10.;
183 origin[2] = origin0[2]+iparticle->Vz()/10.;
184 kf = CheckPDGCode(iparticle->GetPdgCode());
185 iparent = -1;
21391258 186 tof = time0 + kconv*iparticle->T();
0b5dd071 187 if (ks == 1) trackIt = 1;
188 else trackIt = 0;
642f15cf 189 PushTrack(fTrackIt*trackIt,iparent,kf,p,origin,polar,tof,kPPrimary,nt,weight,ks);
0b5dd071 190 KeepTrack(nt);
191
192 if (fDebug == 2)
193 printf("id=%+4d, parent=%3d, ks=%d, p = (%+11.4e,%+11.4e,%+11.4e) GeV\n",
194 kf,iparent,fTrackIt*trackIt,p[0],p[1],p[2]);
195 }
196 fEvent++;
197 fTPHICgen->SetNEVENT(fEvent);
198 if (fDebug == 1 && fEvent%100 == 0) {
199 Info("Generate","Event %d\n",fEvent);
200 fTPHICgen->Finish();
201 }
202}
203
204//____________________________________________________________
205void AliGenTPHIC::SetEventListRange(Int_t eventFirst, Int_t eventLast)
206{
207 // Set a range of event numbers, for which a table
208 // of generated particle will be printed
209 fDebugEventFirst = eventFirst;
210 fDebugEventLast = eventLast;
211 if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
212}
213
214//____________________________________________________________
215void AliGenTPHIC::SetProcess (Int_t proc )
216{
217 // Set process number:
218 // proc=1 - gamma gamma -> X
219 // proc=2 - gamma gamma -> quarkonium
220 // proc=3 - gamma gamma -> fermion+ fermion-
221 // proc=4 - gamma gamma -> W+ W-
222 // proc=5 - not implemented
223 // proc=6 - gamma gamma -> V1 V2 (vector meson pair)
224
225 fTPHICgen->SetIPROC(proc);
226}
227//____________________________________________________________
228 void AliGenTPHIC::SetBeamEnergy (Float_t energy)
229{
230 // Set energy of the beam ion per nucleon in GeV
231 fTPHICgen->SetEBMN(energy);
232}
233//____________________________________________________________
234 void AliGenTPHIC::SetBeamZ (Int_t z )
235{
236 // Set beam ion charge
237 fTPHICgen->SetIZ(z);
238}
239//____________________________________________________________
240 void AliGenTPHIC::SetBeamA (Int_t a )
241{
242 // Set beam ion atomic number
243 fTPHICgen->SetIA(a);
244}
245//____________________________________________________________
246 void AliGenTPHIC::SetYggRange (Float_t ymin, Float_t ymax)
247{
248 // Set rapidity range of 2-photon system for the
249 // luminosity function calculation
250 fTPHICgen->SetYMIN(ymin);
251 fTPHICgen->SetYMAX(ymax);
252}
253//____________________________________________________________
254 void AliGenTPHIC::SetMggRange (Float_t mmin, Float_t mmax)
255{
256 // Set invariant mass range of 2-photon system for the
257 // luminosity function calculation
258 fTPHICgen->SetAMIN(mmin);
259 fTPHICgen->SetAMAX(mmax);
260}
261//____________________________________________________________
262 void AliGenTPHIC::SetNgridY (Int_t ny )
263{
264 // Set number of nodes on the grid along the rapidity axis
265 // to calculate the 2-photon luminosity function
266 fTPHICgen->SetNY(ny);
267}
268//____________________________________________________________
269 void AliGenTPHIC::SetNgridM (Int_t nm )
270{
271 // Set number of nodes on the grid along the mass axis
272 // to calculate the 2-photon luminosity function
273 fTPHICgen->SetNMAS(nm);
274}
275//____________________________________________________________
276 void AliGenTPHIC::SetLumFunName (TString name )
277{
278 // Set filename to store the 2-photon luminosity
279 // function calculated on the grid
280 fTPHICgen->SetLUMFIL(name);
281}
282//____________________________________________________________
283 void AliGenTPHIC::SetLumFunFlag (Int_t flag )
284{
285 // Set flag to calculate the 2-photon luminosity function:
286 // flag=-1 if a new lumimosity function to be calculated
287 // and stored in a file
288 // flag=+1 if a previously calculated function to be read
289 // from a file
290 fTPHICgen->SetILUMF(flag);
291}
292//____________________________________________________________
293 void AliGenTPHIC::SetKfFermion (Int_t kf )
294{
295 // Set a PDG flavour code of a fermion for the process 3,
296 // gamma gamma -> fermion+ fermion-
297 fTPHICgen->SetKFERM(kf);
298}
299//____________________________________________________________
300 void AliGenTPHIC::SetKfOnium (Int_t kf )
301{
302 // Set a PDG flavour code of a quarkonium for the process 2,
303 // gamma gamma -> quarkonium
304 fTPHICgen->SetKFONIUM(kf);
305}
306//____________________________________________________________
307 void AliGenTPHIC::SetMassOnium (Float_t mass )
308{
309 // Set a quarkonium mass [GeV] for the process 2 if it
310 // differes from the Pythia's particle table.
311 // For the well-known quarkonia no need to set this mass
312 // because it will be taken from the Pythia table
313 fTPHICgen->SetXMRES(mass);
314}
315//____________________________________________________________
316 void AliGenTPHIC::SetGGwidthOnium(Float_t width )
317{
318 // Set 2-photon partial width [GeV] of the quarkonium for
319 // the process 2 if it differes fromthe Pythia's particle table.
320 // For the well-known quarkonia no need to set this width
321 // because it will be taken from the Pythia table
322 fTPHICgen->SetXGGRES(width);
323}
324//____________________________________________________________
325 void AliGenTPHIC::SetKfVmesons (Int_t kf1, Int_t kf2)
326{
327 // Set PDG flavour codes of the two vector vesons
328 // for the process 2: gamma gamma -> V1 V2
329 // So far this processes is implemented for the following
330 // mesons and their combinations only:
331 // pho0 (113), omega (223), phi (333), J/psi (443)
332 fTPHICgen->SetKV1(kf1);
333 fTPHICgen->SetKV2(kf2);
334}
335//____________________________________________________________
336 Float_t AliGenTPHIC::GetGGmass ()
337{
338 // Get invariant mass of generated 2-photon system
339 return fTPHICgen->GetWSQ();
340}
341//____________________________________________________________
342 Float_t AliGenTPHIC::GetGGrapidity()
343{
344 // Get rapidity of generated 2-photon system
345 return fTPHICgen->GetYGG();
346}
347//____________________________________________________________
348 Float_t AliGenTPHIC::GetG1mass ()
349{
350 // Get a mass of the first virtual photon of
351 // the 2-photon process (-sqrt(q1^2)).
352 return fTPHICgen->GetXMG1();
353}
354//____________________________________________________________
355 Float_t AliGenTPHIC::GetG2mass ()
356{
357 // Get a mass of the second virtual photon of
358 // the 2-photon process (-sqrt(q2^2)).
359 return fTPHICgen->GetXMG2();
360}
361//____________________________________________________________
362 TClonesArray* AliGenTPHIC::GetParticleList ()
363{
364 // Get array of particles of the event listing
365 return fParticles;
366}
367//____________________________________________________________
368 TLorentzVector AliGenTPHIC::MomentumRecNucl1()
369{
370 // Get 4-momentum of the first recoil nucleus after
371 // the 2-photon process.
2c23f883 372 return TLorentzVector(fTPHICgen->GetPTAG1(0),
373 fTPHICgen->GetPTAG1(1),
0b5dd071 374 fTPHICgen->GetPTAG1(2),
2c23f883 375 fTPHICgen->GetPTAG1(3));
0b5dd071 376}
377//____________________________________________________________
378 TLorentzVector AliGenTPHIC::MomentumRecNucl2()
379{
380 // Get 4-momentum of the first recoil nucleus after
381 // the 2-photon process.
2c23f883 382 return TLorentzVector(fTPHICgen->GetPTAG2(0),
383 fTPHICgen->GetPTAG2(1),
0b5dd071 384 fTPHICgen->GetPTAG2(2),
2c23f883 385 fTPHICgen->GetPTAG2(3));
0b5dd071 386}
387//____________________________________________________________
388 Float_t AliGenTPHIC::GetXSectionCurrent()
389{
390 // Get the cross section of the produced event for the
391 // Monte Carlo integral cross section calculation
392 return fTPHICgen->GetXSCUR();
393}
394//____________________________________________________________
395 Float_t AliGenTPHIC::GetXSection ()
396{
397 // Get the integral cross section of the process
398 // calculated so far
399 return fTPHICgen->GetXSTOT();
400}
401//____________________________________________________________
402 Float_t AliGenTPHIC::GetXSectionError ()
403{
404 // Get the error of the integral cross section of the process
405 // calculated so far
406 return fTPHICgen->GetXSTOTE();
407}