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