1 /**************************************************************************
2 * Copyright(c) 1998-2003, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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. *
15 **************************************************************************/
19 // Event generator of two-photon processes
20 // in ultra-peripheral ion collisions.
21 // 5 two-photon process are implemented, see comments to SetProcess().
23 // The example of the generator initialization for the process
24 // gamma gamma -> X in CaCa collisions at 7000 TeV/nucleon is the following:
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);
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.
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.
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.
53 // Author: Yuri.Kharlov@cern.ch
56 #include <TParticle.h>
57 #include <TParticlePDG.h>
58 #include <TDatabasePDG.h>
59 #include <TClonesArray.h>
61 #include "AliPythia.h"
63 #include <AliGenTPHIC.h>
65 #include <TPHICcommon.h>
69 //------------------------------------------------------------
71 AliGenTPHIC::AliGenTPHIC() :
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
86 SetMC(new TPHICgen());
87 fPythia = AliPythia::Instance();
88 fParticles = new TClonesArray("TParticle",100);
103 //____________________________________________________________
104 AliGenTPHIC::~AliGenTPHIC()
106 // Destroys the object, deletes and disposes all TParticles currently on list.
108 fParticles->Delete();
114 //____________________________________________________________
115 void AliGenTPHIC::Init()
117 // Initialize the generator TPHIC
119 fTPHICgen->Initialize();
123 //____________________________________________________________
124 void AliGenTPHIC::Generate()
126 // Generate one event of two-photon process.
127 // Gaussian smearing on the vertex is done if selected.
128 // All particles from the TPHIC/PYTHIA event listing
129 // are stored in TreeK, and only final-state particles are tracked.
130 // The event differectial cross section is assigned as a weight
131 // to each track of the event.
133 Float_t polar[3] = {0,0,0};
134 Float_t origin0[3] = {0,0,0};
136 Float_t origin[3] = {0,0,0};
140 Int_t ks,kf,iparent,nt, trackIt;
142 const Float_t kconv=0.001/2.999792458e8;
144 fTPHICgen->GenerateEvent();
145 weight = GetXSectionCurrent();
146 if (gAlice->GetEvNumber()>=fDebugEventFirst &&
147 gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(1);
148 fPythia->ImportParticles(fParticles,"All");
151 Info("Generate()","one event is produced");
154 for (j=0;j<3;j++) origin0[j]=fOrigin[j];
156 if(fVertexSmear==kPerEvent) {
159 origin0[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
160 TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
163 time0 += fOsigma[2]/TMath::Ccgs()*
164 TMath::Cos(2*random[0]*TMath::Pi())*
165 TMath::Sqrt(-2*TMath::Log(random[1]));
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());
185 tof = time0 + kconv*iparticle->T();
186 if (ks == 1) trackIt = 1;
188 PushTrack(fTrackIt*trackIt,iparent,kf,p,origin,polar,tof,kPPrimary,nt,weight,ks);
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]);
196 fTPHICgen->SetNEVENT(fEvent);
197 if (fDebug == 1 && fEvent%100 == 0) {
198 Info("Generate","Event %d\n",fEvent);
203 //____________________________________________________________
204 void AliGenTPHIC::SetEventListRange(Int_t eventFirst, Int_t eventLast)
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;
213 //____________________________________________________________
214 void AliGenTPHIC::SetProcess (Int_t proc )
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)
224 fTPHICgen->SetIPROC(proc);
226 //____________________________________________________________
227 void AliGenTPHIC::SetBeamEnergy (Float_t energy)
229 // Set energy of the beam ion per nucleon in GeV
230 fTPHICgen->SetEBMN(energy);
232 //____________________________________________________________
233 void AliGenTPHIC::SetBeamZ (Int_t z )
235 // Set beam ion charge
238 //____________________________________________________________
239 void AliGenTPHIC::SetBeamA (Int_t a )
241 // Set beam ion atomic number
244 //____________________________________________________________
245 void AliGenTPHIC::SetYggRange (Float_t ymin, Float_t ymax)
247 // Set rapidity range of 2-photon system for the
248 // luminosity function calculation
249 fTPHICgen->SetYMIN(ymin);
250 fTPHICgen->SetYMAX(ymax);
252 //____________________________________________________________
253 void AliGenTPHIC::SetMggRange (Float_t mmin, Float_t mmax)
255 // Set invariant mass range of 2-photon system for the
256 // luminosity function calculation
257 fTPHICgen->SetAMIN(mmin);
258 fTPHICgen->SetAMAX(mmax);
260 //____________________________________________________________
261 void AliGenTPHIC::SetNgridY (Int_t ny )
263 // Set number of nodes on the grid along the rapidity axis
264 // to calculate the 2-photon luminosity function
265 fTPHICgen->SetNY(ny);
267 //____________________________________________________________
268 void AliGenTPHIC::SetNgridM (Int_t nm )
270 // Set number of nodes on the grid along the mass axis
271 // to calculate the 2-photon luminosity function
272 fTPHICgen->SetNMAS(nm);
274 //____________________________________________________________
275 void AliGenTPHIC::SetLumFunName (TString name )
277 // Set filename to store the 2-photon luminosity
278 // function calculated on the grid
279 fTPHICgen->SetLUMFIL(name);
281 //____________________________________________________________
282 void AliGenTPHIC::SetLumFunFlag (Int_t flag )
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
289 fTPHICgen->SetILUMF(flag);
291 //____________________________________________________________
292 void AliGenTPHIC::SetKfFermion (Int_t kf )
294 // Set a PDG flavour code of a fermion for the process 3,
295 // gamma gamma -> fermion+ fermion-
296 fTPHICgen->SetKFERM(kf);
298 //____________________________________________________________
299 void AliGenTPHIC::SetKfOnium (Int_t kf )
301 // Set a PDG flavour code of a quarkonium for the process 2,
302 // gamma gamma -> quarkonium
303 fTPHICgen->SetKFONIUM(kf);
305 //____________________________________________________________
306 void AliGenTPHIC::SetMassOnium (Float_t mass )
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);
314 //____________________________________________________________
315 void AliGenTPHIC::SetGGwidthOnium(Float_t width )
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);
323 //____________________________________________________________
324 void AliGenTPHIC::SetKfVmesons (Int_t kf1, Int_t kf2)
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);
334 //____________________________________________________________
335 Float_t AliGenTPHIC::GetGGmass ()
337 // Get invariant mass of generated 2-photon system
338 return fTPHICgen->GetWSQ();
340 //____________________________________________________________
341 Float_t AliGenTPHIC::GetGGrapidity()
343 // Get rapidity of generated 2-photon system
344 return fTPHICgen->GetYGG();
346 //____________________________________________________________
347 Float_t AliGenTPHIC::GetG1mass ()
349 // Get a mass of the first virtual photon of
350 // the 2-photon process (-sqrt(q1^2)).
351 return fTPHICgen->GetXMG1();
353 //____________________________________________________________
354 Float_t AliGenTPHIC::GetG2mass ()
356 // Get a mass of the second virtual photon of
357 // the 2-photon process (-sqrt(q2^2)).
358 return fTPHICgen->GetXMG2();
360 //____________________________________________________________
361 TClonesArray* AliGenTPHIC::GetParticleList ()
363 // Get array of particles of the event listing
366 //____________________________________________________________
367 TLorentzVector AliGenTPHIC::MomentumRecNucl1()
369 // Get 4-momentum of the first recoil nucleus after
370 // the 2-photon process.
371 return TLorentzVector(fTPHICgen->GetPTAG1(0),
372 fTPHICgen->GetPTAG1(1),
373 fTPHICgen->GetPTAG1(2),
374 fTPHICgen->GetPTAG1(3));
376 //____________________________________________________________
377 TLorentzVector AliGenTPHIC::MomentumRecNucl2()
379 // Get 4-momentum of the first recoil nucleus after
380 // the 2-photon process.
381 return TLorentzVector(fTPHICgen->GetPTAG2(0),
382 fTPHICgen->GetPTAG2(1),
383 fTPHICgen->GetPTAG2(2),
384 fTPHICgen->GetPTAG2(3));
386 //____________________________________________________________
387 Float_t AliGenTPHIC::GetXSectionCurrent()
389 // Get the cross section of the produced event for the
390 // Monte Carlo integral cross section calculation
391 return fTPHICgen->GetXSCUR();
393 //____________________________________________________________
394 Float_t AliGenTPHIC::GetXSection ()
396 // Get the integral cross section of the process
398 return fTPHICgen->GetXSTOT();
400 //____________________________________________________________
401 Float_t AliGenTPHIC::GetXSectionError ()
403 // Get the error of the integral cross section of the process
405 return fTPHICgen->GetXSTOTE();