Further fixes to the par file generation
[u/mrichter/AliRoot.git] / macros / Config_PythiaHeavyFlavours.C
1 //*******************************************************************//
2 // Configuration file for charm / beauty generation with PYTHIA      //
3 //                                                                   //
4 // The parameters have been tuned in order to reproduce the inclusive//
5 // heavy quark pt distribution given by the NLO pQCD calculation by  //
6 // Mangano, Nason and Ridolfi.                                       //
7 //                                                                   //
8 // For details and for the NORMALIZATION of the yields see:          //
9 //   N.Carrer and A.Dainese,                                         //
10 //   "Charm and beauty production at the LHC",                       //
11 //   ALICE-INT-2003-019, [arXiv:hep-ph/0311225];                     //
12 //   PPR Chapter 6.6, CERN/LHCC 2005-030 (2005).                     //
13 //*******************************************************************//
14 #if !defined(__CINT__) || defined(__MAKECINT__)
15 #include <Riostream.h>
16 #include <TRandom.h>
17 #include <TDatime.h>
18 #include <TSystem.h>
19 #include <TVirtualMC.h>
20 #include <TGeant3TGeo.h>
21 #include "STEER/AliRunLoader.h"
22 #include "STEER/AliRun.h"
23 #include "STEER/AliConfig.h"
24 #include "PYTHIA6/AliDecayerPythia.h"
25 #include "PYTHIA6/AliGenPythia.h"
26 #include "STEER/AliMagFMaps.h"
27 #include "STRUCT/AliBODY.h"
28 #include "STRUCT/AliMAG.h"
29 #include "STRUCT/AliABSOv3.h"
30 #include "STRUCT/AliDIPOv3.h"
31 #include "STRUCT/AliHALLv3.h"
32 #include "STRUCT/AliFRAMEv2.h"
33 #include "STRUCT/AliSHILv3.h"
34 #include "STRUCT/AliPIPEv3.h"
35 #include "ITS/AliITSv11Hybrid.h"
36 #include "TPC/AliTPCv2.h"
37 #include "TOF/AliTOFv6T0.h"
38 #include "HMPID/AliHMPIDv3.h"
39 #include "ZDC/AliZDCv3.h"
40 #include "TRD/AliTRDv1.h"
41 #include "FMD/AliFMDv0.h"
42 #include "MUON/AliMUONv1.h"
43 #include "PHOS/AliPHOSv1.h"
44 #include "PMD/AliPMDv1.h"
45 #include "T0/AliT0v1.h"
46 #include "ACORDE/AliACORDEv1.h"
47 #endif
48
49 //--- Heavy Flavour Production ---
50 enum ProcessHvFl_t 
51 {
52   kCharmPbPb5500,  kCharmpPb8800,  kCharmpp14000,  kCharmpp14000wmi,
53   kD0PbPb5500,     kD0pPb8800,     kD0pp14000,
54   kDPlusPbPb5500,  kDPluspPb8800,  kDPluspp14000,
55   kDPlusStrangePbPb5500, kDPlusStrangepPb8800, kDPlusStrangepp14000,
56   kBeautyPbPb5500, kBeautypPb8800, kBeautypp14000, kBeautypp14000wmi
57 };
58 //--- Decay Mode ---
59 enum DecayHvFl_t 
60 {
61   kNature,  kHadr, kSemiEl, kSemiMu
62 };
63 //--- Rapidity Cut ---
64 enum YCut_t
65 {
66   kFull, kBarrel, kMuonArm
67 };
68 //--- Magnetic Field ---
69 enum Mag_t
70 {
71     k2kG, k4kG, k5kG
72 };
73 //--- Trigger config ---
74 enum TrigConf_t
75 {
76     kDefaultPPTrig, kDefaultPbPbTrig
77 };
78 const char * TrigConfName[] = {
79     "p-p","Pb-Pb"
80 };
81 enum PprGeo_t
82   {
83     kHoles, kNoHoles
84   };
85
86 //--- Functions ---
87 class AliGenPythia;
88 AliGenPythia *PythiaHVQ(ProcessHvFl_t proc);
89 void          LoadPythia();
90
91
92 // This part for configuration
93 static ProcessHvFl_t procHvFl = kCharmPbPb5500;
94 static DecayHvFl_t   decHvFl  = kNature; 
95 static YCut_t        ycut     = kFull;
96 static Mag_t         mag      = k5kG; 
97 static TrigConf_t    trig     = kDefaultPbPbTrig; // default PbPb trigger configuration
98 static PprGeo_t      geo      = kHoles;
99 // nEvts = -1  : you get 1 QQbar pair and all the fragmentation and 
100 //               decay chain
101 // nEvts = N>0 : you get N charm / beauty Hadrons 
102 Int_t nEvts = -1; 
103 // stars = kTRUE : all heavy resonances and their decay stored
104 //       = kFALSE: only final heavy hadrons and their decays stored
105 Bool_t stars = kTRUE;
106
107 // To be used only with kCharmppMNRwmi and kBeautyppMNRwmi
108 // To get a "reasonable" agreement with MNR results, events have to be 
109 // generated with the minimum ptHard set to 2.76 GeV.
110 // To get a "perfect" agreement with MNR results, events have to be 
111 // generated in four ptHard bins with the following relative 
112 // normalizations:
113 //  CHARM
114 // 2.76-3 GeV: 25%
115 //    3-4 GeV: 40%
116 //    4-8 GeV: 29%
117 //     >8 GeV:  6%
118 //  BEAUTY
119 // 2.76-4 GeV:  5% 
120 //    4-6 GeV: 31%
121 //    6-8 GeV: 28%
122 //     >8 GeV: 36%
123 Float_t ptHardMin =  2.76;
124 Float_t ptHardMax = -1.;
125
126
127 // Comment line
128 static TString comment;
129
130 void Config()
131 {
132  
133   //========================//
134   // Set Random Number seed //
135   //========================//
136   TDatime dt;
137   UInt_t curtime=dt.Get();
138   UInt_t procid=gSystem->GetPid();
139   UInt_t seed=curtime-procid;
140
141   //  gRandom->SetSeed(seed);
142   //  cerr<<"Seed for random number generation= "<<seed<<endl; 
143   gRandom->SetSeed(12345);
144   
145   // Load Pythia libraries
146   LoadPythia();
147   // libraries required by geant321
148 #if defined(__CINT__)
149   gSystem->Load("libgeant321");
150 #endif
151
152   new TGeant3TGeo("C++ Interface to Geant3");
153
154   if(!AliCDBManager::Instance()->IsDefaultStorageSet()){
155     AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT");
156     AliCDBManager::Instance()->SetRun(0);
157   }
158   
159   //=======================================================================
160   //  Create the output file
161
162    
163   AliRunLoader* rl=0x0;
164
165   cout<<"Config.C: Creating Run Loader ..."<<endl;
166   rl = AliRunLoader::Open("galice.root",
167                           AliConfig::GetDefaultEventFolderName(),
168                           "recreate");
169   if (rl == 0x0)
170     {
171       gAlice->Fatal("Config.C","Can not instatiate the Run Loader");
172       return;
173     }
174   rl->SetCompressionLevel(2);
175   rl->SetNumberOfEventsPerFile(3);
176   gAlice->SetRunLoader(rl);
177   // gAlice->SetGeometryFromFile("geometry.root");
178   // gAlice->SetGeometryFromCDB();
179
180   // Set the trigger configuration
181   gAlice->SetTriggerDescriptor(TrigConfName[trig]);
182   cout<<"Trigger configuration is set to  "<<TrigConfName[trig]<<endl;
183
184   //
185   //=======================================================================
186   // ************* STEERING parameters FOR ALICE SIMULATION **************
187   // --- Specify event type to be tracked through the ALICE setup
188   // --- All positions are in cm, angles in degrees, and P and E in GeV
189
190
191     gMC->SetProcess("DCAY",1);
192     gMC->SetProcess("PAIR",1);
193     gMC->SetProcess("COMP",1);
194     gMC->SetProcess("PHOT",1);
195     gMC->SetProcess("PFIS",0);
196     gMC->SetProcess("DRAY",0);
197     gMC->SetProcess("ANNI",1);
198     gMC->SetProcess("BREM",1);
199     gMC->SetProcess("MUNU",1);
200     gMC->SetProcess("CKOV",1);
201     gMC->SetProcess("HADR",1);
202     gMC->SetProcess("LOSS",2);
203     gMC->SetProcess("MULS",1);
204     gMC->SetProcess("RAYL",1);
205
206     Float_t cut = 1.e-3;        // 1MeV cut by default
207     Float_t tofmax = 1.e10;
208
209     gMC->SetCut("CUTGAM", cut);
210     gMC->SetCut("CUTELE", cut);
211     gMC->SetCut("CUTNEU", cut);
212     gMC->SetCut("CUTHAD", cut);
213     gMC->SetCut("CUTMUO", cut);
214     gMC->SetCut("BCUTE",  cut); 
215     gMC->SetCut("BCUTM",  cut); 
216     gMC->SetCut("DCUTE",  cut); 
217     gMC->SetCut("DCUTM",  cut); 
218     gMC->SetCut("PPCUTM", cut);
219     gMC->SetCut("TOFMAX", tofmax); 
220
221
222
223
224   // Set External decayer //
225   //======================//
226   TVirtualMCDecayer* decayer = new AliDecayerPythia();
227   // DECAYS
228   //
229   switch(decHvFl) {
230   case kNature:
231     decayer->SetForceDecay(kAll);
232     break;
233   case kHadr:
234     decayer->SetForceDecay(kHadronicD);
235     break;
236   case kSemiEl:
237     decayer->SetForceDecay(kSemiElectronic);
238     break;
239   case kSemiMu:
240     decayer->SetForceDecay(kSemiMuonic);
241     break;
242   }
243   decayer->Init();
244   gMC->SetExternalDecayer(decayer);
245
246   //=========================//
247   // Generator Configuration //
248   //=========================//
249   AliGenPythia *pythia = PythiaHVQ(procHvFl);
250   // FeedDown option
251   pythia->SetFeedDownHigherFamily(kFALSE);
252   // Stack filling option
253   if(!stars) pythia->SetStackFillOpt(AliGenPythia::kParentSelection);
254   // Set Count mode
255   if(nEvts>0) pythia->SetCountMode(AliGenPythia::kCountParents);
256   // DECAYS
257   //
258   switch(decHvFl) {
259   case kNature:
260     pythia->SetForceDecay(kAll);
261     break;
262   case kHadr:
263     pythia->SetForceDecay(kHadronicD);
264     break;
265   case kSemiEl:
266     pythia->SetForceDecay(kSemiElectronic);
267     break;
268   case kSemiMu:
269     pythia->SetForceDecay(kSemiMuonic);
270     break;
271   }
272   // GEOM & KINE CUTS
273   //
274   pythia->SetMomentumRange(0,99999999);
275   pythia->SetPhiRange(0., 360.);
276   pythia->SetThetaRange(0,180);
277   switch(ycut) {
278   case kFull:
279     pythia->SetYRange(-999,999);
280     break;
281   case kBarrel:
282     pythia->SetYRange(-2,2);
283     break;
284   case kMuonArm:
285     pythia->SetYRange(1,6);
286     break;
287   }
288
289   // PRIMARY VERTEX
290   //
291   pythia->SetOrigin(0, 0, 0);    // vertex position
292   pythia->SetSigma(0, 0, 5.3);   // Sigma in (X,Y,Z) (cm) on IP position
293   pythia->SetCutVertexZ(1.);     // Truncate at 1 sigma
294   pythia->SetVertexSmear(kPerEvent); 
295
296
297
298   pythia->SetTrackingFlag(0);
299   // Specify GEANT tracking limits (Rmax, Zmax)
300   //gAlice->GetMCApp()->TrackingLimits(90.,1.0e10);
301
302
303   pythia->Init();
304
305   // FIELD
306   //    
307   if (mag == k2kG) {
308     comment = comment.Append(" | L3 field 0.2 T");
309   } else if (mag == k4kG) {
310     comment = comment.Append(" | L3 field 0.4 T");
311   } else if (mag == k5kG) {
312     comment = comment.Append(" | L3 field 0.5 T");
313   }
314   printf("\n \n Comment: %s \n \n", comment.Data());
315     
316   AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 2, 1., 10., mag);
317   field->SetL3ConstField(0); //Using const. field in the barrel
318   rl->CdGAFile();
319   gAlice->SetField(field);    
320
321
322   // By default all ALICE is switched off
323   Int_t iABSO=0;
324   Int_t iACORDE=0;
325   Int_t iDIPO=0;
326   Int_t iFMD=0;
327   Int_t iFRAME=0;
328   Int_t iHALL=0;
329   Int_t iITS=0;
330   Int_t iMAG=0;
331   Int_t iMUON=0;
332   Int_t iPHOS=0;
333   Int_t iPIPE=0;
334   Int_t iPMD=0;
335   Int_t iHMPID=0;
336   Int_t iSHIL=0;
337   Int_t iT0=0;
338   Int_t iTOF=0;
339   Int_t iTPC=0;
340   Int_t iTRD=0;
341   Int_t iZDC=0;
342
343   //=================== Alice BODY parameters =============================
344   AliBODY *BODY = new AliBODY("BODY","Alice envelop");
345
346   if(iMAG) {
347     //=================== MAG parameters ============================
348     // --- Start with Magnet since detector layouts may be depending ---
349     // --- on the selected Magnet dimensions ---
350     AliMAG *MAG  = new AliMAG("MAG","Magnet");
351   }
352
353
354   if(iABSO) {
355     //=================== ABSO parameters ============================
356     AliABSO *ABSO  = new AliABSOv3("ABSO","Muon Absorber");
357   }
358
359   if(iDIPO) {
360     //=================== DIPO parameters ============================
361
362     AliDIPO *DIPO  = new AliDIPOv3("DIPO","Dipole version 3");
363   }
364
365   if(iHALL) {
366     //=================== HALL parameters ============================
367
368     AliHALL *HALL  = new AliHALLv3("HALL","Alice Hall");
369   }
370
371
372   if(iFRAME) {
373     //=================== FRAME parameters ============================
374
375     AliFRAME *FRAME  = new AliFRAMEv2("FRAME","Space Frame");
376     if (geo == kHoles) {
377       FRAME->SetHoles(1);
378     } else {
379       FRAME->SetHoles(0);
380     }
381
382   }
383
384   if(iSHIL) {
385     //=================== SHIL parameters ============================
386
387     AliSHIL *SHIL  = new AliSHILv3("SHIL","Shielding");
388   }
389
390
391   if(iPIPE) {
392     //=================== PIPE parameters ============================
393
394     AliPIPE *PIPE  = new AliPIPEv3("PIPE","Beam Pipe");
395   }
396
397
398   if(iITS) {
399     //=================== ITS parameters ============================
400
401     AliITS *ITS  = new AliITSv11Hybrid("ITS","ITS v11Hybrid");
402   }
403   
404
405   if(iTPC) {
406     //============================ TPC parameters ===================
407     AliTPC *TPC  = new AliTPCv2("TPC","Default");
408   }
409
410
411   if(iTOF) {
412     //=================== TOF parameters ============================
413     AliTOF *TOF  = new AliTOFv6T0("TOF","normal TOF");
414   }
415
416   if(iHMPID) {
417     //=================== HMPID parameters ===========================
418     AliHMPID *HMPID  = new AliHMPIDv3("HMPID","normal HMPID");    
419
420   }
421
422
423   if(iZDC) {
424     //=================== ZDC parameters ============================
425
426     AliZDC *ZDC  = new AliZDCv3("ZDC","normal ZDC");
427   }
428
429   if(iACORDE) {
430     //=================== ACORDE parameters ============================
431
432     AliACORDE *ACORDE  = new AliACORDEv1("ACORDE","normal ACORDE");
433   }
434
435   if(iTRD) {
436     //=================== TRD parameters ============================
437   
438     AliTRD *TRD  = new AliTRDv1("TRD","TRD slow simulator");
439   }
440
441   if(iFMD) {
442     //=================== FMD parameters ============================
443
444     AliFMD *FMD  = new AliFMDv0("FMD","normal FMD");
445   }
446
447   if(iMUON) {
448     //=================== MUON parameters ===========================
449     AliMUON *MUON  = new AliMUONv1("MUON","default");
450   }
451  
452   //=================== PHOS parameters ===========================
453
454   if(iPHOS) {
455     AliPHOS *PHOS  = new AliPHOSv1("PHOS","GPS2");
456   }
457
458
459   //=================== ACORDE parameters ===========================
460
461   if(iACORDE) {
462     AliACORDE *ACORDE  = new AliACORDEv1("ACORDE","Normal ACORDEGPS2");
463   }
464
465
466   if(iPMD) {
467     //=================== PMD parameters ============================
468
469     AliPMD *PMD  = new AliPMDv1("PMD","normal PMD");
470     PMD->SetPAR(1., 1., 0.8, 0.02);
471     PMD->SetIN(6., 18., -580., 27., 27.);
472     PMD->SetGEO(0.0, 0.2, 4.);
473     PMD->SetPadSize(0.8, 1.0, 1.0, 1.5);
474
475   }
476
477   if(iT0) {
478     //=================== T0 parameters ============================
479     AliT0 *T0  = new AliT0v1("T0","T0 Detector");
480   }
481
482          
483 }
484 //
485 //           PYTHIA
486 //
487 AliGenPythia *PythiaHVQ(ProcessHvFl_t proc) {
488
489   AliGenPythia * gener = 0x0;
490
491   switch(proc) {
492   case kCharmPbPb5500:
493     comment = comment.Append(" Charm in Pb-Pb at 5.5 TeV");
494     gener = new AliGenPythia(nEvts);
495     gener->SetProcess(kPyCharmPbPbMNR);
496     gener->SetStrucFunc(kCTEQ4L);
497     gener->SetPtHard(2.1,-1.0);
498     gener->SetEnergyCMS(5500.);
499     gener->SetNuclei(208,208);
500     break;
501   case kCharmpPb8800:
502     comment = comment.Append(" Charm in p-Pb at 8.8 TeV");
503     gener = new AliGenPythia(nEvts);
504     gener->SetProcess(kPyCharmpPbMNR);
505     gener->SetStrucFunc(kCTEQ4L);
506     gener->SetPtHard(2.1,-1.0);
507     gener->SetEnergyCMS(8800.);
508     gener->SetProjectile("P",1,1);
509     gener->SetTarget("Pb",208,82);
510     break;
511   case kCharmpp14000:
512     comment = comment.Append(" Charm in pp at 14 TeV");
513     gener = new AliGenPythia(nEvts);
514     gener->SetProcess(kPyCharmppMNR);
515     gener->SetStrucFunc(kCTEQ4L);
516     gener->SetPtHard(2.1,-1.0);
517     gener->SetEnergyCMS(14000.);
518     break;
519   case kCharmpp14000wmi:
520     comment = comment.Append(" Charm in pp at 14 TeV with mult. interactions");
521     gener = new AliGenPythia(-1);
522     gener->SetProcess(kPyCharmppMNRwmi);
523     gener->SetStrucFunc(kCTEQ5L);
524     gener->SetPtHard(ptHardMin,ptHardMax);
525     gener->SetEnergyCMS(14000.);
526     break;
527   case kD0PbPb5500:
528     comment = comment.Append(" D0 in Pb-Pb at 5.5 TeV");
529     gener = new AliGenPythia(nEvts);
530     gener->SetProcess(kPyD0PbPbMNR);
531     gener->SetStrucFunc(kCTEQ4L);
532     gener->SetPtHard(2.1,-1.0);
533     gener->SetEnergyCMS(5500.);
534     gener->SetNuclei(208,208);
535     break;
536   case kD0pPb8800:
537     comment = comment.Append(" D0 in p-Pb at 8.8 TeV");
538     gener = new AliGenPythia(nEvts);
539     gener->SetProcess(kPyD0pPbMNR);
540     gener->SetStrucFunc(kCTEQ4L);
541     gener->SetPtHard(2.1,-1.0);
542     gener->SetEnergyCMS(8800.);
543     gener->SetProjectile("P",1,1);
544     gener->SetTarget("Pb",208,82);
545     break;
546   case kD0pp14000:
547     comment = comment.Append(" D0 in pp at 14 TeV");
548     gener = new AliGenPythia(nEvts);
549     gener->SetProcess(kPyD0ppMNR);
550     gener->SetStrucFunc(kCTEQ4L);
551     gener->SetPtHard(2.1,-1.0);
552     gener->SetEnergyCMS(14000.);
553     break;
554   case kDPlusPbPb5500:
555     comment = comment.Append(" DPlus in Pb-Pb at 5.5 TeV");
556     gener = new AliGenPythia(nEvts);
557     gener->SetProcess(kPyDPlusPbPbMNR);
558     gener->SetStrucFunc(kCTEQ4L);
559     gener->SetPtHard(2.1,-1.0);
560     gener->SetEnergyCMS(5500.);
561     gener->SetNuclei(208,208);
562     break;
563   case kDPluspPb8800:
564     comment = comment.Append(" DPlus in p-Pb at 8.8 TeV");
565     gener = new AliGenPythia(nEvts);
566     gener->SetProcess(kPyDPluspPbMNR);
567     gener->SetStrucFunc(kCTEQ4L);
568     gener->SetPtHard(2.1,-1.0);
569     gener->SetEnergyCMS(8800.);
570     gener->SetProjectile("P",1,1);
571     gener->SetTarget("Pb",208,82);
572     break;
573   case kDPluspp14000:
574     comment = comment.Append(" DPlus in pp at 14 TeV");
575     gener = new AliGenPythia(nEvts);
576     gener->SetProcess(kPyDPlusppMNR);
577     gener->SetStrucFunc(kCTEQ4L);
578     gener->SetPtHard(2.1,-1.0);
579     gener->SetEnergyCMS(14000.);
580     break;
581  case kDPlusStrangePbPb5500:
582     comment = comment.Append(" DPlusStrange in Pb-Pb at 5.5 TeV");
583     gener = new AliGenPythia(nEvts);
584     gener->SetProcess(kPyDPlusStrangePbPbMNR);
585     gener->SetStrucFunc(kCTEQ4L);
586     gener->SetPtHard(2.1,-1.0);
587     gener->SetEnergyCMS(5500.);
588     gener->SetNuclei(208,208);
589     break;
590   case kDPlusStrangepPb8800:
591     comment = comment.Append(" DPlusStrange in p-Pb at 8.8 TeV");
592     gener = new AliGenPythia(nEvts);
593     gener->SetProcess(kPyDPlusStrangepPbMNR);
594     gener->SetStrucFunc(kCTEQ4L);
595     gener->SetPtHard(2.1,-1.0);
596     gener->SetEnergyCMS(8800.);
597     gener->SetProjectile("P",1,1);
598     gener->SetTarget("Pb",208,82);
599     break;
600   case kDPlusStrangepp14000:
601     comment = comment.Append(" DPlusStrange in pp at 14 TeV");
602     gener = new AliGenPythia(nEvts);
603     gener->SetProcess(kPyDPlusStrangeppMNR);
604     gener->SetStrucFunc(kCTEQ4L);
605     gener->SetPtHard(2.1,-1.0);
606     gener->SetEnergyCMS(14000.);
607     break;
608   case kBeautyPbPb5500:
609     comment = comment.Append(" Beauty in Pb-Pb at 5.5 TeV");
610     gener = new AliGenPythia(nEvts);
611     gener->SetProcess(kPyBeautyPbPbMNR);
612     gener->SetStrucFunc(kCTEQ4L);
613     gener->SetPtHard(2.75,-1.0);
614     gener->SetEnergyCMS(5500.);
615     gener->SetNuclei(208,208);
616     break;
617   case kBeautypPb8800:
618     comment = comment.Append(" Beauty in p-Pb at 8.8 TeV");
619     gener = new AliGenPythia(nEvts);
620     gener->SetProcess(kPyBeautypPbMNR);
621     gener->SetStrucFunc(kCTEQ4L);
622     gener->SetPtHard(2.75,-1.0);
623     gener->SetEnergyCMS(8800.);
624     gener->SetProjectile("P",1,1);
625     gener->SetTarget("Pb",208,82);
626     break;
627   case kBeautypp14000:
628     comment = comment.Append(" Beauty in pp at 14 TeV");
629     gener = new AliGenPythia(nEvts);
630     gener->SetProcess(kPyBeautyppMNR);
631     gener->SetStrucFunc(kCTEQ4L);
632     gener->SetPtHard(2.75,-1.0);
633     gener->SetEnergyCMS(14000.);
634     break;
635   case kBeautypp14000wmi:
636     comment = comment.Append(" Beauty in pp at 14 TeV with mult. interactions");
637     gener = new AliGenPythia(-1);
638     gener->SetProcess(kPyBeautyppMNRwmi);
639     gener->SetStrucFunc(kCTEQ5L);
640     gener->SetPtHard(ptHardMin,ptHardMax);
641     gener->SetEnergyCMS(14000.);
642     break;
643   }
644
645   return gener;
646 }
647
648
649
650 void LoadPythia()
651 {
652     // Load Pythia related libraries
653     gSystem->Load("liblhapdf.so");      // Parton density functions
654     gSystem->Load("libEGPythia6.so");   // TGenerator interface
655     gSystem->Load("libpythia6.so");     // Pythia
656     gSystem->Load("libAliPythia6.so");  // ALICE specific implementations
657 }
658
659
660