Line data Source code
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 2 : 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 0 : TPythia6(),
74 0 : AliPythiaBase(),
75 0 : fProcess(kPyMb),
76 0 : fEcms(0.),
77 0 : fStrucFunc(kCTEQ5L),
78 0 : fProjectile("p"),
79 0 : fTarget("p"),
80 0 : fXJet(0.),
81 0 : fYJet(0.),
82 0 : fNGmax(30),
83 0 : fZmax(0.97),
84 0 : fGlauber(0),
85 0 : fQuenchingWeights(0)
86 0 : {
87 : // Default Constructor
88 : //
89 : // Set random number
90 : Int_t i;
91 0 : for (i = 0; i < 501; i++) fDefMDCY[i] = 0;
92 0 : for (i = 0; i < 2001; i++) fDefMDME[i] = 0;
93 0 : for (i = 0; i < 4; i++) fZQuench[i] = 0;
94 :
95 0 : if (!AliPythiaRndm::GetPythiaRandom())
96 0 : AliPythiaRndm::SetPythiaRandom(GetRandom());
97 0 : fGlauber = 0;
98 0 : fQuenchingWeights = 0;
99 0 : }
100 :
101 : AliPythia6::AliPythia6(const AliPythia6& pythia):
102 0 : TPythia6(),
103 0 : AliPythiaBase(),
104 0 : fProcess(kPyMb),
105 0 : fEcms(0.),
106 0 : fStrucFunc(kCTEQ5L),
107 0 : fProjectile("p"),
108 0 : fTarget("p"),
109 0 : fXJet(0.),
110 0 : fYJet(0.),
111 0 : fNGmax(30),
112 0 : fZmax(0.97),
113 0 : fGlauber(0),
114 0 : fQuenchingWeights(0)
115 0 : {
116 : // Copy Constructor
117 : Int_t i;
118 0 : for (i = 0; i < 501; i++) fDefMDCY[i] = 0;
119 0 : for (i = 0; i < 2001; i++) fDefMDME[i] = 0;
120 0 : for (i = 0; i < 4; i++) fZQuench[i] = 0;
121 0 : pythia.Copy(*this);
122 0 : }
123 :
124 0 : void AliPythia6::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfunc, Int_t /*tune*/)
125 : {
126 : // Initialise the process to generate
127 0 : if (!AliPythiaRndm::GetPythiaRandom())
128 0 : AliPythiaRndm::SetPythiaRandom(GetRandom());
129 :
130 0 : fProcess = process;
131 0 : fEcms = energy;
132 0 : fStrucFunc = strucfunc;
133 : //...Switch off decay of pi0, K0S, Lambda, Sigma+-, Xi0-, Omega-.
134 0 : SetMDCY(Pycomp(111) ,1,0);
135 0 : SetMDCY(Pycomp(310) ,1,0);
136 0 : SetMDCY(Pycomp(3122),1,0);
137 0 : SetMDCY(Pycomp(3112),1,0);
138 0 : SetMDCY(Pycomp(3212),1,0);
139 0 : SetMDCY(Pycomp(3222),1,0);
140 0 : SetMDCY(Pycomp(3312),1,0);
141 0 : SetMDCY(Pycomp(3322),1,0);
142 0 : SetMDCY(Pycomp(3334),1,0);
143 : // Select structure function
144 0 : SetMSTP(52,2);
145 0 : 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 0 : SetMSTU(16,2);
149 :
150 : //
151 : // Pythia initialisation for selected processes//
152 : //
153 : // Make MSEL clean
154 : //
155 0 : for (Int_t i=1; i<= 200; i++) {
156 0 : SetMSUB(i,0);
157 : }
158 : // select charm production
159 0 : switch (process)
160 : {
161 : case kPyOldUEQ2ordered: //Old underlying events with Q2 ordered QCD processes
162 : // Multiple interactions on.
163 0 : SetMSTP(81,1);
164 : // Double Gaussian matter distribution.
165 0 : SetMSTP(82,4);
166 0 : SetPARP(83,0.5);
167 0 : SetPARP(84,0.4);
168 : // pT0.
169 0 : SetPARP(82,2.0);
170 : // Reference energy for pT0 and energy rescaling pace.
171 0 : SetPARP(89,1800);
172 0 : SetPARP(90,0.25);
173 : // String drawing almost completely minimizes string length.
174 0 : SetPARP(85,0.9);
175 0 : SetPARP(86,0.95);
176 : // ISR and FSR activity.
177 0 : SetPARP(67,4);
178 0 : SetPARP(71,4);
179 : // Lambda_FSR scale.
180 0 : SetPARJ(81,0.29);
181 0 : break;
182 : case kPyOldUEQ2ordered2:
183 : // Old underlying events with Q2 ordered QCD processes
184 : // Multiple interactions on.
185 0 : SetMSTP(81,1);
186 : // Double Gaussian matter distribution.
187 0 : SetMSTP(82,4);
188 0 : SetPARP(83,0.5);
189 0 : SetPARP(84,0.4);
190 : // pT0.
191 0 : SetPARP(82,2.0);
192 : // Reference energy for pT0 and energy rescaling pace.
193 0 : SetPARP(89,1800);
194 0 : SetPARP(90,0.16); // here is the difference with kPyOldUEQ2ordered
195 : // String drawing almost completely minimizes string length.
196 0 : SetPARP(85,0.9);
197 0 : SetPARP(86,0.95);
198 : // ISR and FSR activity.
199 0 : SetPARP(67,4);
200 0 : SetPARP(71,4);
201 : // Lambda_FSR scale.
202 0 : SetPARJ(81,0.29);
203 0 : break;
204 : case kPyOldPopcorn:
205 : // Old production mechanism: Old Popcorn
206 0 : SetMSEL(1);
207 0 : SetMSTJ(12,3);
208 : // (D=2) Like MSTJ(12)=2 but added prod ofthe 1er rank baryon
209 0 : SetMSTP(88,2);
210 : // (D=1)see can be used to form baryons (BARYON JUNCTION)
211 0 : SetMSTJ(1,1);
212 0 : AtlasTuning();
213 0 : break;
214 : case kPyCharm:
215 0 : SetMSEL(4);
216 : // heavy quark masses
217 :
218 0 : SetPMAS(4,1,1.2);
219 : //
220 : // primordial pT
221 0 : SetMSTP(91,1);
222 0 : SetPARP(91,1.);
223 0 : SetPARP(93,5.);
224 : //
225 0 : break;
226 : case kPyBeauty:
227 0 : SetMSEL(5);
228 0 : SetPMAS(5,1,4.75);
229 0 : break;
230 : case kPyJpsi:
231 0 : SetMSEL(0);
232 : // gg->J/Psi g
233 0 : SetMSUB(86,1);
234 0 : break;
235 : case kPyJpsiChi:
236 0 : SetMSEL(0);
237 : // gg->J/Psi g
238 0 : SetMSUB(86,1);
239 : // gg-> chi_0c g
240 0 : SetMSUB(87,1);
241 : // gg-> chi_1c g
242 0 : SetMSUB(88,1);
243 : // gg-> chi_2c g
244 0 : SetMSUB(89,1);
245 0 : break;
246 : case kPyCharmUnforced:
247 0 : SetMSEL(0);
248 : // gq->qg
249 0 : SetMSUB(28,1);
250 : // gg->qq
251 0 : SetMSUB(53,1);
252 : // gg->gg
253 0 : SetMSUB(68,1);
254 0 : break;
255 : case kPyBeautyUnforced:
256 0 : SetMSEL(0);
257 : // gq->qg
258 0 : SetMSUB(28,1);
259 : // gg->qq
260 0 : SetMSUB(53,1);
261 : // gg->gg
262 0 : SetMSUB(68,1);
263 0 : break;
264 : case kPyMb:
265 : // Minimum Bias pp-Collisions
266 : //
267 : //
268 : // select Pythia min. bias model
269 0 : SetMSEL(0);
270 0 : SetMSUB(92,1); // single diffraction AB-->XB
271 0 : SetMSUB(93,1); // single diffraction AB-->AX
272 0 : SetMSUB(94,1); // double diffraction
273 0 : SetMSUB(95,1); // low pt production
274 :
275 0 : AtlasTuning();
276 0 : break;
277 : case kPyMbAtlasTuneMC09:
278 : // Minimum Bias pp-Collisions
279 : //
280 : //
281 : // select Pythia min. bias model
282 0 : SetMSEL(0);
283 0 : SetMSUB(92,1); // single diffraction AB-->XB
284 0 : SetMSUB(93,1); // single diffraction AB-->AX
285 0 : SetMSUB(94,1); // double diffraction
286 0 : SetMSUB(95,1); // low pt production
287 :
288 0 : AtlasTuningMC09();
289 0 : break;
290 :
291 : case kPyMbWithDirectPhoton:
292 : // Minimum Bias pp-Collisions with direct photon processes added
293 : //
294 : //
295 : // select Pythia min. bias model
296 0 : SetMSEL(0);
297 0 : SetMSUB(92,1); // single diffraction AB-->XB
298 0 : SetMSUB(93,1); // single diffraction AB-->AX
299 0 : SetMSUB(94,1); // double diffraction
300 0 : SetMSUB(95,1); // low pt production
301 :
302 0 : SetMSUB(14,1); //
303 0 : SetMSUB(18,1); //
304 0 : SetMSUB(29,1); //
305 0 : SetMSUB(114,1); //
306 0 : SetMSUB(115,1); //
307 :
308 :
309 0 : AtlasTuning();
310 0 : break;
311 :
312 : case kPyMbDefault:
313 : // Minimum Bias pp-Collisions
314 : //
315 : //
316 : // select Pythia min. bias model
317 0 : SetMSEL(0);
318 0 : SetMSUB(92,1); // single diffraction AB-->XB
319 0 : SetMSUB(93,1); // single diffraction AB-->AX
320 0 : SetMSUB(94,1); // double diffraction
321 0 : SetMSUB(95,1); // low pt production
322 :
323 0 : 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 0 : SetMSEL(0);
329 0 : SetMSUB(92,1); // single diffraction AB-->XB
330 0 : SetMSUB(93,1); // single diffraction AB-->AX
331 0 : SetMSUB(94,1); // double diffraction
332 0 : SetMSUB(95,1); // low pt production
333 0 : SetMSTP(51,AliStructFuncType::PDFsetIndex(kCTEQ6ll));
334 0 : SetMSTP(52,2);
335 0 : SetMSTP(68,1);
336 0 : SetMSTP(70,2);
337 0 : SetMSTP(81,1); // Multiple Interactions ON
338 0 : SetMSTP(82,4); // Double Gaussian Model
339 0 : SetMSTP(88,1);
340 :
341 0 : SetPARP(82,2.3); // [GeV] PT_min at Ref. energy
342 0 : SetPARP(83,0.5); // Core density in proton matter distribution (def.value)
343 0 : SetPARP(84,0.5); // Core radius
344 0 : SetPARP(85,0.9); // Regulates gluon prod. mechanism
345 0 : SetPARP(90,0.2); // 2*epsilon (exponent in power law)
346 :
347 0 : break;
348 : case kPyMbNonDiffr:
349 : // Minimum Bias pp-Collisions
350 : //
351 : //
352 : // select Pythia min. bias model
353 0 : SetMSEL(0);
354 0 : SetMSUB(95,1); // low pt production
355 :
356 0 : AtlasTuning();
357 0 : break;
358 : case kPyMbMSEL1:
359 0 : ConfigHeavyFlavor();
360 : // Intrinsic <kT^2>
361 0 : SetMSTP(91,1);// Width (1=gaussian) primordial kT dist. inside hadrons
362 0 : SetPARP(91,1.); // <kT^2> = PARP(91,1.)^2
363 0 : SetPARP(93,5.); // Upper cut-off
364 : // Set Q-quark mass
365 0 : SetPMAS(4,1,1.2); // Charm quark mass
366 0 : SetPMAS(5,1,4.78); // Beauty quark mass
367 0 : SetPARP(71,4.); // Defaut value
368 : // Atlas Tuning
369 0 : AtlasTuning();
370 0 : break;
371 : case kPyJets:
372 : //
373 : // QCD Jets
374 : //
375 0 : SetMSEL(1);
376 : // Pythia Tune A (CDF)
377 : //
378 0 : SetPARP(67,2.5); // Regulates Initial State Radiation (value from best fit to D0 dijet analysis)
379 0 : SetMSTP(82,4); // Double Gaussian Model
380 0 : SetPARP(82,2.0); // [GeV] PT_min at Ref. energy
381 0 : SetPARP(84,0.4); // Core radius
382 0 : SetPARP(85,0.90) ; // Regulates gluon prod. mechanism
383 0 : SetPARP(86,0.95); // Regulates gluon prod. mechanism
384 0 : SetPARP(89,1800.); // [GeV] Ref. energy
385 0 : SetPARP(90,0.25); // 2*epsilon (exponent in power law)
386 0 : break;
387 : case kPyDirectGamma:
388 0 : SetMSEL(10);
389 0 : 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 0 : ConfigHeavyFlavor();
402 : // Intrinsic <kT>
403 0 : SetMSTP(91,1);
404 0 : SetPARP(91,1.304);
405 0 : SetPARP(93,6.52);
406 : // Set c-quark mass
407 0 : SetPMAS(4,1,1.2);
408 0 : 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 0 : ConfigHeavyFlavor();
421 : // Intrinsic <kT>
422 0 : SetMSTP(91,1);
423 0 : SetPARP(91,1.16);
424 0 : SetPARP(93,5.8);
425 :
426 : // Set c-quark mass
427 0 : SetPMAS(4,1,1.2);
428 0 : 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 0 : ConfigHeavyFlavor();
442 : // Intrinsic <kT^2>
443 0 : SetMSTP(91,1);
444 0 : SetPARP(91,1.);
445 0 : SetPARP(93,5.);
446 :
447 : // Set c-quark mass
448 0 : SetPMAS(4,1,1.2);
449 0 : 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 0 : ConfigHeavyFlavor();
468 : // Intrinsic <kT^2>
469 0 : SetMSTP(91,1);
470 0 : SetPARP(91,1.);
471 0 : SetPARP(93,5.);
472 :
473 : // Set c-quark mass
474 0 : SetPMAS(4,1,1.2);
475 0 : AtlasTuning();
476 0 : 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 0 : ConfigHeavyFlavor();
486 : // QCD scales
487 0 : SetPARP(67,1.0);
488 0 : SetPARP(71,1.0);
489 : // Intrinsic <kT>
490 0 : SetMSTP(91,1);
491 0 : SetPARP(91,2.035);
492 0 : SetPARP(93,10.17);
493 : // Set b-quark mass
494 0 : SetPMAS(5,1,4.75);
495 0 : 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 0 : ConfigHeavyFlavor();
505 : // QCD scales
506 0 : SetPARP(67,1.0);
507 0 : SetPARP(71,1.0);
508 : // Intrinsic <kT>
509 0 : SetMSTP(91,1);
510 0 : SetPARP(91,1.60);
511 0 : SetPARP(93,8.00);
512 : // Set b-quark mass
513 0 : SetPMAS(5,1,4.75);
514 0 : 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 0 : ConfigHeavyFlavor();
524 : // QCD scales
525 0 : SetPARP(67,1.0);
526 0 : SetPARP(71,1.0);
527 :
528 : // Intrinsic <kT>
529 0 : SetMSTP(91,1);
530 0 : SetPARP(91,1.);
531 0 : SetPARP(93,5.);
532 :
533 : // Set b-quark mass
534 0 : SetPMAS(5,1,4.75);
535 0 : 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 0 : ConfigHeavyFlavor();
555 : // QCD scales
556 0 : SetPARP(67,1.0);
557 0 : SetPARP(71,1.0);
558 :
559 : // Intrinsic <kT>
560 0 : SetMSTP(91,1);
561 0 : SetPARP(91,1.);
562 0 : SetPARP(93,5.);
563 :
564 : // Set b-quark mass
565 0 : SetPMAS(5,1,4.75);
566 :
567 0 : AtlasTuning();
568 0 : break;
569 : case kPyW:
570 :
571 : //Inclusive production of W+/-
572 0 : SetMSEL(0);
573 : //f fbar -> W+
574 0 : 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 0 : SetMSTP(61,1); //Initial QCD & QED showers on
587 0 : SetMSTP(71,1); //Final QCD & QED showers on
588 :
589 0 : break;
590 :
591 : case kPyZ:
592 :
593 : //Inclusive production of Z
594 0 : SetMSEL(0);
595 : //f fbar -> Z/gamma
596 0 : 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 0 : SetMSTP(43,2);
609 :
610 : // Initial/final parton shower on (Pythia default)
611 : // With parton showers on we are generating "Z inclusive process"
612 0 : SetMSTP(61,1); //Initial QCD & QED showers on
613 0 : SetMSTP(71,1); //Final QCD & QED showers on
614 0 : break;
615 : case kPyZgamma:
616 : //Inclusive production of Z
617 0 : SetMSEL(0);
618 : //f fbar -> Z/gamma
619 0 : SetMSUB(1,1);
620 : // Initial/final parton shower on (Pythia default)
621 : // With parton showers on we are generating "Z inclusive process"
622 0 : SetMSTP(61,1); //Initial QCD & QED showers on
623 0 : SetMSTP(71,1); //Final QCD & QED showers on
624 0 : 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 0 : SetMSTP(86,1);
638 : // maximum number of errors before pythia aborts (def=10)
639 0 : SetMSTU(22,10);
640 : // number of warnings printed on the shell
641 0 : SetMSTU(26,20);
642 0 : break;
643 : case kPyCharmPWHG:
644 : case kPyBeautyPWHG:
645 : case kPyWPWHG:
646 : // number of warnings printed on the shell
647 0 : SetMSTU(26,20);
648 :
649 0 : }
650 : //
651 : // Initialize PYTHIA
652 0 : SetMSTP(41,1); // all resonance decays switched on
653 0 : if (process == kPyJetsPWHG || process == kPyCharmPWHG || process == kPyBeautyPWHG || process == kPyWPWHG) {
654 0 : Initialize("USER","","",0.);
655 0 : } else {
656 0 : Initialize("CMS",fProjectile,fTarget,fEcms);
657 : }
658 0 : }
659 :
660 0 : Int_t AliPythia6::CheckedLuComp(Int_t kf)
661 : {
662 : // Check Lund particle code (for debugging)
663 0 : Int_t kc=Pycomp(kf);
664 0 : return kc;
665 : }
666 :
667 0 : 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 0 : SetMSTP(52,2);
679 0 : SetMSTP(192, a1);
680 0 : SetMSTP(193, a2);
681 0 : }
682 :
683 :
684 : AliPythia6* AliPythia6::Instance()
685 : {
686 : // Set random number generator
687 0 : if (fgAliPythia) {
688 0 : return fgAliPythia;
689 : } else {
690 0 : fgAliPythia = new AliPythia6();
691 0 : return fgAliPythia;
692 : }
693 0 : }
694 :
695 0 : void AliPythia6::PrintParticles()
696 : {
697 : // Print list of particl properties
698 : Int_t np = 0;
699 0 : char* name = new char[16];
700 0 : for (Int_t kf=0; kf<1000000; kf++) {
701 0 : for (Int_t c = 1; c > -2; c-=2) {
702 0 : Int_t kc = Pycomp(c*kf);
703 0 : if (kc) {
704 0 : Float_t mass = GetPMAS(kc,1);
705 0 : Float_t width = GetPMAS(kc,2);
706 0 : Float_t tau = GetPMAS(kc,4);
707 :
708 0 : Pyname(kf,name);
709 :
710 0 : np++;
711 :
712 0 : printf("\n mass, width, tau: %6d %s %10.3f %10.3e %10.3e",
713 0 : c*kf, name, mass, width, tau);
714 0 : }
715 : }
716 : }
717 0 : printf("\n Number of particles %d \n \n", np);
718 0 : }
719 :
720 0 : void AliPythia6::ResetDecayTable()
721 : {
722 : // Set default values for pythia decay switches
723 : Int_t i;
724 0 : for (i = 1; i < 501; i++) SetMDCY(i,1,fDefMDCY[i]);
725 0 : for (i = 1; i < 2001; i++) SetMDME(i,1,fDefMDME[i]);
726 0 : }
727 :
728 0 : void AliPythia6::SetDecayTable()
729 : {
730 : // Set default values for pythia decay switches
731 : //
732 : Int_t i;
733 0 : for (i = 1; i < 501; i++) fDefMDCY[i] = GetMDCY(i,1);
734 0 : for (i = 1; i < 2001; i++) fDefMDME[i] = GetMDME(i,1);
735 0 : }
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 0 : pyjoin(npart, ipart);
743 0 : }
744 :
745 : void AliPythia6::Pyshowq(Int_t ip1, Int_t ip2, Double_t qmax)
746 : {
747 : // Call qPythia showering
748 : //
749 0 : pyshowq(ip1, ip2, qmax);
750 0 : }
751 :
752 : void AliPythia6::Qpygin0()
753 : {
754 : //position of the hard scattering in the nuclear overlapping area.
755 : //just for qpythia.
756 0 : qpygin0();
757 0 : }
758 :
759 0 : void AliPythia6::Pyclus(Int_t& njet)
760 : {
761 : // Call Pythia clustering algorithm
762 : //
763 0 : pyclus(njet);
764 0 : }
765 :
766 0 : void AliPythia6::Pycell(Int_t& njet)
767 : {
768 : // Call Pythia jet reconstruction algorithm
769 : //
770 0 : pycell(njet);
771 0 : }
772 :
773 0 : 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 0 : Int_t n = GetN();
777 0 : px = GetPyjets()->P[0][n+i];
778 0 : py = GetPyjets()->P[1][n+i];
779 0 : pz = GetPyjets()->P[2][n+i];
780 0 : e = GetPyjets()->P[3][n+i];
781 0 : }
782 :
783 0 : void AliPythia6::Pyshow(Int_t ip1, Int_t ip2, Double_t qmax)
784 : {
785 : // Call Pythia showering
786 : //
787 0 : pyshow(ip1, ip2, qmax);
788 0 : }
789 :
790 0 : 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 0 : pyrobo(imi, ima, the, phi, bex, bey, bez);
793 0 : }
794 :
795 :
796 :
797 0 : 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 0 : fGlauber = AliFastGlauber::Instance();
805 0 : fGlauber->Init(2);
806 0 : fGlauber->SetCentralityClass(cMin, cMax);
807 :
808 0 : fQuenchingWeights = new AliQuenchingWeights();
809 0 : fQuenchingWeights->InitMult();
810 0 : fQuenchingWeights->SetK(k);
811 0 : fQuenchingWeights->SetECMethod(AliQuenchingWeights::kECMethod(iECMethod));
812 0 : fNGmax = ngmax;
813 0 : fZmax = zmax;
814 :
815 0 : }
816 :
817 :
818 0 : 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 0 : Double_t p0[4][5];
841 0 : Double_t p1[4][5];
842 0 : Double_t p2[4][5];
843 0 : Int_t klast[4] = {-1, -1, -1, -1};
844 :
845 0 : 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 0 : Double_t pxq[4], pyq[4], pzq[4], eq[4], yq[4], mq[4], pq[4], phiq[4], thetaq[4], ptq[4];
848 0 : Bool_t quenched[4];
849 0 : Double_t wjtKick[4] = {0., 0., 0., 0.};
850 0 : Int_t nGluon[4];
851 0 : 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 0 : for (Int_t i = 2; i <= 7; i++) {
862 : Int_t j = 0;
863 : // Skip gluons that participate in hard scattering
864 0 : if (i == 4 || i == 5) continue;
865 : // Gluons from hard Scattering
866 0 : if (i == 6 || i == 7) {
867 0 : j = i - 6;
868 0 : pxq[j] = fPyjets->P[0][i];
869 0 : pyq[j] = fPyjets->P[1][i];
870 0 : pzq[j] = fPyjets->P[2][i];
871 0 : eq[j] = fPyjets->P[3][i];
872 0 : mq[j] = fPyjets->P[4][i];
873 0 : } 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 0 : pxq[j] = fPyjets->P[0][i] - fPyjets->P[0][i+2];
885 0 : pyq[j] = fPyjets->P[1][i] - fPyjets->P[1][i+2];
886 0 : pzq[j] = fPyjets->P[2][i] - fPyjets->P[2][i+2];
887 0 : mq[j] = fPyjets->P[4][i];
888 0 : 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 0 : yq[j] = 0.5 * TMath::Log((eq[j] + pzq[j] + 1.e-14) / (eq[j] - pzq[j] + 1.e-14));
894 0 : pq[j] = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j] + pzq[j] * pzq[j]);
895 0 : phiq[j] = TMath::Pi()+TMath::ATan2(-pyq[j], -pxq[j]);
896 0 : ptq[j] = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j]);
897 0 : thetaq[j] = TMath::ATan2(ptq[j], pzq[j]);
898 0 : qPdg[j] = fPyjets->K[1][i];
899 0 : }
900 :
901 0 : Double_t int0[4];
902 0 : Double_t int1[4];
903 :
904 0 : fGlauber->GetI0I1ForPythiaAndXY(4, phiq, int0, int1, fXJet, fYJet, 15.);
905 :
906 0 : for (Int_t j = 0; j < 4; j++) {
907 : //
908 : // Quench only central jets and with E > 10.
909 : //
910 :
911 :
912 0 : Int_t itype = (qPdg[j] == 21) ? 2 : 1;
913 : // Double_t eloss = fQuenchingWeights->GetELossRandomKFast(itype, int0[j], int1[j], eq[j]);
914 0 : Double_t eloss = fQuenchingWeights->GetELossRandomK(itype, int0[j], int1[j], eq[j]);
915 :
916 0 : if (TMath::Abs(yq[j]) > 2.5 || eq[j] < 10.) {
917 0 : fZQuench[j] = 0.;
918 0 : } else {
919 0 : if (eq[j] > 40. && TMath::Abs(yq[j]) < 0.5) {
920 0 : icall ++;
921 0 : eMean += eloss;
922 0 : }
923 : //
924 : // Extra pt
925 0 : Double_t l = fQuenchingWeights->CalcLk(int0[j], int1[j]);
926 0 : wjtKick[j] = TMath::Sqrt(l * fQuenchingWeights->CalcQk(int0[j], int1[j]));
927 : //
928 : // Fractional energy loss
929 0 : fZQuench[j] = eloss / eq[j];
930 : //
931 : // Avoid complete loss
932 : //
933 0 : 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 0 : quenched[j] = (fZQuench[j] > 0.01);
946 : } // primary partons
947 :
948 :
949 :
950 0 : Double_t pNew[1000][4];
951 0 : Int_t kNew[1000];
952 : Int_t icount = 0;
953 0 : Double_t zquench[4];
954 :
955 : //
956 : // System Loop
957 0 : for (Int_t isys = 0; isys < 4; isys++) {
958 : // Skip to next system if not quenched.
959 0 : if (!quenched[isys]) continue;
960 :
961 0 : nGluon[isys] = 1 + Int_t(fZQuench[isys] / (1. - fZQuench[isys]));
962 0 : if (nGluon[isys] > fNGmax) nGluon[isys] = fNGmax;
963 0 : zquench[isys] = 1. - TMath::Power(1. - fZQuench[isys], 1./Double_t(nGluon[isys]));
964 0 : 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 0 : for (Int_t iglu = 0; iglu < nGluon[isys]; iglu++) {
976 : while (1) {
977 : icount = 0;
978 0 : for (Int_t k = 0; k < 4; k++)
979 : {
980 0 : p0[isys][k] = 0.;
981 0 : p1[isys][k] = 0.;
982 0 : p2[isys][k] = 0.;
983 : }
984 : // Loop over partons
985 0 : for (Int_t i = 0; i < numpart; i++)
986 : {
987 0 : imo = fPyjets->K[2][i];
988 0 : kst = fPyjets->K[0][i];
989 0 : pdg = fPyjets->K[1][i];
990 :
991 :
992 :
993 : // Quarks and gluons only
994 0 : if (pdg != 21 && TMath::Abs(pdg) > 6) continue;
995 : // Particles from hard scattering only
996 :
997 0 : if (imo > 8 && imo < 1000) imo = fPyjets->K[2][imo - 1];
998 0 : Int_t imom = imo % 1000;
999 0 : if ((isys == 0 || isys == 1) && ((imom != (isys + 7)))) continue;
1000 0 : if ((isys == 2 || isys == 3) && ((imom != (isys + 1)))) continue;
1001 :
1002 :
1003 : // Skip comment lines
1004 0 : if (kst != 1 && kst != 2) continue;
1005 : //
1006 : // Parton kinematic
1007 0 : px = fPyjets->P[0][i];
1008 0 : py = fPyjets->P[1][i];
1009 0 : pz = fPyjets->P[2][i];
1010 0 : e = fPyjets->P[3][i];
1011 0 : m = fPyjets->P[4][i];
1012 0 : pt = TMath::Sqrt(px * px + py * py);
1013 0 : p = TMath::Sqrt(px * px + py * py + pz * pz);
1014 0 : phi = TMath::Pi() + TMath::ATan2(-py, -px);
1015 0 : theta = TMath::ATan2(pt, pz);
1016 :
1017 : //
1018 : // Save 4-momentum sum for balancing
1019 : Int_t index = isys;
1020 :
1021 0 : p0[index][0] += px;
1022 0 : p0[index][1] += py;
1023 0 : p0[index][2] += pz;
1024 0 : p0[index][3] += e;
1025 :
1026 0 : klast[index] = i;
1027 :
1028 : //
1029 : // Fractional energy loss
1030 0 : Double_t z = zquench[index];
1031 :
1032 :
1033 : // Don't fully quench radiated gluons
1034 : //
1035 0 : 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 0 : TVector3 v(px, py, pz);
1051 0 : v.RotateZ(-phiq[index]); v.RotateY(-thetaq[index]);
1052 0 : Double_t pxs = v.X(); Double_t pys = v.Y(); Double_t pl = v.Z();
1053 :
1054 0 : Double_t jt = TMath::Sqrt(pxs * pxs + pys * pys);
1055 0 : Double_t mt2 = jt * jt + m * m;
1056 : Double_t zmax = 1.;
1057 : //
1058 : // Kinematic limit on z
1059 : //
1060 0 : if (m > 0.) zmax = 1. - m / TMath::Sqrt(m * m + jt * jt);
1061 : //
1062 : // Change light-cone kinematics rel. to initial parton
1063 : //
1064 0 : Double_t eppzOld = e + pl;
1065 0 : Double_t empzOld = e - pl;
1066 :
1067 0 : Double_t eppzNew = (1. - z) * eppzOld;
1068 0 : Double_t empzNew = empzOld - mt2 * z / eppzOld;
1069 0 : Double_t eNew = 0.5 * (eppzNew + empzNew);
1070 0 : 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 0 : Double_t mt2New = eppzNew * empzNew;
1076 0 : if (mt2New < 1.e-8) mt2New = 0.;
1077 0 : if (z < zmax) {
1078 0 : if (m * m > mt2New) {
1079 : //
1080 : // This should not happen
1081 : //
1082 0 : Fatal("Quench()", "This should never happen %e %e %e!", m, eppzNew, empzNew);
1083 : jtNew = 0;
1084 0 : } else {
1085 0 : 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 0 : 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 0 : if (jt > 0.) {
1102 0 : pxNew = jtNew / jt * pxs;
1103 0 : pyNew = jtNew / jt * pys;
1104 0 : }
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 0 : TVector3 w(pxNew, pyNew, plNew);
1117 0 : w.RotateY(thetaq[index]); w.RotateZ(phiq[index]);
1118 0 : pxNew = w.X(); pyNew = w.Y(); plNew = w.Z();
1119 :
1120 0 : p1[index][0] += pxNew;
1121 0 : p1[index][1] += pyNew;
1122 0 : p1[index][2] += plNew;
1123 0 : p1[index][3] += eNew;
1124 : //
1125 : // Updated 4-momentum vectors
1126 : //
1127 0 : pNew[icount][0] = pxNew;
1128 0 : pNew[icount][1] = pyNew;
1129 0 : pNew[icount][2] = plNew;
1130 0 : pNew[icount][3] = eNew;
1131 0 : kNew[icount] = i;
1132 0 : icount++;
1133 0 : } // parton loop
1134 : //
1135 : // Check if there was phase-space for quenching
1136 : //
1137 :
1138 0 : if (icount == 0) quenched[isys] = kFALSE;
1139 0 : if (!quenched[isys]) break;
1140 :
1141 0 : for (Int_t j = 0; j < 4; j++)
1142 : {
1143 0 : p2[isys][j] = p0[isys][j] - p1[isys][j];
1144 : }
1145 0 : 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 0 : if (p2[isys][4] > 0.) {
1147 0 : p2[isys][4] = TMath::Sqrt(p2[isys][4]);
1148 0 : break;
1149 : } else {
1150 0 : printf("Warning negative mass squared in system %d %f ! \n", isys, zquench[isys]);
1151 0 : 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 0 : if (p2[isys][4] < -0.01) {
1153 0 : 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 0 : p2[isys][4] = 0.;
1157 0 : p2[isys][3] = TMath::Sqrt(p2[isys][0] * p2[isys][0] + p2[isys][1] * p2[isys][1] + p2[isys][2] * p2[isys][2]);
1158 0 : break;
1159 : } else {
1160 0 : p2[isys][4] = 0.;
1161 0 : 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 0 : 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 0 : fPyjets->P[0][kNew[k]] = pNew[k][0];
1180 0 : fPyjets->P[1][kNew[k]] = pNew[k][1];
1181 0 : fPyjets->P[2][kNew[k]] = pNew[k][2];
1182 0 : 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 0 : if (!quenched[isys]) continue;
1190 : //
1191 : // Last parton from shower i
1192 0 : Int_t in = klast[isys];
1193 : //
1194 : // Continue if no parton in shower i selected
1195 0 : 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 0 : 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 0 : if (p2[isys][4] > 0.05) ish = 2;
1206 : //
1207 : // Position of gluons
1208 : iGlu = numpart;
1209 0 : if (iglu == 0) igMin = iGlu;
1210 : igMax = iGlu;
1211 0 : numpart += ish;
1212 0 : (fPyjets->N) += ish;
1213 :
1214 0 : if (ish == 1) {
1215 0 : fPyjets->P[0][iGlu] = p2[isys][0];
1216 0 : fPyjets->P[1][iGlu] = p2[isys][1];
1217 0 : fPyjets->P[2][iGlu] = p2[isys][2];
1218 0 : fPyjets->P[3][iGlu] = p2[isys][3];
1219 0 : fPyjets->P[4][iGlu] = p2[isys][4];
1220 :
1221 0 : fPyjets->K[0][iGlu] = 1;
1222 0 : if (iglu == nGluon[isys] - 1) fPyjets->K[0][iGlu] = 1;
1223 0 : fPyjets->K[1][iGlu] = 21;
1224 0 : fPyjets->K[2][iGlu] = fPyjets->K[2][in] + 1000;
1225 0 : fPyjets->K[3][iGlu] = -1;
1226 0 : fPyjets->K[4][iGlu] = -1;
1227 :
1228 0 : pg[0] += p2[isys][0];
1229 0 : pg[1] += p2[isys][1];
1230 0 : pg[2] += p2[isys][2];
1231 0 : pg[3] += p2[isys][3];
1232 0 : } else {
1233 : //
1234 : // Split gluon in rest frame.
1235 : //
1236 0 : Double_t bx = p2[isys][0] / p2[isys][3];
1237 0 : Double_t by = p2[isys][1] / p2[isys][3];
1238 0 : Double_t bz = p2[isys][2] / p2[isys][3];
1239 0 : Double_t pst = p2[isys][4] / 2.;
1240 : //
1241 : // Isotropic decay ????
1242 0 : Double_t cost = 2. * gRandom->Rndm() - 1.;
1243 0 : Double_t sint = TMath::Sqrt((1.-cost)*(1.+cost));
1244 0 : Double_t phis = 2. * TMath::Pi() * gRandom->Rndm();
1245 :
1246 0 : Double_t pz1 = pst * cost;
1247 0 : Double_t pz2 = -pst * cost;
1248 0 : Double_t pt1 = pst * sint;
1249 0 : Double_t pt2 = -pst * sint;
1250 0 : Double_t px1 = pt1 * TMath::Cos(phis);
1251 0 : Double_t py1 = pt1 * TMath::Sin(phis);
1252 0 : Double_t px2 = pt2 * TMath::Cos(phis);
1253 0 : Double_t py2 = pt2 * TMath::Sin(phis);
1254 :
1255 0 : fPyjets->P[0][iGlu] = px1;
1256 0 : fPyjets->P[1][iGlu] = py1;
1257 0 : fPyjets->P[2][iGlu] = pz1;
1258 0 : fPyjets->P[3][iGlu] = pst;
1259 0 : fPyjets->P[4][iGlu] = 0.;
1260 :
1261 0 : fPyjets->K[0][iGlu] = 1 ;
1262 0 : fPyjets->K[1][iGlu] = 21;
1263 0 : fPyjets->K[2][iGlu] = fPyjets->K[2][in] + 1000;
1264 0 : fPyjets->K[3][iGlu] = -1;
1265 0 : fPyjets->K[4][iGlu] = -1;
1266 :
1267 0 : fPyjets->P[0][iGlu+1] = px2;
1268 0 : fPyjets->P[1][iGlu+1] = py2;
1269 0 : fPyjets->P[2][iGlu+1] = pz2;
1270 0 : fPyjets->P[3][iGlu+1] = pst;
1271 0 : fPyjets->P[4][iGlu+1] = 0.;
1272 :
1273 0 : fPyjets->K[0][iGlu+1] = 1;
1274 0 : if (iglu == nGluon[isys] - 1) fPyjets->K[0][iGlu+1] = 1;
1275 0 : fPyjets->K[1][iGlu+1] = 21;
1276 0 : fPyjets->K[2][iGlu+1] = fPyjets->K[2][in] + 1000;
1277 0 : fPyjets->K[3][iGlu+1] = -1;
1278 0 : fPyjets->K[4][iGlu+1] = -1;
1279 0 : SetMSTU(1,0);
1280 0 : SetMSTU(2,0);
1281 : //
1282 : // Boost back
1283 : //
1284 0 : 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 0 : } // 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 0 : for (Int_t i = 0; i < numpart; i++)
1321 : {
1322 0 : kst = fPyjets->K[0][i];
1323 0 : if (kst != 1 && kst != 2) continue;
1324 0 : pxs += fPyjets->P[0][i];
1325 0 : pys += fPyjets->P[1][i];
1326 0 : pzs += fPyjets->P[2][i];
1327 0 : es -= fPyjets->P[3][i];
1328 0 : }
1329 0 : if (TMath::Abs(pxs) > 1.e-2 ||
1330 0 : TMath::Abs(pys) > 1.e-2 ||
1331 0 : TMath::Abs(pzs) > 1.e-1) {
1332 0 : printf("%e %e %e %e\n", pxs, pys, pzs, es);
1333 : // Fatal("Quench()", "4-Momentum non-conservation");
1334 0 : }
1335 :
1336 0 : } // end quenching loop (systems)
1337 : // Clean-up
1338 0 : for (Int_t i = 0; i < numpart; i++)
1339 : {
1340 0 : imo = fPyjets->K[2][i];
1341 0 : if (imo > 1000) {
1342 0 : fPyjets->K[2][i] = fPyjets->K[2][i] % 1000;
1343 0 : }
1344 : }
1345 : // this->Pylist(1);
1346 0 : } // end quench
1347 :
1348 :
1349 0 : 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 0 : pyquen(a, ibf, b);
1355 0 : }
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 0 : PYQPAR.t0 = t0;
1364 0 : PYQPAR.tau0 = tau0;
1365 0 : PYQPAR.nf = nf;
1366 0 : PYQPAR.iengl = iengl;
1367 0 : PYQPAR.iangl = iangl;
1368 0 : }
1369 :
1370 0 : void AliPythia6::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
1371 : {
1372 : //
1373 : // Load event into Pythia Common Block
1374 : //
1375 :
1376 0 : Int_t npart = stack -> GetNprimary();
1377 : Int_t n0 = 0;
1378 :
1379 0 : if (!flag) {
1380 0 : GetPyjets()->N = npart;
1381 0 : } else {
1382 0 : n0 = GetPyjets()->N;
1383 0 : GetPyjets()->N = n0 + npart;
1384 : }
1385 :
1386 :
1387 0 : for (Int_t part = 0; part < npart; part++) {
1388 0 : TParticle *mPart = stack->Particle(part);
1389 :
1390 0 : Int_t kf = mPart->GetPdgCode();
1391 0 : Int_t ks = mPart->GetStatusCode();
1392 0 : Int_t idf = mPart->GetFirstDaughter();
1393 0 : Int_t idl = mPart->GetLastDaughter();
1394 :
1395 0 : if (reHadr) {
1396 0 : if (ks == 11 || ks == 12) {
1397 0 : ks -= 10;
1398 : idf = -1;
1399 : idl = -1;
1400 0 : }
1401 : }
1402 :
1403 0 : Float_t px = mPart->Px();
1404 0 : Float_t py = mPart->Py();
1405 0 : Float_t pz = mPart->Pz();
1406 0 : Float_t e = mPart->Energy();
1407 0 : Float_t m = mPart->GetCalcMass();
1408 :
1409 :
1410 0 : (GetPyjets())->P[0][part+n0] = px;
1411 0 : (GetPyjets())->P[1][part+n0] = py;
1412 0 : (GetPyjets())->P[2][part+n0] = pz;
1413 0 : (GetPyjets())->P[3][part+n0] = e;
1414 0 : (GetPyjets())->P[4][part+n0] = m;
1415 :
1416 0 : (GetPyjets())->K[1][part+n0] = kf;
1417 0 : (GetPyjets())->K[0][part+n0] = ks;
1418 0 : (GetPyjets())->K[3][part+n0] = idf + 1;
1419 0 : (GetPyjets())->K[4][part+n0] = idl + 1;
1420 0 : (GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1421 : }
1422 0 : }
1423 :
1424 :
1425 : void AliPythia6::Pyevnw()
1426 : {
1427 : // New multiple interaction scenario
1428 0 : pyevnw();
1429 0 : }
1430 :
1431 0 : void AliPythia6::GetQuenchingParameters(Double_t& xp, Double_t& yp, Double_t z[4])
1432 : {
1433 : // Return event specific quenching parameters
1434 0 : xp = fXJet;
1435 0 : yp = fYJet;
1436 0 : for (Int_t i = 0; i < 4; i++) z[i] = fZQuench[i];
1437 :
1438 0 : }
1439 :
1440 0 : void AliPythia6::ConfigHeavyFlavor()
1441 : {
1442 : //
1443 : // Default configuration for Heavy Flavor production
1444 : //
1445 : // All QCD processes
1446 : //
1447 0 : SetMSEL(1);
1448 :
1449 : // No multiple interactions
1450 0 : SetMSTP(81,0);
1451 0 : SetPARP(81, 0.);
1452 0 : SetPARP(82, 0.);
1453 : // Initial/final parton shower on (Pythia default)
1454 0 : SetMSTP(61,1);
1455 0 : SetMSTP(71,1);
1456 :
1457 : // 2nd order alpha_s
1458 0 : SetMSTP(2,2);
1459 :
1460 : // QCD scales
1461 0 : SetMSTP(32,2);
1462 0 : SetPARP(34,1.0);
1463 0 : }
1464 :
1465 0 : void AliPythia6::AtlasTuning()
1466 : {
1467 : //
1468 : // Configuration for the ATLAS tuning
1469 0 : SetMSTP(51,AliStructFuncType::PDFsetIndex(kCTEQ5L));
1470 0 : SetMSTP(81,1); // Multiple Interactions ON
1471 0 : SetMSTP(82,4); // Double Gaussian Model
1472 0 : SetPARP(81,1.9); // Min. pt for multiple interactions (default in 6.2-14)
1473 0 : SetPARP(82,1.8); // [GeV] PT_min at Ref. energy
1474 0 : SetPARP(89,1000.); // [GeV] Ref. energy
1475 0 : SetPARP(90,0.16); // 2*epsilon (exponent in power law)
1476 0 : SetPARP(83,0.5); // Core density in proton matter distribution (def.value)
1477 0 : SetPARP(84,0.5); // Core radius
1478 0 : SetPARP(85,0.33); // Regulates gluon prod. mechanism
1479 0 : SetPARP(86,0.66); // Regulates gluon prod. mechanism
1480 0 : SetPARP(67,1); // Regulates Initial State Radiation
1481 0 : }
1482 :
1483 : void AliPythia6::AtlasTuningMC09()
1484 : {
1485 : //
1486 : // Configuration for the ATLAS tuning
1487 0 : printf("ATLAS New TUNE MC09\n");
1488 0 : SetMSTP(81,21); // treatment for MI, ISR, FSR and beam remnants: MI on, new model
1489 0 : SetMSTP(82, 4); // Double Gaussian Model
1490 0 : SetMSTP(52, 2); // External PDF
1491 0 : SetMSTP(51, 20650); // MRST LO*
1492 :
1493 :
1494 0 : SetMSTP(70, 0); // (was 2: def manual 1, def code 0) virtuality scale for ISR
1495 0 : SetMSTP(72, 1); // (was 0: def 1) maximum scale for FSR
1496 0 : SetMSTP(88, 1); // (was 0: def 1) strategy for qq junction to di-quark or baryon in beam remnant
1497 0 : SetMSTP(90, 0); // (was 1: def 0) strategy of compensate the primordial kT
1498 :
1499 0 : SetPARP(78, 0.3); // the amount of color reconnection in the final state
1500 0 : SetPARP(80, 0.1); // probability of color partons kicked out from beam remnant
1501 0 : SetPARP(82, 2.3); // [GeV] PT_min at Ref. energy
1502 0 : SetPARP(83, 0.8); // Core density in proton matter distribution (def.value)
1503 0 : SetPARP(84, 0.7); // Core radius
1504 0 : SetPARP(90, 0.25); // 2*epsilon (exponent in power law)
1505 0 : SetPARJ(81, 0.29); // (was 0.14: def 0.29) Labmda value in running alpha_s for parton showers
1506 :
1507 0 : SetMSTP(95, 6);
1508 0 : SetPARJ(41, 0.3); // a and b parameters of the symmm. Lund FF
1509 0 : SetPARJ(42, 0.58);
1510 0 : SetPARJ(46, 0.75); // mod. of the Lund FF for heavy end-point quarks
1511 0 : SetPARP(89,1800.); // [GeV] Ref. energy
1512 0 : }
1513 :
1514 : void AliPythia6::SetWeightPower(Double_t pow)
1515 : {
1516 0 : setpowwght(pow);
1517 0 : SetMSTP(142, 1); // Tell Pythia to use pyevwt to calculate event wghts
1518 0 : }
1519 :
1520 0 : void AliPythia6::SetPtHardRange(Float_t ptmin, Float_t ptmax)
1521 : {
1522 : // Set the pt hard range
1523 0 : SetCKIN(3, ptmin);
1524 0 : SetCKIN(4, ptmax);
1525 0 : }
1526 :
1527 0 : void AliPythia6::SetYHardRange(Float_t ymin, Float_t ymax)
1528 : {
1529 : // Set the y hard range
1530 0 : SetCKIN(7, ymin);
1531 0 : SetCKIN(8, ymax);
1532 0 : }
1533 :
1534 :
1535 0 : void AliPythia6::SetFragmentation(Int_t flag)
1536 : {
1537 : // Switch fragmentation on/off
1538 0 : SetMSTP(111, flag);
1539 0 : }
1540 :
1541 0 : void AliPythia6::SetInitialAndFinalStateRadiation(Int_t flag1, Int_t flag2)
1542 : {
1543 : // initial state radiation
1544 0 : SetMSTP(61, flag1);
1545 : // final state radiation
1546 0 : SetMSTP(71, flag2);
1547 0 : }
1548 :
1549 0 : void AliPythia6::SetIntrinsicKt(Float_t kt)
1550 : {
1551 : // Set the inreinsic kt
1552 0 : if (kt > 0.) {
1553 0 : SetMSTP(91,1);
1554 0 : SetPARP(91,kt);
1555 0 : SetPARP(93, 4. * kt);
1556 0 : } else {
1557 0 : SetMSTP(91,0);
1558 : }
1559 0 : }
1560 :
1561 0 : void AliPythia6::SwitchHFOff()
1562 : {
1563 : // Switch off heavy flavor
1564 : // Maximum number of quark flavours used in pdf
1565 0 : SetMSTP(58, 3);
1566 : // Maximum number of flavors that can be used in showers
1567 0 : SetMSTJ(45, 3);
1568 0 : }
1569 :
1570 0 : 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 0 : SetPARU(51, etamax);
1575 0 : SetMSTU(51, neta);
1576 0 : SetMSTU(52, nphi);
1577 0 : SetPARU(58, thresh);
1578 0 : SetPARU(52, etseed);
1579 0 : SetPARU(53, minet);
1580 0 : SetPARU(54, r);
1581 0 : SetMSTU(54, 2);
1582 0 : }
1583 :
1584 0 : void AliPythia6::ModifiedSplitting()
1585 : {
1586 : // Modified splitting probability as a model for quenching
1587 0 : SetPARJ(200, 0.8);
1588 0 : SetMSTJ(41, 1); // QCD radiation only
1589 0 : SetMSTJ(42, 2); // angular ordering
1590 0 : SetMSTJ(44, 2); // option to run alpha_s
1591 0 : SetMSTJ(47, 0); // No correction back to hard scattering element
1592 0 : SetMSTJ(50, 0); // No coherence in first branching
1593 0 : SetPARJ(82, 1.); // Cut off for parton showers
1594 0 : }
1595 :
1596 0 : void AliPythia6::SwitchHadronisationOff()
1597 : {
1598 : // Switch off hadronisarion
1599 0 : SetMSTJ(1, 0);
1600 0 : }
1601 :
1602 0 : void AliPythia6::SwitchHadronisationOn()
1603 : {
1604 : // Switch on hadronisarion
1605 0 : SetMSTJ(1, 1);
1606 0 : }
1607 :
1608 :
1609 0 : void AliPythia6::GetXandQ(Float_t& x1, Float_t& x2, Float_t& q)
1610 : {
1611 : // Get x1, x2 and Q for this event
1612 :
1613 0 : q = GetVINT(51);
1614 0 : x1 = GetVINT(41);
1615 0 : x2 = GetVINT(42);
1616 0 : }
1617 :
1618 0 : Float_t AliPythia6::GetXSection()
1619 : {
1620 : // Get the total cross-section
1621 0 : return (GetPARI(1));
1622 : }
1623 :
1624 0 : Float_t AliPythia6::GetPtHard()
1625 : {
1626 : // Get the pT hard for this event
1627 0 : return GetVINT(47);
1628 : }
1629 :
1630 0 : Int_t AliPythia6::ProcessCode()
1631 : {
1632 : // Get the subprocess code
1633 0 : return GetMSTI(1);
1634 : }
1635 :
1636 0 : void AliPythia6::PrintStatistics()
1637 : {
1638 : // End of run statistics
1639 0 : Pystat(1);
1640 0 : }
1641 :
1642 0 : void AliPythia6::EventListing()
1643 : {
1644 : // End of run statistics
1645 0 : Pylist(2);
1646 0 : }
1647 :
1648 : AliPythia6& AliPythia6::operator=(const AliPythia6& rhs)
1649 : {
1650 : // Assignment operator
1651 0 : rhs.Copy(*this);
1652 0 : return *this;
1653 : }
1654 :
1655 : void AliPythia6::Copy(TObject&) const
1656 : {
1657 : //
1658 : // Copy
1659 : //
1660 0 : Fatal("Copy","Not implemented!\n");
1661 0 : }
1662 :
1663 0 : Int_t AliPythia6::GetNMPI()
1664 : {
1665 0 : return (GetMSTI(31));
1666 : }
|
