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