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