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