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