New classe needed for PDC06 muon generation. (Nicole Bastid)
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 20 Mar 2006 09:12:10 +0000 (09:12 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 20 Mar 2006 09:12:10 +0000 (09:12 +0000)
EVGEN/AliGenMUONCocktailpp.cxx [new file with mode: 0644]
EVGEN/AliGenMUONCocktailpp.h [new file with mode: 0644]
EVGEN/EVGENLinkDef.h
EVGEN/libEVGEN.pkg

diff --git a/EVGEN/AliGenMUONCocktailpp.cxx b/EVGEN/AliGenMUONCocktailpp.cxx
new file mode 100644 (file)
index 0000000..c2ce74a
--- /dev/null
@@ -0,0 +1,304 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id$ */
+
+//
+// Classe to create the coktail for physics with muons for pp at 14 TeV
+// The followoing sources: 
+// jpsi,psiP, upsilon, upsilonP, upsilonPP are added to Pythia
+// The free parameeters are :
+//      pp reaction cross-section
+//      production cross-sections in pp collisions 
+// 
+
+#include <TObjArray.h>
+#include <TParticle.h>
+#include <TF1.h>
+
+#include "AliGenCocktailEventHeader.h"
+
+#include "AliGenCocktailEntry.h"
+#include "AliGenMUONCocktailpp.h"
+#include "AliGenMUONlib.h"
+#include "AliGenParam.h"
+#include "AliMC.h"
+#include "AliRun.h"
+#include "AliStack.h"
+#include "AliLog.h"
+#include "AliGenPythia.h"
+
+
+ClassImp(AliGenMUONCocktailpp)  
+  
+//________________________________________________________________________
+AliGenMUONCocktailpp::AliGenMUONCocktailpp()
+  :AliGenCocktail()
+{
+// 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()
+{
+// MinBias NN cross section @ pp 14 TeV -PR  Vol II p:473
+    Double_t sigmaReaction =   0.070; 
+    
+// These limits are only used for renormalization of quarkonia cross section
+// Pythia events are generated in 4pi  
+    Double_t ptMin  = fPtMin;
+    Double_t ptMax  = fPtMax;
+    Double_t yMin   = fYMin;;
+    Double_t yMax   = fYMax;;
+    Double_t phiMin = fPhiMin*180./TMath::Pi();
+    Double_t phiMax = fPhiMax*180./TMath::Pi();
+    AliInfo(Form("Ranges pT:%4.1f : %4.1f GeV/c, y:%4.2f : %4.2f, Phi:%5.1f : %5.1f degres",ptMin,ptMax,yMin,yMax,phiMin,phiMax));
+    
+// Ratios with respect to the reaction cross-section in the 
+// kinematics limit of the MUONCocktail
+    Double_t ratiojpsi;
+    Double_t ratiopsiP;
+    Double_t ratioupsilon;
+    Double_t ratioupsilonP;
+    Double_t ratioupsilonPP;
+    
+// Cross sections in barns (from PPR Vol. II p: 552) pp - 14 TeV
+
+    Double_t sigmajpsi = 53.85e-6;  
+    Double_t sigmapsiP = 7.67e-6;  
+    Double_t sigmaupsilon = 1.15e-6;  
+    Double_t sigmaupsilonP = 0.526e-6;  
+    Double_t sigmaupsilonPP = 0.228e-6;
+    
+// Generation using CDF scaled distribution
+//(See http://clrwww.in2p3.fr/DIMUON2004/talks/sgrigoryan.pdf )
+    
+//------------------------------------------------------------------
+// Jpsi generator
+    AliGenParam * genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF scaled", "Jpsi");
+// first step: generation in 4pi
+    genjpsi->SetPtRange(0.,100.);
+    genjpsi->SetYRange(-8.,8.);
+    genjpsi->SetPhiRange(0.,360.);
+    genjpsi->SetForceDecay(kAll);
+    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));
+// 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 scaled", "PsiP");
+// first step: generation in 4pi
+    genpsiP->SetPtRange(0.,100.);
+    genpsiP->SetYRange(-8.,8.);
+    genpsiP->SetPhiRange(0.,360.);
+    genpsiP->SetForceDecay(kAll);
+    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));
+// 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 scaled", "Upsilon");
+// first step: generation in 4pi
+    genupsilon->SetPtRange(0.,100.); 
+    genupsilon->SetYRange(-8.,8.);
+    genupsilon->SetPhiRange(0.,360.);
+    genupsilon->SetForceDecay(kAll);
+    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));
+// 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 scaled", "UpsilonP");
+// first step: generation in 4pi
+    genupsilonP->SetPtRange(0.,100.);
+    genupsilonP->SetYRange(-8.,8.);
+    genupsilonP->SetPhiRange(0.,360.);
+    genupsilonP->SetForceDecay(kAll);
+    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));
+// 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
+    AliGenParam * genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF scaled", "UpsilonPP");
+// first step: generation in 4pi
+    genupsilonPP->SetPtRange(0.,100.); 
+    genupsilonPP->SetYRange(-8.,8.);
+    genupsilonPP->SetPhiRange(0.,360.);
+    genupsilonPP->SetForceDecay(kAll);
+    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));
+// 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;
+    
+//------------------------------------------------------------------
+// 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.;
+
+}
+
+//_________________________________________________________________________
+void AliGenMUONCocktailpp::Generate()
+{
+// Generate event 
+    TIter next(fEntries);
+    AliGenCocktailEntry *entry = 0;
+    AliGenCocktailEntry *preventry = 0;
+    AliGenerator* gen = 0;
+
+    TObjArray *partArray = gAlice->GetMCApp()->Particles();
+    
+// Generate the vertex position used by all generators    
+    if(fVertexSmear == kPerEvent) Vertex();
+
+// Loop on primordialTrigger: 
+// minimum muon multiplicity above a pt cut in a theta acceptance region
+
+    Bool_t primordialTrigger = kFALSE;
+    while(!primordialTrigger) {
+       //Reseting stack
+       AliRunLoader * runloader = gAlice->GetRunLoader();
+       if (runloader)
+           if (runloader->Stack())
+               runloader->Stack()->Reset();
+       // Loop over generators and generate events
+       Int_t igen = 0;
+       Int_t npart = 0;        
+       const char* genName = 0;
+       while((entry = (AliGenCocktailEntry*)next())) {
+           gen = entry->Generator();
+           genName = entry->GetName();
+           gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
+
+           npart = (strcmp(genName,"Pythia") == 0) ? 1 :
+               gRandom->Poisson(entry->Rate());
+
+           if (npart > 0) {
+               igen++; 
+               if (igen == 1) entry->SetFirst(0);              
+               else  entry->SetFirst((partArray->GetEntriesFast())+1);
+               
+               gen->SetNumberParticles(npart);         
+               gen->Generate();
+               entry->SetLast(partArray->GetEntriesFast());
+               preventry = entry;
+           }
+       }  
+       next.Reset();
+
+// Testing primordial trigger: Single muons or dimuons with Pt above a Pt cut 
+// in the muon spectrometer acceptance
+       Int_t iPart;
+       fNGenerated++;
+       Int_t numberOfMuons=0;  
+       
+       for(iPart=0; iPart<partArray->GetEntriesFast(); 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++;
+           }
+       }
+       
+       if (numberOfMuons >= fMuonMultiplicity) primordialTrigger = kTRUE;
+    }
+    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));
+}
+
+    
diff --git a/EVGEN/AliGenMUONCocktailpp.h b/EVGEN/AliGenMUONCocktailpp.h
new file mode 100644 (file)
index 0000000..cc545f3
--- /dev/null
@@ -0,0 +1,52 @@
+#ifndef ALIGENMUONCOCKTAILPP_H
+#define ALIGENMUONCOCKTAILPP_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+#include "AliGenCocktail.h"
+
+class AliGenCocktailEntry;
+
+class AliGenMUONCocktailpp : public AliGenCocktail
+{
+ public:
+
+    AliGenMUONCocktailpp();
+    AliGenMUONCocktailpp(const AliGenMUONCocktailpp &cocktail); 
+    virtual ~AliGenMUONCocktailpp();    
+    virtual void Init();
+    virtual void Generate();    
+    Int_t   GetNSucceded()         const {return fNSucceded;}    
+    Int_t   GetNGenerated()        const {return fNGenerated;}
+    Int_t   GetMuonMultiplicity()  const {return fMuonMultiplicity;}
+    Float_t GetMuonPtCut()         const {return fMuonPtCut;}
+    Float_t GetMuonThetaMin()      const {return fMuonThetaMinCut;}
+    Float_t GetMuonThetaMax()      const {return fMuonThetaMaxCut;}        
+    
+    void    SetMuonMultiplicity(Int_t MuonMultiplicity) { fMuonMultiplicity= MuonMultiplicity;}
+    void    SetMuonPtCut(Float_t PtCut) { fMuonPtCut = PtCut;}
+    void    SetMuonThetaRange(Float_t ThetaMin, Float_t ThetaMax){
+       fMuonThetaMinCut=ThetaMin;
+       fMuonThetaMaxCut=ThetaMax; }    
+
+ protected:
+
+    //
+ private:
+    Float_t fTotalRate;// Total rate of the full cocktail processes
+    Int_t   fMuonMultiplicity; // Muon multiplicity for the primordial trigger
+    Float_t fMuonPtCut;// Transverse momentum cut for muons
+    Float_t fMuonThetaMinCut;// Minimum theta cut for muons
+    Float_t fMuonThetaMaxCut; // Maximum theta cut for muons
+    Int_t   fNSucceded;// Number of Succes in the (di)-muon generation in the acceptance
+    Int_t   fNGenerated;// Number of generated cocktails
+
+    ClassDef(AliGenMUONCocktailpp,1)  //  cocktail for physics in the Alice
+};
+
+#endif
+
+
+
+
+
index 9fe8a8b..4b25369 100644 (file)
@@ -21,6 +21,7 @@
 #pragma link C++ class  AliGenParam+;
 #pragma link C++ class  AliGenCocktail+;
 #pragma link C++ class  AliGenMUONCocktail+;
+#pragma link C++ class  AliGenMUONCocktailpp+;
 #pragma link C++ class  AliGenCocktailAfterBurner+;
 #pragma link C++ class  AliGenCocktailEntry+;
 #pragma link C++ class  AliGenExtFile+;
index 96ca807..c3e21af 100644 (file)
@@ -16,7 +16,7 @@ SRCS          = AliGenHIJINGpara.cxx AliGenBox.cxx AliGenFixed.cxx \
                AliStructFuncType.cxx AliGenSlowNucleons.cxx \
                AliGenGeVSimEventHeader.cxx\
                AliSlowNucleonModel.cxx AliSlowNucleonModelExp.cxx \
-               AliGenMUONCocktail.cxx AliGenHBTosl.cxx AliGenCocktailEventHeader.cxx \
+               AliGenMUONCocktail.cxx AliGenMUONCocktailpp.cxx AliGenHBTosl.cxx AliGenCocktailEventHeader.cxx \
                AliGenReaderEMD.cxx
 
 # Headerfiles for this particular package (Path respect to own directory)
@@ -28,6 +28,6 @@ DHDR:=EVGENLinkDef.h
 
 EXPORT:=AliDecayer.h AliGenMC.h AliStructFuncType.h AliGenCocktailAfterBurner.h AliGenCocktail.h
 
-EINCLUDE:=FASTSIM
+EINCLUDE:=FASTSIM PYTHIA6