Updated from the TPC-PreRelease branch
[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$
eb342d2d 18Revision 1.1 2000/06/09 20:47:27 morsch
19AliGenerator interface class to HIJING using THijing (test version)
20
110410f9 21*/
22
23#include "AliGenHijing.h"
24#include "AliRun.h"
25
26#include <TArrayI.h>
27#include <TParticle.h>
28#include <THijing.h>
29
30 ClassImp(AliGenHijing)
31
32AliGenHijing::AliGenHijing()
33 :AliGenerator()
34{
35// Constructor
36}
37
38AliGenHijing::AliGenHijing(Int_t npart)
39 :AliGenerator(npart)
40{
41// Default PbPb collisions at 5. 5 TeV
42//
43 SetEnergyCMS();
44 SetImpactParameterRange();
45 SetTarget();
46 SetProjectile();
47 fKeep=0;
48 fQuench=1;
49 fShadowing=1;
50 fTrigger=0;
51 fDecaysOff=1;
52 fEvaluate=0;
53 fSelectAll=0;
eb342d2d 54 fFlavor=0;
110410f9 55}
56
57AliGenHijing::AliGenHijing(const AliGenHijing & Hijing)
58{
59// copy constructor
60}
61
62
63AliGenHijing::~AliGenHijing()
64{
65// Destructor
66}
67
68void AliGenHijing::Init()
69{
70// Initialisation
71 SetMC(new THijing(fEnergyCMS, fFrame, fProjectile, fTarget,
72 fAProjectile, fZProjectile, fATarget, fZTarget,
73 fMinImpactParam, fMaxImpactParam));
74
75 fHijing=(THijing*) fgMCEvGen;
76 fHijing->Initialize();
77 fHijing->SetIHPR2(3, fTrigger);
78 fHijing->SetIHPR2(4, fQuench);
79 fHijing->SetIHPR2(6, fShadowing);
80 fHijing->SetIHPR2(12, fDecaysOff);
81 fHijing->SetIHPR2(21, fKeep);
82//
83 if (fEvaluate) EvaluateCrossSections();
84}
85
86void AliGenHijing::Generate()
87{
88// Generate one event
89
90 Float_t polar[3] = {0,0,0};
91 Float_t origin[3]= {0,0,0};
92 Float_t origin0[3]= {0,0,0};
93 Float_t p[3], random[6];
94 Float_t tof;
95
96 static TClonesArray *particles;
97// converts from mm/c to s
98 const Float_t kconv=0.001/2.999792458e8;
99//
100 Int_t nt=0;
101 Int_t jev=0;
102 Int_t j, kf, ks, imo;
eb342d2d 103 kf=0;
104
110410f9 105 if(!particles) particles=new TClonesArray("TParticle",10000);
106
107 fTrials=0;
108 for (j=0;j<3;j++) origin0[j]=fOrigin[j];
109 if(fVertexSmear==perEvent) {
110 gMC->Rndm(random,6);
111 for (j=0;j<3;j++) {
112 origin0[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
113 TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
114// fHijing->SetMSTP(151,0);
115 }
116 } else if (fVertexSmear==perTrack) {
117// fHijing->SetMSTP(151,0);
118 for (j=0;j<3;j++) {
119// fHijing->SetPARP(151+j, fOsigma[j]*10.);
120 }
121 }
122 while(1)
123 {
124
125 fHijing->GenerateEvent();
126 fTrials++;
127 fHijing->ImportParticles(particles,"All");
128 Int_t np = particles->GetEntriesFast();
129 printf("\n **************************************************%d\n",np);
130 Int_t nc=0;
131 if (np == 0 ) continue;
132 Int_t i;
133 Int_t * newPos = new Int_t[np];
eb342d2d 134 for (i = 0; i<np; i++) *(newPos+i)=i;
110410f9 135
eb342d2d 136 for (i = 0; i<np; i++) {
110410f9 137 TParticle * iparticle = (TParticle *) particles->At(i);
138
139 Bool_t hasMother = (iparticle->GetFirstMother() >=0);
140 Bool_t hasDaughter = (iparticle->GetFirstDaughter() >=0);
141 Bool_t selected = kTRUE;
142 Bool_t hasSelectedDaughters = kFALSE;
143
110410f9 144
eb342d2d 145 kf = iparticle->GetPdgCode();
146 if (!fSelectAll) selected = KinematicSelection(iparticle)&&SelectFlavor(kf);
110410f9 147 if (hasDaughter && !selected) hasSelectedDaughters = DaughtersSelection(iparticle, particles);
148//
149// Put particle on the stack if it is either selected or it is the mother of at least one seleted particle
150//
151
152 if (selected || hasSelectedDaughters) {
153 nc++;
154 ks = iparticle->GetStatusCode();
155 p[0]=iparticle->Px();
156 p[1]=iparticle->Py();
157 p[2]=iparticle->Pz();
158 origin[0]=origin0[0]+iparticle->Vx()/10;
159 origin[1]=origin0[1]+iparticle->Vy()/10;
160 origin[2]=origin0[2]+iparticle->Vz()/10;
161 tof=kconv*iparticle->T();
162 imo=-1;
163 if (hasMother) {
164 imo=iparticle->GetFirstMother();
165 imo=*(newPos+imo);
166 }
167
168// printf("\n selected iparent %d %d %d \n",i, kf, imo);
169 if (hasDaughter) {
170 gAlice->SetTrack(0,imo,kf,p,origin,polar,
171 tof,"Primary",nt);
172 } else {
173 gAlice->SetTrack(fTrackIt,imo,kf,p,origin,polar,
174 tof,"Secondary",nt);
175 }
176 *(newPos+i)=nt;
177 } // selected
178 } // particle loop
179 delete newPos;
180 printf("\n I've put %i particles on the stack \n",nc);
181 if (nc > 0) {
182 jev+=nc;
183 if (jev >= fNpart || fNpart == -1) {
184 fKineBias=Float_t(fNpart)/Float_t(fTrials);
185 printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
186 break;
187 }
188 }
189 } // event loop
190}
191
192Bool_t AliGenHijing::KinematicSelection(TParticle *particle)
193{
194// Perform kinematic selection
195 Float_t px=particle->Px();
196 Float_t py=particle->Py();
197 Float_t pz=particle->Pz();
198 Float_t e=particle->Energy();
199
200//
201// transverse momentum cut
202 Float_t pt=TMath::Sqrt(px*px+py*py);
203 if (pt > fPtMax || pt < fPtMin)
204 {
205// printf("\n failed pt cut %f %f %f \n",pt,fPtMin,fPtMax);
206 return kFALSE;
207 }
208//
209// momentum cut
210 Float_t p=TMath::Sqrt(px*px+py*py+pz*pz);
211 if (p > fPMax || p < fPMin)
212 {
213// printf("\n failed p cut %f %f %f \n",p,fPMin,fPMax);
214 return kFALSE;
215 }
216
217//
218// theta cut
219 Float_t theta = Float_t(TMath::ATan2(Double_t(pt),Double_t(pz)));
220 if (theta > fThetaMax || theta < fThetaMin)
221 {
222
223// printf("\n failed theta cut %f %f %f \n",theta,fThetaMin,fThetaMax);
224 return kFALSE;
225 }
226
227//
228// rapidity cut
229 Float_t y = 0.5*TMath::Log((e+pz)/(e-pz));
230 if (y > fYMax || y < fYMin)
231 {
232// printf("\n failed y cut %f %f %f \n",y,fYMin,fYMax);
233 return kFALSE;
234 }
235
236//
237// phi cut
238 Float_t phi=Float_t(TMath::ATan2(Double_t(py),Double_t(px)));
239 if (phi > fPhiMax || phi < fPhiMin)
240 {
241// printf("\n failed phi cut %f %f %f \n",phi,fPhiMin,fPhiMax);
242 return kFALSE;
243 }
244
245 return kTRUE;
246}
247
248void AliGenHijing::KeepFullEvent()
249{
250 fKeep=1;
251}
252
253void AliGenHijing::EvaluateCrossSections()
254{
255// Glauber Calculation of geometrical x-section
256//
257 Float_t xTot=0.; // barn
258 Float_t xTotHard=0.; // barn
259 Float_t xPart=0.; // barn
260 Float_t xPartHard=0.; // barn
261 Float_t sigmaHard=0.1; // mbarn
262 Float_t bMin=0.;
263 Float_t bMax=fHijing->GetHIPR1(34)+fHijing->GetHIPR1(35);
264 const Float_t kdib=0.2;
265 Int_t kMax=Int_t((bMax-bMin)/kdib)+1;
266
267
268 printf("\n Projectile Radius (fm): %f \n",fHijing->GetHIPR1(34));
269 printf("\n Target Radius (fm): %f \n",fHijing->GetHIPR1(35));
270 Int_t i;
271 Float_t oldvalue=0.;
272
273 for (i=0; i<kMax; i++)
274 {
275 Float_t xb=bMin+i*kdib;
276 Float_t ov;
277 ov=fHijing->Profile(xb);
278 Float_t gb = 2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*(1.-TMath::Exp(-fHijing->GetHINT1(12)*ov));
279 Float_t gbh = 2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*sigmaHard*ov;
280 xTot+=gb;
281 xTotHard+=gbh;
282 if (xb > fMinImpactParam && xb < fMaxImpactParam)
283 {
284 xPart+=gb;
285 xPartHard+=gbh;
286 }
287
288 if ((xTot-oldvalue)/oldvalue<0.0001) break;
289 oldvalue=xTot;
290 printf("\n Total cross section (barn): %d %f %f \n",i, xb, xTot);
291 printf("\n Hard cross section (barn): %d %f %f \n\n",i, xb, xTotHard);
292 }
293 printf("\n Total cross section (barn): %f \n",xTot);
294 printf("\n Hard cross section (barn): %f \n \n",xTotHard);
295 printf("\n Partial cross section (barn): %f %f \n",xPart, xPart/xTot*100.);
296 printf("\n Partial hard cross section (barn): %f %f \n",xPartHard, xPartHard/xTotHard*100.);
297}
298
299Bool_t AliGenHijing::DaughtersSelection(TParticle* iparticle, TClonesArray* particles)
300{
301//
302// Looks recursively if one of the daughters has been selected
303//
304// printf("\n Consider daughters %d:",iparticle->GetPdgCode());
305 Int_t imin=-1;
306 Int_t imax=-1;
307 Int_t i;
308 Bool_t hasDaughters= (iparticle->GetFirstDaughter() >=0);
309 Bool_t selected=kFALSE;
310 if (hasDaughters) {
311 imin=iparticle->GetFirstDaughter();
312 imax=iparticle->GetLastDaughter();
110410f9 313 for (i=imin; i<= imax; i++){
314 TParticle * jparticle = (TParticle *) particles->At(i);
eb342d2d 315 Int_t ip=jparticle->GetPdgCode();
316 if (KinematicSelection(jparticle)&&SelectFlavor(ip)) {selected=kTRUE; break;}
110410f9 317 if (DaughtersSelection(jparticle, particles)) {selected=kTRUE; break; }
318 }
319 } else {
320 return kFALSE;
321 }
322
323 return selected;
324}
325
326
eb342d2d 327Bool_t AliGenHijing::SelectFlavor(Int_t pid)
328{
329// Select flavor of particle
330// 0: all
331// 4: charm and beauty
332// 5: beauty
333 if (fFlavor == 0) return kTRUE;
334
335 Int_t ifl=TMath::Abs(pid/100);
336 if (ifl > 10) ifl/=10;
337 return ((fFlavor==4 && (ifl==4 || ifl==5)) ||
338 (fFlavor==5 && ifl==5));
339
340}
110410f9 341
342AliGenHijing& AliGenHijing::operator=(const AliGenHijing& rhs)
343{
344// Assignment operator
345 return *this;
346}