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