X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;ds=sidebyside;f=EVGEN%2FAliGenMUONCocktailpp.cxx;h=f71a6b1137e73536217a6ae26bd5873f61049db7;hb=33c3c91a040995417ea28c5ae034899b77cca3a6;hp=b3aef2c439789d5101de6478e4f34ee7c0931355;hpb=93a2041b6acdbf31d3a3eb48708a2a54a759543d;p=u%2Fmrichter%2FAliRoot.git diff --git a/EVGEN/AliGenMUONCocktailpp.cxx b/EVGEN/AliGenMUONCocktailpp.cxx index b3aef2c4397..f71a6b1137e 100644 --- a/EVGEN/AliGenMUONCocktailpp.cxx +++ b/EVGEN/AliGenMUONCocktailpp.cxx @@ -22,11 +22,14 @@ // 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 #include #include #include +#include #include "AliGenCocktailEventHeader.h" @@ -37,20 +40,24 @@ #include "AliMC.h" #include "AliRun.h" #include "AliStack.h" +#include "AliDecayer.h" #include "AliLog.h" -#include "AliGenPythia.h" - +#include "AliGenCorrHF.h" ClassImp(AliGenMUONCocktailpp) //________________________________________________________________________ AliGenMUONCocktailpp::AliGenMUONCocktailpp() :AliGenCocktail(), - fTotalRate(0), + fDecayer(0), + fDecayModeResonance(kAll), + fDecayModePythia(kAll), fMuonMultiplicity(0), fMuonPtCut(0.), + fMuonPCut(0.), fMuonThetaMinCut(0.), - fMuonThetaMaxCut(0.), + fMuonThetaMaxCut(180.), + fMuonOriginCut(-999.), fNSucceded(0), fNGenerated(0) { @@ -65,7 +72,7 @@ AliGenMUONCocktailpp::~AliGenMUONCocktailpp() } //_________________________________________________________________________ -void AliGenMUONCocktailpp::Init() +void AliGenMUONCocktailpp::CreateCocktail() { // MinBias NN cross section @ pp 14 TeV -PR Vol II p:473 Double_t sigmaReaction = 0.070; @@ -87,7 +94,9 @@ void AliGenMUONCocktailpp::Init() Double_t ratioupsilon; Double_t ratioupsilonP; Double_t ratioupsilonPP; - + Double_t ratioccbar; + Double_t ratiobbbar; + // Cross sections in barns (from PPR Vol. II p: 552) pp - 14 TeV and // corrected from feed down of higher resonances @@ -96,7 +105,54 @@ void AliGenMUONCocktailpp::Init() 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; + + Bool_t showSIGJPSI = kTRUE; + if (gSystem->Getenv("COCKTAIL_SIGJPSI")) { + sigmajpsi = atof(gSystem->Getenv("COCKTAIL_SIGJPSI")); + AliInfo("Cross-section for JPsi set from external value"); + if (!gSystem->Getenv("COCKTAIL_SIGSHOW")) showSIGJPSI = kFALSE; + } + Bool_t showSIGPSIP = kTRUE; + if (gSystem->Getenv("COCKTAIL_SIGPSIP")) { + sigmapsiP = atof(gSystem->Getenv("COCKTAIL_SIGPSIP")); + AliInfo("Cross-section for Psi-prime set from external value"); + if (!gSystem->Getenv("COCKTAIL_SIGSHOW")) showSIGPSIP = kFALSE; + } + Bool_t showSIGUPSI = kTRUE; + if (gSystem->Getenv("COCKTAIL_SIGUPSI")) { + sigmaupsilon = atof(gSystem->Getenv("COCKTAIL_SIGUPSI")); + AliInfo("Cross-section for Upsilon 1S set from external value"); + if (!gSystem->Getenv("COCKTAIL_SIGSHOW")) showSIGUPSI = kFALSE; + } + Bool_t showSIGUPSIP = kTRUE; + if (gSystem->Getenv("COCKTAIL_SIGUPSIP")) { + sigmaupsilonP = atof(gSystem->Getenv("COCKTAIL_SIGUPSIP")); + AliInfo("Cross-section for Upsilon 2S set from external value"); + if (!gSystem->Getenv("COCKTAIL_SIGSHOW")) showSIGUPSIP = kFALSE; + } + Bool_t showSIGUPSIPP = kTRUE; + if (gSystem->Getenv("COCKTAIL_SIGUPSIPP")) { + sigmaupsilonPP = atof(gSystem->Getenv("COCKTAIL_SIGUPSIPP")); + AliInfo("Cross-section for Upsilon 3S set from external value"); + if (!gSystem->Getenv("COCKTAIL_SIGSHOW")) showSIGUPSIPP = kFALSE; + } + Bool_t showSIGCCBAR = kTRUE; + if (gSystem->Getenv("COCKTAIL_SIGCCBAR")) { + sigmaccbar = atof(gSystem->Getenv("COCKTAIL_SIGCCBAR")); + AliInfo("Cross-section for c-cbar set from external value"); + if (!gSystem->Getenv("COCKTAIL_SIGSHOW")) showSIGCCBAR = kFALSE; + } + Bool_t showSIGBBBAR = kTRUE; + if (gSystem->Getenv("COCKTAIL_SIGBBBAR")) { + sigmabbbar = atof(gSystem->Getenv("COCKTAIL_SIGBBBAR")); + AliInfo("Cross-section for b-bbar set from external value"); + if (!gSystem->Getenv("COCKTAIL_SIGSHOW")) showSIGBBBAR = kFALSE; + } + + AliInfo(Form("the parametrised resonances uses the decay mode %d",fDecayModeResonance)); + // Generation using CDF scaled distribution for pp@14TeV (from D. Stocco) //---------------------------------------------------------------------- // Jpsi generator @@ -105,19 +161,21 @@ void AliGenMUONCocktailpp::Init() genjpsi->SetPtRange(0.,100.); genjpsi->SetYRange(-8.,8.); genjpsi->SetPhiRange(0.,360.); - genjpsi->SetForceDecay(kAll); + 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)); + if (showSIGJPSI) { + 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; - //------------------------------------------------------------------ // Psi prime generator AliGenParam * genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF pp", "PsiP"); @@ -125,19 +183,21 @@ void AliGenMUONCocktailpp::Init() genpsiP->SetPtRange(0.,100.); genpsiP->SetYRange(-8.,8.); genpsiP->SetPhiRange(0.,360.); - genpsiP->SetForceDecay(kAll); + 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)); + if (showSIGPSIP) { + 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; - //------------------------------------------------------------------ // Upsilon 1S generator AliGenParam * genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "CDF pp", "Upsilon"); @@ -145,19 +205,20 @@ void AliGenMUONCocktailpp::Init() genupsilon->SetPtRange(0.,100.); genupsilon->SetYRange(-8.,8.); genupsilon->SetPhiRange(0.,360.); - genupsilon->SetForceDecay(kAll); + 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)); + if (showSIGUPSI) { + 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)); + } // 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"); @@ -165,18 +226,20 @@ void AliGenMUONCocktailpp::Init() genupsilonP->SetPtRange(0.,100.); genupsilonP->SetYRange(-8.,8.); genupsilonP->SetPhiRange(0.,360.); - genupsilonP->SetForceDecay(kAll); + 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)); + if (showSIGUPSIP) { + 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; //------------------------------------------------------------------ // Upsilon 3S generator @@ -185,34 +248,83 @@ void AliGenMUONCocktailpp::Init() genupsilonPP->SetPtRange(0.,100.); genupsilonPP->SetYRange(-8.,8.); genupsilonPP->SetPhiRange(0.,360.); - genupsilonPP->SetForceDecay(kAll); + 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)); + if (showSIGUPSIPP) { + 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; + +//------------------------------------------------------------------ +// Generator of charm + AliGenCorrHF *gencharm = new AliGenCorrHF(1, 4); + gencharm->SetMomentumRange(0,9999); + gencharm->SetForceDecay(kAll); + ratioccbar = sigmaccbar/sigmaReaction; + if (!gMC) gencharm->SetDecayer(fDecayer); + gencharm->Init(); + if (showSIGCCBAR) { + AliInfo(Form("c-cbar prod. cross-section in pp(14 TeV) %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; + if (!gMC) genbeauty->SetDecayer(fDecayer); + genbeauty->Init(); + if (showSIGBBBAR) { + AliInfo(Form("b-bbar prod. cross-section in pp(14 TeV) %5.3g b",sigmabbbar)); + AliInfo(Form("b-bbar prod. probability per collision in acceptance %5.3g",ratiobbbar)); + } + AddGenerator(genbeauty,"CorrHFBeauty",ratiobbbar); + +//------------------------------------------------------------------- // Pythia generator - AliGenPythia *pythia = new AliGenPythia(1); - pythia->SetProcess(kPyMbMSEL1); - pythia->SetStrucFunc(kCTEQ5L); - pythia->SetEnergyCMS(14000.); - pythia->SetForceDecay(kAll); - 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.; +// +// 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) { + while((entry = (AliGenCocktailEntry*)next())) { + entry->Generator()->SetStack(fStack); + } + } } //_________________________________________________________________________ @@ -224,7 +336,7 @@ void AliGenMUONCocktailpp::Generate() AliGenCocktailEntry *preventry = 0; AliGenerator* gen = 0; - TObjArray *partArray = gAlice->GetMCApp()->Particles(); + const TObjArray *partArray = gAlice->GetMCApp()->Particles(); // Generate the vertex position used by all generators if(fVertexSmear == kPerEvent) Vertex(); @@ -235,10 +347,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; @@ -268,24 +380,25 @@ void AliGenMUONCocktailpp::Generate() // in the muon spectrometer acceptance Int_t iPart; fNGenerated++; - Int_t numberOfMuons=0; - - for(iPart=0; iPartGetEntriesFast(); iPart++){ + Int_t numberOfMuons=0;Int_t maxPart = partArray->GetEntriesFast(); + for(iPart=0; iPartGetMCApp()->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) ) { - - numberOfMuons++; + TParticle *part = gAlice->GetMCApp()->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; } fNSucceded++; - AliInfo(Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded)); +// 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)); }