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