New options added for correlated dimuons. By A.Morsch.

diff --git a/EVGEN/DimuonCombinator.cxx b/EVGEN/DimuonCombinator.cxx
index e19d8e7f266196847eb14aeeb57882322fd68098..707bc1b09541db6e124ff78da2b29f3db37d5f32 100644 (file)
--- a/EVGEN/DimuonCombinator.cxx
+++ b/EVGEN/DimuonCombinator.cxx
@@ -159,7 +159,7 @@ void DimuonCombinator::SetSecondRange(Int_t from, Int_t to)
 //                       Selection
 //
 
-Int_t DimuonCombinator::Selected(GParticle* part)
+Bool_t DimuonCombinator::Selected(GParticle* part)
 {
 // 
 //
@@ -174,7 +174,7 @@ Int_t DimuonCombinator::Selected(GParticle* part)
     
 }
 
-Int_t DimuonCombinator::Selected(GParticle* part1, GParticle* part2)
+Bool_t DimuonCombinator::Selected(GParticle* part1, GParticle* part2)
 {
      return Selected(part1)*Selected(part2);
 }
@@ -225,17 +225,82 @@ void DimuonCombinator::SmearGauss(Float_t width, Float_t & value)
 //              Weighting
 // 
 
-Float_t DimuonCombinator::Weight(GParticle* part1, GParticle* part2)
+Float_t DimuonCombinator::Decay_Prob(GParticle* part)
 {
-    if (part1->GetParent() == part2->GetParent()) {
-       return (part1->GetWgt())*fRate1;
+    Float_t d, h, theta, CTau;
+    GParticle* parent = Parent(part);
+    Int_t ipar=Type(parent);
+    if (ipar==8 || ipar==9) {
+       CTau=780.4;
+    } else if (ipar==11 || ipar==12) {
+       CTau=370.9;
+    } else {
+       CTau=0;
+    }
+    
+    
+    Float_t GammaBeta=(parent->GetMomentum())/(parent->GetMass());
+//
+// this part is still very ALICE muon-arm specific
+//
+    theta=parent->GetTheta();
+    h=90*TMath::Tan(theta);
+    
+    if (h<4) {
+       d=4/TMath::Sin(theta);
     } else {
-       return (part1->GetWgt())*(part2->GetWgt())*fRate1*fRate2;
+       d=90/TMath::Cos(theta);
+    }
+    
+    if (CTau > 0) {
+       return 1-TMath::Exp(-d/CTau/GammaBeta);
+    } else {
+       return 1;
     }
 }
 
+Float_t DimuonCombinator::Weight(GParticle* part1, GParticle* part2)
+{
+    Float_t wgt=(part1->GetWgt())*(part2->GetWgt());
+    
+    if (Correlated(part1, part2)) {
+       return wgt/(Parent(part1)->GetWgt())*fRate1;
+    } else {
+       return wgt*fRate1*fRate2;
+    }
+} 
+
+
 Float_t DimuonCombinator::Weight(GParticle* part)
 {
-    return part->GetWgt()*fRate1;
+    return (part->GetWgt())*(Parent(part)->GetWgt())*fRate1;
+}
+Bool_t  DimuonCombinator::Correlated(GParticle* part1, GParticle* part2)
+{
+    if (Origin(part1) == Origin(part2)) {
+       return kTRUE;
+    } else {
+       return kFALSE;
+    }
+}
+GParticle* DimuonCombinator::Parent(GParticle* part)
+{
+    return (GParticle*) (fPartArray->UncheckedAt(part->GetParent()));
+}
+
+Int_t DimuonCombinator::Origin(GParticle* part)
+{
+    Int_t iparent= part->GetParent();
+    if (iparent < 0) return iparent;
+    Int_t ip;
+    while(1) {
+       ip=((GParticle*) fPartArray->UncheckedAt(iparent))->GetParent();
+       if (ip < 0) {
+           break;
+       } else {
+           iparent=ip;
+       }
+    }
+    return iparent;
 }
 

diff --git a/EVGEN/DimuonCombinator.h b/EVGEN/DimuonCombinator.h
index 30c31da6d722684987549475d9fcb413fafcd4c9..569353fdf5180e2d29202133407c57e01aef609a 100644 (file)
--- a/EVGEN/DimuonCombinator.h
+++ b/EVGEN/DimuonCombinator.h
@@ -30,24 +30,28 @@ public:
        fRate2=1.;
     }
 //    
-//  Iterators    
+//  Iterators
+//  Single muons
     GParticle* FirstMuon();
     GParticle* NextMuon();
+//  Single muons selected
     GParticle* FirstMuonSelected();
     GParticle* NextMuonSelected();
-    
+//  Dimuons    
     void FirstMuonPair(GParticle* & muon1, GParticle* & muon2);
     void NextMuonPair(GParticle* & muon1, GParticle* & muon2);
+//  Dimuons selected    
     void FirstMuonPairSelected(GParticle* & muon1, GParticle* & muon2);
     void NextMuonPairSelected(GParticle* & muon1, GParticle* & muon2);
+//  Loop over all prticles    
     void ResetRange();
+//  Set two ranges for dimuon loop    
     void SetFirstRange (Int_t from, Int_t to);
     void SetSecondRange(Int_t from, Int_t to);    
 //  Cuts
-    void SetPtMin(Float_t ptmin){fPtMin=ptmin;}
-    void SetEtaCut(Float_t etamin, Float_t etamax){fEtaMin=etamin; fEtaMax=etamax;}    
-    Int_t Selected(GParticle* part);
-    Int_t Selected(GParticle* part1, GParticle* part2);
+    void SetPtMin(Float_t ptmin) {fPtMin=ptmin;}
+    void SetEtaCut(Float_t etamin, Float_t etamax){fEtaMin=etamin; fEtaMax=etamax;}      Bool_t Selected(GParticle* part);
+    Bool_t Selected(GParticle* part1, GParticle* part2);
 // Kinematics
     Float_t Mass(GParticle* part1, GParticle* part);
     Float_t PT(GParticle* part1, GParticle* part);
@@ -55,19 +59,21 @@ public:
     Float_t Y(GParticle* part1, GParticle* part);
 // Response
     void SmearGauss(Float_t width, Float_t & value);
-    
 // Weight
+    Bool_t  Correlated(GParticle* part1, GParticle* part2);
     void    SetRate(Float_t rate){fRate1=rate;}
     void    SetRate(Float_t rate1, Float_t rate2 ){fRate1=rate1; fRate2=rate2;}
     Float_t Weight(GParticle* part);
     Float_t Weight(GParticle* part1, GParticle* part);
-
+    Float_t Decay_Prob(GParticle* part);
+    
  private:
     void FirstPartner();
     void NextPartner();
     void FirstPartnerSelected();
     void NextPartnerSelected();
-
+    Int_t Origin(GParticle* part);
+    GParticle* Parent(GParticle* part);
     GParticle* Partner();
     Int_t Type(GParticle *part){return part->GetKF();}
 private: