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