Corrected library names and paths to macros
[u/mrichter/AliRoot.git] / test / genkine / gen / fastMcProduction.C
1 /***************************************************************************
2  *  Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                         *
4  *  Author: ALICE OFFLICE.                                                 *
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      $Satyajit Jena || alien:sjena Sun Apr 21 14:05:19 CEST 2013$
17                                                                           
18        Sterring macro for fast production of MC events - Followed from 
19        the original macro of fastGenAmpt.C     
20      
21      Implemented Generators: (Version 1.0: Sun Apr 21 14:05:19 CEST 2013)
22      ----------------------------------------------------------------
23       kPythia6,            kPythia8,               kPythia6D6T,       
24       kPythiaPerugia0,     kPythia6ATLAS,          
25       kPythiaJets,         
26       kPhojet,                
27       kDPMjet,             kDPMjet_pA, 
28       kAmptDefault,        kAmptStringMelting,          kAmptStringMeltingNoART,
29       kAmptpA,                kAmptFlow,
30       kAmptReducedFlow,
31
32      FIXME: 
33      kPythia6ATLAS_Flat, 
34      kHijing,             kHijing2000,            kHijing_pA,             
35
36                                                                 
37  ***************************************************************************/
38
39 #if !defined(__CINT__) || defined(__MAKECINT__)
40 #include <Riostream.h>
41 #include <TStopwatch.h>
42 #include <TDatime.h>
43 #include <TRandom.h>
44 #include <TFile.h>
45 #include <TTree.h>
46 #include <TMath.h>
47 #include <TParticle.h>
48 #include <TPDGCode.h>
49 #include <TDatabasePDG.h>
50 #include <TRandom3.h>
51 #include <TChain.h>
52 #include "AliRun.h"
53 #include "AliRunLoader.h"
54 #include "AliHeader.h"
55 #include "AliStack.h"
56 #include "AliPDG.h"
57 #include "AliGenAmpt.h"
58 #include "TAmpt.h"
59 #include "PYTHIA6/AliDecayerPythia.h"
60 #include <TVirtualMC.h>
61 #include <TF1.h>
62 #include "STEER/AliConfig.h"
63 #include "STEER/AliGenerator.h"
64 #include "STEER/AliLog.h"
65 #include "EVGEN/AliGenHIJINGpara.h"
66 #include "THijing/AliGenHijing.h"
67 #include "EVGEN/AliGenCocktail.h"
68 #include "EVGEN/AliGenSlowNucleons.h"
69 #include "EVGEN/AliSlowNucleonModelExp.h"
70 #include "EVGEN/AliGenParam.h"
71 #include "EVGEN/AliGenMUONlib.h"
72 #include "EVGEN/AliGenSTRANGElib.h"
73 #include "EVGEN/AliGenMUONCocktail.h"
74 #include "EVGEN/AliGenCocktail.h"
75 #include "EVGEN/AliGenGeVSim.h"
76 #include "EVGEN/AliGeVSimParticle.h"
77 #include "PYTHIA6/AliGenPythia.h"
78 #endif
79
80 //___________________________________________________________________________
81 enum PDC06Proc_t { 
82   kPythia6,
83   kPythia8, 
84   kPythia6D6T, 
85   kPythiaPerugia0, 
86   kPythia6ATLAS, 
87   kPythia6ATLAS_Flat, 
88   kPythiaJets,
89   kPhojet, 
90   kHijing, 
91   kHijing2000, 
92   kHijing_pA,
93   kDPMjet, 
94   kDPMjet_pA,
95   kAmptDefault,
96   kAmpt, 
97   kAmptpA,
98   kAmptFlowStringMelting,
99   kAmptStringMeltingNoART,
100   kAmptFlow,
101   kAmptReducedFlow,
102   kRunMax
103 };
104
105 //___________________________________________________________________________
106 const char * pprRunName[] = {
107   "kPythia6", 
108   "kPythia8",
109   "kPythia6D6T", 
110   "kPythiaPerugia0", 
111   "kPythia6ATLAS", 
112   "kPythia6ATLAS_Flat", 
113   "kPythiaJets", 
114   "kPhojet",  
115   "kHijing", 
116   "kHijing2000", 
117   "kHijing_pA",
118   "kDPMjet", 
119   "kDPMjet_pA", 
120   "kAmpt",
121   "kAmptpA", 
122   "kAmptFlow",
123   "kAmptReducedFlow"  
124 };
125
126 enum PprTrigConf_t {kDefaultPPTrig, kDefaultPbPbTrig };
127 const char * pprTrigConfName[] = {"p-p","Pb-Pb"};
128
129 //___________________________________________________________________________
130 void ProcessEnvironmentVars();
131 class AliGenPythia;
132 AliGenerator *MbPythia();
133 AliGenerator *Pythia8();
134 AliGenerator *MbPythiaTuneD6T();
135 AliGenerator *MbPythiaTunePerugia0();
136 AliGenerator *MbPythiaTuneATLAS();
137 AliGenerator *MbPythiaTuneATLAS_Flat();
138 AliGenerator *PythiaJets();
139 AliGenerator *MbPhojet();
140 AliGenerator *Hijing();
141 AliGenerator *Hijing2000();
142 AliGenerator *Hijing_pA(Bool_t kSlowN);
143 AliGenerator *DPMjet();
144 AliGenerator *DPMjet_pA(Bool_t fragments);
145 AliGenerator *Ampt();
146 AliGenerator *AmptpA();
147 AliGenerator* AmptFlow();
148 AliGenerator *AmptReducedFlow();
149
150
151 //_________________________________________________________________________
152 //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
153 // Geterator, field, beam energy
154
155 //static Double_t    pBeamEnergy = 4000.0;  // Used during pA runs
156 //static Double_t  energy  = 2.*pBeamEnergy*TMath::Sqrt(82./208.); //energy in CMS 
157
158 static Double_t      energy    = 7000.;
159 static PDC06Proc_t   proc      = kPythia6;
160 static Float_t       bMin      = 0.;
161 static Float_t       bMax      = 100.;
162 static PprTrigConf_t strig     = kDefaultPPTrig; // default pp trigger configuration
163
164
165 static Double_t  JpsiPol      = 0; // Jpsi polarisation
166 static Bool_t    JpsiHarderPt = kFALSE; // Jpsi harder pt spectrum (8.8 TeV)
167 static TString comment;
168 //static PprTrigConf_t strig    = kDefaultPbPbTrig; // default pp trigger configuration
169 TDatime dt; 
170 static UInt_t seed    = dt.Get();
171
172 //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
173 //_________________________________________________________________________
174 void fastMcProduction(Int_t nev = 300) {
175
176   ProcessEnvironmentVars();
177
178   gRandom->SetSeed(seed);
179   cerr<<" Seed for random number generation = "<<seed<<endl; 
180   
181 #if defined(__CINT__)
182   gSystem->Load("liblhapdf");  
183   gSystem->Load("libEGPythia6"); 
184   
185   if (proc == kPythia6 || proc == kPhojet || proc == kDPMjet || proc==kDPMjet_pA) {
186     gSystem->Load("libpythia6");        // Pythia 6.2
187     gSystem->Load("libAliPythia6");     // ALICE specific implementations
188   }
189     
190   if (proc == kHijing || proc == kHijing2000 || proc == kHijing_pA ) {
191     gSystem->Load("libHIJING"); 
192     gSystem->Load("libTHijing");
193   } 
194   
195   else if ( proc == kDPMjet || proc== kDPMjet_pA ) {
196     gSystem->Load("libDPMJET"); 
197     gSystem->Load("libTDPMjet");
198   } 
199   
200   else if (proc == kAmptDefault || kAmptFlowStringMelting || proc ==  kAmptStringMeltingNoART || proc == kAmptpA || proc == kAmptReducedFlow) {
201     gSystem->Load("libampt");  
202     gSystem->Load("libTAmpt");
203     gSystem->Load("libpythia6.so");
204     gSystem->Load("libAliPythia6.so");
205   } 
206
207   if (proc == kPythia8) {
208     gSystem->Load("libpythia8.so");
209     gSystem->Load("libAliPythia8.so");
210     gSystem->Setenv("PYTHIA8DATA", gSystem->ExpandPathName("$ALICE_ROOT/PYTHIA8/pythia8145/xmldoc"));
211     gSystem->Setenv("LHAPDF",      gSystem->ExpandPathName("$ALICE_ROOT/LHAPDF"));
212     gSystem->Setenv("LHAPATH",     gSystem->ExpandPathName("$ALICE_ROOT/LHAPDF/PDFsets"));
213   }
214 #endif
215
216
217  AliGenerator* gener = 0x0;
218  
219  cout<<"Run type set to ------------- "<<pprRunName[proc]<<"   " << proc << "    " << kDPMjet_pA<< endl;
220
221  if (proc == kPythia6) {
222    Printf("<<<<<<<<<<<<<<<<<<< Processing Mb. PYTHIA >>>>>>>>>>>>>>>>>>>>"); 
223    gener = MbPythia();
224  } 
225  
226  else if (proc == kPythia8) {
227    Printf("<<<<<<<<<<<<<<<<<<< Processing Mb. Pythia 8 >>>>>>>>>>>>>>>>>>>>"); 
228    gener = Pythia8();
229  }  
230  
231  else if (proc == kPythia6D6T) {
232    Printf("<<<<<<<<<<<<<<<<<<< Processing Mb. PYTHIA D6T >>>>>>>>>>>>>>>>>>>>"); 
233    gener = MbPythiaTuneD6T();
234  } 
235
236  else if (proc == kPythiaPerugia0) {
237    Printf("<<<<<<<<<<<<<<<<<<< Processing Mb. PYTHIA Perugia0 >>>>>>>>>>>>>>>>>>>>"); 
238    gener = MbPythiaTunePerugia0();
239  } 
240  
241  else if (proc == kPythia6ATLAS) {
242    Printf("<<<<<<<<<<<<<<<<<<< Processing Mb. PYTHIA ATLAS>>>>>>>>>>>>>>>>>>>>"); 
243    gener = MbPythiaTuneATLAS();
244  } 
245  
246  else if (proc == kPythia6ATLAS_Flat) {
247    Printf("<<<<<<<<<<<<<<<<<<< Processing Mb. PYTHIA ATLAS_FLAT >>>>>>>>>>>>>>>>>>>>"); 
248    gener = MbPythiaTuneATLAS_Flat();
249  } 
250  
251  else if (proc == kPythiaJets ) {
252    Printf("<<<<<<<<<<<<<<<<<<< Processing Pythia Jets >>>>>>>>>>>>>>>>>>>>"); 
253    gener = PythiaJets();
254  } 
255  
256  else if (proc == kPhojet) {
257    Printf("<<<<<<<<<<<<<<<<<<< Processing Mb. PHOJET >>>>>>>>>>>>>>>>>>>>"); 
258    gener = MbPhojet();
259  } 
260
261  else if (proc == kHijing) {
262    Printf("<<<<<<<<<<<<<<<<<<< Processing Mb. HIJING >>>>>>>>>>>>>>>>>>>>"); 
263    gener = Hijing();    
264  } 
265  
266  else if (proc == kHijing2000) {
267    Printf("<<<<<<<<<<<<<<<<<<< Processing Mb. HIJING 2000 >>>>>>>>>>>>>>>>>>>>"); 
268    gener = Hijing2000();        
269  }
270  
271  else if (proc ==kHijing_pA) {
272    Printf("<<<<<<<<<<<<<<<<<<< Processing Mb. pA Hijing >>>>>>>>>>>>>>>>>>>>"); 
273    gener = Hijing_pA(kTRUE);
274  } 
275  
276  else if (proc == kDPMjet) {
277    Printf("<<<<<<<<<<<<<<<<<<< Processing  DMPJet  >>>>>>>>>>>>>>>>>>>>");
278    gener = DPMjet();    
279  } 
280
281  else if (proc == kDPMjet_pA) {
282    Printf("<<<<<<<<<<<<<<<<<<< Processing  DMPJet pA >>>>>>>>>>>>>>>>>>>>");
283    gener = DPMjet_pA(kFALSE);   
284  } 
285  
286  else if (proc == kAmptDefault) {
287    Printf("<<<<<<<<<<<<<<<<<<< Processing  AMPT Default >>>>>>>>>>>>>>>>>>>>");
288    gener = AmptDefault();
289  } 
290
291  else if (proc == kAmptStringMelting) {
292      Printf("<<<<<<<<<<<<<<<<<<< Processing  AMPT With Flow  >>>>>>>>>>>>>>>>>>>>");
293      gener = AmptStringMelting();
294  }
295
296  else if (proc == kAmptStringMeltingNoART) {
297      Printf("<<<<<<<<<<<<<<<<<<< Processing  AMPT With Flow  >>>>>>>>>>>>>>>>>>>>");
298      gener = AmptStringMeltingNoART();
299  }
300
301  else if (proc == kAmptpA) {
302    Printf("<<<<<<<<<<<<<<<<<<< Processing  AMPT pA  >>>>>>>>>>>>>>>>>>>>");
303    gener = AmptpA();
304  } 
305  
306
307  else if (proc == kAmptReducedFlow) {
308    // Specific Fastgen
309    Printf("<<<<<<<<<<<<<<<<<<< Processing  AMPT With Reduced Flow >>>>>>>>>>>>>>>>>>>>");
310    gener = AmptReducedFlow();
311  } 
312
313  else {
314    cout << "ERROR : Wrong Procss Selcted !!!" << endl;
315    return;
316  }
317
318
319  AliPDG::AddParticlesToPdgDataBase();
320  TDatabasePDG::Instance();
321  
322  const char* filename = "galice.root";
323  AliRunLoader* rl = AliRunLoader::Open(filename,"FASTRUN","recreate");
324  
325  rl->SetCompressionLevel(2);
326  rl->SetNumberOfEventsPerFile(nev);
327  rl->LoadKinematics("RECREATE");
328  rl->MakeTree("E");
329  gAlice->SetRunLoader(rl);
330  rl->MakeStack();
331  AliStack* stack = rl->Stack();
332  
333  AliHeader* header = rl->GetHeader();
334  
335  /*
336    Float_t sigmaz  = 5.4 / TMath::Sqrt(2.); // [cm]
337    Float_t betast  = 3.5;                      // beta* [m]
338    Float_t eps     = 3.75e-6;                   // emittance [m]
339    Float_t gamma   = energy / 2.0 / 0.938272;  // relativistic gamma [1]
340    Float_t sigmaxy = TMath::Sqrt(eps * betast / gamma) / TMath::Sqrt(2.) * 100.;  // [cm]
341    printf("\n \n Diamond size x-y: %10.3e z: %10.3e\n \n", sigmaxy, sigmaz);
342    gener->SetSigma(sigmaxy, sigmaxy, sigmaz);      // Sigma in (X,Y,Z) (cm) on IP position
343    gener->SetVertexSmear(kPerEvent);
344  */
345  
346
347  gener->Init();
348  gener->SetStack(stack);
349  
350  rl->CdGAFile();
351  
352  TStopwatch timer;
353  timer.Start();
354  for (Int_t iev = 0; iev < nev; iev++) {
355    cout <<"============================================= Event number "<< iev << endl;
356    //  Initialize event
357    header->Reset(0,iev);
358    rl->SetEventNumber(iev);
359    stack->Reset();
360    rl->MakeTree("K");
361    Int_t nprim = 0;
362    Int_t ntrial = 0;
363    //  Int_t ndstar = 0;
364    stack->Reset();
365    stack->ConnectTree(rl->TreeK());
366    gener->Generate();
367    ntrial++;
368    nprim = stack->GetNprimary();
369    cout << "Number of particles " << nprim << endl;
370    cout << "Number of trials " << ntrial << endl;
371    header->SetNprimary(stack->GetNprimary());
372    header->SetNtrack(stack->GetNtrack());  
373    stack->FinishEvent();
374     header->SetStack(stack);
375     rl->TreeE()->Fill();
376     rl->WriteKinematics("OVERWRITE");
377  } // event loop
378   timer.Stop();
379   timer.Print();
380   gener->FinishRun();
381   rl->WriteHeader("OVERWRITE");
382   gener->Write();
383   rl->Write();
384 }
385
386 //___________________________________________________//
387 void ProcessEnvironmentVars() {
388     // Run type
389     if (gSystem->Getenv("CONFIG_RUN_TYPE")) {
390       for (Int_t iRun = 0; iRun < kRunMax; iRun++) {
391         if (strcmp(gSystem->Getenv("CONFIG_RUN_TYPE"), pprRunName[iRun]) == 0) {
392           proc = (PDC06Proc_t)iRun;
393           cout<<"Run type set to "<<pprRunName[iRun]<<endl;
394         }
395       }
396     }
397
398     
399     // Energy
400     if (gSystem->Getenv("CONFIG_ENERGY")) {
401       energy = atoi(gSystem->Getenv("CONFIG_ENERGY"));
402       cout<<"Energy set to "<<energy<<" GeV"<<endl;
403     }
404
405     // Random Number seed
406     if (gSystem->Getenv("CONFIG_SEED")) {
407       seed = atoi(gSystem->Getenv("CONFIG_SEED"));
408     }
409
410     // Impact param
411     if (gSystem->Getenv("CONFIG_BMIN")) {
412       bMin = atof(gSystem->Getenv("CONFIG_BMIN"));
413     }
414
415     if (gSystem->Getenv("CONFIG_BMAX")) {
416       bMax = atof(gSystem->Getenv("CONFIG_BMAX"));
417     }
418     cout<<"Impact parameter in ["<<bMin<<","<<bMax<<"]"<<endl;
419 }
420
421
422
423 //______________________________________________________________________
424 AliGenerator* MbPythia() // Mb Pythia
425 {
426       comment = comment.Append(" pp: Pythia low-pt");
427       AliGenPythia* pythia = new AliGenPythia(-1); 
428       /* pythia->SetMomentumRange(0, 999999.);
429       pythia->SetThetaRange(0., 180.);
430       pythia->SetYRange(-12.,12.);
431       pythia->SetPtRange(0,1000.);*/
432       pythia->SetProcess(kPyMb);
433       pythia->SetEnergyCMS(energy);
434       
435       return pythia;
436 }
437
438
439 //______________________________________________________________________
440 AliGenerator* Pythia8()
441 {
442   AliGenPythiaPlus* gener = new AliGenPythiaPlus(AliPythia8::Instance());
443   gener->SetProcess(kPyMbDefault);
444   gener->SetEnergyCMS(energy);
445   gener->SetEventListRange(-1, 2);
446   return gener;
447 }
448
449
450
451 //______________________________________________________________________
452 AliGenerator* MbPythiaTuneD6T()
453 {
454       comment = comment.Append(" pp: Pythia low-pt");
455       AliGenPythia* pythia = new AliGenPythia(-1); 
456       pythia->SetMomentumRange(0, 999999.);
457       pythia->SetThetaRange(0., 180.);
458       pythia->SetYRange(-12.,12.);
459       pythia->SetPtRange(0,1000.);
460       pythia->SetProcess(kPyMb);
461       pythia->SetEnergyCMS(energy);
462 //    Tune
463 //    109     D6T : Rick Field's CDF Tune D6T (NB: needs CTEQ6L pdfs externally)
464       pythia->SetTune(109); // F I X 
465       pythia->SetStrucFunc(kCTEQ6l);
466 //
467       return pythia;
468 }
469
470 //______________________________________________________________________
471 AliGenerator* MbPythiaTunePerugia0()
472 {
473       comment = comment.Append(" pp: Pythia low-pt (Perugia0)");
474 //
475 //    Pythia
476       AliGenPythia* pythia = new AliGenPythia(-1); 
477       /* pythia->SetMomentumRange(0, 999999.);
478       pythia->SetThetaRange(0., 180.);
479       pythia->SetYRange(-12.,12.);
480       pythia->SetPtRange(0,1000.);*/
481       pythia->SetProcess(kPyMb);
482       pythia->SetEnergyCMS(energy);
483 //    Tune
484 //    320     Perugia 0
485       pythia->SetTune(320); 
486       pythia->UseNewMultipleInteractionsScenario();
487 //
488       return pythia;
489 }
490
491 //______________________________________________________________________
492 AliGenerator* MbPythiaTuneATLAS()
493 {
494       comment = comment.Append(" pp: Pythia low-pt");
495 //
496 //    Pythia
497       AliGenPythia* pythia = new AliGenPythia(-1); 
498       /*   pythia->SetMomentumRange(0, 999999.);
499       pythia->SetThetaRange(0., 180.);
500       pythia->SetYRange(-12.,12.);
501       pythia->SetPtRange(0,1000.);*/
502       pythia->SetProcess(kPyMb);
503       pythia->SetEnergyCMS(energy);
504 //    Tune
505 //    C   306 ATLAS-CSC: Arthur Moraes' (new) ATLAS tune (needs CTEQ6L externally)
506       pythia->SetTune(306);
507       pythia->SetStrucFunc(kCTEQ6l);
508 //
509       return pythia;
510 }
511
512
513 //______________________________________________________________________
514 AliGenerator* MbPythiaTuneATLAS_Flat()
515 {
516   AliGenPythia* pythia = MbPythiaTuneATLAS();
517   
518   comment = comment.Append("; flat multiplicity distribution");
519   
520   // set high multiplicity trigger
521   // this weight achieves a flat multiplicity distribution
522   TH1 *weight = new TH1D("weight","weight",201,-0.5,200.5);
523   weight->SetBinContent(1,5.49443);
524   weight->SetBinContent(2,8.770816);
525   weight->SetBinContent(6,0.4568624);
526   weight->SetBinContent(7,0.2919915);
527   weight->SetBinContent(8,0.6674189);
528   weight->SetBinContent(9,0.364737);
529   weight->SetBinContent(10,0.8818444);
530   weight->SetBinContent(11,0.531885);
531   weight->SetBinContent(12,1.035197);
532   weight->SetBinContent(13,0.9394057);
533   weight->SetBinContent(14,0.9643193);
534   weight->SetBinContent(15,0.94543);
535   weight->SetBinContent(16,0.9426507);
536   weight->SetBinContent(17,0.9423649);
537   weight->SetBinContent(18,0.789456);
538   weight->SetBinContent(19,1.149026);
539   weight->SetBinContent(20,1.100491);
540   weight->SetBinContent(21,0.6350525);
541   weight->SetBinContent(22,1.351941);
542   weight->SetBinContent(23,0.03233504);
543   weight->SetBinContent(24,0.9574557);
544   weight->SetBinContent(25,0.868133);
545   weight->SetBinContent(26,1.030998);
546   weight->SetBinContent(27,1.08897);
547   weight->SetBinContent(28,1.251382);
548   weight->SetBinContent(29,0.1391099);
549   weight->SetBinContent(30,1.192876);
550   weight->SetBinContent(31,0.448944);
551   weight->SetBinContent(32,1);
552   weight->SetBinContent(33,1);
553   weight->SetBinContent(34,1);
554   weight->SetBinContent(35,1);
555   weight->SetBinContent(36,0.9999997);
556   weight->SetBinContent(37,0.9999997);
557   weight->SetBinContent(38,0.9999996);
558   weight->SetBinContent(39,0.9999996);
559   weight->SetBinContent(40,0.9999995);
560   weight->SetBinContent(41,0.9999993);
561   weight->SetBinContent(42,1);
562   weight->SetBinContent(43,1);
563   weight->SetBinContent(44,1);
564   weight->SetBinContent(45,1);
565   weight->SetBinContent(46,1);
566   weight->SetBinContent(47,0.9999999);
567   weight->SetBinContent(48,0.9999998);
568   weight->SetBinContent(49,0.9999998);
569   weight->SetBinContent(50,0.9999999);
570   weight->SetBinContent(51,0.9999999);
571   weight->SetBinContent(52,0.9999999);
572   weight->SetBinContent(53,0.9999999);
573   weight->SetBinContent(54,0.9999998);
574   weight->SetBinContent(55,0.9999998);
575   weight->SetBinContent(56,0.9999998);
576   weight->SetBinContent(57,0.9999997);
577   weight->SetBinContent(58,0.9999996);
578   weight->SetBinContent(59,0.9999995);
579   weight->SetBinContent(60,1);
580   weight->SetBinContent(61,1);
581   weight->SetBinContent(62,1);
582   weight->SetBinContent(63,1);
583   weight->SetBinContent(64,1);
584   weight->SetBinContent(65,0.9999999);
585   weight->SetBinContent(66,0.9999998);
586   weight->SetBinContent(67,0.9999998);
587   weight->SetBinContent(68,0.9999999);
588   weight->SetBinContent(69,1);
589   weight->SetBinContent(70,1);
590   weight->SetBinContent(71,0.9999997);
591   weight->SetBinContent(72,0.9999995);
592   weight->SetBinContent(73,0.9999994);
593   weight->SetBinContent(74,1);
594   weight->SetBinContent(75,1);
595   weight->SetBinContent(76,1);
596   weight->SetBinContent(77,1);
597   weight->SetBinContent(78,0.9999999);
598   weight->SetBinContent(79,1);
599   weight->SetBinContent(80,1);
600   weight->SetEntries(526);
601   Int_t limit = weight->GetRandom();
602   pythia->SetTriggerChargedMultiplicity(limit, 1.4);
603   comment = comment.Append(Form("; multiplicity threshold set to %d in |eta| < 1.4", limit));
604   return pythia;
605 }
606
607 //______________________________________________________________________
608 AliGenerator* PythiaJets()
609 {
610       comment = comment.Append(" pp: Pythia low-pt");
611 //
612 //    Pythia
613       AliGenPythia* pythia = new AliGenPythia(-1); 
614       /*   pythia->SetMomentumRange(0, 999999.);
615       pythia->SetThetaRange(0., 180.);
616       pythia->SetYRange(-12.,12.);
617       pythia->SetPtRange(0,1000.);*/
618       pythia->SetProcess(kPyJets);
619       pythia->SetEnergyCMS(energy);
620       pythia->SetStrucFunc(kCTEQ6l);
621       //  pythia->SetPtHard(50, 1000);
622 //
623       return pythia;
624 }
625
626
627
628 //___________________________________________________________________
629 AliGenerator* MbPhojet()
630 {
631   comment = comment.Append(" pp: Pythia low-pt");
632 #if defined(__CINT__)
633   gSystem->Load("libDPMJET");      // Parton density functions
634   gSystem->Load("libTDPMjet");      // Parton density functions
635 #endif
636   AliGenDPMjet* dpmjet = new AliGenDPMjet(-1); 
637
638   /*  dpmjet->SetMomentumRange(0, 999999.);
639   dpmjet->SetThetaRange(0., 180.);
640   dpmjet->SetYRange(-12.,12.);
641   dpmjet->SetPtRange(0,1000.);*/
642
643   dpmjet->SetProcess(kDpmMb);
644   dpmjet->SetEnergyCMS(energy);
645   return dpmjet;
646 }
647
648 //__________________________________________________________________
649 AliGenerator* Hijing()
650 {
651     AliGenHijing *gener = new AliGenHijing(-1);
652     // centre of mass energy 
653     gener->SetEnergyCMS(energy);
654     gener->SetImpactParameterRange(bMin, bMax); 
655     // reference frame
656     gener->SetReferenceFrame("CMS");
657     // projectile
658     gener->SetProjectile("A", 208, 82);
659     gener->SetTarget    ("A", 208, 82);
660     // tell hijing to keep the full parent child chain
661     gener->KeepFullEvent();
662     // enable jet quenching
663     gener->SetJetQuenching(1);
664     // enable shadowing
665     gener->SetDecaysOff(0);
666     gener->SetShadowing(1);
667     // Don't track spectators
668     gener->SetSpectators(0);
669     // kinematic selection
670     gener->SetSelectAll(0);
671     return gener;
672 }
673
674
675 //__________________________________________________________________
676 AliGenerator* Hijing2000()
677 {
678   AliGenHijing *gener = (AliGenHijing*) Hijing();
679   gener->SetJetQuenching(0);    
680   gener->SetPtHardMin (2.3);
681   return gener;
682 }
683
684
685
686 //_____________________________________________________________
687 AliGenerator* Hijing_pA(Bool_t kSlowN) {
688   AliGenCocktail *gener = 0x0;
689   if (kSlowN) {
690     gener  = new AliGenCocktail();
691     gener->SetProjectile("A", 208, 82);
692     gener->SetTarget    ("P",   1,  1);
693     gener->SetEnergyCMS(energy);
694   }
695   
696   AliGenHijing   *hijing = new AliGenHijing(-1);
697   // centre of mass energy 
698   hijing->SetEnergyCMS(energy);
699   // impact parameter range
700   hijing->SetImpactParameterRange(0., 100.);
701   // reference frame
702   hijing->SetReferenceFrame("CMS");
703   hijing->SetBoostLHC(1);
704   // projectile
705   hijing->SetProjectile("A", 208, 82);
706   hijing->SetTarget    ("P", 1, 1);
707   // tell hijing to keep the full parent child chain
708   hijing->KeepFullEvent();
709   // enable jet quenching
710   hijing->SetJetQuenching(0);
711   // enable shadowing
712   hijing->SetShadowing(1);
713   // kinematic selection
714   hijing->SetSelectAll(0);
715   
716   if (!kSlowN) {
717     // DO track spectators
718     hijing->SetSpectators(1);
719     return hijing;
720   }
721   else {
722     // Cocktail with slow nucleons generator
723     // DO NOT track spectators
724     hijing->SetSpectators(0);
725     AliGenSlowNucleons* gray   = new AliGenSlowNucleons(1);
726     AliCollisionGeometry* coll = hijing->CollisionGeometry();
727     AliSlowNucleonModelExp* model = new AliSlowNucleonModelExp();
728     //  Not yet in the release...
729     //    model->SetSaturation(kTRUE);
730     gray->SetSlowNucleonModel(model);
731     gray->SetTarget(208,82);
732     gray->SetThetaDist(1);
733     gray->SetProtonDirection(1);
734     //    gray->SetDebug(1);
735     gray->SetNominalCmsEnergy(2.*pBeamEnergy);
736     gray->NeedsCollisionGeometry();
737     gray->SetCollisionGeometry(coll);
738     
739     gener->AddGenerator(hijing, "Hijing pPb", 1);
740     gener->AddGenerator(gray, "Gray Particles", 1);
741     
742     return gener;
743   }
744 }
745
746
747
748 //__________________________________________________________________
749 AliGenerator* DPMjet()
750 {
751   AliGenDPMjet* dpmjet = new AliGenDPMjet(-1); 
752   dpmjet->SetEnergyCMS(energy);
753   dpmjet->SetProjectile("A", 208, 82);
754   dpmjet->SetTarget    ("A", 208, 82);
755   dpmjet->SetImpactParameterRange(bMin, bMax);
756   dpmjet->SetPi0Decay(0);
757   return dpmjet;
758 }
759
760
761 //____________________________________________
762 AliGenerator* DPMjet_pA(Bool_t fragments)
763 {
764   AliGenDPMjet *gener = new AliGenDPMjet(-1);
765   //  A-p
766   gener->SetProjectile("A", 208, 82);
767   gener->SetTarget("P", 1, 1);
768   //
769   gener->SetEnergyCMS(energy);
770   gener->SetProjectileBeamEnergy(82.*pBeamEnergy/208.);
771
772   //
773   gener->SetProcess(kDpmMb);
774   gener->SetImpactParameterRange(0., 100.);
775   // DO NOT BOOST !
776   //  gener->SetBoostLHC(1);
777   //
778   gener->SetFragmentProd(fragments);
779   
780   return gener;
781
782 }
783
784 //__________________________________________________________________
785 AliGenerator* AmptDefault()
786 {
787   AliGenAmpt *genHi = new AliGenAmpt(-1);
788   //=========================================================================
789   // THE DECAYER
790   AliDecayer *decayer = new AliDecayerPythia();
791   cout << "*****************************************" << endl;
792   genHi->SetForceDecay( kHadronicD );
793   genHi->SetDecayer( decayer );
794   //=========================================================================
795   genHi->SetEnergyCMS(energy);
796   genHi->SetReferenceFrame("CMS");
797   genHi->SetProjectile("A",208,82);
798   genHi->SetTarget("A",208,82);
799     
800   genHi->SetIsoft(1); //1: defaul - 4: string melting
801   genHi->SetStringFrag(0.5,0.9); //Lund string frgmentation parameters
802   genHi->SetPtHardMin (3);
803   //genHi->SetImpactParameterRange(9.,9.5);
804   genHi->SetImpactParameterRange(bMin,bMax);
805     
806   // Xmu = 3.2 fm^-1 and as = 0.33 ==> sigma_{partonic} = 1.5mb
807   // Ntmax = 150
808   // v_2{2} = 0.102105 +/- 0.000266894
809   // v_2{4} = 0.0829477 +/- 0.00106158
810     
811   genHi->SetNtMax(150);        // NTMAX: number of timesteps (D=150)
812   genHi->SetXmu(3.2264);        // parton screening mass in fm^(-1) (D=3.2264d0)
813     
814   genHi->SetJetQuenching(0);  // enable jet quenching
815   genHi->SetShadowing(1);     // enable shadowing
816   genHi->SetDecaysOff(1);     // neutral pion and heavy particle decays switched off
817   genHi->SetSpectators(0);    // track spectators
818   //Boost into LHC lab frame
819   genHi->SetBoostLHC(1);
820   //  genHi->Init();
821   genHi->SetRandomReactionPlane(kTRUE);
822  
823   return genHi;
824 }
825
826 //__________________________________________________________________
827 AliGenerator* AmptStringMelting()
828 {
829   AliGenAmpt *genHi = new AliGenAmpt(-1);
830   //=========================================================================
831   // THE DECAYER
832   AliDecayer *decayer = new AliDecayerPythia();
833   cout << "*****************************************" << endl;
834   genHi->SetForceDecay( kHadronicD );
835   genHi->SetDecayer( decayer );
836   //=========================================================================
837   genHi->SetEnergyCMS(energy);
838   genHi->SetReferenceFrame("CMS");
839   genHi->SetProjectile("A",208,82);
840   genHi->SetTarget("A",208,82);
841     
842   genHi->SetIsoft(4); //1: defaul - 4: string melting
843   genHi->SetStringFrag(0.5,0.9); //Lund string frgmentation parameters
844   genHi->SetPtHardMin (3);
845   //genHi->SetImpactParameterRange(9.,9.5);
846   genHi->SetImpactParameterRange(bMin,bMax);
847
848   // Xmu = 3.2 fm^-1 and as = 0.33 ==> sigma_{partonic} = 1.5mb
849   // Ntmax = 150
850   // v_2{2} = 0.102105 +/- 0.000266894
851   // v_2{4} = 0.0829477 +/- 0.00106158
852
853   genHi->SetNtMax(150);        // NTMAX: number of timesteps (D=150)
854   genHi->SetXmu(3.2264);        // parton screening mass in fm^(-1) (D=3.2264d0)
855
856   genHi->SetJetQuenching(0);  // enable jet quenching
857   genHi->SetShadowing(1);     // enable shadowing
858   genHi->SetDecaysOff(1);     // neutral pion and heavy particle decays switched off
859   genHi->SetSpectators(0);    // track spectators
860   //Boost into LHC lab frame
861   genHi->SetBoostLHC(1);
862 //  genHi->Init();
863   genHi->SetRandomReactionPlane(kTRUE);
864   return genHi;
865
866 }
867
868 //__________________________________________________________________
869 AliGenerator* AmptStringMeltingNoART()
870 {
871     AliGenAmpt *genHi = new AliGenAmpt(-1);
872     //=========================================================================
873     // THE DECAYER
874     AliDecayer *decayer = new AliDecayerPythia();
875     cout << "*****************************************" << endl;
876     genHi->SetForceDecay( kHadronicD );
877     genHi->SetDecayer( decayer );
878     //=========================================================================
879     genHi->SetEnergyCMS(energy);
880     genHi->SetReferenceFrame("CMS");
881     genHi->SetProjectile("A",208,82);
882     genHi->SetTarget("A",208,82);
883     
884     genHi->SetIsoft(4); //1: defaul - 4: string melting
885     genHi->SetStringFrag(0.5,0.9); //Lund string frgmentation parameters
886     genHi->SetPtHardMin (3);
887     //genHi->SetImpactParameterRange(9.,9.5);
888     genHi->SetImpactParameterRange(bMin,bMax);
889     
890     // Xmu = 3.2 fm^-1 and as = 0.33 ==> sigma_{partonic} = 1.5mb
891     // Ntmax = 150
892     // v_2{2} = 0.102105 +/- 0.000266894
893     // v_2{4} = 0.0829477 +/- 0.00106158
894     
895     genHi->SetNtMax(3);        // NTMAX: number of timesteps (D=150)
896     genHi->SetXmu(3.2264);        // parton screening mass in fm^(-1) (D=3.2264d0)
897     
898     genHi->SetJetQuenching(0);  // enable jet quenching
899     genHi->SetShadowing(1);     // enable shadowing
900     genHi->SetDecaysOff(1);     // neutral pion and heavy particle decays switched off
901     genHi->SetSpectators(0);    // track spectators
902     //Boost into LHC lab frame
903     genHi->SetBoostLHC(1);
904     //  genHi->Init();
905     genHi->SetRandomReactionPlane(kTRUE);
906     return genHi;
907     
908 }
909
910
911 //__________________________________________________________________
912 AliGenerator* AmptReducedFlow()
913 {
914   AliGenAmpt *genHi = new AliGenAmpt(-1);
915   //=========================================================================
916   // THE DECAYER
917   AliDecayer *decayer = new AliDecayerPythia();
918   cout << "*****************************************" << endl;
919   genHi->SetForceDecay( kHadronicD );
920   genHi->SetDecayer( decayer );
921   //=========================================================================
922   genHi->SetEnergyCMS(energy);
923   genHi->SetReferenceFrame("CMS");
924   genHi->SetProjectile("A",208,82);
925   genHi->SetTarget("A",208,82);
926     
927   genHi->SetIsoft(4); //1: defaul - 4: string melting
928   genHi->SetStringFrag(0.5,0.9); //Lund string frgmentation parameters
929   genHi->SetPtHardMin (3);
930   //genHi->SetImpactParameterRange(9.,9.5);
931   genHi->SetImpactParameterRange(bMin,bMax);
932
933   // Xmu = 12.4 fm^-1 and as = 0.33 ==> sigma_{partonic} = 0.1mb
934   // Ntmax = 20
935   // flow estimates from Q-cumulants
936   // (POI, without weights)
937   // v_2{2} = 0.0549735 +/- 0.000270249
938   // v_2{4} = 0.0421905 +/- 0.00189449
939
940   genHi->SetNtMax(20);        // NTMAX: number of timesteps (D=150)
941   genHi->SetXmu(12.4);        // parton screening mass in fm^(-1) (D=3.2264d0)
942
943   genHi->SetJetQuenching(0);  // enable jet quenching
944   genHi->SetShadowing(1);     // enable shadowing
945   genHi->SetDecaysOff(1);     // neutral pion and heavy particle decays switched off
946   genHi->SetSpectators(0);    // track spectators
947   //Boost into LHC lab frame
948   genHi->SetBoostLHC(1);
949  // genHi->Init();
950   genHi->SetRandomReactionPlane(kTRUE);
951   return genHi;
952
953 }
954
955 //__________________________________________________________________
956 AliGenerator* AmptpA()
957 {
958     AliGenAmpt *genHi = new AliGenAmpt(-1);
959     //=========================================================================
960     // THE DECAYER
961     AliDecayer *decayer = new AliDecayerPythia();
962     cout << "*****************************************" << endl;
963     genHi->SetForceDecay( kHadronicD );
964     genHi->SetDecayer( decayer );
965     //=========================================================================
966     genHi->SetEnergyCMS(energy);
967     genHi->SetReferenceFrame("CMS");
968     genHi->SetProjectile("A", 208, 82);
969     genHi->SetTarget    ("P", 1, 1);
970     genHi->SetIsoft(4); //1: defaul - 4: string melting
971     genHi->SetStringFrag(0.5,0.9); //Lund string frgmentation parameters
972     genHi->SetPtHardMin (3);
973     //genHi->SetImpactParameterRange(9.,9.5);
974     genHi->SetImpactParameterRange(bMin,bMax);
975     genHi->SetNtMax(1500); //NTMAX: number of timesteps (D=150)
976     genHi->SetXmu(3.2264); //parton screening mass in fm^(-1) (D=3.2264d0)
977     genHi->SetJetQuenching(0); // enable jet quenching
978     genHi->SetShadowing(1);    // enable shadowing
979     genHi->SetDecaysOff(1);    // neutral pion and heavy particle decays switched off
980     genHi->SetSpectators(0);   // track spectators
981     //Boost into LHC lab frame
982     genHi->SetBoostLHC(1);
983     //  genHi->Init();
984     genHi->SetRandomReactionPlane(kTRUE);
985     return genHi;
986     
987 }