/************************************************************************** * 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. * **************************************************************************/ /* $Id$ */ // // 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" #include "AliGenMUONCocktailpp.h" #include "AliGenMUONlib.h" #include "AliGenParam.h" #include "AliMC.h" #include "AliRun.h" #include "AliStack.h" #include "AliDecayer.h" #include "AliLog.h" #include "AliGenCorrHF.h" #include "AliDecayerPolarized.h" ClassImp(AliGenMUONCocktailpp) //________________________________________________________________________ AliGenMUONCocktailpp::AliGenMUONCocktailpp() :AliGenCocktail(), fDecayer(0), fDecayModeResonance(kAll), fDecayModePythia(kAll), fMuonMultiplicity(0), fMuonPtCut(0.), fMuonPCut(0.), fMuonThetaMinCut(0.), fMuonThetaMaxCut(180.), fMuonOriginCut(-999.), fNSucceded(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::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() { // 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; Double_t ptMax = fPtMax; Double_t yMin = fYMin;; Double_t yMax = fYMax;; Double_t phiMin = fPhiMin*180./TMath::Pi(); 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)); // 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)); // 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; iSetMomentumRange(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); //------------------------------------------------------------------ // 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 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); } } AddGenerator(genReso,nameReso,ratioReso); // Adding Generator } //------------------------------------------------------------------- void AliGenMUONCocktailpp::Init() { // Initialisation TIter next(fEntries); AliGenCocktailEntry *entry; if (fStack) { while((entry = (AliGenCocktailEntry*)next())) { entry->Generator()->SetStack(fStack); } } } //_________________________________________________________________________ void AliGenMUONCocktailpp::Generate() { // Generate event TIter next(fEntries); AliGenCocktailEntry *entry = 0; AliGenCocktailEntry *preventry = 0; AliGenerator* gen = 0; 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(); // Loop on primordialTrigger: // minimum muon multiplicity above a pt cut in a theta acceptance region Bool_t primordialTrigger = kFALSE; while(!primordialTrigger) { //Reseting stack (if fMuonMultiplicity > 0 : S. Grigoryan, June 2012) AliRunLoader * runloader = AliRunLoader::Instance(); // if (runloader) if (runloader && fMuonMultiplicity > 0) if (runloader->Stack()) runloader->Stack()->Clean(); // Loop over generators and generate events Int_t igen = 0; Int_t npart = 0; const char* genName = 0; while((entry = (AliGenCocktailEntry*)next())) { 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()); if (npart > 0) { igen++; if (igen == 1) entry->SetFirst(0); else entry->SetFirst((partArray->GetEntriesFast())+1); 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(); for(iPart=0; iPartGetMCApp()->Particle(iPart); if ( TMath::Abs(part->GetPdgCode()) == 13 ){ 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->P()>fMuonPCut)) { numberOfMuons++; } } } 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)); }