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