Introducing the interaction time into the aliroot generators. In case of gaussian...
[u/mrichter/AliRoot.git] / EVGEN / AliGenMUONCocktail.cxx
index e31db39..46dd2c8 100644 (file)
@@ -47,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()
@@ -99,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);
@@ -382,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();
@@ -390,56 +373,81 @@ 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);
-         // 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();
-      // 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));
+                       gen->SetTime(fTime);
+                       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), fTime);   
+                                       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));
 }
-
-
-
-
-
-