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