PDC'06 configurations for charm and beauty added.
[u/mrichter/AliRoot.git] / PYTHIA6 / AliPythia.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 #include "AliPythia.h"
19 #include "AliPythiaRndm.h"
20 #include "../FASTSIM/AliFastGlauber.h"
21 #include "../FASTSIM/AliQuenchingWeights.h"
22 #include "TVector3.h"
23
24 ClassImp(AliPythia)
25
26 #ifndef WIN32
27 # define pyclus pyclus_
28 # define pycell pycell_
29 # define pyshow pyshow_
30 # define pyrobo pyrobo_
31 # define pyquen pyquen_
32 # define pyevnw pyevnw_
33 # define type_of_call
34 #else
35 # define pyclus PYCLUS
36 # define pycell PYCELL
37 # define pyrobo PYROBO
38 # define pyquen PYQUEN
39 # define pyevnw PYEVNW
40 # define type_of_call _stdcall
41 #endif
42
43 extern "C" void type_of_call pyclus(Int_t & );
44 extern "C" void type_of_call pycell(Int_t & );
45 extern "C" void type_of_call pyshow(Int_t &, Int_t &, Double_t &);
46 extern "C" void type_of_call pyrobo(Int_t &, Int_t &, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &);
47 extern "C" void type_of_call pyquen(Double_t &, Int_t &, Double_t &);
48 extern "C" void type_of_call pyevnw(){;}
49
50 //_____________________________________________________________________________
51
52 AliPythia* AliPythia::fgAliPythia=NULL;
53
54 AliPythia::AliPythia()
55 {
56 // Default Constructor
57 //
58 //  Set random number
59     if (!AliPythiaRndm::GetPythiaRandom()) 
60       AliPythiaRndm::SetPythiaRandom(GetRandom());
61     fGlauber          = 0;
62     fQuenchingWeights = 0;
63 }
64
65 void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfunc)
66 {
67 // Initialise the process to generate 
68     if (!AliPythiaRndm::GetPythiaRandom()) 
69       AliPythiaRndm::SetPythiaRandom(GetRandom());
70     
71     fProcess = process;
72     fEcms = energy;
73     fStrucFunc = strucfunc;
74 //...Switch off decay of pi0, K0S, Lambda, Sigma+-, Xi0-, Omega-.
75     SetMDCY(Pycomp(111) ,1,0);
76     SetMDCY(Pycomp(310) ,1,0);
77     SetMDCY(Pycomp(3122),1,0);
78     SetMDCY(Pycomp(3112),1,0);
79     SetMDCY(Pycomp(3212),1,0);
80     SetMDCY(Pycomp(3222),1,0);
81     SetMDCY(Pycomp(3312),1,0);
82     SetMDCY(Pycomp(3322),1,0);
83     SetMDCY(Pycomp(3334),1,0);
84     //  select structure function 
85     SetMSTP(52,2);
86     SetMSTP(51,strucfunc);
87 //
88 // Pythia initialisation for selected processes//
89 //
90 // Make MSEL clean
91 //
92     for (Int_t i=1; i<= 200; i++) {
93         SetMSUB(i,0);
94     }
95 //  select charm production
96     switch (process) 
97     {
98     case kPyOldUEQ2ordered:  //Old underlying events with Q2 ordered QCD processes
99 //        Multiple interactions on.
100         SetMSTP(81,1);
101 // Double Gaussian matter distribution.
102         SetMSTP(82,4);
103         SetPARP(83,0.5);
104         SetPARP(84,0.4);
105 //  pT0.
106         SetPARP(82,2.0);
107 //  Reference energy for pT0 and energy rescaling pace.
108         SetPARP(89,1800);
109         SetPARP(90,0.25);
110 //  String drawing almost completely minimizes string length.
111         SetPARP(85,0.9);
112         SetPARP(86,0.95);
113 // ISR and FSR activity.
114         SetPARP(67,4);
115         SetPARP(71,4);
116 // Lambda_FSR scale.
117         SetPARJ(81,0.29);
118         break;
119     case kPyOldUEQ2ordered2:   
120 // Old underlying events with Q2 ordered QCD processes
121 // Multiple interactions on.
122         SetMSTP(81,1);
123 // Double Gaussian matter distribution.
124         SetMSTP(82,4);
125         SetPARP(83,0.5);
126         SetPARP(84,0.4);
127 // pT0.
128         SetPARP(82,2.0);
129 // Reference energy for pT0 and energy rescaling pace.
130         SetPARP(89,1800);
131         SetPARP(90,0.16);  // here is the difference with  kPyOldUEQ2ordered
132 // String drawing almost completely minimizes string length.
133         SetPARP(85,0.9);
134         SetPARP(86,0.95);
135 // ISR and FSR activity.
136         SetPARP(67,4);
137         SetPARP(71,4);
138 // Lambda_FSR scale.
139         SetPARJ(81,0.29);       
140         break;
141     case kPyOldPopcorn:  
142 // Old production mechanism: Old Popcorn
143         SetMSEL(1);
144         SetMSTJ(12,3); 
145 // (D=2) Like MSTJ(12)=2 but added prod ofthe 1er rank baryon
146         SetMSTP(88,2); 
147 // (D=1)see can be used to form  baryons (BARYON JUNCTION)
148         SetMSTJ(1,1);  
149         AtlasTuning();
150         break;
151     case kPyCharm:
152         SetMSEL(4);
153 //  heavy quark masses
154
155         SetPMAS(4,1,1.2);
156         SetMSTU(16,2);
157 //
158 //    primordial pT
159         SetMSTP(91,1);
160         SetPARP(91,1.);
161         SetPARP(93,5.);
162 //
163         break;
164     case kPyBeauty:
165         SetMSEL(5);
166         SetPMAS(5,1,4.75);
167         SetMSTU(16,2);
168         break;
169     case kPyJpsi:
170         SetMSEL(0);
171 // gg->J/Psi g
172         SetMSUB(86,1);
173         break;
174     case kPyJpsiChi:
175         SetMSEL(0);
176 // gg->J/Psi g
177         SetMSUB(86,1);
178 // gg-> chi_0c g
179         SetMSUB(87,1);
180 // gg-> chi_1c g
181         SetMSUB(88,1);
182 // gg-> chi_2c g
183         SetMSUB(89,1);  
184         break;
185     case kPyCharmUnforced:
186         SetMSEL(0);
187 // gq->qg   
188         SetMSUB(28,1);
189 // gg->qq
190         SetMSUB(53,1);
191 // gg->gg
192         SetMSUB(68,1);
193         break;
194     case kPyBeautyUnforced:
195         SetMSEL(0);
196 // gq->qg   
197         SetMSUB(28,1);
198 // gg->qq
199         SetMSUB(53,1);
200 // gg->gg
201         SetMSUB(68,1);
202         break;
203     case kPyMb:
204 // Minimum Bias pp-Collisions
205 //
206 //   
207 //      select Pythia min. bias model
208         SetMSEL(0);
209         SetMSUB(92,1);             // single diffraction AB-->XB
210         SetMSUB(93,1);             // single diffraction AB-->AX
211         SetMSUB(94,1);             // double diffraction
212         SetMSUB(95,1);             // low pt production
213
214         AtlasTuning();
215         break;
216     case kPyMbNonDiffr:
217 // Minimum Bias pp-Collisions
218 //
219 //   
220 //      select Pythia min. bias model
221         SetMSEL(0);
222         SetMSUB(95,1);             // low pt production
223
224         AtlasTuning();
225         break;
226     case kPyJets:
227 //
228 //  QCD Jets
229 //
230         SetMSEL(1);
231  // Pythia Tune A (CDF)
232  //
233        SetPARP(67,4.);            // Regulates Initial State Radiation
234        SetMSTP(82,4);             // Double Gaussian Model
235        SetPARP(82,2.0);           // [GeV]    PT_min at Ref. energy
236        SetPARP(84,0.4);           // Core radius
237        SetPARP(85,0.90) ;         // Regulates gluon prod. mechanism
238        SetPARP(86,0.95);          // Regulates gluon prod. mechanism
239        SetPARP(89,1800.);         // [GeV]   Ref. energy
240        SetPARP(90,0.25);          // 2*epsilon (exponent in power law)
241        break;
242     case kPyDirectGamma:
243         SetMSEL(10);
244         break;
245     case kPyCharmPbPbMNR:
246     case kPyD0PbPbMNR:
247     case kPyDPlusPbPbMNR:
248     case kPyDPlusStrangePbPbMNR:
249       // Tuning of Pythia parameters aimed to get a resonable agreement
250       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
251       // c-cbar single inclusive and double differential distributions.
252       // This parameter settings are meant to work with Pb-Pb collisions
253       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
254       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
255       // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
256         ConfigHeavyFlavor();
257       // Intrinsic <kT>
258       SetMSTP(91,1);
259       SetPARP(91,1.304);
260       SetPARP(93,6.52);
261       // Set c-quark mass
262       SetPMAS(4,1,1.2);
263       break;
264     case kPyCharmpPbMNR:
265     case kPyD0pPbMNR:
266     case kPyDPluspPbMNR:
267     case kPyDPlusStrangepPbMNR:
268       // Tuning of Pythia parameters aimed to get a resonable agreement
269       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
270       // c-cbar single inclusive and double differential distributions.
271       // This parameter settings are meant to work with p-Pb collisions
272       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
273       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
274       // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
275         ConfigHeavyFlavor();
276       // Intrinsic <kT>
277         SetMSTP(91,1);
278         SetPARP(91,1.16);
279         SetPARP(93,5.8);
280         
281       // Set c-quark mass
282         SetPMAS(4,1,1.2);
283       break;
284     case kPyCharmppMNR:
285     case kPyD0ppMNR:
286     case kPyDPlusppMNR:
287     case kPyDPlusStrangeppMNR:
288       // Tuning of Pythia parameters aimed to get a resonable agreement
289       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
290       // c-cbar single inclusive and double differential distributions.
291       // This parameter settings are meant to work with pp collisions
292       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
293       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
294       // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
295         ConfigHeavyFlavor();
296       // Intrinsic <kT^2>
297         SetMSTP(91,1);
298         SetPARP(91,1.);
299         SetPARP(93,5.);
300         
301       // Set c-quark mass
302         SetPMAS(4,1,1.2);
303       break;
304     case kPyCharmppMNRwmi:
305       // Tuning of Pythia parameters aimed to get a resonable agreement
306       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
307       // c-cbar single inclusive and double differential distributions.
308       // This parameter settings are meant to work with pp collisions
309       // and with kCTEQ5L PDFs.
310       // Added multiple interactions according to ATLAS tune settings.
311       // To get a "reasonable" agreement with MNR results, events have to be 
312       // generated with the minimum ptHard (AliGenPythia::SetPtHard)
313       // set to 2.76 GeV.
314       // To get a "perfect" agreement with MNR results, events have to be 
315       // generated in four ptHard bins with the following relative 
316       // normalizations:
317       // 2.76-3 GeV: 25%
318       //    3-4 GeV: 40%
319       //    4-8 GeV: 29%
320       //     >8 GeV:  6%
321         ConfigHeavyFlavor();
322       // Intrinsic <kT^2>
323         SetMSTP(91,1);
324         SetPARP(91,1.);
325         SetPARP(93,5.);
326
327       // Set c-quark mass
328         SetPMAS(4,1,1.2);
329         AtlasTuning();
330         break;
331     case kPyBeautyPbPbMNR:
332       // Tuning of Pythia parameters aimed to get a resonable agreement
333       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
334       // b-bbar single inclusive and double differential distributions.
335       // This parameter settings are meant to work with Pb-Pb collisions
336       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
337       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
338       // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
339         ConfigHeavyFlavor();
340       // QCD scales
341         SetPARP(67,1.0);
342         SetPARP(71,1.0);
343       // Intrinsic <kT>
344         SetMSTP(91,1);
345         SetPARP(91,2.035);
346         SetPARP(93,10.17);
347       // Set b-quark mass
348         SetPMAS(5,1,4.75);
349       break;
350     case kPyBeautypPbMNR:
351       // Tuning of Pythia parameters aimed to get a resonable agreement
352       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
353       // b-bbar single inclusive and double differential distributions.
354       // This parameter settings are meant to work with p-Pb collisions
355       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
356       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
357       // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
358         ConfigHeavyFlavor();
359       // QCD scales
360         SetPARP(67,1.0);
361         SetPARP(71,1.0);
362       // Intrinsic <kT>
363         SetMSTP(91,1);
364         SetPARP(91,1.60);
365         SetPARP(93,8.00);
366       // Set b-quark mass
367         SetPMAS(5,1,4.75);
368       break;
369     case kPyBeautyppMNR:
370       // Tuning of Pythia parameters aimed to get a resonable agreement
371       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
372       // b-bbar single inclusive and double differential distributions.
373       // This parameter settings are meant to work with pp collisions
374       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
375       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
376       // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
377         ConfigHeavyFlavor();
378       // QCD scales
379         SetPARP(67,1.0);
380         SetPARP(71,1.0);
381         
382         // Intrinsic <kT>
383         SetMSTP(91,1);
384         SetPARP(91,1.);
385         SetPARP(93,5.);
386         
387         // Set b-quark mass
388         SetPMAS(5,1,4.75);
389       break;
390      case kPyBeautyppMNRwmi:
391       // Tuning of Pythia parameters aimed to get a resonable agreement
392       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
393       // b-bbar single inclusive and double differential distributions.
394       // This parameter settings are meant to work with pp collisions
395       // and with kCTEQ5L PDFs.
396       // Added multiple interactions according to ATLAS tune settings.
397       // To get a "reasonable" agreement with MNR results, events have to be 
398       // generated with the minimum ptHard (AliGenPythia::SetPtHard)
399       // set to 2.76 GeV.
400       // To get a "perfect" agreement with MNR results, events have to be 
401       // generated in four ptHard bins with the following relative 
402       // normalizations:
403       // 2.76-4 GeV:  5% 
404       //    4-6 GeV: 31%
405       //    6-8 GeV: 28%
406       //     >8 GeV: 36%
407          ConfigHeavyFlavor();
408       // QCD scales
409          SetPARP(67,1.0);
410          SetPARP(71,1.0);
411          
412          // Intrinsic <kT>
413          SetMSTP(91,1);
414          SetPARP(91,1.);
415          SetPARP(93,5.);
416
417       // Set b-quark mass
418          SetPMAS(5,1,4.75);
419
420          AtlasTuning();
421          break; 
422     case kPyW:
423
424       //Inclusive production of W+/-
425       SetMSEL(0);
426       //f fbar -> W+ 
427       SetMSUB(2,1);
428       //        //f fbar -> g W+
429       //        SetMSUB(16,1);
430       //        //f fbar -> gamma W+
431       //        SetMSUB(20,1);
432       //        //f g -> f W+  
433       //        SetMSUB(31,1);
434       //        //f gamma -> f W+
435       //        SetMSUB(36,1);
436       
437       // Initial/final parton shower on (Pythia default)
438       // With parton showers on we are generating "W inclusive process"
439       SetMSTP(61,1); //Initial QCD & QED showers on
440       SetMSTP(71,1); //Final QCD & QED showers on
441       
442       break;  
443
444     case kPyZ:
445
446       //Inclusive production of Z
447       SetMSEL(0);
448       //f fbar -> Z/gamma
449       SetMSUB(1,1);
450       
451       //       // f fbar -> g Z/gamma
452       //       SetMSUB(15,1);
453       //       // f fbar -> gamma Z/gamma
454       //       SetMSUB(19,1);
455       //       // f g -> f Z/gamma
456       //       SetMSUB(30,1);
457       //       // f gamma -> f Z/gamma
458       //       SetMSUB(35,1);
459       
460       //only Z included, not gamma
461       SetMSTP(43,2);
462       
463       // Initial/final parton shower on (Pythia default)
464       // With parton showers on we are generating "Z inclusive process"
465       SetMSTP(61,1); //Initial QCD & QED showers on
466       SetMSTP(71,1); //Final QCD & QED showers on
467       
468       break;  
469
470     }
471 //
472 //  Initialize PYTHIA
473     SetMSTP(41,1);   // all resonance decays switched on
474
475     Initialize("CMS","p","p",fEcms);
476
477 }
478
479 Int_t AliPythia::CheckedLuComp(Int_t kf)
480 {
481 // Check Lund particle code (for debugging)
482     Int_t kc=Pycomp(kf);
483     printf("\n Lucomp kf,kc %d %d",kf,kc);
484     return kc;
485 }
486
487 void AliPythia::SetNuclei(Int_t a1, Int_t a2)
488 {
489 // Treat protons as inside nuclei with mass numbers a1 and a2  
490 //    The MSTP array in the PYPARS common block is used to enable and 
491 //    select the nuclear structure functions. 
492 //    MSTP(52)  : (D=1) choice of proton and nuclear structure-function library
493 //            =1: internal PYTHIA acording to MSTP(51) 
494 //            =2: PDFLIB proton  s.f., with MSTP(51)  = 1000xNGROUP+NSET
495 //    If the following mass number both not equal zero, nuclear corrections of the stf are used.
496 //    MSTP(192) : Mass number of nucleus side 1
497 //    MSTP(193) : Mass number of nucleus side 2
498     SetMSTP(52,2);
499     SetMSTP(192, a1);
500     SetMSTP(193, a2);  
501 }
502         
503
504 AliPythia* AliPythia::Instance()
505
506 // Set random number generator 
507     if (fgAliPythia) {
508         return fgAliPythia;
509     } else {
510         fgAliPythia = new AliPythia();
511         return fgAliPythia;
512     }
513 }
514
515 void AliPythia::PrintParticles()
516
517 // Print list of particl properties
518     Int_t np = 0;
519     char*   name = new char[16];    
520     for (Int_t kf=0; kf<1000000; kf++) {
521         for (Int_t c = 1;  c > -2; c-=2) {
522             Int_t kc = Pycomp(c*kf);
523             if (kc) {
524                 Float_t mass  = GetPMAS(kc,1);
525                 Float_t width = GetPMAS(kc,2);  
526                 Float_t tau   = GetPMAS(kc,4);
527
528                 Pyname(kf,name);
529         
530                 np++;
531                 
532                 printf("\n mass, width, tau: %6d %s %10.3f %10.3e %10.3e", 
533                        c*kf, name, mass, width, tau);
534             }
535         }
536     }
537     printf("\n Number of particles %d \n \n", np);
538 }
539
540 void  AliPythia::ResetDecayTable()
541 {
542 //  Set default values for pythia decay switches
543     Int_t i;
544     for (i = 1; i <  501; i++) SetMDCY(i,1,fDefMDCY[i]);
545     for (i = 1; i < 2001; i++) SetMDME(i,1,fDefMDME[i]);
546 }
547
548 void  AliPythia::SetDecayTable()
549 {
550 //  Set default values for pythia decay switches
551 //
552     Int_t i;
553     for (i = 1; i <  501; i++) fDefMDCY[i] = GetMDCY(i,1);
554     for (i = 1; i < 2001; i++) fDefMDME[i] = GetMDME(i,1);
555 }
556
557 void  AliPythia::Pyclus(Int_t& njet)
558 {
559 //  Call Pythia clustering algorithm
560 //
561     pyclus(njet);
562 }
563
564 void  AliPythia::Pycell(Int_t& njet)
565 {
566 //  Call Pythia jet reconstruction algorithm
567 //
568     pycell(njet);
569 }
570
571 void  AliPythia::Pyshow(Int_t ip1, Int_t ip2, Double_t qmax)
572 {
573 //  Call Pythia jet reconstruction algorithm
574 //
575     pyshow(ip1, ip2, qmax);
576 }
577
578 void AliPythia::Pyrobo(Int_t imi, Int_t ima, Double_t the, Double_t phi, Double_t bex, Double_t bey, Double_t bez)
579 {
580     pyrobo(imi, ima, the, phi, bex, bey, bez);
581 }
582
583
584
585 void AliPythia::InitQuenching(Float_t cMin, Float_t cMax, Float_t k, Int_t iECMethod)
586 {
587 // Initializes 
588 // (1) The quenching model using quenching weights according to C. Salgado and U. Wiedemann
589 // (2) The nuclear geometry using the Glauber Model
590 //     
591
592
593     fGlauber = new AliFastGlauber();
594     fGlauber->Init(2);
595     fGlauber->SetCentralityClass(cMin, cMax); 
596
597     fQuenchingWeights = new AliQuenchingWeights();
598     fQuenchingWeights->InitMult();
599     fQuenchingWeights->SetK(k);
600     fQuenchingWeights->SetECMethod(AliQuenchingWeights::kECMethod(iECMethod));
601 }
602
603
604 void  AliPythia::Quench()
605 {
606 //
607 //
608 //  Simple Jet Quenching routine:
609 //  =============================
610 //  The jet formed by all final state partons radiated by the parton created 
611 //  in the hard collisions is quenched by a factor (1-z) using light cone variables in 
612 //  the initial parton reference frame:
613 //  (E + p_z)new = (1-z) (E + p_z)old
614 //
615 //
616 //
617 //
618 //   The lost momentum is first balanced by one gluon with virtuality > 0.   
619 //   Subsequently the gluon splits to yield two gluons with E = p.
620 //
621 //
622 // 
623     static Float_t eMean = 0.;
624     static Int_t   icall = 0;
625     
626     Double_t p0[4][5];
627     Double_t p1[4][5];
628     Double_t p2[4][5];
629     Int_t   klast[4] = {-1, -1, -1, -1};
630
631     Int_t numpart   = fPyjets->N;
632     Double_t px = 0., py = 0., pz = 0., e = 0., m = 0., p = 0., pt = 0., theta = 0., phi = 0.;
633     Double_t pxq[4], pyq[4], pzq[4], eq[4], yq[4], mq[4], pq[4], phiq[4], thetaq[4], ptq[4];
634     Bool_t  quenched[4];
635     Double_t wjtKick[4];
636     Int_t nGluon[4];
637     Int_t qPdg[4];
638     Int_t   imo, kst, pdg;
639     
640 //
641 //  Sore information about Primary partons
642 //
643 //  j =
644 //  0, 1 partons from hard scattering
645 //  2, 3 partons from initial state radiation
646 // 
647     for (Int_t i = 2; i <= 7; i++) {
648         Int_t j = 0;
649         // Skip gluons that participate in hard scattering
650         if (i == 4 || i == 5) continue;
651         // Gluons from hard Scattering
652         if (i == 6 || i == 7) {
653             j = i - 6;
654             pxq[j]    = fPyjets->P[0][i];
655             pyq[j]    = fPyjets->P[1][i];
656             pzq[j]    = fPyjets->P[2][i];
657             eq[j]     = fPyjets->P[3][i];
658             mq[j]     = fPyjets->P[4][i];
659         } else {
660             // Gluons from initial state radiation
661             //
662             // Obtain 4-momentum vector from difference between original parton and parton after gluon 
663             // radiation. Energy is calculated independently because initial state radition does not 
664             // conserve strictly momentum and energy for each partonic system independently.
665             //
666             // Not very clean. Should be improved !
667             //
668             //
669             j = i;
670             pxq[j]    = fPyjets->P[0][i] - fPyjets->P[0][i+2];
671             pyq[j]    = fPyjets->P[1][i] - fPyjets->P[1][i+2];
672             pzq[j]    = fPyjets->P[2][i] - fPyjets->P[2][i+2];
673             mq[j]     = fPyjets->P[4][i];
674             eq[j]     = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j] + pzq[j] * pzq[j] + mq[j] * mq[j]);
675         }
676 //
677 //  Calculate some kinematic variables
678 //
679         yq[j]     = 0.5 * TMath::Log((eq[j] + pzq[j] + 1.e-14) / (eq[j] - pzq[j] + 1.e-14));
680         pq[j]     = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j] + pzq[j] * pzq[j]);
681         phiq[j]   = TMath::Pi()+TMath::ATan2(-pyq[j], -pxq[j]);
682         ptq[j]    = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j]);
683         thetaq[j] = TMath::ATan2(ptq[j], pzq[j]);
684         qPdg[j]   =  fPyjets->K[1][i];
685     }
686   
687     Double_t int0[4];
688     Double_t int1[4];
689     
690     fGlauber->GetI0I1ForPythiaAndXY(4, phiq, int0, int1, fXJet, fYJet, 15.);
691
692     for (Int_t j = 0; j < 4; j++) {
693         //
694         // Quench only central jets and with E > 10.
695         //
696
697
698         Int_t itype = (qPdg[j] == 21) ? 2 : 1;
699         Double_t eloss = fQuenchingWeights->GetELossRandomKFast(itype, int0[j], int1[j], eq[j]);
700
701         if (TMath::Abs(yq[j]) > 2.5 || eq[j] < 10.) {
702             fZQuench[j] = 0.;
703         } else {
704             if (eq[j] > 40. && TMath::Abs(yq[j]) < 0.5) {
705                 icall ++;
706                 eMean += eloss;
707             }
708             //
709             // Extra pt
710             Double_t l =   fQuenchingWeights->CalcLk(int0[j], int1[j]);     
711             wjtKick[j] = TMath::Sqrt(l *  fQuenchingWeights->CalcQk(int0[j], int1[j]));
712             //
713             // Fractional energy loss
714             fZQuench[j] = eloss / eq[j];
715             //
716             // Avoid complete loss
717             //
718             if (fZQuench[j] == 1.) fZQuench[j] = 0.95;
719             //
720             // Some debug printing
721
722             
723 //          printf("Initial parton # %3d, Type %3d Energy %10.3f Phi %10.3f Length %10.3f Loss %10.3f Kick %10.3f Mean: %10.3f %10.3f\n", 
724 //                 j, itype, eq[j], phiq[j], l, eloss, wjtKick[j], eMean / Float_t(icall+1), yq[j]);
725             
726 //          fZQuench[j] = 0.8;
727 //          while (fZQuench[j] >= 0.95)  fZQuench[j] = gRandom->Exp(0.2);
728         }
729         
730         quenched[j] = (fZQuench[j] > 0.01);
731     } // primary partons
732     
733     
734
735     Double_t pNew[1000][4];
736     Int_t    kNew[1000];
737     Int_t icount = 0;
738     Double_t zquench[4];
739     
740 //
741 //  System Loop    
742     for (Int_t isys = 0; isys < 4; isys++) {
743 //      Skip to next system if not quenched.
744         if (!quenched[isys]) continue;
745         
746         nGluon[isys]   = 1 + Int_t(fZQuench[isys] / (1. - fZQuench[isys]));
747         if (nGluon[isys] > 6) nGluon[isys] = 6;
748         zquench[isys] = 1. - TMath::Power(1. - fZQuench[isys], 1./Double_t(nGluon[isys]));
749         wjtKick[isys]  = wjtKick[isys] / TMath::Sqrt(Double_t(nGluon[isys]));
750
751
752         
753         Int_t igMin = -1;
754         Int_t igMax = -1;
755         Double_t pg[4] = {0., 0., 0., 0.};
756         
757 //
758 // Loop on radiation events
759
760         for (Int_t iglu = 0; iglu < nGluon[isys]; iglu++) {
761             while (1) {
762                 icount = 0;
763                 for (Int_t k = 0; k < 4; k++)
764                 {
765                     p0[isys][k] = 0.;
766                     p1[isys][k] = 0.;
767                     p2[isys][k] = 0.;
768                 }
769 //      Loop over partons
770                 for (Int_t i = 0; i < numpart; i++)
771                 {
772                     imo =  fPyjets->K[2][i];
773                     kst =  fPyjets->K[0][i];
774                     pdg =  fPyjets->K[1][i];
775                     
776                 
777                 
778 //      Quarks and gluons only
779                     if (pdg != 21 && TMath::Abs(pdg) > 6) continue;
780 //      Particles from hard scattering only
781                     
782                     if (imo > 8 && imo < 1000) imo = fPyjets->K[2][imo - 1];
783                     Int_t imom = imo % 1000;
784                     if ((isys == 0 || isys == 1) && ((imom != (isys + 7)))) continue;
785                     if ((isys == 2 || isys == 3) && ((imom != (isys + 1)))) continue;               
786                     
787                     
788 //      Skip comment lines
789                     if (kst != 1 && kst != 2) continue;
790 //
791 //      Parton kinematic
792                     px    = fPyjets->P[0][i];
793                     py    = fPyjets->P[1][i];
794                     pz    = fPyjets->P[2][i];
795                     e     = fPyjets->P[3][i];
796                     m     = fPyjets->P[4][i];
797                     pt    = TMath::Sqrt(px * px + py * py);
798                     p     = TMath::Sqrt(px * px + py * py + pz * pz); 
799                     phi   = TMath::Pi() + TMath::ATan2(-py, -px);
800                     theta = TMath::ATan2(pt, pz);
801                 
802 //
803 //      Save 4-momentum sum for balancing
804                     Int_t index = isys;
805                     
806                     p0[index][0] += px;
807                     p0[index][1] += py;
808                     p0[index][2] += pz;
809                     p0[index][3] += e;
810                 
811                     klast[index] = i;
812                     
813 //
814 //      Fractional energy loss
815                     Double_t z = zquench[index];
816                     
817                     
818 //      Don't fully quench radiated gluons
819 //
820                     if (imo > 1000) {
821 //      This small factor makes sure that the gluons are not too close in phase space to avoid recombination
822 //
823
824                         z = 0.02;
825                     }
826 //                  printf("z: %d %f\n", imo, z);
827                     
828
829 //
830                     
831                     //
832                     //
833                     //      Transform into frame in which initial parton is along z-axis
834                     //
835                     TVector3 v(px, py, pz);
836                     v.RotateZ(-phiq[index]);  v.RotateY(-thetaq[index]);
837                     Double_t pxs = v.X(); Double_t pys = v.Y(); Double_t pl  = v.Z();
838
839                     Double_t jt  = TMath::Sqrt(pxs * pxs + pys * pys);
840                     Double_t mt2 = jt * jt + m * m;
841                     Double_t zmax = 1.;     
842                     //
843                     // Kinematic limit on z
844                     //
845                     if (m > 0.) zmax = 1. - m / TMath::Sqrt(m * m + jt * jt);
846                     //
847                     // Change light-cone kinematics rel. to initial parton
848                     //  
849                     Double_t eppzOld = e + pl;
850                     Double_t empzOld = e - pl;
851                     
852                     Double_t eppzNew = (1. - z) * eppzOld;
853                     Double_t empzNew = empzOld - mt2 * z / eppzOld;
854                     Double_t eNew    = 0.5 * (eppzNew + empzNew);
855                     Double_t plNew   = 0.5 * (eppzNew - empzNew);
856                     
857                     Double_t jtNew;
858                     //
859                     // if mt very small (or sometimes even < 0 for numerical reasons) set it to 0
860                     Double_t mt2New = eppzNew * empzNew;
861                     if (mt2New < 1.e-8) mt2New = 0.;
862                     if (z < zmax) {
863                         if (m * m > mt2New) {
864                             //
865                             // This should not happen 
866                             //
867                             Fatal("Quench()", "This should never happen %e %e %e!", m, eppzNew, empzNew);
868                             jtNew = 0;
869                         } else {
870                             jtNew    = TMath::Sqrt(mt2New - m * m);
871                         }
872                     } else {
873                         // If pT is to small (probably a leading massive particle) we scale only the energy
874                         // This can cause negative masses of the radiated gluon
875                         // Let's hope for the best ...
876                         jtNew = jt;
877                         eNew  = TMath::Sqrt(plNew * plNew + mt2);
878                         
879                     }
880                     //
881                     //     Calculate new px, py
882                     //
883                     Double_t pxNew   = jtNew / jt * pxs;
884                     Double_t pyNew   = jtNew / jt * pys;        
885                     
886 //                  Double_t dpx = pxs - pxNew;
887 //                  Double_t dpy = pys - pyNew;
888 //                  Double_t dpz = pl  - plNew;
889 //                  Double_t de  = e   - eNew;
890 //                  Double_t dmass2 = de * de  - dpx * dpx - dpy * dpy - dpz * dpz;
891 //                  printf("New mass (1) %e %e %e %e %e %e %e \n", dmass2, jt, jtNew, pl, plNew, e, eNew);
892 //                  printf("New mass (2) %e %e \n", pxNew, pyNew);
893                     //
894                     //      Rotate back
895                     //  
896                     TVector3 w(pxNew, pyNew, plNew);
897                     w.RotateY(thetaq[index]); w.RotateZ(phiq[index]);
898                     pxNew = w.X(); pyNew = w.Y(); plNew = w.Z();
899                 
900                     p1[index][0] += pxNew;
901                     p1[index][1] += pyNew;
902                     p1[index][2] += plNew;
903                     p1[index][3] += eNew;       
904                     //
905                     // Updated 4-momentum vectors
906                     //
907                     pNew[icount][0]  = pxNew;
908                     pNew[icount][1]  = pyNew;
909                     pNew[icount][2]  = plNew;
910                     pNew[icount][3]  = eNew;
911                     kNew[icount]     = i;
912                     icount++;
913                 } // parton loop
914                 //
915                 // Check if there was phase-space for quenching
916                 //
917
918                 if (icount == 0) quenched[isys] = kFALSE;
919                 if (!quenched[isys]) break;
920                 
921                 for (Int_t j = 0; j < 4; j++) 
922                 {
923                     p2[isys][j] = p0[isys][j] - p1[isys][j];
924                 }
925                 p2[isys][4] = p2[isys][3] * p2[isys][3] - p2[isys][0] * p2[isys][0] - p2[isys][1] * p2[isys][1] - p2[isys][2] * p2[isys][2];
926                 if (p2[isys][4] > 0.) {
927                     p2[isys][4] = TMath::Sqrt(p2[isys][4]);
928                     break;
929                 } else {
930                     printf("Warning negative mass squared in system %d %f ! \n", isys, zquench[isys]);
931                     printf("4-Momentum: %10.3e %10.3e %10.3e %10.3e %10.3e \n", p2[isys][0], p2[isys][1], p2[isys][2], p2[isys][3], p2[isys][4]);
932                     if (p2[isys][4] < -0.01) {
933                         printf("Negative mass squared !\n");
934                         // Here we have to put the gluon back to mass shell
935                         // This will lead to a small energy imbalance
936                         p2[isys][4]  = 0.;
937                         p2[isys][3]  = TMath::Sqrt(p2[isys][0] * p2[isys][0] + p2[isys][1] * p2[isys][1] + p2[isys][2] * p2[isys][2]);
938                         break;
939                     } else {
940                         p2[isys][4] = 0.;
941                         break;
942                     }
943                 }
944                 /*
945                 zHeavy *= 0.98;
946                 printf("zHeavy lowered to %f\n", zHeavy);
947                 if (zHeavy < 0.01) {
948                     printf("No success ! \n");
949                     icount = 0;
950                     quenched[isys] = kFALSE;
951                     break;
952                 }
953                 */
954             } // iteration on z (while)
955             
956 //          Update  event record
957             for (Int_t k = 0; k < icount; k++) {
958 //              printf("%6d %6d %10.3e %10.3e %10.3e %10.3e\n", k, kNew[k], pNew[k][0],pNew[k][1], pNew[k][2], pNew[k][3] );
959                 fPyjets->P[0][kNew[k]] = pNew[k][0];
960                 fPyjets->P[1][kNew[k]] = pNew[k][1];
961                 fPyjets->P[2][kNew[k]] = pNew[k][2];
962                 fPyjets->P[3][kNew[k]] = pNew[k][3];
963             }
964             //
965             // Add the gluons
966             //
967             Int_t ish = 0;    
968             Int_t iGlu;
969             if (!quenched[isys]) continue;
970 //
971 //      Last parton from shower i
972             Int_t in = klast[isys];
973 //
974 //      Continue if no parton in shower i selected
975             if (in == -1) continue;
976 //  
977 //      If this is the second initial parton and it is behind the first move pointer by previous ish
978             if (isys == 1 && klast[1] > klast[0]) in += ish;
979 //
980 //      Starting index
981             
982 //          jmin = in - 1;
983 // How many additional gluons will be generated
984             ish  = 1;
985             if (p2[isys][4] > 0.05) ish = 2;
986 //
987 //      Position of gluons
988             iGlu = numpart;
989             if (iglu == 0) igMin = iGlu;
990             igMax = iGlu;
991             numpart += ish;
992             (fPyjets->N) += ish;
993             
994             if (ish == 1) {
995                 fPyjets->P[0][iGlu] = p2[isys][0];
996                 fPyjets->P[1][iGlu] = p2[isys][1];
997                 fPyjets->P[2][iGlu] = p2[isys][2];
998                 fPyjets->P[3][iGlu] = p2[isys][3];
999                 fPyjets->P[4][iGlu] = p2[isys][4];
1000                 
1001                 fPyjets->K[0][iGlu] = 1;
1002                 if (iglu == nGluon[isys] - 1) fPyjets->K[0][iGlu] = 1;
1003                 fPyjets->K[1][iGlu] = 21;       
1004                 fPyjets->K[2][iGlu] = fPyjets->K[2][in] + 1000;
1005                 fPyjets->K[3][iGlu] = -1;       
1006                 fPyjets->K[4][iGlu] = -1;
1007                 
1008                 pg[0] += p2[isys][0];
1009                 pg[1] += p2[isys][1];
1010                 pg[2] += p2[isys][2];
1011                 pg[3] += p2[isys][3];
1012             } else {
1013                 //
1014                 // Split gluon in rest frame.
1015                 //
1016                 Double_t bx   =  p2[isys][0] / p2[isys][3];
1017                 Double_t by   =  p2[isys][1] / p2[isys][3];
1018                 Double_t bz   =  p2[isys][2] / p2[isys][3];
1019                 Double_t pst  =  p2[isys][4] / 2.;
1020                 //
1021                 // Isotropic decay ????
1022                 Double_t cost = 2. * gRandom->Rndm() - 1.;
1023                 Double_t sint = TMath::Sqrt(1. - cost * cost);
1024                 Double_t phi =  2. * TMath::Pi() * gRandom->Rndm();
1025                 
1026                 Double_t pz1 =   pst * cost;
1027                 Double_t pz2 =  -pst * cost;
1028                 Double_t pt1 =   pst * sint;
1029                 Double_t pt2 =  -pst * sint;
1030                 Double_t px1 =   pt1 * TMath::Cos(phi);
1031                 Double_t py1 =   pt1 * TMath::Sin(phi);     
1032                 Double_t px2 =   pt2 * TMath::Cos(phi);
1033                 Double_t py2 =   pt2 * TMath::Sin(phi);     
1034                 
1035                 fPyjets->P[0][iGlu] = px1;
1036                 fPyjets->P[1][iGlu] = py1;
1037                 fPyjets->P[2][iGlu] = pz1;
1038                 fPyjets->P[3][iGlu] = pst;
1039                 fPyjets->P[4][iGlu] = 0.;
1040                 
1041                 fPyjets->K[0][iGlu] = 1 ;
1042                 fPyjets->K[1][iGlu] = 21;       
1043                 fPyjets->K[2][iGlu] = fPyjets->K[2][in] + 1000;
1044                 fPyjets->K[3][iGlu] = -1;       
1045                 fPyjets->K[4][iGlu] = -1;
1046                 
1047                 fPyjets->P[0][iGlu+1] = px2;
1048                 fPyjets->P[1][iGlu+1] = py2;
1049                 fPyjets->P[2][iGlu+1] = pz2;
1050                 fPyjets->P[3][iGlu+1] = pst;
1051                 fPyjets->P[4][iGlu+1] = 0.;
1052                 
1053                 fPyjets->K[0][iGlu+1] = 1;
1054                 if (iglu == nGluon[isys] - 1) fPyjets->K[0][iGlu+1] = 1;
1055                 fPyjets->K[1][iGlu+1] = 21;     
1056                 fPyjets->K[2][iGlu+1] = fPyjets->K[2][in] + 1000;
1057                 fPyjets->K[3][iGlu+1] = -1;     
1058                 fPyjets->K[4][iGlu+1] = -1;
1059                 SetMSTU(1,0);
1060                 SetMSTU(2,0);
1061                 //
1062                 // Boost back
1063                 //
1064                 Pyrobo(iGlu + 1, iGlu + 2, 0., 0., bx, by, bz);
1065             }
1066 /*    
1067             for (Int_t ig = iGlu; ig < iGlu+ish; ig++) {
1068                 Double_t px, py, pz;
1069                 px = fPyjets->P[0][ig]; 
1070                 py = fPyjets->P[1][ig]; 
1071                 pz = fPyjets->P[2][ig]; 
1072                 TVector3 v(px, py, pz);
1073                 v.RotateZ(-phiq[isys]);
1074                 v.RotateY(-thetaq[isys]);
1075                 Double_t pxs     = v.X(); Double_t pys = v.Y(); Double_t pzs  = v.Z();     
1076                 Double_t r       = AliPythiaRndm::GetPythiaRandom()->Rndm();
1077                 Double_t jtKick  = 0.3 * TMath::Sqrt(-TMath::Log(r));
1078                 if (ish == 2)   jtKick  = wjtKick[i] * TMath::Sqrt(-TMath::Log(r)) / TMath::Sqrt(2.);
1079                 Double_t phiKick = 2. * TMath::Pi() * AliPythiaRndm::GetPythiaRandom()->Rndm();
1080                 pxs += jtKick * TMath::Cos(phiKick);
1081                 pys += jtKick * TMath::Sin(phiKick);
1082                 TVector3 w(pxs, pys, pzs);
1083                 w.RotateY(thetaq[isys]);
1084                 w.RotateZ(phiq[isys]);
1085                 fPyjets->P[0][ig] = w.X(); 
1086                 fPyjets->P[1][ig] = w.Y(); 
1087                 fPyjets->P[2][ig] = w.Z(); 
1088                 fPyjets->P[2][ig] = w.Mag();
1089             }
1090 */
1091         } // kGluon         
1092         
1093         
1094     // Check energy conservation
1095         Double_t pxs = 0.;
1096         Double_t pys = 0.;
1097         Double_t pzs = 0.;      
1098         Double_t es  = 14000.;
1099         
1100         for (Int_t i = 0; i < numpart; i++)
1101         {
1102             kst =  fPyjets->K[0][i];
1103             if (kst != 1 && kst != 2) continue;
1104             pxs += fPyjets->P[0][i];
1105             pys += fPyjets->P[1][i];
1106             pzs += fPyjets->P[2][i];        
1107             es  -= fPyjets->P[3][i];        
1108         }
1109         if (TMath::Abs(pxs) > 1.e-2 ||
1110             TMath::Abs(pys) > 1.e-2 ||
1111             TMath::Abs(pzs) > 1.e-1) {
1112             printf("%e %e %e %e\n", pxs, pys, pzs, es);
1113 //              Fatal("Quench()", "4-Momentum non-conservation");
1114         }
1115         
1116     } // end quenching loop (systems)
1117 // Clean-up
1118     for (Int_t i = 0; i < numpart; i++)
1119     {
1120         imo =  fPyjets->K[2][i];
1121         if (imo > 1000) {
1122             fPyjets->K[2][i] = fPyjets->K[2][i] % 1000;
1123         }
1124     }
1125 //      this->Pylist(1);
1126 } // end quench
1127
1128
1129 void AliPythia::Pyquen(Double_t a, Int_t ibf, Double_t b)
1130 {
1131     // Igor Lokthine's quenching routine
1132     pyquen(a, ibf, b);
1133 }
1134
1135 void AliPythia::Pyevnw()
1136 {
1137     // New multiple interaction scenario
1138     pyevnw();
1139 }
1140
1141 void AliPythia::GetQuenchingParameters(Double_t& xp, Double_t& yp, Double_t z[4])
1142 {
1143     // Return event specific quenching parameters
1144     xp = fXJet;
1145     yp = fYJet;
1146     for (Int_t i = 0; i < 4; i++) z[i] = fZQuench[i];
1147
1148 }
1149
1150 void AliPythia::ConfigHeavyFlavor()
1151 {
1152     //
1153     // Default configuration for Heavy Flavor production
1154     //
1155     // All QCD processes
1156     //
1157     SetMSEL(1);
1158     
1159     // No multiple interactions
1160     SetMSTP(81,0);
1161     SetPARP(81,0.0);
1162     SetPARP(82,0.0);
1163     
1164     // Initial/final parton shower on (Pythia default)
1165     SetMSTP(61,1);
1166     SetMSTP(71,1);
1167     
1168     // 2nd order alpha_s
1169     SetMSTP(2,2);
1170     
1171     // QCD scales
1172     SetMSTP(32,2);
1173     SetPARP(34,1.0);
1174 }
1175
1176 void AliPythia::AtlasTuning()
1177 {
1178     //
1179     // Configuration for the ATLAS tuning
1180         SetMSTP(51, kCTEQ5L);      // CTEQ5L pdf
1181         SetMSTP(81,1);             // Multiple Interactions ON
1182         SetMSTP(82,4);             // Double Gaussian Model
1183         SetPARP(82,1.8);           // [GeV]    PT_min at Ref. energy
1184         SetPARP(89,1000.);         // [GeV]   Ref. energy
1185         SetPARP(90,0.16);          // 2*epsilon (exponent in power law)
1186         SetPARP(83,0.5);           // Core density in proton matter distribution (def.value)
1187         SetPARP(84,0.5);           // Core radius
1188         SetPARP(85,0.33);          // Regulates gluon prod. mechanism
1189         SetPARP(86,0.66);          // Regulates gluon prod. mechanism
1190         SetPARP(67,1);             // Regulates Initial State Radiation
1191 }