]> 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
ESD track based on AliExternalTrackParam. Class redesign and clean-up (Yu.Belikov)
[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"
0f482ae4 20#include "../FASTSIM/AliFastGlauber.h"
21#include "../FASTSIM/AliQuenchingWeights.h"
22#include "TVector3.h"
8d2cd130 23
24ClassImp(AliPythia)
25
26#ifndef WIN32
27# define pyclus pyclus_
28# define pycell pycell_
452af8c7 29# define pyshow pyshow_
30# define pyrobo pyrobo_
992f2843 31# define pyquen pyquen_
8d2cd130 32# define type_of_call
33#else
34# define pyclus PYCLUS
35# define pycell PYCELL
452af8c7 36# define pyrobo PYROBO
992f2843 37# define pyquen PYQUEN
8d2cd130 38# define type_of_call _stdcall
39#endif
40
41extern "C" void type_of_call pyclus(Int_t & );
42extern "C" void type_of_call pycell(Int_t & );
452af8c7 43extern "C" void type_of_call pyshow(Int_t &, Int_t &, Double_t &);
44extern "C" void type_of_call pyrobo(Int_t &, Int_t &, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &);
992f2843 45extern "C" void type_of_call pyquen(Double_t &, Int_t &, Double_t &);
8d2cd130 46
47//_____________________________________________________________________________
48
49AliPythia* AliPythia::fgAliPythia=NULL;
50
51AliPythia::AliPythia()
52{
53// Default Constructor
54//
55// Set random number
7cdba479 56 if (!AliPythiaRndm::GetPythiaRandom())
57 AliPythiaRndm::SetPythiaRandom(GetRandom());
0f482ae4 58 fGlauber = 0;
59 fQuenchingWeights = 0;
8d2cd130 60}
61
62void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfunc)
63{
64// Initialise the process to generate
7cdba479 65 if (!AliPythiaRndm::GetPythiaRandom())
66 AliPythiaRndm::SetPythiaRandom(GetRandom());
8d2cd130 67
68 fProcess = process;
69 fEcms = energy;
70 fStrucFunc = strucfunc;
1d5b1b20 71//...Switch off decay of pi0, K0S, Lambda, Sigma+-, Xi0-, Omega-.
72 SetMDCY(Pycomp(111) ,1,0);
73 SetMDCY(Pycomp(310) ,1,0);
74 SetMDCY(Pycomp(3122),1,0);
75 SetMDCY(Pycomp(3112),1,0);
76 SetMDCY(Pycomp(3212),1,0);
77 SetMDCY(Pycomp(3222),1,0);
78 SetMDCY(Pycomp(3312),1,0);
79 SetMDCY(Pycomp(3322),1,0);
80 SetMDCY(Pycomp(3334),1,0);
81 // select structure function
8d2cd130 82 SetMSTP(52,2);
83 SetMSTP(51,strucfunc);
84//
85// Pythia initialisation for selected processes//
86//
87// Make MSEL clean
88//
89 for (Int_t i=1; i<= 200; i++) {
90 SetMSUB(i,0);
91 }
92// select charm production
93 switch (process)
94 {
65f2626c 95 case kPyOldUEQ2ordered: //Old underlying events with Q2 ordered QCD processes
96// Multiple interactions on.
97 SetMSTP(81,1);
98// Double Gaussian matter distribution.
99 SetMSTP(82,4);
100 SetPARP(83,0.5);
101 SetPARP(84,0.4);
102// pT0.
103 SetPARP(82,2.0);
104// Reference energy for pT0 and energy rescaling pace.
105 SetPARP(89,1800);
106 SetPARP(90,0.25);
107// String drawing almost completely minimizes string length.
108 SetPARP(85,0.9);
109 SetPARP(86,0.95);
110// ISR and FSR activity.
111 SetPARP(67,4);
112 SetPARP(71,4);
113// Lambda_FSR scale.
114 SetPARJ(81,0.29);
115 break;
116 case kPyOldUEQ2ordered2:
117// Old underlying events with Q2 ordered QCD processes
118// Multiple interactions on.
119 SetMSTP(81,1);
120// Double Gaussian matter distribution.
121 SetMSTP(82,4);
122 SetPARP(83,0.5);
123 SetPARP(84,0.4);
124// pT0.
125 SetPARP(82,2.0);
126// Reference energy for pT0 and energy rescaling pace.
127 SetPARP(89,1800);
128 SetPARP(90,0.16); // here is the difference with kPyOldUEQ2ordered
129// String drawing almost completely minimizes string length.
130 SetPARP(85,0.9);
131 SetPARP(86,0.95);
132// ISR and FSR activity.
133 SetPARP(67,4);
134 SetPARP(71,4);
135// Lambda_FSR scale.
136 SetPARJ(81,0.29);
137 break;
138 case kPyOldPopcorn:
139// Old production mechanism: Old Popcorn
140 SetMSEL(1);
141 SetMSTJ(12,3);
142// (D=2) Like MSTJ(12)=2 but added prod ofthe 1er rank baryon
143 SetMSTP(88,2);
144// (D=1)see can be used to form baryons (BARYON JUNCTION)
145 SetMSTJ(1,1);
146 SetMSTP(51,kCTEQ5L);// CTEQ 5L ! CTEQ5L pdf
147 SetMSTP(81,1); // Multiple Interactions ON
148 SetMSTP(82,4); // Double Gaussian Model
149 SetPARP(82,1.8); // [GeV] PT_min at Ref. energy
150 SetPARP(89,1000.); // [GeV] Ref. energy
151 SetPARP(90,0.16); // 2*epsilon (exponent in power law)
152 SetPARP(83,0.5); // Core density in proton matter dist. (def.value)
153 SetPARP(84,0.5); // Core radius
154 SetPARP(85,0.33); // Regulates gluon prod. mechanism
155 SetPARP(86,0.66); // Regulates gluon prod. mechanism
156 SetPARP(67,1); // Regulate gluon prod. mechanism
157 break;
8d2cd130 158 case kPyCharm:
159 SetMSEL(4);
8d2cd130 160// heavy quark masses
161
162 SetPMAS(4,1,1.2);
163 SetMSTU(16,2);
164//
165// primordial pT
166 SetMSTP(91,1);
167 SetPARP(91,1.);
168 SetPARP(93,5.);
169//
170 break;
171 case kPyBeauty:
172 SetMSEL(5);
173 SetPMAS(5,1,4.75);
174 SetMSTU(16,2);
175 break;
176 case kPyJpsi:
177 SetMSEL(0);
178// gg->J/Psi g
179 SetMSUB(86,1);
180 break;
181 case kPyJpsiChi:
182 SetMSEL(0);
183// gg->J/Psi g
184 SetMSUB(86,1);
185// gg-> chi_0c g
186 SetMSUB(87,1);
187// gg-> chi_1c g
188 SetMSUB(88,1);
189// gg-> chi_2c g
190 SetMSUB(89,1);
191 break;
192 case kPyCharmUnforced:
193 SetMSEL(0);
194// gq->qg
195 SetMSUB(28,1);
196// gg->qq
197 SetMSUB(53,1);
198// gg->gg
199 SetMSUB(68,1);
200 break;
201 case kPyBeautyUnforced:
202 SetMSEL(0);
203// gq->qg
204 SetMSUB(28,1);
205// gg->qq
206 SetMSUB(53,1);
207// gg->gg
208 SetMSUB(68,1);
209 break;
210 case kPyMb:
211// Minimum Bias pp-Collisions
212//
213//
214// select Pythia min. bias model
215 SetMSEL(0);
511db649 216 SetMSUB(92,1); // single diffraction AB-->XB
217 SetMSUB(93,1); // single diffraction AB-->AX
218 SetMSUB(94,1); // double diffraction
219 SetMSUB(95,1); // low pt production
220
221//
222// ATLAS Tuning
223//
c01c4118 224
225 SetMSTP(51, kCTEQ5L); // CTEQ5L pdf
511db649 226 SetMSTP(81,1); // Multiple Interactions ON
227 SetMSTP(82,4); // Double Gaussian Model
228
229 SetPARP(82,1.8); // [GeV] PT_min at Ref. energy
230 SetPARP(89,1000.); // [GeV] Ref. energy
231 SetPARP(90,0.16); // 2*epsilon (exponent in power law)
232 SetPARP(83,0.5); // Core density in proton matter distribution (def.value)
233 SetPARP(84,0.5); // Core radius
234 SetPARP(85,0.33); // Regulates gluon prod. mechanism
235 SetPARP(86,0.66); // Regulates gluon prod. mechanism
236 SetPARP(67,1); // Regulates Initial State Radiation
237 break;
8d2cd130 238 case kPyMbNonDiffr:
239// Minimum Bias pp-Collisions
240//
241//
242// select Pythia min. bias model
243 SetMSEL(0);
511db649 244 SetMSUB(95,1); // low pt production
0f482ae4 245
246//
247// ATLAS Tuning
248//
511db649 249
d682afd1 250 SetMSTP(51,kCTEQ5L); // CTEQ5L pdf
511db649 251 SetMSTP(81,1); // Multiple Interactions ON
252 SetMSTP(82,4); // Double Gaussian Model
253
254 SetPARP(82,1.8); // [GeV] PT_min at Ref. energy
255 SetPARP(89,1000.); // [GeV] Ref. energy
256 SetPARP(90,0.16); // 2*epsilon (exponent in power law)
257 SetPARP(83,0.5); // Core density in proton matter distribution (def.value)
258 SetPARP(84,0.5); // Core radius
259 SetPARP(85,0.33); // Regulates gluon prod. mechanism
260 SetPARP(86,0.66); // Regulates gluon prod. mechanism
261 SetPARP(67,1); // Regulates Initial State Radiation
8d2cd130 262 break;
263 case kPyJets:
264//
265// QCD Jets
266//
267 SetMSEL(1);
65f2626c 268 // Pythia Tune A (CDF)
269 //
270 SetPARP(67,4.); // Regulates Initial State Radiation
271 SetMSTP(82,4); // Double Gaussian Model
272 SetPARP(82,2.0); // [GeV] PT_min at Ref. energy
273 SetPARP(84,0.4); // Core radius
274 SetPARP(85,0.90) ; // Regulates gluon prod. mechanism
275 SetPARP(86,0.95); // Regulates gluon prod. mechanism
276 SetPARP(89,1800.); // [GeV] Ref. energy
277 SetPARP(90,0.25); // 2*epsilon (exponent in power law)
278 break;
8d2cd130 279 case kPyDirectGamma:
280 SetMSEL(10);
281 break;
adf4d898 282 case kPyCharmPbPbMNR:
283 case kPyD0PbPbMNR:
8d2cd130 284 // Tuning of Pythia parameters aimed to get a resonable agreement
285 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
286 // c-cbar single inclusive and double differential distributions.
287 // This parameter settings are meant to work with Pb-Pb collisions
adf4d898 288 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
8d2cd130 289 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
290 // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
291
292 // All QCD processes
293 SetMSEL(1);
294
295 // No multiple interactions
296 SetMSTP(81,0);
297 SetPARP(81,0.0);
298 SetPARP(82,0.0);
299
300 // Initial/final parton shower on (Pythia default)
301 SetMSTP(61,1);
302 SetMSTP(71,1);
303
304 // 2nd order alpha_s
305 SetMSTP(2,2);
306
307 // QCD scales
308 SetMSTP(32,2);
309 SetPARP(34,1.0);
310
adf4d898 311 // Intrinsic <kT>
8d2cd130 312 SetMSTP(91,1);
313 SetPARP(91,1.304);
314 SetPARP(93,6.52);
315
316 // Set c-quark mass
317 SetPMAS(4,1,1.2);
318
90d7b703 319 break;
320 case kPyDPlusPbPbMNR:
321 // Tuning of Pythia parameters aimed to get a resonable agreement
322 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
323 // c-cbar single inclusive and double differential distributions.
324 // This parameter settings are meant to work with Pb-Pb collisions
325 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
326 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
327 // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
328
329 // All QCD processes
330 SetMSEL(1);
331
332 // No multiple interactions
333 SetMSTP(81,0);
334 SetPARP(81,0.0);
335 SetPARP(82,0.0);
336
337 // Initial/final parton shower on (Pythia default)
338 SetMSTP(61,1);
339 SetMSTP(71,1);
340
341 // 2nd order alpha_s
342 SetMSTP(2,2);
343
344 // QCD scales
345 SetMSTP(32,2);
346 SetPARP(34,1.0);
347
348 // Intrinsic <kT>
349 SetMSTP(91,1);
350 SetPARP(91,1.304);
351 SetPARP(93,6.52);
352
353 // Set c-quark mass
354 SetPMAS(4,1,1.2);
355
8d2cd130 356 break;
adf4d898 357 case kPyCharmpPbMNR:
358 case kPyD0pPbMNR:
359 // Tuning of Pythia parameters aimed to get a resonable agreement
360 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
361 // c-cbar single inclusive and double differential distributions.
362 // This parameter settings are meant to work with p-Pb collisions
363 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
364 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
365 // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
366
367 // All QCD processes
368 SetMSEL(1);
369
370 // No multiple interactions
371 SetMSTP(81,0);
372 SetPARP(81,0.0);
373 SetPARP(82,0.0);
374
375 // Initial/final parton shower on (Pythia default)
376 SetMSTP(61,1);
377 SetMSTP(71,1);
378
379 // 2nd order alpha_s
380 SetMSTP(2,2);
381
382 // QCD scales
383 SetMSTP(32,2);
384 SetPARP(34,1.0);
385
386 // Intrinsic <kT>
387 SetMSTP(91,1);
388 SetPARP(91,1.16);
389 SetPARP(93,5.8);
390
391 // Set c-quark mass
392 SetPMAS(4,1,1.2);
393
90d7b703 394 break;
395 case kPyDPluspPbMNR:
396 // Tuning of Pythia parameters aimed to get a resonable agreement
397 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
398 // c-cbar single inclusive and double differential distributions.
399 // This parameter settings are meant to work with p-Pb collisions
400 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
401 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
402 // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
403
404 // All QCD processes
405 SetMSEL(1);
406
407 // No multiple interactions
408 SetMSTP(81,0);
409 SetPARP(81,0.0);
410 SetPARP(82,0.0);
411
412 // Initial/final parton shower on (Pythia default)
413 SetMSTP(61,1);
414 SetMSTP(71,1);
415
416 // 2nd order alpha_s
417 SetMSTP(2,2);
418
419 // QCD scales
420 SetMSTP(32,2);
421 SetPARP(34,1.0);
422
423 // Intrinsic <kT>
424 SetMSTP(91,1);
425 SetPARP(91,1.16);
426 SetPARP(93,5.8);
427
428 // Set c-quark mass
429 SetPMAS(4,1,1.2);
430
adf4d898 431 break;
432 case kPyCharmppMNR:
433 case kPyD0ppMNR:
434 // Tuning of Pythia parameters aimed to get a resonable agreement
435 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
436 // c-cbar single inclusive and double differential distributions.
437 // This parameter settings are meant to work with pp collisions
438 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
439 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
440 // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
441
442 // All QCD processes
443 SetMSEL(1);
444
445 // No multiple interactions
446 SetMSTP(81,0);
447 SetPARP(81,0.0);
448 SetPARP(82,0.0);
449
450 // Initial/final parton shower on (Pythia default)
451 SetMSTP(61,1);
452 SetMSTP(71,1);
453
454 // 2nd order alpha_s
455 SetMSTP(2,2);
456
457 // QCD scales
458 SetMSTP(32,2);
459 SetPARP(34,1.0);
460
461 // Intrinsic <kT^2>
462 SetMSTP(91,1);
463 SetPARP(91,1.);
464 SetPARP(93,5.);
465
466 // Set c-quark mass
467 SetPMAS(4,1,1.2);
468
90d7b703 469 break;
470 case kPyDPlusppMNR:
471 // Tuning of Pythia parameters aimed to get a resonable agreement
472 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
473 // c-cbar single inclusive and double differential distributions.
474 // This parameter settings are meant to work with pp collisions
475 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
476 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
477 // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
478
479 // All QCD processes
480 SetMSEL(1);
481
482 // No multiple interactions
483 SetMSTP(81,0);
484 SetPARP(81,0.0);
485 SetPARP(82,0.0);
486
487 // Initial/final parton shower on (Pythia default)
488 SetMSTP(61,1);
489 SetMSTP(71,1);
490
491 // 2nd order alpha_s
492 SetMSTP(2,2);
493
494 // QCD scales
495 SetMSTP(32,2);
496 SetPARP(34,1.0);
497
498 // Intrinsic <kT^2>
499 SetMSTP(91,1);
500 SetPARP(91,1.);
501 SetPARP(93,5.);
502
503 // Set c-quark mass
504 SetPMAS(4,1,1.2);
505
adf4d898 506 break;
507 case kPyBeautyPbPbMNR:
8d2cd130 508 // Tuning of Pythia parameters aimed to get a resonable agreement
509 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
510 // b-bbar single inclusive and double differential distributions.
511 // This parameter settings are meant to work with Pb-Pb collisions
512 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
513 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
514 // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
515
516 // All QCD processes
517 SetMSEL(1);
518
519 // No multiple interactions
520 SetMSTP(81,0);
521 SetPARP(81,0.0);
522 SetPARP(82,0.0);
523
524 // Initial/final parton shower on (Pythia default)
525 SetMSTP(61,1);
526 SetMSTP(71,1);
527
528 // 2nd order alpha_s
529 SetMSTP(2,2);
530
531 // QCD scales
532 SetMSTP(32,2);
533 SetPARP(34,1.0);
534 SetPARP(67,1.0);
535 SetPARP(71,1.0);
536
adf4d898 537 // Intrinsic <kT>
8d2cd130 538 SetMSTP(91,1);
539 SetPARP(91,2.035);
540 SetPARP(93,10.17);
541
542 // Set b-quark mass
543 SetPMAS(5,1,4.75);
544
adf4d898 545 break;
546 case kPyBeautypPbMNR:
547 // Tuning of Pythia parameters aimed to get a resonable agreement
548 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
549 // b-bbar single inclusive and double differential distributions.
550 // This parameter settings are meant to work with p-Pb collisions
551 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
552 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
553 // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
554
555 // All QCD processes
556 SetMSEL(1);
557
558 // No multiple interactions
559 SetMSTP(81,0);
560 SetPARP(81,0.0);
561 SetPARP(82,0.0);
562
563 // Initial/final parton shower on (Pythia default)
564 SetMSTP(61,1);
565 SetMSTP(71,1);
566
567 // 2nd order alpha_s
568 SetMSTP(2,2);
569
570 // QCD scales
571 SetMSTP(32,2);
572 SetPARP(34,1.0);
573 SetPARP(67,1.0);
574 SetPARP(71,1.0);
575
576 // Intrinsic <kT>
577 SetMSTP(91,1);
578 SetPARP(91,1.60);
579 SetPARP(93,8.00);
580
581 // Set b-quark mass
582 SetPMAS(5,1,4.75);
583
584 break;
585 case kPyBeautyppMNR:
586 // Tuning of Pythia parameters aimed to get a resonable agreement
587 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
588 // b-bbar single inclusive and double differential distributions.
589 // This parameter settings are meant to work with pp collisions
590 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
591 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
592 // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
593
594 // All QCD processes
595 SetMSEL(1);
596
597 // No multiple interactions
598 SetMSTP(81,0);
599 SetPARP(81,0.0);
600 SetPARP(82,0.0);
601
602 // Initial/final parton shower on (Pythia default)
603 SetMSTP(61,1);
604 SetMSTP(71,1);
605
606 // 2nd order alpha_s
607 SetMSTP(2,2);
608
609 // QCD scales
610 SetMSTP(32,2);
611 SetPARP(34,1.0);
612 SetPARP(67,1.0);
613 SetPARP(71,1.0);
614
615 // Intrinsic <kT>
616 SetMSTP(91,1);
617 SetPARP(91,1.);
618 SetPARP(93,5.);
619
620 // Set b-quark mass
621 SetPMAS(5,1,4.75);
622
8d2cd130 623 break;
589380c6 624
625 case kPyW:
626
627 //Inclusive production of W+/-
628 SetMSEL(0);
629 //f fbar -> W+
630 SetMSUB(2,1);
631 // //f fbar -> g W+
632 // SetMSUB(16,1);
633 // //f fbar -> gamma W+
634 // SetMSUB(20,1);
635 // //f g -> f W+
636 // SetMSUB(31,1);
637 // //f gamma -> f W+
638 // SetMSUB(36,1);
639
640 // Initial/final parton shower on (Pythia default)
641 // With parton showers on we are generating "W inclusive process"
642 SetMSTP(61,1); //Initial QCD & QED showers on
643 SetMSTP(71,1); //Final QCD & QED showers on
644
645 break;
8d2cd130 646 }
647//
648// Initialize PYTHIA
649 SetMSTP(41,1); // all resonance decays switched on
650
651 Initialize("CMS","p","p",fEcms);
652
653}
654
655Int_t AliPythia::CheckedLuComp(Int_t kf)
656{
657// Check Lund particle code (for debugging)
658 Int_t kc=Pycomp(kf);
659 printf("\n Lucomp kf,kc %d %d",kf,kc);
660 return kc;
661}
662
663void AliPythia::SetNuclei(Int_t a1, Int_t a2)
664{
665// Treat protons as inside nuclei with mass numbers a1 and a2
666// The MSTP array in the PYPARS common block is used to enable and
667// select the nuclear structure functions.
668// MSTP(52) : (D=1) choice of proton and nuclear structure-function library
669// =1: internal PYTHIA acording to MSTP(51)
670// =2: PDFLIB proton s.f., with MSTP(51) = 1000xNGROUP+NSET
671// If the following mass number both not equal zero, nuclear corrections of the stf are used.
672// MSTP(192) : Mass number of nucleus side 1
673// MSTP(193) : Mass number of nucleus side 2
674 SetMSTP(52,2);
675 SetMSTP(192, a1);
676 SetMSTP(193, a2);
677}
678
679
680AliPythia* AliPythia::Instance()
681{
682// Set random number generator
683 if (fgAliPythia) {
684 return fgAliPythia;
685 } else {
686 fgAliPythia = new AliPythia();
687 return fgAliPythia;
688 }
689}
690
691void AliPythia::PrintParticles()
692{
693// Print list of particl properties
694 Int_t np = 0;
c31f1d37 695 char* name = new char[16];
8d2cd130 696 for (Int_t kf=0; kf<1000000; kf++) {
697 for (Int_t c = 1; c > -2; c-=2) {
8d2cd130 698 Int_t kc = Pycomp(c*kf);
699 if (kc) {
700 Float_t mass = GetPMAS(kc,1);
701 Float_t width = GetPMAS(kc,2);
702 Float_t tau = GetPMAS(kc,4);
c31f1d37 703
8d2cd130 704 Pyname(kf,name);
705
706 np++;
707
708 printf("\n mass, width, tau: %6d %s %10.3f %10.3e %10.3e",
709 c*kf, name, mass, width, tau);
710 }
711 }
712 }
713 printf("\n Number of particles %d \n \n", np);
714}
715
716void AliPythia::ResetDecayTable()
717{
718// Set default values for pythia decay switches
719 Int_t i;
720 for (i = 1; i < 501; i++) SetMDCY(i,1,fDefMDCY[i]);
721 for (i = 1; i < 2001; i++) SetMDME(i,1,fDefMDME[i]);
722}
723
724void AliPythia::SetDecayTable()
725{
726// Set default values for pythia decay switches
727//
728 Int_t i;
729 for (i = 1; i < 501; i++) fDefMDCY[i] = GetMDCY(i,1);
730 for (i = 1; i < 2001; i++) fDefMDME[i] = GetMDME(i,1);
731}
732
733void AliPythia::Pyclus(Int_t& njet)
734{
735// Call Pythia clustering algorithm
736//
737 pyclus(njet);
738}
739
740void AliPythia::Pycell(Int_t& njet)
741{
742// Call Pythia jet reconstruction algorithm
743//
744 pycell(njet);
745}
746
452af8c7 747void AliPythia::Pyshow(Int_t ip1, Int_t ip2, Double_t qmax)
748{
749// Call Pythia jet reconstruction algorithm
750//
452af8c7 751 pyshow(ip1, ip2, qmax);
752}
753
754void AliPythia::Pyrobo(Int_t imi, Int_t ima, Double_t the, Double_t phi, Double_t bex, Double_t bey, Double_t bez)
755{
756 pyrobo(imi, ima, the, phi, bex, bey, bez);
757}
758
759
760
86b6ad68 761void AliPythia::InitQuenching(Float_t cMin, Float_t cMax, Float_t k, Int_t iECMethod)
0f482ae4 762{
763// Initializes
764// (1) The quenching model using quenching weights according to C. Salgado and U. Wiedemann
765// (2) The nuclear geometry using the Glauber Model
766//
767
768
769 fGlauber = new AliFastGlauber();
770 fGlauber->Init(2);
771 fGlauber->SetCentralityClass(cMin, cMax);
772
773 fQuenchingWeights = new AliQuenchingWeights();
774 fQuenchingWeights->InitMult();
86b6ad68 775 fQuenchingWeights->SetK(k);
0f482ae4 776 fQuenchingWeights->SetECMethod(AliQuenchingWeights::kECMethod(iECMethod));
0f482ae4 777}
778
779
452af8c7 780void AliPythia::Quench()
781{
782//
783//
784// Simple Jet Quenching routine:
785// =============================
786// The jet formed by all final state partons radiated by the parton created
0f482ae4 787// in the hard collisions is quenched by a factor (1-z) using light cone variables in
788// the initial parton reference frame:
452af8c7 789// (E + p_z)new = (1-z) (E + p_z)old
790//
0f482ae4 791//
792//
793//
452af8c7 794// The lost momentum is first balanced by one gluon with virtuality > 0.
795// Subsequently the gluon splits to yield two gluons with E = p.
796//
0f482ae4 797//
798//
4e383037 799 static Float_t eMean = 0.;
800 static Int_t icall = 0;
0f482ae4 801
c2c598a3 802 Double_t p0[4][5];
803 Double_t p1[4][5];
804 Double_t p2[4][5];
805 Int_t klast[4] = {-1, -1, -1, -1};
452af8c7 806
807 Int_t numpart = fPyjets->N;
86b6ad68 808 Double_t px = 0., py = 0., pz = 0., e = 0., m = 0., p = 0., pt = 0., theta = 0., phi = 0.;
c2c598a3 809 Double_t pxq[4], pyq[4], pzq[4], eq[4], yq[4], mq[4], pq[4], phiq[4], thetaq[4], ptq[4];
810 Bool_t quenched[4];
b280c4cc 811 Double_t wjtKick[4];
c2c598a3 812 Int_t nGluon[4];
86b6ad68 813 Int_t qPdg[4];
0f482ae4 814 Int_t imo, kst, pdg;
b280c4cc 815
511db649 816//
c2c598a3 817// Sore information about Primary partons
818//
819// j =
820// 0, 1 partons from hard scattering
821// 2, 3 partons from initial state radiation
822//
823 for (Int_t i = 2; i <= 7; i++) {
824 Int_t j = 0;
825 // Skip gluons that participate in hard scattering
826 if (i == 4 || i == 5) continue;
827 // Gluons from hard Scattering
828 if (i == 6 || i == 7) {
829 j = i - 6;
830 pxq[j] = fPyjets->P[0][i];
831 pyq[j] = fPyjets->P[1][i];
832 pzq[j] = fPyjets->P[2][i];
833 eq[j] = fPyjets->P[3][i];
834 mq[j] = fPyjets->P[4][i];
835 } else {
836 // Gluons from initial state radiation
837 //
838 // Obtain 4-momentum vector from difference between original parton and parton after gluon
839 // radiation. Energy is calculated independently because initial state radition does not
840 // conserve strictly momentum and energy for each partonic system independently.
841 //
842 // Not very clean. Should be improved !
843 //
844 //
845 j = i;
846 pxq[j] = fPyjets->P[0][i] - fPyjets->P[0][i+2];
847 pyq[j] = fPyjets->P[1][i] - fPyjets->P[1][i+2];
848 pzq[j] = fPyjets->P[2][i] - fPyjets->P[2][i+2];
849 mq[j] = fPyjets->P[4][i];
850 eq[j] = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j] + pzq[j] * pzq[j] + mq[j] * mq[j]);
851 }
852//
853// Calculate some kinematic variables
511db649 854//
4e383037 855 yq[j] = 0.5 * TMath::Log((eq[j] + pzq[j] + 1.e-14) / (eq[j] - pzq[j] + 1.e-14));
0f482ae4 856 pq[j] = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j] + pzq[j] * pzq[j]);
857 phiq[j] = TMath::Pi()+TMath::ATan2(-pyq[j], -pxq[j]);
858 ptq[j] = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j]);
859 thetaq[j] = TMath::ATan2(ptq[j], pzq[j]);
86b6ad68 860 qPdg[j] = fPyjets->K[1][i];
861 }
862
863 Double_t int0[4];
864 Double_t int1[4];
86b6ad68 865
b280c4cc 866 fGlauber->GetI0I1ForPythiaAndXY(4, phiq, int0, int1, fXJet, fYJet, 15.);
867
86b6ad68 868 for (Int_t j = 0; j < 4; j++) {
c2c598a3 869 //
870 // Quench only central jets and with E > 10.
871 //
86b6ad68 872
873
874 Int_t itype = (qPdg[j] == 21) ? 2 : 1;
875 Double_t eloss = fQuenchingWeights->GetELossRandomKFast(itype, int0[j], int1[j], eq[j]);
876
c2c598a3 877 if (TMath::Abs(yq[j]) > 2.5 || eq[j] < 10.) {
b280c4cc 878 fZQuench[j] = 0.;
0f482ae4 879 } else {
c2c598a3 880 if (eq[j] > 40. && TMath::Abs(yq[j]) < 0.5) {
4e383037 881 icall ++;
882 eMean += eloss;
883 }
0f482ae4 884 //
885 // Extra pt
86b6ad68 886 Double_t l = fQuenchingWeights->CalcLk(int0[j], int1[j]);
887 wjtKick[j] = TMath::Sqrt(l * fQuenchingWeights->CalcQk(int0[j], int1[j]));
0f482ae4 888 //
889 // Fractional energy loss
b280c4cc 890 fZQuench[j] = eloss / eq[j];
0f482ae4 891 //
892 // Avoid complete loss
893 //
b280c4cc 894 if (fZQuench[j] == 1.) fZQuench[j] = 0.95;
0f482ae4 895 //
896 // Some debug printing
86b6ad68 897
898
bf9bb016 899// 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",
900// j, itype, eq[j], phiq[j], l, eloss, wjtKick[j], eMean / Float_t(icall+1), yq[j]);
4e383037 901
b280c4cc 902// fZQuench[j] = 0.8;
903// while (fZQuench[j] >= 0.95) fZQuench[j] = gRandom->Exp(0.2);
0f482ae4 904 }
4e383037 905
b280c4cc 906 quenched[j] = (fZQuench[j] > 0.01);
4e383037 907 } // primary partons
c2c598a3 908
b280c4cc 909
910
6e90ad26 911 Double_t pNew[1000][4];
912 Int_t kNew[1000];
913 Int_t icount = 0;
b280c4cc 914 Double_t zquench[4];
915
6e90ad26 916//
4e383037 917// System Loop
c2c598a3 918 for (Int_t isys = 0; isys < 4; isys++) {
6e90ad26 919// Skip to next system if not quenched.
4e383037 920 if (!quenched[isys]) continue;
921
b280c4cc 922 nGluon[isys] = 1 + Int_t(fZQuench[isys] / (1. - fZQuench[isys]));
4e383037 923 if (nGluon[isys] > 6) nGluon[isys] = 6;
b280c4cc 924 zquench[isys] = 1. - TMath::Power(1. - fZQuench[isys], 1./Double_t(nGluon[isys]));
4e383037 925 wjtKick[isys] = wjtKick[isys] / TMath::Sqrt(Double_t(nGluon[isys]));
0f482ae4 926
4e383037 927
928
929 Int_t igMin = -1;
930 Int_t igMax = -1;
931 Double_t pg[4] = {0., 0., 0., 0.};
932
933//
934// Loop on radiation events
935
936 for (Int_t iglu = 0; iglu < nGluon[isys]; iglu++) {
6e90ad26 937 while (1) {
938 icount = 0;
939 for (Int_t k = 0; k < 4; k++)
940 {
941 p0[isys][k] = 0.;
942 p1[isys][k] = 0.;
943 p2[isys][k] = 0.;
944 }
945// Loop over partons
946 for (Int_t i = 0; i < numpart; i++)
947 {
948 imo = fPyjets->K[2][i];
949 kst = fPyjets->K[0][i];
950 pdg = fPyjets->K[1][i];
951
952
953
0f482ae4 954// Quarks and gluons only
6e90ad26 955 if (pdg != 21 && TMath::Abs(pdg) > 6) continue;
0f482ae4 956// Particles from hard scattering only
c2c598a3 957
6e90ad26 958 if (imo > 8 && imo < 1000) imo = fPyjets->K[2][imo - 1];
c2c598a3 959 Int_t imom = imo % 1000;
960 if ((isys == 0 || isys == 1) && ((imom != (isys + 7)))) continue;
961 if ((isys == 2 || isys == 3) && ((imom != (isys + 1)))) continue;
962
6e90ad26 963
0f482ae4 964// Skip comment lines
6e90ad26 965 if (kst != 1 && kst != 2) continue;
0f482ae4 966//
967// Parton kinematic
6e90ad26 968 px = fPyjets->P[0][i];
969 py = fPyjets->P[1][i];
970 pz = fPyjets->P[2][i];
971 e = fPyjets->P[3][i];
972 m = fPyjets->P[4][i];
973 pt = TMath::Sqrt(px * px + py * py);
974 p = TMath::Sqrt(px * px + py * py + pz * pz);
975 phi = TMath::Pi() + TMath::ATan2(-py, -px);
976 theta = TMath::ATan2(pt, pz);
977
0f482ae4 978//
c2c598a3 979// Save 4-momentum sum for balancing
980 Int_t index = isys;
6e90ad26 981
982 p0[index][0] += px;
983 p0[index][1] += py;
984 p0[index][2] += pz;
985 p0[index][3] += e;
6e90ad26 986
987 klast[index] = i;
988
0f482ae4 989//
990// Fractional energy loss
b280c4cc 991 Double_t z = zquench[index];
4e383037 992
c2c598a3 993
4e383037 994// Don't fully quench radiated gluons
995//
996 if (imo > 1000) {
997// This small factor makes sure that the gluons are not too close in phase space to avoid recombination
998//
999
c2c598a3 1000 z = 0.02;
4e383037 1001 }
c2c598a3 1002// printf("z: %d %f\n", imo, z);
1003
4e383037 1004
1005//
6e90ad26 1006
1007 //
1008 //
1009 // Transform into frame in which initial parton is along z-axis
1010 //
1011 TVector3 v(px, py, pz);
1012 v.RotateZ(-phiq[index]); v.RotateY(-thetaq[index]);
1013 Double_t pxs = v.X(); Double_t pys = v.Y(); Double_t pl = v.Z();
1014
1015 Double_t jt = TMath::Sqrt(pxs * pxs + pys * pys);
1016 Double_t mt2 = jt * jt + m * m;
1017 Double_t zmax = 1.;
1018 //
1019 // Kinematic limit on z
1020 //
4e383037 1021 if (m > 0.) zmax = 1. - m / TMath::Sqrt(m * m + jt * jt);
6e90ad26 1022 //
1023 // Change light-cone kinematics rel. to initial parton
1024 //
1025 Double_t eppzOld = e + pl;
1026 Double_t empzOld = e - pl;
1027
1028 Double_t eppzNew = (1. - z) * eppzOld;
1029 Double_t empzNew = empzOld - mt2 * z / eppzOld;
1030 Double_t eNew = 0.5 * (eppzNew + empzNew);
1031 Double_t plNew = 0.5 * (eppzNew - empzNew);
1032
1033 Double_t jtNew;
1034 //
1035 // if mt very small (or sometimes even < 0 for numerical reasons) set it to 0
1036 Double_t mt2New = eppzNew * empzNew;
1037 if (mt2New < 1.e-8) mt2New = 0.;
4e383037 1038 if (z < zmax) {
1039 if (m * m > mt2New) {
1040 //
1041 // This should not happen
1042 //
1043 Fatal("Quench()", "This should never happen %e %e %e!", m, eppzNew, empzNew);
1044 jtNew = 0;
1045 } else {
1046 jtNew = TMath::Sqrt(mt2New - m * m);
1047 }
6e90ad26 1048 } else {
4e383037 1049 // If pT is to small (probably a leading massive particle) we scale only the energy
1050 // This can cause negative masses of the radiated gluon
1051 // Let's hope for the best ...
1052 jtNew = jt;
1053 eNew = TMath::Sqrt(plNew * plNew + mt2);
1054
6e90ad26 1055 }
6e90ad26 1056 //
1057 // Calculate new px, py
1058 //
1059 Double_t pxNew = jtNew / jt * pxs;
1060 Double_t pyNew = jtNew / jt * pys;
1061
1062// Double_t dpx = pxs - pxNew;
1063// Double_t dpy = pys - pyNew;
1064// Double_t dpz = pl - plNew;
1065// Double_t de = e - eNew;
1066// Double_t dmass2 = de * de - dpx * dpx - dpy * dpy - dpz * dpz;
1067// printf("New mass (1) %e %e %e %e %e %e %e \n", dmass2, jt, jtNew, pl, plNew, e, eNew);
1068// printf("New mass (2) %e %e \n", pxNew, pyNew);
1069 //
1070 // Rotate back
1071 //
1072 TVector3 w(pxNew, pyNew, plNew);
1073 w.RotateY(thetaq[index]); w.RotateZ(phiq[index]);
1074 pxNew = w.X(); pyNew = w.Y(); plNew = w.Z();
1075
1076 p1[index][0] += pxNew;
1077 p1[index][1] += pyNew;
1078 p1[index][2] += plNew;
1079 p1[index][3] += eNew;
1080 //
1081 // Updated 4-momentum vectors
1082 //
1083 pNew[icount][0] = pxNew;
1084 pNew[icount][1] = pyNew;
1085 pNew[icount][2] = plNew;
1086 pNew[icount][3] = eNew;
1087 kNew[icount] = i;
1088 icount++;
1089 } // parton loop
0f482ae4 1090 //
6e90ad26 1091 // Check if there was phase-space for quenching
0f482ae4 1092 //
0f482ae4 1093
6e90ad26 1094 if (icount == 0) quenched[isys] = kFALSE;
1095 if (!quenched[isys]) break;
1096
1097 for (Int_t j = 0; j < 4; j++)
1098 {
1099 p2[isys][j] = p0[isys][j] - p1[isys][j];
1100 }
1101 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 1102 if (p2[isys][4] > 0.) {
1103 p2[isys][4] = TMath::Sqrt(p2[isys][4]);
1104 break;
1105 } else {
b280c4cc 1106 printf("Warning negative mass squared in system %d %f ! \n", isys, zquench[isys]);
4e383037 1107 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 1108 if (p2[isys][4] < -0.01) {
4e383037 1109 printf("Negative mass squared !\n");
1110 // Here we have to put the gluon back to mass shell
1111 // This will lead to a small energy imbalance
1112 p2[isys][4] = 0.;
1113 p2[isys][3] = TMath::Sqrt(p2[isys][0] * p2[isys][0] + p2[isys][1] * p2[isys][1] + p2[isys][2] * p2[isys][2]);
1114 break;
6e90ad26 1115 } else {
1116 p2[isys][4] = 0.;
1117 break;
1118 }
1119 }
6e90ad26 1120 /*
6e90ad26 1121 zHeavy *= 0.98;
1122 printf("zHeavy lowered to %f\n", zHeavy);
1123 if (zHeavy < 0.01) {
1124 printf("No success ! \n");
1125 icount = 0;
1126 quenched[isys] = kFALSE;
1127 break;
1128 }
4e383037 1129 */
1130 } // iteration on z (while)
1131
6e90ad26 1132// Update event record
1133 for (Int_t k = 0; k < icount; k++) {
1134// 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] );
1135 fPyjets->P[0][kNew[k]] = pNew[k][0];
1136 fPyjets->P[1][kNew[k]] = pNew[k][1];
1137 fPyjets->P[2][kNew[k]] = pNew[k][2];
1138 fPyjets->P[3][kNew[k]] = pNew[k][3];
0f482ae4 1139 }
4e383037 1140 //
1141 // Add the gluons
1142 //
1143 Int_t ish = 0;
1837e95c 1144 Int_t iGlu;
4e383037 1145 if (!quenched[isys]) continue;
0f482ae4 1146//
1147// Last parton from shower i
4e383037 1148 Int_t in = klast[isys];
0f482ae4 1149//
1150// Continue if no parton in shower i selected
1151 if (in == -1) continue;
1152//
1153// If this is the second initial parton and it is behind the first move pointer by previous ish
4e383037 1154 if (isys == 1 && klast[1] > klast[0]) in += ish;
0f482ae4 1155//
1156// Starting index
452af8c7 1157
4e383037 1158// jmin = in - 1;
0f482ae4 1159// How many additional gluons will be generated
1160 ish = 1;
4e383037 1161 if (p2[isys][4] > 0.05) ish = 2;
0f482ae4 1162//
1163// Position of gluons
4e383037 1164 iGlu = numpart;
1165 if (iglu == 0) igMin = iGlu;
1166 igMax = iGlu;
0f482ae4 1167 numpart += ish;
1168 (fPyjets->N) += ish;
4e383037 1169
0f482ae4 1170 if (ish == 1) {
4e383037 1171 fPyjets->P[0][iGlu] = p2[isys][0];
1172 fPyjets->P[1][iGlu] = p2[isys][1];
1173 fPyjets->P[2][iGlu] = p2[isys][2];
1174 fPyjets->P[3][iGlu] = p2[isys][3];
1175 fPyjets->P[4][iGlu] = p2[isys][4];
0f482ae4 1176
4e383037 1177 fPyjets->K[0][iGlu] = 1;
1178 if (iglu == nGluon[isys] - 1) fPyjets->K[0][iGlu] = 1;
0f482ae4 1179 fPyjets->K[1][iGlu] = 21;
4e383037 1180 fPyjets->K[2][iGlu] = fPyjets->K[2][in] + 1000;
0f482ae4 1181 fPyjets->K[3][iGlu] = -1;
1182 fPyjets->K[4][iGlu] = -1;
4e383037 1183
1184 pg[0] += p2[isys][0];
1185 pg[1] += p2[isys][1];
1186 pg[2] += p2[isys][2];
1187 pg[3] += p2[isys][3];
0f482ae4 1188 } else {
1189 //
1190 // Split gluon in rest frame.
1191 //
4e383037 1192 Double_t bx = p2[isys][0] / p2[isys][3];
1193 Double_t by = p2[isys][1] / p2[isys][3];
1194 Double_t bz = p2[isys][2] / p2[isys][3];
1195 Double_t pst = p2[isys][4] / 2.;
0f482ae4 1196 //
1197 // Isotropic decay ????
1198 Double_t cost = 2. * gRandom->Rndm() - 1.;
1199 Double_t sint = TMath::Sqrt(1. - cost * cost);
1200 Double_t phi = 2. * TMath::Pi() * gRandom->Rndm();
1201
1202 Double_t pz1 = pst * cost;
1203 Double_t pz2 = -pst * cost;
1204 Double_t pt1 = pst * sint;
1205 Double_t pt2 = -pst * sint;
1206 Double_t px1 = pt1 * TMath::Cos(phi);
1207 Double_t py1 = pt1 * TMath::Sin(phi);
1208 Double_t px2 = pt2 * TMath::Cos(phi);
1209 Double_t py2 = pt2 * TMath::Sin(phi);
1210
1211 fPyjets->P[0][iGlu] = px1;
1212 fPyjets->P[1][iGlu] = py1;
1213 fPyjets->P[2][iGlu] = pz1;
1214 fPyjets->P[3][iGlu] = pst;
1215 fPyjets->P[4][iGlu] = 0.;
1216
4e383037 1217 fPyjets->K[0][iGlu] = 1 ;
0f482ae4 1218 fPyjets->K[1][iGlu] = 21;
4e383037 1219 fPyjets->K[2][iGlu] = fPyjets->K[2][in] + 1000;
0f482ae4 1220 fPyjets->K[3][iGlu] = -1;
1221 fPyjets->K[4][iGlu] = -1;
1222
1223 fPyjets->P[0][iGlu+1] = px2;
1224 fPyjets->P[1][iGlu+1] = py2;
1225 fPyjets->P[2][iGlu+1] = pz2;
1226 fPyjets->P[3][iGlu+1] = pst;
1227 fPyjets->P[4][iGlu+1] = 0.;
1228
4e383037 1229 fPyjets->K[0][iGlu+1] = 1;
1230 if (iglu == nGluon[isys] - 1) fPyjets->K[0][iGlu+1] = 1;
0f482ae4 1231 fPyjets->K[1][iGlu+1] = 21;
4e383037 1232 fPyjets->K[2][iGlu+1] = fPyjets->K[2][in] + 1000;
0f482ae4 1233 fPyjets->K[3][iGlu+1] = -1;
1234 fPyjets->K[4][iGlu+1] = -1;
1235 SetMSTU(1,0);
1236 SetMSTU(2,0);
1237 //
1238 // Boost back
1239 //
1240 Pyrobo(iGlu + 1, iGlu + 2, 0., 0., bx, by, bz);
1241 }
4e383037 1242/*
1243 for (Int_t ig = iGlu; ig < iGlu+ish; ig++) {
1244 Double_t px, py, pz;
1245 px = fPyjets->P[0][ig];
1246 py = fPyjets->P[1][ig];
1247 pz = fPyjets->P[2][ig];
1248 TVector3 v(px, py, pz);
1249 v.RotateZ(-phiq[isys]);
1250 v.RotateY(-thetaq[isys]);
1251 Double_t pxs = v.X(); Double_t pys = v.Y(); Double_t pzs = v.Z();
1252 Double_t r = AliPythiaRndm::GetPythiaRandom()->Rndm();
1253 Double_t jtKick = 0.3 * TMath::Sqrt(-TMath::Log(r));
1254 if (ish == 2) jtKick = wjtKick[i] * TMath::Sqrt(-TMath::Log(r)) / TMath::Sqrt(2.);
1255 Double_t phiKick = 2. * TMath::Pi() * AliPythiaRndm::GetPythiaRandom()->Rndm();
1256 pxs += jtKick * TMath::Cos(phiKick);
1257 pys += jtKick * TMath::Sin(phiKick);
1258 TVector3 w(pxs, pys, pzs);
1259 w.RotateY(thetaq[isys]);
1260 w.RotateZ(phiq[isys]);
1261 fPyjets->P[0][ig] = w.X();
1262 fPyjets->P[1][ig] = w.Y();
1263 fPyjets->P[2][ig] = w.Z();
1264 fPyjets->P[2][ig] = w.Mag();
1265 }
1266*/
1267 } // kGluon
1268
6e90ad26 1269
4e383037 1270 // Check energy conservation
0f482ae4 1271 Double_t pxs = 0.;
1272 Double_t pys = 0.;
1273 Double_t pzs = 0.;
1274 Double_t es = 14000.;
1275
1276 for (Int_t i = 0; i < numpart; i++)
1277 {
1278 kst = fPyjets->K[0][i];
1279 if (kst != 1 && kst != 2) continue;
1280 pxs += fPyjets->P[0][i];
1281 pys += fPyjets->P[1][i];
1282 pzs += fPyjets->P[2][i];
1283 es -= fPyjets->P[3][i];
1284 }
1285 if (TMath::Abs(pxs) > 1.e-2 ||
1286 TMath::Abs(pys) > 1.e-2 ||
1287 TMath::Abs(pzs) > 1.e-1) {
1288 printf("%e %e %e %e\n", pxs, pys, pzs, es);
4e383037 1289// Fatal("Quench()", "4-Momentum non-conservation");
452af8c7 1290 }
4e383037 1291
1292 } // end quenching loop (systems)
6e90ad26 1293// Clean-up
0f482ae4 1294 for (Int_t i = 0; i < numpart; i++)
1295 {
4e383037 1296 imo = fPyjets->K[2][i];
1297 if (imo > 1000) {
1298 fPyjets->K[2][i] = fPyjets->K[2][i] % 1000;
1299 }
0f482ae4 1300 }
4e383037 1301// this->Pylist(1);
0f482ae4 1302} // end quench
90d7b703 1303
992f2843 1304
1305void AliPythia::Pyquen(Double_t a, Int_t ibf, Double_t b)
1306{
1307 // Igor Lokthine's quenching routine
1308 pyquen(a, ibf, b);
1309}
b280c4cc 1310
1311void AliPythia::GetQuenchingParameters(Double_t& xp, Double_t& yp, Double_t z[4])
1312{
1313 // Return event specific quenching parameters
1314 xp = fXJet;
1315 yp = fYJet;
1316 for (Int_t i = 0; i < 4; i++) z[i] = fZQuench[i];
1317
1318}
1319
65f2626c 1320