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