]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EVGEN/AliGenMUONCocktailpp.cxx
First cleanup of Loaders
[u/mrichter/AliRoot.git] / EVGEN / AliGenMUONCocktailpp.cxx
index c2ce74af54db1fcc42e56968bf5f02e5b62e92c3..476c3c67432064da7d64ffaec1b69810262d781c 100644 (file)
 // 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
@@ -77,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; 
@@ -99,26 +94,33 @@ void AliGenMUONCocktailpp::Init()
     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));
@@ -129,16 +131,16 @@ void AliGenMUONCocktailpp::Init()
     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));
@@ -149,16 +151,15 @@ void AliGenMUONCocktailpp::Init()
     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));
@@ -169,16 +170,15 @@ void AliGenMUONCocktailpp::Init()
     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));
@@ -189,16 +189,16 @@ void AliGenMUONCocktailpp::Init()
     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));
@@ -209,22 +209,61 @@ void AliGenMUONCocktailpp::Init()
     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);
+       }
+    }
 }
 
 //_________________________________________________________________________
@@ -236,7 +275,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();
@@ -247,10 +286,10 @@ void AliGenMUONCocktailpp::Generate()
     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;        
@@ -280,24 +319,25 @@ void AliGenMUONCocktailpp::Generate()
 // 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));
 }