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