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