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