]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EVGEN/AliGenMUONCocktail.cxx
Warning fixed
[u/mrichter/AliRoot.git] / EVGEN / AliGenMUONCocktail.cxx
index 8e95ac2467caca93719db7932e1aee866c0efe3e..89fe4de562b6fd9a6fb432be1040a52418734fc7 100644 (file)
 // There is a primordial trigger which requires :
 //      a minimum muon multiplicity above a pT cut in a theta acceptance cone
 //
-// Gines Martinez, jan 2004, Nantes  martinez@in2p3.fr
-// Gines Martinez  sep 2004, Nantes  martinez@in2p3.fr
-//
-//
 
 #include <TObjArray.h>
 #include <TParticle.h>
@@ -51,43 +47,26 @@ ClassImp(AliGenMUONCocktail)
   
 //________________________________________________________________________
 AliGenMUONCocktail::AliGenMUONCocktail()
-  :AliGenCocktail()
+    :AliGenCocktail(),
+     fFastGlauber (0x0),
+     fTotalRate(0),  
+     fMuonMultiplicity(1),
+     fMuonPtCut(1.),
+     fMuonThetaMinCut(171.), 
+     fMuonThetaMaxCut(178.),
+     fNSucceded(0), 
+     fNGenerated(0), 
+     fLowImpactParameter(0.),
+     fHighImpactParameter(5.),
+     fAverageImpactParameter(0.),
+     fNumberOfCollisions(0.), 
+     fNumberOfParticipants(0.),
+     fHadronicMuons(kTRUE),
+     fInvMassCut (kFALSE),
+     fInvMassMinCut (0.),
+     fInvMassMaxCut (100.)
 {
 // Constructor
-  fFastGlauber =0x0;
-  fTotalRate =0;  
-  fNSucceded=0; 
-  fNGenerated=0; 
-  fMuonMultiplicity=1;
-  fMuonPtCut= 1.;
-  fMuonThetaMinCut=171.; 
-  fMuonThetaMaxCut=178.;
-  fLowImpactParameter = 0.;
-  fHighImpactParameter = 5.;
-  fAverageImpactParameter =0.;
-  fNumberOfCollisions   = 0.; 
-  fNumberOfParticipants = 0.;
-  fHadronicMuons = kTRUE;
-}
-//_________________________________________________________________________
-AliGenMUONCocktail::AliGenMUONCocktail(const AliGenMUONCocktail & cocktail):
-    AliGenCocktail(cocktail)
-{
-// Copy constructor
-  fFastGlauber =0x0;
-  fTotalRate =0; 
-  fNSucceded=0; 
-  fNGenerated=0; 
-  fMuonMultiplicity=1;
-  fMuonPtCut= 1.;
-  fMuonThetaMinCut=171.; 
-  fMuonThetaMaxCut=178.;  
-  fLowImpactParameter = 0.;
-  fHighImpactParameter = 5.;
-  fAverageImpactParameter =0.;
-  fNumberOfCollisions = 0.; 
-  fNumberOfParticipants = 0.;
-  fHadronicMuons = kTRUE;
 }
 //_________________________________________________________________________
 AliGenMUONCocktail::~AliGenMUONCocktail()
@@ -103,7 +82,7 @@ void AliGenMUONCocktail::Init()
   Double_t sigmaReaction =   0.072;   // MinBias NN cross section for PbPb LHC energies  http://arxiv.org/pdf/nucl-ex/0302016
 
   // Initialising Fast Glauber object
-  fFastGlauber = new AliFastGlauber();
+  fFastGlauber = AliFastGlauber::Instance();
   fFastGlauber->SetPbPbLHC();
   fFastGlauber->SetNNCrossSection(sigmaReaction*1000.); //Expected NN cross-section in mb at LHC with diffractive part http://arxiv.org/pdf/nucl-ex/0302016 )
   fFastGlauber->Init(1);
@@ -164,8 +143,8 @@ void AliGenMUONCocktail::Init()
   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));
 
   // Generating J/Psi Physics 
-  // Using Ramona Vogt distribution form CERN Yelow report and Andreas MORSCH private communication 
-  AliGenParam * genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "Vogt", "Jpsi");
+  // Using CFD scaled distribution (see http://clrwww.in2p3.fr/DIMUON2004/talks/sgrigoryan.pdf )
+  AliGenParam * genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF scaled", "Jpsi");
   genjpsi->SetPtRange(0,100.); // 4pi generation
   genjpsi->SetYRange(-8.,8);
   genjpsi->SetPhiRange(0.,360.);
@@ -189,8 +168,8 @@ void AliGenMUONCocktail::Init()
   fTotalRate+=ratiojpsi;
 
  // Generating Psi prime Physics
- // Using Ramona Vogt distribution form CERN Yelow report and Andreas MORSCH private communication 
-  AliGenParam * genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "Vogt", "PsiP");
+ // Using CFD scaled distribution (see http://clrwww.in2p3.fr/DIMUON2004/talks/sgrigoryan.pdf )
+  AliGenParam * genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF scaled", "PsiP");
   genpsiP->SetPtRange(0,100.);// 4pi generation
   genpsiP->SetYRange(-8.,8);
   genpsiP->SetPhiRange(0.,360.);
@@ -214,8 +193,8 @@ void AliGenMUONCocktail::Init()
   fTotalRate+=ratiopsiP;
 
   // Generating Upsilon Physics
-  // Using Ramona Vogt distribution form CERN Yelow report and Andreas MORSCH private communication 
-  AliGenParam * genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "Vogt", "Upsilon");  
+  // Using CFD scaled distribution (see http://clrwww.in2p3.fr/DIMUON2004/talks/sgrigoryan.pdf )
+  AliGenParam * genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "CDF scaled", "Upsilon");  
   genupsilon->SetPtRange(0,100.);  
   genupsilon->SetYRange(-8.,8);
   genupsilon->SetPhiRange(0.,360.);
@@ -236,8 +215,8 @@ void AliGenMUONCocktail::Init()
   fTotalRate+=ratioupsilon;
 
   // Generating UpsilonP Physics
-  // Using Ramona Vogt distribution form CERN Yelow report and Andreas MORSCH private communication 
-  AliGenParam * genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "Vogt", "UpsilonP");  
+  // Using CFD scaled distribution (see http://clrwww.in2p3.fr/DIMUON2004/talks/sgrigoryan.pdf )
+  AliGenParam * genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "CDF Scaled", "UpsilonP");  
   genupsilonP->SetPtRange(0,100.);  
   genupsilonP->SetYRange(-8.,8);
   genupsilonP->SetPhiRange(0.,360.);
@@ -258,8 +237,8 @@ void AliGenMUONCocktail::Init()
   fTotalRate+=ratioupsilonP;
 
   // Generating UpsilonPP Physics
-  // Using Ramona Vogt distribution form CERN Yelow report and Andreas MORSCH private communication 
-  AliGenParam * genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "Vogt", "UpsilonPP");  
+  // Using CFD scaled distribution (see http://clrwww.in2p3.fr/DIMUON2004/talks/sgrigoryan.pdf )
+  AliGenParam * genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF Scaled", "UpsilonPP");  
   genupsilonPP->SetPtRange(0,100.);  
   genupsilonPP->SetYRange(-8.,8);
   genupsilonPP->SetPhiRange(0.,360.);
@@ -280,7 +259,7 @@ void AliGenMUONCocktail::Init()
   fTotalRate+=ratioupsilonPP;
 
 // Generating non-correlated Charm Physics 
-  AliGenParam * gencharm = new AliGenParam(1, AliGenMUONlib::kCharm, "Vogt", "Charm");  
+  AliGenParam * gencharm = new AliGenParam(1, AliGenMUONlib::kCharm, "pp", "Charm");  
   gencharm->SetPtRange(0,100.);  
   gencharm->SetYRange(-8.,8);
   gencharm->SetPhiRange(0.,360.);
@@ -301,7 +280,7 @@ void AliGenMUONCocktail::Init()
   fTotalRate+=ratiocharm;
 
 // Generating non-correlated Beauty Physics 
-  AliGenParam * genbeauty = new AliGenParam(1, AliGenMUONlib::kBeauty, "Vogt", "Beauty");  
+  AliGenParam * genbeauty = new AliGenParam(1, AliGenMUONlib::kBeauty, "pp", "Beauty");  
   genbeauty->SetPtRange(0,100.);  
   genbeauty->SetYRange(-8.,8);
   genbeauty->SetPhiRange(0.,360.);
@@ -334,7 +313,7 @@ void AliGenMUONCocktail::Init()
     //    9 - 12       188    71    0.655
     //   12 - 16        15    10    0.086
     // We found the hadronic muons scales quite well with the number of participants  
-    AliGenParam * genpion = new AliGenParam(1, AliGenMUONlib::kPion, "Vogt", "Pion");  
+    AliGenParam * genpion = new AliGenParam(1, AliGenMUONlib::kPion, "default", "Pion");  
     genpion->SetPtRange(0,100.);  
     genpion->SetYRange(-8.,8);
     genpion->SetPhiRange(0.,360.);
@@ -355,7 +334,7 @@ void AliGenMUONCocktail::Init()
     fTotalRate+=ratiopion;
     
     // Generating Kaon Physics
-    AliGenParam * genkaon = new AliGenParam(1, AliGenMUONlib::kKaon, "Vogt", "Kaon");  
+    AliGenParam * genkaon = new AliGenParam(1, AliGenMUONlib::kKaon, "default", "Kaon");  
     genkaon->SetPtRange(0,100.);  
     genkaon->SetYRange(-8.,8);
     genkaon->SetPhiRange(0.,360.);
@@ -386,7 +365,7 @@ void AliGenMUONCocktail::Generate()
     AliGenCocktailEntry *entry = 0;
     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();
@@ -394,50 +373,80 @@ void AliGenMUONCocktail::Generate()
     // Loop on primordialTrigger
     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;
-      while((entry = (AliGenCocktailEntry*)next())) {
-       gen = entry->Generator();
-       gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
-       if ( (npart = gRandom->Poisson(entry->Rate())) >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();
-      // Tesitng primordial trigger : Muon  pair in the MUON spectrometer acceptance and pTCut
-      Int_t iPart;
-      fNGenerated++;
-      Int_t numberOfMuons=0;
-      //      printf(">>>fNGenerated is %d\n",fNGenerated);
-      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)                             ) { 
-         gAlice->GetMCApp()->Particle(iPart)->SetProductionVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2), 0.);   
-         numberOfMuons++;
+               //Reseting stack
+               AliRunLoader * runloader = AliRunLoader::Instance();
+               if (runloader)
+                       if (runloader->Stack())
+                               runloader->Stack()->Clean();
+               // Loop over generators and generate events
+               Int_t igen=0;
+               Int_t npart =0;
+               while((entry = (AliGenCocktailEntry*)next())) {
+                       gen = entry->Generator();
+                       gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
+                       if ( (npart = gRandom->Poisson(entry->Rate())) >0 ) {
+                               igen++; 
+                               if (igen == 1) entry->SetFirst(0);
+                               else  entry->SetFirst((partArray->GetEntriesFast())+1);
+                               // if ( (fHadronicMuons == kFALSE) && ( (gen->GetName() == "Pions") || (gen->GetName() == "Kaons") ) )
+                               //  { AliInfo(Form("This generator %s is finally not generated. This is option for hadronic muons.",gen->GetName() ) ); }
+                               // else {
+                               gen->SetNumberParticles(npart);
+                               gen->Generate();
+                               entry->SetLast(partArray->GetEntriesFast());
+                               preventry = entry;
+                         // }
+                       }
+               }  
+               next.Reset();
+               // Testing primordial trigger : Muon  pair in the MUON spectrometer acceptance and pTCut
+               Int_t iPart;
+               fNGenerated++;
+               Int_t numberOfMuons=0;
+               //      printf(">>>fNGenerated is %d\n",fNGenerated);
+               
+               TObjArray GoodMuons; // Used in the Invariant Mass selection cut
+               
+               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)                             ) { 
+                                       gAlice->GetMCApp()->Particle(iPart)->SetProductionVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2), 0.);   
+                                       GoodMuons.AddLast(gAlice->GetMCApp()->Particle(iPart));
+                                       numberOfMuons++;
+                       }                       
+               }
+               
+               // Test the invariant mass of each pair (if cut on Invariant mass is required)
+               Bool_t InvMassRangeOK = kTRUE;
+               if(fInvMassCut && (numberOfMuons>=2) ){
+                       TLorentzVector fV1, fV2, fVtot;
+                       InvMassRangeOK = kFALSE;
+                       for(iPart=0; iPart<GoodMuons.GetEntriesFast(); iPart++){      
+                       TParticle * mu1 = ((TParticle *)GoodMuons.At(iPart));
+
+                               for(int iPart2=iPart+1; iPart2<GoodMuons.GetEntriesFast(); iPart2++){      
+                                       TParticle * mu2 = ((TParticle *)GoodMuons.At(iPart2));
+                                               
+                                               fV1.SetPxPyPzE(mu1->Px() ,mu1->Py() ,mu1->Pz() ,mu1->Energy() );
+                                               fV2.SetPxPyPzE(mu2->Px() ,mu2->Py() ,mu2->Pz() ,mu2->Energy() );
+                                               fVtot = fV1 + fV2;
+                                               
+                                               if(fVtot.M()>fInvMassMinCut && fVtot.M()<fInvMassMaxCut) {
+                                                       InvMassRangeOK = kTRUE;
+                                                       break;
+                                               }
+                               }
+                               if(InvMassRangeOK) break; // Invariant Mass Cut pass as soon as one pair satisfy the criterion
+                       }       
+               }
+               
+               
+               if ((numberOfMuons >= fMuonMultiplicity) &&  InvMassRangeOK ) primordialTrigger = kTRUE;
        }
-      }
-      //  printf(">>> Number of Muons is %d \n", numberOfMuons);
-      if (numberOfMuons >= fMuonMultiplicity ) primordialTrigger = kTRUE;
-    }
     fNSucceded++;
-}
-
-
-
-
-
 
+    AliDebug(5,Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
+}