X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=EVGEN%2FAliGenMUONCocktailpp.cxx;h=6cd7ce84cd7239751a3efbbc491194d60a47afb8;hp=e96d1e4a69c2544d4b6f7211218008a015b36d4e;hb=787342d6495a9e0bc314ed69c310cf7a4961f4a0;hpb=aaa95f22c7d5b29a10d6a410a4797c40e79ec380 diff --git a/EVGEN/AliGenMUONCocktailpp.cxx b/EVGEN/AliGenMUONCocktailpp.cxx index e96d1e4a69c..6cd7ce84cd7 100644 --- a/EVGEN/AliGenMUONCocktailpp.cxx +++ b/EVGEN/AliGenMUONCocktailpp.cxx @@ -16,19 +16,79 @@ /* $Id$ */ // -// Classe to create the coktail for physics with muons for pp at 14 TeV -// The followoing sources: -// jpsi,psiP, upsilon, upsilonP, upsilonPP are added to Pythia +// Class to create the coktail for physics with muons for pp collisions +// using the The followoing sources: +// jpsi, psiP, upsilon, upsilonP, upsilonPP, open charm and open beauty // The free parameeters are : // pp reaction cross-section // 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 +//------------------------ +// 04/2010: +// - CMS energy passed via parameter +// i.e. gener->SetCMSEnergy(AliGenMUONCocktailpp::kCMS07TeV) in Config.C +// - resonances now added to the cocktail via AddReso2Generator +// - cleaning +// B.Vulpescu & P.Crochet +//----------------------- +// 10/2011: +// - added the cocktail for p-Pb & Pb-p @ 8.8 TeV with 4 centrality bins and +// for Pb-Pb @ 2.76 TeV with 11 centrality bins. Bins should be defined also +// in the Config.C with one AliGenMUONCocktailpp per bin. These generators +// included in a AliGenCocktail together with an event generator (e.g. Hijing) +// providing the underlying event and collision centrality. The bin number n +// passed via AliGenMUONCocktailpp::SetCentralityBin(n). +// See details in my presentation at the PWG3-Muon meeting (05.10.2011): +// https://indico.cern.ch/conferenceDisplay.py?confId=157367 +// - simplifications and bug fix in CreateCocktail() +// S. Grigoryan +//----------------------- +// 06/2012: +// - added the cocktail for p-Pb & Pb-p @ 2.76, 4.4 & 5.03 TeV with 4 centrality +// bins, using the EPS09-LO shadowing computed for 5.03 TeV. Energies are set by +// AliGenMUONCocktailpp::SetCMSEnergy(int CMSEnergyCode), CMSEnergy codes are +// defined in AliGenMUONCocktailpp.h. +// - added functions to scale x-section of JPsi, Charmonia, Bottomonia, CCbar & BBbar +// in Config.C to manage the statistics. Example of usage (in a cocktail with Hijing): +/* + AliGenCocktail *cocktail = new AliGenCocktail(); + cocktail->AddGenerator(hijing,"hijing",1); + TFormula *form[nb]; // nb - number of centrality bins with impact params bBins[i] + AliGenMUONCocktailpp *gen[nb]; + for (Int_t i=0; iSetParameters(i+1, bBins[i], bBins[i+1]); + gen[i] = MuonCocktail(); + gen[i]->SetCentralityBin(i+1); + gen[i]->ScaleJPsi(100.); + gen[i]->CreateCocktail(); + cocktail->AddGenerator(gen[i], Form("g%d",i), 101+i, form[i]); + } +AliGenMUONCocktailpp* MuonCocktail() { + AliGenMUONCocktailpp *mc = new AliGenMUONCocktailpp(); + .................................... + return mc; +} +*/ +// - a bug fixed in the function Generate() +// S. Grigoryan + #include #include #include #include - #include "AliGenCocktailEventHeader.h" #include "AliGenCocktailEntry.h" @@ -40,8 +100,8 @@ #include "AliStack.h" #include "AliDecayer.h" #include "AliLog.h" -#include "AliGenPythia.h" - +#include "AliGenCorrHF.h" +#include "AliDecayerPolarized.h" ClassImp(AliGenMUONCocktailpp) @@ -51,31 +111,294 @@ 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), + fCentralityBin(0), + fScaleJPsi(1), + fScaleCharmonia(1), + fScaleBottomonia(1), + fScaleCCbar(1), + fScaleBBbar(1), + + fJpsiPol(0), + fChic1Pol(0), + fChic2Pol(0), + fPsiPPol(0), + fUpsPol(0), + fUpsPPol(0), + fUpsPPPol(0), + fPolFrame(0), + + fCMSEnergyTeV(0), + fCMSEnergyTeVArray(), + fSigmaReaction(0), + fSigmaReactionArray(), + fSigmaJPsi(0), + fSigmaJPsiArray(), + fSigmaChic1(0), + fSigmaChic1Array(), + fSigmaChic2(0), + fSigmaChic2Array(), + fSigmaPsiP(0), + fSigmaPsiPArray(), + fSigmaUpsilon(0), + fSigmaUpsilonArray(), + fSigmaUpsilonP(0), + fSigmaUpsilonPArray(), + fSigmaUpsilonPP(0), + fSigmaUpsilonPPArray(), + fSigmaCCbar(0), + fSigmaCCbarArray(), + fSigmaBBbar(0), + fSigmaBBbarArray(), + + fSigmaSilent(kFALSE) { // Constructor + +// x-sections for pp @ 7 TeV: +// -charmonia: 4pi integral of fit function for inclusive J/psi dsigma/dy LHC data +// gives 60 mub; so sigma_prompt = 54 mub, while Ref = R.Vogt_arXiv:1003.3497 (Table 2) +// gives 35 mub. Below we use sigma_direct from the Ref scaled by the factor 54/35. +// -bottomonia: 4pi integral of fit function for inclusive Upsilon1S dsigma/dy LHC data +// gives 0.56 mub, sigmas for 2S & 3S obtained using LHCb data for ratios 2S/1S & 3S/1S +// -ccbar & bbbar: NLO pQCD computations - http://www-alice.gsi.de/ana/MNR/results.html + fCMSEnergyTeVArray[0] = 7.00; + fSigmaReactionArray[0] = 0.070; + fSigmaJPsiArray[0] = 33.6e-6; + fSigmaChic1Array[0] = 32.6e-6; + fSigmaChic2Array[0] = 53.8e-6; + fSigmaPsiPArray[0] = 7.6e-6; + fSigmaUpsilonArray[0] = 0.56e-6; + fSigmaUpsilonPArray[0] = 0.18e-6; + fSigmaUpsilonPPArray[0] = 0.08e-6; + fSigmaCCbarArray[0] = 6.91e-3; + fSigmaBBbarArray[0] = 0.232e-3; + +//x-sections for pp @ 10 TeV: charmonia and bottomonia from 14 TeV numbers +// scaled down according to ccbar and bbbar cross-sections + fCMSEnergyTeVArray[1] = 10.00; + fSigmaReactionArray[1] = 0.070; + fSigmaJPsiArray[1] = 26.06e-6; + fSigmaChic1Array[1] = 25.18e-6; + fSigmaChic2Array[1] = 41.58e-6; + fSigmaPsiPArray[1] = 5.88e-6; + fSigmaUpsilonArray[1] = 0.658e-6; + fSigmaUpsilonPArray[1] = 0.218e-6; + fSigmaUpsilonPPArray[1] = 0.122e-6; + fSigmaCCbarArray[1] = 8.9e-3; + fSigmaBBbarArray[1] = 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 into account that +// feed-down from chib is included + fCMSEnergyTeVArray[2] = 14.00; + fSigmaReactionArray[2] = 0.070; + fSigmaJPsiArray[2] = 32.9e-6; + fSigmaChic1Array[2] = 31.8e-6; + fSigmaChic2Array[2] = 52.5e-6; + fSigmaPsiPArray[2] = 7.43e-6; + fSigmaUpsilonArray[2] = 0.989e-6; + fSigmaUpsilonPArray[2] = 0.502e-6; + fSigmaUpsilonPPArray[2] = 0.228e-6; + fSigmaCCbarArray[2] = 11.2e-3; + fSigmaBBbarArray[2] = 0.445e-3; + +// x-sections for Min. Bias p-Pb & Pb-p @ 2.76 TeV: charmonia and bottomonia +// from 7 TeV numbers scaled according to pQCD ccbar and bbbar x-sections +// and with Glauber scaling + fCMSEnergyTeVArray[3] = 2.00; // for 2.76 TeV + fSigmaReactionArray[3] = 2.10; + fSigmaJPsiArray[3] = 3.54e-3; // 208*0.507*33.6e-6 + fSigmaChic1Array[3] = 3.44e-3; + fSigmaChic2Array[3] = 5.66e-3; + fSigmaPsiPArray[3] = 0.80e-3; + fSigmaUpsilonArray[3] = 0.045e-3; // 208*0.384*0.56e-6 + fSigmaUpsilonPArray[3] = 0.014e-3; + fSigmaUpsilonPPArray[3] = 0.006e-3; + fSigmaCCbarArray[3] = 0.73; // 208*3.50e-3 + fSigmaBBbarArray[3] = 0.019; // 208*0.089e-3 + + fCMSEnergyTeVArray[4] = -fCMSEnergyTeVArray[3]; + fSigmaReactionArray[4] = fSigmaReactionArray[3]; + fSigmaJPsiArray[4] = fSigmaJPsiArray[3]; + fSigmaChic1Array[4] = fSigmaChic1Array[3]; + fSigmaChic2Array[4] = fSigmaChic2Array[3]; + fSigmaPsiPArray[4] = fSigmaPsiPArray[3]; + fSigmaUpsilonArray[4] = fSigmaUpsilonArray[3]; + fSigmaUpsilonPArray[4] = fSigmaUpsilonPArray[3]; + fSigmaUpsilonPPArray[4] = fSigmaUpsilonPPArray[3]; + fSigmaCCbarArray[4] = fSigmaCCbarArray[3]; + fSigmaBBbarArray[4] = fSigmaBBbarArray[3]; + +// x-sections for Min. Bias p-Pb & Pb-p @ 4.4 TeV: charmonia and bottomonia +// from 7 TeV numbers scaled according to pQCD ccbar and bbbar x-sections +// and with Glauber scaling + fCMSEnergyTeVArray[5] = 4.00; // for 4.4 TeV + fSigmaReactionArray[5] = 2.10; + fSigmaJPsiArray[5] = 5.00e-3; // 208*0.715*33.6e-6 + fSigmaChic1Array[5] = 4.86e-3; + fSigmaChic2Array[5] = 7.99e-3; + fSigmaPsiPArray[5] = 1.12e-3; + fSigmaUpsilonArray[5] = 0.074e-3; // 208*0.629*0.56e-6 + fSigmaUpsilonPArray[5] = 0.023e-3; + fSigmaUpsilonPPArray[5] = 0.010e-3; + fSigmaCCbarArray[5] = 1.03; // 208*4.94e-3 + fSigmaBBbarArray[5] = 0.030; // 208*0.146e-3 + + fCMSEnergyTeVArray[6] = -fCMSEnergyTeVArray[5]; + fSigmaReactionArray[6] = fSigmaReactionArray[5]; + fSigmaJPsiArray[6] = fSigmaJPsiArray[5]; + fSigmaChic1Array[6] = fSigmaChic1Array[5]; + fSigmaChic2Array[6] = fSigmaChic2Array[5]; + fSigmaPsiPArray[6] = fSigmaPsiPArray[5]; + fSigmaUpsilonArray[6] = fSigmaUpsilonArray[5]; + fSigmaUpsilonPArray[6] = fSigmaUpsilonPArray[5]; + fSigmaUpsilonPPArray[6] = fSigmaUpsilonPPArray[5]; + fSigmaCCbarArray[6] = fSigmaCCbarArray[5]; + fSigmaBBbarArray[6] = fSigmaBBbarArray[5]; + +// x-sections for Min. Bias p-Pb & Pb-p @ 5.03 TeV: charmonia and bottomonia +// from 7 TeV numbers scaled according to pQCD ccbar and bbbar x-sections +// and with Glauber scaling + fCMSEnergyTeVArray[7] = 5.00; // for 5.03 TeV + fSigmaReactionArray[7] = 2.10; + fSigmaJPsiArray[7] = 5.50e-3; // 208*0.787*33.6e-6 + fSigmaChic1Array[7] = 5.35e-3; + fSigmaChic2Array[7] = 8.79e-3; + fSigmaPsiPArray[7] = 1.23e-3; + fSigmaUpsilonArray[7] = 0.083e-3; // 208*0.716*0.56e-6 + fSigmaUpsilonPArray[7] = 0.026e-3; + fSigmaUpsilonPPArray[7] = 0.011e-3; + fSigmaCCbarArray[7] = 1.13; // 208*5.44e-3 + fSigmaBBbarArray[7] = 0.035; // 208*0.166e-3 + + fCMSEnergyTeVArray[8] = -fCMSEnergyTeVArray[7]; + fSigmaReactionArray[8] = fSigmaReactionArray[7]; + fSigmaJPsiArray[8] = fSigmaJPsiArray[7]; + fSigmaChic1Array[8] = fSigmaChic1Array[7]; + fSigmaChic2Array[8] = fSigmaChic2Array[7]; + fSigmaPsiPArray[8] = fSigmaPsiPArray[7]; + fSigmaUpsilonArray[8] = fSigmaUpsilonArray[7]; + fSigmaUpsilonPArray[8] = fSigmaUpsilonPArray[7]; + fSigmaUpsilonPPArray[8] = fSigmaUpsilonPPArray[7]; + fSigmaCCbarArray[8] = fSigmaCCbarArray[7]; + fSigmaBBbarArray[8] = fSigmaBBbarArray[7]; + +// x-sections for Min. Bias p-Pb & Pb-p @ 8.8 TeV: charmonia and bottomonia +// from 7 TeV numbers scaled according to pQCD ccbar and bbbar x-sections +// and with Glauber scaling + fCMSEnergyTeVArray[9] = 9.00; // for 8.8 TeV + fSigmaReactionArray[9] = 2.10; + fSigmaJPsiArray[9] = 8.19e-3; // 208*1.172*33.6e-6 + fSigmaChic1Array[9] = 7.95e-3; + fSigmaChic2Array[9] = 13.1e-3; + fSigmaPsiPArray[9] = 1.85e-3; + fSigmaUpsilonArray[9] = 0.146e-3; // 208*1.25*0.56e-6 + fSigmaUpsilonPArray[9] = 0.047e-3; + fSigmaUpsilonPPArray[9] = 0.021e-3; + fSigmaCCbarArray[9] = 1.68; // 208*8.1e-3 + fSigmaBBbarArray[9] = 0.061; // 208*0.29e-3 + + fCMSEnergyTeVArray[10] = -fCMSEnergyTeVArray[9]; + fSigmaReactionArray[10] = fSigmaReactionArray[9]; + fSigmaJPsiArray[10] = fSigmaJPsiArray[9]; + fSigmaChic1Array[10] = fSigmaChic1Array[9]; + fSigmaChic2Array[10] = fSigmaChic2Array[9]; + fSigmaPsiPArray[10] = fSigmaPsiPArray[9]; + fSigmaUpsilonArray[10] = fSigmaUpsilonArray[9]; + fSigmaUpsilonPArray[10] = fSigmaUpsilonPArray[9]; + fSigmaUpsilonPPArray[10] = fSigmaUpsilonPPArray[9]; + fSigmaCCbarArray[10] = fSigmaCCbarArray[9]; + fSigmaBBbarArray[10] = fSigmaBBbarArray[9]; + +// x-sections for Min. Bias Pb-Pb @ 2.76 TeV: charmonia and bottomonia +// from 7 TeV numbers scaled according to pQCD ccbar and bbbar x-sections +// and with Glauber scaling + fCMSEnergyTeVArray[11] = 3.00; // for 2.76 TeV + fSigmaReactionArray[11] = 7.65; + fSigmaJPsiArray[11] = 0.737; // 208*208*0.507*33.6e-6 + fSigmaChic1Array[11] = 0.715; + fSigmaChic2Array[11] = 1.179; + fSigmaPsiPArray[11] = 0.166; + fSigmaUpsilonArray[11] = 0.0093; // 208*208*0.384*0.56e-6 + fSigmaUpsilonPArray[11] = 0.0030; + fSigmaUpsilonPPArray[11] = 0.0013; + fSigmaCCbarArray[11] = 151.; // 208*208*3.50e-3 + fSigmaBBbarArray[11] = 3.8; // 208*208*0.089e-3 + } //_________________________________________________________________________ AliGenMUONCocktailpp::~AliGenMUONCocktailpp() +{ // Destructor + +} + +//_________________________________________________________________________ +void AliGenMUONCocktailpp::SetCMSEnergy(CMSEnergyCode cmsEnergy) { - +// setter for CMSEnergy and corresponding cross-sections + fCMSEnergyTeV = fCMSEnergyTeVArray[cmsEnergy]; + fSigmaReaction = fSigmaReactionArray[cmsEnergy]; + fSigmaJPsi = fSigmaJPsiArray[cmsEnergy]; + fSigmaChic1 = fSigmaChic1Array[cmsEnergy]; + fSigmaChic2 = fSigmaChic2Array[cmsEnergy]; + fSigmaPsiP = fSigmaPsiPArray[cmsEnergy]; + fSigmaUpsilon = fSigmaUpsilonArray[cmsEnergy]; + fSigmaUpsilonP = fSigmaUpsilonPArray[cmsEnergy]; + fSigmaUpsilonPP = fSigmaUpsilonPPArray[cmsEnergy]; + fSigmaCCbar = fSigmaCCbarArray[cmsEnergy]; + fSigmaBBbar = fSigmaBBbarArray[cmsEnergy]; } //_________________________________________________________________________ -void AliGenMUONCocktailpp::Init() +void AliGenMUONCocktailpp::SetResPolarization(Double_t JpsiPol, Double_t PsiPPol, Double_t UpsPol, + Double_t UpsPPol, Double_t UpsPPPol, char *PolFrame){ +// setter for resonances polarization + 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; - +// create and add resonances and open HF to the coctail + Int_t cmsEnergy = Int_t(fCMSEnergyTeV); + + // For temporary use of p-Pb & Pb-p shadowing at 5.03 TeV for lower energies + if (cmsEnergy == 2 || cmsEnergy == 4) cmsEnergy = 5; + if (cmsEnergy == -2 || cmsEnergy == -4) cmsEnergy = -5; + // These limits are only used for renormalization of quarkonia cross section // Pythia events are generated in 4pi Double_t ptMin = fPtMin; @@ -86,148 +409,257 @@ void AliGenMUONCocktailpp::Init() Double_t phiMax = fPhiMax*180./TMath::Pi(); AliInfo(Form("Ranges pT:%4.1f : %4.1f GeV/c, y:%4.2f : %4.2f, Phi:%5.1f : %5.1f degres",ptMin,ptMax,yMin,yMax,phiMin,phiMax)); -// Ratios with respect to the reaction cross-section in the -// kinematics limit of the MUONCocktail - Double_t ratiojpsi; - Double_t ratiopsiP; - Double_t ratioupsilon; - Double_t ratioupsilonP; - Double_t ratioupsilonPP; - -// 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; - +// Cross sections in barns + + 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;} + +// Cross sections scaled to manage the statistics + + sigmajpsi *= fScaleJPsi*fScaleCharmonia; + sigmachic1 *= fScaleCharmonia; + sigmachic2 *= fScaleCharmonia; + sigmapsiP *= fScaleCharmonia; + sigmaupsilon *= fScaleBottomonia; + sigmaupsilonP *= fScaleBottomonia; + sigmaupsilonPP *= fScaleBottomonia; + sigmaccbar *= fScaleCCbar; + sigmabbbar *= fScaleBBbar; + 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"); -// 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); - - 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)); -// 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 - AddGenerator(genjpsi, "Jpsi", ratiojpsi); // Adding Generator - fTotalRate+=ratiojpsi; +// Create and add resonances to the generator + AliGenParam * genjpsi=0; + AliGenParam * genchic1=0; + AliGenParam * genchic2=0; + AliGenParam * genpsiP=0; + AliGenParam * genupsilon=0; + AliGenParam * genupsilonP=0; + AliGenParam * genupsilonPP=0; + Char_t nameJpsi[10]; + Char_t nameChic1[10]; + Char_t nameChic2[10]; + Char_t namePsiP[10]; + Char_t nameUps[10]; + Char_t nameUpsP[10]; + Char_t nameUpsPP[10]; + + snprintf(nameJpsi,10, "Jpsi"); + snprintf(nameChic1,10, "Chic1"); + snprintf(nameChic2,10, "Chic2"); + snprintf(namePsiP,10, "PsiP"); + snprintf(nameUps,10, "Ups"); + snprintf(nameUpsP,10, "UpsP"); + snprintf(nameUpsPP,10, "UpsPP"); + + Char_t tname[40] = ""; + if(cmsEnergy == 10) {snprintf(tname, 40, "CDF pp 10"); + } else if (cmsEnergy == 14){snprintf(tname, 40, "CDF pp"); + } else if (cmsEnergy == 7) {snprintf(tname, 40, "pp 7"); + // } else if (cmsEnergy == 2) {snprintf(tname, 40, "pp 2.76"); + } else if (cmsEnergy == 5) {snprintf(tname, 40, "pPb 5.03"); + if (fCentralityBin > 0) snprintf(tname, 40, "pPb 5.03c%d",fCentralityBin); + } else if (cmsEnergy == -5){snprintf(tname, 40, "Pbp 5.03"); + if (fCentralityBin > 0) snprintf(tname, 40, "Pbp 5.03c%d",fCentralityBin); + } else if (cmsEnergy == 9) {snprintf(tname, 40, "pPb 8.8"); + if (fCentralityBin > 0) snprintf(tname, 40, "pPb 8.8c%d",fCentralityBin); + } else if (cmsEnergy == -9){snprintf(tname, 40, "Pbp 8.8"); + if (fCentralityBin > 0) snprintf(tname, 40, "Pbp 8.8c%d",fCentralityBin); + } else if (cmsEnergy == 3) {snprintf(tname, 40, "PbPb 2.76"); + if (fCentralityBin > 0) snprintf(tname, 40, "PbPb 2.76c%d",fCentralityBin); + } else { + AliError("Initialisation failed, wrong cmsEnergy"); + return; + } + genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, tname, "Jpsi"); + genchic1 = new AliGenParam(1, AliGenMUONlib::kChic1, tname, "Chic1"); + genchic2 = new AliGenParam(1, AliGenMUONlib::kChic2, tname, "Chic2"); + genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, tname, "PsiP"); + genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, tname, "Upsilon"); + genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, tname, "UpsilonP"); + genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, tname, "UpsilonPP"); + +// Hard process yield per pA or AA collision for i-th centrality bin is R*r[i]*shad[i] +// where R is the ratio of hard and geometrical x-sections, r[i] is the ratio of these +// x-section fractions for given centrality and shad[i] is the shadowing factor (in 4pi). +// The latter is assumed to be the same for HF-hadrons & quarkonia of the same flavour. + Int_t i = 0; + Double_t chard[20] = {0}; // charm & beauty shadowing factors are different + Double_t bhard[20] = {0}; + chard[0] = 1; // 1st element for pp and min. bias (MB) collisions + bhard[0] = 1; + +// 4 centrality bins for p-Pb & Pb-p at 5.03 TeV: 0-20-40-60-100 % + if (cmsEnergy == 5 || cmsEnergy == -5) { + const Int_t n5 = 5; // 1st element for MB collisions + Double_t r5[n5] = {1, 1.936, 1.473, 0.914, 0.333}; // ratio of hard-over-geo fractions + Double_t cshad5[n5] = {0.806, 0.742, 0.796, 0.870, 0.955};// EPS09-LO shadowing factors + Double_t bshad5[n5] = {0.917, 0.889, 0.913, 0.944, 0.981}; + for(i=0; iSetPtRange(0.,100.); - genpsiP->SetYRange(-8.,8.); - genpsiP->SetPhiRange(0.,360.); - genpsiP->SetForceDecay(fDecayModeResonance); - 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)); -// 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 - AddGenerator(genpsiP, "PsiP", ratiopsiP); // Adding Generator - fTotalRate+=ratiopsiP; - +// Generator of charm + AliGenCorrHF *gencharm = new AliGenCorrHF(1, 4, cmsEnergy); + gencharm->SetMomentumRange(0,9999); + gencharm->SetForceDecay(kAll); + Double_t ratioccbar = chard[fCentralityBin]*sigmaccbar/fSigmaReaction; + if (!TVirtualMC::GetMC()) 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); //------------------------------------------------------------------ -// Upsilon 1S generator - AliGenParam * 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); - 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)); +// Generator of beauty + AliGenCorrHF *genbeauty = new AliGenCorrHF(1, 5, cmsEnergy); + genbeauty->SetMomentumRange(0,9999); + genbeauty->SetForceDecay(kAll); + Double_t ratiobbbar = bhard[fCentralityBin]*sigmabbbar/fSigmaReaction; + if (!TVirtualMC::GetMC()) 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); + +//------------------------------------------------------------------- +// Pythia generator +// +// 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::AddReso2Generator(Char_t* nameReso, + AliGenParam* const genReso, + Double_t sigmaReso, + Double_t polReso) +{ +// add resonances to the cocktail + Double_t phiMin = fPhiMin*180./TMath::Pi(); + Double_t phiMax = fPhiMax*180./TMath::Pi(); + + // first step: generation in 4pi + genReso->SetPtRange(0.,100.); + genReso->SetYRange(-8.,8.); + genReso->SetPhiRange(0.,360.); + genReso->SetForceDecay(fDecayModeResonance); + if (!TVirtualMC::GetMC()) genReso->SetDecayer(fDecayer); + genReso->Init(); // generation in 4pi +// Ratios with respect to the reaction cross-section in the +// kinematics limit of the MUONCocktail + Double_t ratioReso = sigmaReso / fSigmaReaction * genReso->GetRelativeArea(fPtMin,fPtMax,fYMin,fYMax,phiMin,phiMax); + if (!fSigmaSilent) { + AliInfo(Form("%s prod. cross-section in pp %5.3g b",nameReso,sigmaReso)); + AliInfo(Form("%s prod. probability per collision in acceptance %5.3g",nameReso,ratioReso)); + } // 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 - AddGenerator(genupsilon,"Upsilon", ratioupsilon); // Adding Generator - fTotalRate+=ratioupsilon; - -//------------------------------------------------------------------ -// Upsilon 2S generator - AliGenParam * 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); - 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)); -// 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 - AddGenerator(genupsilonP,"UpsilonP", ratioupsilonP); // Adding Generator - fTotalRate+=ratioupsilonP; + genReso->SetPtRange(fPtMin, fPtMax); + genReso->SetYRange(fYMin, fYMax); + genReso->SetPhiRange(phiMin, phiMax); + genReso->Init(); // generation in selected kinematical range + + TVirtualMCDecayer *decReso = 0; + if(TMath::Abs(polReso) > 1.e-30){ + AliInfo(Form("******Setting polarized decayer for %s''",nameReso)); + if(fPolFrame==0) { + decReso = new AliDecayerPolarized(polReso,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon); + AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",polReso)); + } + if(fPolFrame==1) { + decReso = new AliDecayerPolarized(polReso,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon); + AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",polReso)); + } + if (decReso) { + decReso->SetForceDecay(kAll); + decReso->Init(); + genReso->SetDecayer(decReso); + } + } -//------------------------------------------------------------------ -// Upsilon 3S generator - AliGenParam * 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); - 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)); -// 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 - AddGenerator(genupsilonPP,"UpsilonPP", ratioupsilonPP); // Adding Generator - fTotalRate+=ratioupsilonPP; -//------------------------------------------------------------------ -// 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->Init(); - AddGenerator(pythia,"Pythia",1); - fTotalRate+=1.; + AddGenerator(genReso,nameReso,ratioReso); // Adding Generator +} +//------------------------------------------------------------------- +void AliGenMUONCocktailpp::Init() +{ + // Initialisation TIter next(fEntries); AliGenCocktailEntry *entry; if (fStack) { @@ -243,10 +675,12 @@ void AliGenMUONCocktailpp::Generate() // Generate event TIter next(fEntries); AliGenCocktailEntry *entry = 0; - 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(); @@ -256,11 +690,12 @@ void AliGenMUONCocktailpp::Generate() Bool_t primordialTrigger = kFALSE; while(!primordialTrigger) { - //Reseting stack - AliRunLoader * runloader = gAlice->GetRunLoader(); - if (runloader) + //Reseting stack (if fMuonMultiplicity > 0 : S. Grigoryan, June 2012) + AliRunLoader * runloader = AliRunLoader::Instance(); + // if (runloader) + if (runloader && fMuonMultiplicity > 0) if (runloader->Stack()) - runloader->Stack()->Reset(); + runloader->Stack()->Clean(); // Loop over generators and generate events Int_t igen = 0; Int_t npart = 0; @@ -269,6 +704,7 @@ void AliGenMUONCocktailpp::Generate() gen = entry->Generator(); genName = entry->GetName(); gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2)); + gen->SetTime(fTime); npart = (strcmp(genName,"Pythia") == 0) ? 1 : gRandom->Poisson(entry->Rate()); @@ -281,18 +717,16 @@ void AliGenMUONCocktailpp::Generate() gen->SetNumberParticles(npart); gen->Generate(); entry->SetLast(partArray->GetEntriesFast()); - preventry = entry; } } next.Reset(); + // Testing primordial trigger: Single muons or dimuons with Pt above a Pt cut // in the muon spectrometer acceptance Int_t iPart; fNGenerated++; - Int_t numberOfMuons=0; - - Int_t maxPart = partArray->GetEntriesFast(); + Int_t numberOfMuons=0;Int_t maxPart = partArray->GetEntriesFast(); for(iPart=0; iPartGetMCApp()->Particle(iPart); @@ -300,15 +734,29 @@ 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); + fHeader->SetInteractionTime(fTime); + + 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)); }