// 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 <TObjArray.h>
#include <TParticle.h>
#include <TF1.h>
+#include <TVirtualMC.h>
#include "AliGenCocktailEventHeader.h"
#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()
+ :AliGenCocktail(),
+ fDecayer(0),
+ fDecayModeResonance(kAll),
+ fDecayModePythia(kAll),
+ fMuonMultiplicity(0),
+ fMuonPtCut(0.),
+ fMuonPCut(0.),
+ fMuonThetaMinCut(0.),
+ fMuonThetaMaxCut(180.),
+ fMuonOriginCut(-999.),
+ fNSucceded(0),
+ fNGenerated(0)
{
// Constructor
- fTotalRate =0;
- fNSucceded=0;
- fNGenerated=0;
- fMuonMultiplicity=0;
- fMuonPtCut= 0.;
- fMuonThetaMinCut=0.;
- fMuonThetaMaxCut=0.;
-}
-//_________________________________________________________________________
-AliGenMUONCocktailpp::AliGenMUONCocktailpp(const AliGenMUONCocktailpp & cocktail):
- AliGenCocktail(cocktail)
-{
-// Copy constructor
- fTotalRate =0;
- fNSucceded=0;
- fNGenerated=0;
- fMuonMultiplicity=0;
- fMuonPtCut=0.;
- fMuonThetaMinCut=0.;
- fMuonThetaMaxCut=0.;
}
+
//_________________________________________________________________________
AliGenMUONCocktailpp::~AliGenMUONCocktailpp()
// Destructor
}
//_________________________________________________________________________
-void AliGenMUONCocktailpp::Init()
+void AliGenMUONCocktailpp::CreateCocktail()
{
// MinBias NN cross section @ pp 14 TeV -PR Vol II p:473
Double_t sigmaReaction = 0.070;
Double_t ratioupsilon;
Double_t ratioupsilonP;
Double_t ratioupsilonPP;
-
-// Cross sections in barns (from PPR Vol. II p: 552) pp - 14 TeV
+ 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
- Double_t sigmajpsi = 53.85e-6;
+ Double_t sigmajpsi = 49.44e-6;
Double_t sigmapsiP = 7.67e-6;
- Double_t sigmaupsilon = 1.15e-6;
- Double_t sigmaupsilonP = 0.526e-6;
+ Double_t sigmaupsilon = 0.989e-6;
+ Double_t sigmaupsilonP = 0.502e-6;
Double_t sigmaupsilonPP = 0.228e-6;
+ Double_t sigmaccbar = 11.2e-3;
+ Double_t sigmabbbar = 0.51e-3;
-// Generation using CDF scaled distribution
-//(See http://clrwww.in2p3.fr/DIMUON2004/talks/sgrigoryan.pdf )
-
-//------------------------------------------------------------------
+ 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 scaled", "Jpsi");
+ 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(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));
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 scaled", "PsiP");
+ AliGenParam * genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF pp", "PsiP");
// first step: generation in 4pi
genpsiP->SetPtRange(0.,100.);
genpsiP->SetYRange(-8.,8.);
genpsiP->SetPhiRange(0.,360.);
- genpsiP->SetForceDecay(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));
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 scaled", "Upsilon");
+ 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(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));
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 scaled", "UpsilonP");
+ 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(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));
genupsilonP->SetPhiRange(phiMin, phiMax);
genupsilonP->Init(); // generation in selected kinematical range
AddGenerator(genupsilonP,"UpsilonP", ratioupsilonP); // Adding Generator
- fTotalRate+=ratioupsilonP;
//------------------------------------------------------------------
// Upsilon 3S generator
- AliGenParam * genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF scaled", "UpsilonPP");
+ 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(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));
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();
+ 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();
+ 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->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);
+ }
+ }
}
//_________________________________________________________________________
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();
Bool_t primordialTrigger = kFALSE;
while(!primordialTrigger) {
//Reseting stack
- AliRunLoader * runloader = gAlice->GetRunLoader();
+ AliRunLoader * runloader = AliRunLoader::GetRunLoader();
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;
// in the muon spectrometer acceptance
Int_t iPart;
fNGenerated++;
- Int_t numberOfMuons=0;
-
- for(iPart=0; iPart<partArray->GetEntriesFast(); iPart++){
+ Int_t numberOfMuons=0;Int_t maxPart = partArray->GetEntriesFast();
+ for(iPart=0; iPart<maxPart; 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()<fMuonThetaMaxCut) &&
- (gAlice->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()<fMuonThetaMaxCut) &&
+ (part->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));
}