X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=EVGEN%2FAliGenMUONCocktail.cxx;h=89fe4de562b6fd9a6fb432be1040a52418734fc7;hb=9902ce1e57248fe067265ced3ac5171115d0269e;hp=028a25c54f629c3c24c7403559d99efa0ee5c3ab;hpb=84954c474d5959c7aa0260ac1f7eb905bb5cc086;p=u%2Fmrichter%2FAliRoot.git diff --git a/EVGEN/AliGenMUONCocktail.cxx b/EVGEN/AliGenMUONCocktail.cxx index 028a25c54f6..89fe4de562b 100644 --- a/EVGEN/AliGenMUONCocktail.cxx +++ b/EVGEN/AliGenMUONCocktail.cxx @@ -15,209 +15,345 @@ /* $Id$ */ +// // Classe to create the MUON coktail for physics in the Alice muon spectrometer -// Gines Martinez, jan 2004, Nantes martinez@in2p3.fr - - +// The followoing muons sources are included in this cocktail: +// jpsi, upsilon, non-correlated open and beauty. +// The free parameeters are : +// pp reaction cross-section +// production cross-sections in pp collisions and +// branching ratios in the muon channel and shadowing +// Hard probes are supposed to scale with Ncoll and hadronic muon production with (0.8Ncoll+0.2*Npart) +// There is a primordial trigger which requires : +// a minimum muon multiplicity above a pT cut in a theta acceptance cone // -#include #include -#include #include +#include -#include "AliGenParam.h" -#include "AliGenMUONlib.h" -#include "AliGenMUONCocktail.h" #include "AliGenCocktailEntry.h" -#include "AliCollisionGeometry.h" -#include "AliRun.h" +#include "AliGenMUONCocktail.h" +#include "AliGenMUONlib.h" +#include "AliFastGlauber.h" +#include "AliGenParam.h" #include "AliMC.h" +#include "AliRun.h" #include "AliStack.h" +#include "AliLog.h" ClassImp(AliGenMUONCocktail) //________________________________________________________________________ AliGenMUONCocktail::AliGenMUONCocktail() - :AliGenCocktail() + :AliGenCocktail(), + fFastGlauber (0x0), + fTotalRate(0), + fMuonMultiplicity(1), + fMuonPtCut(1.), + fMuonThetaMinCut(171.), + fMuonThetaMaxCut(178.), + fNSucceded(0), + fNGenerated(0), + fLowImpactParameter(0.), + fHighImpactParameter(5.), + fAverageImpactParameter(0.), + fNumberOfCollisions(0.), + fNumberOfParticipants(0.), + fHadronicMuons(kTRUE), + fInvMassCut (kFALSE), + fInvMassMinCut (0.), + fInvMassMaxCut (100.) { // Constructor - fTotalRate =0; - fNSucceded=0; - fNGenerated=0; - fMuonMultiplicity=1; - fMuonPtCut= 1.; - fMuonThetaMinCut=171.; - fMuonThetaMaxCut=178.; - fNumberOfCollisions = 400; // Minimum bias Pb+Pb collisions - // -} -//_________________________________________________________________________ -AliGenMUONCocktail::AliGenMUONCocktail(const AliGenMUONCocktail & cocktail): - AliGenCocktail(cocktail) -{ -// Copy constructor - fTotalRate =0; - fNSucceded=0; - fNGenerated=0; - fMuonMultiplicity=1; - fMuonPtCut= 1.; - fMuonThetaMinCut=171.; - fMuonThetaMaxCut=178.; - fNumberOfCollisions = 400; // Minimum bias Pb+Pb collisions - // } //_________________________________________________________________________ AliGenMUONCocktail::~AliGenMUONCocktail() { // Destructor - delete fEntries; + if (fFastGlauber) delete fFastGlauber; } //_________________________________________________________________________ void AliGenMUONCocktail::Init() { - // Defining MUON physics cocktail + // NN cross section + Double_t sigmaReaction = 0.072; // MinBias NN cross section for PbPb LHC energies http://arxiv.org/pdf/nucl-ex/0302016 + + // Initialising Fast Glauber object + fFastGlauber = AliFastGlauber::Instance(); + fFastGlauber->SetPbPbLHC(); + fFastGlauber->SetNNCrossSection(sigmaReaction*1000.); //Expected NN cross-section in mb at LHC with diffractive part http://arxiv.org/pdf/nucl-ex/0302016 ) + fFastGlauber->Init(1); + + // Calculating average number of collisions + Int_t ib=0; + Int_t ibmax=10000; + Double_t b = 0.; + fAverageImpactParameter=0.; + fNumberOfCollisions = 0.; + fNumberOfParticipants = 0.; + for(ib=0; ibGetRandomImpactParameter(fLowImpactParameter,fHighImpactParameter); + fAverageImpactParameter+=b; + fNumberOfCollisions += fFastGlauber->GetNumberOfCollisions( b )/(1.-TMath::Exp(-fFastGlauber->GetNumberOfCollisions(b))); + fNumberOfParticipants += fFastGlauber->GetNumberOfParticipants( b ); + } + fAverageImpactParameter/= ((Double_t) ibmax); + fNumberOfCollisions /= ((Double_t) ibmax); + fNumberOfParticipants /= ((Double_t) ibmax);; + AliInfo(Form("=%4.2f, =%5.1f and and =%5.1f",b, fNumberOfCollisions, fNumberOfParticipants)); + // Estimating shadowing on charm a beaty production + // ----------------------------------------------------- + // Extrapolation of the cross sections from $p-p$ to \mbox{Pb--Pb} + // interactions + // is done by means of the Glauber model. For the impact parameter dependence + // of the shadowing factor we use a simple formula: + // $C_{sh}(b) = C_{sh}(0) + (1 - C_{sh}(0))(b/16~fm)4$, + // motivated by the theoretical predictions (see e.g. + // V. Emelyanov et al., Phys. Rev. C61, 044904 (2000)) and HIJING + // simulations showing an almost flat behaviour + // up to 10~$fm$ and a rapid increase to 1 for larger impact parameters. + // C_{sh}(0) = 0.60 for Psi and 0.76 for Upsilon (Smba communication). + // for open charm and beauty is 0.65 and 0.84 + // ----------------------------------------------------- + Double_t charmshadowing = 0.65 + (1.0-0.65)*TMath::Power(fAverageImpactParameter/16.,4); + Double_t beautyshadowing = 0.84 + (1.0-0.84)*TMath::Power(fAverageImpactParameter/16.,4); + Double_t charmoniumshadowing = 0.60 + (1.0-0.60)*TMath::Power(fAverageImpactParameter/16.,4); + Double_t beautoniumshadowing = 0.76 + (1.0-0.76)*TMath::Power(fAverageImpactParameter/16.,4); + if (fAverageImpactParameter>16.) { + charmoniumshadowing = 1.0; + beautoniumshadowing = 1.0; + charmshadowing = 1.0; + beautyshadowing= 1.0; + } + AliInfo(Form("Shadowing for charmonium and beautonium production are %4.2f and %4.2f respectively",charmoniumshadowing,beautoniumshadowing)); + AliInfo(Form("Shadowing for charm and beauty production are %4.2f and %4.2f respectively",charmshadowing,beautyshadowing)); + + // Defining MUON physics cocktail // Kinematical limits for particle generation - Float_t ptMin = fPtMin; - Float_t ptMax = fPtMax; - Float_t yMin = fYMin;; - Float_t yMax = fYMax;; - Float_t phiMin = fPhiMin*180./TMath::Pi(); - Float_t phiMax = fPhiMax*180./TMath::Pi(); - printf(">>> Kinematical range pT:%f:%f y:%f:%f Phi:%f:%f\n",ptMin,ptMax,yMin,yMax,phiMin,phiMax); + 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)); - Float_t sigma_reaction = 0.072; // MINB pp at LHC energies 72 mb - - // Generating J/Psi Physics - AliGenParam * gen_jpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "Vogt", "Jpsi"); - // 4pi Generation - gen_jpsi->SetPtRange(0,100.); - gen_jpsi->SetYRange(-8.,8); - gen_jpsi->SetPhiRange(0.,360.); - gen_jpsi->SetForceDecay(kDiMuon); - gen_jpsi->SetTrackingFlag(1); - // Calculation of the paritcle multiplicity per event in the muonic channel - Float_t ratio_jpsi; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail - Float_t sigma_jpsi = 31.0e-6 * 0.437; // section "6.7 Quarkonia Production" table 6.5 for pp times shadowing - Float_t br_jpsi = 0.0588; // Branching Ratio for JPsi - gen_jpsi->Init(); // Generating pT and Y parametrsation for the 4pi - ratio_jpsi = sigma_jpsi * br_jpsi * fNumberOfCollisions / sigma_reaction * gen_jpsi->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); - printf(">>> ratio jpsi %g et %g Ncol %g sigma %g\n",ratio_jpsi,gen_jpsi->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax),fNumberOfCollisions, sigma_jpsi ); + // Generating J/Psi Physics + // Using CFD scaled distribution (see http://clrwww.in2p3.fr/DIMUON2004/talks/sgrigoryan.pdf ) + AliGenParam * genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF scaled", "Jpsi"); + genjpsi->SetPtRange(0,100.); // 4pi generation + genjpsi->SetYRange(-8.,8); + genjpsi->SetPhiRange(0.,360.); + genjpsi->SetForceDecay(kDiMuon); + genjpsi->SetTrackingFlag(1); + // Calculation of the particle multiplicity per event in the muonic channel + Double_t ratiojpsi; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail + Double_t sigmajpsi = 31.0e-6 * charmoniumshadowing; // section "6.7 Quarkonia Production" table 6.5 for pp times shadowing + Double_t brjpsi = 0.0588; // Branching Ratio for JPsi PDG PRC15 (200) + genjpsi->Init(); // Generating pT and Y parametrsation for the 4pi + ratiojpsi = sigmajpsi * brjpsi * fNumberOfCollisions / sigmaReaction * genjpsi->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); + AliInfo(Form("Jpsi production cross-section in pp with shadowing %5.3g barns",sigmajpsi)); + AliInfo(Form("Jpsi production probability per collisions in acceptance via the muonic channel %5.3g",ratiojpsi)); // Generation in the kinematical limits of AliGenMUONCocktail - gen_jpsi->SetPtRange(ptMin, ptMax); - gen_jpsi->SetYRange(yMin, yMax); - gen_jpsi->SetPhiRange(phiMin, phiMax); - gen_jpsi->Init(); // Generating pT and Y parametrsation in the desired kinematic range + genjpsi->SetPtRange(ptMin, ptMax); + genjpsi->SetYRange(yMin, yMax); + genjpsi->SetPhiRange(phiMin, phiMax); + genjpsi->Init(); // Generating pT and Y parametrsation in the desired kinematic range // Adding Generator - AddGenerator(gen_jpsi, "Jpsi", ratio_jpsi); - fTotalRate+=ratio_jpsi; + AddGenerator(genjpsi, "Jpsi", ratiojpsi); + fTotalRate+=ratiojpsi; -// Generating Upsilon Physics - AliGenParam * gen_upsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "Vogt", "Upsilon"); - gen_upsilon->SetPtRange(0,100.); - gen_upsilon->SetYRange(-8.,8); - gen_upsilon->SetPhiRange(0.,360.); - gen_upsilon->SetForceDecay(kDiMuon); - gen_upsilon->SetTrackingFlag(1); - Float_t ratio_upsilon; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail - Float_t sigma_upsilon = 0.501e-6 * 0.674; // section "6.7 Quarkonia Production" table 6.5 for pp times shadowing - Float_t br_upsilon = 0.0248; // Branching Ratio for Upsilon - gen_upsilon->Init(); // Generating pT and Y parametrsation for the 4pi - ratio_upsilon = sigma_upsilon * br_upsilon * fNumberOfCollisions / sigma_reaction * gen_upsilon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); - printf(">>> ratio upsilon %g et %g\n",ratio_upsilon, gen_upsilon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax)); - gen_upsilon->SetPtRange(ptMin, ptMax); - gen_upsilon->SetYRange(yMin, yMax); - gen_upsilon->SetPhiRange(phiMin, phiMax); - gen_upsilon->Init(); // Generating pT and Y parametrsation in the desired kinematic range - AddGenerator(gen_upsilon,"Upsilon", ratio_upsilon); - fTotalRate+=ratio_upsilon; + // Generating Psi prime Physics + // Using CFD scaled distribution (see http://clrwww.in2p3.fr/DIMUON2004/talks/sgrigoryan.pdf ) + AliGenParam * genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF scaled", "PsiP"); + genpsiP->SetPtRange(0,100.);// 4pi generation + genpsiP->SetYRange(-8.,8); + genpsiP->SetPhiRange(0.,360.); + genpsiP->SetForceDecay(kDiMuon); + genpsiP->SetTrackingFlag(1); + // Calculation of the paritcle multiplicity per event in the muonic channel + Double_t ratiopsiP; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail + Double_t sigmapsiP = 4.68e-6 * charmoniumshadowing; // section "6.7 Quarkonia Production" table 6.5 for pp times shadowing + Double_t brpsiP = 0.0103; // Branching Ratio for PsiP + genpsiP->Init(); // Generating pT and Y parametrsation for the 4pi + ratiopsiP = sigmapsiP * brpsiP * fNumberOfCollisions / sigmaReaction * genpsiP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); + AliInfo(Form("Psi prime production cross-section in pp with shadowing %5.3g barns",sigmapsiP)); + AliInfo(Form("Psi prime production probability per collisions in acceptance via the muonic channel %5.3g",ratiopsiP)); + // Generation in the kinematical limits of AliGenMUONCocktail + genpsiP->SetPtRange(ptMin, ptMax); + genpsiP->SetYRange(yMin, yMax); + genpsiP->SetPhiRange(phiMin, phiMax); + genpsiP->Init(); // Generating pT and Y parametrsation in the desired kinematic range + // Adding Generator + AddGenerator(genpsiP, "PsiP", ratiopsiP); + fTotalRate+=ratiopsiP; -// Generating Charm Physics - AliGenParam * gen_charm = new AliGenParam(1, AliGenMUONlib::kCharm, "Vogt", "Charm"); - gen_charm->SetPtRange(0,100.); - gen_charm->SetYRange(-8.,8); - gen_charm->SetPhiRange(0.,360.); - gen_charm->SetForceDecay(kSemiMuonic); - gen_charm->SetTrackingFlag(1); - Float_t ratio_charm; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail - Float_t sigma_charm = 2. * 6.64e-3 * 0.65 ; // - Float_t br_charm = 0.12; // Branching Ratio for Charm - gen_charm->Init(); // Generating pT and Y parametrsation for the 4pi - ratio_charm = sigma_charm * br_charm * fNumberOfCollisions / sigma_reaction * - gen_charm->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); - gen_charm->SetPtRange(ptMin, ptMax); - gen_charm->SetYRange(yMin, yMax); - gen_charm->SetPhiRange(phiMin, phiMax); - gen_charm->Init(); // Generating pT and Y parametrsation in the desired kinematic range + // Generating Upsilon Physics + // Using CFD scaled distribution (see http://clrwww.in2p3.fr/DIMUON2004/talks/sgrigoryan.pdf ) + AliGenParam * genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "CDF scaled", "Upsilon"); + genupsilon->SetPtRange(0,100.); + genupsilon->SetYRange(-8.,8); + genupsilon->SetPhiRange(0.,360.); + genupsilon->SetForceDecay(kDiMuon); + genupsilon->SetTrackingFlag(1); + Double_t ratioupsilon; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail + Double_t sigmaupsilon = 0.501e-6 * beautoniumshadowing; // section "6.7 Quarkonia Production" table 6.5 for pp times shadowing + Double_t brupsilon = 0.0248; // Branching Ratio for Upsilon + genupsilon->Init(); // Generating pT and Y parametrsation for the 4pi + ratioupsilon = sigmaupsilon * brupsilon * fNumberOfCollisions / sigmaReaction * genupsilon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); + AliInfo(Form("Upsilon 1S production cross-section in pp with shadowing %5.3g barns",sigmaupsilon)); + AliInfo(Form("Upsilon 1S production probability per collisions in acceptance via the muonic channel %5.3g",ratioupsilon)); + genupsilon->SetPtRange(ptMin, ptMax); + genupsilon->SetYRange(yMin, yMax); + genupsilon->SetPhiRange(phiMin, phiMax); + genupsilon->Init(); // Generating pT and Y parametrsation in the desired kinematic range + AddGenerator(genupsilon,"Upsilon", ratioupsilon); + fTotalRate+=ratioupsilon; - printf(">>> ratio charm %f\n",ratio_charm); - AddGenerator(gen_charm,"Charm", ratio_charm); - fTotalRate+=ratio_charm; + // Generating UpsilonP Physics + // Using CFD scaled distribution (see http://clrwww.in2p3.fr/DIMUON2004/talks/sgrigoryan.pdf ) + AliGenParam * genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "CDF Scaled", "UpsilonP"); + genupsilonP->SetPtRange(0,100.); + genupsilonP->SetYRange(-8.,8); + genupsilonP->SetPhiRange(0.,360.); + genupsilonP->SetForceDecay(kDiMuon); + genupsilonP->SetTrackingFlag(1); + Double_t ratioupsilonP; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail + Double_t sigmaupsilonP = 0.246e-6 * beautoniumshadowing; // section "6.7 Quarkonia Production" table 6.5 for pp times shadowing + Double_t brupsilonP = 0.0131; // Branching Ratio for UpsilonP + genupsilonP->Init(); // Generating pT and Y parametrsation for the 4pi + ratioupsilonP = sigmaupsilonP * brupsilonP * fNumberOfCollisions / sigmaReaction * genupsilonP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); + AliInfo(Form("Upsilon 2S production cross-section in pp with shadowing %5.3g barns",sigmaupsilonP)); + AliInfo(Form("Upsilon 2S production probability per collisions in acceptance via the muonic channel %5.3g",ratioupsilonP)); + genupsilonP->SetPtRange(ptMin, ptMax); + genupsilonP->SetYRange(yMin, yMax); + genupsilonP->SetPhiRange(phiMin, phiMax); + genupsilonP->Init(); // Generating pT and Y parametrsation in the desired kinematic range + AddGenerator(genupsilonP,"UpsilonP", ratioupsilonP); + fTotalRate+=ratioupsilonP; -// Generating Beauty Physics "Correlated Pairs" - AliGenParam * gen_beauty = new AliGenParam(2, AliGenMUONlib::kBeauty, "Vogt", "Beauty"); - gen_beauty->SetPtRange(0,100.); - gen_beauty->SetYRange(-8.,8); - gen_beauty->SetPhiRange(0.,360.); - gen_beauty->SetForceDecay(kSemiMuonic); - gen_beauty->SetTrackingFlag(1); - Float_t ratio_beauty; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail - Float_t sigma_beauty = 2. * 0.210e-3 * 0.84; // - Float_t br_beauty = 0.15; // Branching Ratio for Beauty - gen_beauty->Init(); // Generating pT and Y parametrsation for the 4pi - ratio_beauty = sigma_beauty * br_beauty * fNumberOfCollisions / sigma_reaction * - gen_beauty->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); - gen_beauty->SetPtRange(ptMin, ptMax); - gen_beauty->SetYRange(yMin, yMax); - gen_beauty->SetPhiRange(phiMin, phiMax); - gen_beauty->Init(); // Generating pT and Y parametrisation in the desired kinematic range + // Generating UpsilonPP Physics + // Using CFD scaled distribution (see http://clrwww.in2p3.fr/DIMUON2004/talks/sgrigoryan.pdf ) + AliGenParam * genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF Scaled", "UpsilonPP"); + genupsilonPP->SetPtRange(0,100.); + genupsilonPP->SetYRange(-8.,8); + genupsilonPP->SetPhiRange(0.,360.); + genupsilonPP->SetForceDecay(kDiMuon); + genupsilonPP->SetTrackingFlag(1); + Double_t ratioupsilonPP; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail + Double_t sigmaupsilonPP = 0.100e-6 * beautoniumshadowing; // section "6.7 Quarkonia Production" table 6.5 for pp times shadowing + Double_t brupsilonPP = 0.0181; // Branching Ratio for UpsilonPP + genupsilonPP->Init(); // Generating pT and Y parametrsation for the 4pi + ratioupsilonPP = sigmaupsilonPP * brupsilonPP * fNumberOfCollisions / sigmaReaction * genupsilonPP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); + AliInfo(Form("Upsilon 3S production cross-section in pp with shadowing %5.3g barns",sigmaupsilonPP)); + AliInfo(Form("Upsilon 3S production probability per collisions in acceptance via the muonic channel %5.3g",ratioupsilonPP)); + genupsilonPP->SetPtRange(ptMin, ptMax); + genupsilonPP->SetYRange(yMin, yMax); + genupsilonPP->SetPhiRange(phiMin, phiMax); + genupsilonPP->Init(); // Generating pT and Y parametrsation in the desired kinematic range + AddGenerator(genupsilonPP,"UpsilonPP", ratioupsilonPP); + fTotalRate+=ratioupsilonPP; - printf(">>> ratio beauty %f\n",ratio_beauty); - AddGenerator(gen_beauty,"Beauty", ratio_beauty); - fTotalRate+=ratio_beauty; +// Generating non-correlated Charm Physics + AliGenParam * gencharm = new AliGenParam(1, AliGenMUONlib::kCharm, "pp", "Charm"); + gencharm->SetPtRange(0,100.); + gencharm->SetYRange(-8.,8); + gencharm->SetPhiRange(0.,360.); + gencharm->SetForceDecay(kSemiMuonic); + gencharm->SetTrackingFlag(1); + Double_t ratiocharm; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail + Double_t sigmacharm = 2. * 6.64e-3 * charmshadowing ; // + Double_t brcharm = 0.12; // Branching Ratio for Charm + gencharm->Init(); // Generating pT and Y parametrsation for the 4pi + ratiocharm = sigmacharm * brcharm * fNumberOfCollisions / sigmaReaction * gencharm->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); + AliInfo(Form("Charm production cross-section in pp with shadowing %5.3g barns",sigmacharm)); + AliInfo(Form("Charm production probability per collisions in acceptance via the semi-muonic channel %5.3g",ratiocharm)); + gencharm->SetPtRange(ptMin, ptMax); + gencharm->SetYRange(yMin, yMax); + gencharm->SetPhiRange(phiMin, phiMax); + gencharm->Init(); // Generating pT and Y parametrsation in the desired kinematic range + AddGenerator(gencharm,"Charm", ratiocharm); + fTotalRate+=ratiocharm; -// Generating Pion Physics - AliGenParam * gen_pion = new AliGenParam(1, AliGenMUONlib::kPion, "Vogt", "Pion"); - gen_pion->SetPtRange(0,100.); - gen_pion->SetYRange(-8.,8); - gen_pion->SetPhiRange(0.,360.); - gen_pion->SetForceDecay(kPiToMu); - gen_pion->SetTrackingFlag(1); - Float_t ratio_pion; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail - Float_t sigma_pion = 1.6e-2; // A ojo TO be studied in detail. - Float_t br_pion = 0.9999; // Branching Ratio for Pion - gen_pion->Init(); // Generating pT and Y parametrsation for the 4pi - ratio_pion = sigma_pion * br_pion * fNumberOfParticipants / sigma_reaction * gen_pion->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); - gen_pion->SetPtRange(ptMin, ptMax); - gen_pion->SetYRange(yMin, yMax); - gen_pion->SetPhiRange(phiMin, phiMax); - gen_pion->Init(); // Generating pT and Y parametrsation in the desired kinematic range - printf(">>> ratio pion %f\n",ratio_pion); - AddGenerator(gen_pion,"Pion", ratio_pion); - fTotalRate+=ratio_pion; +// Generating non-correlated Beauty Physics + AliGenParam * genbeauty = new AliGenParam(1, AliGenMUONlib::kBeauty, "pp", "Beauty"); + genbeauty->SetPtRange(0,100.); + genbeauty->SetYRange(-8.,8); + genbeauty->SetPhiRange(0.,360.); + genbeauty->SetForceDecay(kSemiMuonic); + genbeauty->SetTrackingFlag(1); + Double_t ratiobeauty; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail + Double_t sigmabeauty = 2. * 0.210e-3 * beautyshadowing; // + Double_t brbeauty = 0.15; // Branching Ratio for Beauty + genbeauty->Init(); // Generating pT and Y parametrsation for the 4pi + ratiobeauty = sigmabeauty * brbeauty * fNumberOfCollisions / sigmaReaction * + genbeauty->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); + AliInfo(Form("Beauty production cross-section in pp with shadowing %5.3g barns",sigmabeauty)); + AliInfo(Form("Beauty production probability per collisions in acceptance via the semi-muonic channel %5.3g",ratiobeauty)); + genbeauty->SetPtRange(ptMin, ptMax); + genbeauty->SetYRange(yMin, yMax); + genbeauty->SetPhiRange(phiMin, phiMax); + genbeauty->Init(); // Generating pT and Y parametrisation in the desired kinematic range + AddGenerator(genbeauty,"Beauty", ratiobeauty); + fTotalRate+=ratiobeauty; -// Generating Kaon Physics - AliGenParam * gen_kaon = new AliGenParam(1, AliGenMUONlib::kKaon, "Vogt", "Kaon"); - gen_kaon->SetPtRange(0,100.); - gen_kaon->SetYRange(-8.,8); - gen_kaon->SetPhiRange(0.,360.); - gen_kaon->SetForceDecay(kKaToMu); - gen_kaon->SetTrackingFlag(1); - Float_t ratio_kaon; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail - Float_t sigma_kaon = 1.8e-4; // OJO - Float_t br_kaon = 0.6351 ; // Branching Ratio for Kaon - gen_kaon->Init(); // Generating pT and Y parametrsation for the 4pi - ratio_kaon = sigma_kaon * br_kaon * fNumberOfParticipants/ sigma_reaction * gen_kaon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); - gen_kaon->SetPtRange(ptMin, ptMax); - gen_kaon->SetYRange(yMin, yMax); - gen_kaon->SetPhiRange(phiMin, phiMax); - gen_kaon->Init(); // Generating pT and Y parametrsation in the desired kinematic range - printf(">>> ratio kaon %f\n",ratio_kaon); - AddGenerator(gen_kaon,"Kaon", ratio_kaon); - fTotalRate+=ratio_kaon; + // Only if hadronic muons are included in the cocktail + if(fHadronicMuons) { + // Generating Pion Physics + // The scaling with Npart and Ncoll has been obtained to reproduced tha values presented by Valeri lors de presentatation + // a Clermont Ferrand http://pcrochet.home.cern.ch/pcrochet/files/valerie.pdf + // b range(fm) Ncoll Npart N_mu pT>0.4 GeV/c + // 0 - 3 1982 381 3.62 + // 3 - 6 1388 297 3.07 + // 6 - 9 674 177 1.76 + // 9 - 12 188 71 0.655 + // 12 - 16 15 10 0.086 + // We found the hadronic muons scales quite well with the number of participants + AliGenParam * genpion = new AliGenParam(1, AliGenMUONlib::kPion, "default", "Pion"); + genpion->SetPtRange(0,100.); + genpion->SetYRange(-8.,8); + genpion->SetPhiRange(0.,360.); + genpion->SetForceDecay(kPiToMu); + genpion->SetTrackingFlag(1); + Double_t ratiopion; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail + Double_t sigmapion = 1.80e-2; // Just for reproducing Valeries's data + Double_t brpion = 0.9999; // Branching Ratio for Pion + genpion->Init(); // Generating pT and Y parametrsation for the 4pi + ratiopion = sigmapion * brpion * (0.93*fNumberOfParticipants+0.07*fNumberOfCollisions) / sigmaReaction * genpion->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); + AliInfo(Form("Pseudo-Pion production cross-section in pp with shadowing %5.3g barns",sigmapion)); + AliInfo(Form("Pion production probability per collisions in acceptance via the muonic channel %5.3g",ratiopion)); + genpion->SetPtRange(ptMin, ptMax); + genpion->SetYRange(yMin, yMax); + genpion->SetPhiRange(phiMin, phiMax); + genpion->Init(); // Generating pT and Y parametrsation in the desired kinematic range + AddGenerator(genpion,"Pion", ratiopion); + fTotalRate+=ratiopion; + + // Generating Kaon Physics + AliGenParam * genkaon = new AliGenParam(1, AliGenMUONlib::kKaon, "default", "Kaon"); + genkaon->SetPtRange(0,100.); + genkaon->SetYRange(-8.,8); + genkaon->SetPhiRange(0.,360.); + genkaon->SetForceDecay(kKaToMu); + genkaon->SetTrackingFlag(1); + Double_t ratiokaon; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail + Double_t sigmakaon = 2.40e-4; // Valerie presentation Clermont-16-jan-2004 and Alice-int-2002-06 + Double_t brkaon = 0.6351 ; // Branching Ratio for Kaon + genkaon->Init(); // Generating pT and Y parametrsation for the 4pi + ratiokaon = sigmakaon * brkaon * (0.93*fNumberOfParticipants+0.07*fNumberOfCollisions)/ sigmaReaction * genkaon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); + AliInfo(Form("Pseudo-kaon production cross-section in pp with shadowing %5.3g barns",sigmakaon)); + AliInfo(Form("Kaon production probability per collisions in acceptance via the muonic channel %5.3g",ratiokaon)); + genkaon->SetPtRange(ptMin, ptMax); + genkaon->SetYRange(yMin, yMax); + genkaon->SetPhiRange(phiMin, phiMax); + genkaon->Init(); // Generating pT and Y parametrsation in the desired kinematic range + AddGenerator(genkaon,"Kaon", ratiokaon); + fTotalRate+=ratiokaon; + } } //_________________________________________________________________________ @@ -229,68 +365,88 @@ void AliGenMUONCocktail::Generate() AliGenCocktailEntry *entry = 0; AliGenCocktailEntry *preventry = 0; AliGenerator* gen = 0; + const TObjArray *partArray = gAlice->GetMCApp()->Particles(); - TObjArray *partArray = gAlice->GetMCApp()->Particles(); - -// -// Generate the vertex position used by all generators -// +// Generate the vertex position used by all generators if(fVertexSmear == kPerEvent) Vertex(); - Bool_t PrimordialTrigger = kFALSE; - while(!PrimordialTrigger) { + // Loop on primordialTrigger + Bool_t primordialTrigger = kFALSE; + while(!primordialTrigger) { + //Reseting stack + AliRunLoader * runloader = AliRunLoader::Instance(); + if (runloader) + if (runloader->Stack()) + runloader->Stack()->Clean(); + // Loop over generators and generate events + Int_t igen=0; + Int_t npart =0; + while((entry = (AliGenCocktailEntry*)next())) { + gen = entry->Generator(); + gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2)); + if ( (npart = gRandom->Poisson(entry->Rate())) >0 ) { + igen++; + if (igen == 1) entry->SetFirst(0); + else entry->SetFirst((partArray->GetEntriesFast())+1); + // if ( (fHadronicMuons == kFALSE) && ( (gen->GetName() == "Pions") || (gen->GetName() == "Kaons") ) ) + // { AliInfo(Form("This generator %s is finally not generated. This is option for hadronic muons.",gen->GetName() ) ); } + // else { + gen->SetNumberParticles(npart); + gen->Generate(); + entry->SetLast(partArray->GetEntriesFast()); + preventry = entry; + // } + } + } + next.Reset(); + // Testing primordial trigger : Muon pair in the MUON spectrometer acceptance and pTCut + Int_t iPart; + fNGenerated++; + Int_t numberOfMuons=0; + // printf(">>>fNGenerated is %d\n",fNGenerated); + + TObjArray GoodMuons; // Used in the Invariant Mass selection cut + + for(iPart=0; iPartGetEntriesFast(); iPart++){ + + if ( (TMath::Abs(gAlice->GetMCApp()->Particle(iPart)->GetPdgCode())==13) && + (gAlice->GetMCApp()->Particle(iPart)->Theta()*180./TMath::Pi()>fMuonThetaMinCut) && + (gAlice->GetMCApp()->Particle(iPart)->Theta()*180./TMath::Pi()GetMCApp()->Particle(iPart)->Pt()>fMuonPtCut) ) { + gAlice->GetMCApp()->Particle(iPart)->SetProductionVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2), 0.); + GoodMuons.AddLast(gAlice->GetMCApp()->Particle(iPart)); + numberOfMuons++; + } + } + + // Test the invariant mass of each pair (if cut on Invariant mass is required) + Bool_t InvMassRangeOK = kTRUE; + if(fInvMassCut && (numberOfMuons>=2) ){ + TLorentzVector fV1, fV2, fVtot; + InvMassRangeOK = kFALSE; + for(iPart=0; iPartGetRunLoader(); - if (runloader) - if (runloader->Stack()) - runloader->Stack()->Reset(); - // - // Loop over generators and generate events - Int_t igen=0; - Int_t npart =0; - - while((entry = (AliGenCocktailEntry*)next())) { - gen = entry->Generator(); - gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2)); - if ( (npart = gRandom->Poisson(entry->Rate())) >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; + for(int iPart2=iPart+1; iPart2Px() ,mu1->Py() ,mu1->Pz() ,mu1->Energy() ); + fV2.SetPxPyPzE(mu2->Px() ,mu2->Py() ,mu2->Pz() ,mu2->Energy() ); + fVtot = fV1 + fV2; + + if(fVtot.M()>fInvMassMinCut && fVtot.M()= fMuonMultiplicity) && InvMassRangeOK ) primordialTrigger = kTRUE; } - } - next.Reset(); - - // Tesitng primordial trigger : Muon pair in the MUON spectrometer acceptance 171,178 and pTCut - Int_t iPart; - fNGenerated++; - Int_t numberOfMuons=0; - // printf(">>>fNGenerated is %d\n",fNGenerated); - for(iPart=0; iPartGetEntriesFast(); iPart++){ - // gAlice->GetMCApp()->Particle(iPart)->Print(); - if ( (TMath::Abs(gAlice->GetMCApp()->Particle(iPart)->GetPdgCode())==13) && - (gAlice->GetMCApp()->Particle(iPart)->Theta()*180./TMath::Pi()>fMuonThetaMinCut) && - (gAlice->GetMCApp()->Particle(iPart)->Theta()*180./TMath::Pi()GetMCApp()->Particle(iPart)->Pt()>fMuonPtCut) ) { - gAlice->GetMCApp()->Particle(iPart)->SetProductionVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2), 0.); - numberOfMuons++; - } - } - // printf(">>> Number of Muons is %d \n", numberOfMuons); - if (numberOfMuons >= fMuonMultiplicity ) PrimordialTrigger = kTRUE; - } - //printf(">>> Trigger Accepted!!! \n"); fNSucceded++; - // Float_t Ratio = (Float_t) fNSucceded/fNGenerated; - // printf("Generated %d, Succeded %d and Ratio %f\n",fNGenerated, fNSucceded,Ratio); -} - - - - - + AliDebug(5,Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded)); +}