X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=PYTHIA6%2FAliPythia.cxx;h=9173c056a66f04a4b56d8566ac17ec5194b70794;hp=e1e055a76b8d59c95c5010e5088692d20f951e40;hb=2509792614ce13054e26ba4996fa3aaa88c3d525;hpb=1837e95c797311b0513dae0c3763ea0bb8cba9b6 diff --git a/PYTHIA6/AliPythia.cxx b/PYTHIA6/AliPythia.cxx index e1e055a76b8..9173c056a66 100644 --- a/PYTHIA6/AliPythia.cxx +++ b/PYTHIA6/AliPythia.cxx @@ -1,3 +1,4 @@ + /************************************************************************** * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * * @@ -17,9 +18,10 @@ #include "AliPythia.h" #include "AliPythiaRndm.h" -#include "../FASTSIM/AliFastGlauber.h" -#include "../FASTSIM/AliQuenchingWeights.h" +#include "AliFastGlauber.h" +#include "AliQuenchingWeights.h" #include "TVector3.h" +#include "PyquenCommon.h" ClassImp(AliPythia) @@ -28,11 +30,15 @@ ClassImp(AliPythia) # define pycell pycell_ # define pyshow pyshow_ # define pyrobo pyrobo_ +# define pyquen pyquen_ +# define pyevnw pyevnw_ # define type_of_call #else # define pyclus PYCLUS # define pycell PYCELL # define pyrobo PYROBO +# define pyquen PYQUEN +# define pyevnw PYEVNW # define type_of_call _stdcall #endif @@ -40,12 +46,23 @@ extern "C" void type_of_call pyclus(Int_t & ); extern "C" void type_of_call pycell(Int_t & ); extern "C" void type_of_call pyshow(Int_t &, Int_t &, Double_t &); extern "C" void type_of_call pyrobo(Int_t &, Int_t &, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &); +extern "C" void type_of_call pyquen(Double_t &, Int_t &, Double_t &); +extern "C" void type_of_call pyevnw(){;} //_____________________________________________________________________________ AliPythia* AliPythia::fgAliPythia=NULL; -AliPythia::AliPythia() +AliPythia::AliPythia(): + fProcess(kPyMb), + fEcms(0.), + fStrucFunc(kCTEQ5L), + fXJet(0.), + fYJet(0.), + fNGmax(30), + fZmax(0.97), + fGlauber(0), + fQuenchingWeights(0) { // Default Constructor // @@ -56,6 +73,23 @@ AliPythia::AliPythia() fQuenchingWeights = 0; } +AliPythia::AliPythia(const AliPythia& pythia): + TPythia6(pythia), + AliRndm(pythia), + fProcess(kPyMb), + fEcms(0.), + fStrucFunc(kCTEQ5L), + fXJet(0.), + fYJet(0.), + fNGmax(30), + fZmax(0.97), + fGlauber(0), + fQuenchingWeights(0) +{ + // Copy Constructor + pythia.Copy(*this); +} + void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfunc) { // Initialise the process to generate @@ -65,11 +99,23 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun fProcess = process; fEcms = energy; fStrucFunc = strucfunc; -// don't decay p0 - SetMDCY(Pycomp(111),1,0); -// select structure function +//...Switch off decay of pi0, K0S, Lambda, Sigma+-, Xi0-, Omega-. + SetMDCY(Pycomp(111) ,1,0); + SetMDCY(Pycomp(310) ,1,0); + SetMDCY(Pycomp(3122),1,0); + SetMDCY(Pycomp(3112),1,0); + SetMDCY(Pycomp(3212),1,0); + SetMDCY(Pycomp(3222),1,0); + SetMDCY(Pycomp(3312),1,0); + SetMDCY(Pycomp(3322),1,0); + SetMDCY(Pycomp(3334),1,0); + // Select structure function SetMSTP(52,2); SetMSTP(51,strucfunc); + // Particles produced in string fragmentation point directly to either of the two endpoints + // of the string (depending in the side they were generated from). + SetMSTU(16,2); + // // Pythia initialisation for selected processes// // @@ -81,13 +127,64 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun // select charm production switch (process) { + case kPyOldUEQ2ordered: //Old underlying events with Q2 ordered QCD processes +// Multiple interactions on. + SetMSTP(81,1); +// Double Gaussian matter distribution. + SetMSTP(82,4); + SetPARP(83,0.5); + SetPARP(84,0.4); +// pT0. + SetPARP(82,2.0); +// Reference energy for pT0 and energy rescaling pace. + SetPARP(89,1800); + SetPARP(90,0.25); +// String drawing almost completely minimizes string length. + SetPARP(85,0.9); + SetPARP(86,0.95); +// ISR and FSR activity. + SetPARP(67,4); + SetPARP(71,4); +// Lambda_FSR scale. + SetPARJ(81,0.29); + break; + case kPyOldUEQ2ordered2: +// Old underlying events with Q2 ordered QCD processes +// Multiple interactions on. + SetMSTP(81,1); +// Double Gaussian matter distribution. + SetMSTP(82,4); + SetPARP(83,0.5); + SetPARP(84,0.4); +// pT0. + SetPARP(82,2.0); +// Reference energy for pT0 and energy rescaling pace. + SetPARP(89,1800); + SetPARP(90,0.16); // here is the difference with kPyOldUEQ2ordered +// String drawing almost completely minimizes string length. + SetPARP(85,0.9); + SetPARP(86,0.95); +// ISR and FSR activity. + SetPARP(67,4); + SetPARP(71,4); +// Lambda_FSR scale. + SetPARJ(81,0.29); + break; + case kPyOldPopcorn: +// Old production mechanism: Old Popcorn + SetMSEL(1); + SetMSTJ(12,3); +// (D=2) Like MSTJ(12)=2 but added prod ofthe 1er rank baryon + SetMSTP(88,2); +// (D=1)see can be used to form baryons (BARYON JUNCTION) + SetMSTJ(1,1); + AtlasTuning(); + break; case kPyCharm: SetMSEL(4); -// // heavy quark masses SetPMAS(4,1,1.2); - SetMSTU(16,2); // // primordial pT SetMSTP(91,1); @@ -98,7 +195,6 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun case kPyBeauty: SetMSEL(5); SetPMAS(5,1,4.75); - SetMSTU(16,2); break; case kPyJpsi: SetMSEL(0); @@ -145,21 +241,44 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun SetMSUB(94,1); // double diffraction SetMSUB(95,1); // low pt production + AtlasTuning(); + break; + case kPyMbDefault: +// Minimum Bias pp-Collisions // -// ATLAS Tuning -// - SetMSTP(51,7); // CTEQ5L pdf +// +// select Pythia min. bias model + SetMSEL(0); + SetMSUB(92,1); // single diffraction AB-->XB + SetMSUB(93,1); // single diffraction AB-->AX + SetMSUB(94,1); // double diffraction + SetMSUB(95,1); // low pt production + + break; + case kPyLhwgMb: +// Les Houches Working Group 05 Minimum Bias pp-Collisions: hep-ph/0604120 +// -> Pythia 6.3 or above is needed +// + SetMSEL(0); + SetMSUB(92,1); // single diffraction AB-->XB + SetMSUB(93,1); // single diffraction AB-->AX + SetMSUB(94,1); // double diffraction + SetMSUB(95,1); // low pt production + + SetMSTP(51,kCTEQ6ll); // CTEQ6ll pdf + SetMSTP(52,2); + SetMSTP(68,1); + SetMSTP(70,2); SetMSTP(81,1); // Multiple Interactions ON SetMSTP(82,4); // Double Gaussian Model + SetMSTP(88,1); - SetPARP(82,1.8); // [GeV] PT_min at Ref. energy - SetPARP(89,1000.); // [GeV] Ref. energy - SetPARP(90,0.16); // 2*epsilon (exponent in power law) + SetPARP(82,2.3); // [GeV] PT_min at Ref. energy SetPARP(83,0.5); // Core density in proton matter distribution (def.value) SetPARP(84,0.5); // Core radius - SetPARP(85,0.33); // Regulates gluon prod. mechanism - SetPARP(86,0.66); // Regulates gluon prod. mechanism - SetPARP(67,1); // Regulates Initial State Radiation + SetPARP(85,0.9); // Regulates gluon prod. mechanism + SetPARP(90,0.2); // 2*epsilon (exponent in power law) + break; case kPyMbNonDiffr: // Minimum Bias pp-Collisions @@ -169,34 +288,44 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun SetMSEL(0); SetMSUB(95,1); // low pt production -// -// ATLAS Tuning -// - - SetMSTP(51,7); // CTEQ5L pdf - SetMSTP(81,1); // Multiple Interactions ON - SetMSTP(82,4); // Double Gaussian Model - - SetPARP(82,1.8); // [GeV] PT_min at Ref. energy - SetPARP(89,1000.); // [GeV] Ref. energy - SetPARP(90,0.16); // 2*epsilon (exponent in power law) - SetPARP(83,0.5); // Core density in proton matter distribution (def.value) - SetPARP(84,0.5); // Core radius - SetPARP(85,0.33); // Regulates gluon prod. mechanism - SetPARP(86,0.66); // Regulates gluon prod. mechanism - SetPARP(67,1); // Regulates Initial State Radiation + AtlasTuning(); + break; + case kPyMbMSEL1: + ConfigHeavyFlavor(); +// Intrinsic + SetMSTP(91,1);// Width (1=gaussian) primordial kT dist. inside hadrons + SetPARP(91,1.); // = PARP(91,1.)^2 + SetPARP(93,5.); // Upper cut-off +// Set Q-quark mass + SetPMAS(4,1,1.2); // Charm quark mass + SetPMAS(5,1,4.78); // Beauty quark mass + SetPARP(71,4.); // Defaut value +// Atlas Tuning + AtlasTuning(); break; case kPyJets: // // QCD Jets // SetMSEL(1); - break; + // Pythia Tune A (CDF) + // + SetPARP(67,2.5); // Regulates Initial State Radiation (value from best fit to D0 dijet analysis) + SetMSTP(82,4); // Double Gaussian Model + SetPARP(82,2.0); // [GeV] PT_min at Ref. energy + SetPARP(84,0.4); // Core radius + SetPARP(85,0.90) ; // Regulates gluon prod. mechanism + SetPARP(86,0.95); // Regulates gluon prod. mechanism + SetPARP(89,1800.); // [GeV] Ref. energy + SetPARP(90,0.25); // 2*epsilon (exponent in power law) + break; case kPyDirectGamma: SetMSEL(10); break; case kPyCharmPbPbMNR: case kPyD0PbPbMNR: + case kPyDPlusPbPbMNR: + case kPyDPlusStrangePbPbMNR: // Tuning of Pythia parameters aimed to get a resonable agreement // between with the NLO calculation by Mangano, Nason, Ridolfi for the // c-cbar single inclusive and double differential distributions. @@ -204,37 +333,18 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs. // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard) // has to be set to 2.1GeV. Example in ConfigCharmPPR.C. - - // All QCD processes - SetMSEL(1); - - // No multiple interactions - SetMSTP(81,0); - SetPARP(81,0.0); - SetPARP(82,0.0); - - // Initial/final parton shower on (Pythia default) - SetMSTP(61,1); - SetMSTP(71,1); - - // 2nd order alpha_s - SetMSTP(2,2); - - // QCD scales - SetMSTP(32,2); - SetPARP(34,1.0); - + ConfigHeavyFlavor(); // Intrinsic SetMSTP(91,1); SetPARP(91,1.304); SetPARP(93,6.52); - // Set c-quark mass SetPMAS(4,1,1.2); - break; case kPyCharmpPbMNR: case kPyD0pPbMNR: + case kPyDPluspPbMNR: + case kPyDPlusStrangepPbMNR: // Tuning of Pythia parameters aimed to get a resonable agreement // between with the NLO calculation by Mangano, Nason, Ridolfi for the // c-cbar single inclusive and double differential distributions. @@ -242,37 +352,19 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs. // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard) // has to be set to 2.1GeV. Example in ConfigCharmPPR.C. - - // All QCD processes - SetMSEL(1); - - // No multiple interactions - SetMSTP(81,0); - SetPARP(81,0.0); - SetPARP(82,0.0); - - // Initial/final parton shower on (Pythia default) - SetMSTP(61,1); - SetMSTP(71,1); - - // 2nd order alpha_s - SetMSTP(2,2); - - // QCD scales - SetMSTP(32,2); - SetPARP(34,1.0); - + ConfigHeavyFlavor(); // Intrinsic - SetMSTP(91,1); - SetPARP(91,1.16); - SetPARP(93,5.8); - + SetMSTP(91,1); + SetPARP(91,1.16); + SetPARP(93,5.8); + // Set c-quark mass - SetPMAS(4,1,1.2); - + SetPMAS(4,1,1.2); break; case kPyCharmppMNR: case kPyD0ppMNR: + case kPyDPlusppMNR: + case kPyDPlusStrangeppMNR: // Tuning of Pythia parameters aimed to get a resonable agreement // between with the NLO calculation by Mangano, Nason, Ridolfi for the // c-cbar single inclusive and double differential distributions. @@ -280,35 +372,42 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs. // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard) // has to be set to 2.1GeV. Example in ConfigCharmPPR.C. - - // All QCD processes - SetMSEL(1); - - // No multiple interactions - SetMSTP(81,0); - SetPARP(81,0.0); - SetPARP(82,0.0); - - // Initial/final parton shower on (Pythia default) - SetMSTP(61,1); - SetMSTP(71,1); - - // 2nd order alpha_s - SetMSTP(2,2); - - // QCD scales - SetMSTP(32,2); - SetPARP(34,1.0); - + ConfigHeavyFlavor(); // Intrinsic - SetMSTP(91,1); - SetPARP(91,1.); - SetPARP(93,5.); - + SetMSTP(91,1); + SetPARP(91,1.); + SetPARP(93,5.); + // Set c-quark mass - SetPMAS(4,1,1.2); - + SetPMAS(4,1,1.2); break; + case kPyCharmppMNRwmi: + // Tuning of Pythia parameters aimed to get a resonable agreement + // between with the NLO calculation by Mangano, Nason, Ridolfi for the + // c-cbar single inclusive and double differential distributions. + // This parameter settings are meant to work with pp collisions + // and with kCTEQ5L PDFs. + // Added multiple interactions according to ATLAS tune settings. + // To get a "reasonable" agreement with MNR results, events have to be + // generated with the minimum ptHard (AliGenPythia::SetPtHard) + // set to 2.76 GeV. + // To get a "perfect" agreement with MNR results, events have to be + // generated in four ptHard bins with the following relative + // normalizations: + // 2.76-3 GeV: 25% + // 3-4 GeV: 40% + // 4-8 GeV: 29% + // >8 GeV: 6% + ConfigHeavyFlavor(); + // Intrinsic + SetMSTP(91,1); + SetPARP(91,1.); + SetPARP(93,5.); + + // Set c-quark mass + SetPMAS(4,1,1.2); + AtlasTuning(); + break; case kPyBeautyPbPbMNR: // Tuning of Pythia parameters aimed to get a resonable agreement // between with the NLO calculation by Mangano, Nason, Ridolfi for the @@ -317,36 +416,16 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs. // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard) // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C. - - // All QCD processes - SetMSEL(1); - - // No multiple interactions - SetMSTP(81,0); - SetPARP(81,0.0); - SetPARP(82,0.0); - - // Initial/final parton shower on (Pythia default) - SetMSTP(61,1); - SetMSTP(71,1); - - // 2nd order alpha_s - SetMSTP(2,2); - + ConfigHeavyFlavor(); // QCD scales - SetMSTP(32,2); - SetPARP(34,1.0); - SetPARP(67,1.0); - SetPARP(71,1.0); - + SetPARP(67,1.0); + SetPARP(71,1.0); // Intrinsic - SetMSTP(91,1); - SetPARP(91,2.035); - SetPARP(93,10.17); - + SetMSTP(91,1); + SetPARP(91,2.035); + SetPARP(93,10.17); // Set b-quark mass - SetPMAS(5,1,4.75); - + SetPMAS(5,1,4.75); break; case kPyBeautypPbMNR: // Tuning of Pythia parameters aimed to get a resonable agreement @@ -356,36 +435,16 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs. // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard) // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C. - - // All QCD processes - SetMSEL(1); - - // No multiple interactions - SetMSTP(81,0); - SetPARP(81,0.0); - SetPARP(82,0.0); - - // Initial/final parton shower on (Pythia default) - SetMSTP(61,1); - SetMSTP(71,1); - - // 2nd order alpha_s - SetMSTP(2,2); - + ConfigHeavyFlavor(); // QCD scales - SetMSTP(32,2); - SetPARP(34,1.0); - SetPARP(67,1.0); - SetPARP(71,1.0); - + SetPARP(67,1.0); + SetPARP(71,1.0); // Intrinsic - SetMSTP(91,1); - SetPARP(91,1.60); - SetPARP(93,8.00); - + SetMSTP(91,1); + SetPARP(91,1.60); + SetPARP(93,8.00); // Set b-quark mass - SetPMAS(5,1,4.75); - + SetPMAS(5,1,4.75); break; case kPyBeautyppMNR: // Tuning of Pythia parameters aimed to get a resonable agreement @@ -395,44 +454,105 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs. // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard) // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C. - - // All QCD processes - SetMSEL(1); - - // No multiple interactions - SetMSTP(81,0); - SetPARP(81,0.0); - SetPARP(82,0.0); - - // Initial/final parton shower on (Pythia default) - SetMSTP(61,1); - SetMSTP(71,1); - - // 2nd order alpha_s - SetMSTP(2,2); - + ConfigHeavyFlavor(); // QCD scales - SetMSTP(32,2); - SetPARP(34,1.0); - SetPARP(67,1.0); - SetPARP(71,1.0); - - // Intrinsic - SetMSTP(91,1); - SetPARP(91,1.); - SetPARP(93,5.); + SetPARP(67,1.0); + SetPARP(71,1.0); + + // Intrinsic + SetMSTP(91,1); + SetPARP(91,1.); + SetPARP(93,5.); + + // Set b-quark mass + SetPMAS(5,1,4.75); + break; + case kPyBeautyppMNRwmi: + // Tuning of Pythia parameters aimed to get a resonable agreement + // between with the NLO calculation by Mangano, Nason, Ridolfi for the + // b-bbar single inclusive and double differential distributions. + // This parameter settings are meant to work with pp collisions + // and with kCTEQ5L PDFs. + // Added multiple interactions according to ATLAS tune settings. + // To get a "reasonable" agreement with MNR results, events have to be + // generated with the minimum ptHard (AliGenPythia::SetPtHard) + // set to 2.76 GeV. + // To get a "perfect" agreement with MNR results, events have to be + // generated in four ptHard bins with the following relative + // normalizations: + // 2.76-4 GeV: 5% + // 4-6 GeV: 31% + // 6-8 GeV: 28% + // >8 GeV: 36% + ConfigHeavyFlavor(); + // QCD scales + SetPARP(67,1.0); + SetPARP(71,1.0); + + // Intrinsic + SetMSTP(91,1); + SetPARP(91,1.); + SetPARP(93,5.); // Set b-quark mass - SetPMAS(5,1,4.75); + SetPMAS(5,1,4.75); + + AtlasTuning(); + break; + case kPyW: + + //Inclusive production of W+/- + SetMSEL(0); + //f fbar -> W+ + SetMSUB(2,1); + // //f fbar -> g W+ + // SetMSUB(16,1); + // //f fbar -> gamma W+ + // SetMSUB(20,1); + // //f g -> f W+ + // SetMSUB(31,1); + // //f gamma -> f W+ + // SetMSUB(36,1); + + // Initial/final parton shower on (Pythia default) + // With parton showers on we are generating "W inclusive process" + SetMSTP(61,1); //Initial QCD & QED showers on + SetMSTP(71,1); //Final QCD & QED showers on + + break; + + case kPyZ: + + //Inclusive production of Z + SetMSEL(0); + //f fbar -> Z/gamma + SetMSUB(1,1); + + // // f fbar -> g Z/gamma + // SetMSUB(15,1); + // // f fbar -> gamma Z/gamma + // SetMSUB(19,1); + // // f g -> f Z/gamma + // SetMSUB(30,1); + // // f gamma -> f Z/gamma + // SetMSUB(35,1); + + //only Z included, not gamma + SetMSTP(43,2); + + // Initial/final parton shower on (Pythia default) + // With parton showers on we are generating "Z inclusive process" + SetMSTP(61,1); //Initial QCD & QED showers on + SetMSTP(71,1); //Final QCD & QED showers on + + break; - break; } // // Initialize PYTHIA SetMSTP(41,1); // all resonance decays switched on - Initialize("CMS","p","p",fEcms); - + } Int_t AliPythia::CheckedLuComp(Int_t kf) @@ -541,24 +661,23 @@ void AliPythia::Pyrobo(Int_t imi, Int_t ima, Double_t the, Double_t phi, Double_ -void AliPythia::InitQuenching(Float_t cMin, Float_t cMax, Float_t qTransport, Float_t maxLength, Int_t iECMethod) +void AliPythia::InitQuenching(Float_t cMin, Float_t cMax, Float_t k, Int_t iECMethod, Float_t zmax, Int_t ngmax) { // Initializes // (1) The quenching model using quenching weights according to C. Salgado and U. Wiedemann // (2) The nuclear geometry using the Glauber Model // - - + fGlauber = new AliFastGlauber(); fGlauber->Init(2); fGlauber->SetCentralityClass(cMin, cMax); fQuenchingWeights = new AliQuenchingWeights(); fQuenchingWeights->InitMult(); - fQuenchingWeights->SetQTransport(qTransport); + fQuenchingWeights->SetK(k); fQuenchingWeights->SetECMethod(AliQuenchingWeights::kECMethod(iECMethod)); - fQuenchingWeights->SetLengthMax(Int_t(maxLength)); - fQuenchingWeights->SampleEnergyLoss(); + fNGmax = ngmax; + fZmax = zmax; } @@ -585,95 +704,129 @@ void AliPythia::Quench() static Float_t eMean = 0.; static Int_t icall = 0; - Double_t p0[2][5]; - Double_t p1[2][5]; - Double_t p2[2][5]; - Int_t klast[2] = {-1, -1}; + Double_t p0[4][5]; + Double_t p1[4][5]; + Double_t p2[4][5]; + Int_t klast[4] = {-1, -1, -1, -1}; Int_t numpart = fPyjets->N; - Double_t px = 0., py = 0., pz = 0., e = 0., m = 0., p = 0., pt = 0., theta = 0.; - Double_t pxq[2], pyq[2], pzq[2], eq[2], yq[2], mq[2], pq[2], phiq[2], thetaq[2], ptq[2]; - Bool_t quenched[2]; - Double_t phi; - Double_t zInitial[2], wjtKick[2]; - Int_t nGluon[2]; - + Double_t px = 0., py = 0., pz = 0., e = 0., m = 0., p = 0., pt = 0., theta = 0., phi = 0.; + Double_t pxq[4], pyq[4], pzq[4], eq[4], yq[4], mq[4], pq[4], phiq[4], thetaq[4], ptq[4]; + Bool_t quenched[4]; + Double_t wjtKick[4]; + Int_t nGluon[4]; + Int_t qPdg[4]; Int_t imo, kst, pdg; + // -// Primary partons +// Sore information about Primary partons +// +// j = +// 0, 1 partons from hard scattering +// 2, 3 partons from initial state radiation +// + for (Int_t i = 2; i <= 7; i++) { + Int_t j = 0; + // Skip gluons that participate in hard scattering + if (i == 4 || i == 5) continue; + // Gluons from hard Scattering + if (i == 6 || i == 7) { + j = i - 6; + pxq[j] = fPyjets->P[0][i]; + pyq[j] = fPyjets->P[1][i]; + pzq[j] = fPyjets->P[2][i]; + eq[j] = fPyjets->P[3][i]; + mq[j] = fPyjets->P[4][i]; + } else { + // Gluons from initial state radiation + // + // Obtain 4-momentum vector from difference between original parton and parton after gluon + // radiation. Energy is calculated independently because initial state radition does not + // conserve strictly momentum and energy for each partonic system independently. + // + // Not very clean. Should be improved ! + // + // + j = i; + pxq[j] = fPyjets->P[0][i] - fPyjets->P[0][i+2]; + pyq[j] = fPyjets->P[1][i] - fPyjets->P[1][i+2]; + pzq[j] = fPyjets->P[2][i] - fPyjets->P[2][i+2]; + mq[j] = fPyjets->P[4][i]; + eq[j] = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j] + pzq[j] * pzq[j] + mq[j] * mq[j]); + } +// +// Calculate some kinematic variables // - - - - for (Int_t i = 6; i <= 7; i++) { - Int_t j = i - 6; - - pxq[j] = fPyjets->P[0][i]; - pyq[j] = fPyjets->P[1][i]; - pzq[j] = fPyjets->P[2][i]; - eq[j] = fPyjets->P[3][i]; - mq[j] = fPyjets->P[4][i]; yq[j] = 0.5 * TMath::Log((eq[j] + pzq[j] + 1.e-14) / (eq[j] - pzq[j] + 1.e-14)); pq[j] = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j] + pzq[j] * pzq[j]); phiq[j] = TMath::Pi()+TMath::ATan2(-pyq[j], -pxq[j]); ptq[j] = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j]); thetaq[j] = TMath::ATan2(ptq[j], pzq[j]); - phi = phiq[j]; - - // Quench only central jets - if (TMath::Abs(yq[j]) > 2.5) { - zInitial[j] = 0.; + qPdg[j] = fPyjets->K[1][i]; + } + + Double_t int0[4]; + Double_t int1[4]; + + fGlauber->GetI0I1ForPythiaAndXY(4, phiq, int0, int1, fXJet, fYJet, 15.); + + for (Int_t j = 0; j < 4; j++) { + // + // Quench only central jets and with E > 10. + // + + + Int_t itype = (qPdg[j] == 21) ? 2 : 1; + Double_t eloss = fQuenchingWeights->GetELossRandomKFast(itype, int0[j], int1[j], eq[j]); + + if (TMath::Abs(yq[j]) > 2.5 || eq[j] < 10.) { + fZQuench[j] = 0.; } else { - pdg = fPyjets->K[1][i]; - - // Get length in nucleus - Double_t l; - fGlauber->GetLengthsForPythia(1, &phi, &l, -1.); - // - // Energy loss for given length and parton typr - Int_t itype = (pdg == 21) ? 2 : 1; - - Double_t eloss = fQuenchingWeights->GetELossRandom(itype, l, eq[j]); - if (eq[j] > 80. && TMath::Abs(yq[j]) < 0.5) { + if (eq[j] > 40. && TMath::Abs(yq[j]) < 0.5) { icall ++; eMean += eloss; } - // // Extra pt - - wjtKick[j] = TMath::Sqrt(l * fQuenchingWeights->GetQTransport()); + Double_t l = fQuenchingWeights->CalcLk(int0[j], int1[j]); + wjtKick[j] = TMath::Sqrt(l * fQuenchingWeights->CalcQk(int0[j], int1[j])); // // Fractional energy loss - zInitial[j] = eloss / eq[j]; + fZQuench[j] = eloss / eq[j]; // // Avoid complete loss // - if (zInitial[j] == 1.) zInitial[j] = 0.95; + if (fZQuench[j] == 1.) fZQuench[j] = fZmax; // // Some debug printing - 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", - j, itype, eq[j], phi, l, eloss, wjtKick[j], eMean / Float_t(icall+1), yq[j]); + + +// 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", +// j, itype, eq[j], phiq[j], l, eloss, wjtKick[j], eMean / Float_t(icall+1), yq[j]); - zInitial[j] = 1.; - while (zInitial[j] >= 0.95) zInitial[j] = gRandom->Exp(0.2); +// fZQuench[j] = 0.8; +// while (fZQuench[j] >= 0.95) fZQuench[j] = gRandom->Exp(0.2); } - quenched[j] = (zInitial[j] > 0.01); + quenched[j] = (fZQuench[j] > 0.01); } // primary partons - + + + Double_t pNew[1000][4]; Int_t kNew[1000]; Int_t icount = 0; + Double_t zquench[4]; + // // System Loop - for (Int_t isys = 0; isys < 2; isys++) { + for (Int_t isys = 0; isys < 4; isys++) { // Skip to next system if not quenched. if (!quenched[isys]) continue; - nGluon[isys] = 1 + Int_t(zInitial[isys] / (1. - zInitial[isys])); - if (nGluon[isys] > 6) nGluon[isys] = 6; - zInitial[isys] = 1. - TMath::Power(1. - zInitial[isys], 1./Double_t(nGluon[isys])); + nGluon[isys] = 1 + Int_t(fZQuench[isys] / (1. - fZQuench[isys])); + if (nGluon[isys] > fNGmax) nGluon[isys] = fNGmax; + zquench[isys] = 1. - TMath::Power(1. - fZQuench[isys], 1./Double_t(nGluon[isys])); wjtKick[isys] = wjtKick[isys] / TMath::Sqrt(Double_t(nGluon[isys])); @@ -686,9 +839,6 @@ void AliPythia::Quench() // Loop on radiation events for (Int_t iglu = 0; iglu < nGluon[isys]; iglu++) { - Double_t zHeavy = zInitial[isys]; -// - while (1) { icount = 0; for (Int_t k = 0; k < 4; k++) @@ -709,8 +859,12 @@ void AliPythia::Quench() // Quarks and gluons only if (pdg != 21 && TMath::Abs(pdg) > 6) continue; // Particles from hard scattering only + if (imo > 8 && imo < 1000) imo = fPyjets->K[2][imo - 1]; - if (imo != (isys + 7) && (imo % 1000) != (isys + 7)) continue; + Int_t imom = imo % 1000; + if ((isys == 0 || isys == 1) && ((imom != (isys + 7)))) continue; + if ((isys == 2 || isys == 3) && ((imom != (isys + 1)))) continue; + // Skip comment lines if (kst != 1 && kst != 2) continue; @@ -727,9 +881,8 @@ void AliPythia::Quench() theta = TMath::ATan2(pt, pz); // -// Save 4-momentum sum for balancing - Int_t index = imo - 7; - if (index >= 1000) index = imo % 1000 - 7; +// Save 4-momentum sum for balancing + Int_t index = isys; p0[index][0] += px; p0[index][1] += py; @@ -740,7 +893,8 @@ void AliPythia::Quench() // // Fractional energy loss - Double_t z = zInitial[index]; + Double_t z = zquench[index]; + // Don't fully quench radiated gluons // @@ -748,12 +902,12 @@ void AliPythia::Quench() // This small factor makes sure that the gluons are not too close in phase space to avoid recombination // - z = 0.05; + z = 0.02; } +// printf("z: %d %f\n", imo, z); + // - - if (m > 0.) z = zHeavy; // // @@ -854,7 +1008,7 @@ void AliPythia::Quench() p2[isys][4] = TMath::Sqrt(p2[isys][4]); break; } else { - printf("Warning negative mass squared in system %d %f ! \n", isys, zInitial[isys]); + printf("Warning negative mass squared in system %d %f ! \n", isys, zquench[isys]); 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]); if (p2[isys][4] < -0.01) { printf("Negative mass squared !\n"); @@ -1051,3 +1205,99 @@ void AliPythia::Quench() } // this->Pylist(1); } // end quench + + +void AliPythia::Pyquen(Double_t a, Int_t ibf, Double_t b) +{ + // Igor Lokthine's quenching routine + // http://lokhtin.web.cern.ch/lokhtin/pyquen/pyquen.txt + + pyquen(a, ibf, b); +} + +void AliPythia::SetPyquenParameters(Double_t t0, Double_t tau0, Int_t nf, Int_t iengl, Int_t iangl) +{ + // Set the parameters for the PYQUEN package. + // See comments in PyquenCommon.h + + + PYQPAR.t0 = t0; + PYQPAR.tau0 = tau0; + PYQPAR.nf = nf; + PYQPAR.iengl = iengl; + PYQPAR.iangl = iangl; +} + + +void AliPythia::Pyevnw() +{ + // New multiple interaction scenario + pyevnw(); +} + +void AliPythia::GetQuenchingParameters(Double_t& xp, Double_t& yp, Double_t z[4]) +{ + // Return event specific quenching parameters + xp = fXJet; + yp = fYJet; + for (Int_t i = 0; i < 4; i++) z[i] = fZQuench[i]; + +} + +void AliPythia::ConfigHeavyFlavor() +{ + // + // Default configuration for Heavy Flavor production + // + // All QCD processes + // + SetMSEL(1); + + // No multiple interactions + SetMSTP(81,0); + SetPARP(81, 0.); + SetPARP(82, 0.); + // Initial/final parton shower on (Pythia default) + SetMSTP(61,1); + SetMSTP(71,1); + + // 2nd order alpha_s + SetMSTP(2,2); + + // QCD scales + SetMSTP(32,2); + SetPARP(34,1.0); +} + +void AliPythia::AtlasTuning() +{ + // + // Configuration for the ATLAS tuning + SetMSTP(51, kCTEQ5L); // CTEQ5L pdf + SetMSTP(81,1); // Multiple Interactions ON + SetMSTP(82,4); // Double Gaussian Model + SetPARP(81,1.9); // Min. pt for multiple interactions (default in 6.2-14) + SetPARP(82,1.8); // [GeV] PT_min at Ref. energy + SetPARP(89,1000.); // [GeV] Ref. energy + SetPARP(90,0.16); // 2*epsilon (exponent in power law) + SetPARP(83,0.5); // Core density in proton matter distribution (def.value) + SetPARP(84,0.5); // Core radius + SetPARP(85,0.33); // Regulates gluon prod. mechanism + SetPARP(86,0.66); // Regulates gluon prod. mechanism + SetPARP(67,1); // Regulates Initial State Radiation +} + +AliPythia& AliPythia::operator=(const AliPythia& rhs) +{ +// Assignment operator + rhs.Copy(*this); + return *this; +} + + void AliPythia::Copy(TObject&) const +{ + // + // Copy + // + Fatal("Copy","Not implemented!\n"); +}