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