Protect EvaluateCrossSections() against division by zero.
[u/mrichter/AliRoot.git] / EVGEN / AliGenHijing.cxx
CommitLineData
110410f9 1/**************************************************************************
2 * Copyright(c) 1998-1999, 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$Log$
73bc6a56 18Revision 1.7 2000/10/02 21:28:06 fca
19Removal of useless dependecies via forward declarations
20
94de3818 21Revision 1.6 2000/09/11 13:23:37 morsch
22Write last seed to file (fortran lun 50) and reed back from same lun using calls to
23luget_hijing and luset_hijing.
24
b78ff279 25Revision 1.5 2000/09/07 16:55:40 morsch
26fHijing->Initialize(); after change of parameters. (Dmitri Yurevitch Peressounko)
27
ba244bac 28Revision 1.4 2000/07/11 18:24:56 fca
29Coding convention corrections + few minor bug fixes
30
aee8290b 31Revision 1.3 2000/06/30 12:08:36 morsch
32In member data: char* replaced by TString, Init takes care of resizing the strings to
338 characters required by Hijing.
34
3b8bd63f 35Revision 1.2 2000/06/15 14:15:05 morsch
36Add possibility for heavy flavor selection: charm and beauty.
37
eb342d2d 38Revision 1.1 2000/06/09 20:47:27 morsch
39AliGenerator interface class to HIJING using THijing (test version)
40
110410f9 41*/
42
43#include "AliGenHijing.h"
3b8bd63f 44#include "AliGenHijingEventHeader.h"
110410f9 45#include "AliRun.h"
94de3818 46#include "AliMC.h"
110410f9 47
48#include <TArrayI.h>
49#include <TParticle.h>
50#include <THijing.h>
51
52 ClassImp(AliGenHijing)
53
54AliGenHijing::AliGenHijing()
55 :AliGenerator()
56{
57// Constructor
58}
59
60AliGenHijing::AliGenHijing(Int_t npart)
61 :AliGenerator(npart)
62{
63// Default PbPb collisions at 5. 5 TeV
64//
65 SetEnergyCMS();
66 SetImpactParameterRange();
67 SetTarget();
68 SetProjectile();
69 fKeep=0;
70 fQuench=1;
71 fShadowing=1;
72 fTrigger=0;
73 fDecaysOff=1;
74 fEvaluate=0;
75 fSelectAll=0;
eb342d2d 76 fFlavor=0;
110410f9 77}
78
79AliGenHijing::AliGenHijing(const AliGenHijing & Hijing)
80{
81// copy constructor
82}
83
84
85AliGenHijing::~AliGenHijing()
86{
87// Destructor
88}
89
90void AliGenHijing::Init()
91{
92// Initialisation
3b8bd63f 93 fFrame.Resize(8);
94 fTarget.Resize(8);
95 fProjectile.Resize(8);
96
110410f9 97 SetMC(new THijing(fEnergyCMS, fFrame, fProjectile, fTarget,
98 fAProjectile, fZProjectile, fATarget, fZTarget,
99 fMinImpactParam, fMaxImpactParam));
100
101 fHijing=(THijing*) fgMCEvGen;
ba244bac 102
110410f9 103 fHijing->SetIHPR2(3, fTrigger);
104 fHijing->SetIHPR2(4, fQuench);
105 fHijing->SetIHPR2(6, fShadowing);
106 fHijing->SetIHPR2(12, fDecaysOff);
107 fHijing->SetIHPR2(21, fKeep);
b78ff279 108 fHijing->Rluset(50,0);
ba244bac 109 fHijing->Initialize();
b78ff279 110
111
110410f9 112//
113 if (fEvaluate) EvaluateCrossSections();
114}
115
116void AliGenHijing::Generate()
117{
118// Generate one event
119
120 Float_t polar[3] = {0,0,0};
121 Float_t origin[3]= {0,0,0};
122 Float_t origin0[3]= {0,0,0};
123 Float_t p[3], random[6];
124 Float_t tof;
125
126 static TClonesArray *particles;
127// converts from mm/c to s
128 const Float_t kconv=0.001/2.999792458e8;
129//
130 Int_t nt=0;
131 Int_t jev=0;
132 Int_t j, kf, ks, imo;
eb342d2d 133 kf=0;
134
110410f9 135 if(!particles) particles=new TClonesArray("TParticle",10000);
136
137 fTrials=0;
138 for (j=0;j<3;j++) origin0[j]=fOrigin[j];
aee8290b 139 if(fVertexSmear==kPerEvent) {
110410f9 140 gMC->Rndm(random,6);
141 for (j=0;j<3;j++) {
142 origin0[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
143 TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
144// fHijing->SetMSTP(151,0);
145 }
aee8290b 146 } else if (fVertexSmear==kPerTrack) {
110410f9 147// fHijing->SetMSTP(151,0);
148 for (j=0;j<3;j++) {
149// fHijing->SetPARP(151+j, fOsigma[j]*10.);
150 }
151 }
152 while(1)
153 {
154
155 fHijing->GenerateEvent();
156 fTrials++;
157 fHijing->ImportParticles(particles,"All");
158 Int_t np = particles->GetEntriesFast();
159 printf("\n **************************************************%d\n",np);
160 Int_t nc=0;
161 if (np == 0 ) continue;
162 Int_t i;
163 Int_t * newPos = new Int_t[np];
eb342d2d 164 for (i = 0; i<np; i++) *(newPos+i)=i;
110410f9 165
eb342d2d 166 for (i = 0; i<np; i++) {
110410f9 167 TParticle * iparticle = (TParticle *) particles->At(i);
168
169 Bool_t hasMother = (iparticle->GetFirstMother() >=0);
170 Bool_t hasDaughter = (iparticle->GetFirstDaughter() >=0);
171 Bool_t selected = kTRUE;
172 Bool_t hasSelectedDaughters = kFALSE;
173
110410f9 174
eb342d2d 175 kf = iparticle->GetPdgCode();
176 if (!fSelectAll) selected = KinematicSelection(iparticle)&&SelectFlavor(kf);
110410f9 177 if (hasDaughter && !selected) hasSelectedDaughters = DaughtersSelection(iparticle, particles);
178//
179// Put particle on the stack if it is either selected or it is the mother of at least one seleted particle
180//
181
182 if (selected || hasSelectedDaughters) {
183 nc++;
184 ks = iparticle->GetStatusCode();
185 p[0]=iparticle->Px();
186 p[1]=iparticle->Py();
187 p[2]=iparticle->Pz();
188 origin[0]=origin0[0]+iparticle->Vx()/10;
189 origin[1]=origin0[1]+iparticle->Vy()/10;
190 origin[2]=origin0[2]+iparticle->Vz()/10;
191 tof=kconv*iparticle->T();
192 imo=-1;
193 if (hasMother) {
194 imo=iparticle->GetFirstMother();
195 imo=*(newPos+imo);
196 }
197
198// printf("\n selected iparent %d %d %d \n",i, kf, imo);
199 if (hasDaughter) {
200 gAlice->SetTrack(0,imo,kf,p,origin,polar,
201 tof,"Primary",nt);
202 } else {
203 gAlice->SetTrack(fTrackIt,imo,kf,p,origin,polar,
204 tof,"Secondary",nt);
205 }
206 *(newPos+i)=nt;
207 } // selected
208 } // particle loop
209 delete newPos;
210 printf("\n I've put %i particles on the stack \n",nc);
211 if (nc > 0) {
212 jev+=nc;
213 if (jev >= fNpart || fNpart == -1) {
214 fKineBias=Float_t(fNpart)/Float_t(fTrials);
215 printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
216 break;
217 }
218 }
219 } // event loop
b78ff279 220 fHijing->Rluget(50,-1);
110410f9 221}
222
223Bool_t AliGenHijing::KinematicSelection(TParticle *particle)
224{
225// Perform kinematic selection
226 Float_t px=particle->Px();
227 Float_t py=particle->Py();
228 Float_t pz=particle->Pz();
229 Float_t e=particle->Energy();
230
231//
232// transverse momentum cut
233 Float_t pt=TMath::Sqrt(px*px+py*py);
234 if (pt > fPtMax || pt < fPtMin)
235 {
236// printf("\n failed pt cut %f %f %f \n",pt,fPtMin,fPtMax);
237 return kFALSE;
238 }
239//
240// momentum cut
241 Float_t p=TMath::Sqrt(px*px+py*py+pz*pz);
242 if (p > fPMax || p < fPMin)
243 {
244// printf("\n failed p cut %f %f %f \n",p,fPMin,fPMax);
245 return kFALSE;
246 }
247
248//
249// theta cut
250 Float_t theta = Float_t(TMath::ATan2(Double_t(pt),Double_t(pz)));
251 if (theta > fThetaMax || theta < fThetaMin)
252 {
253
254// printf("\n failed theta cut %f %f %f \n",theta,fThetaMin,fThetaMax);
255 return kFALSE;
256 }
257
258//
259// rapidity cut
260 Float_t y = 0.5*TMath::Log((e+pz)/(e-pz));
261 if (y > fYMax || y < fYMin)
262 {
263// printf("\n failed y cut %f %f %f \n",y,fYMin,fYMax);
264 return kFALSE;
265 }
266
267//
268// phi cut
269 Float_t phi=Float_t(TMath::ATan2(Double_t(py),Double_t(px)));
270 if (phi > fPhiMax || phi < fPhiMin)
271 {
272// printf("\n failed phi cut %f %f %f \n",phi,fPhiMin,fPhiMax);
273 return kFALSE;
274 }
275
276 return kTRUE;
277}
278
279void AliGenHijing::KeepFullEvent()
280{
281 fKeep=1;
282}
283
284void AliGenHijing::EvaluateCrossSections()
285{
286// Glauber Calculation of geometrical x-section
287//
288 Float_t xTot=0.; // barn
289 Float_t xTotHard=0.; // barn
290 Float_t xPart=0.; // barn
291 Float_t xPartHard=0.; // barn
292 Float_t sigmaHard=0.1; // mbarn
293 Float_t bMin=0.;
294 Float_t bMax=fHijing->GetHIPR1(34)+fHijing->GetHIPR1(35);
295 const Float_t kdib=0.2;
296 Int_t kMax=Int_t((bMax-bMin)/kdib)+1;
297
298
299 printf("\n Projectile Radius (fm): %f \n",fHijing->GetHIPR1(34));
300 printf("\n Target Radius (fm): %f \n",fHijing->GetHIPR1(35));
301 Int_t i;
302 Float_t oldvalue=0.;
303
304 for (i=0; i<kMax; i++)
305 {
306 Float_t xb=bMin+i*kdib;
307 Float_t ov;
308 ov=fHijing->Profile(xb);
309 Float_t gb = 2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*(1.-TMath::Exp(-fHijing->GetHINT1(12)*ov));
310 Float_t gbh = 2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*sigmaHard*ov;
311 xTot+=gb;
312 xTotHard+=gbh;
313 if (xb > fMinImpactParam && xb < fMaxImpactParam)
314 {
315 xPart+=gb;
316 xPartHard+=gbh;
317 }
73bc6a56 318 if (xTot == 0.) continue;
110410f9 319
73bc6a56 320 if ((xTot-oldvalue)/xTot<0.0001) break;
110410f9 321 oldvalue=xTot;
322 printf("\n Total cross section (barn): %d %f %f \n",i, xb, xTot);
323 printf("\n Hard cross section (barn): %d %f %f \n\n",i, xb, xTotHard);
324 }
325 printf("\n Total cross section (barn): %f \n",xTot);
326 printf("\n Hard cross section (barn): %f \n \n",xTotHard);
327 printf("\n Partial cross section (barn): %f %f \n",xPart, xPart/xTot*100.);
328 printf("\n Partial hard cross section (barn): %f %f \n",xPartHard, xPartHard/xTotHard*100.);
329}
330
331Bool_t AliGenHijing::DaughtersSelection(TParticle* iparticle, TClonesArray* particles)
332{
333//
334// Looks recursively if one of the daughters has been selected
335//
336// printf("\n Consider daughters %d:",iparticle->GetPdgCode());
337 Int_t imin=-1;
338 Int_t imax=-1;
339 Int_t i;
340 Bool_t hasDaughters= (iparticle->GetFirstDaughter() >=0);
341 Bool_t selected=kFALSE;
342 if (hasDaughters) {
343 imin=iparticle->GetFirstDaughter();
344 imax=iparticle->GetLastDaughter();
110410f9 345 for (i=imin; i<= imax; i++){
346 TParticle * jparticle = (TParticle *) particles->At(i);
eb342d2d 347 Int_t ip=jparticle->GetPdgCode();
348 if (KinematicSelection(jparticle)&&SelectFlavor(ip)) {selected=kTRUE; break;}
110410f9 349 if (DaughtersSelection(jparticle, particles)) {selected=kTRUE; break; }
350 }
351 } else {
352 return kFALSE;
353 }
354
355 return selected;
356}
357
358
eb342d2d 359Bool_t AliGenHijing::SelectFlavor(Int_t pid)
360{
361// Select flavor of particle
362// 0: all
363// 4: charm and beauty
364// 5: beauty
365 if (fFlavor == 0) return kTRUE;
366
367 Int_t ifl=TMath::Abs(pid/100);
368 if (ifl > 10) ifl/=10;
369 return ((fFlavor==4 && (ifl==4 || ifl==5)) ||
370 (fFlavor==5 && ifl==5));
371
372}
110410f9 373
3b8bd63f 374void AliGenHijing::MakeHeader()
375{
376// Builds the event header, to be called after each event
377 AliGenHijingEventHeader* header = new AliGenHijingEventHeader("Hijing");
378// header->SetDate(date);
379// header->SetRunNumber(run);
380// header->SetEventNumber(event);
381 header->SetNProduced(fHijing->GetNATT());
382 header->SetImpactParameter(fHijing->GetHINT1(19));
383 header->SetTotalEnergy(fHijing->GetEATT());
384 header->SetHardScatters(fHijing->GetJATT());
385 header->SetParticipants(fHijing->GetNP(), fHijing->GetNT());
386 header->SetCollisions(fHijing->GetN0(),
387 fHijing->GetN01(),
388 fHijing->GetN10(),
389 fHijing->GetN11());
390}
391
110410f9 392AliGenHijing& AliGenHijing::operator=(const AliGenHijing& rhs)
393{
394// Assignment operator
395 return *this;
396}