- Use Vertex() method
[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 //
22 // andreas.morsch@cern.ch
23
24 #include <TArrayI.h>
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 {
90 // copy constructor
91 }
92
93
94 AliGenHijing::~AliGenHijing()
95 {
96 // Destructor
97     if ( fDsigmaDb) delete  fDsigmaDb;  
98     if ( fDnDb)     delete  fDnDb;  
99     delete fParticles;
100 }
101
102 void AliGenHijing::Init()
103 {
104 // Initialisation
105     fFrame.Resize(8);
106     fTarget.Resize(8);
107     fProjectile.Resize(8);
108     
109     SetMC(new THijing(fEnergyCMS, fFrame, fProjectile, fTarget, 
110                       fAProjectile, fZProjectile, fATarget, fZTarget, 
111                       fMinImpactParam, fMaxImpactParam));
112
113     fHijing=(THijing*) fgMCEvGen;
114     fHijing->SetIHPR2(2,  fRadiation);
115     fHijing->SetIHPR2(3,  fTrigger);
116     fHijing->SetIHPR2(6,  fShadowing);
117     fHijing->SetIHPR2(12, fDecaysOff);    
118     fHijing->SetIHPR2(21, fKeep);
119     fHijing->SetHIPR1(10, fPtMinJet);   
120     fHijing->SetHIPR1(50, fSimpleJet);
121 //
122 //  Quenching
123 //
124 //
125 //  fQuench = 0:  no quenching
126 //  fQuench = 1:  hijing default
127 //  fQuench = 2:  new LHC  parameters for HIPR1(11) and HIPR1(14)
128 //  fQuench = 3:  new RHIC parameters for HIPR1(11) and HIPR1(14)
129 //  fQuench = 4:  new LHC  parameters with log(e) dependence
130 //  fQuench = 5:  new RHIC parameters with log(e) dependence
131     fHijing->SetIHPR2(50, 0);
132     if (fQuench > 0) 
133         fHijing->SetIHPR2(4,  1);
134     else
135         fHijing->SetIHPR2(4,  0);
136 // New LHC parameters from Xin-Nian Wang
137     if (fQuench == 2) {
138         fHijing->SetHIPR1(14, 1.1);
139         fHijing->SetHIPR1(11, 3.7);
140     } else if (fQuench == 3) {
141         fHijing->SetHIPR1(14, 0.20);
142         fHijing->SetHIPR1(11, 2.5);
143     } else if (fQuench == 4) {
144         fHijing->SetIHPR2(50, 1);
145         fHijing->SetHIPR1(14, 4.*0.34);
146         fHijing->SetHIPR1(11, 3.7);
147     } else if (fQuench == 5) {
148         fHijing->SetIHPR2(50, 1);
149         fHijing->SetHIPR1(14, 0.34);
150         fHijing->SetHIPR1(11, 2.5);
151     }
152     
153     
154     AliGenMC::Init();
155     
156 //
157 //  Initialize Hijing  
158 //    
159     fHijing->Initialize();
160 //
161     if (fEvaluate) EvaluateCrossSections();
162 //
163 }
164
165 void AliGenHijing::Generate()
166 {
167 // Generate one event
168
169   Float_t polar[3]    =   {0,0,0};
170   Float_t origin[3]   =   {0,0,0};
171   Float_t origin0[3]  =   {0,0,0};
172   Float_t p[3];
173   Float_t tof;
174
175 //  converts from mm/c to s
176   const Float_t kconv = 0.001/2.999792458e8;
177 //
178   Int_t nt  = 0;
179   Int_t jev = 0;
180   Int_t j, kf, ks, imo;
181   kf = 0;
182     
183
184     
185   fTrials = 0;
186   for (j = 0;j < 3; j++) origin0[j] = fOrigin[j];
187   if(fVertexSmear == kPerEvent) {
188       Vertex();
189       for (j=0; j < 3; j++) origin0[j] = fVertex[j];
190   } 
191
192   while(1)
193   {
194 //    Generate one event
195 // --------------------------------------------------------------------------
196       fSpecn    = 0;  
197       fSpecp    = 0;
198 // --------------------------------------------------------------------------
199       fHijing->GenerateEvent();
200       fTrials++;
201       fHijing->ImportParticles(fParticles,"All");
202       if (fTrigger != kNoTrigger) {
203           if (!CheckTrigger()) continue;
204       }
205       if (fLHC) Boost();
206       
207       
208       Int_t np = fParticles->GetEntriesFast();
209       printf("\n **************************************************%d\n",np);
210       Int_t nc = 0;
211       if (np == 0 ) continue;
212       Int_t i;
213       Int_t* newPos     = new Int_t[np];
214       Int_t* pSelected  = new Int_t[np];
215
216       for (i = 0; i < np; i++) {
217           newPos[i]    = i;
218           pSelected[i] = 0;
219       }
220       
221 //      Get event vertex
222 //
223       TParticle *  iparticle = (TParticle *) fParticles->At(0);
224       fVertex[0] = origin0[0];
225       fVertex[1] = origin0[1];  
226       fVertex[2] = origin0[2];
227       
228 //
229 //      First select parent particles
230 //
231
232       for (i = 0; i < np; i++) {
233           iparticle = (TParticle *) fParticles->At(i);
234
235 // Is this a parent particle ?
236           if (Stable(iparticle)) continue;
237 //
238           Bool_t  selected             =  kTRUE;
239           Bool_t  hasSelectedDaughters =  kFALSE;
240           
241           
242           kf        = iparticle->GetPdgCode();
243           ks        = iparticle->GetStatusCode();
244           if (kf == 92) continue;
245             
246           if (!fSelectAll) selected = KinematicSelection(iparticle, 0) && 
247                                SelectFlavor(kf);
248           hasSelectedDaughters = DaughtersSelection(iparticle);
249 //
250 // Put particle on the stack if it is either selected or 
251 // it is the mother of at least one seleted particle
252 //
253           if (selected || hasSelectedDaughters) {
254               nc++;
255               pSelected[i] = 1;
256           } // selected
257       } // particle loop parents
258 //
259 // Now select the final state particles
260 //
261
262       for (i = 0; i<np; i++) {
263           TParticle *  iparticle = (TParticle *) fParticles->At(i);
264 // Is this a final state particle ?
265           if (!Stable(iparticle)) continue;
266       
267           Bool_t  selected             =  kTRUE;
268           kf        = iparticle->GetPdgCode();
269           ks        = iparticle->GetStatusCode();
270           
271 // --------------------------------------------------------------------------
272 // Count spectator neutrons and protons
273           if(ks == 0 || ks == 1 || ks == 10 || ks == 11){
274               if(kf == kNeutron) fSpecn += 1;
275               if(kf == kProton)  fSpecp += 1;
276           }
277 // --------------------------------------------------------------------------
278 //          
279           if (!fSelectAll) {
280               selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
281               if (!fSpectators && selected) selected = (ks != 0 && ks != 1 && ks != 10
282                                                         && ks != 11);
283           }
284 //
285 // Put particle on the stack if selected
286 //
287           if (selected) {
288               nc++;
289               pSelected[i] = 1;
290           } // selected
291       } // particle loop final state
292 //
293 // Write particles to stack
294 //
295       for (i = 0; i<np; i++) {
296           TParticle *  iparticle = (TParticle *) fParticles->At(i);
297           Bool_t  hasMother   = (iparticle->GetFirstMother()     >=0);
298           Bool_t  hasDaughter = (iparticle->GetFirstDaughter()   >=0);
299
300           if (pSelected[i]) {
301               kf   = iparticle->GetPdgCode();
302               ks   = iparticle->GetStatusCode();
303               p[0] = iparticle->Px();
304               p[1] = iparticle->Py();
305               p[2] = iparticle->Pz();
306               origin[0] = origin0[0]+iparticle->Vx()/10;
307               origin[1] = origin0[1]+iparticle->Vy()/10;
308               origin[2] = origin0[2]+iparticle->Vz()/10;
309               tof = kconv*iparticle->T();
310               imo = -1;
311               TParticle* mother = 0;
312               if (hasMother) {
313                   imo = iparticle->GetFirstMother();
314                   mother = (TParticle *) fParticles->At(imo);
315                   imo = (mother->GetPdgCode() != 92) ? imo = newPos[imo] : -1;
316               } // if has mother   
317               Bool_t tFlag = (fTrackIt && !hasDaughter);
318               SetTrack(tFlag,imo,kf,p,origin,polar,
319                        tof,kPNoProcess,nt, 1., ks);
320               KeepTrack(nt);
321               newPos[i] = nt;
322           } // if selected
323       } // particle loop
324       delete[] newPos;
325       delete[] pSelected;
326       
327       printf("\n I've put %i particles on the stack \n",nc);
328       if (nc > 0) {
329           jev += nc;
330           if (jev >= fNpart || fNpart == -1) {
331               fKineBias = Float_t(fNpart)/Float_t(fTrials);
332               printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
333               break;
334           }
335       }
336   } // event loop
337   MakeHeader();
338   SetHighWaterMark(nt);
339 }
340
341 void AliGenHijing::KeepFullEvent()
342 {
343     fKeep=1;
344 }
345
346 void AliGenHijing::EvaluateCrossSections()
347 {
348 //     Glauber Calculation of geometrical x-section
349 //
350     Float_t xTot       = 0.;          // barn
351     Float_t xTotHard   = 0.;          // barn 
352     Float_t xPart      = 0.;          // barn
353     Float_t xPartHard  = 0.;          // barn 
354     Float_t sigmaHard  = 0.1;         // mbarn
355     Float_t bMin       = 0.;
356     Float_t bMax       = fHijing->GetHIPR1(34)+fHijing->GetHIPR1(35);
357     const Float_t kdib = 0.2;
358     Int_t   kMax       = Int_t((bMax-bMin)/kdib)+1;
359
360
361     printf("\n Projectile Radius (fm): %f \n",fHijing->GetHIPR1(34));
362     printf("\n Target     Radius (fm): %f \n",fHijing->GetHIPR1(35));    
363     Int_t i;
364     Float_t oldvalue= 0.;
365
366     Float_t* b   = new Float_t[kMax];
367     Float_t* si1 = new Float_t[kMax];    
368     Float_t* si2 = new Float_t[kMax];    
369     
370     for (i = 0; i < kMax; i++)
371     {
372         Float_t xb  = bMin+i*kdib;
373         Float_t ov;
374         ov=fHijing->Profile(xb);
375         Float_t gb  =  2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*(1.-TMath::Exp(-fHijing->GetHINT1(12)*ov));
376         Float_t gbh =  2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*sigmaHard*ov;
377         xTot+=gb;
378         xTotHard += gbh;
379         printf("profile %f %f %f\n", xb, ov, fHijing->GetHINT1(12));
380         
381         if (xb > fMinImpactParam && xb < fMaxImpactParam)
382         {
383             xPart += gb;
384             xPartHard += gbh;
385         }
386         
387         if(oldvalue) if ((xTot-oldvalue)/oldvalue<0.0001) break;
388         oldvalue = xTot;
389         printf("\n Total cross section (barn): %d %f %f \n",i, xb, xTot);
390         printf("\n Hard  cross section (barn): %d %f %f \n\n",i, xb, xTotHard);
391         if (i>0) {
392             si1[i] = gb/kdib;
393             si2[i] = gbh/gb;
394             b[i]  = xb;
395         }
396     }
397
398     printf("\n Total cross section (barn): %f \n",xTot);
399     printf("\n Hard  cross section (barn): %f \n \n",xTotHard);
400     printf("\n Partial       cross section (barn): %f %f \n",xPart, xPart/xTot*100.);
401     printf("\n Partial  hard cross section (barn): %f %f \n",xPartHard, xPartHard/xTotHard*100.);
402
403 //  Store result as a graph
404     b[0] = 0;
405     si1[0] = 0;
406     si2[0]=si2[1];
407     
408     fDsigmaDb  = new TGraph(i, b, si1);
409     fDnDb      = new TGraph(i, b, si2);
410 }
411
412 Bool_t AliGenHijing::DaughtersSelection(TParticle* iparticle)
413 {
414 //
415 // Looks recursively if one of the daughters has been selected
416 //
417 //    printf("\n Consider daughters %d:",iparticle->GetPdgCode());
418     Int_t imin = -1;
419     Int_t imax = -1;
420     Int_t i;
421     Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
422     Bool_t selected = kFALSE;
423     if (hasDaughters) {
424         imin = iparticle->GetFirstDaughter();
425         imax = iparticle->GetLastDaughter();       
426         for (i = imin; i <= imax; i++){
427             TParticle *  jparticle = (TParticle *) fParticles->At(i);   
428             Int_t ip = jparticle->GetPdgCode();
429             if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
430                 selected=kTRUE; break;
431             }
432             if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
433         }
434     } else {
435         return kFALSE;
436     }
437     return selected;
438 }
439
440
441 Bool_t AliGenHijing::SelectFlavor(Int_t pid)
442 {
443 // Select flavor of particle
444 // 0: all
445 // 4: charm and beauty
446 // 5: beauty
447     Bool_t res = 0;
448     
449     if (fFlavor == 0) {
450         res = kTRUE;
451     } else {
452         Int_t ifl = TMath::Abs(pid/100);
453         if (ifl > 10) ifl/=10;
454         res = (fFlavor == ifl);
455     }
456 //
457 //  This part if gamma writing is inhibited
458     if (fNoGammas) 
459         res = res && (pid != kGamma && pid != kPi0);
460 //
461     return res;
462 }
463
464 Bool_t AliGenHijing::Stable(TParticle*  particle)
465 {
466 // Return true for a stable particle
467 //
468     
469     if (particle->GetFirstDaughter() < 0 )
470     {
471         return kTRUE;
472     } else {
473         return kFALSE;
474     }
475 }
476
477
478
479 void AliGenHijing::MakeHeader()
480 {
481 // Builds the event header, to be called after each event
482     AliGenEventHeader* header = new AliGenHijingEventHeader("Hijing");
483     ((AliGenHijingEventHeader*) header)->SetNProduced(fHijing->GetNATT());
484     ((AliGenHijingEventHeader*) header)->SetImpactParameter(fHijing->GetHINT1(19));
485     ((AliGenHijingEventHeader*) header)->SetTotalEnergy(fHijing->GetEATT());
486     ((AliGenHijingEventHeader*) header)->SetHardScatters(fHijing->GetJATT());
487     ((AliGenHijingEventHeader*) header)->SetParticipants(fHijing->GetNP(), fHijing->GetNT());
488     ((AliGenHijingEventHeader*) header)->SetCollisions(fHijing->GetN0(),
489                                                        fHijing->GetN01(),
490                                                        fHijing->GetN10(),
491                                                        fHijing->GetN11());
492     ((AliGenHijingEventHeader*) header)->SetSpectators(fSpecn, fSpecp);
493
494 // 4-momentum vectors of the triggered jets.
495 //
496 // Before final state gluon radiation.
497     TLorentzVector* jet1 = new TLorentzVector(fHijing->GetHINT1(21), 
498                                               fHijing->GetHINT1(22),
499                                               fHijing->GetHINT1(23),
500                                               fHijing->GetHINT1(24));
501
502     TLorentzVector* jet2 = new TLorentzVector(fHijing->GetHINT1(31), 
503                                               fHijing->GetHINT1(32),
504                                               fHijing->GetHINT1(33),
505                                               fHijing->GetHINT1(34));
506 // After final state gluon radiation.
507     TLorentzVector* jet3 = new TLorentzVector(fHijing->GetHINT1(26), 
508                                               fHijing->GetHINT1(27),
509                                               fHijing->GetHINT1(28),
510                                               fHijing->GetHINT1(29));
511
512     TLorentzVector* jet4 = new TLorentzVector(fHijing->GetHINT1(36), 
513                                               fHijing->GetHINT1(37),
514                                               fHijing->GetHINT1(38),
515                                               fHijing->GetHINT1(39));
516     ((AliGenHijingEventHeader*) header)->SetJets(jet1, jet2, jet3, jet4);
517 // Bookkeeping for kinematic bias
518     ((AliGenHijingEventHeader*) header)->SetTrials(fTrials);
519 // Event Vertex
520     header->SetPrimaryVertex(fVertex);
521     gAlice->SetGenEventHeader(header);   
522     fCollisionGeometry = (AliGenHijingEventHeader*)  header;
523 }
524
525 Bool_t AliGenHijing::CheckTrigger()
526 {
527 // Check the kinematic trigger condition
528 //
529     Bool_t   triggered = kFALSE;
530  
531     if (fTrigger == 1) {
532 //
533 //  jet-jet Trigger     
534         
535         TLorentzVector* jet1 = new TLorentzVector(fHijing->GetHINT1(26), 
536                                                   fHijing->GetHINT1(27),
537                                                   fHijing->GetHINT1(28),
538                                                   fHijing->GetHINT1(29));
539         
540         TLorentzVector* jet2 = new TLorentzVector(fHijing->GetHINT1(36), 
541                                                   fHijing->GetHINT1(37),
542                                                   fHijing->GetHINT1(38),
543                                                   fHijing->GetHINT1(39));
544         Double_t eta1      = jet1->Eta();
545         Double_t eta2      = jet2->Eta();
546         Double_t phi1      = jet1->Phi();
547         Double_t phi2      = jet2->Phi();
548 //    printf("\n Trigger: %f %f %f %f",
549 //         fEtaMinJet, fEtaMaxJet, fPhiMinJet, fPhiMaxJet);
550         if (
551             (eta1 < fEtaMaxJet && eta1 > fEtaMinJet &&  
552              phi1 < fPhiMaxJet && phi1 > fPhiMinJet) 
553             ||
554             (eta2 < fEtaMaxJet && eta2 > fEtaMinJet &&  
555              phi2 < fPhiMaxJet && phi2 > fPhiMinJet)
556             ) 
557             triggered = kTRUE;
558     } else if (fTrigger == 2) {
559 //  Gamma Jet
560 //
561         Int_t np = fParticles->GetEntriesFast();
562         for (Int_t i = 0; i < np; i++) {
563             TParticle* part = (TParticle*) fParticles->At(i);
564             Int_t kf = part->GetPdgCode();
565             Int_t ks = part->GetStatusCode();
566             if (kf == 22 && ks == 40) {
567                 Float_t phi = part->Phi();
568                 Float_t eta = part->Eta();
569                 if  (eta < fEtaMaxJet && 
570                      eta > fEtaMinJet &&
571                      phi < fPhiMaxJet && 
572                      phi > fPhiMinJet) {
573                     triggered = 1;
574                     break;
575                 } // check phi,eta within limits
576             } // direct gamma ? 
577         } // particle loop
578     } // fTrigger == 2
579     return triggered;
580 }
581
582
583
584
585 AliGenHijing& AliGenHijing::operator=(const  AliGenHijing& rhs)
586 {
587 // Assignment operator
588     return *this;
589 }
590