]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EVGEN/AliGenMUONCocktailpp.cxx
Load pythia libraries.
[u/mrichter/AliRoot.git] / EVGEN / AliGenMUONCocktailpp.cxx
index c159b60d908c8e155db1c751eab38f192c8b5c43..709eecab428ffaaff3c855fea893c10564b02b39 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(),
-     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)
 {
 // Constructor
 }
-//_________________________________________________________________________
-AliGenMUONCocktailpp::AliGenMUONCocktailpp(const AliGenMUONCocktailpp & cocktail):
-    AliGenCocktail(cocktail),
-    fTotalRate(0),
-    fMuonMultiplicity(0),
-    fMuonPtCut(0.),
-    fMuonThetaMinCut(0.),
-    fMuonThetaMaxCut(0.),
-    fNSucceded(0),
-    fNGenerated(0)
-{
-// Copy constructor
-  fTotalRate =0;
-  fNSucceded=0; 
-  fNGenerated=0;
-  fMuonMultiplicity=0;
-  fMuonPtCut=0.;
-  fMuonThetaMinCut=0.; 
-  fMuonThetaMaxCut=0.;
-}
+
 //_________________________________________________________________________
 AliGenMUONCocktailpp::~AliGenMUONCocktailpp()
 // Destructor
@@ -84,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; 
@@ -106,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 
 
@@ -115,7 +105,11 @@ 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;
     
+    AliInfo(Form("the parametrised resonances uses the decay mode %d",fDecayModeResonance));
+
 // Generation using CDF scaled distribution for pp@14TeV (from D. Stocco)
 //----------------------------------------------------------------------
 // Jpsi generator
@@ -124,7 +118,9 @@ 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));
@@ -135,8 +131,6 @@ 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 pp", "PsiP");
@@ -144,7 +138,9 @@ 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));
@@ -155,8 +151,6 @@ 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 pp", "Upsilon");
@@ -164,7 +158,8 @@ 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));
@@ -175,8 +170,6 @@ 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 pp", "UpsilonP");
@@ -184,7 +177,8 @@ 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));
@@ -195,7 +189,6 @@ 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
@@ -204,7 +197,8 @@ 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));
@@ -215,23 +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->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);
+       }
+    }
 }
 
 //_________________________________________________________________________
@@ -257,7 +289,7 @@ void AliGenMUONCocktailpp::Generate()
        AliRunLoader * runloader = gAlice->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;        
@@ -287,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));
 }