932737d13a6598874abc12fa7a4f5b5a4047f0b4
[u/mrichter/AliRoot.git] / TPHIC / AliGenTPHIC.cxx
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
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
65 ClassImp(AliGenTPHIC)
66
67 //------------------------------------------------------------
68
69 AliGenTPHIC::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 //____________________________________________________________
94 AliGenTPHIC::AliGenTPHIC(const AliGenTPHIC & gen):
95     AliGenMC(gen)
96 {
97   // copy constructor
98   gen.Copy(*this);
99 }
100
101 //____________________________________________________________
102 AliGenTPHIC::~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 //____________________________________________________________
113 void AliGenTPHIC::Init()
114 {
115   // Initialize the generator TPHIC
116
117   fTPHICgen->Initialize();
118   fEvent = 0;
119 }
120
121 //____________________________________________________________
122 void 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;
179     PushTrack(fTrackIt*trackIt,iparent,kf,p,origin,polar,tof,kPPrimary,nt,weight,ks);
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 //____________________________________________________________
195 void 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 //____________________________________________________________
205 void 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 }
398
399 void AliGenTPHIC::Copy(AliGenTPHIC&) const
400 {
401     //
402     // Copy 
403     //
404     Fatal("Copy","Not implemented!\n");
405 }
406