Fix in LoadCTPTimeParamsCDB() (Plamen)
[u/mrichter/AliRoot.git] / PYTHIA6 / AliPythia.cxx
CommitLineData
6b435cde 1
8d2cd130 2/**************************************************************************
3 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * *
5 * Author: The ALICE Off-line Project. *
6 * Contributors are mentioned in the code where appropriate. *
7 * *
8 * Permission to use, copy, modify and distribute this software and its *
9 * documentation strictly for non-commercial purposes is hereby granted *
10 * without fee, provided that the above copyright notice appears in all *
11 * copies and that both the copyright notice and this permission notice *
12 * appear in the supporting documentation. The authors make no claims *
13 * about the suitability of this software for any purpose. It is *
14 * provided "as is" without express or implied warranty. *
15 **************************************************************************/
16
7cdba479 17/* $Id$ */
8d2cd130 18
19#include "AliPythia.h"
7cdba479 20#include "AliPythiaRndm.h"
12cb0bc0 21#include "AliFastGlauber.h"
22#include "AliQuenchingWeights.h"
0f482ae4 23#include "TVector3.h"
12cb0bc0 24#include "PyquenCommon.h"
8d2cd130 25
26ClassImp(AliPythia)
27
28#ifndef WIN32
29# define pyclus pyclus_
30# define pycell pycell_
452af8c7 31# define pyshow pyshow_
32# define pyrobo pyrobo_
992f2843 33# define pyquen pyquen_
16a82508 34# define pyevnw pyevnw_
cd07c39b 35# define pyshowq pyshowq_
694b39f9 36# define pytune pytune_
9b61ba2a 37# define py2ent py2ent_
8d2cd130 38# define type_of_call
39#else
40# define pyclus PYCLUS
41# define pycell PYCELL
452af8c7 42# define pyrobo PYROBO
992f2843 43# define pyquen PYQUEN
16a82508 44# define pyevnw PYEVNW
cd07c39b 45# define pyshowq PYSHOWQ
9b61ba2a 46# define pytune PYTUNE
47# define py2ent PY2ENT
8d2cd130 48# define type_of_call _stdcall
49#endif
50
51extern "C" void type_of_call pyclus(Int_t & );
52extern "C" void type_of_call pycell(Int_t & );
452af8c7 53extern "C" void type_of_call pyshow(Int_t &, Int_t &, Double_t &);
54extern "C" void type_of_call pyrobo(Int_t &, Int_t &, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &);
992f2843 55extern "C" void type_of_call pyquen(Double_t &, Int_t &, Double_t &);
0a2cfc0a 56extern "C" void type_of_call pyevnw(){;}
cd07c39b 57extern "C" void type_of_call pyshowq(Int_t &, Int_t &, Double_t &);
694b39f9 58extern "C" void type_of_call pytune(Int_t &);
9b61ba2a 59extern "C" void type_of_call py2ent(Int_t &, Int_t&, Int_t&, Double_t&);
8d2cd130 60
61//_____________________________________________________________________________
62
63AliPythia* AliPythia::fgAliPythia=NULL;
64
e8a8adcd 65AliPythia::AliPythia():
66 fProcess(kPyMb),
67 fEcms(0.),
68 fStrucFunc(kCTEQ5L),
69 fXJet(0.),
70 fYJet(0.),
32c8e463 71 fNGmax(30),
72 fZmax(0.97),
e8a8adcd 73 fGlauber(0),
74 fQuenchingWeights(0)
8d2cd130 75{
76// Default Constructor
77//
78// Set random number
7cdba479 79 if (!AliPythiaRndm::GetPythiaRandom())
80 AliPythiaRndm::SetPythiaRandom(GetRandom());
0f482ae4 81 fGlauber = 0;
82 fQuenchingWeights = 0;
8d2cd130 83}
84
e8a8adcd 85AliPythia::AliPythia(const AliPythia& pythia):
6b435cde 86 TPythia6(pythia),
87 AliRndm(pythia),
e8a8adcd 88 fProcess(kPyMb),
89 fEcms(0.),
90 fStrucFunc(kCTEQ5L),
91 fXJet(0.),
92 fYJet(0.),
32c8e463 93 fNGmax(30),
94 fZmax(0.97),
e8a8adcd 95 fGlauber(0),
96 fQuenchingWeights(0)
97{
98 // Copy Constructor
99 pythia.Copy(*this);
100}
101
efe3b1cd 102void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfunc, Int_t itune)
8d2cd130 103{
104// Initialise the process to generate
7cdba479 105 if (!AliPythiaRndm::GetPythiaRandom())
106 AliPythiaRndm::SetPythiaRandom(GetRandom());
8d2cd130 107
108 fProcess = process;
109 fEcms = energy;
110 fStrucFunc = strucfunc;
1d5b1b20 111//...Switch off decay of pi0, K0S, Lambda, Sigma+-, Xi0-, Omega-.
e2de0ce1 112 SetMDCY(Pycomp(111) ,1,0); // pi0
113 SetMDCY(Pycomp(310) ,1,0); // K0S
114 SetMDCY(Pycomp(3122),1,0); // kLambda
115 SetMDCY(Pycomp(3112),1,0); // sigma -
116 SetMDCY(Pycomp(3212),1,0); // sigma 0
117 SetMDCY(Pycomp(3222),1,0); // sigma +
118 SetMDCY(Pycomp(3312),1,0); // xi -
119 SetMDCY(Pycomp(3322),1,0); // xi 0
120 SetMDCY(Pycomp(3334),1,0); // omega-
1c50ec12 121 // Select structure function
8d2cd130 122 SetMSTP(52,2);
e2de0ce1 123 SetMSTP(51, AliStructFuncType::PDFsetIndex(strucfunc));
1c50ec12 124 // Particles produced in string fragmentation point directly to either of the two endpoints
125 // of the string (depending in the side they were generated from).
126 SetMSTU(16,2);
127
8d2cd130 128//
129// Pythia initialisation for selected processes//
130//
131// Make MSEL clean
132//
133 for (Int_t i=1; i<= 200; i++) {
134 SetMSUB(i,0);
135 }
136// select charm production
137 switch (process)
138 {
65f2626c 139 case kPyOldUEQ2ordered: //Old underlying events with Q2 ordered QCD processes
140// Multiple interactions on.
141 SetMSTP(81,1);
142// Double Gaussian matter distribution.
143 SetMSTP(82,4);
144 SetPARP(83,0.5);
145 SetPARP(84,0.4);
146// pT0.
147 SetPARP(82,2.0);
148// Reference energy for pT0 and energy rescaling pace.
149 SetPARP(89,1800);
150 SetPARP(90,0.25);
151// String drawing almost completely minimizes string length.
152 SetPARP(85,0.9);
153 SetPARP(86,0.95);
154// ISR and FSR activity.
155 SetPARP(67,4);
156 SetPARP(71,4);
157// Lambda_FSR scale.
158 SetPARJ(81,0.29);
159 break;
160 case kPyOldUEQ2ordered2:
161// 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.16); // here is the difference with kPyOldUEQ2ordered
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 kPyOldPopcorn:
183// Old production mechanism: Old Popcorn
184 SetMSEL(1);
185 SetMSTJ(12,3);
186// (D=2) Like MSTJ(12)=2 but added prod ofthe 1er rank baryon
187 SetMSTP(88,2);
188// (D=1)see can be used to form baryons (BARYON JUNCTION)
189 SetMSTJ(1,1);
e0e89f40 190 AtlasTuning();
65f2626c 191 break;
8d2cd130 192 case kPyCharm:
193 SetMSEL(4);
8d2cd130 194// heavy quark masses
195
196 SetPMAS(4,1,1.2);
8d2cd130 197//
198// primordial pT
199 SetMSTP(91,1);
200 SetPARP(91,1.);
201 SetPARP(93,5.);
202//
203 break;
204 case kPyBeauty:
205 SetMSEL(5);
206 SetPMAS(5,1,4.75);
8d2cd130 207 break;
208 case kPyJpsi:
209 SetMSEL(0);
210// gg->J/Psi g
211 SetMSUB(86,1);
212 break;
213 case kPyJpsiChi:
214 SetMSEL(0);
215// gg->J/Psi g
216 SetMSUB(86,1);
217// gg-> chi_0c g
218 SetMSUB(87,1);
219// gg-> chi_1c g
220 SetMSUB(88,1);
221// gg-> chi_2c g
222 SetMSUB(89,1);
223 break;
224 case kPyCharmUnforced:
225 SetMSEL(0);
226// gq->qg
227 SetMSUB(28,1);
228// gg->qq
229 SetMSUB(53,1);
230// gg->gg
231 SetMSUB(68,1);
232 break;
233 case kPyBeautyUnforced:
234 SetMSEL(0);
235// gq->qg
236 SetMSUB(28,1);
237// gg->qq
238 SetMSUB(53,1);
239// gg->gg
240 SetMSUB(68,1);
241 break;
242 case kPyMb:
243// Minimum Bias pp-Collisions
244//
245//
246// select Pythia min. bias model
247 SetMSEL(0);
511db649 248 SetMSUB(92,1); // single diffraction AB-->XB
249 SetMSUB(93,1); // single diffraction AB-->AX
250 SetMSUB(94,1); // double diffraction
251 SetMSUB(95,1); // low pt production
252
e0e89f40 253 AtlasTuning();
511db649 254 break;
04081a91 255
256 case kPyMbWithDirectPhoton:
257// Minimum Bias pp-Collisions with direct photon processes added
258//
259//
260// select Pythia min. bias model
261 SetMSEL(0);
262 SetMSUB(92,1); // single diffraction AB-->XB
263 SetMSUB(93,1); // single diffraction AB-->AX
264 SetMSUB(94,1); // double diffraction
265 SetMSUB(95,1); // low pt production
266
267 SetMSUB(14,1); //
268 SetMSUB(18,1); //
269 SetMSUB(29,1); //
270 SetMSUB(114,1); //
271 SetMSUB(115,1); //
272
273
274 AtlasTuning();
275 break;
276
f529e69b 277 case kPyMbDefault:
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 break;
0a0cbcfd 289 case kPyLhwgMb:
290// Les Houches Working Group 05 Minimum Bias pp-Collisions: hep-ph/0604120
291// -> Pythia 6.3 or above is needed
292//
293 SetMSEL(0);
294 SetMSUB(92,1); // single diffraction AB-->XB
295 SetMSUB(93,1); // single diffraction AB-->AX
296 SetMSUB(94,1); // double diffraction
297 SetMSUB(95,1); // low pt production
298
e2de0ce1 299 SetMSTP(51,AliStructFuncType::PDFsetIndex(kCTEQ6ll)); // CTEQ6ll pdf
0a0cbcfd 300 SetMSTP(52,2);
301 SetMSTP(68,1);
302 SetMSTP(70,2);
303 SetMSTP(81,1); // Multiple Interactions ON
304 SetMSTP(82,4); // Double Gaussian Model
305 SetMSTP(88,1);
306
307 SetPARP(82,2.3); // [GeV] PT_min at Ref. energy
308 SetPARP(83,0.5); // Core density in proton matter distribution (def.value)
309 SetPARP(84,0.5); // Core radius
310 SetPARP(85,0.9); // Regulates gluon prod. mechanism
311 SetPARP(90,0.2); // 2*epsilon (exponent in power law)
312
313 break;
8d2cd130 314 case kPyMbNonDiffr:
315// Minimum Bias pp-Collisions
316//
317//
318// select Pythia min. bias model
319 SetMSEL(0);
511db649 320 SetMSUB(95,1); // low pt production
0f482ae4 321
e0e89f40 322 AtlasTuning();
8d2cd130 323 break;
d7de4a67 324 case kPyMbMSEL1:
325 ConfigHeavyFlavor();
326// Intrinsic <kT^2>
327 SetMSTP(91,1);// Width (1=gaussian) primordial kT dist. inside hadrons
328 SetPARP(91,1.); // <kT^2> = PARP(91,1.)^2
329 SetPARP(93,5.); // Upper cut-off
330// Set Q-quark mass
331 SetPMAS(4,1,1.2); // Charm quark mass
332 SetPMAS(5,1,4.78); // Beauty quark mass
333 SetPARP(71,4.); // Defaut value
334// Atlas Tuning
335 AtlasTuning();
336 break;
8d2cd130 337 case kPyJets:
338//
339// QCD Jets
340//
341 SetMSEL(1);
65f2626c 342 // Pythia Tune A (CDF)
343 //
4167b79f 344 SetPARP(67,2.5); // Regulates Initial State Radiation (value from best fit to D0 dijet analysis)
65f2626c 345 SetMSTP(82,4); // Double Gaussian Model
346 SetPARP(82,2.0); // [GeV] PT_min at Ref. energy
347 SetPARP(84,0.4); // Core radius
348 SetPARP(85,0.90) ; // Regulates gluon prod. mechanism
349 SetPARP(86,0.95); // Regulates gluon prod. mechanism
350 SetPARP(89,1800.); // [GeV] Ref. energy
351 SetPARP(90,0.25); // 2*epsilon (exponent in power law)
352 break;
8d2cd130 353 case kPyDirectGamma:
354 SetMSEL(10);
355 break;
adf4d898 356 case kPyCharmPbPbMNR:
357 case kPyD0PbPbMNR:
90d7b703 358 case kPyDPlusPbPbMNR:
e0e89f40 359 case kPyDPlusStrangePbPbMNR:
90d7b703 360 // Tuning of Pythia parameters aimed to get a resonable agreement
361 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
362 // c-cbar single inclusive and double differential distributions.
363 // This parameter settings are meant to work with Pb-Pb collisions
364 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
365 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
366 // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
3dc3ec94 367 ConfigHeavyFlavor();
90d7b703 368 // Intrinsic <kT>
369 SetMSTP(91,1);
370 SetPARP(91,1.304);
371 SetPARP(93,6.52);
90d7b703 372 // Set c-quark mass
373 SetPMAS(4,1,1.2);
90d7b703 374 break;
adf4d898 375 case kPyCharmpPbMNR:
376 case kPyD0pPbMNR:
90d7b703 377 case kPyDPluspPbMNR:
e0e89f40 378 case kPyDPlusStrangepPbMNR:
90d7b703 379 // Tuning of Pythia parameters aimed to get a resonable agreement
380 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
381 // c-cbar single inclusive and double differential distributions.
382 // This parameter settings are meant to work with p-Pb collisions
383 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
384 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
385 // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
3dc3ec94 386 ConfigHeavyFlavor();
90d7b703 387 // Intrinsic <kT>
3dc3ec94 388 SetMSTP(91,1);
389 SetPARP(91,1.16);
390 SetPARP(93,5.8);
391
90d7b703 392 // Set c-quark mass
3dc3ec94 393 SetPMAS(4,1,1.2);
90d7b703 394 break;
adf4d898 395 case kPyCharmppMNR:
396 case kPyD0ppMNR:
90d7b703 397 case kPyDPlusppMNR:
e0e89f40 398 case kPyDPlusStrangeppMNR:
90d7b703 399 // Tuning of Pythia parameters aimed to get a resonable agreement
400 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
401 // c-cbar single inclusive and double differential distributions.
402 // This parameter settings are meant to work with pp collisions
403 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
404 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
405 // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
3dc3ec94 406 ConfigHeavyFlavor();
90d7b703 407 // Intrinsic <kT^2>
3dc3ec94 408 SetMSTP(91,1);
409 SetPARP(91,1.);
410 SetPARP(93,5.);
411
90d7b703 412 // Set c-quark mass
3dc3ec94 413 SetPMAS(4,1,1.2);
90d7b703 414 break;
e0e89f40 415 case kPyCharmppMNRwmi:
416 // Tuning of Pythia parameters aimed to get a resonable agreement
417 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
418 // c-cbar single inclusive and double differential distributions.
419 // This parameter settings are meant to work with pp collisions
420 // and with kCTEQ5L PDFs.
421 // Added multiple interactions according to ATLAS tune settings.
422 // To get a "reasonable" agreement with MNR results, events have to be
423 // generated with the minimum ptHard (AliGenPythia::SetPtHard)
424 // set to 2.76 GeV.
425 // To get a "perfect" agreement with MNR results, events have to be
426 // generated in four ptHard bins with the following relative
427 // normalizations:
428 // 2.76-3 GeV: 25%
429 // 3-4 GeV: 40%
430 // 4-8 GeV: 29%
431 // >8 GeV: 6%
432 ConfigHeavyFlavor();
433 // Intrinsic <kT^2>
434 SetMSTP(91,1);
435 SetPARP(91,1.);
436 SetPARP(93,5.);
437
438 // Set c-quark mass
439 SetPMAS(4,1,1.2);
440 AtlasTuning();
441 break;
adf4d898 442 case kPyBeautyPbPbMNR:
8d2cd130 443 // Tuning of Pythia parameters aimed to get a resonable agreement
444 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
445 // b-bbar single inclusive and double differential distributions.
446 // This parameter settings are meant to work with Pb-Pb collisions
447 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
448 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
449 // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
3dc3ec94 450 ConfigHeavyFlavor();
8d2cd130 451 // QCD scales
3dc3ec94 452 SetPARP(67,1.0);
453 SetPARP(71,1.0);
adf4d898 454 // Intrinsic <kT>
3dc3ec94 455 SetMSTP(91,1);
456 SetPARP(91,2.035);
457 SetPARP(93,10.17);
8d2cd130 458 // Set b-quark mass
3dc3ec94 459 SetPMAS(5,1,4.75);
8d2cd130 460 break;
adf4d898 461 case kPyBeautypPbMNR:
462 // Tuning of Pythia parameters aimed to get a resonable agreement
463 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
464 // b-bbar single inclusive and double differential distributions.
465 // This parameter settings are meant to work with p-Pb collisions
466 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
467 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
468 // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
3dc3ec94 469 ConfigHeavyFlavor();
adf4d898 470 // QCD scales
3dc3ec94 471 SetPARP(67,1.0);
472 SetPARP(71,1.0);
adf4d898 473 // Intrinsic <kT>
3dc3ec94 474 SetMSTP(91,1);
475 SetPARP(91,1.60);
476 SetPARP(93,8.00);
adf4d898 477 // Set b-quark mass
3dc3ec94 478 SetPMAS(5,1,4.75);
adf4d898 479 break;
480 case kPyBeautyppMNR:
481 // Tuning of Pythia parameters aimed to get a resonable agreement
482 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
483 // b-bbar single inclusive and double differential distributions.
484 // This parameter settings are meant to work with pp collisions
485 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
486 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
487 // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
3dc3ec94 488 ConfigHeavyFlavor();
adf4d898 489 // QCD scales
3dc3ec94 490 SetPARP(67,1.0);
491 SetPARP(71,1.0);
492
493 // Intrinsic <kT>
494 SetMSTP(91,1);
495 SetPARP(91,1.);
496 SetPARP(93,5.);
497
498 // Set b-quark mass
499 SetPMAS(5,1,4.75);
adf4d898 500 break;
9dfe63b3 501 case kPyBeautyJets:
e0e89f40 502 case kPyBeautyppMNRwmi:
503 // Tuning of Pythia parameters aimed to get a resonable agreement
504 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
505 // b-bbar single inclusive and double differential distributions.
506 // This parameter settings are meant to work with pp collisions
507 // and with kCTEQ5L PDFs.
508 // Added multiple interactions according to ATLAS tune settings.
509 // To get a "reasonable" agreement with MNR results, events have to be
510 // generated with the minimum ptHard (AliGenPythia::SetPtHard)
511 // set to 2.76 GeV.
512 // To get a "perfect" agreement with MNR results, events have to be
513 // generated in four ptHard bins with the following relative
514 // normalizations:
515 // 2.76-4 GeV: 5%
516 // 4-6 GeV: 31%
517 // 6-8 GeV: 28%
518 // >8 GeV: 36%
519 ConfigHeavyFlavor();
520 // QCD scales
521 SetPARP(67,1.0);
522 SetPARP(71,1.0);
523
524 // Intrinsic <kT>
525 SetMSTP(91,1);
526 SetPARP(91,1.);
527 SetPARP(93,5.);
528
529 // Set b-quark mass
530 SetPMAS(5,1,4.75);
531
532 AtlasTuning();
533 break;
589380c6 534 case kPyW:
535
536 //Inclusive production of W+/-
537 SetMSEL(0);
538 //f fbar -> W+
539 SetMSUB(2,1);
540 // //f fbar -> g W+
541 // SetMSUB(16,1);
542 // //f fbar -> gamma W+
543 // SetMSUB(20,1);
544 // //f g -> f W+
545 // SetMSUB(31,1);
546 // //f gamma -> f W+
547 // SetMSUB(36,1);
548
549 // Initial/final parton shower on (Pythia default)
550 // With parton showers on we are generating "W inclusive process"
551 SetMSTP(61,1); //Initial QCD & QED showers on
552 SetMSTP(71,1); //Final QCD & QED showers on
553
554 break;
0f6ee828 555
556 case kPyZ:
557
558 //Inclusive production of Z
559 SetMSEL(0);
560 //f fbar -> Z/gamma
561 SetMSUB(1,1);
562
563 // // f fbar -> g Z/gamma
564 // SetMSUB(15,1);
565 // // f fbar -> gamma Z/gamma
566 // SetMSUB(19,1);
567 // // f g -> f Z/gamma
568 // SetMSUB(30,1);
569 // // f gamma -> f Z/gamma
570 // SetMSUB(35,1);
571
572 //only Z included, not gamma
573 SetMSTP(43,2);
574
575 // Initial/final parton shower on (Pythia default)
576 // With parton showers on we are generating "Z inclusive process"
577 SetMSTP(61,1); //Initial QCD & QED showers on
578 SetMSTP(71,1); //Final QCD & QED showers on
579
580 break;
581
8d2cd130 582 }
583//
584// Initialize PYTHIA
efe3b1cd 585//
586// Select the tune
587 if (itune > -1) Pytune(itune);
588
589//
8d2cd130 590 SetMSTP(41,1); // all resonance decays switched on
8d2cd130 591 Initialize("CMS","p","p",fEcms);
4167b79f 592
8d2cd130 593}
594
595Int_t AliPythia::CheckedLuComp(Int_t kf)
596{
597// Check Lund particle code (for debugging)
598 Int_t kc=Pycomp(kf);
599 printf("\n Lucomp kf,kc %d %d",kf,kc);
600 return kc;
601}
602
20e47f08 603void AliPythia::SetNuclei(Int_t a1, Int_t a2, Int_t pdf)
8d2cd130 604{
605// Treat protons as inside nuclei with mass numbers a1 and a2
606// The MSTP array in the PYPARS common block is used to enable and
607// select the nuclear structure functions.
608// MSTP(52) : (D=1) choice of proton and nuclear structure-function library
609// =1: internal PYTHIA acording to MSTP(51)
610// =2: PDFLIB proton s.f., with MSTP(51) = 1000xNGROUP+NSET
611// If the following mass number both not equal zero, nuclear corrections of the stf are used.
612// MSTP(192) : Mass number of nucleus side 1
613// MSTP(193) : Mass number of nucleus side 2
20e47f08 614// MSTP(194) : Nuclear structure function: 0: EKS98 1:EPS08
8d2cd130 615 SetMSTP(52,2);
616 SetMSTP(192, a1);
20e47f08 617 SetMSTP(193, a2);
618 SetMSTP(194, pdf);
8d2cd130 619}
620
621
622AliPythia* AliPythia::Instance()
623{
624// Set random number generator
625 if (fgAliPythia) {
626 return fgAliPythia;
627 } else {
628 fgAliPythia = new AliPythia();
629 return fgAliPythia;
630 }
631}
632
633void AliPythia::PrintParticles()
634{
635// Print list of particl properties
636 Int_t np = 0;
c31f1d37 637 char* name = new char[16];
8d2cd130 638 for (Int_t kf=0; kf<1000000; kf++) {
639 for (Int_t c = 1; c > -2; c-=2) {
8d2cd130 640 Int_t kc = Pycomp(c*kf);
641 if (kc) {
642 Float_t mass = GetPMAS(kc,1);
643 Float_t width = GetPMAS(kc,2);
644 Float_t tau = GetPMAS(kc,4);
c31f1d37 645
8d2cd130 646 Pyname(kf,name);
647
648 np++;
649
650 printf("\n mass, width, tau: %6d %s %10.3f %10.3e %10.3e",
651 c*kf, name, mass, width, tau);
652 }
653 }
654 }
655 printf("\n Number of particles %d \n \n", np);
656}
657
658void AliPythia::ResetDecayTable()
659{
660// Set default values for pythia decay switches
661 Int_t i;
662 for (i = 1; i < 501; i++) SetMDCY(i,1,fDefMDCY[i]);
663 for (i = 1; i < 2001; i++) SetMDME(i,1,fDefMDME[i]);
664}
665
666void AliPythia::SetDecayTable()
667{
668// Set default values for pythia decay switches
669//
670 Int_t i;
671 for (i = 1; i < 501; i++) fDefMDCY[i] = GetMDCY(i,1);
672 for (i = 1; i < 2001; i++) fDefMDME[i] = GetMDME(i,1);
673}
674
675void AliPythia::Pyclus(Int_t& njet)
676{
677// Call Pythia clustering algorithm
678//
679 pyclus(njet);
680}
681
682void AliPythia::Pycell(Int_t& njet)
683{
684// Call Pythia jet reconstruction algorithm
685//
686 pycell(njet);
687}
688
452af8c7 689void AliPythia::Pyshow(Int_t ip1, Int_t ip2, Double_t qmax)
690{
691// Call Pythia jet reconstruction algorithm
692//
452af8c7 693 pyshow(ip1, ip2, qmax);
694}
695
696void AliPythia::Pyrobo(Int_t imi, Int_t ima, Double_t the, Double_t phi, Double_t bex, Double_t bey, Double_t bez)
697{
698 pyrobo(imi, ima, the, phi, bex, bey, bez);
699}
700
694b39f9 701void AliPythia::Pytune(Int_t itune)
702{
c5e2801a 703/*
704C
705C ITUNE NAME (detailed descriptions below)
706C 0 Default : No settings changed => linked Pythia version's defaults.
707C ====== Old UE, Q2-ordered showers ==========================================
708C 100 A : Rick Field's CDF Tune A
709C 101 AW : Rick Field's CDF Tune AW
710C 102 BW : Rick Field's CDF Tune BW
711C 103 DW : Rick Field's CDF Tune DW
712C 104 DWT : Rick Field's CDF Tune DW with slower UE energy scaling
713C 105 QW : Rick Field's CDF Tune QW (NB: needs CTEQ6.1M pdfs externally)
714C 106 ATLAS-DC2: Arthur Moraes' (old) ATLAS tune (ATLAS DC2 / Rome)
715C 107 ACR : Tune A modified with annealing CR
716C 108 D6 : Rick Field's CDF Tune D6 (NB: needs CTEQ6L pdfs externally)
717C 109 D6T : Rick Field's CDF Tune D6T (NB: needs CTEQ6L pdfs externally)
718C ====== Intermediate Models =================================================
719C 200 IM 1 : Intermediate model: new UE, Q2-ordered showers, annealing CR
720C 201 APT : Tune A modified to use pT-ordered final-state showers
721C ====== New UE, interleaved pT-ordered showers, annealing CR ================
722C 300 S0 : Sandhoff-Skands Tune 0
723C 301 S1 : Sandhoff-Skands Tune 1
724C 302 S2 : Sandhoff-Skands Tune 2
725C 303 S0A : S0 with "Tune A" UE energy scaling
726C 304 NOCR : New UE "best try" without colour reconnections
727C 305 Old : New UE, original (primitive) colour reconnections
728C 306 ATLAS-CSC: Arthur Moraes' (new) ATLAS tune (needs CTEQ6L externally)
729C ======= The Uppsala models =================================================
730C ( NB! must be run with special modified Pythia 6.215 version )
731C ( available from http://www.isv.uu.se/thep/MC/scigal/ )
732C 400 GAL 0 : Generalized area-law model. Old parameters
733C 401 SCI 0 : Soft-Colour-Interaction model. Old parameters
734C 402 GAL 1 : Generalized area-law model. Tevatron MB retuned (Skands)
735*/
694b39f9 736 pytune(itune);
737}
738
9b61ba2a 739void AliPythia::Py2ent(Int_t idx, Int_t pdg1, Int_t pdg2, Double_t p){
740 // Inset 2-parton system at line idx
741 py2ent(idx, pdg1, pdg2, p);
742}
452af8c7 743
744
32c8e463 745void AliPythia::InitQuenching(Float_t cMin, Float_t cMax, Float_t k, Int_t iECMethod, Float_t zmax, Int_t ngmax)
0f482ae4 746{
747// Initializes
748// (1) The quenching model using quenching weights according to C. Salgado and U. Wiedemann
749// (2) The nuclear geometry using the Glauber Model
750//
6b435cde 751
18b7a4a1 752 fGlauber = AliFastGlauber::Instance();
0f482ae4 753 fGlauber->Init(2);
754 fGlauber->SetCentralityClass(cMin, cMax);
755
756 fQuenchingWeights = new AliQuenchingWeights();
757 fQuenchingWeights->InitMult();
86b6ad68 758 fQuenchingWeights->SetK(k);
0f482ae4 759 fQuenchingWeights->SetECMethod(AliQuenchingWeights::kECMethod(iECMethod));
32c8e463 760 fNGmax = ngmax;
761 fZmax = zmax;
762
0f482ae4 763}
764
765
452af8c7 766void AliPythia::Quench()
767{
768//
769//
770// Simple Jet Quenching routine:
771// =============================
772// The jet formed by all final state partons radiated by the parton created
0f482ae4 773// in the hard collisions is quenched by a factor (1-z) using light cone variables in
774// the initial parton reference frame:
452af8c7 775// (E + p_z)new = (1-z) (E + p_z)old
776//
0f482ae4 777//
778//
779//
452af8c7 780// The lost momentum is first balanced by one gluon with virtuality > 0.
781// Subsequently the gluon splits to yield two gluons with E = p.
782//
0f482ae4 783//
784//
4e383037 785 static Float_t eMean = 0.;
786 static Int_t icall = 0;
0f482ae4 787
c2c598a3 788 Double_t p0[4][5];
789 Double_t p1[4][5];
790 Double_t p2[4][5];
791 Int_t klast[4] = {-1, -1, -1, -1};
452af8c7 792
793 Int_t numpart = fPyjets->N;
86b6ad68 794 Double_t px = 0., py = 0., pz = 0., e = 0., m = 0., p = 0., pt = 0., theta = 0., phi = 0.;
c2c598a3 795 Double_t pxq[4], pyq[4], pzq[4], eq[4], yq[4], mq[4], pq[4], phiq[4], thetaq[4], ptq[4];
796 Bool_t quenched[4];
b280c4cc 797 Double_t wjtKick[4];
c2c598a3 798 Int_t nGluon[4];
86b6ad68 799 Int_t qPdg[4];
0f482ae4 800 Int_t imo, kst, pdg;
b280c4cc 801
511db649 802//
c2c598a3 803// Sore information about Primary partons
804//
805// j =
806// 0, 1 partons from hard scattering
807// 2, 3 partons from initial state radiation
808//
809 for (Int_t i = 2; i <= 7; i++) {
810 Int_t j = 0;
811 // Skip gluons that participate in hard scattering
812 if (i == 4 || i == 5) continue;
813 // Gluons from hard Scattering
814 if (i == 6 || i == 7) {
815 j = i - 6;
816 pxq[j] = fPyjets->P[0][i];
817 pyq[j] = fPyjets->P[1][i];
818 pzq[j] = fPyjets->P[2][i];
819 eq[j] = fPyjets->P[3][i];
820 mq[j] = fPyjets->P[4][i];
821 } else {
822 // Gluons from initial state radiation
823 //
824 // Obtain 4-momentum vector from difference between original parton and parton after gluon
825 // radiation. Energy is calculated independently because initial state radition does not
826 // conserve strictly momentum and energy for each partonic system independently.
827 //
828 // Not very clean. Should be improved !
829 //
830 //
831 j = i;
832 pxq[j] = fPyjets->P[0][i] - fPyjets->P[0][i+2];
833 pyq[j] = fPyjets->P[1][i] - fPyjets->P[1][i+2];
834 pzq[j] = fPyjets->P[2][i] - fPyjets->P[2][i+2];
835 mq[j] = fPyjets->P[4][i];
836 eq[j] = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j] + pzq[j] * pzq[j] + mq[j] * mq[j]);
837 }
838//
839// Calculate some kinematic variables
511db649 840//
4e383037 841 yq[j] = 0.5 * TMath::Log((eq[j] + pzq[j] + 1.e-14) / (eq[j] - pzq[j] + 1.e-14));
0f482ae4 842 pq[j] = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j] + pzq[j] * pzq[j]);
843 phiq[j] = TMath::Pi()+TMath::ATan2(-pyq[j], -pxq[j]);
844 ptq[j] = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j]);
845 thetaq[j] = TMath::ATan2(ptq[j], pzq[j]);
86b6ad68 846 qPdg[j] = fPyjets->K[1][i];
847 }
848
849 Double_t int0[4];
850 Double_t int1[4];
86b6ad68 851
b280c4cc 852 fGlauber->GetI0I1ForPythiaAndXY(4, phiq, int0, int1, fXJet, fYJet, 15.);
853
86b6ad68 854 for (Int_t j = 0; j < 4; j++) {
c2c598a3 855 //
856 // Quench only central jets and with E > 10.
857 //
86b6ad68 858
859
860 Int_t itype = (qPdg[j] == 21) ? 2 : 1;
861 Double_t eloss = fQuenchingWeights->GetELossRandomKFast(itype, int0[j], int1[j], eq[j]);
862
c2c598a3 863 if (TMath::Abs(yq[j]) > 2.5 || eq[j] < 10.) {
b280c4cc 864 fZQuench[j] = 0.;
0f482ae4 865 } else {
c2c598a3 866 if (eq[j] > 40. && TMath::Abs(yq[j]) < 0.5) {
4e383037 867 icall ++;
868 eMean += eloss;
869 }
0f482ae4 870 //
871 // Extra pt
86b6ad68 872 Double_t l = fQuenchingWeights->CalcLk(int0[j], int1[j]);
873 wjtKick[j] = TMath::Sqrt(l * fQuenchingWeights->CalcQk(int0[j], int1[j]));
0f482ae4 874 //
875 // Fractional energy loss
b280c4cc 876 fZQuench[j] = eloss / eq[j];
0f482ae4 877 //
878 // Avoid complete loss
879 //
1044c4d8 880 if (fZQuench[j] > fZmax) fZQuench[j] = fZmax;
0f482ae4 881 //
882 // Some debug printing
86b6ad68 883
884
bf9bb016 885// 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",
886// j, itype, eq[j], phiq[j], l, eloss, wjtKick[j], eMean / Float_t(icall+1), yq[j]);
4e383037 887
b280c4cc 888// fZQuench[j] = 0.8;
889// while (fZQuench[j] >= 0.95) fZQuench[j] = gRandom->Exp(0.2);
0f482ae4 890 }
4e383037 891
b280c4cc 892 quenched[j] = (fZQuench[j] > 0.01);
4e383037 893 } // primary partons
c2c598a3 894
b280c4cc 895
896
6e90ad26 897 Double_t pNew[1000][4];
898 Int_t kNew[1000];
899 Int_t icount = 0;
b280c4cc 900 Double_t zquench[4];
901
6e90ad26 902//
4e383037 903// System Loop
c2c598a3 904 for (Int_t isys = 0; isys < 4; isys++) {
6e90ad26 905// Skip to next system if not quenched.
4e383037 906 if (!quenched[isys]) continue;
907
b280c4cc 908 nGluon[isys] = 1 + Int_t(fZQuench[isys] / (1. - fZQuench[isys]));
32c8e463 909 if (nGluon[isys] > fNGmax) nGluon[isys] = fNGmax;
b280c4cc 910 zquench[isys] = 1. - TMath::Power(1. - fZQuench[isys], 1./Double_t(nGluon[isys]));
4e383037 911 wjtKick[isys] = wjtKick[isys] / TMath::Sqrt(Double_t(nGluon[isys]));
0f482ae4 912
4e383037 913
914
915 Int_t igMin = -1;
916 Int_t igMax = -1;
917 Double_t pg[4] = {0., 0., 0., 0.};
918
919//
920// Loop on radiation events
921
922 for (Int_t iglu = 0; iglu < nGluon[isys]; iglu++) {
6e90ad26 923 while (1) {
924 icount = 0;
925 for (Int_t k = 0; k < 4; k++)
926 {
927 p0[isys][k] = 0.;
928 p1[isys][k] = 0.;
929 p2[isys][k] = 0.;
930 }
931// Loop over partons
932 for (Int_t i = 0; i < numpart; i++)
933 {
934 imo = fPyjets->K[2][i];
935 kst = fPyjets->K[0][i];
936 pdg = fPyjets->K[1][i];
937
938
939
0f482ae4 940// Quarks and gluons only
6e90ad26 941 if (pdg != 21 && TMath::Abs(pdg) > 6) continue;
0f482ae4 942// Particles from hard scattering only
c2c598a3 943
6e90ad26 944 if (imo > 8 && imo < 1000) imo = fPyjets->K[2][imo - 1];
c2c598a3 945 Int_t imom = imo % 1000;
946 if ((isys == 0 || isys == 1) && ((imom != (isys + 7)))) continue;
947 if ((isys == 2 || isys == 3) && ((imom != (isys + 1)))) continue;
948
6e90ad26 949
0f482ae4 950// Skip comment lines
6e90ad26 951 if (kst != 1 && kst != 2) continue;
0f482ae4 952//
953// Parton kinematic
6e90ad26 954 px = fPyjets->P[0][i];
955 py = fPyjets->P[1][i];
956 pz = fPyjets->P[2][i];
957 e = fPyjets->P[3][i];
958 m = fPyjets->P[4][i];
959 pt = TMath::Sqrt(px * px + py * py);
960 p = TMath::Sqrt(px * px + py * py + pz * pz);
961 phi = TMath::Pi() + TMath::ATan2(-py, -px);
962 theta = TMath::ATan2(pt, pz);
963
0f482ae4 964//
c2c598a3 965// Save 4-momentum sum for balancing
966 Int_t index = isys;
6e90ad26 967
968 p0[index][0] += px;
969 p0[index][1] += py;
970 p0[index][2] += pz;
971 p0[index][3] += e;
6e90ad26 972
973 klast[index] = i;
974
0f482ae4 975//
976// Fractional energy loss
b280c4cc 977 Double_t z = zquench[index];
4e383037 978
c2c598a3 979
4e383037 980// Don't fully quench radiated gluons
981//
982 if (imo > 1000) {
983// This small factor makes sure that the gluons are not too close in phase space to avoid recombination
984//
985
c2c598a3 986 z = 0.02;
4e383037 987 }
c2c598a3 988// printf("z: %d %f\n", imo, z);
989
4e383037 990
991//
6e90ad26 992
993 //
994 //
995 // Transform into frame in which initial parton is along z-axis
996 //
997 TVector3 v(px, py, pz);
998 v.RotateZ(-phiq[index]); v.RotateY(-thetaq[index]);
999 Double_t pxs = v.X(); Double_t pys = v.Y(); Double_t pl = v.Z();
1000
1001 Double_t jt = TMath::Sqrt(pxs * pxs + pys * pys);
1002 Double_t mt2 = jt * jt + m * m;
1003 Double_t zmax = 1.;
1004 //
1005 // Kinematic limit on z
1006 //
4e383037 1007 if (m > 0.) zmax = 1. - m / TMath::Sqrt(m * m + jt * jt);
6e90ad26 1008 //
1009 // Change light-cone kinematics rel. to initial parton
1010 //
1011 Double_t eppzOld = e + pl;
1012 Double_t empzOld = e - pl;
1013
1014 Double_t eppzNew = (1. - z) * eppzOld;
1015 Double_t empzNew = empzOld - mt2 * z / eppzOld;
1016 Double_t eNew = 0.5 * (eppzNew + empzNew);
1017 Double_t plNew = 0.5 * (eppzNew - empzNew);
1018
1019 Double_t jtNew;
1020 //
1021 // if mt very small (or sometimes even < 0 for numerical reasons) set it to 0
1022 Double_t mt2New = eppzNew * empzNew;
1023 if (mt2New < 1.e-8) mt2New = 0.;
4e383037 1024 if (z < zmax) {
1025 if (m * m > mt2New) {
1026 //
1027 // This should not happen
1028 //
1029 Fatal("Quench()", "This should never happen %e %e %e!", m, eppzNew, empzNew);
1030 jtNew = 0;
1031 } else {
1032 jtNew = TMath::Sqrt(mt2New - m * m);
1033 }
6e90ad26 1034 } else {
4e383037 1035 // If pT is to small (probably a leading massive particle) we scale only the energy
1036 // This can cause negative masses of the radiated gluon
1037 // Let's hope for the best ...
1038 jtNew = jt;
1039 eNew = TMath::Sqrt(plNew * plNew + mt2);
1040
6e90ad26 1041 }
6e90ad26 1042 //
1043 // Calculate new px, py
1044 //
b07be423 1045 Double_t pxNew = 0;
1046 Double_t pyNew = 0;
6e90ad26 1047
b07be423 1048 if (jt>0) {
6b118b3c 1049 pxNew = jtNew / jt * pxs;
1050 pyNew = jtNew / jt * pys;
b07be423 1051 }
6e90ad26 1052// Double_t dpx = pxs - pxNew;
1053// Double_t dpy = pys - pyNew;
1054// Double_t dpz = pl - plNew;
1055// Double_t de = e - eNew;
1056// Double_t dmass2 = de * de - dpx * dpx - dpy * dpy - dpz * dpz;
1057// printf("New mass (1) %e %e %e %e %e %e %e \n", dmass2, jt, jtNew, pl, plNew, e, eNew);
1058// printf("New mass (2) %e %e \n", pxNew, pyNew);
1059 //
1060 // Rotate back
1061 //
1062 TVector3 w(pxNew, pyNew, plNew);
1063 w.RotateY(thetaq[index]); w.RotateZ(phiq[index]);
1064 pxNew = w.X(); pyNew = w.Y(); plNew = w.Z();
1065
1066 p1[index][0] += pxNew;
1067 p1[index][1] += pyNew;
1068 p1[index][2] += plNew;
1069 p1[index][3] += eNew;
1070 //
1071 // Updated 4-momentum vectors
1072 //
1073 pNew[icount][0] = pxNew;
1074 pNew[icount][1] = pyNew;
1075 pNew[icount][2] = plNew;
1076 pNew[icount][3] = eNew;
1077 kNew[icount] = i;
1078 icount++;
1079 } // parton loop
0f482ae4 1080 //
6e90ad26 1081 // Check if there was phase-space for quenching
0f482ae4 1082 //
0f482ae4 1083
6e90ad26 1084 if (icount == 0) quenched[isys] = kFALSE;
1085 if (!quenched[isys]) break;
1086
1087 for (Int_t j = 0; j < 4; j++)
1088 {
1089 p2[isys][j] = p0[isys][j] - p1[isys][j];
1090 }
1091 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];
6e90ad26 1092 if (p2[isys][4] > 0.) {
1093 p2[isys][4] = TMath::Sqrt(p2[isys][4]);
1094 break;
1095 } else {
b280c4cc 1096 printf("Warning negative mass squared in system %d %f ! \n", isys, zquench[isys]);
4e383037 1097 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]);
6e90ad26 1098 if (p2[isys][4] < -0.01) {
4e383037 1099 printf("Negative mass squared !\n");
1100 // Here we have to put the gluon back to mass shell
1101 // This will lead to a small energy imbalance
1102 p2[isys][4] = 0.;
1103 p2[isys][3] = TMath::Sqrt(p2[isys][0] * p2[isys][0] + p2[isys][1] * p2[isys][1] + p2[isys][2] * p2[isys][2]);
1104 break;
6e90ad26 1105 } else {
1106 p2[isys][4] = 0.;
1107 break;
1108 }
1109 }
6e90ad26 1110 /*
6e90ad26 1111 zHeavy *= 0.98;
1112 printf("zHeavy lowered to %f\n", zHeavy);
1113 if (zHeavy < 0.01) {
1114 printf("No success ! \n");
1115 icount = 0;
1116 quenched[isys] = kFALSE;
1117 break;
1118 }
4e383037 1119 */
1120 } // iteration on z (while)
1121
6e90ad26 1122// Update event record
1123 for (Int_t k = 0; k < icount; k++) {
1124// 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] );
1125 fPyjets->P[0][kNew[k]] = pNew[k][0];
1126 fPyjets->P[1][kNew[k]] = pNew[k][1];
1127 fPyjets->P[2][kNew[k]] = pNew[k][2];
1128 fPyjets->P[3][kNew[k]] = pNew[k][3];
0f482ae4 1129 }
4e383037 1130 //
1131 // Add the gluons
1132 //
1133 Int_t ish = 0;
1837e95c 1134 Int_t iGlu;
4e383037 1135 if (!quenched[isys]) continue;
0f482ae4 1136//
1137// Last parton from shower i
4e383037 1138 Int_t in = klast[isys];
0f482ae4 1139//
1140// Continue if no parton in shower i selected
1141 if (in == -1) continue;
1142//
1143// If this is the second initial parton and it is behind the first move pointer by previous ish
4e383037 1144 if (isys == 1 && klast[1] > klast[0]) in += ish;
0f482ae4 1145//
1146// Starting index
452af8c7 1147
4e383037 1148// jmin = in - 1;
0f482ae4 1149// How many additional gluons will be generated
1150 ish = 1;
4e383037 1151 if (p2[isys][4] > 0.05) ish = 2;
0f482ae4 1152//
1153// Position of gluons
4e383037 1154 iGlu = numpart;
1155 if (iglu == 0) igMin = iGlu;
1156 igMax = iGlu;
0f482ae4 1157 numpart += ish;
1158 (fPyjets->N) += ish;
4e383037 1159
0f482ae4 1160 if (ish == 1) {
4e383037 1161 fPyjets->P[0][iGlu] = p2[isys][0];
1162 fPyjets->P[1][iGlu] = p2[isys][1];
1163 fPyjets->P[2][iGlu] = p2[isys][2];
1164 fPyjets->P[3][iGlu] = p2[isys][3];
1165 fPyjets->P[4][iGlu] = p2[isys][4];
0f482ae4 1166
4e383037 1167 fPyjets->K[0][iGlu] = 1;
1168 if (iglu == nGluon[isys] - 1) fPyjets->K[0][iGlu] = 1;
0f482ae4 1169 fPyjets->K[1][iGlu] = 21;
4e383037 1170 fPyjets->K[2][iGlu] = fPyjets->K[2][in] + 1000;
0f482ae4 1171 fPyjets->K[3][iGlu] = -1;
1172 fPyjets->K[4][iGlu] = -1;
4e383037 1173
1174 pg[0] += p2[isys][0];
1175 pg[1] += p2[isys][1];
1176 pg[2] += p2[isys][2];
1177 pg[3] += p2[isys][3];
0f482ae4 1178 } else {
1179 //
1180 // Split gluon in rest frame.
1181 //
4e383037 1182 Double_t bx = p2[isys][0] / p2[isys][3];
1183 Double_t by = p2[isys][1] / p2[isys][3];
1184 Double_t bz = p2[isys][2] / p2[isys][3];
1185 Double_t pst = p2[isys][4] / 2.;
0f482ae4 1186 //
1187 // Isotropic decay ????
1188 Double_t cost = 2. * gRandom->Rndm() - 1.;
60e55aee 1189 Double_t sint = TMath::Sqrt((1.-cost)*(1.+cost));
2ab330c9 1190 Double_t phis = 2. * TMath::Pi() * gRandom->Rndm();
0f482ae4 1191
1192 Double_t pz1 = pst * cost;
1193 Double_t pz2 = -pst * cost;
1194 Double_t pt1 = pst * sint;
1195 Double_t pt2 = -pst * sint;
2ab330c9 1196 Double_t px1 = pt1 * TMath::Cos(phis);
1197 Double_t py1 = pt1 * TMath::Sin(phis);
1198 Double_t px2 = pt2 * TMath::Cos(phis);
1199 Double_t py2 = pt2 * TMath::Sin(phis);
0f482ae4 1200
1201 fPyjets->P[0][iGlu] = px1;
1202 fPyjets->P[1][iGlu] = py1;
1203 fPyjets->P[2][iGlu] = pz1;
1204 fPyjets->P[3][iGlu] = pst;
1205 fPyjets->P[4][iGlu] = 0.;
1206
4e383037 1207 fPyjets->K[0][iGlu] = 1 ;
0f482ae4 1208 fPyjets->K[1][iGlu] = 21;
4e383037 1209 fPyjets->K[2][iGlu] = fPyjets->K[2][in] + 1000;
0f482ae4 1210 fPyjets->K[3][iGlu] = -1;
1211 fPyjets->K[4][iGlu] = -1;
1212
1213 fPyjets->P[0][iGlu+1] = px2;
1214 fPyjets->P[1][iGlu+1] = py2;
1215 fPyjets->P[2][iGlu+1] = pz2;
1216 fPyjets->P[3][iGlu+1] = pst;
1217 fPyjets->P[4][iGlu+1] = 0.;
1218
4e383037 1219 fPyjets->K[0][iGlu+1] = 1;
1220 if (iglu == nGluon[isys] - 1) fPyjets->K[0][iGlu+1] = 1;
0f482ae4 1221 fPyjets->K[1][iGlu+1] = 21;
4e383037 1222 fPyjets->K[2][iGlu+1] = fPyjets->K[2][in] + 1000;
0f482ae4 1223 fPyjets->K[3][iGlu+1] = -1;
1224 fPyjets->K[4][iGlu+1] = -1;
1225 SetMSTU(1,0);
1226 SetMSTU(2,0);
1227 //
1228 // Boost back
1229 //
1230 Pyrobo(iGlu + 1, iGlu + 2, 0., 0., bx, by, bz);
1231 }
4e383037 1232/*
1233 for (Int_t ig = iGlu; ig < iGlu+ish; ig++) {
1234 Double_t px, py, pz;
1235 px = fPyjets->P[0][ig];
1236 py = fPyjets->P[1][ig];
1237 pz = fPyjets->P[2][ig];
1238 TVector3 v(px, py, pz);
1239 v.RotateZ(-phiq[isys]);
1240 v.RotateY(-thetaq[isys]);
1241 Double_t pxs = v.X(); Double_t pys = v.Y(); Double_t pzs = v.Z();
1242 Double_t r = AliPythiaRndm::GetPythiaRandom()->Rndm();
1243 Double_t jtKick = 0.3 * TMath::Sqrt(-TMath::Log(r));
1244 if (ish == 2) jtKick = wjtKick[i] * TMath::Sqrt(-TMath::Log(r)) / TMath::Sqrt(2.);
1245 Double_t phiKick = 2. * TMath::Pi() * AliPythiaRndm::GetPythiaRandom()->Rndm();
1246 pxs += jtKick * TMath::Cos(phiKick);
1247 pys += jtKick * TMath::Sin(phiKick);
1248 TVector3 w(pxs, pys, pzs);
1249 w.RotateY(thetaq[isys]);
1250 w.RotateZ(phiq[isys]);
1251 fPyjets->P[0][ig] = w.X();
1252 fPyjets->P[1][ig] = w.Y();
1253 fPyjets->P[2][ig] = w.Z();
1254 fPyjets->P[2][ig] = w.Mag();
1255 }
1256*/
1257 } // kGluon
1258
6e90ad26 1259
4e383037 1260 // Check energy conservation
0f482ae4 1261 Double_t pxs = 0.;
1262 Double_t pys = 0.;
1263 Double_t pzs = 0.;
1264 Double_t es = 14000.;
1265
1266 for (Int_t i = 0; i < numpart; i++)
1267 {
1268 kst = fPyjets->K[0][i];
1269 if (kst != 1 && kst != 2) continue;
1270 pxs += fPyjets->P[0][i];
1271 pys += fPyjets->P[1][i];
1272 pzs += fPyjets->P[2][i];
1273 es -= fPyjets->P[3][i];
1274 }
1275 if (TMath::Abs(pxs) > 1.e-2 ||
1276 TMath::Abs(pys) > 1.e-2 ||
1277 TMath::Abs(pzs) > 1.e-1) {
1278 printf("%e %e %e %e\n", pxs, pys, pzs, es);
4e383037 1279// Fatal("Quench()", "4-Momentum non-conservation");
452af8c7 1280 }
4e383037 1281
1282 } // end quenching loop (systems)
6e90ad26 1283// Clean-up
0f482ae4 1284 for (Int_t i = 0; i < numpart; i++)
1285 {
4e383037 1286 imo = fPyjets->K[2][i];
1287 if (imo > 1000) {
1288 fPyjets->K[2][i] = fPyjets->K[2][i] % 1000;
1289 }
0f482ae4 1290 }
4e383037 1291// this->Pylist(1);
0f482ae4 1292} // end quench
90d7b703 1293
992f2843 1294
1295void AliPythia::Pyquen(Double_t a, Int_t ibf, Double_t b)
1296{
1297 // Igor Lokthine's quenching routine
12cb0bc0 1298 // http://lokhtin.web.cern.ch/lokhtin/pyquen/pyquen.txt
1299
992f2843 1300 pyquen(a, ibf, b);
1301}
b280c4cc 1302
12cb0bc0 1303void AliPythia::SetPyquenParameters(Double_t t0, Double_t tau0, Int_t nf, Int_t iengl, Int_t iangl)
1304{
1305 // Set the parameters for the PYQUEN package.
1306 // See comments in PyquenCommon.h
1307
1308
1309 PYQPAR.t0 = t0;
1310 PYQPAR.tau0 = tau0;
1311 PYQPAR.nf = nf;
1312 PYQPAR.iengl = iengl;
1313 PYQPAR.iangl = iangl;
1314}
1315
1316
16a82508 1317void AliPythia::Pyevnw()
1318{
1319 // New multiple interaction scenario
1320 pyevnw();
1321}
1322
cd07c39b 1323void AliPythia::Pyshowq(Int_t ip1, Int_t ip2, Double_t qmax)
1324{
1325 // Call medium-modified Pythia jet reconstruction algorithm
1326 //
1327 pyshowq(ip1, ip2, qmax);
1328}
1329
b280c4cc 1330void AliPythia::GetQuenchingParameters(Double_t& xp, Double_t& yp, Double_t z[4])
1331{
1332 // Return event specific quenching parameters
1333 xp = fXJet;
1334 yp = fYJet;
1335 for (Int_t i = 0; i < 4; i++) z[i] = fZQuench[i];
1336
1337}
1338
3dc3ec94 1339void AliPythia::ConfigHeavyFlavor()
1340{
1341 //
1342 // Default configuration for Heavy Flavor production
1343 //
1344 // All QCD processes
1345 //
1346 SetMSEL(1);
1347
1348 // No multiple interactions
1349 SetMSTP(81,0);
39c2e610 1350 SetPARP(81, 0.);
1351 SetPARP(82, 0.);
3dc3ec94 1352 // Initial/final parton shower on (Pythia default)
1353 SetMSTP(61,1);
1354 SetMSTP(71,1);
1355
1356 // 2nd order alpha_s
1357 SetMSTP(2,2);
1358
1359 // QCD scales
1360 SetMSTP(32,2);
1361 SetPARP(34,1.0);
1362}
e0e89f40 1363
1364void AliPythia::AtlasTuning()
1365{
1366 //
1367 // Configuration for the ATLAS tuning
e2de0ce1 1368 SetMSTP(51, AliStructFuncType::PDFsetIndex(kCTEQ5L)); // CTEQ5L pdf
e0e89f40 1369 SetMSTP(81,1); // Multiple Interactions ON
1370 SetMSTP(82,4); // Double Gaussian Model
39c2e610 1371 SetPARP(81,1.9); // Min. pt for multiple interactions (default in 6.2-14)
e0e89f40 1372 SetPARP(82,1.8); // [GeV] PT_min at Ref. energy
1373 SetPARP(89,1000.); // [GeV] Ref. energy
1374 SetPARP(90,0.16); // 2*epsilon (exponent in power law)
1375 SetPARP(83,0.5); // Core density in proton matter distribution (def.value)
1376 SetPARP(84,0.5); // Core radius
1377 SetPARP(85,0.33); // Regulates gluon prod. mechanism
1378 SetPARP(86,0.66); // Regulates gluon prod. mechanism
1379 SetPARP(67,1); // Regulates Initial State Radiation
1380}
e8a8adcd 1381
1382AliPythia& AliPythia::operator=(const AliPythia& rhs)
1383{
1384// Assignment operator
1385 rhs.Copy(*this);
1386 return *this;
1387}
1388
1389 void AliPythia::Copy(TObject&) const
1390{
1391 //
1392 // Copy
1393 //
1394 Fatal("Copy","Not implemented!\n");
1395}
cd07c39b 1396