MUON trigger classes to collaborate with CTP (E. Lopez Torres)
[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>
59#include "AliPythia.h"
60#include "AliRun.h"
61#include <AliGenTPHIC.h>
62#include <TPHICgen.h>
63#include <TPHICcommon.h>
64
65ClassImp(AliGenTPHIC)
66
67//------------------------------------------------------------
68
69AliGenTPHIC::AliGenTPHIC()
70{
71 // Constructor: create generator instance,
72 // create particle array
73 // Set TPHIC parameters to default values:
74 // eta_b production in Ca-Ca collisions at 7 A*TeV
75
76 SetMC(new TPHICgen());
77 fPythia = AliPythia::Instance();
78 fParticles = new TClonesArray("TParticle",100);
79
80 SetProcess ();
81 SetBeamEnergy();
82 SetBeamZ ();
83 SetBeamA ();
84 SetYggRange ();
85 SetMggRange ();
86 SetNgridY ();
87 SetNgridM ();
88 SetLumFunName();
89 SetLumFunFlag();
90 SetKfOnium ();
91}
92
93//____________________________________________________________
a06bb94f 94AliGenTPHIC::AliGenTPHIC(const AliGenTPHIC & gen):
95 AliGenMC(gen)
0b5dd071 96{
97 // copy constructor
98 gen.Copy(*this);
99}
100
101//____________________________________________________________
102AliGenTPHIC::~AliGenTPHIC()
103{
104 // Destroys the object, deletes and disposes all TParticles currently on list.
105 if (fParticles) {
106 fParticles->Delete();
107 delete fParticles;
108 fParticles = 0;
109 }
110}
111
112//____________________________________________________________
113void AliGenTPHIC::Init()
114{
115 // Initialize the generator TPHIC
116
117 fTPHICgen->Initialize();
118 fEvent = 0;
119}
120
121//____________________________________________________________
122void AliGenTPHIC::Generate()
123{
124 // Generate one event of two-photon process.
125 // Gaussian smearing on the vertex is done if selected.
126 // All particles from the TPHIC/PYTHIA event listing
127 // are stored in TreeK, and only final-state particles are tracked.
128 // The event differectial cross section is assigned as a weight
129 // to each track of the event.
130
131 Float_t polar[3]= {0,0,0};
132 Float_t origin0[3],origin[3];
133 Float_t p[3], tof;
134 Double_t weight;
135
136 Int_t ks,kf,iparent,nt, trackIt;
137 Float_t random[6];
138 const Float_t kconv=0.001/2.999792458e8;
139
140 fTPHICgen->GenerateEvent();
141 weight = GetXSectionCurrent();
142 if (gAlice->GetEvNumber()>=fDebugEventFirst &&
143 gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(1);
144 fPythia->ImportParticles(fParticles,"All");
145
146 if (fDebug == 1)
147 Info("Generate()","one event is produced");
148
149 Int_t j;
150 for (j=0;j<3;j++) origin[j]=fOrigin[j];
151 if(fVertexSmear==kPerEvent) {
152 Rndm(random,6);
153 for (j=0;j<3;j++) {
154 origin0[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
155 TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
156 }
157 }
158
159 Int_t ip;
160 Int_t np = fParticles->GetEntriesFast();
161 TParticle *iparticle;
162// Int_t* pParent = new Int_t[np];
163 for (ip=0; ip<np; ip++) {
164 iparticle = (TParticle *) fParticles->At(ip);
165 ks = iparticle->GetStatusCode();
166// // No initial state partons
167// if (ks==21) continue;
168 p[0] = iparticle->Px();
169 p[1] = iparticle->Py();
170 p[2] = iparticle->Pz();
171 origin[0] = origin0[0]+iparticle->Vx()/10.;
172 origin[1] = origin0[1]+iparticle->Vy()/10.;
173 origin[2] = origin0[2]+iparticle->Vz()/10.;
174 kf = CheckPDGCode(iparticle->GetPdgCode());
175 iparent = -1;
176 tof = kconv*iparticle->T();
177 if (ks == 1) trackIt = 1;
178 else trackIt = 0;
642f15cf 179 PushTrack(fTrackIt*trackIt,iparent,kf,p,origin,polar,tof,kPPrimary,nt,weight,ks);
0b5dd071 180 KeepTrack(nt);
181
182 if (fDebug == 2)
183 printf("id=%+4d, parent=%3d, ks=%d, p = (%+11.4e,%+11.4e,%+11.4e) GeV\n",
184 kf,iparent,fTrackIt*trackIt,p[0],p[1],p[2]);
185 }
186 fEvent++;
187 fTPHICgen->SetNEVENT(fEvent);
188 if (fDebug == 1 && fEvent%100 == 0) {
189 Info("Generate","Event %d\n",fEvent);
190 fTPHICgen->Finish();
191 }
192}
193
194//____________________________________________________________
195void AliGenTPHIC::SetEventListRange(Int_t eventFirst, Int_t eventLast)
196{
197 // Set a range of event numbers, for which a table
198 // of generated particle will be printed
199 fDebugEventFirst = eventFirst;
200 fDebugEventLast = eventLast;
201 if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
202}
203
204//____________________________________________________________
205void AliGenTPHIC::SetProcess (Int_t proc )
206{
207 // Set process number:
208 // proc=1 - gamma gamma -> X
209 // proc=2 - gamma gamma -> quarkonium
210 // proc=3 - gamma gamma -> fermion+ fermion-
211 // proc=4 - gamma gamma -> W+ W-
212 // proc=5 - not implemented
213 // proc=6 - gamma gamma -> V1 V2 (vector meson pair)
214
215 fTPHICgen->SetIPROC(proc);
216}
217//____________________________________________________________
218 void AliGenTPHIC::SetBeamEnergy (Float_t energy)
219{
220 // Set energy of the beam ion per nucleon in GeV
221 fTPHICgen->SetEBMN(energy);
222}
223//____________________________________________________________
224 void AliGenTPHIC::SetBeamZ (Int_t z )
225{
226 // Set beam ion charge
227 fTPHICgen->SetIZ(z);
228}
229//____________________________________________________________
230 void AliGenTPHIC::SetBeamA (Int_t a )
231{
232 // Set beam ion atomic number
233 fTPHICgen->SetIA(a);
234}
235//____________________________________________________________
236 void AliGenTPHIC::SetYggRange (Float_t ymin, Float_t ymax)
237{
238 // Set rapidity range of 2-photon system for the
239 // luminosity function calculation
240 fTPHICgen->SetYMIN(ymin);
241 fTPHICgen->SetYMAX(ymax);
242}
243//____________________________________________________________
244 void AliGenTPHIC::SetMggRange (Float_t mmin, Float_t mmax)
245{
246 // Set invariant mass range of 2-photon system for the
247 // luminosity function calculation
248 fTPHICgen->SetAMIN(mmin);
249 fTPHICgen->SetAMAX(mmax);
250}
251//____________________________________________________________
252 void AliGenTPHIC::SetNgridY (Int_t ny )
253{
254 // Set number of nodes on the grid along the rapidity axis
255 // to calculate the 2-photon luminosity function
256 fTPHICgen->SetNY(ny);
257}
258//____________________________________________________________
259 void AliGenTPHIC::SetNgridM (Int_t nm )
260{
261 // Set number of nodes on the grid along the mass axis
262 // to calculate the 2-photon luminosity function
263 fTPHICgen->SetNMAS(nm);
264}
265//____________________________________________________________
266 void AliGenTPHIC::SetLumFunName (TString name )
267{
268 // Set filename to store the 2-photon luminosity
269 // function calculated on the grid
270 fTPHICgen->SetLUMFIL(name);
271}
272//____________________________________________________________
273 void AliGenTPHIC::SetLumFunFlag (Int_t flag )
274{
275 // Set flag to calculate the 2-photon luminosity function:
276 // flag=-1 if a new lumimosity function to be calculated
277 // and stored in a file
278 // flag=+1 if a previously calculated function to be read
279 // from a file
280 fTPHICgen->SetILUMF(flag);
281}
282//____________________________________________________________
283 void AliGenTPHIC::SetKfFermion (Int_t kf )
284{
285 // Set a PDG flavour code of a fermion for the process 3,
286 // gamma gamma -> fermion+ fermion-
287 fTPHICgen->SetKFERM(kf);
288}
289//____________________________________________________________
290 void AliGenTPHIC::SetKfOnium (Int_t kf )
291{
292 // Set a PDG flavour code of a quarkonium for the process 2,
293 // gamma gamma -> quarkonium
294 fTPHICgen->SetKFONIUM(kf);
295}
296//____________________________________________________________
297 void AliGenTPHIC::SetMassOnium (Float_t mass )
298{
299 // Set a quarkonium mass [GeV] for the process 2 if it
300 // differes from the Pythia's particle table.
301 // For the well-known quarkonia no need to set this mass
302 // because it will be taken from the Pythia table
303 fTPHICgen->SetXMRES(mass);
304}
305//____________________________________________________________
306 void AliGenTPHIC::SetGGwidthOnium(Float_t width )
307{
308 // Set 2-photon partial width [GeV] of the quarkonium for
309 // the process 2 if it differes fromthe Pythia's particle table.
310 // For the well-known quarkonia no need to set this width
311 // because it will be taken from the Pythia table
312 fTPHICgen->SetXGGRES(width);
313}
314//____________________________________________________________
315 void AliGenTPHIC::SetKfVmesons (Int_t kf1, Int_t kf2)
316{
317 // Set PDG flavour codes of the two vector vesons
318 // for the process 2: gamma gamma -> V1 V2
319 // So far this processes is implemented for the following
320 // mesons and their combinations only:
321 // pho0 (113), omega (223), phi (333), J/psi (443)
322 fTPHICgen->SetKV1(kf1);
323 fTPHICgen->SetKV2(kf2);
324}
325//____________________________________________________________
326 Float_t AliGenTPHIC::GetGGmass ()
327{
328 // Get invariant mass of generated 2-photon system
329 return fTPHICgen->GetWSQ();
330}
331//____________________________________________________________
332 Float_t AliGenTPHIC::GetGGrapidity()
333{
334 // Get rapidity of generated 2-photon system
335 return fTPHICgen->GetYGG();
336}
337//____________________________________________________________
338 Float_t AliGenTPHIC::GetG1mass ()
339{
340 // Get a mass of the first virtual photon of
341 // the 2-photon process (-sqrt(q1^2)).
342 return fTPHICgen->GetXMG1();
343}
344//____________________________________________________________
345 Float_t AliGenTPHIC::GetG2mass ()
346{
347 // Get a mass of the second virtual photon of
348 // the 2-photon process (-sqrt(q2^2)).
349 return fTPHICgen->GetXMG2();
350}
351//____________________________________________________________
352 TClonesArray* AliGenTPHIC::GetParticleList ()
353{
354 // Get array of particles of the event listing
355 return fParticles;
356}
357//____________________________________________________________
358 TLorentzVector AliGenTPHIC::MomentumRecNucl1()
359{
360 // Get 4-momentum of the first recoil nucleus after
361 // the 2-photon process.
362 return TLorentzVector(fTPHICgen->GetPTAG1(1),
363 fTPHICgen->GetPTAG1(2),
364 fTPHICgen->GetPTAG1(3),
365 fTPHICgen->GetPTAG1(4));
366}
367//____________________________________________________________
368 TLorentzVector AliGenTPHIC::MomentumRecNucl2()
369{
370 // Get 4-momentum of the first recoil nucleus after
371 // the 2-photon process.
372 return TLorentzVector(fTPHICgen->GetPTAG2(1),
373 fTPHICgen->GetPTAG2(2),
374 fTPHICgen->GetPTAG2(3),
375 fTPHICgen->GetPTAG2(4));
376}
377//____________________________________________________________
378 Float_t AliGenTPHIC::GetXSectionCurrent()
379{
380 // Get the cross section of the produced event for the
381 // Monte Carlo integral cross section calculation
382 return fTPHICgen->GetXSCUR();
383}
384//____________________________________________________________
385 Float_t AliGenTPHIC::GetXSection ()
386{
387 // Get the integral cross section of the process
388 // calculated so far
389 return fTPHICgen->GetXSTOT();
390}
391//____________________________________________________________
392 Float_t AliGenTPHIC::GetXSectionError ()
393{
394 // Get the error of the integral cross section of the process
395 // calculated so far
396 return fTPHICgen->GetXSTOTE();
397}
a06bb94f 398
82764a9c 399void AliGenTPHIC::Copy(TObject&) const
a06bb94f 400{
401 //
402 // Copy
403 //
404 Fatal("Copy","Not implemented!\n");
405}
406