X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=EVGEN%2FAliGenMUONCocktailpp.cxx;h=e8fe9c3fd1987a1a1bc34ff1caf73ff5156c0bac;hb=1f0f39e23ac3a7fcbf0f995a6e8b41f636df8483;hp=b8f95538fc60eb668b30cafe5ef5640c4aa0ffac;hpb=b44c390145ebd4fb7c11ba9e29dc5152ce09f1ea;p=u%2Fmrichter%2FAliRoot.git diff --git a/EVGEN/AliGenMUONCocktailpp.cxx b/EVGEN/AliGenMUONCocktailpp.cxx index b8f95538fc6..e8fe9c3fd19 100644 --- a/EVGEN/AliGenMUONCocktailpp.cxx +++ b/EVGEN/AliGenMUONCocktailpp.cxx @@ -24,6 +24,17 @@ // production cross-sections in pp collisions // July 07:added heavy quark production from AliGenCorrHF and heavy quark // production switched off in Pythia +// Aug. 07: added trigger cut on total momentum +// 2009: added possibility to hide x-sections (B. Vulpescu) +// 2009: added possibility to have the cocktail (fast generator and param.) +// for pp @ 10 TeV or pp @ 14 TeV (N. Bastid) +//----------------- +// 2009: added polarization (L. Bianchi) +//------------------ +// 11/2009: added chi_c1 & chi_c2 (P.Crochet & N.Bastid). +// Cross-sections for charmonia are now directly taken from the Yellow Report +// (hep-ph/0311048) Tab.9, page 19. See below for details w.r.t. beam energy. +// usage: see example of Config in $ALICE_ROOT/prod/LHC09a10/Config.C #include #include @@ -41,8 +52,8 @@ #include "AliStack.h" #include "AliDecayer.h" #include "AliLog.h" -#include "AliGenPythia.h" #include "AliGenCorrHF.h" +#include "AliDecayerPolarized.h" ClassImp(AliGenMUONCocktailpp) @@ -52,14 +63,70 @@ AliGenMUONCocktailpp::AliGenMUONCocktailpp() fDecayer(0), fDecayModeResonance(kAll), fDecayModePythia(kAll), - fTotalRate(0), fMuonMultiplicity(0), fMuonPtCut(0.), + fMuonPCut(0.), fMuonThetaMinCut(0.), - fMuonThetaMaxCut(0.), + fMuonThetaMaxCut(180.), fMuonOriginCut(-999.), fNSucceded(0), - fNGenerated(0) + fNGenerated(0), + + fJpsiPol(0), + fChic1Pol(0), + fChic2Pol(0), + fPsiPPol(0), + fUpsPol(0), + fUpsPPol(0), + fUpsPPPol(0), + fPolFrame(0), + +// x-sections for pp @ 7 TeV: charmonia from hep-ph/0311048 Tab.9, page 19, +// bottomnium as for 10 TeV +/* fCMSEnergy(7), + fSigmaReaction(0.0695), + fSigmaJPsi(21.8e-6), + fSigmaChic1(21.1e-6), + fSigmaChic2(34.9e-6), + fSigmaPsiP(4.93e-6), + fSigmaUpsilon(0.463e-6), + fSigmaUpsilonP(0.154e-6), + fSigmaUpsilonPP(0.0886e-6), + fSigmaCCbar(6.91e-3), + fSigmaBBbar(0.232e-3), +*/ + +//x-sections for pp @ 10 TeV: charmonia and bottomonia from 14 TeV numbers +// scaled down according to ccbar and bbar cross-sections + fCMSEnergy(10), + fSigmaReaction(0.0695), + fSigmaJPsi(26.06e-6), + fSigmaChic1(25.18e-6), + fSigmaChic2(41.58e-6), + fSigmaPsiP(5.88e-6), + fSigmaUpsilon(0.658e-6), + fSigmaUpsilonP(0.218e-6), + fSigmaUpsilonPP(0.122e-6), + fSigmaCCbar(8.9e-3), + fSigmaBBbar(0.33e-3), + + +//x-sections for pp @ 14 TeV: charmonia from hep-ph/0311048 Tab.9, page 19, +// bottomonium from hep-ph/0311048 Tab.9, page 19 taken inton account that +// feed-down from chib is included +/* fCMSEnergy(14), + fSigmaReaction(0.070), + fSigmaJPsi(32.9e-6), + fSigmaChic1(31.8e-6), + fSigmaChic2(52.5e-6), + fSigmaPsiP(7.43e-6), + fSigmaUpsilon(0.989e-6), + fSigmaUpsilonP(0.502e-6), + fSigmaUpsilonPP(0.228e-6), + fSigmaCCbar(11.2e-3), + fSigmaBBbar(0.51e-3), +*/ + fSigmaSilent(kFALSE) { // Constructor } @@ -72,10 +139,38 @@ AliGenMUONCocktailpp::~AliGenMUONCocktailpp() } //_________________________________________________________________________ -void AliGenMUONCocktailpp::Init() +void AliGenMUONCocktailpp::SetResPolarization(Double_t JpsiPol, Double_t PsiPPol, Double_t UpsPol, + Double_t UpsPPol, Double_t UpsPPPol, char *PolFrame){ + + if (strcmp(PolFrame,"kColSop")==0){ + fJpsiPol = (JpsiPol>=-1 && JpsiPol<=1) ? JpsiPol : 0; + fPsiPPol = (PsiPPol>=-1 && PsiPPol<=1) ? PsiPPol : 0; + fUpsPol = (UpsPol>=-1 && UpsPol<=1) ? UpsPol : 0; + fUpsPPol = (UpsPPol>=-1 && UpsPPol<=1) ? UpsPPol : 0; + fUpsPPPol = (UpsPPPol>=-1 && UpsPPPol<=1) ? UpsPPPol : 0; + fPolFrame = 0; + } else if (strcmp(PolFrame,"kHelicity")==0){ + fJpsiPol = (JpsiPol>=-1 && JpsiPol<=1) ? JpsiPol : 0; + fPsiPPol = (PsiPPol>=-1 && PsiPPol<=1) ? PsiPPol : 0; + fUpsPol = (UpsPol>=-1 && UpsPol<=1) ? UpsPol : 0; + fUpsPPol = (UpsPPol>=-1 && UpsPPol<=1) ? UpsPPol : 0; + fUpsPPPol = (UpsPPPol>=-1 && UpsPPPol<=1) ? UpsPPPol : 0; + fPolFrame = 1; + + } else { + AliInfo(Form("The polarization frame is not valid")); + AliInfo(Form("No polarization will be set")); + fJpsiPol=0.; + fPsiPPol=0.; + fUpsPol=0.; + fUpsPPol=0.; + fUpsPPPol=0.; + } +} + +//_________________________________________________________________________ +void AliGenMUONCocktailpp::CreateCocktail() { -// MinBias NN cross section @ pp 14 TeV -PR Vol II p:473 - Double_t sigmaReaction = 0.070; // These limits are only used for renormalization of quarkonia cross section // Pythia events are generated in 4pi @@ -90,6 +185,8 @@ void AliGenMUONCocktailpp::Init() // Ratios with respect to the reaction cross-section in the // kinematics limit of the MUONCocktail Double_t ratiojpsi; + Double_t ratiochic1; + Double_t ratiochic2; Double_t ratiopsiP; Double_t ratioupsilon; Double_t ratioupsilonP; @@ -97,166 +194,412 @@ void AliGenMUONCocktailpp::Init() Double_t ratioccbar; Double_t ratiobbbar; +// Beam energy + Int_t cmsEnergy = fCMSEnergy; + // Cross sections in barns (from PPR Vol. II p: 552) pp - 14 TeV and // corrected from feed down of higher resonances - Double_t sigmajpsi = 49.44e-6; - Double_t sigmapsiP = 7.67e-6; - Double_t sigmaupsilon = 0.989e-6; - Double_t sigmaupsilonP = 0.502e-6; - Double_t sigmaupsilonPP = 0.228e-6; - Double_t sigmaccbar = 11.2e-3; - Double_t sigmabbbar = 0.51e-3; - + Double_t sigmajpsi = fSigmaJPsi; + Double_t sigmachic1 = fSigmaChic1; + Double_t sigmachic2 = fSigmaChic2; + Double_t sigmapsiP = fSigmaPsiP; + Double_t sigmaupsilon = fSigmaUpsilon; + Double_t sigmaupsilonP = fSigmaUpsilonP; + Double_t sigmaupsilonPP = fSigmaUpsilonPP; + Double_t sigmaccbar = fSigmaCCbar; + Double_t sigmabbbar = fSigmaBBbar; + +// Cross sections corrected with the BR in mu+mu- +// (only in case of use of AliDecayerPolarized) + + if(TMath::Abs(fJpsiPol) > 1.e-30) {sigmajpsi = fSigmaJPsi*0.0593;} + if(TMath::Abs(fChic1Pol) > 1.e-30) {sigmachic1 = fSigmaChic1*0.;} // tb consistent + if(TMath::Abs(fChic2Pol) > 1.e-30) {sigmachic2 = fSigmaChic2*0.;} // tb consistent + if(TMath::Abs(fPsiPPol) > 1.e-30) {sigmapsiP = fSigmaPsiP*0.0075;} + if(TMath::Abs(fUpsPol) > 1.e-30) {sigmaupsilon = fSigmaUpsilon*0.0248;} + if(TMath::Abs(fUpsPPol) > 1.e-30) {sigmaupsilonP = fSigmaUpsilonP*0.0193;} + if(TMath::Abs(fUpsPPPol) > 1.e-30) {sigmaupsilonPP = fSigmaUpsilonPP*0.0218;} + +// MinBias NN cross section @ pp 14 TeV -PR Vol II p:473 + Double_t sigmaReaction = fSigmaReaction; + + Int_t eincStart = 10; + AliInfo(Form("the parametrised resonances uses the decay mode %d",fDecayModeResonance)); // Generation using CDF scaled distribution for pp@14TeV (from D. Stocco) //---------------------------------------------------------------------- // Jpsi generator - AliGenParam * genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp", "Jpsi"); + AliGenParam * genjpsi; + if(cmsEnergy == eincStart){ + genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp 10", "Jpsi"); + } else { + genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp ", "Jpsi"); + } // first step: generation in 4pi genjpsi->SetPtRange(0.,100.); genjpsi->SetYRange(-8.,8.); genjpsi->SetPhiRange(0.,360.); genjpsi->SetForceDecay(fDecayModeResonance); - if (!gMC) genjpsi->SetDecayer(fDecayer); + //if (!gMC) genjpsi->SetDecayer(fDecayer); genjpsi->Init(); // generation in 4pi ratiojpsi = sigmajpsi / sigmaReaction * genjpsi->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); // get weight - AliInfo(Form("jpsi prod. cross-section in pp(14 TeV) %5.3g b",sigmajpsi)); - AliInfo(Form("jpsi prod. probability per collision in acceptance %5.3g",ratiojpsi)); + if (!fSigmaSilent) { + AliInfo(Form("jpsi prod. cross-section in pp %5.3g b",sigmajpsi)); + AliInfo(Form("jpsi prod. probability per collision in acceptance %5.3g",ratiojpsi)); + } // second step: generation in selected kinematical range genjpsi->SetPtRange(ptMin, ptMax); genjpsi->SetYRange(yMin, yMax); genjpsi->SetPhiRange(phiMin, phiMax); genjpsi->Init(); // generation in selected kinematic range + + TVirtualMCDecayer *jpsiDec = 0; + if(TMath::Abs(fJpsiPol) > 1.e-30){ + AliInfo(Form("******Setting polarized decayer for J/psi")); + if(fPolFrame==0) { + jpsiDec = new AliDecayerPolarized(fJpsiPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon); + AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fJpsiPol)); + } + if(fPolFrame==1) { + jpsiDec = new AliDecayerPolarized(fJpsiPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon); + AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fJpsiPol)); + } + jpsiDec->SetForceDecay(kAll); + jpsiDec->Init(); + genjpsi->SetDecayer(jpsiDec); + } + AddGenerator(genjpsi, "Jpsi", ratiojpsi); // Adding Generator - fTotalRate+=ratiojpsi; - + + +//---------------------------------------------------------------------- +// Chic1 generator + AliGenParam * genchic1; + if(cmsEnergy == eincStart){ + genchic1 = new AliGenParam(1, AliGenMUONlib::kChic1, "CDF pp 10", "Chic1"); + } else { + genchic1 = new AliGenParam(1, AliGenMUONlib::kChic1, "CDF pp ", "Chic1"); + } +// first step: generation in 4pi + genchic1->SetPtRange(0.,100.); + genchic1->SetYRange(-8.,8.); + genchic1->SetPhiRange(0.,360.); + genchic1->SetForceDecay(fDecayModeResonance); + //if (!gMC) genchic1->SetDecayer(fDecayer); + + genchic1->Init(); // generation in 4pi + ratiochic1 = sigmachic1 / sigmaReaction * genchic1->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); // get weight + if (!fSigmaSilent) { + AliInfo(Form("chic1 prod. cross-section in pp %5.3g b",sigmachic1)); + AliInfo(Form("chic1 prod. probability per collision in acceptance %5.3g",ratiochic1)); + } +// second step: generation in selected kinematical range + genchic1->SetPtRange(ptMin, ptMax); + genchic1->SetYRange(yMin, yMax); + genchic1->SetPhiRange(phiMin, phiMax); + genchic1->Init(); // generation in selected kinematic range + + TVirtualMCDecayer *chic1Dec = 0; + if(fChic1Pol != 0){ + AliInfo(Form("******Setting polarized decayer for Chic1")); + if(fPolFrame==0) { + chic1Dec = new AliDecayerPolarized(fChic1Pol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon); + AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fChic1Pol)); + } + if(fPolFrame==1) { + chic1Dec = new AliDecayerPolarized(fChic1Pol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon); + AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fChic1Pol)); + } + chic1Dec->SetForceDecay(kAll); + chic1Dec->Init(); + genchic1->SetDecayer(chic1Dec); + } + + AddGenerator(genchic1, "Chic1", ratiochic1); // Adding Generator + +//---------------------------------------------------------------------- +// Chic2 generator + AliGenParam * genchic2; + if(cmsEnergy == eincStart){ + genchic2 = new AliGenParam(1, AliGenMUONlib::kChic2, "CDF pp 10", "Chic2"); + } else { + genchic2 = new AliGenParam(1, AliGenMUONlib::kChic2, "CDF pp ", "Chic2"); + } +// first step: generation in 4pi + genchic2->SetPtRange(0.,100.); + genchic2->SetYRange(-8.,8.); + genchic2->SetPhiRange(0.,360.); + genchic2->SetForceDecay(fDecayModeResonance); + //if (!gMC) genchic1->SetDecayer(fDecayer); + + genchic2->Init(); // generation in 4pi + ratiochic2 = sigmachic2 / sigmaReaction * genchic2->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); // get weight + if (!fSigmaSilent) { + AliInfo(Form("chic2 prod. cross-section in pp %5.3g b",sigmachic2)); + AliInfo(Form("chic2 prod. probability per collision in acceptance %5.3g",ratiochic2)); + } +// second step: generation in selected kinematical range + genchic2->SetPtRange(ptMin, ptMax); + genchic2->SetYRange(yMin, yMax); + genchic2->SetPhiRange(phiMin, phiMax); + genchic2->Init(); // generation in selected kinematic range + + TVirtualMCDecayer *chic2Dec = 0; + if(fChic2Pol != 0){ + AliInfo(Form("******Setting polarized decayer for Chic2")); + if(fPolFrame==0) { + chic2Dec = new AliDecayerPolarized(fChic2Pol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon); + AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fChic2Pol)); + } + if(fPolFrame==1) { + chic2Dec = new AliDecayerPolarized(fChic2Pol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon); + AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fChic2Pol)); + } + chic2Dec->SetForceDecay(kAll); + chic2Dec->Init(); + genchic2->SetDecayer(chic2Dec); + } + + AddGenerator(genchic2, "Chic2", ratiochic2); // Adding Generator + //------------------------------------------------------------------ // Psi prime generator - AliGenParam * genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF pp", "PsiP"); + AliGenParam * genpsiP; + if(cmsEnergy == eincStart){ + genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF pp 10", "PsiP"); + } else { + genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF pp", "PsiP"); + } // first step: generation in 4pi genpsiP->SetPtRange(0.,100.); genpsiP->SetYRange(-8.,8.); genpsiP->SetPhiRange(0.,360.); genpsiP->SetForceDecay(fDecayModeResonance); - if (!gMC) genpsiP->SetDecayer(fDecayer); + //if (!gMC) genpsiP->SetDecayer(fDecayer); genpsiP->Init(); // generation in 4pi ratiopsiP = sigmapsiP / sigmaReaction * genpsiP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); - AliInfo(Form("psiP prod. cross-section in pp(14 TeV) %5.3g b",sigmapsiP)); - AliInfo(Form("psiP prod. probability per collision in acceptance %5.3g",ratiopsiP)); + if (!fSigmaSilent) { + AliInfo(Form("psiP prod. cross-section in pp %5.3g b",sigmapsiP)); + AliInfo(Form("psiP prod. probability per collision in acceptance %5.3g",ratiopsiP)); + } // second step: generation in selected kinematical range genpsiP->SetPtRange(ptMin, ptMax); genpsiP->SetYRange(yMin, yMax); genpsiP->SetPhiRange(phiMin, phiMax); genpsiP->Init(); // generation in selected kinematic range + + TVirtualMCDecayer *psipDec = 0; + if(TMath::Abs(fPsiPPol) > 1.e-30){ + AliInfo(Form("******Setting polarized decayer for psi'")); + if(fPolFrame==0) { + psipDec = new AliDecayerPolarized(fPsiPPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon); + AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fPsiPPol)); + } + if(fPolFrame==1) { + psipDec = new AliDecayerPolarized(fPsiPPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon); + AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fPsiPPol)); + } + psipDec->SetForceDecay(kAll); + psipDec->Init(); + genpsiP->SetDecayer(psipDec); + } + AddGenerator(genpsiP, "PsiP", ratiopsiP); // Adding Generator - fTotalRate+=ratiopsiP; - + //------------------------------------------------------------------ // Upsilon 1S generator - AliGenParam * genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "CDF pp", "Upsilon"); + AliGenParam * genupsilon; + if(cmsEnergy == eincStart) { + genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "CDF pp 10", "Upsilon"); + } else { + genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "CDF pp", "Upsilon"); + } // first step: generation in 4pi genupsilon->SetPtRange(0.,100.); genupsilon->SetYRange(-8.,8.); genupsilon->SetPhiRange(0.,360.); genupsilon->SetForceDecay(fDecayModeResonance); - if (!gMC) genupsilon->SetDecayer(fDecayer); + //if (!gMC) genupsilon->SetDecayer(fDecayer); genupsilon->Init(); // generation in 4pi ratioupsilon = sigmaupsilon / sigmaReaction * genupsilon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); - AliInfo(Form("Upsilon 1S prod. cross-section in pp(14 TeV) %5.3g b",sigmaupsilon)); - AliInfo(Form("Upsilon 1S prod. probability per collision in acceptance %5.3g",ratioupsilon)); + if (!fSigmaSilent) { + AliInfo(Form("Upsilon 1S prod. cross-section in pp %5.3g b",sigmaupsilon)); + AliInfo(Form("Upsilon 1S prod. probability per collision in acceptance %5.3g",ratioupsilon)); + } // second step: generation in selected kinematical range genupsilon->SetPtRange(ptMin, ptMax); genupsilon->SetYRange(yMin, yMax); genupsilon->SetPhiRange(phiMin, phiMax); genupsilon->Init(); // generation in selected kinematical range + + TVirtualMCDecayer *upsDec = 0; + if(TMath::Abs(fUpsPol) > 1.e-30){ + AliInfo(Form("******Setting polarized decayer for Upsilon")); + if(fPolFrame==0) { + upsDec = new AliDecayerPolarized(fUpsPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon); + AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fUpsPol)); + } + if(fPolFrame==1) { + upsDec = new AliDecayerPolarized(fUpsPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon); + AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fUpsPol)); + } + upsDec->SetForceDecay(kAll); + upsDec->Init(); + genupsilon->SetDecayer(upsDec); + } + AddGenerator(genupsilon,"Upsilon", ratioupsilon); // Adding Generator - fTotalRate+=ratioupsilon; - //------------------------------------------------------------------ // Upsilon 2S generator - AliGenParam * genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "CDF pp", "UpsilonP"); + AliGenParam * genupsilonP; + if(cmsEnergy == eincStart){ + genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "CDF pp 10", "UpsilonP"); + } else { + genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "CDF pp", "UpsilonP"); + } + // first step: generation in 4pi genupsilonP->SetPtRange(0.,100.); genupsilonP->SetYRange(-8.,8.); genupsilonP->SetPhiRange(0.,360.); genupsilonP->SetForceDecay(fDecayModeResonance); - if (!gMC) genupsilonP->SetDecayer(fDecayer); + //if (!gMC) genupsilonP->SetDecayer(fDecayer); genupsilonP->Init(); // generation in 4pi ratioupsilonP = sigmaupsilonP / sigmaReaction * genupsilonP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); - AliInfo(Form("Upsilon 2S prod. cross-section in pp(14 TeV) %5.3g b",sigmaupsilonP)); - AliInfo(Form("Upsilon 2S prod. probability per collision in acceptance %5.3g",ratioupsilonP)); + if (!fSigmaSilent) { + AliInfo(Form("Upsilon 2S prod. cross-section in pp %5.3g b",sigmaupsilonP)); + AliInfo(Form("Upsilon 2S prod. probability per collision in acceptance %5.3g",ratioupsilonP)); + } // second step: generation in the kinematical range genupsilonP->SetPtRange(ptMin, ptMax); genupsilonP->SetYRange(yMin, yMax); genupsilonP->SetPhiRange(phiMin, phiMax); genupsilonP->Init(); // generation in selected kinematical range + + TVirtualMCDecayer *upspDec = 0; + if(TMath::Abs(fUpsPPol) > 1.e-30){ + AliInfo(Form("******Setting polarized decayer for Upsilon'")); + if(fPolFrame==0) { + upspDec = new AliDecayerPolarized(fUpsPPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon); + AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fUpsPPol)); + } + if(fPolFrame==1) { + upspDec = new AliDecayerPolarized(fUpsPPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon); + AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fUpsPPol)); + } + upspDec->SetForceDecay(kAll); + upspDec->Init(); + genupsilonP->SetDecayer(upspDec); + } + AddGenerator(genupsilonP,"UpsilonP", ratioupsilonP); // Adding Generator - fTotalRate+=ratioupsilonP; //------------------------------------------------------------------ // Upsilon 3S generator - AliGenParam * genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF pp", "UpsilonPP"); + AliGenParam * genupsilonPP; + if(cmsEnergy == eincStart){ + genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF pp 10", "UpsilonPP"); + } + else { + genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF pp", "UpsilonPP"); + } + // first step: generation in 4pi genupsilonPP->SetPtRange(0.,100.); genupsilonPP->SetYRange(-8.,8.); genupsilonPP->SetPhiRange(0.,360.); genupsilonPP->SetForceDecay(fDecayModeResonance); - if (!gMC) genupsilonPP->SetDecayer(fDecayer); + //if (!gMC) genupsilonPP->SetDecayer(fDecayer); genupsilonPP->Init(); // generation in 4pi ratioupsilonPP = sigmaupsilonPP / sigmaReaction * genupsilonPP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); - AliInfo(Form("Upsilon 3S prod. cross-section in pp(14 TeV) %5.3g b",sigmaupsilonPP)); - AliInfo(Form("Upsilon 3S prod. probability per collision in acceptance %5.3g",ratioupsilonPP)); + if (!fSigmaSilent) { + AliInfo(Form("Upsilon 3S prod. cross-section in pp %5.3g b",sigmaupsilonPP)); + AliInfo(Form("Upsilon 3S prod. probability per collision in acceptance %5.3g",ratioupsilonPP)); + } // second step: generation in selected kinematical range genupsilonPP->SetPtRange(ptMin, ptMax); genupsilonPP->SetYRange(yMin, yMax); genupsilonPP->SetPhiRange(phiMin, phiMax); genupsilonPP->Init(); // generation in selected kinematical range + + TVirtualMCDecayer *upsppDec = 0; + if(fUpsPPPol != 0){ + AliInfo(Form("******Setting polarized decayer for Upsilon''")); + if(fPolFrame==0) { + upsppDec = new AliDecayerPolarized(fUpsPPPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon); + AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fUpsPPPol)); + } + if(fPolFrame==1) { + upsppDec = new AliDecayerPolarized(fUpsPPPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon); + AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fUpsPPPol)); + } + upsppDec->SetForceDecay(kAll); + upsppDec->Init(); + genupsilonPP->SetDecayer(upsppDec); + } + AddGenerator(genupsilonPP,"UpsilonPP", ratioupsilonPP); // Adding Generator - fTotalRate+=ratioupsilonPP; //------------------------------------------------------------------ // Generator of charm - AliGenCorrHF *gencharm = new AliGenCorrHF(1, 4); - gencharm->SetMomentumRange(0,9999); - gencharm->SetForceDecay(kAll); - ratioccbar = sigmaccbar/sigmaReaction; - gencharm->Init(); - AddGenerator(gencharm,"CorrHFCharm",ratioccbar); - fTotalRate+=ratioccbar; + AliGenCorrHF *gencharm = new AliGenCorrHF(1, 4, cmsEnergy); + gencharm->SetMomentumRange(0,9999); + gencharm->SetForceDecay(kAll); + ratioccbar = sigmaccbar/sigmaReaction; + //if (!gMC) gencharm->SetDecayer(fDecayer); + gencharm->Init(); + if (!fSigmaSilent) { + AliInfo(Form("c-cbar prod. cross-section in pp %5.3g b",sigmaccbar)); + AliInfo(Form("c-cbar prod. probability per collision in acceptance %5.3g",ratioccbar)); + } + AddGenerator(gencharm,"CorrHFCharm",ratioccbar); //------------------------------------------------------------------ // Generator of beauty - AliGenCorrHF *genbeauty = new AliGenCorrHF(1, 5); - genbeauty->SetMomentumRange(0,9999); - genbeauty->SetForceDecay(kAll); - ratiobbbar = sigmabbbar/sigmaReaction; - genbeauty->Init(); - AddGenerator(genbeauty,"CorrHFBeauty",ratiobbbar); - fTotalRate+=ratiobbbar; + AliGenCorrHF *genbeauty = new AliGenCorrHF(1, 5, cmsEnergy); + genbeauty->SetMomentumRange(0,9999); + genbeauty->SetForceDecay(kAll); + ratiobbbar = sigmabbbar/sigmaReaction; + //if (!gMC) genbeauty->SetDecayer(fDecayer); + genbeauty->Init(); + if (!fSigmaSilent) { + AliInfo(Form("b-bbar prod. cross-section in pp %5.3g b",sigmabbbar)); + AliInfo(Form("b-bbar prod. probability per collision in acceptance %5.3g",ratiobbbar)); + } + AddGenerator(genbeauty,"CorrHFBeauty",ratiobbbar); -//--------------------------t---------------------------------------- +//------------------------------------------------------------------- // Pythia generator - AliGenPythia *pythia = new AliGenPythia(1); - pythia->SetProcess(kPyMbMSEL1); - pythia->SetStrucFunc(kCTEQ5L); - pythia->SetEnergyCMS(14000.); - AliInfo(Form("\n\npythia uses the decay mode %d",fDecayModePythia)); - pythia->SetForceDecay(fDecayModePythia); - pythia->SetPtRange(0.,100.); - pythia->SetYRange(-8.,8.); - pythia->SetPhiRange(0.,360.); - pythia->SetPtHard(2.76,-1.0); - pythia->SwitchHFOff(); - pythia->Init(); - AddGenerator(pythia,"Pythia",1); - fTotalRate+=1.; +// +// This has to go into the Config.C +// +// AliGenPythia *pythia = new AliGenPythia(1); +// pythia->SetProcess(kPyMbMSEL1); +// pythia->SetStrucFunc(kCTEQ5L); +// pythia->SetEnergyCMS(14000.); +// AliInfo(Form("\n\npythia uses the decay mode %d", GetDecayModePythia())); +// Decay_t dt = gener->GetDecayModePythia(); +// pythia->SetForceDecay(dt); +// pythia->SetPtRange(0.,100.); +// pythia->SetYRange(-8.,8.); +// pythia->SetPhiRange(0.,360.); +// pythia->SetPtHard(2.76,-1.0); +// pythia->SwitchHFOff(); +// pythia->Init(); +// AddGenerator(pythia,"Pythia",1); + +} +void AliGenMUONCocktailpp::Init() +{ + // + // Initialisation TIter next(fEntries); AliGenCocktailEntry *entry; if (fStack) { @@ -275,7 +618,10 @@ void AliGenMUONCocktailpp::Generate() AliGenCocktailEntry *preventry = 0; AliGenerator* gen = 0; - TObjArray *partArray = gAlice->GetMCApp()->Particles(); + if (fHeader) delete fHeader; + fHeader = new AliGenCocktailEventHeader("MUON Cocktail Header"); + + const TObjArray *partArray = gAlice->GetMCApp()->Particles(); // Generate the vertex position used by all generators if(fVertexSmear == kPerEvent) Vertex(); @@ -286,10 +632,10 @@ void AliGenMUONCocktailpp::Generate() Bool_t primordialTrigger = kFALSE; while(!primordialTrigger) { //Reseting stack - AliRunLoader * runloader = gAlice->GetRunLoader(); + AliRunLoader * runloader = AliRunLoader::Instance(); if (runloader) if (runloader->Stack()) - runloader->Stack()->Reset(); + runloader->Stack()->Clean(); // Loop over generators and generate events Int_t igen = 0; Int_t npart = 0; @@ -327,15 +673,28 @@ void AliGenMUONCocktailpp::Generate() if((part->Vz() > fMuonOriginCut) && //take only the muons that decayed before the abs + 1 int. length in C abs (part->Theta()*180./TMath::Pi()>fMuonThetaMinCut) && (part->Theta()*180./TMath::Pi()Pt()>fMuonPtCut)) { + (part->Pt()>fMuonPtCut) && + (part->P()>fMuonPCut)) { numberOfMuons++; } } } - if (numberOfMuons >= fMuonMultiplicity) primordialTrigger = kTRUE; + if (numberOfMuons >= fMuonMultiplicity) { + primordialTrigger = kTRUE; + fHeader->SetNProduced(maxPart); + } + } fNSucceded++; + TArrayF eventVertex; + eventVertex.Set(3); + for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j]; + + fHeader->SetPrimaryVertex(eventVertex); + + gAlice->SetGenEventHeader(fHeader); + // AliInfo(Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded)); AliDebug(5,Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded)); }