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
147 SetMSTP(51, 7); // CTEQ5L pdf
148 SetMSTP(81,1); // Multiple Interactions ON
149 SetMSTP(82,4); // Double Gaussian Model
151 SetPARP(82,1.8); // [GeV] PT_min at Ref. energy
152 SetPARP(89,1000.); // [GeV] Ref. energy
153 SetPARP(90,0.16); // 2*epsilon (exponent in power law)
154 SetPARP(83,0.5); // Core density in proton matter distribution (def.value)
155 SetPARP(84,0.5); // Core radius
156 SetPARP(85,0.33); // Regulates gluon prod. mechanism
157 SetPARP(86,0.66); // Regulates gluon prod. mechanism
158 SetPARP(67,1); // Regulates Initial State Radiation
161 // Minimum Bias pp-Collisions
164 // select Pythia min. bias model
166 SetMSUB(95,1); // low pt production
168 SetMSTP(51, 7); // CTEQ5L pdf
169 SetMSTP(81,1); // Multiple Interactions ON
170 SetMSTP(82,4); // Double Gaussian Model
172 SetPARP(82,1.8); // [GeV] PT_min at Ref. energy
173 SetPARP(89,1000.); // [GeV] Ref. energy
174 SetPARP(90,0.16); // 2*epsilon (exponent in power law)
175 SetPARP(83,0.5); // Core density in proton matter distribution (def.value)
176 SetPARP(84,0.5); // Core radius
177 SetPARP(85,0.33); // Regulates gluon prod. mechanism
178 SetPARP(86,0.66); // Regulates gluon prod. mechanism
179 SetPARP(67,1); // Regulates Initial State Radiation
190 case kPyCharmPbPbMNR:
192 // Tuning of Pythia parameters aimed to get a resonable agreement
193 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
194 // c-cbar single inclusive and double differential distributions.
195 // This parameter settings are meant to work with Pb-Pb collisions
196 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
197 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
198 // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
203 // No multiple interactions
208 // Initial/final parton shower on (Pythia default)
230 // Tuning of Pythia parameters aimed to get a resonable agreement
231 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
232 // c-cbar single inclusive and double differential distributions.
233 // This parameter settings are meant to work with p-Pb collisions
234 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
235 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
236 // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
241 // No multiple interactions
246 // Initial/final parton shower on (Pythia default)
268 // Tuning of Pythia parameters aimed to get a resonable agreement
269 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
270 // c-cbar single inclusive and double differential distributions.
271 // This parameter settings are meant to work with pp collisions
272 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
273 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
274 // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
279 // No multiple interactions
284 // Initial/final parton shower on (Pythia default)
304 case kPyBeautyPbPbMNR:
305 // Tuning of Pythia parameters aimed to get a resonable agreement
306 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
307 // b-bbar single inclusive and double differential distributions.
308 // This parameter settings are meant to work with Pb-Pb collisions
309 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
310 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
311 // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
316 // No multiple interactions
321 // Initial/final parton shower on (Pythia default)
343 case kPyBeautypPbMNR:
344 // Tuning of Pythia parameters aimed to get a resonable agreement
345 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
346 // b-bbar single inclusive and double differential distributions.
347 // This parameter settings are meant to work with p-Pb collisions
348 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
349 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
350 // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
355 // No multiple interactions
360 // Initial/final parton shower on (Pythia default)
383 // Tuning of Pythia parameters aimed to get a resonable agreement
384 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
385 // b-bbar single inclusive and double differential distributions.
386 // This parameter settings are meant to work with pp collisions
387 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
388 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
389 // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
394 // No multiple interactions
399 // Initial/final parton shower on (Pythia default)
424 SetMSTP(41,1); // all resonance decays switched on
426 Initialize("CMS","p","p",fEcms);
430 Int_t AliPythia::CheckedLuComp(Int_t kf)
432 // Check Lund particle code (for debugging)
434 printf("\n Lucomp kf,kc %d %d",kf,kc);
438 void AliPythia::SetNuclei(Int_t a1, Int_t a2)
440 // Treat protons as inside nuclei with mass numbers a1 and a2
441 // The MSTP array in the PYPARS common block is used to enable and
442 // select the nuclear structure functions.
443 // MSTP(52) : (D=1) choice of proton and nuclear structure-function library
444 // =1: internal PYTHIA acording to MSTP(51)
445 // =2: PDFLIB proton s.f., with MSTP(51) = 1000xNGROUP+NSET
446 // If the following mass number both not equal zero, nuclear corrections of the stf are used.
447 // MSTP(192) : Mass number of nucleus side 1
448 // MSTP(193) : Mass number of nucleus side 2
455 AliPythia* AliPythia::Instance()
457 // Set random number generator
461 fgAliPythia = new AliPythia();
466 void AliPythia::PrintParticles()
468 // Print list of particl properties
470 char* name = new char[16];
471 for (Int_t kf=0; kf<1000000; kf++) {
472 for (Int_t c = 1; c > -2; c-=2) {
473 Int_t kc = Pycomp(c*kf);
475 Float_t mass = GetPMAS(kc,1);
476 Float_t width = GetPMAS(kc,2);
477 Float_t tau = GetPMAS(kc,4);
483 printf("\n mass, width, tau: %6d %s %10.3f %10.3e %10.3e",
484 c*kf, name, mass, width, tau);
488 printf("\n Number of particles %d \n \n", np);
491 void AliPythia::ResetDecayTable()
493 // Set default values for pythia decay switches
495 for (i = 1; i < 501; i++) SetMDCY(i,1,fDefMDCY[i]);
496 for (i = 1; i < 2001; i++) SetMDME(i,1,fDefMDME[i]);
499 void AliPythia::SetDecayTable()
501 // Set default values for pythia decay switches
504 for (i = 1; i < 501; i++) fDefMDCY[i] = GetMDCY(i,1);
505 for (i = 1; i < 2001; i++) fDefMDME[i] = GetMDME(i,1);
508 void AliPythia::Pyclus(Int_t& njet)
510 // Call Pythia clustering algorithm
515 void AliPythia::Pycell(Int_t& njet)
517 // Call Pythia jet reconstruction algorithm
522 void AliPythia::Pyshow(Int_t ip1, Int_t ip2, Double_t qmax)
524 // Call Pythia jet reconstruction algorithm
526 Int_t numpart = fPyjets->N;
527 for (Int_t i = 0; i < numpart; i++)
529 if (fPyjets->K[2][i] == 7) ip1 = i+1;
530 if (fPyjets->K[2][i] == 8) ip2 = i+1;
534 qmax = 2. * GetVINT(51);
535 printf("Pyshow %d %d %f", ip1, ip2, qmax);
537 pyshow(ip1, ip2, qmax);
540 void AliPythia::Pyrobo(Int_t imi, Int_t ima, Double_t the, Double_t phi, Double_t bex, Double_t bey, Double_t bez)
542 pyrobo(imi, ima, the, phi, bex, bey, bez);
547 void AliPythia::Quench()
551 // Simple Jet Quenching routine:
552 // =============================
553 // The jet formed by all final state partons radiated by the parton created
554 // in the hard collisions is quenched by a factor z using:
555 // (E + p_z)new = (1-z) (E + p_z)old
557 // The lost momentum is first balanced by one gluon with virtuality > 0.
558 // Subsequently the gluon splits to yield two gluons with E = p.
563 Int_t klast[2] = {-1, -1};
565 for (Int_t i = 0; i < 4; i++)
575 Int_t numpart = fPyjets->N;
577 for (Int_t i = 0; i < numpart; i++)
579 Int_t imo = fPyjets->K[2][i];
580 Int_t kst = fPyjets->K[0][i];
581 Int_t pdg = fPyjets->K[1][i];
583 // Quarks and gluons only
584 if (pdg != 21 && TMath::Abs(pdg) > 6) continue;
586 // Particles from hard scattering only
589 Float_t px = fPyjets->P[0][i];
590 Float_t py = fPyjets->P[1][i];
591 Float_t pz = fPyjets->P[2][i];
592 Float_t e = fPyjets->P[3][i];
593 Float_t m = fPyjets->P[4][i];
594 Float_t pt = TMath::Sqrt(px * px + py * py);
595 // Skip comment lines
596 if (kst != 1 && kst != 2) continue;
598 Float_t mt = TMath::Sqrt(px * px + py * py + m * m);
601 // Some cuts to be in a save kinematic region
603 if (imo != 7 && imo != 8) continue;
604 Int_t index = imo - 7;
616 Float_t eppzOld = e + pz;
617 Float_t empzOld = e - pz;
621 // Kinematics of the original parton
624 Float_t eppzNew = (1. - z) * eppzOld;
625 Float_t empzNew = empzOld - mt * mt * z / eppzOld;
626 Float_t eNew0 = 0.5 * (eppzNew + empzNew);
627 Float_t pzNew0 = 0.5 * (eppzNew - empzNew);
629 // Skip if pt too small
631 if (m * m > eppzNew * empzNew) continue;
632 Float_t ptNew = TMath::Sqrt(eppzNew * empzNew - m * m);
633 Float_t pxNew0 = ptNew / pt * px;
634 Float_t pyNew0 = ptNew / pt * py;
636 p1[index][0] += pxNew0;
637 p1[index][1] += pyNew0;
638 p1[index][2] += pzNew0;
639 p1[index][3] += eNew0;
641 // Update event record
643 fPyjets->P[0][i] = pxNew0;
644 fPyjets->P[1][i] = pyNew0;
645 fPyjets->P[2][i] = pzNew0;
646 fPyjets->P[3][i] = eNew0;
654 for (Int_t k = 0; k < 2; k++)
656 for (Int_t j = 0; j < 4; j++)
658 p2[k][j] = p0[k][j] - p1[k][j];
660 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];
666 // Bring gluon back to mass shell via momentum scaling
667 // (momentum will not be conserved, but energy)
671 Float_t psq = p2[k][0] * p2[k][0] + p2[k][1] * p2[k][1] + p2[k][2] * p2[k][2];
672 Float_t fact = TMath::Sqrt(1. + p2[k][4] / psq);
676 p2[k][3] = TMath::Sqrt(psq) * fact;
683 p2[0][4] = TMath::Sqrt(p2[0][4]);
685 printf("Warning negative mass squared ! \n");
689 p2[1][4] = TMath::Sqrt(p2[1][4]);
691 printf("Warning negative mass squared ! \n");
699 for (Int_t i = 0; i < 2; i++) {
700 Int_t ish, jmin, jmax, iGlu, iNew;
704 if (in == -1) continue;
705 if (i == 1 && klast[1] > klast[0]) in += ish;
710 if (p2[i][4] > 0) ish = 2;
714 jmax = numpart + ish - 1;
716 if (fPyjets->K[0][in-1] == 1 || fPyjets->K[0][in-1] == 21 || fPyjets->K[0][in-1] == 11) {
724 for (Int_t j = jmax; j > jmin; j--)
726 for (Int_t k = 0; k < 5; k++) {
727 fPyjets->K[k][j] = fPyjets->K[k][j-ish];
728 fPyjets->P[k][j] = fPyjets->P[k][j-ish];
729 fPyjets->V[k][j] = fPyjets->V[k][j-ish];
736 fPyjets->P[0][iGlu] = p2[i][0];
737 fPyjets->P[1][iGlu] = p2[i][1];
738 fPyjets->P[2][iGlu] = p2[i][2];
739 fPyjets->P[3][iGlu] = p2[i][3];
740 fPyjets->P[4][iGlu] = p2[i][4];
742 fPyjets->K[0][iGlu] = 2;
743 fPyjets->K[1][iGlu] = 21;
744 fPyjets->K[2][iGlu] = fPyjets->K[2][iNew];
745 fPyjets->K[3][iGlu] = -1;
746 fPyjets->K[4][iGlu] = -1;
749 // Split gluon in rest frame.
751 Double_t bx = p2[i][0] / p2[i][3];
752 Double_t by = p2[i][1] / p2[i][3];
753 Double_t bz = p2[i][2] / p2[i][3];
755 Float_t pst = p2[i][4] / 2.;
757 Float_t cost = 2. * gRandom->Rndm() - 1.;
758 Float_t sint = TMath::Sqrt(1. - cost * cost);
759 Float_t phi = 2. * TMath::Pi() * gRandom->Rndm();
761 Float_t pz1 = pst * cost;
762 Float_t pz2 = -pst * cost;
763 Float_t pt1 = pst * sint;
764 Float_t pt2 = -pst * sint;
765 Float_t px1 = pt1 * TMath::Cos(phi);
766 Float_t py1 = pt1 * TMath::Sin(phi);
767 Float_t px2 = pt2 * TMath::Cos(phi);
768 Float_t py2 = pt2 * TMath::Sin(phi);
770 fPyjets->P[0][iGlu] = px1;
771 fPyjets->P[1][iGlu] = py1;
772 fPyjets->P[2][iGlu] = pz1;
773 fPyjets->P[3][iGlu] = pst;
774 fPyjets->P[4][iGlu] = 0.;
776 fPyjets->K[0][iGlu] = 2;
777 fPyjets->K[1][iGlu] = 21;
778 fPyjets->K[2][iGlu] = fPyjets->K[2][iNew];
779 fPyjets->K[3][iGlu] = -1;
780 fPyjets->K[4][iGlu] = -1;
782 fPyjets->P[0][iGlu+1] = px2;
783 fPyjets->P[1][iGlu+1] = py2;
784 fPyjets->P[2][iGlu+1] = pz2;
785 fPyjets->P[3][iGlu+1] = pst;
786 fPyjets->P[4][iGlu+1] = 0.;
788 fPyjets->K[0][iGlu+1] = 2;
789 fPyjets->K[1][iGlu+1] = 21;
790 fPyjets->K[2][iGlu+1] = fPyjets->K[2][iNew];
791 fPyjets->K[3][iGlu+1] = -1;
792 fPyjets->K[4][iGlu+1] = -1;
799 Pyrobo(iGlu + 1, iGlu + 2, 0., 0., bx, by, bz);
802 } // end adding gluons