1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 #include "AliPythia.h"
19 #include "AliPythiaRndm.h"
24 # define pyclus pyclus_
25 # define pycell pycell_
26 # define pyshow pyshow_
27 # define pyrobo pyrobo_
30 # define pyclus PYCLUS
31 # define pycell PYCELL
32 # define pyrobo PYROBO
33 # define type_of_call _stdcall
36 extern "C" void type_of_call pyclus(Int_t & );
37 extern "C" void type_of_call pycell(Int_t & );
38 extern "C" void type_of_call pyshow(Int_t &, Int_t &, Double_t &);
39 extern "C" void type_of_call pyrobo(Int_t &, Int_t &, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &);
41 //_____________________________________________________________________________
43 AliPythia* AliPythia::fgAliPythia=NULL;
45 AliPythia::AliPythia()
47 // Default Constructor
50 if (!AliPythiaRndm::GetPythiaRandom())
51 AliPythiaRndm::SetPythiaRandom(GetRandom());
55 void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfunc)
57 // Initialise the process to generate
58 if (!AliPythiaRndm::GetPythiaRandom())
59 AliPythiaRndm::SetPythiaRandom(GetRandom());
63 fStrucFunc = strucfunc;
65 SetMDCY(Pycomp(111),1,0);
66 // select structure function
68 SetMSTP(51,strucfunc);
70 // Pythia initialisation for selected processes//
74 for (Int_t i=1; i<= 200; i++) {
77 // select charm production
115 case kPyCharmUnforced:
124 case kPyBeautyUnforced:
134 // Minimum Bias pp-Collisions
137 // select Pythia min. bias model
139 SetMSUB(92,1); // single diffraction AB-->XB
140 SetMSUB(93,1); // single diffraction AB-->AX
141 SetMSUB(94,1); // double diffraction
142 SetMSUB(95,1); // low pt production
143 SetMSTP(81,1); // multiple interactions switched on
144 SetMSTP(82,3); // model with varying impact param. & a single Gaussian
145 SetPARP(82,3.47); // set value pT_0 for turn-off of the cross section of
146 // multiple interaction at a reference energy = 14000 GeV
147 SetPARP(89,14000.); // reference energy for the above parameter
148 SetPARP(90,0.174); // set exponent for energy dependence of pT_0
150 // Minimum Bias pp-Collisions
153 // select Pythia min. bias model
155 SetMSUB(95,1); // low pt production
156 SetMSTP(81,1); // multiple interactions switched on
157 SetMSTP(82,3); // model with varying impact param. & a single Gaussian
158 SetPARP(82,3.47); // set value pT_0 for turn-off of the cross section of
159 // multiple interaction at a reference energy = 14000 GeV
160 SetPARP(89,14000.); // reference energy for the above parameter
161 SetPARP(90,0.174); // set exponent for energy dependence of pT_0
173 case kPyCharmPbPbMNR:
175 // Tuning of Pythia parameters aimed to get a resonable agreement
176 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
177 // c-cbar single inclusive and double differential distributions.
178 // This parameter settings are meant to work with Pb-Pb collisions
179 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
180 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
181 // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
186 // No multiple interactions
191 // Initial/final parton shower on (Pythia default)
213 // Tuning of Pythia parameters aimed to get a resonable agreement
214 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
215 // c-cbar single inclusive and double differential distributions.
216 // This parameter settings are meant to work with p-Pb collisions
217 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
218 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
219 // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
224 // No multiple interactions
229 // Initial/final parton shower on (Pythia default)
251 // Tuning of Pythia parameters aimed to get a resonable agreement
252 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
253 // c-cbar single inclusive and double differential distributions.
254 // This parameter settings are meant to work with pp collisions
255 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
256 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
257 // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
262 // No multiple interactions
267 // Initial/final parton shower on (Pythia default)
287 case kPyBeautyPbPbMNR:
288 // Tuning of Pythia parameters aimed to get a resonable agreement
289 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
290 // b-bbar single inclusive and double differential distributions.
291 // This parameter settings are meant to work with Pb-Pb collisions
292 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
293 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
294 // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
299 // No multiple interactions
304 // Initial/final parton shower on (Pythia default)
326 case kPyBeautypPbMNR:
327 // Tuning of Pythia parameters aimed to get a resonable agreement
328 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
329 // b-bbar single inclusive and double differential distributions.
330 // This parameter settings are meant to work with p-Pb collisions
331 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
332 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
333 // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
338 // No multiple interactions
343 // Initial/final parton shower on (Pythia default)
366 // Tuning of Pythia parameters aimed to get a resonable agreement
367 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
368 // b-bbar single inclusive and double differential distributions.
369 // This parameter settings are meant to work with pp collisions
370 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
371 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
372 // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
377 // No multiple interactions
382 // Initial/final parton shower on (Pythia default)
407 SetMSTP(41,1); // all resonance decays switched on
409 Initialize("CMS","p","p",fEcms);
413 Int_t AliPythia::CheckedLuComp(Int_t kf)
415 // Check Lund particle code (for debugging)
417 printf("\n Lucomp kf,kc %d %d",kf,kc);
421 void AliPythia::SetNuclei(Int_t a1, Int_t a2)
423 // Treat protons as inside nuclei with mass numbers a1 and a2
424 // The MSTP array in the PYPARS common block is used to enable and
425 // select the nuclear structure functions.
426 // MSTP(52) : (D=1) choice of proton and nuclear structure-function library
427 // =1: internal PYTHIA acording to MSTP(51)
428 // =2: PDFLIB proton s.f., with MSTP(51) = 1000xNGROUP+NSET
429 // If the following mass number both not equal zero, nuclear corrections of the stf are used.
430 // MSTP(192) : Mass number of nucleus side 1
431 // MSTP(193) : Mass number of nucleus side 2
438 AliPythia* AliPythia::Instance()
440 // Set random number generator
444 fgAliPythia = new AliPythia();
449 void AliPythia::PrintParticles()
451 // Print list of particl properties
453 char* name = new char[16];
454 for (Int_t kf=0; kf<1000000; kf++) {
455 for (Int_t c = 1; c > -2; c-=2) {
456 Int_t kc = Pycomp(c*kf);
458 Float_t mass = GetPMAS(kc,1);
459 Float_t width = GetPMAS(kc,2);
460 Float_t tau = GetPMAS(kc,4);
466 printf("\n mass, width, tau: %6d %s %10.3f %10.3e %10.3e",
467 c*kf, name, mass, width, tau);
471 printf("\n Number of particles %d \n \n", np);
474 void AliPythia::ResetDecayTable()
476 // Set default values for pythia decay switches
478 for (i = 1; i < 501; i++) SetMDCY(i,1,fDefMDCY[i]);
479 for (i = 1; i < 2001; i++) SetMDME(i,1,fDefMDME[i]);
482 void AliPythia::SetDecayTable()
484 // Set default values for pythia decay switches
487 for (i = 1; i < 501; i++) fDefMDCY[i] = GetMDCY(i,1);
488 for (i = 1; i < 2001; i++) fDefMDME[i] = GetMDME(i,1);
491 void AliPythia::Pyclus(Int_t& njet)
493 // Call Pythia clustering algorithm
498 void AliPythia::Pycell(Int_t& njet)
500 // Call Pythia jet reconstruction algorithm
505 void AliPythia::Pyshow(Int_t ip1, Int_t ip2, Double_t qmax)
507 // Call Pythia jet reconstruction algorithm
509 Int_t numpart = fPyjets->N;
510 for (Int_t i = 0; i < numpart; i++)
512 if (fPyjets->K[2][i] == 7) ip1 = i+1;
513 if (fPyjets->K[2][i] == 8) ip2 = i+1;
517 qmax = 2. * GetVINT(51);
518 printf("Pyshow %d %d %f", ip1, ip2, qmax);
520 pyshow(ip1, ip2, qmax);
523 void AliPythia::Pyrobo(Int_t imi, Int_t ima, Double_t the, Double_t phi, Double_t bex, Double_t bey, Double_t bez)
525 pyrobo(imi, ima, the, phi, bex, bey, bez);
530 void AliPythia::Quench()
534 // Simple Jet Quenching routine:
535 // =============================
536 // The jet formed by all final state partons radiated by the parton created
537 // in the hard collisions is quenched by a factor z using:
538 // (E + p_z)new = (1-z) (E + p_z)old
540 // The lost momentum is first balanced by one gluon with virtuality > 0.
541 // Subsequently the gluon splits to yield two gluons with E = p.
546 Int_t klast[2] = {-1, -1};
548 for (Int_t i = 0; i < 4; i++)
558 Int_t numpart = fPyjets->N;
560 for (Int_t i = 0; i < numpart; i++)
562 Int_t imo = fPyjets->K[2][i];
563 Int_t kst = fPyjets->K[0][i];
564 Int_t pdg = fPyjets->K[1][i];
566 // Quarks and gluons only
567 if (pdg != 21 && TMath::Abs(pdg) > 6) continue;
569 // Particles from hard scattering only
572 Float_t px = fPyjets->P[0][i];
573 Float_t py = fPyjets->P[1][i];
574 Float_t pz = fPyjets->P[2][i];
575 Float_t e = fPyjets->P[3][i];
576 Float_t m = fPyjets->P[4][i];
577 Float_t pt = TMath::Sqrt(px * px + py * py);
578 // Skip comment lines
579 if (kst != 1 && kst != 2) continue;
581 Float_t mt = TMath::Sqrt(px * px + py * py + m * m);
584 // Some cuts to be in a save kinematic region
586 if (imo != 7 && imo != 8) continue;
587 Int_t index = imo - 7;
600 Float_t eppzOld = e + pz;
601 Float_t empzOld = e - pz;
605 // Kinematics of the original parton
608 Float_t eppzNew = (1. - z) * eppzOld;
609 Float_t empzNew = empzOld - mt * mt * z / eppzOld;
610 Float_t eNew0 = 0.5 * (eppzNew + empzNew);
611 Float_t pzNew0 = 0.5 * (eppzNew - empzNew);
613 // Skip if pt too small
615 if (m * m > eppzNew * empzNew) continue;
616 Float_t ptNew = TMath::Sqrt(eppzNew * empzNew - m * m);
617 Float_t pxNew0 = ptNew / pt * px;
618 Float_t pyNew0 = ptNew / pt * py;
620 p1[index][0] += pxNew0;
621 p1[index][1] += pyNew0;
622 p1[index][2] += pzNew0;
623 p1[index][3] += eNew0;
625 // Update event record
627 fPyjets->P[0][i] = pxNew0;
628 fPyjets->P[1][i] = pyNew0;
629 fPyjets->P[2][i] = pzNew0;
630 fPyjets->P[3][i] = eNew0;
638 for (Int_t k = 0; k < 2; k++)
640 for (Int_t j = 0; j < 4; j++)
642 p2[k][j] = p0[k][j] - p1[k][j];
644 p2[k][4] = p2[k][3] * p2[k][3] - p2[k][0] * p2[k][0] - p2[k][1] * p2[k][1] - p2[k][2] * p2[k][2];
650 // Bring gluon back to mass shell via momentum scaling
651 // (momentum will not be conserved, but energy)
655 Float_t psq = p2[k][0] * p2[k][0] + p2[k][1] * p2[k][1] + p2[k][2] * p2[k][2];
656 Float_t fact = TMath::Sqrt(1. + p2[k][4] / psq);
660 p2[k][3] = TMath::Sqrt(psq) * fact;
667 p2[0][4] = TMath::Sqrt(p2[0][4]);
669 printf("Warning negative mass squared ! \n");
673 p2[1][4] = TMath::Sqrt(p2[1][4]);
675 printf("Warning negative mass squared ! \n");
683 for (Int_t i = 0; i < 2; i++) {
684 Int_t ish, jmin, jmax, iGlu, iNew;
688 if (in == -1) continue;
689 if (i == 1 && klast[1] > klast[0]) in += ish;
694 if (p2[i][4] > 0) ish = 2;
698 jmax = numpart + ish - 1;
700 if (fPyjets->K[0][in-1] == 1 || fPyjets->K[0][in-1] == 21 || fPyjets->K[0][in-1] == 11) {
708 for (Int_t j = jmax; j > jmin; j--)
710 for (Int_t k = 0; k < 5; k++) {
711 fPyjets->K[k][j] = fPyjets->K[k][j-ish];
712 fPyjets->P[k][j] = fPyjets->P[k][j-ish];
713 fPyjets->V[k][j] = fPyjets->V[k][j-ish];
720 fPyjets->P[0][iGlu] = p2[i][0];
721 fPyjets->P[1][iGlu] = p2[i][1];
722 fPyjets->P[2][iGlu] = p2[i][2];
723 fPyjets->P[3][iGlu] = p2[i][3];
724 fPyjets->P[4][iGlu] = p2[i][4];
726 fPyjets->K[0][iGlu] = 2;
727 fPyjets->K[1][iGlu] = 21;
728 fPyjets->K[2][iGlu] = fPyjets->K[2][iNew];
729 fPyjets->K[3][iGlu] = -1;
730 fPyjets->K[4][iGlu] = -1;
733 // Split gluon in rest frame.
735 Double_t bx = p2[i][0] / p2[i][3];
736 Double_t by = p2[i][1] / p2[i][3];
737 Double_t bz = p2[i][2] / p2[i][3];
739 Float_t pst = p2[i][4] / 2.;
741 Float_t cost = 2. * gRandom->Rndm() - 1.;
742 Float_t sint = TMath::Sqrt(1. - cost * cost);
743 Float_t phi = 2. * TMath::Pi() * gRandom->Rndm();
745 Float_t pz1 = pst * cost;
746 Float_t pz2 = -pst * cost;
747 Float_t pt1 = pst * sint;
748 Float_t pt2 = -pst * sint;
749 Float_t px1 = pt1 * TMath::Cos(phi);
750 Float_t py1 = pt1 * TMath::Sin(phi);
751 Float_t px2 = pt2 * TMath::Cos(phi);
752 Float_t py2 = pt2 * TMath::Sin(phi);
754 fPyjets->P[0][iGlu] = px1;
755 fPyjets->P[1][iGlu] = py1;
756 fPyjets->P[2][iGlu] = pz1;
757 fPyjets->P[3][iGlu] = pst;
758 fPyjets->P[4][iGlu] = 0.;
760 fPyjets->K[0][iGlu] = 2;
761 fPyjets->K[1][iGlu] = 21;
762 fPyjets->K[2][iGlu] = fPyjets->K[2][iNew];
763 fPyjets->K[3][iGlu] = -1;
764 fPyjets->K[4][iGlu] = -1;
766 fPyjets->P[0][iGlu+1] = px2;
767 fPyjets->P[1][iGlu+1] = py2;
768 fPyjets->P[2][iGlu+1] = pz2;
769 fPyjets->P[3][iGlu+1] = pst;
770 fPyjets->P[4][iGlu+1] = 0.;
772 fPyjets->K[0][iGlu+1] = 2;
773 fPyjets->K[1][iGlu+1] = 21;
774 fPyjets->K[2][iGlu+1] = fPyjets->K[2][iNew];
775 fPyjets->K[3][iGlu+1] = -1;
776 fPyjets->K[4][iGlu+1] = -1;
783 Pyrobo(iGlu + 1, iGlu + 2, 0., 0., bx, by, bz);
786 } // end adding gluons