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