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