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