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