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