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