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