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