Implement setters for weighted event generation and pyevwt
[u/mrichter/AliRoot.git] / PYTHIA6 / AliPythia6.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: AliPythia.cxx,v 1.40 2007/10/09 08:43:24 morsch Exp $ */
17
18 #include "AliPythia6.h"
19 #include "AliStack.h"
20 #include "AliPythiaRndm.h"
21 #include "AliFastGlauber.h"
22 #include "AliQuenchingWeights.h"
23
24 #include "TVector3.h"
25 #include "TParticle.h"
26 #include "PyquenCommon.h"
27
28 ClassImp(AliPythia6)
29
30 #ifndef WIN32
31 # define pyclus pyclus_
32 # define pycell pycell_
33 # define pyshow pyshow_
34 # define pyshowq pyshowq_
35 # define pyrobo pyrobo_
36 # define pyquen pyquen_
37 # define pyevnw pyevnw_
38 # define pyjoin pyjoin_
39 # define qpygin0 qpygin0_
40 # define setpowwght setpowwght_
41 # define type_of_call
42 #else
43 # define pyclus PYCLUS
44 # define pycell PYCELL
45 # define pyshow PYSHOW
46 # define pyshowq PYSHOWQ
47 # define pyrobo PYROBO
48 # define pyquen PYQUEN
49 # define pyevnw PYEVNW
50 # define pyjoin PYJOIN
51 # define qpygin0 QPYGIN0
52 # define setpowwght SETPOWWGHT
53 # define type_of_call _stdcall
54 #endif
55
56 extern "C" void type_of_call pyjoin(Int_t &, Int_t * );
57 extern "C" void type_of_call pyclus(Int_t & );
58 extern "C" void type_of_call pycell(Int_t & );
59 extern "C" void type_of_call pyshow(Int_t &, Int_t &, Double_t &);
60 extern "C" void type_of_call pyshowq(Int_t &, Int_t &, Double_t &);
61 extern "C" void type_of_call qpygin0();
62 extern "C" void type_of_call pyrobo(Int_t &, Int_t &, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &);
63 extern "C" void type_of_call pyquen(Double_t &, Int_t &, Double_t &);
64 extern "C" void type_of_call pyevnw();
65 extern "C" void type_of_call setpowwght(Double_t &);
66
67
68 //_____________________________________________________________________________
69
70 AliPythia6* AliPythia6::fgAliPythia=NULL;
71
72 AliPythia6::AliPythia6():
73     TPythia6(),
74     AliPythiaBase(),
75     fProcess(kPyMb),
76     fEcms(0.),
77     fStrucFunc(kCTEQ5L),
78     fProjectile("p"),
79     fTarget("p"),
80     fXJet(0.),
81     fYJet(0.),
82     fNGmax(30),
83     fZmax(0.97),
84     fGlauber(0),
85     fQuenchingWeights(0)
86 {
87 // Default Constructor
88 //
89 //  Set random number
90     Int_t i;
91     for (i = 0; i <  501; i++) fDefMDCY[i] = 0;
92     for (i = 0; i < 2001; i++) fDefMDME[i] = 0;
93     for (i = 0; i <    4; i++) fZQuench[i] = 0;
94
95     if (!AliPythiaRndm::GetPythiaRandom()) 
96       AliPythiaRndm::SetPythiaRandom(GetRandom());
97     fGlauber          = 0;
98     fQuenchingWeights = 0;
99 }
100
101 AliPythia6::AliPythia6(const AliPythia6& pythia):
102     TPythia6(),
103     AliPythiaBase(),
104     fProcess(kPyMb),
105     fEcms(0.),
106     fStrucFunc(kCTEQ5L),
107     fProjectile("p"),
108     fTarget("p"),
109     fXJet(0.),
110     fYJet(0.),
111     fNGmax(30),
112     fZmax(0.97),
113     fGlauber(0),
114     fQuenchingWeights(0)
115 {
116     // Copy Constructor
117     Int_t i;
118     for (i = 0; i <  501; i++) fDefMDCY[i] = 0;
119     for (i = 0; i < 2001; i++) fDefMDME[i] = 0;
120     for (i = 0; i <    4; i++) fZQuench[i] = 0;
121     pythia.Copy(*this);
122 }
123
124 void AliPythia6::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfunc, Int_t /*tune*/)
125 {
126 // Initialise the process to generate 
127     if (!AliPythiaRndm::GetPythiaRandom()) 
128       AliPythiaRndm::SetPythiaRandom(GetRandom());
129     
130     fProcess = process;
131     fEcms = energy;
132     fStrucFunc = strucfunc;
133 //...Switch off decay of pi0, K0S, Lambda, Sigma+-, Xi0-, Omega-.
134     SetMDCY(Pycomp(111) ,1,0);
135     SetMDCY(Pycomp(310) ,1,0);
136     SetMDCY(Pycomp(3122),1,0);
137     SetMDCY(Pycomp(3112),1,0);
138     SetMDCY(Pycomp(3212),1,0);
139     SetMDCY(Pycomp(3222),1,0);
140     SetMDCY(Pycomp(3312),1,0);
141     SetMDCY(Pycomp(3322),1,0);
142     SetMDCY(Pycomp(3334),1,0);
143     // Select structure function 
144     SetMSTP(52,2);
145     SetMSTP(51,AliStructFuncType::PDFsetIndex(strucfunc));
146     // Particles produced in string fragmentation point directly to either of the two endpoints
147     // of the string (depending in the side they were generated from).
148     SetMSTU(16,2);
149
150 //
151 // Pythia initialisation for selected processes//
152 //
153 // Make MSEL clean
154 //
155     for (Int_t i=1; i<= 200; i++) {
156         SetMSUB(i,0);
157     }
158 //  select charm production
159     switch (process) 
160     {
161     case kPyOldUEQ2ordered:  //Old underlying events with Q2 ordered QCD processes
162 //        Multiple interactions on.
163         SetMSTP(81,1);
164 // Double Gaussian matter distribution.
165         SetMSTP(82,4);
166         SetPARP(83,0.5);
167         SetPARP(84,0.4);
168 //  pT0.
169         SetPARP(82,2.0);
170 //  Reference energy for pT0 and energy rescaling pace.
171         SetPARP(89,1800);
172         SetPARP(90,0.25);
173 //  String drawing almost completely minimizes string length.
174         SetPARP(85,0.9);
175         SetPARP(86,0.95);
176 // ISR and FSR activity.
177         SetPARP(67,4);
178         SetPARP(71,4);
179 // Lambda_FSR scale.
180         SetPARJ(81,0.29);
181         break;
182     case kPyOldUEQ2ordered2:   
183 // Old underlying events with Q2 ordered QCD processes
184 // Multiple interactions on.
185         SetMSTP(81,1);
186 // Double Gaussian matter distribution.
187         SetMSTP(82,4);
188         SetPARP(83,0.5);
189         SetPARP(84,0.4);
190 // pT0.
191         SetPARP(82,2.0);
192 // Reference energy for pT0 and energy rescaling pace.
193         SetPARP(89,1800);
194         SetPARP(90,0.16);  // here is the difference with  kPyOldUEQ2ordered
195 // String drawing almost completely minimizes string length.
196         SetPARP(85,0.9);
197         SetPARP(86,0.95);
198 // ISR and FSR activity.
199         SetPARP(67,4);
200         SetPARP(71,4);
201 // Lambda_FSR scale.
202         SetPARJ(81,0.29);       
203         break;
204     case kPyOldPopcorn:  
205 // Old production mechanism: Old Popcorn
206         SetMSEL(1);
207         SetMSTJ(12,3); 
208 // (D=2) Like MSTJ(12)=2 but added prod ofthe 1er rank baryon
209         SetMSTP(88,2); 
210 // (D=1)see can be used to form  baryons (BARYON JUNCTION)
211         SetMSTJ(1,1);  
212         AtlasTuning();
213         break;
214     case kPyCharm:
215         SetMSEL(4);
216 //  heavy quark masses
217
218         SetPMAS(4,1,1.2);
219 //
220 //    primordial pT
221         SetMSTP(91,1);
222         SetPARP(91,1.);
223         SetPARP(93,5.);
224 //
225         break;
226     case kPyBeauty:
227         SetMSEL(5);
228         SetPMAS(5,1,4.75);
229         break;
230     case kPyJpsi:
231         SetMSEL(0);
232 // gg->J/Psi g
233         SetMSUB(86,1);
234         break;
235     case kPyJpsiChi:
236         SetMSEL(0);
237 // gg->J/Psi g
238         SetMSUB(86,1);
239 // gg-> chi_0c g
240         SetMSUB(87,1);
241 // gg-> chi_1c g
242         SetMSUB(88,1);
243 // gg-> chi_2c g
244         SetMSUB(89,1);  
245         break;
246     case kPyCharmUnforced:
247         SetMSEL(0);
248 // gq->qg   
249         SetMSUB(28,1);
250 // gg->qq
251         SetMSUB(53,1);
252 // gg->gg
253         SetMSUB(68,1);
254         break;
255     case kPyBeautyUnforced:
256         SetMSEL(0);
257 // gq->qg   
258         SetMSUB(28,1);
259 // gg->qq
260         SetMSUB(53,1);
261 // gg->gg
262         SetMSUB(68,1);
263         break;
264     case kPyMb:
265 // Minimum Bias pp-Collisions
266 //
267 //   
268 //      select Pythia min. bias model
269         SetMSEL(0);
270         SetMSUB(92,1);             // single diffraction AB-->XB
271         SetMSUB(93,1);             // single diffraction AB-->AX
272         SetMSUB(94,1);             // double diffraction
273         SetMSUB(95,1);             // low pt production
274
275         AtlasTuning();
276         break;
277     case kPyMbAtlasTuneMC09:
278 // Minimum Bias pp-Collisions
279 //
280 //   
281 //      select Pythia min. bias model
282         SetMSEL(0);
283         SetMSUB(92,1);             // single diffraction AB-->XB
284         SetMSUB(93,1);             // single diffraction AB-->AX
285         SetMSUB(94,1);             // double diffraction
286         SetMSUB(95,1);             // low pt production
287
288         AtlasTuningMC09();
289         break;
290
291     case kPyMbWithDirectPhoton:
292 // Minimum Bias pp-Collisions with direct photon processes added 
293 //
294 //   
295 //      select Pythia min. bias model
296         SetMSEL(0);
297         SetMSUB(92,1);             // single diffraction AB-->XB
298         SetMSUB(93,1);             // single diffraction AB-->AX
299         SetMSUB(94,1);             // double diffraction
300         SetMSUB(95,1);             // low pt production
301
302         SetMSUB(14,1);             //
303         SetMSUB(18,1);             //
304         SetMSUB(29,1);             //
305         SetMSUB(114,1);            //
306         SetMSUB(115,1);            //
307
308
309         AtlasTuning();
310         break;
311
312     case kPyMbDefault:
313 // Minimum Bias pp-Collisions
314 //
315 //   
316 //      select Pythia min. bias model
317         SetMSEL(0);
318         SetMSUB(92,1);             // single diffraction AB-->XB
319         SetMSUB(93,1);             // single diffraction AB-->AX
320         SetMSUB(94,1);             // double diffraction
321         SetMSUB(95,1);             // low pt production
322
323         break;
324     case kPyLhwgMb:
325 // Les Houches Working Group 05 Minimum Bias pp-Collisions: hep-ph/0604120
326 //  -> Pythia 6.3 or above is needed
327 //   
328         SetMSEL(0);
329         SetMSUB(92,1);             // single diffraction AB-->XB
330         SetMSUB(93,1);             // single diffraction AB-->AX
331         SetMSUB(94,1);             // double diffraction
332         SetMSUB(95,1);             // low pt production
333         SetMSTP(51,AliStructFuncType::PDFsetIndex(kCTEQ6ll));
334         SetMSTP(52,2);
335         SetMSTP(68,1);
336         SetMSTP(70,2);
337         SetMSTP(81,1);             // Multiple Interactions ON
338         SetMSTP(82,4);             // Double Gaussian Model
339         SetMSTP(88,1);
340
341         SetPARP(82,2.3);           // [GeV]    PT_min at Ref. energy
342         SetPARP(83,0.5);           // Core density in proton matter distribution (def.value)
343         SetPARP(84,0.5);           // Core radius
344         SetPARP(85,0.9);           // Regulates gluon prod. mechanism
345         SetPARP(90,0.2);           // 2*epsilon (exponent in power law)
346
347         break;
348     case kPyMbNonDiffr:
349 // Minimum Bias pp-Collisions
350 //
351 //   
352 //      select Pythia min. bias model
353         SetMSEL(0);
354         SetMSUB(95,1);             // low pt production
355
356         AtlasTuning();
357         break;
358     case kPyMbMSEL1:
359         ConfigHeavyFlavor();
360 // Intrinsic <kT^2>
361         SetMSTP(91,1);// Width (1=gaussian) primordial kT dist. inside hadrons
362         SetPARP(91,1.);     // <kT^2> = PARP(91,1.)^2
363         SetPARP(93,5.);     // Upper cut-off
364 // Set Q-quark mass
365         SetPMAS(4,1,1.2);   // Charm quark mass
366         SetPMAS(5,1,4.78);  // Beauty quark mass
367         SetPARP(71,4.);     // Defaut value
368 // Atlas Tuning
369         AtlasTuning();
370         break;
371     case kPyJets:
372 //
373 //  QCD Jets
374 //
375         SetMSEL(1);
376  // Pythia Tune A (CDF)
377  //
378        SetPARP(67,2.5);           // Regulates Initial State Radiation (value from best fit to D0 dijet analysis)
379        SetMSTP(82,4);             // Double Gaussian Model
380        SetPARP(82,2.0);           // [GeV]    PT_min at Ref. energy
381        SetPARP(84,0.4);           // Core radius
382        SetPARP(85,0.90) ;         // Regulates gluon prod. mechanism
383        SetPARP(86,0.95);          // Regulates gluon prod. mechanism
384        SetPARP(89,1800.);         // [GeV]   Ref. energy
385        SetPARP(90,0.25);          // 2*epsilon (exponent in power law)
386        break;
387     case kPyDirectGamma:
388         SetMSEL(10);
389         break;
390     case kPyCharmPbPbMNR:
391     case kPyD0PbPbMNR:
392     case kPyDPlusPbPbMNR:
393     case kPyDPlusStrangePbPbMNR:
394       // Tuning of Pythia parameters aimed to get a resonable agreement
395       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
396       // c-cbar single inclusive and double differential distributions.
397       // This parameter settings are meant to work with Pb-Pb collisions
398       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
399       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
400       // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
401         ConfigHeavyFlavor();
402       // Intrinsic <kT>
403       SetMSTP(91,1);
404       SetPARP(91,1.304);
405       SetPARP(93,6.52);
406       // Set c-quark mass
407       SetPMAS(4,1,1.2);
408       break;
409     case kPyCharmpPbMNR:
410     case kPyD0pPbMNR:
411     case kPyDPluspPbMNR:
412     case kPyDPlusStrangepPbMNR:
413       // Tuning of Pythia parameters aimed to get a resonable agreement
414       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
415       // c-cbar single inclusive and double differential distributions.
416       // This parameter settings are meant to work with p-Pb collisions
417       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
418       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
419       // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
420         ConfigHeavyFlavor();
421       // Intrinsic <kT>
422         SetMSTP(91,1);
423         SetPARP(91,1.16);
424         SetPARP(93,5.8);
425         
426       // Set c-quark mass
427         SetPMAS(4,1,1.2);
428       break;
429     case kPyCharmppMNR:
430     case kPyD0ppMNR:
431     case kPyDPlusppMNR:
432     case kPyDPlusStrangeppMNR:
433     case kPyLambdacppMNR:
434       // Tuning of Pythia parameters aimed to get a resonable agreement
435       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
436       // c-cbar single inclusive and double differential distributions.
437       // This parameter settings are meant to work with pp collisions
438       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
439       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
440       // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
441         ConfigHeavyFlavor();
442       // Intrinsic <kT^2>
443         SetMSTP(91,1);
444         SetPARP(91,1.);
445         SetPARP(93,5.);
446         
447       // Set c-quark mass
448         SetPMAS(4,1,1.2);
449       break;
450     case kPyCharmppMNRwmi:
451       // Tuning of Pythia parameters aimed to get a resonable agreement
452       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
453       // c-cbar single inclusive and double differential distributions.
454       // This parameter settings are meant to work with pp collisions
455       // and with kCTEQ5L PDFs.
456       // Added multiple interactions according to ATLAS tune settings.
457       // To get a "reasonable" agreement with MNR results, events have to be 
458       // generated with the minimum ptHard (AliGenPythia::SetPtHard)
459       // set to 2.76 GeV.
460       // To get a "perfect" agreement with MNR results, events have to be 
461       // generated in four ptHard bins with the following relative 
462       // normalizations:
463       // 2.76-3 GeV: 25%
464       //    3-4 GeV: 40%
465       //    4-8 GeV: 29%
466       //     >8 GeV:  6%
467         ConfigHeavyFlavor();
468       // Intrinsic <kT^2>
469         SetMSTP(91,1);
470         SetPARP(91,1.);
471         SetPARP(93,5.);
472
473       // Set c-quark mass
474         SetPMAS(4,1,1.2);
475         AtlasTuning();
476         break;
477     case kPyBeautyPbPbMNR:
478       // Tuning of Pythia parameters aimed to get a resonable agreement
479       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
480       // b-bbar single inclusive and double differential distributions.
481       // This parameter settings are meant to work with Pb-Pb collisions
482       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
483       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
484       // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
485         ConfigHeavyFlavor();
486       // QCD scales
487         SetPARP(67,1.0);
488         SetPARP(71,1.0);
489       // Intrinsic <kT>
490         SetMSTP(91,1);
491         SetPARP(91,2.035);
492         SetPARP(93,10.17);
493       // Set b-quark mass
494         SetPMAS(5,1,4.75);
495       break;
496     case kPyBeautypPbMNR:
497       // Tuning of Pythia parameters aimed to get a resonable agreement
498       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
499       // b-bbar single inclusive and double differential distributions.
500       // This parameter settings are meant to work with p-Pb collisions
501       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
502       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
503       // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
504         ConfigHeavyFlavor();
505       // QCD scales
506         SetPARP(67,1.0);
507         SetPARP(71,1.0);
508       // Intrinsic <kT>
509         SetMSTP(91,1);
510         SetPARP(91,1.60);
511         SetPARP(93,8.00);
512       // Set b-quark mass
513         SetPMAS(5,1,4.75);
514       break;
515     case kPyBeautyppMNR:
516       // Tuning of Pythia parameters aimed to get a resonable agreement
517       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
518       // b-bbar single inclusive and double differential distributions.
519       // This parameter settings are meant to work with pp collisions
520       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
521       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
522       // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
523         ConfigHeavyFlavor();
524       // QCD scales
525         SetPARP(67,1.0);
526         SetPARP(71,1.0);
527         
528         // Intrinsic <kT>
529         SetMSTP(91,1);
530         SetPARP(91,1.);
531         SetPARP(93,5.);
532         
533         // Set b-quark mass
534         SetPMAS(5,1,4.75);
535       break;
536      case kPyBeautyJets: 
537      case kPyBeautyppMNRwmi:
538       // Tuning of Pythia parameters aimed to get a resonable agreement
539       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
540       // b-bbar single inclusive and double differential distributions.
541       // This parameter settings are meant to work with pp collisions
542       // and with kCTEQ5L PDFs.
543       // Added multiple interactions according to ATLAS tune settings.
544       // To get a "reasonable" agreement with MNR results, events have to be 
545       // generated with the minimum ptHard (AliGenPythia::SetPtHard)
546       // set to 2.76 GeV.
547       // To get a "perfect" agreement with MNR results, events have to be 
548       // generated in four ptHard bins with the following relative 
549       // normalizations:
550       // 2.76-4 GeV:  5% 
551       //    4-6 GeV: 31%
552       //    6-8 GeV: 28%
553       //     >8 GeV: 36%
554          ConfigHeavyFlavor();
555       // QCD scales
556          SetPARP(67,1.0);
557          SetPARP(71,1.0);
558          
559          // Intrinsic <kT>
560          SetMSTP(91,1);
561          SetPARP(91,1.);
562          SetPARP(93,5.);
563
564       // Set b-quark mass
565          SetPMAS(5,1,4.75);
566
567          AtlasTuning();
568          break; 
569     case kPyW:
570
571       //Inclusive production of W+/-
572       SetMSEL(0);
573       //f fbar -> W+ 
574       SetMSUB(2,1);
575       //        //f fbar -> g W+
576       //        SetMSUB(16,1);
577       //        //f fbar -> gamma W+
578       //        SetMSUB(20,1);
579       //        //f g -> f W+  
580       //        SetMSUB(31,1);
581       //        //f gamma -> f W+
582       //        SetMSUB(36,1);
583       
584       // Initial/final parton shower on (Pythia default)
585       // With parton showers on we are generating "W inclusive process"
586       SetMSTP(61,1); //Initial QCD & QED showers on
587       SetMSTP(71,1); //Final QCD & QED showers on
588       
589       break;  
590
591     case kPyZ:
592
593       //Inclusive production of Z
594       SetMSEL(0);
595       //f fbar -> Z/gamma
596       SetMSUB(1,1);
597       
598       //       // f fbar -> g Z/gamma
599       //       SetMSUB(15,1);
600       //       // f fbar -> gamma Z/gamma
601       //       SetMSUB(19,1);
602       //       // f g -> f Z/gamma
603       //       SetMSUB(30,1);
604       //       // f gamma -> f Z/gamma
605       //       SetMSUB(35,1);
606       
607       //only Z included, not gamma
608       SetMSTP(43,2);
609       
610       // Initial/final parton shower on (Pythia default)
611       // With parton showers on we are generating "Z inclusive process"
612       SetMSTP(61,1); //Initial QCD & QED showers on
613       SetMSTP(71,1); //Final QCD & QED showers on
614       break;
615     case kPyZgamma:
616       //Inclusive production of Z
617       SetMSEL(0);
618       //f fbar -> Z/gamma
619       SetMSUB(1,1);
620       // Initial/final parton shower on (Pythia default)
621       // With parton showers on we are generating "Z inclusive process"
622       SetMSTP(61,1); //Initial QCD & QED showers on
623       SetMSTP(71,1); //Final QCD & QED showers on
624       break;
625       case kPyMBRSingleDiffraction:
626       case kPyMBRDoubleDiffraction:
627       case kPyMBRCentralDiffraction:
628       break;  
629       case kPyJetsPWHG:
630       //    N.B.
631       //    ====
632       //    For the case of jet production the following parameter setting
633       //    limits the transverse momentum of secondary scatterings, due
634       //    to multiple parton interactions, to be less than that of the
635       //    primary interaction (see POWHEG Dijet paper arXiv:1012.3380
636       //    [hep-ph] sec. 4.1 and also the PYTHIA Manual).
637       SetMSTP(86,1);
638       //    maximum number of errors before pythia aborts (def=10)
639       SetMSTU(22,10);
640       //    number of warnings printed on the shell
641       SetMSTU(26,20);
642       break;
643       case kPyCharmPWHG:
644       case kPyBeautyPWHG:
645       case kPyWPWHG:
646       //    number of warnings printed on the shell
647       SetMSTU(26,20);
648
649     }
650 //
651 //  Initialize PYTHIA
652     SetMSTP(41,1);   // all resonance decays switched on
653     if (process == kPyJetsPWHG || process == kPyCharmPWHG || process == kPyBeautyPWHG || process == kPyWPWHG) {
654       Initialize("USER","","",0.);
655     } else {    
656       Initialize("CMS",fProjectile,fTarget,fEcms);
657     }
658 }
659
660 Int_t AliPythia6::CheckedLuComp(Int_t kf)
661 {
662 // Check Lund particle code (for debugging)
663     Int_t kc=Pycomp(kf);
664     return kc;
665 }
666
667 void AliPythia6::SetNuclei(Int_t a1, Int_t a2)
668 {
669 // Treat protons as inside nuclei with mass numbers a1 and a2  
670 //    The MSTP array in the PYPARS common block is used to enable and 
671 //    select the nuclear structure functions. 
672 //    MSTP(52)  : (D=1) choice of proton and nuclear structure-function library
673 //            =1: internal PYTHIA acording to MSTP(51) 
674 //            =2: PDFLIB proton  s.f., with MSTP(51)  = 1000xNGROUP+NSET
675 //    If the following mass number both not equal zero, nuclear corrections of the stf are used.
676 //    MSTP(192) : Mass number of nucleus side 1
677 //    MSTP(193) : Mass number of nucleus side 2
678     SetMSTP(52,2);
679     SetMSTP(192, a1);
680     SetMSTP(193, a2);  
681 }
682         
683
684 AliPythia6* AliPythia6::Instance()
685
686 // Set random number generator 
687     if (fgAliPythia) {
688         return fgAliPythia;
689     } else {
690         fgAliPythia = new AliPythia6();
691         return fgAliPythia;
692     }
693 }
694
695 void AliPythia6::PrintParticles()
696
697 // Print list of particl properties
698     Int_t np = 0;
699     char*   name = new char[16];    
700     for (Int_t kf=0; kf<1000000; kf++) {
701         for (Int_t c = 1;  c > -2; c-=2) {
702             Int_t kc = Pycomp(c*kf);
703             if (kc) {
704                 Float_t mass  = GetPMAS(kc,1);
705                 Float_t width = GetPMAS(kc,2);  
706                 Float_t tau   = GetPMAS(kc,4);
707
708                 Pyname(kf,name);
709         
710                 np++;
711                 
712                 printf("\n mass, width, tau: %6d %s %10.3f %10.3e %10.3e", 
713                        c*kf, name, mass, width, tau);
714             }
715         }
716     }
717     printf("\n Number of particles %d \n \n", np);
718 }
719
720 void  AliPythia6::ResetDecayTable()
721 {
722 //  Set default values for pythia decay switches
723     Int_t i;
724     for (i = 1; i <  501; i++) SetMDCY(i,1,fDefMDCY[i]);
725     for (i = 1; i < 2001; i++) SetMDME(i,1,fDefMDME[i]);
726 }
727
728 void  AliPythia6::SetDecayTable()
729 {
730 //  Set default values for pythia decay switches
731 //
732     Int_t i;
733     for (i = 1; i <  501; i++) fDefMDCY[i] = GetMDCY(i,1);
734     for (i = 1; i < 2001; i++) fDefMDME[i] = GetMDME(i,1);
735 }
736
737 void  AliPythia6::Pyjoin(Int_t& npart, Int_t *ipart)
738 {
739 //  Call Pythia join alogorithm to set up a string between
740 //  npart partons, given by indices in array ipart[npart]
741 //
742     pyjoin(npart, ipart);
743 }
744
745 void  AliPythia6::Pyshowq(Int_t ip1, Int_t ip2, Double_t qmax)
746 {
747 //  Call qPythia showering
748 //
749     pyshowq(ip1, ip2, qmax);
750 }
751
752 void AliPythia6::Qpygin0()
753 {
754     //position of the hard scattering in the nuclear overlapping area.
755     //just for qpythia.
756     qpygin0();
757 }
758
759 void  AliPythia6::Pyclus(Int_t& njet)
760 {
761 //  Call Pythia clustering algorithm
762 //
763     pyclus(njet);
764 }
765
766 void  AliPythia6::Pycell(Int_t& njet)
767 {
768 //  Call Pythia jet reconstruction algorithm
769 //
770     pycell(njet);
771 }
772
773 void AliPythia6::GetJet(Int_t i, Float_t& px, Float_t& py, Float_t& pz, Float_t& e)
774 {
775     // Get jet number i
776     Int_t n = GetN();
777     px    = GetPyjets()->P[0][n+i];
778     py    = GetPyjets()->P[1][n+i];
779     pz    = GetPyjets()->P[2][n+i];
780     e     = GetPyjets()->P[3][n+i];
781 }
782
783 void  AliPythia6::Pyshow(Int_t ip1, Int_t ip2, Double_t qmax)
784 {
785 //  Call Pythia showering
786 //
787     pyshow(ip1, ip2, qmax);
788 }
789
790 void AliPythia6::Pyrobo(Int_t imi, Int_t ima, Double_t the, Double_t phi, Double_t bex, Double_t bey, Double_t bez)
791 {
792     pyrobo(imi, ima, the, phi, bex, bey, bez);
793 }
794
795
796
797 void AliPythia6::InitQuenching(Float_t cMin, Float_t cMax, Float_t k, Int_t iECMethod, Float_t zmax, Int_t ngmax)
798 {
799 // Initializes 
800 // (1) The quenching model using quenching weights according to C. Salgado and U. Wiedemann
801 // (2) The nuclear geometry using the Glauber Model
802 //     
803     
804     fGlauber = AliFastGlauber::Instance();
805     fGlauber->Init(2);
806     fGlauber->SetCentralityClass(cMin, cMax); 
807
808     fQuenchingWeights = new AliQuenchingWeights();
809     fQuenchingWeights->InitMult();
810     fQuenchingWeights->SetK(k);
811     fQuenchingWeights->SetECMethod(AliQuenchingWeights::kECMethod(iECMethod));
812     fNGmax = ngmax;
813     fZmax  = zmax;
814     
815 }
816
817
818 void  AliPythia6::Quench()
819 {
820 //
821 //
822 //  Simple Jet Quenching routine:
823 //  =============================
824 //  The jet formed by all final state partons radiated by the parton created 
825 //  in the hard collisions is quenched by a factor (1-z) using light cone variables in 
826 //  the initial parton reference frame:
827 //  (E + p_z)new = (1-z) (E + p_z)old
828 //
829 //
830 //
831 //
832 //   The lost momentum is first balanced by one gluon with virtuality > 0.   
833 //   Subsequently the gluon splits to yield two gluons with E = p.
834 //
835 //
836 // 
837     static Float_t eMean = 0.;
838     static Int_t   icall = 0;
839     
840     Double_t p0[4][5];
841     Double_t p1[4][5];
842     Double_t p2[4][5];
843     Int_t   klast[4] = {-1, -1, -1, -1};
844
845     Int_t numpart   = fPyjets->N;
846     Double_t px = 0., py = 0., pz = 0., e = 0., m = 0., p = 0., pt = 0., theta = 0., phi = 0.;
847     Double_t pxq[4], pyq[4], pzq[4], eq[4], yq[4], mq[4], pq[4], phiq[4], thetaq[4], ptq[4];
848     Bool_t  quenched[4];
849     Double_t wjtKick[4] = {0., 0., 0., 0.};
850     Int_t nGluon[4];
851     Int_t qPdg[4];
852     Int_t   imo, kst, pdg;
853     
854 //
855 //  Sore information about Primary partons
856 //
857 //  j =
858 //  0, 1 partons from hard scattering
859 //  2, 3 partons from initial state radiation
860 // 
861     for (Int_t i = 2; i <= 7; i++) {
862         Int_t j = 0;
863         // Skip gluons that participate in hard scattering
864         if (i == 4 || i == 5) continue;
865         // Gluons from hard Scattering
866         if (i == 6 || i == 7) {
867             j = i - 6;
868             pxq[j]    = fPyjets->P[0][i];
869             pyq[j]    = fPyjets->P[1][i];
870             pzq[j]    = fPyjets->P[2][i];
871             eq[j]     = fPyjets->P[3][i];
872             mq[j]     = fPyjets->P[4][i];
873         } else {
874             // Gluons from initial state radiation
875             //
876             // Obtain 4-momentum vector from difference between original parton and parton after gluon 
877             // radiation. Energy is calculated independently because initial state radition does not 
878             // conserve strictly momentum and energy for each partonic system independently.
879             //
880             // Not very clean. Should be improved !
881             //
882             //
883             j = i;
884             pxq[j]    = fPyjets->P[0][i] - fPyjets->P[0][i+2];
885             pyq[j]    = fPyjets->P[1][i] - fPyjets->P[1][i+2];
886             pzq[j]    = fPyjets->P[2][i] - fPyjets->P[2][i+2];
887             mq[j]     = fPyjets->P[4][i];
888             eq[j]     = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j] + pzq[j] * pzq[j] + mq[j] * mq[j]);
889         }
890 //
891 //  Calculate some kinematic variables
892 //
893         yq[j]     = 0.5 * TMath::Log((eq[j] + pzq[j] + 1.e-14) / (eq[j] - pzq[j] + 1.e-14));
894         pq[j]     = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j] + pzq[j] * pzq[j]);
895         phiq[j]   = TMath::Pi()+TMath::ATan2(-pyq[j], -pxq[j]);
896         ptq[j]    = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j]);
897         thetaq[j] = TMath::ATan2(ptq[j], pzq[j]);
898         qPdg[j]   =  fPyjets->K[1][i];
899     }
900   
901     Double_t int0[4];
902     Double_t int1[4];
903     
904     fGlauber->GetI0I1ForPythiaAndXY(4, phiq, int0, int1, fXJet, fYJet, 15.);
905
906     for (Int_t j = 0; j < 4; j++) {
907         //
908         // Quench only central jets and with E > 10.
909         //
910
911
912         Int_t itype = (qPdg[j] == 21) ? 2 : 1;
913         //      Double_t eloss = fQuenchingWeights->GetELossRandomKFast(itype, int0[j], int1[j], eq[j]);
914         Double_t eloss = fQuenchingWeights->GetELossRandomK(itype, int0[j], int1[j], eq[j]);
915
916         if (TMath::Abs(yq[j]) > 2.5 || eq[j] < 10.) {
917             fZQuench[j] = 0.;
918         } else {
919             if (eq[j] > 40. && TMath::Abs(yq[j]) < 0.5) {
920                 icall ++;
921                 eMean += eloss;
922             }
923             //
924             // Extra pt
925             Double_t l =   fQuenchingWeights->CalcLk(int0[j], int1[j]);     
926             wjtKick[j] = TMath::Sqrt(l *  fQuenchingWeights->CalcQk(int0[j], int1[j]));
927             //
928             // Fractional energy loss
929             fZQuench[j] = eloss / eq[j];
930             //
931             // Avoid complete loss
932             //
933             if (fZQuench[j] > fZmax) fZQuench[j] = fZmax;
934             //
935             // Some debug printing
936
937             
938 //          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", 
939 //                 j, itype, eq[j], phiq[j], l, eloss, wjtKick[j], eMean / Float_t(icall+1), yq[j]);
940             
941 //          fZQuench[j] = 0.8;
942 //          while (fZQuench[j] >= 0.95)  fZQuench[j] = gRandom->Exp(0.2);
943         }
944         
945         quenched[j] = (fZQuench[j] > 0.01);
946     } // primary partons
947     
948     
949
950     Double_t pNew[1000][4];
951     Int_t    kNew[1000];
952     Int_t icount = 0;
953     Double_t zquench[4];
954     
955 //
956 //  System Loop    
957     for (Int_t isys = 0; isys < 4; isys++) {
958 //      Skip to next system if not quenched.
959         if (!quenched[isys]) continue;
960         
961         nGluon[isys]   = 1 + Int_t(fZQuench[isys] / (1. - fZQuench[isys]));
962         if (nGluon[isys] > fNGmax) nGluon[isys] = fNGmax;
963         zquench[isys] = 1. - TMath::Power(1. - fZQuench[isys], 1./Double_t(nGluon[isys]));
964         wjtKick[isys]  = wjtKick[isys] / TMath::Sqrt(Double_t(nGluon[isys]));
965
966
967         
968         Int_t igMin = -1;
969         Int_t igMax = -1;
970         Double_t pg[4] = {0., 0., 0., 0.};
971         
972 //
973 // Loop on radiation events
974
975         for (Int_t iglu = 0; iglu < nGluon[isys]; iglu++) {
976             while (1) {
977                 icount = 0;
978                 for (Int_t k = 0; k < 4; k++)
979                 {
980                     p0[isys][k] = 0.;
981                     p1[isys][k] = 0.;
982                     p2[isys][k] = 0.;
983                 }
984 //      Loop over partons
985                 for (Int_t i = 0; i < numpart; i++)
986                 {
987                     imo =  fPyjets->K[2][i];
988                     kst =  fPyjets->K[0][i];
989                     pdg =  fPyjets->K[1][i];
990                     
991                 
992                 
993 //      Quarks and gluons only
994                     if (pdg != 21 && TMath::Abs(pdg) > 6) continue;
995 //      Particles from hard scattering only
996                     
997                     if (imo > 8 && imo < 1000) imo = fPyjets->K[2][imo - 1];
998                     Int_t imom = imo % 1000;
999                     if ((isys == 0 || isys == 1) && ((imom != (isys + 7)))) continue;
1000                     if ((isys == 2 || isys == 3) && ((imom != (isys + 1)))) continue;               
1001                     
1002                     
1003 //      Skip comment lines
1004                     if (kst != 1 && kst != 2) continue;
1005 //
1006 //      Parton kinematic
1007                     px    = fPyjets->P[0][i];
1008                     py    = fPyjets->P[1][i];
1009                     pz    = fPyjets->P[2][i];
1010                     e     = fPyjets->P[3][i];
1011                     m     = fPyjets->P[4][i];
1012                     pt    = TMath::Sqrt(px * px + py * py);
1013                     p     = TMath::Sqrt(px * px + py * py + pz * pz); 
1014                     phi   = TMath::Pi() + TMath::ATan2(-py, -px);
1015                     theta = TMath::ATan2(pt, pz);
1016                 
1017 //
1018 //      Save 4-momentum sum for balancing
1019                     Int_t index = isys;
1020                     
1021                     p0[index][0] += px;
1022                     p0[index][1] += py;
1023                     p0[index][2] += pz;
1024                     p0[index][3] += e;
1025                 
1026                     klast[index] = i;
1027                     
1028 //
1029 //      Fractional energy loss
1030                     Double_t z = zquench[index];
1031                     
1032                     
1033 //      Don't fully quench radiated gluons
1034 //
1035                     if (imo > 1000) {
1036 //      This small factor makes sure that the gluons are not too close in phase space to avoid recombination
1037 //
1038
1039                         z = 0.02;
1040                     }
1041 //                  printf("z: %d %f\n", imo, z);
1042                     
1043
1044 //
1045                     
1046                     //
1047                     //
1048                     //      Transform into frame in which initial parton is along z-axis
1049                     //
1050                     TVector3 v(px, py, pz);
1051                     v.RotateZ(-phiq[index]);  v.RotateY(-thetaq[index]);
1052                     Double_t pxs = v.X(); Double_t pys = v.Y(); Double_t pl  = v.Z();
1053
1054                     Double_t jt  = TMath::Sqrt(pxs * pxs + pys * pys);
1055                     Double_t mt2 = jt * jt + m * m;
1056                     Double_t zmax = 1.;     
1057                     //
1058                     // Kinematic limit on z
1059                     //
1060                     if (m > 0.) zmax = 1. - m / TMath::Sqrt(m * m + jt * jt);
1061                     //
1062                     // Change light-cone kinematics rel. to initial parton
1063                     //  
1064                     Double_t eppzOld = e + pl;
1065                     Double_t empzOld = e - pl;
1066                     
1067                     Double_t eppzNew = (1. - z) * eppzOld;
1068                     Double_t empzNew = empzOld - mt2 * z / eppzOld;
1069                     Double_t eNew    = 0.5 * (eppzNew + empzNew);
1070                     Double_t plNew   = 0.5 * (eppzNew - empzNew);
1071                     
1072                     Double_t jtNew;
1073                     //
1074                     // if mt very small (or sometimes even < 0 for numerical reasons) set it to 0
1075                     Double_t mt2New = eppzNew * empzNew;
1076                     if (mt2New < 1.e-8) mt2New = 0.;
1077                     if (z < zmax) {
1078                         if (m * m > mt2New) {
1079                             //
1080                             // This should not happen 
1081                             //
1082                             Fatal("Quench()", "This should never happen %e %e %e!", m, eppzNew, empzNew);
1083                             jtNew = 0;
1084                         } else {
1085                             jtNew    = TMath::Sqrt(mt2New - m * m);
1086                         }
1087                     } else {
1088                         // If pT is to small (probably a leading massive particle) we scale only the energy
1089                         // This can cause negative masses of the radiated gluon
1090                         // Let's hope for the best ...
1091                         jtNew = jt;
1092                         eNew  = TMath::Sqrt(plNew * plNew + mt2);
1093                         
1094                     }
1095                     //
1096                     //     Calculate new px, py
1097                     //
1098                     Double_t pxNew = 0;
1099                     Double_t pyNew = 0;
1100
1101                     if (jt > 0.) {
1102                         pxNew   = jtNew / jt * pxs;
1103                         pyNew   = jtNew / jt * pys;     
1104                     }
1105                     
1106 //                  Double_t dpx = pxs - pxNew;
1107 //                  Double_t dpy = pys - pyNew;
1108 //                  Double_t dpz = pl  - plNew;
1109 //                  Double_t de  = e   - eNew;
1110 //                  Double_t dmass2 = de * de  - dpx * dpx - dpy * dpy - dpz * dpz;
1111 //                  printf("New mass (1) %e %e %e %e %e %e %e \n", dmass2, jt, jtNew, pl, plNew, e, eNew);
1112 //                  printf("New mass (2) %e %e \n", pxNew, pyNew);
1113                     //
1114                     //      Rotate back
1115                     //  
1116                     TVector3 w(pxNew, pyNew, plNew);
1117                     w.RotateY(thetaq[index]); w.RotateZ(phiq[index]);
1118                     pxNew = w.X(); pyNew = w.Y(); plNew = w.Z();
1119                 
1120                     p1[index][0] += pxNew;
1121                     p1[index][1] += pyNew;
1122                     p1[index][2] += plNew;
1123                     p1[index][3] += eNew;       
1124                     //
1125                     // Updated 4-momentum vectors
1126                     //
1127                     pNew[icount][0]  = pxNew;
1128                     pNew[icount][1]  = pyNew;
1129                     pNew[icount][2]  = plNew;
1130                     pNew[icount][3]  = eNew;
1131                     kNew[icount]     = i;
1132                     icount++;
1133                 } // parton loop
1134                 //
1135                 // Check if there was phase-space for quenching
1136                 //
1137
1138                 if (icount == 0) quenched[isys] = kFALSE;
1139                 if (!quenched[isys]) break;
1140                 
1141                 for (Int_t j = 0; j < 4; j++) 
1142                 {
1143                     p2[isys][j] = p0[isys][j] - p1[isys][j];
1144                 }
1145                 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];
1146                 if (p2[isys][4] > 0.) {
1147                     p2[isys][4] = TMath::Sqrt(p2[isys][4]);
1148                     break;
1149                 } else {
1150                     printf("Warning negative mass squared in system %d %f ! \n", isys, zquench[isys]);
1151                     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]);
1152                     if (p2[isys][4] < -0.01) {
1153                         printf("Negative mass squared !\n");
1154                         // Here we have to put the gluon back to mass shell
1155                         // This will lead to a small energy imbalance
1156                         p2[isys][4]  = 0.;
1157                         p2[isys][3]  = TMath::Sqrt(p2[isys][0] * p2[isys][0] + p2[isys][1] * p2[isys][1] + p2[isys][2] * p2[isys][2]);
1158                         break;
1159                     } else {
1160                         p2[isys][4] = 0.;
1161                         break;
1162                     }
1163                 }
1164                 /*
1165                 zHeavy *= 0.98;
1166                 printf("zHeavy lowered to %f\n", zHeavy);
1167                 if (zHeavy < 0.01) {
1168                     printf("No success ! \n");
1169                     icount = 0;
1170                     quenched[isys] = kFALSE;
1171                     break;
1172                 }
1173                 */
1174             } // iteration on z (while)
1175             
1176 //          Update  event record
1177             for (Int_t k = 0; k < icount; k++) {
1178 //              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] );
1179                 fPyjets->P[0][kNew[k]] = pNew[k][0];
1180                 fPyjets->P[1][kNew[k]] = pNew[k][1];
1181                 fPyjets->P[2][kNew[k]] = pNew[k][2];
1182                 fPyjets->P[3][kNew[k]] = pNew[k][3];
1183             }
1184             //
1185             // Add the gluons
1186             //
1187             Int_t ish = 0;    
1188             Int_t iGlu;
1189             if (!quenched[isys]) continue;
1190 //
1191 //      Last parton from shower i
1192             Int_t in = klast[isys];
1193 //
1194 //      Continue if no parton in shower i selected
1195             if (in == -1) continue;
1196 //  
1197 //      If this is the second initial parton and it is behind the first move pointer by previous ish
1198             if (isys == 1 && klast[1] > klast[0]) in += ish;
1199 //
1200 //      Starting index
1201             
1202 //          jmin = in - 1;
1203 // How many additional gluons will be generated
1204             ish  = 1;
1205             if (p2[isys][4] > 0.05) ish = 2;
1206 //
1207 //      Position of gluons
1208             iGlu = numpart;
1209             if (iglu == 0) igMin = iGlu;
1210             igMax = iGlu;
1211             numpart += ish;
1212             (fPyjets->N) += ish;
1213             
1214             if (ish == 1) {
1215                 fPyjets->P[0][iGlu] = p2[isys][0];
1216                 fPyjets->P[1][iGlu] = p2[isys][1];
1217                 fPyjets->P[2][iGlu] = p2[isys][2];
1218                 fPyjets->P[3][iGlu] = p2[isys][3];
1219                 fPyjets->P[4][iGlu] = p2[isys][4];
1220                 
1221                 fPyjets->K[0][iGlu] = 1;
1222                 if (iglu == nGluon[isys] - 1) fPyjets->K[0][iGlu] = 1;
1223                 fPyjets->K[1][iGlu] = 21;       
1224                 fPyjets->K[2][iGlu] = fPyjets->K[2][in] + 1000;
1225                 fPyjets->K[3][iGlu] = -1;       
1226                 fPyjets->K[4][iGlu] = -1;
1227                 
1228                 pg[0] += p2[isys][0];
1229                 pg[1] += p2[isys][1];
1230                 pg[2] += p2[isys][2];
1231                 pg[3] += p2[isys][3];
1232             } else {
1233                 //
1234                 // Split gluon in rest frame.
1235                 //
1236                 Double_t bx   =  p2[isys][0] / p2[isys][3];
1237                 Double_t by   =  p2[isys][1] / p2[isys][3];
1238                 Double_t bz   =  p2[isys][2] / p2[isys][3];
1239                 Double_t pst  =  p2[isys][4] / 2.;
1240                 //
1241                 // Isotropic decay ????
1242                 Double_t cost = 2. * gRandom->Rndm() - 1.;
1243                 Double_t sint = TMath::Sqrt((1.-cost)*(1.+cost));
1244                 Double_t phis =  2. * TMath::Pi() * gRandom->Rndm();
1245                 
1246                 Double_t pz1 =   pst * cost;
1247                 Double_t pz2 =  -pst * cost;
1248                 Double_t pt1 =   pst * sint;
1249                 Double_t pt2 =  -pst * sint;
1250                 Double_t px1 =   pt1 * TMath::Cos(phis);
1251                 Double_t py1 =   pt1 * TMath::Sin(phis);            
1252                 Double_t px2 =   pt2 * TMath::Cos(phis);
1253                 Double_t py2 =   pt2 * TMath::Sin(phis);            
1254                 
1255                 fPyjets->P[0][iGlu] = px1;
1256                 fPyjets->P[1][iGlu] = py1;
1257                 fPyjets->P[2][iGlu] = pz1;
1258                 fPyjets->P[3][iGlu] = pst;
1259                 fPyjets->P[4][iGlu] = 0.;
1260                 
1261                 fPyjets->K[0][iGlu] = 1 ;
1262                 fPyjets->K[1][iGlu] = 21;       
1263                 fPyjets->K[2][iGlu] = fPyjets->K[2][in] + 1000;
1264                 fPyjets->K[3][iGlu] = -1;       
1265                 fPyjets->K[4][iGlu] = -1;
1266                 
1267                 fPyjets->P[0][iGlu+1] = px2;
1268                 fPyjets->P[1][iGlu+1] = py2;
1269                 fPyjets->P[2][iGlu+1] = pz2;
1270                 fPyjets->P[3][iGlu+1] = pst;
1271                 fPyjets->P[4][iGlu+1] = 0.;
1272                 
1273                 fPyjets->K[0][iGlu+1] = 1;
1274                 if (iglu == nGluon[isys] - 1) fPyjets->K[0][iGlu+1] = 1;
1275                 fPyjets->K[1][iGlu+1] = 21;     
1276                 fPyjets->K[2][iGlu+1] = fPyjets->K[2][in] + 1000;
1277                 fPyjets->K[3][iGlu+1] = -1;     
1278                 fPyjets->K[4][iGlu+1] = -1;
1279                 SetMSTU(1,0);
1280                 SetMSTU(2,0);
1281                 //
1282                 // Boost back
1283                 //
1284                 Pyrobo(iGlu + 1, iGlu + 2, 0., 0., bx, by, bz);
1285             }
1286 /*    
1287             for (Int_t ig = iGlu; ig < iGlu+ish; ig++) {
1288                 Double_t px, py, pz;
1289                 px = fPyjets->P[0][ig]; 
1290                 py = fPyjets->P[1][ig]; 
1291                 pz = fPyjets->P[2][ig]; 
1292                 TVector3 v(px, py, pz);
1293                 v.RotateZ(-phiq[isys]);
1294                 v.RotateY(-thetaq[isys]);
1295                 Double_t pxs     = v.X(); Double_t pys = v.Y(); Double_t pzs  = v.Z();     
1296                 Double_t r       = AliPythiaRndm::GetPythiaRandom()->Rndm();
1297                 Double_t jtKick  = 0.3 * TMath::Sqrt(-TMath::Log(r));
1298                 if (ish == 2)   jtKick  = wjtKick[i] * TMath::Sqrt(-TMath::Log(r)) / TMath::Sqrt(2.);
1299                 Double_t phiKick = 2. * TMath::Pi() * AliPythiaRndm::GetPythiaRandom()->Rndm();
1300                 pxs += jtKick * TMath::Cos(phiKick);
1301                 pys += jtKick * TMath::Sin(phiKick);
1302                 TVector3 w(pxs, pys, pzs);
1303                 w.RotateY(thetaq[isys]);
1304                 w.RotateZ(phiq[isys]);
1305                 fPyjets->P[0][ig] = w.X(); 
1306                 fPyjets->P[1][ig] = w.Y(); 
1307                 fPyjets->P[2][ig] = w.Z(); 
1308                 fPyjets->P[2][ig] = w.Mag();
1309             }
1310 */
1311         } // kGluon         
1312         
1313         
1314     // Check energy conservation
1315         Double_t pxs = 0.;
1316         Double_t pys = 0.;
1317         Double_t pzs = 0.;      
1318         Double_t es  = 14000.;
1319         
1320         for (Int_t i = 0; i < numpart; i++)
1321         {
1322             kst =  fPyjets->K[0][i];
1323             if (kst != 1 && kst != 2) continue;
1324             pxs += fPyjets->P[0][i];
1325             pys += fPyjets->P[1][i];
1326             pzs += fPyjets->P[2][i];        
1327             es  -= fPyjets->P[3][i];        
1328         }
1329         if (TMath::Abs(pxs) > 1.e-2 ||
1330             TMath::Abs(pys) > 1.e-2 ||
1331             TMath::Abs(pzs) > 1.e-1) {
1332             printf("%e %e %e %e\n", pxs, pys, pzs, es);
1333 //              Fatal("Quench()", "4-Momentum non-conservation");
1334         }
1335         
1336     } // end quenching loop (systems)
1337 // Clean-up
1338     for (Int_t i = 0; i < numpart; i++)
1339     {
1340         imo =  fPyjets->K[2][i];
1341         if (imo > 1000) {
1342             fPyjets->K[2][i] = fPyjets->K[2][i] % 1000;
1343         }
1344     }
1345 //      this->Pylist(1);
1346 } // end quench
1347
1348
1349 void AliPythia6::Pyquen(Double_t a, Int_t ibf, Double_t b)
1350 {
1351     // Igor Lokthine's quenching routine
1352     // http://lokhtin.web.cern.ch/lokhtin/pyquen/pyquen.txt
1353
1354     pyquen(a, ibf, b);
1355 }
1356
1357 void AliPythia6::SetPyquenParameters(Double_t t0, Double_t tau0, Int_t nf, Int_t iengl, Int_t iangl)
1358 {
1359     // Set the parameters for the PYQUEN package.
1360     // See comments in PyquenCommon.h
1361     
1362     
1363     PYQPAR.t0    = t0;
1364     PYQPAR.tau0  = tau0;
1365     PYQPAR.nf    = nf;
1366     PYQPAR.iengl = iengl;
1367     PYQPAR.iangl = iangl;
1368 }
1369
1370 void  AliPythia6::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
1371 {
1372 //
1373 // Load event into Pythia Common Block
1374 //
1375
1376     Int_t npart = stack -> GetNprimary();
1377     Int_t n0 = 0;
1378     
1379     if (!flag) {
1380         GetPyjets()->N = npart;
1381     } else {
1382         n0 = GetPyjets()->N;
1383         GetPyjets()->N = n0 + npart;
1384     }
1385     
1386     
1387     for (Int_t part = 0; part < npart; part++) {
1388         TParticle *mPart = stack->Particle(part);
1389         
1390         Int_t kf     =  mPart->GetPdgCode();
1391         Int_t ks     =  mPart->GetStatusCode();
1392         Int_t idf    =  mPart->GetFirstDaughter();
1393         Int_t idl    =  mPart->GetLastDaughter();
1394         
1395         if (reHadr) {
1396             if (ks == 11 || ks == 12) {
1397                 ks  -= 10;
1398                 idf  = -1;
1399                 idl  = -1;
1400             }
1401         }
1402         
1403         Float_t px = mPart->Px();
1404         Float_t py = mPart->Py();
1405         Float_t pz = mPart->Pz();
1406         Float_t e  = mPart->Energy();
1407         Float_t m  = mPart->GetCalcMass();
1408         
1409         
1410         (GetPyjets())->P[0][part+n0] = px;
1411         (GetPyjets())->P[1][part+n0] = py;
1412         (GetPyjets())->P[2][part+n0] = pz;
1413         (GetPyjets())->P[3][part+n0] = e;
1414         (GetPyjets())->P[4][part+n0] = m;
1415         
1416         (GetPyjets())->K[1][part+n0] = kf;
1417         (GetPyjets())->K[0][part+n0] = ks;
1418         (GetPyjets())->K[3][part+n0] = idf + 1;
1419         (GetPyjets())->K[4][part+n0] = idl + 1;
1420         (GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1421     }
1422 }
1423
1424
1425 void AliPythia6::Pyevnw()
1426 {
1427     // New multiple interaction scenario
1428     pyevnw();
1429 }
1430
1431 void AliPythia6::GetQuenchingParameters(Double_t& xp, Double_t& yp, Double_t z[4])
1432 {
1433     // Return event specific quenching parameters
1434     xp = fXJet;
1435     yp = fYJet;
1436     for (Int_t i = 0; i < 4; i++) z[i] = fZQuench[i];
1437
1438 }
1439
1440 void AliPythia6::ConfigHeavyFlavor()
1441 {
1442     //
1443     // Default configuration for Heavy Flavor production
1444     //
1445     // All QCD processes
1446     //
1447     SetMSEL(1);
1448     
1449     // No multiple interactions
1450     SetMSTP(81,0);
1451     SetPARP(81, 0.);
1452     SetPARP(82, 0.);    
1453     // Initial/final parton shower on (Pythia default)
1454     SetMSTP(61,1);
1455     SetMSTP(71,1);
1456     
1457     // 2nd order alpha_s
1458     SetMSTP(2,2);
1459     
1460     // QCD scales
1461     SetMSTP(32,2);
1462     SetPARP(34,1.0);
1463 }
1464
1465 void AliPythia6::AtlasTuning()
1466 {
1467     //
1468     // Configuration for the ATLAS tuning
1469         SetMSTP(51,AliStructFuncType::PDFsetIndex(kCTEQ5L));
1470         SetMSTP(81,1);             // Multiple Interactions ON
1471         SetMSTP(82,4);             // Double Gaussian Model
1472         SetPARP(81,1.9);           // Min. pt for multiple interactions (default in 6.2-14) 
1473         SetPARP(82,1.8);           // [GeV]    PT_min at Ref. energy
1474         SetPARP(89,1000.);         // [GeV]   Ref. energy
1475         SetPARP(90,0.16);          // 2*epsilon (exponent in power law)
1476         SetPARP(83,0.5);           // Core density in proton matter distribution (def.value)
1477         SetPARP(84,0.5);           // Core radius
1478         SetPARP(85,0.33);          // Regulates gluon prod. mechanism
1479         SetPARP(86,0.66);          // Regulates gluon prod. mechanism
1480         SetPARP(67,1);             // Regulates Initial State Radiation
1481 }
1482
1483 void AliPythia6::AtlasTuningMC09()
1484 {
1485     //
1486     // Configuration for the ATLAS tuning
1487     printf("ATLAS New TUNE MC09\n");
1488     SetMSTP(81,21);             // treatment for MI, ISR, FSR and beam remnants: MI on, new model
1489     SetMSTP(82, 4);             // Double Gaussian Model
1490     SetMSTP(52, 2);             // External PDF
1491     SetMSTP(51, 20650);         // MRST LO*
1492   
1493     
1494     SetMSTP(70, 0);             // (was 2: def manual 1, def code 0) virtuality scale for ISR 
1495     SetMSTP(72, 1);             // (was 0: def 1) maximum scale for FSR
1496     SetMSTP(88, 1);             // (was 0: def 1) strategy for qq junction to di-quark or baryon in beam remnant
1497     SetMSTP(90, 0);             // (was 1: def 0) strategy of compensate the primordial kT
1498
1499     SetPARP(78, 0.3);           // the amount of color reconnection in the final state
1500     SetPARP(80, 0.1);           // probability of color partons kicked out from beam remnant
1501     SetPARP(82, 2.3);           // [GeV]    PT_min at Ref. energy    
1502     SetPARP(83, 0.8);           // Core density in proton matter distribution (def.value)    
1503     SetPARP(84, 0.7);           // Core radius
1504     SetPARP(90, 0.25);          //  2*epsilon (exponent in power law)
1505     SetPARJ(81, 0.29);          // (was 0.14: def 0.29) Labmda value in running alpha_s for parton showers
1506
1507     SetMSTP(95, 6);
1508     SetPARJ(41, 0.3);           // a and b parameters of the symmm. Lund FF
1509     SetPARJ(42, 0.58);
1510     SetPARJ(46, 0.75);          // mod. of the Lund FF for heavy end-point quarks
1511     SetPARP(89,1800.);         // [GeV]   Ref. energy
1512 }
1513
1514 void AliPythia6::SetWeightPower(Double_t pow)
1515 {
1516     setpowwght(pow);
1517     SetMSTP(142, 1); // Tell Pythia to use pyevwt to calculate event wghts
1518 }
1519
1520 void AliPythia6::SetPtHardRange(Float_t ptmin, Float_t ptmax)
1521 {
1522     // Set the pt hard range
1523     SetCKIN(3, ptmin);
1524     SetCKIN(4, ptmax);
1525 }
1526
1527 void AliPythia6::SetYHardRange(Float_t ymin, Float_t ymax)
1528 {
1529     // Set the y hard range
1530     SetCKIN(7, ymin);
1531     SetCKIN(8, ymax);
1532 }
1533
1534
1535 void AliPythia6::SetFragmentation(Int_t flag)
1536 {
1537     // Switch fragmentation on/off
1538     SetMSTP(111, flag);
1539 }
1540
1541 void AliPythia6::SetInitialAndFinalStateRadiation(Int_t flag1, Int_t flag2)
1542 {
1543 //  initial state radiation    
1544     SetMSTP(61, flag1);
1545 //  final state radiation
1546     SetMSTP(71, flag2);
1547 }
1548
1549 void AliPythia6::SetIntrinsicKt(Float_t kt)
1550 {
1551     // Set the inreinsic kt
1552     if (kt > 0.) {
1553         SetMSTP(91,1);
1554         SetPARP(91,kt); 
1555         SetPARP(93, 4. * kt);
1556     } else {
1557         SetMSTP(91,0);
1558     }
1559 }
1560
1561 void AliPythia6::SwitchHFOff()
1562 {
1563     // Switch off heavy flavor
1564     // Maximum number of quark flavours used in pdf 
1565     SetMSTP(58, 3);
1566     // Maximum number of flavors that can be used in showers
1567     SetMSTJ(45, 3);     
1568 }
1569
1570 void AliPythia6::SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
1571                                        Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
1572 {
1573 // Set pycell parameters
1574     SetPARU(51, etamax);
1575     SetMSTU(51, neta);
1576     SetMSTU(52, nphi);
1577     SetPARU(58, thresh);
1578     SetPARU(52, etseed);
1579     SetPARU(53, minet);
1580     SetPARU(54, r);
1581     SetMSTU(54,  2);
1582 }
1583
1584 void AliPythia6::ModifiedSplitting()
1585 {
1586     // Modified splitting probability as a model for quenching
1587     SetPARJ(200, 0.8);
1588     SetMSTJ(41, 1);  // QCD radiation only
1589     SetMSTJ(42, 2);  // angular ordering
1590     SetMSTJ(44, 2);  // option to run alpha_s
1591     SetMSTJ(47, 0);  // No correction back to hard scattering element
1592     SetMSTJ(50, 0);  // No coherence in first branching
1593     SetPARJ(82, 1.); // Cut off for parton showers
1594 }
1595
1596 void AliPythia6::SwitchHadronisationOff()
1597 {
1598     // Switch off hadronisarion
1599     SetMSTJ(1, 0);
1600 }
1601
1602 void AliPythia6::SwitchHadronisationOn()
1603 {
1604     // Switch on hadronisarion
1605     SetMSTJ(1, 1);
1606 }
1607
1608
1609 void AliPythia6::GetXandQ(Float_t& x1, Float_t& x2, Float_t& q)
1610 {
1611     // Get x1, x2 and Q for this event
1612     
1613     q  = GetVINT(51);
1614     x1 = GetVINT(41);
1615     x2 = GetVINT(42);
1616 }
1617
1618 Float_t AliPythia6::GetXSection()
1619 {
1620     // Get the total cross-section
1621     return (GetPARI(1));
1622 }
1623
1624 Float_t AliPythia6::GetPtHard()
1625 {
1626     // Get the pT hard for this event
1627     return GetVINT(47);
1628 }
1629
1630 Int_t AliPythia6::ProcessCode()
1631 {
1632     // Get the subprocess code
1633     return GetMSTI(1);
1634 }
1635
1636 void AliPythia6::PrintStatistics()
1637 {
1638     // End of run statistics
1639     Pystat(1);
1640 }
1641
1642 void AliPythia6::EventListing()
1643 {
1644     // End of run statistics
1645     Pylist(2);
1646 }
1647
1648 AliPythia6& AliPythia6::operator=(const  AliPythia6& rhs)
1649 {
1650 // Assignment operator
1651     rhs.Copy(*this);
1652     return *this;
1653 }
1654
1655  void AliPythia6::Copy(TObject&) const
1656 {
1657     //
1658     // Copy 
1659     //
1660     Fatal("Copy","Not implemented!\n");
1661 }