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