/*
$Log$
+Revision 1.14 2001/03/09 13:03:40 morsch
+Process_t and Struc_Func_t moved to AliPythia.h
+
+Revision 1.13 2000/12/18 08:55:35 morsch
+Make AliPythia dependent generartors work with new scheme of random number generation
+
+Revision 1.12 2000/11/30 07:12:50 alibrary
+Introducing new Rndm and QA classes
+
+Revision 1.11 2000/10/20 06:30:06 fca
+Use version 0 to avoid streamer generation
+
+Revision 1.10 2000/10/06 14:18:44 morsch
+Upper cut of prim. pT distribution set to 5. GeV
+
+Revision 1.9 2000/09/18 10:41:35 morsch
+Add possibility to use nuclear structure functions from PDF library V8.
+
+Revision 1.8 2000/09/06 14:26:24 morsch
+Decayer functionality of AliPythia has been moved to AliDecayerPythia.
+Class is now a singleton.
+
+Revision 1.7 2000/06/09 20:34:50 morsch
+All coding rule violations except RS3 corrected
+
+Revision 1.6 1999/11/09 07:38:48 fca
+Changes for compatibility with version 2.23 of ROOT
+
Revision 1.5 1999/11/03 17:43:20 fca
New version from G.Martinez & A.Morsch
#include "AliPythia.h"
-#include "AliMC.h"
-ClassImp(AliPythia)
-
-#ifndef WIN32
-# define lu1ent lu1ent_
-# define type_of_call
-#else
-# define lu1ent LU1ENT
-# define type_of_call _stdcall
-#endif
-
-extern "C" void type_of_call
- lu1ent(Int_t&, Int_t&, Float_t&, Float_t&, Float_t&);
-
+ClassImp(AliPythia)
//_____________________________________________________________________________
-Int_t AliPythia::fgInit=0;
+AliPythia* AliPythia::fgAliPythia=NULL;
AliPythia::AliPythia()
{
- for (Int_t i=0; i<501; i++) {
- fGPCode[i][0]=0;
- fGPCode[i][1]=0;
- }
-}
-
-void AliPythia::Lu1Ent(Int_t flag, Int_t idpart,
- Float_t mom, Float_t theta,Float_t phi)
-{
- printf("%d %d %f %f %f\n",flag, idpart, mom, theta, phi);
- lu1ent(flag, idpart, mom, theta, phi);
+// Default Constructor
+//
+// Set random number
+ if (!sRandom) sRandom=fRandom;
}
-void AliPythia::DecayParticle(Int_t idpart,
- Float_t mom, Float_t theta,Float_t phi)
-{
- Lu1Ent(0, idpart, mom, theta, phi);
- GetPrimaries();
-}
void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfunc)
{
+// Initialise the process to generate
fProcess = process;
fEcms = energy;
fStrucFunc = strucfunc;
// don't decay p0
- SetMDCY(Lucomp(111),1,0);
+ SetMDCY(Pycomp(111),1,0);
// select structure function
SetMSTP(52,2);
SetMSTP(51,strucfunc);
// select charm production
switch (process)
{
- case charm:
+ case kPyCharm:
SetMSEL(4);
//
// heavy quark masses
//
// primordial pT
SetMSTP(91,1);
- SetPARP(91,1);
- SetPARP(93,3);
+ SetPARP(91,1.);
+ SetPARP(93,5.);
//
break;
- case beauty:
+ case kPyBeauty:
SetMSEL(5);
SetPMAS(5,1,4.75);
break;
- case jpsi:
+ case kPyJpsi:
SetMSEL(0);
// gg->J/Psi g
SetMSUB(86,1);
break;
- case jpsi_chi:
+ case kPyJpsiChi:
SetMSEL(0);
// gg->J/Psi g
SetMSUB(86,1);
SetMSUB(88,1);
// gg-> chi_2c g
SetMSUB(89,1);
- case charm_unforced:
+ case kPyCharmUnforced:
SetMSEL(0);
// gq->qg
SetMSUB(28,1);
SetMSUB(53,1);
// gg->gg
SetMSUB(68,1);
- case beauty_unforced:
+ case kPyBeautyUnforced:
SetMSEL(0);
// gq->qg
SetMSUB(28,1);
// gg->gg
SetMSUB(68,1);
break;
- case mb:
+ case kPyMb:
// Minimum Bias pp-Collisions
//
// Tuning of parameters descibed in G. Ciapetti and A. Di Ciaccio
// Proc. of the LHC Workshop, Aachen 1990, Vol. II p. 155
//
// select Pythia min. bias model
- SetMSEL(2);
+ SetMSEL(0);
SetMSUB(92,1);
SetMSUB(93,1);
SetMSUB(94,1);
SetPARP(86,0.9);
// 90% of gluon interactions have minimum string length
SetPARP(85,0.9);
+ break;
+ case kPyJets:
+ SetMSEL(1);
+ break;
+ case kPyDirectGamma:
+ SetMSEL(10);
+ break;
}
//
// Initialize PYTHIA
+ SetMSTP(41,1);
+
Initialize("CMS","p","p",fEcms);
-}
-Int_t AliPythia::CountProducts(Int_t channel, Int_t particle)
-{
- Int_t np=0;
- for (Int_t i=1; i<=5; i++) {
- if (TMath::Abs(GetKFDP(channel,i)) == particle) np++;
- }
- return np;
}
-void AliPythia::AllowAllDecays()
+Int_t AliPythia::CheckedLuComp(Int_t kf)
{
- Int_t i;
- for (i=1; i<= 2000; i++) {
- SetMDME(i,1,1);
- }
-//
- for (i=0; i<501; i++){
- fBraPart[i]=1;
- }
+// Check Lund particle code (for debugging)
+ Int_t kc=Pycomp(kf);
+ printf("\n Lucomp kf,kc %d %d",kf,kc);
+ return kc;
}
-void AliPythia::ForceParticleDecay(Int_t particle, Int_t product, Int_t mult)
+void AliPythia::SetNuclei(Int_t a1, Int_t a2)
{
-//
-// force decay of particle into products with multiplicity mult
-
- Int_t kc=Lucomp(particle);
- SetMDCY(kc,1,1);
- Int_t ifirst=GetMDCY(kc,2);
- Int_t ilast=ifirst+GetMDCY(kc,3)-1;
- fBraPart[kc] = 1;
-//
-// Loop over decay channels
- for (Int_t channel=ifirst; channel<=ilast;channel++) {
- if (CountProducts(channel,product) >= mult) {
- SetMDME(channel,1,1);
- } else {
- SetMDME(channel,1,0);
- fBraPart[kc]-=GetBRAT(channel);
- }
- }
+// Treat protons as inside nuclei with mass numbers a1 and a2
+// The MSTP array in the PYPARS common block is used to enable and
+// select the nuclear structure functions.
+// MSTP(52) : (D=1) choice of proton and nuclear structure-function library
+// =1: internal PYTHIA acording to MSTP(51)
+// =2: PDFLIB proton s.f., with MSTP(51) = 1000xNGROUP+NSET
+// =3: PDFLIB proton s.f. with nuclar correction:
+// MSTP( 51) = 1000xNPGROUP+NPSET
+// MSTP(151) = 1000xNAGROUP+NASET
+// MSTP(192) : Mass number of nucleus side 1
+// MSTP(193) : Mass number of nucleus side 2
+
+ SetMSTP(52,3);
+ SetMSTP(191, 1001);
+ SetMSTP(192, a1);
+ SetMSTP(193, a2);
}
-
-void AliPythia::ForceDecay(Decay_t decay)
-{
- fDecay=decay;
-//
-// Make clean
-// AllowAllDecays();
-//
-// select mode
-
- switch (decay)
- {
- case semimuonic:
- if (fProcess==charm || fProcess == charm_unforced) {
- ForceParticleDecay( 411,13,1); // D+/-
- ForceParticleDecay( 421,13,1); // D0
- ForceParticleDecay( 431,13,1); // D_s
- ForceParticleDecay( 4122,13,1); // Lambda_c
- }
- if (fProcess==beauty || fProcess == beauty_unforced) {
- ForceParticleDecay( 511,13,1); // B0
- ForceParticleDecay( 521,13,1); // B+/-
- ForceParticleDecay( 531,13,1); // B_s
- ForceParticleDecay( 5122,13,1); // Lambda_b
- }
- break;
- case dimuon:
- ForceParticleDecay( 41,13,2); // phi
- ForceParticleDecay( 443,13,2); // J/Psi
- ForceParticleDecay(30443,13,2); // Psi'
- ForceParticleDecay( 553,13,2); // Upsilon
- ForceParticleDecay(30553,13,2); // Upsilon'
- break;
- case semielectronic:
-
- ForceParticleDecay( 411,11,1); // D+/-
- ForceParticleDecay( 421,11,1); // D0
- ForceParticleDecay( 431,11,1); // D_s
- ForceParticleDecay( 4122,11,1); // Lambda_c
- ForceParticleDecay( 511,11,1); // B0
- ForceParticleDecay( 521,11,1); // B+/-
- ForceParticleDecay( 531,11,1); // B_s
- ForceParticleDecay( 5122,11,1); // Lambda_b
- break;
- case dielectron:
- ForceParticleDecay( 41,11,2); // phi
- ForceParticleDecay( 443,11,2); // J/Psi
- ForceParticleDecay(30443,11,2); // Psi'
- ForceParticleDecay( 553,11,2); // Upsilon
- ForceParticleDecay(30553,11,2); // Upsilon'
- break;
- case b_jpsi_dimuon:
- ForceParticleDecay( 511,443,1); // B0
- ForceParticleDecay( 521,443,1); // B+/-
- ForceParticleDecay( 531,443,1); // B_s
- ForceParticleDecay( 5122,443,1); // Lambda_b
- ForceParticleDecay( 443,13,2); // J/Psi
- break;
- case b_psip_dimuon:
- ForceParticleDecay( 511,30443,1); // B0
- ForceParticleDecay( 521,30443,1); // B+/-
- ForceParticleDecay( 531,30443,1); // B_s
- ForceParticleDecay( 5122,30443,1); // Lambda_b
- ForceParticleDecay(30443,13,2); // Psi'
- break;
- case b_jpsi_dielectron:
- ForceParticleDecay( 511,443,1); // B0
- ForceParticleDecay( 521,443,1); // B+/-
- ForceParticleDecay( 531,443,1); // B_s
- ForceParticleDecay( 5122,443,1); // Lambda_b
- ForceParticleDecay( 443,11,2); // J/Psi
- break;
- case b_psip_dielectron:
- ForceParticleDecay( 511,30443,1); // B0
- ForceParticleDecay( 521,30443,1); // B+/-
- ForceParticleDecay( 531,30443,1); // B_s
- ForceParticleDecay( 5122,30443,1); // Lambda_b
- ForceParticleDecay(30443,11,2); // Psi'
- break;
- case pitomu:
- ForceParticleDecay(211,13,1); // pi->mu
- break;
- case katomu:
- ForceParticleDecay(321,13,1); // K->mu
- break;
- case all:
- case nodecay:
- break;
+AliPythia* AliPythia::Instance()
+{
+// Set random number generator
+ if (fgAliPythia) {
+ return fgAliPythia;
+ } else {
+ fgAliPythia = new AliPythia();
+ return fgAliPythia;
}
}
-
- void AliPythia::DefineParticles()
-{
- if (fgInit) return;
- fgInit=1;
-
- Float_t mass;
- Float_t tlife;
- Int_t kc, nkc, i;
-//
-//
-// Some particles cloned for rare decays
-//
-// phi-> mu+mu- and phi -> e+e-
-// clone the original phi
- kc = Lucomp(333);
- nkc = 41;
-
- for (i=1;i<=3;i++) {
- SetKCHG(nkc,i,GetKCHG(kc,i));
- }
+void AliPythia::PrintParticles()
+{
+// Print list of particl properties
+ Int_t np = 0;
- for (i=1;i<=4;i++) {
- SetPMAS(nkc,i,GetPMAS(kc,i));
+ for (Int_t kf=0; kf<1000000; kf++) {
+ for (Int_t c = 1; c > -2; c-=2) {
+
+ Int_t kc = Pycomp(c*kf);
+ if (kc) {
+ Float_t mass = GetPMAS(kc,1);
+ Float_t width = GetPMAS(kc,2);
+ Float_t tau = GetPMAS(kc,4);
+
+ char* name = new char[8];
+ Pyname(kf,name);
+
+ np++;
+
+ printf("\n mass, width, tau: %6d %s %10.3f %10.3e %10.3e",
+ c*kf, name, mass, width, tau);
+ }
+ }
}
- SetCHAF(nkc,GetCHAF(kc));
- fBraPart[kc]=1;
-//
-// decay
- SetMDCY(nkc,1,1);
- SetMDCY(nkc,2,993);
- SetMDCY(nkc,3,2);
-//
-// phi-> e+e-
- SetMDME(993,1,1);
- SetMDME(993,2,0);
- SetBRAT(993,2.99e-4);
- SetKFDP(993,1,+11);
- SetKFDP(993,2,-11);
- SetKFDP(993,3,0);
- SetKFDP(993,4,0);
- SetKFDP(993,5,0);
-//
-// phi-> mu+mu-
- SetMDME(994,1,1);
- SetMDME(994,2,0);
- SetBRAT(994,2.5e-4);
- SetKFDP(994,1,+13);
- SetKFDP(994,2,-13);
- SetKFDP(994,3,0);
- SetKFDP(994,4,0);
- SetKFDP(994,5,0);
-//
-// Vector mesons
-//
-// phi clone for dilepton decay-channel
- kc=Lucomp(41);
- mass =GetPMAS(kc,1);
- tlife=GetPMAS(kc,4);
- gMC->Gspart(113,"Phi",3,mass,0,tlife);
- fGPCode[kc][0]=113;
- fGPCode[kc][1]=113;
- // J/Psi
- kc=Lucomp(443);
- mass =GetPMAS(kc,1);
- tlife=GetPMAS(kc,4);
- gMC->Gspart(114,"J/Psi",3,mass,0,tlife);
- fGPCode[kc][0]=114;
- fGPCode[kc][1]=114;
- // psi prime
- kc=Lucomp(30443);
- mass =GetPMAS(kc,1);
- tlife=GetPMAS(kc,4);
- gMC->Gspart(115,"Psi'",3,mass,0,tlife);
- fGPCode[kc][0]=115;
- fGPCode[kc][1]=115;
- // upsilon(1s)
- kc=Lucomp(553);
- mass =GetPMAS(kc,1);
- tlife=GetPMAS(kc,4);
- gMC->Gspart(116,"Upsilon",3,mass,0,tlife);
- fGPCode[kc][0]=116;
- fGPCode[kc][1]=116;
- // upsilon(2s)
- kc=Lucomp(30553);
- mass =GetPMAS(kc,1);
- tlife=GetPMAS(kc,4);
- gMC->Gspart(117,"Upsilon'",3,mass,0,tlife);
- fGPCode[kc][0]=117;
- fGPCode[kc][1]=117;
- // upsilon(3s)
- kc=Lucomp(30553);
- mass =GetPMAS(kc,1);
- tlife=GetPMAS(kc,4);
- gMC->Gspart(118,"Upsilon''",3,mass,0,tlife);
- fGPCode[kc][0]=118;
- fGPCode[kc][1]=118;
-//
-// charmed mesons
-//
- // D^+/-
- kc=Lucomp(411);
- mass =GetPMAS(kc,1);
- tlife=GetPMAS(kc,4);
- gMC->Gspart(119,"D^+",3,mass, 1,tlife);
- gMC->Gspart(120,"D^-",3,mass,-1,tlife);
- fGPCode[kc][0]=119;
- fGPCode[kc][1]=120;
- // D^0
- kc=Lucomp(421);
- mass =GetPMAS(kc,1);
- tlife=GetPMAS(kc,4);
- gMC->Gspart(121,"D^0",3,mass,0,tlife);
- gMC->Gspart(122,"D^0bar",3,mass,0,tlife);
- fGPCode[kc][0]=121;
- fGPCode[kc][1]=122;
- // D_s
- kc=Lucomp(431);
- mass =GetPMAS(kc,1);
- tlife=GetPMAS(kc,4);
- gMC->Gspart(123,"D_s^+",3,mass, 1,tlife);
- gMC->Gspart(124,"D_s^-",3,mass,-1,tlife);
- fGPCode[kc][0]=123;
- fGPCode[kc][1]=124;
- // Lambda_c
- kc=Lucomp(4122);
- mass =GetPMAS(kc,1);
- tlife=GetPMAS(kc,4);
- gMC->Gspart(125,"Lambda_c+",3,mass, 1,tlife);
- gMC->Gspart(126,"Lambda_c-",3,mass,-1,tlife);
- fGPCode[kc][0]=125;
- fGPCode[kc][1]=126;
- //
- // beauty mesons
- // B_0
- kc=Lucomp(511);
- mass =GetPMAS(kc,1);
- tlife=GetPMAS(kc,4);
- gMC->Gspart(127,"B^0",3,mass, 0,tlife);
- gMC->Gspart(128,"B^0bar",3,mass, 0,tlife);
- fGPCode[kc][0]=127;
- fGPCode[kc][1]=128;
- // B^+-
- kc=Lucomp(521);
- mass =GetPMAS(kc,1);
- tlife=GetPMAS(kc,4);
- gMC->Gspart(129,"B^+",3,mass, 1,tlife);
- gMC->Gspart(130,"B^-",3,mass,-1,tlife);
- fGPCode[kc][0]=129;
- fGPCode[kc][1]=130;
- // B_s
- kc=Lucomp(531);
- mass =GetPMAS(kc,1);
- tlife=GetPMAS(kc,4);
- gMC->Gspart(131,"B_s",3,mass, 0,tlife);
- gMC->Gspart(132,"B_s^bar",3,mass,0,tlife);
- fGPCode[kc][0]=131;
- fGPCode[kc][1]=132;
- // Lambda_b
- kc=Lucomp(5122);
- mass =GetPMAS(kc,1);
- tlife=GetPMAS(kc,4);
- gMC->Gspart(133,"Lambda_b",3,mass, 0,tlife);
- gMC->Gspart(134,"Lambda_b^bar",3,mass,0,tlife);
- fGPCode[kc][0]=133;
- fGPCode[kc][1]=134;
-//
-// set up correspondance between standard GEANT particle codes
-// and PYTHIA kf
-
- kc=Lucomp(22); // gamma
- fGPCode[kc][0]=1;
- fGPCode[kc][1]=1;
-
- kc=Lucomp(11); // positron
- fGPCode[kc][0]=2;
- fGPCode[kc][1]=3;
-
- kc=Lucomp(12); // neutrino
- fGPCode[kc][0]=4;
- fGPCode[kc][1]=4;
-
- kc=Lucomp(13); // muon
- fGPCode[kc][0]=5;
- fGPCode[kc][1]=6;
-
- kc=Lucomp(111); // pi0
- fGPCode[kc][0]=7;
- fGPCode[kc][1]=7;
-
- kc=Lucomp(211); // pi+
- fGPCode[kc][0]=8;
- fGPCode[kc][1]=9;
-
- kc=Lucomp(130); // K0 short
- fGPCode[kc][0]=10;
- fGPCode[kc][1]=10;
-
- kc=Lucomp(321); // K+/-
- fGPCode[kc][0]=11;
- fGPCode[kc][1]=12;
-
- kc=Lucomp(2112); // neutron/anti-neutron
- fGPCode[kc][0]=13;
- fGPCode[kc][1]=25;
-
- kc=Lucomp(2212); // proton/anti-proton
- fGPCode[kc][0]=14;
- fGPCode[kc][1]=15;
-
- kc=Lucomp(310); // K0 short
- fGPCode[kc][0]=16;
- fGPCode[kc][1]=16;
-
- kc=Lucomp(221); // eta
- fGPCode[kc][0]=17;
- fGPCode[kc][1]=17;
-
- kc=Lucomp(3122); // lambda
- fGPCode[kc][0]=18;
- fGPCode[kc][1]=18;
-
- kc=Lucomp(3222); // sigma+/antisigma+
- fGPCode[kc][0]=19;
- fGPCode[kc][1]=29;
-
- kc=Lucomp(3212); // sigma0/antisigma0
- fGPCode[kc][0]=20;
- fGPCode[kc][1]=28;
-
- kc=Lucomp(3112); // sigma-/antisigma-
- fGPCode[kc][0]=21;
- fGPCode[kc][1]=27;
-
- kc=Lucomp(3322); // xsi0-/antixsi0
- fGPCode[kc][0]=22;
- fGPCode[kc][1]=30;
-
- kc=Lucomp(3312); // xsi-/antixsi+
- fGPCode[kc][0]=23;
- fGPCode[kc][1]=31;
-
- kc=Lucomp(3334); // omega/antiomega
- fGPCode[kc][0]=24;
- fGPCode[kc][1]=32;
+ printf("\n Number of particles %d \n \n", np);
}
-
-
-Int_t AliPythia::GetGeantCode(Int_t kf)
+void AliPythia::ResetDecayTable()
{
- Int_t kc=Lucomp(TMath::Abs(kf));
- return (kf > 0) ? fGPCode[kc][0] : fGPCode[kc][1];
+// Set default values for pythia decay switches
+ Int_t i;
+ for (i = 1; i < 501; i++) SetMDCY(i,1,fDefMDCY[i]);
+ for (i = 1; i < 2001; i++) SetMDME(i,1,fDefMDME[i]);
}
-
-Float_t AliPythia::GetBraPart(Int_t kf)
+
+void AliPythia::SetDecayTable()
{
- Int_t kc=Lucomp(TMath::Abs(kf));
- return fBraPart[kc];
+// Set default values for pythia decay switches
+//
+ Int_t i;
+ for (i = 1; i < 501; i++) fDefMDCY[i] = GetMDCY(i,1);
+ for (i = 1; i < 2001; i++) fDefMDME[i] = GetMDME(i,1);
}
-
+#ifndef WIN32
+#define pyr pyr_
+#define pyrset pyrset_
+#define pyrget pyrget_
+#else
+#define pyr PYR
+#define pyrset PYRSET
+#define pyrget PYRGET
+#endif
+extern "C" {
+ Double_t pyr(Int_t*) {return sRandom->Rndm();}
+ void pyrset(Int_t*,Int_t*) {}
+ void pyrget(Int_t*,Int_t*) {}
+}