Coding rule violations corrected.
[u/mrichter/AliRoot.git] / THijing / 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 /* $Id$ */
17
18 // Generator using HIJING as an external generator
19 // The main HIJING options are accessable for the user through this interface.
20 // Uses the THijing implementation of TGenerator.
21 // Author:
22 // Andreas Morsch    (andreas.morsch@cern.ch)
23 //
24
25 #include <TGraph.h>
26 #include <THijing.h>
27 #include <TLorentzVector.h>
28 #include <TPDGCode.h>
29 #include <TParticle.h>
30
31 #include "AliGenHijing.h"
32 #include "AliGenHijingEventHeader.h"
33 #include "AliRun.h"
34 #include "AliHijingRndm.h"
35
36 ClassImp(AliGenHijing)
37
38 AliGenHijing::AliGenHijing()
39                  :AliGenMC()
40 {
41 // Constructor
42     fParticles = 0;
43     fHijing    = 0;
44     fDsigmaDb  = 0;
45     fDnDb      = 0;
46     AliHijingRndm::SetHijingRandom(GetRandom());
47 }
48
49 AliGenHijing::AliGenHijing(Int_t npart)
50     :AliGenMC(npart)
51 {
52 // Default PbPb collisions at 5. 5 TeV
53 //
54     fName = "Hijing";
55     fTitle= "Particle Generator using HIJING";
56
57     SetEnergyCMS();
58     SetImpactParameterRange();
59     SetBoostLHC();
60     SetJetEtaRange();
61     SetJetPhiRange();
62     
63     fKeep       =  0;
64     fQuench     =  1;
65     fShadowing  =  1;
66     fTrigger    =  0;
67     fDecaysOff  =  1;
68     fEvaluate   =  0;
69     fSelectAll  =  0;
70     fFlavor     =  0;
71     fSpectators =  1;
72     fDsigmaDb   =  0;
73     fDnDb       =  0;
74     fPtMinJet   = -2.5;         
75     fRadiation  =  3;
76     //
77     SetSimpleJets();
78     SetNoGammas();
79 //
80     fParticles = new TClonesArray("TParticle",10000);    
81 //
82 // Set random number generator   
83     AliHijingRndm::SetHijingRandom(GetRandom());
84     fHijing = 0;
85
86 }
87
88 AliGenHijing::AliGenHijing(const AliGenHijing & hijing):
89     AliGenMC(hijing)
90 {
91 // copy constructor
92 }
93
94
95 AliGenHijing::~AliGenHijing()
96 {
97 // Destructor
98     if ( fDsigmaDb) delete  fDsigmaDb;  
99     if ( fDnDb)     delete  fDnDb;  
100     delete fParticles;
101 }
102
103 void AliGenHijing::Init()
104 {
105 // Initialisation
106     fFrame.Resize(8);
107     fTarget.Resize(8);
108     fProjectile.Resize(8);
109     
110     SetMC(new THijing(fEnergyCMS, fFrame, fProjectile, fTarget, 
111                       fAProjectile, fZProjectile, fATarget, fZTarget, 
112                       fMinImpactParam, fMaxImpactParam));
113
114     fHijing=(THijing*) fMCEvGen;
115     fHijing->SetIHPR2(2,  fRadiation);
116     fHijing->SetIHPR2(3,  fTrigger);
117     fHijing->SetIHPR2(6,  fShadowing);
118     fHijing->SetIHPR2(12, fDecaysOff);    
119     fHijing->SetIHPR2(21, fKeep);
120     fHijing->SetHIPR1(10, fPtMinJet);   
121     fHijing->SetHIPR1(50, fSimpleJet);
122 //
123 //  Quenching
124 //
125 //
126 //  fQuench = 0:  no quenching
127 //  fQuench = 1:  hijing default
128 //  fQuench = 2:  new LHC  parameters for HIPR1(11) and HIPR1(14)
129 //  fQuench = 3:  new RHIC parameters for HIPR1(11) and HIPR1(14)
130 //  fQuench = 4:  new LHC  parameters with log(e) dependence
131 //  fQuench = 5:  new RHIC parameters with log(e) dependence
132     fHijing->SetIHPR2(50, 0);
133     if (fQuench > 0) 
134         fHijing->SetIHPR2(4,  1);
135     else
136         fHijing->SetIHPR2(4,  0);
137 // New LHC parameters from Xin-Nian Wang
138     if (fQuench == 2) {
139         fHijing->SetHIPR1(14, 1.1);
140         fHijing->SetHIPR1(11, 3.7);
141     } else if (fQuench == 3) {
142         fHijing->SetHIPR1(14, 0.20);
143         fHijing->SetHIPR1(11, 2.5);
144     } else if (fQuench == 4) {
145         fHijing->SetIHPR2(50, 1);
146         fHijing->SetHIPR1(14, 4.*0.34);
147         fHijing->SetHIPR1(11, 3.7);
148     } else if (fQuench == 5) {
149         fHijing->SetIHPR2(50, 1);
150         fHijing->SetHIPR1(14, 0.34);
151         fHijing->SetHIPR1(11, 2.5);
152     }
153     
154     
155     AliGenMC::Init();
156     
157 //
158 //  Initialize Hijing  
159 //    
160     fHijing->Initialize();
161 //
162     if (fEvaluate) EvaluateCrossSections();
163 //
164 }
165
166 void AliGenHijing::Generate()
167 {
168 // Generate one event
169
170   Float_t polar[3]    =   {0,0,0};
171   Float_t origin[3]   =   {0,0,0};
172   Float_t origin0[3]  =   {0,0,0};
173   Float_t p[3];
174   Float_t tof;
175
176 //  converts from mm/c to s
177   const Float_t kconv = 0.001/2.999792458e8;
178 //
179   Int_t nt  = 0;
180   Int_t jev = 0;
181   Int_t j, kf, ks, imo;
182   kf = 0;
183     
184
185     
186   fTrials = 0;
187   for (j = 0;j < 3; j++) origin0[j] = fOrigin[j];
188   if(fVertexSmear == kPerEvent) {
189       Vertex();
190       for (j=0; j < 3; j++) origin0[j] = fVertex[j];
191   } 
192
193   while(1)
194   {
195 //    Generate one event
196 // --------------------------------------------------------------------------
197       fSpecn    = 0;  
198       fSpecp    = 0;
199 // --------------------------------------------------------------------------
200       fHijing->GenerateEvent();
201       fTrials++;
202       fHijing->ImportParticles(fParticles,"All");
203       if (fTrigger != kNoTrigger) {
204           if (!CheckTrigger()) continue;
205       }
206       if (fLHC) Boost();
207       
208       
209       Int_t np = fParticles->GetEntriesFast();
210       printf("\n **************************************************%d\n",np);
211       Int_t nc = 0;
212       if (np == 0 ) continue;
213       Int_t i;
214       Int_t* newPos     = new Int_t[np];
215       Int_t* pSelected  = new Int_t[np];
216
217       for (i = 0; i < np; i++) {
218           newPos[i]    = i;
219           pSelected[i] = 0;
220       }
221       
222 //      Get event vertex
223 //
224       TParticle *  iparticle = (TParticle *) fParticles->At(0);
225       fVertex[0] = origin0[0];
226       fVertex[1] = origin0[1];  
227       fVertex[2] = origin0[2];
228       
229 //
230 //      First select parent particles
231 //
232
233       for (i = 0; i < np; i++) {
234           iparticle = (TParticle *) fParticles->At(i);
235
236 // Is this a parent particle ?
237           if (Stable(iparticle)) continue;
238 //
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, 0) && 
248                                SelectFlavor(kf);
249           hasSelectedDaughters = DaughtersSelection(iparticle);
250 //
251 // Put particle on the stack if it is either selected or 
252 // it is the mother of at least one seleted particle
253 //
254           if (selected || hasSelectedDaughters) {
255               nc++;
256               pSelected[i] = 1;
257           } // selected
258       } // particle loop parents
259 //
260 // Now select the final state particles
261 //
262
263       for (i = 0; i<np; i++) {
264           TParticle *  iparticle = (TParticle *) fParticles->At(i);
265 // Is this a final state particle ?
266           if (!Stable(iparticle)) continue;
267       
268           Bool_t  selected             =  kTRUE;
269           kf        = iparticle->GetPdgCode();
270           ks        = iparticle->GetStatusCode();
271           
272 // --------------------------------------------------------------------------
273 // Count spectator neutrons and protons
274           if(ks == 0 || ks == 1 || ks == 10 || ks == 11){
275               if(kf == kNeutron) fSpecn += 1;
276               if(kf == kProton)  fSpecp += 1;
277           }
278 // --------------------------------------------------------------------------
279 //          
280           if (!fSelectAll) {
281               selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
282               if (!fSpectators && selected) selected = (ks != 0 && ks != 1 && ks != 10
283                                                         && ks != 11);
284           }
285 //
286 // Put particle on the stack if selected
287 //
288           if (selected) {
289               nc++;
290               pSelected[i] = 1;
291           } // selected
292       } // particle loop final state
293 //
294 // Write particles to stack
295 //
296       for (i = 0; i<np; i++) {
297           TParticle *  iparticle = (TParticle *) fParticles->At(i);
298           Bool_t  hasMother   = (iparticle->GetFirstMother()     >=0);
299           Bool_t  hasDaughter = (iparticle->GetFirstDaughter()   >=0);
300
301           if (pSelected[i]) {
302               kf   = iparticle->GetPdgCode();
303               ks   = iparticle->GetStatusCode();
304               p[0] = iparticle->Px();
305               p[1] = iparticle->Py();
306               p[2] = iparticle->Pz();
307               origin[0] = origin0[0]+iparticle->Vx()/10;
308               origin[1] = origin0[1]+iparticle->Vy()/10;
309               origin[2] = origin0[2]+iparticle->Vz()/10;
310               tof = kconv*iparticle->T();
311               imo = -1;
312               TParticle* mother = 0;
313               if (hasMother) {
314                   imo = iparticle->GetFirstMother();
315                   mother = (TParticle *) fParticles->At(imo);
316                   imo = (mother->GetPdgCode() != 92) ? imo = newPos[imo] : -1;
317               } // if has mother   
318               Bool_t tFlag = (fTrackIt && !hasDaughter);
319               PushTrack(tFlag,imo,kf,p,origin,polar,
320                        tof,kPNoProcess,nt, 1., ks);
321               KeepTrack(nt);
322               newPos[i] = nt;
323           } // if selected
324       } // particle loop
325       delete[] newPos;
326       delete[] pSelected;
327       
328       printf("\n I've put %i particles on the stack \n",nc);
329       if (nc > 0) {
330           jev += nc;
331           if (jev >= fNpart || fNpart == -1) {
332               fKineBias = Float_t(fNpart)/Float_t(fTrials);
333               printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
334               break;
335           }
336       }
337   } // event loop
338   MakeHeader();
339   SetHighWaterMark(nt);
340 }
341
342 void AliGenHijing::KeepFullEvent()
343 {
344     fKeep=1;
345 }
346
347 void AliGenHijing::EvaluateCrossSections()
348 {
349 //     Glauber Calculation of geometrical x-section
350 //
351     Float_t xTot       = 0.;          // barn
352     Float_t xTotHard   = 0.;          // barn 
353     Float_t xPart      = 0.;          // barn
354     Float_t xPartHard  = 0.;          // barn 
355     Float_t sigmaHard  = 0.1;         // mbarn
356     Float_t bMin       = 0.;
357     Float_t bMax       = fHijing->GetHIPR1(34)+fHijing->GetHIPR1(35);
358     const Float_t kdib = 0.2;
359     Int_t   kMax       = Int_t((bMax-bMin)/kdib)+1;
360
361
362     printf("\n Projectile Radius (fm): %f \n",fHijing->GetHIPR1(34));
363     printf("\n Target     Radius (fm): %f \n",fHijing->GetHIPR1(35));    
364     Int_t i;
365     Float_t oldvalue= 0.;
366
367     Float_t* b   = new Float_t[kMax];
368     Float_t* si1 = new Float_t[kMax];    
369     Float_t* si2 = new Float_t[kMax];    
370     
371     for (i = 0; i < kMax; i++)
372     {
373         Float_t xb  = bMin+i*kdib;
374         Float_t ov;
375         ov=fHijing->Profile(xb);
376         Float_t gb  =  2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*(1.-TMath::Exp(-fHijing->GetHINT1(12)*ov));
377         Float_t gbh =  2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*sigmaHard*ov;
378         xTot+=gb;
379         xTotHard += gbh;
380         printf("profile %f %f %f\n", xb, ov, fHijing->GetHINT1(12));
381         
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)
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 *) fParticles->At(i);   
429             Int_t ip = jparticle->GetPdgCode();
430             if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
431                 selected=kTRUE; break;
432             }
433             if (DaughtersSelection(jparticle)) {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     Bool_t res = 0;
449     
450     if (fFlavor == 0) {
451         res = kTRUE;
452     } else {
453         Int_t ifl = TMath::Abs(pid/100);
454         if (ifl > 10) ifl/=10;
455         res = (fFlavor == ifl);
456     }
457 //
458 //  This part if gamma writing is inhibited
459     if (fNoGammas) 
460         res = res && (pid != kGamma && pid != kPi0);
461 //
462     return res;
463 }
464
465 Bool_t AliGenHijing::Stable(TParticle*  particle) const
466 {
467 // Return true for a stable particle
468 //
469     
470     if (particle->GetFirstDaughter() < 0 )
471     {
472         return kTRUE;
473     } else {
474         return kFALSE;
475     }
476 }
477
478
479
480 void AliGenHijing::MakeHeader()
481 {
482 // Builds the event header, to be called after each event
483     AliGenEventHeader* header = new AliGenHijingEventHeader("Hijing");
484     ((AliGenHijingEventHeader*) header)->SetNProduced(fHijing->GetNATT());
485     ((AliGenHijingEventHeader*) header)->SetImpactParameter(fHijing->GetHINT1(19));
486     ((AliGenHijingEventHeader*) header)->SetTotalEnergy(fHijing->GetEATT());
487     ((AliGenHijingEventHeader*) header)->SetHardScatters(fHijing->GetJATT());
488     ((AliGenHijingEventHeader*) header)->SetParticipants(fHijing->GetNP(), fHijing->GetNT());
489     ((AliGenHijingEventHeader*) header)->SetCollisions(fHijing->GetN0(),
490                                                        fHijing->GetN01(),
491                                                        fHijing->GetN10(),
492                                                        fHijing->GetN11());
493     ((AliGenHijingEventHeader*) header)->SetSpectators(fSpecn, fSpecp);
494
495 // 4-momentum vectors of the triggered jets.
496 //
497 // Before final state gluon radiation.
498     TLorentzVector* jet1 = new TLorentzVector(fHijing->GetHINT1(21), 
499                                               fHijing->GetHINT1(22),
500                                               fHijing->GetHINT1(23),
501                                               fHijing->GetHINT1(24));
502
503     TLorentzVector* jet2 = new TLorentzVector(fHijing->GetHINT1(31), 
504                                               fHijing->GetHINT1(32),
505                                               fHijing->GetHINT1(33),
506                                               fHijing->GetHINT1(34));
507 // After final state gluon radiation.
508     TLorentzVector* jet3 = new TLorentzVector(fHijing->GetHINT1(26), 
509                                               fHijing->GetHINT1(27),
510                                               fHijing->GetHINT1(28),
511                                               fHijing->GetHINT1(29));
512
513     TLorentzVector* jet4 = new TLorentzVector(fHijing->GetHINT1(36), 
514                                               fHijing->GetHINT1(37),
515                                               fHijing->GetHINT1(38),
516                                               fHijing->GetHINT1(39));
517     ((AliGenHijingEventHeader*) header)->SetJets(jet1, jet2, jet3, jet4);
518 // Bookkeeping for kinematic bias
519     ((AliGenHijingEventHeader*) header)->SetTrials(fTrials);
520 // Event Vertex
521     header->SetPrimaryVertex(fVertex);
522     if (gAlice) gAlice->SetGenEventHeader(header);   
523     fCollisionGeometry = (AliGenHijingEventHeader*)  header;
524 }
525
526 Bool_t AliGenHijing::CheckTrigger()
527 {
528 // Check the kinematic trigger condition
529 //
530     Bool_t   triggered = kFALSE;
531  
532     if (fTrigger == 1) {
533 //
534 //  jet-jet Trigger     
535         
536         TLorentzVector* jet1 = new TLorentzVector(fHijing->GetHINT1(26), 
537                                                   fHijing->GetHINT1(27),
538                                                   fHijing->GetHINT1(28),
539                                                   fHijing->GetHINT1(29));
540         
541         TLorentzVector* jet2 = new TLorentzVector(fHijing->GetHINT1(36), 
542                                                   fHijing->GetHINT1(37),
543                                                   fHijing->GetHINT1(38),
544                                                   fHijing->GetHINT1(39));
545         Double_t eta1      = jet1->Eta();
546         Double_t eta2      = jet2->Eta();
547         Double_t phi1      = jet1->Phi();
548         Double_t phi2      = jet2->Phi();
549 //    printf("\n Trigger: %f %f %f %f",
550 //         fEtaMinJet, fEtaMaxJet, fPhiMinJet, fPhiMaxJet);
551         if (
552             (eta1 < fEtaMaxJet && eta1 > fEtaMinJet &&  
553              phi1 < fPhiMaxJet && phi1 > fPhiMinJet) 
554             ||
555             (eta2 < fEtaMaxJet && eta2 > fEtaMinJet &&  
556              phi2 < fPhiMaxJet && phi2 > fPhiMinJet)
557             ) 
558             triggered = kTRUE;
559     } else if (fTrigger == 2) {
560 //  Gamma Jet
561 //
562         Int_t np = fParticles->GetEntriesFast();
563         for (Int_t i = 0; i < np; i++) {
564             TParticle* part = (TParticle*) fParticles->At(i);
565             Int_t kf = part->GetPdgCode();
566             Int_t ks = part->GetStatusCode();
567             if (kf == 22 && ks == 40) {
568                 Float_t phi = part->Phi();
569                 Float_t eta = part->Eta();
570                 if  (eta < fEtaMaxJet && 
571                      eta > fEtaMinJet &&
572                      phi < fPhiMaxJet && 
573                      phi > fPhiMinJet) {
574                     triggered = 1;
575                     break;
576                 } // check phi,eta within limits
577             } // direct gamma ? 
578         } // particle loop
579     } // fTrigger == 2
580     return triggered;
581 }
582
583
584 void AliGenHijing::Copy(TObject &) const
585 {
586   Fatal("Copy","Not implemented!\n");
587 }
588
589 AliGenHijing& AliGenHijing::operator=(const  AliGenHijing& rhs)
590 {
591     rhs.Copy(*this); 
592     return (*this);
593 }
594