X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=EVGEN%2FAliPythia.cxx;h=c39fa5787696ecee168178bc10c88ccdedab5db5;hb=8d29edbe2698055fed6dff133ad33b5fed426e12;hp=5e060345d133813da89b7df9228561b3358b8848;hpb=fe4da5cc22f890b04843f1aebec0f1bf4f9c3fc9;p=u%2Fmrichter%2FAliRoot.git diff --git a/EVGEN/AliPythia.cxx b/EVGEN/AliPythia.cxx index 5e060345d13..c39fa578769 100644 --- a/EVGEN/AliPythia.cxx +++ b/EVGEN/AliPythia.cxx @@ -1,45 +1,78 @@ -#include "AliPythia.h" -#include "AliMC.h" -ClassImp(AliPythia) +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ -#ifndef WIN32 -# define lu1ent lu1ent_ -# define type_of_call -#else -# define lu1ent LU1ENT -# define type_of_call _stdcall -#endif +/* +$Log$ +Revision 1.13 2000/12/18 08:55:35 morsch +Make AliPythia dependent generartors work with new scheme of random number generation -extern "C" void type_of_call - lu1ent(Int_t&, Int_t&, Float_t&, Float_t&, Float_t&); +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. -Int_t AliPythia::fgInit=0; +Revision 1.8 2000/09/06 14:26:24 morsch +Decayer functionality of AliPythia has been moved to AliDecayerPythia. +Class is now a singleton. -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); +Revision 1.7 2000/06/09 20:34:50 morsch +All coding rule violations except RS3 corrected -} -void AliPythia::DecayParticle(Int_t idpart, - Float_t mom, Float_t theta,Float_t phi) +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 + +Revision 1.4 1999/09/29 09:24:14 fca +Introduction of the Copyright and cvs Log + +*/ + + +#include "AliPythia.h" + +ClassImp(AliPythia) + +//_____________________________________________________________________________ + +AliPythia* AliPythia::fgAliPythia=NULL; + +AliPythia::AliPythia() { - Lu1Ent(0, idpart, mom, theta, phi); - GetPrimaries(); +// Default Constructor +// +// Set random number + if (!sRandom) sRandom=fRandom; } 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); @@ -54,7 +87,7 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun // select charm production switch (process) { - case charm: + case kPyCharm: SetMSEL(4); // // heavy quark masses @@ -64,20 +97,20 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun // // 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); @@ -87,7 +120,7 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun SetMSUB(88,1); // gg-> chi_2c g SetMSUB(89,1); - case charm_unforced: + case kPyCharmUnforced: SetMSEL(0); // gq->qg SetMSUB(28,1); @@ -95,7 +128,7 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun SetMSUB(53,1); // gg->gg SetMSUB(68,1); - case beauty_unforced: + case kPyBeautyUnforced: SetMSEL(0); // gq->qg SetMSUB(28,1); @@ -104,398 +137,100 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun // gg->gg SetMSUB(68,1); break; + 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); + SetMSUB(92,1); + SetMSUB(93,1); + SetMSUB(94,1); + SetMSUB(95,1); +// Multiple interactions switched on + SetMSTP(81,1); + SetMSTP(82,1); +// Low-pT cut-off for hard scattering + SetPARP(81,1.9); +// model for subsequent non-hardest interaction +// 90% gg->gg 10% gg->qq + 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 - Initialize("CMS","p","p",fEcms); -} + SetMSTP(41,1); -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; + Initialize("CMS","p","p",fEcms); } -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); - 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: - ForceParticleDecay( 411,13,1); // D+/- - ForceParticleDecay( 421,13,1); // D0 - ForceParticleDecay( 431,13,1); // D_s - ForceParticleDecay( 4122,13,1); // Lambda_c - 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; +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; - AliMC *pMC=AliMC::GetMC(); - - 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)); - } - - for (i=1;i<=4;i++) { - SetPMAS(nkc,i,GetPMAS(kc,i)); - } - 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 dileptin decay-channel - kc=LuComp(41); - mass =GetPMAS(kc,1); - tlife=GetPMAS(kc,4); - pMC->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); - pMC->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); - pMC->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); - pMC->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); - pMC->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); - pMC->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); - pMC->Gspart(119,"D^+",3,mass, 1,tlife); - pMC->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); - pMC->Gspart(121,"D^0",3,mass,0,tlife); - pMC->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); - pMC->Gspart(123,"D_s^+",3,mass, 1,tlife); - pMC->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); - pMC->Gspart(125,"Lambda_c+",3,mass, 1,tlife); - pMC->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); - pMC->Gspart(127,"B^0",3,mass, 0,tlife); - pMC->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); - pMC->Gspart(129,"B^+",3,mass, 1,tlife); - pMC->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); - pMC->Gspart(131,"B_s",3,mass, 0,tlife); - pMC->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); - pMC->Gspart(133,"Lambda_b",3,mass, 0,tlife); - pMC->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; -} +#ifndef WIN32 +#define pyr pyr_ +#define pyrset pyrset_ +#define pyrget pyrget_ +#else +#define pyr PYR +#define pyrset PYRSET +#define pyrget PYRGET +#endif - -Int_t AliPythia::GetGeantCode(Int_t kf) -{ - Int_t kc=LuComp(TMath::Abs(kf)); - return (kf > 0) ? fGPCode[kc][0] : fGPCode[kc][1]; +extern "C" { + Double_t pyr(Int_t*) {return sRandom->Rndm();} + void pyrset(Int_t*,Int_t*) {} + void pyrget(Int_t*,Int_t*) {} } - -Float_t AliPythia::GetBraPart(Int_t kf) -{ - Int_t kc=LuComp(TMath::Abs(kf)); - return fBraPart[kc]; -} - - - -