More protections in case of null eloss

[u/mrichter/AliRoot.git] / EVGEN / AliGenMUONCocktailpp.cxx

diff --git a/EVGEN/AliGenMUONCocktailpp.cxx b/EVGEN/AliGenMUONCocktailpp.cxx
index e96d1e4a69c2544d4b6f7211218008a015b36d4e..709eecab428ffaaff3c855fea893c10564b02b39 100644 (file)
--- a/EVGEN/AliGenMUONCocktailpp.cxx
+++ b/EVGEN/AliGenMUONCocktailpp.cxx
@@ -22,7 +22,9 @@
 // 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>
@@ -40,8 +42,7 @@
 #include "AliStack.h"
 #include "AliDecayer.h"
 #include "AliLog.h"
-#include "AliGenPythia.h"
-
+#include "AliGenCorrHF.h"
 
 ClassImp(AliGenMUONCocktailpp)  
   
@@ -51,11 +52,11 @@ AliGenMUONCocktailpp::AliGenMUONCocktailpp()
      fDecayer(0),
      fDecayModeResonance(kAll),
      fDecayModePythia(kAll),
-     fTotalRate(0),
      fMuonMultiplicity(0),
      fMuonPtCut(0.),
+     fMuonPCut(0.),
      fMuonThetaMinCut(0.),
-     fMuonThetaMaxCut(0.),
+     fMuonThetaMaxCut(180.),
      fMuonOriginCut(-999.),
      fNSucceded(0),
      fNGenerated(0)
@@ -71,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; 
@@ -93,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 
 
@@ -102,6 +105,8 @@ 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));
 
@@ -126,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");
@@ -148,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");
@@ -169,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");
@@ -190,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
@@ -211,23 +209,54 @@ 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.);
-    AliInfo(Form("\n\npythia uses the decay mode %d",fDecayModePythia));
-    pythia->SetForceDecay(fDecayModePythia);
-    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) {
@@ -260,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;        
@@ -290,9 +319,7 @@ void AliGenMUONCocktailpp::Generate()
 // in the muon spectrometer acceptance
        Int_t iPart;
        fNGenerated++;
-       Int_t numberOfMuons=0;  
-
-       Int_t maxPart = partArray->GetEntriesFast();
+       Int_t numberOfMuons=0;Int_t maxPart = partArray->GetEntriesFast();
        for(iPart=0; iPart<maxPart; iPart++){      
            
          TParticle *part = gAlice->GetMCApp()->Particle(iPart);
@@ -300,7 +327,8 @@ void AliGenMUONCocktailpp::Generate()
            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->Pt()>fMuonPtCut) &&
+              (part->P()>fMuonPCut)) {
              numberOfMuons++;
            }
          }