Add new global configuration for prod shuttle.
[u/mrichter/AliRoot.git] / MFT / AliMuonForwardTrackAnalysis.cxx
index 42dd32a..ad3cfd1 100644 (file)
@@ -13,6 +13,8 @@
 #include "TMatrixD.h"
 #include "TTree.h"
 #include "TH1D.h"
+#include "TH2D.h"
+#include "TH3D.h"
 #include "AliLog.h"
 #include "TFile.h"
 #include "TParticle.h"
@@ -24,6 +26,7 @@
 #include "TDatabasePDG.h"
 #include "TGraph.h"
 #include "AliMuonForwardTrackAnalysis.h"
+#include "TSystem.h"
 
 ClassImp(AliMuonForwardTrackAnalysis)
 
@@ -33,10 +36,19 @@ AliMuonForwardTrackAnalysis::AliMuonForwardTrackAnalysis():
   TObject(),
   fInputDir(0),
   fOutputDir(0),
-  fInputTree(0),
-  fMuonForwardTracks(0),
-  fMuonForwardTrackPairs(0),
+  fInputTreeWithBranson(0x0),
+  fInputTreeWithoutBranson(0x0),
+  fMuonForwardTracksWithBranson(0),
+  fMuonForwardTracksWithBransonMix(0),
+  fMuonForwardTrackPairsWithBranson(0),
+  fMuonForwardTracksWithoutBranson(0),
+  fMuonForwardTracksWithoutBransonMix(0),
+  fMuonForwardTrackPairsWithoutBranson(0),
+  fMFTTrackWithBranson(0),
+  fMFTTrackWithoutBranson(0),
   fMFTTrack(0),
+  fMFTTrackPairWithBranson(0),
+  fMFTTrackPairWithoutBranson(0),
   fMFTTrackPair(0),
   fMCRefTrack(0),
   fEv(0),
@@ -47,43 +59,73 @@ AliMuonForwardTrackAnalysis::AliMuonForwardTrackAnalysis():
   fNTracksAnalyzed(0),
   fNPairsOfEvent(0),
   fNPairsAnalyzedOfEvent(0),
-  fHistOffsetSingleMuonsX(0x0),
-  fHistOffsetSingleMuonsY(0x0),
-  fHistOffsetSingleMuons(0x0),
-  fHistWOffsetSingleMuons(0x0),
-  fHistErrorSingleMuonsX(0x0),
-  fHistErrorSingleMuonsY(0x0),
-  fHistOffsetSingleMuonsX_vsPtRapidity(0x0),
-  fHistOffsetSingleMuonsY_vsPtRapidity(0x0),
-  fHistSingleMuonsPtRapidity(0x0),
+  fNTracksAnalyzedOfEventAfterCut(0),
+  fNPairsAnalyzedOfEventAfterCut(0),
+  fHistXOffsetSingleMuonsVsEtaVsP(0x0),
+  fHistYOffsetSingleMuonsVsEtaVsP(0x0),
+  fHistOffsetSingleMuonsVsEtaVsP(0x0),
+  fHistWOffsetSingleMuonsVsEtaVsP(0x0),
+  fHistXOffsetSingleMuonsVsEtaVsPt(0x0),
+  fHistYOffsetSingleMuonsVsEtaVsPt(0x0),
+  fHistOffsetSingleMuonsVsEtaVsPt(0x0),
+  fHistWOffsetSingleMuonsVsEtaVsPt(0x0),
+  fHistXErrorSingleMuonsVsEtaVsP(0x0),
+  fHistYErrorSingleMuonsVsEtaVsP(0x0),
+  fHistXErrorSingleMuonsVsEtaVsPt(0x0),
+  fHistYErrorSingleMuonsVsEtaVsPt(0x0),
+  fHistZOriginSingleMuonsMC(0x0),
+  fHistZROriginSingleMuonsMC(0x0), 
+  fHistSingleMuonsPtRapidity(0x0), 
   fHistSingleMuonsOffsetChi2(0x0),
-  fHistWOffsetMuonPairs(0x0),
-  fHistMassMuonPairs(0x0),
-  fHistMassMuonPairsWithoutMFT(0x0),
-  fHistMassMuonPairsMC(0x0),
-  fHistRapidityPtMuonPairsMC(0x0),
-  fGraphSingleMuonsOffsetChi2(0x0),
-  fNMassBins(1000),
+  fHistMassMuonPairsMCVsPt(0x0),
+  fHistDimuonVtxResolutionXVsPt(0x0),
+  fHistDimuonVtxResolutionYVsPt(0x0),
+  fHistDimuonVtxResolutionZVsPt(0x0),
+  fTrueMass(0.),
   fMassMin(0),
-  fMassMax(10),
+  fMassMax(9999),
   fSingleMuonAnalysis(1),
   fMuonPairAnalysis(1),
-  fMatchTrigger(0),
   fOption(0),
+  fTriggerLevel(0),
   fXVertResMC(50.e-4),
   fYVertResMC(50.e-4),
   fZVertResMC(50.e-4),
+  fPrimaryVtxX(0.),
+  fPrimaryVtxY(0.),
+  fPrimaryVtxZ(0.),
   fMaxNWrongClustersMC(999),
-  fPtMinSingleMuons(0)
+  fMinPtSingleMuons(0),
+  fMinEtaSingleMuons(-99999.),
+  fMaxEtaSingleMuons(+99999.),
+  fUseBransonForCut(kFALSE),
+  fUseBransonForKinematics(kFALSE),
+  fCorrelateCutOnOffsetChi2(kFALSE),
+  fMaxChi2SingleMuons(1.e9), 
+  fMaxOffsetSingleMuons(1.e9),
+  fMaxWOffsetMuonPairsAtPrimaryVtx(1.e9), 
+  fMaxWOffsetMuonPairsAtPCA(1.e9), 
+  fMaxDistancePrimaryVtxPCA(1.e9), 
+  fMinPCAQuality(0.),
+  fMixing(kFALSE),
+  fNEventsToMix(10)
 {
 
   // default constructor
 
-  for (Int_t rapBin=0; rapBin<fNRapBinsOffsetSingleMuons; rapBin++) {
-    for (Int_t ptBin=0; ptBin<fNPtBinsOffsetSingleMuons; ptBin++) {
-      fHistOffsetSingleMuonsX_tmp[rapBin][ptBin] = 0x0;
-      fHistOffsetSingleMuonsY_tmp[rapBin][ptBin] = 0x0;
-    }
+  for (Int_t i=0; i<2; i++) {
+    fHistRapidityPtMuonPairs[i]              = 0;
+    fHistMassMuonPairsVsPt[i]                = 0;
+    fHistMassMuonPairsWithoutMFTVsPt[i]      = 0;
+    fHistMassMuonPairsVsPtLSp[i]             = 0;
+    fHistMassMuonPairsWithoutMFTVsPtLSp[i]   = 0;
+    fHistMassMuonPairsVsPtLSm[i]             = 0;
+    fHistMassMuonPairsWithoutMFTVsPtLSm[i]   = 0;
+    fHistWOffsetMuonPairsAtPrimaryVtxVsPt[i] = 0;
+    fHistWOffsetMuonPairsAtPCAVsPt[i]        = 0;
+    fHistDistancePrimaryVtxPCAVsPt[i]        = 0;
+    fHistPCAQualityVsPt[i]                   = 0;
+    fHistPseudoProperDecayLengthVsPt[i]      = 0;
   }
 
 }
@@ -93,67 +135,151 @@ AliMuonForwardTrackAnalysis::AliMuonForwardTrackAnalysis():
 Bool_t AliMuonForwardTrackAnalysis::Init(Char_t *inputFileName) {
 
   BookHistos();
-
+  
   TFile *inputFile = new TFile(Form("%s/%s",fInputDir.Data(),inputFileName));
   if (!inputFile || !inputFile->IsOpen()) {
     AliError(Form("Error opening file %s", inputFileName));
     return kFALSE;
   }
-  fInputTree = (TTree*) inputFile->Get("AliMuonForwardTracks");
-  if (!fInputTree) {
+  
+  fInputTreeWithBranson    = (TTree*) inputFile->Get("AliMuonForwardTracksWithBranson");
+  if (!fInputTreeWithBranson) {
     AliError("Error reading input tree");
     return kFALSE;
   }
-
-  if (fFirstEvent<0 || fLastEvent<0 || fFirstEvent>fLastEvent || fFirstEvent>=fInputTree->GetEntries()) {
+  
+  fInputTreeWithoutBranson    = (TTree*) inputFile->Get("AliMuonForwardTracksWithoutBranson");
+  if (!fInputTreeWithoutBranson) {
+    AliError("Error reading input tree");
+    return kFALSE;
+  }
+  
+  if (fFirstEvent>=fInputTreeWithBranson->GetEntriesFast()) return kFALSE;
+  else if (fFirstEvent<0 || fLastEvent<0 || fFirstEvent>fLastEvent || fFirstEvent>=fInputTreeWithBranson->GetEntriesFast()) {
     fFirstEvent = 0;
-    fLastEvent  = fInputTree->GetEntries()-1;
+    fLastEvent  = fInputTreeWithBranson->GetEntriesFast()-1;
   }
   else {
-    fLastEvent = TMath::Min(fLastEvent, Int_t(fInputTree->GetEntries()-1));
+    fLastEvent = TMath::Min(fLastEvent, Int_t(fInputTreeWithBranson->GetEntriesFast()-1));
   }
-
+  
   AliInfo(Form("Analysing events %d to %d", fFirstEvent, fLastEvent));
-
-  fMuonForwardTracks = new TClonesArray("AliMuonForwardTrack");
-  fInputTree->SetBranchAddress("tracks", &fMuonForwardTracks);  
-
+  
+  fMuonForwardTracksWithBranson = new TClonesArray("AliMuonForwardTrack",30);
+  fInputTreeWithBranson->SetBranchAddress("tracks", &fMuonForwardTracksWithBranson);  
+  
+  fMuonForwardTracksWithoutBranson = new TClonesArray("AliMuonForwardTrack",30);
+  fInputTreeWithoutBranson->SetBranchAddress("tracks", &fMuonForwardTracksWithoutBranson);  
+  
   TGeoManager::Import(Form("%s/geometry.root",fInputDir.Data()));
-
+  
   AliMUONTrackExtrap::SetField();
-
-  fMuonForwardTrackPairs = new TClonesArray("AliMuonForwardTrackPair");
-
+  
+  fMuonForwardTrackPairsWithBranson    = new TClonesArray("AliMuonForwardTrackPair",10);
+  fMuonForwardTrackPairsWithoutBranson = new TClonesArray("AliMuonForwardTrackPair",10);
+  fMuonForwardTrackPairsWithBranson    -> SetOwner(kTRUE);
+  fMuonForwardTrackPairsWithoutBranson -> SetOwner(kTRUE);
+  
   return kTRUE;
-
+  
 }
-
+  
 //====================================================================================================================================================
 
 Bool_t AliMuonForwardTrackAnalysis::LoadNextEvent() {
 
   if (fEv>fLastEvent) return kFALSE;
   if (fEv<fFirstEvent) { fEv++; return kTRUE; }
-  fMuonForwardTracks -> Clear();
-  fInputTree->GetEvent(fEv);
-  AliInfo(Form("**** analyzing event # %4d (%3d tracks) ****", fEv, fMuonForwardTracks->GetEntries()));
+
+  fMuonForwardTracksWithBranson    -> Clear("C");
+  fMuonForwardTracksWithoutBranson -> Clear("C");
+  fMuonForwardTracksWithBranson    -> Delete();
+  fMuonForwardTracksWithoutBranson -> Delete();
+  fInputTreeWithBranson    -> GetEvent(fEv);
+  fInputTreeWithoutBranson -> GetEvent(fEv);
+
+  AliDebug(2,Form("**** analyzing event # %4d (%3d tracks) ****", fEv, fMuonForwardTracksWithBranson->GetEntriesFast()));
+
+  AliInfo(Form("**** analyzing event # %6d of %6d ****", fEv, fLastEvent));
+
+  fPrimaryVtxX = gRandom->Gaus(0., fXVertResMC);
+  fPrimaryVtxY = gRandom->Gaus(0., fYVertResMC);
+  fPrimaryVtxZ = gRandom->Gaus(0., fZVertResMC);
 
   if (fSingleMuonAnalysis) {
     fNTracksAnalyzedOfEvent = 0;
-    fNTracksOfEvent = fMuonForwardTracks->GetEntries();
+    fNTracksAnalyzedOfEventAfterCut = 0;
+    fNTracksOfEvent = fMuonForwardTracksWithBranson->GetEntriesFast();
     while (AnalyzeSingleMuon()) continue;
   }
   
   if (fMuonPairAnalysis) {
-    if (fMuonForwardTrackPairs) {
-      fMuonForwardTrackPairs->Clear();
+    if (fMuonForwardTrackPairsWithBranson || fMuonForwardTrackPairsWithoutBranson) {
+      fMuonForwardTrackPairsWithBranson    -> Clear("C");
+      fMuonForwardTrackPairsWithoutBranson -> Clear("C");
+      fMuonForwardTrackPairsWithBranson    -> Delete();
+      fMuonForwardTrackPairsWithoutBranson -> Delete();
     }
     BuildMuonPairs();
     fNPairsAnalyzedOfEvent = 0;
-    fNPairsOfEvent = fMuonForwardTrackPairs->GetEntries();
-    while (AnalyzeMuonPair()) continue;
+    fNPairsAnalyzedOfEventAfterCut = 0;
+    fNPairsOfEvent = fMuonForwardTrackPairsWithBranson->GetEntriesFast();
+    while (AnalyzeMuonPair(kSingleEvent)) continue;
   }
 
+  AliDebug(2,Form("**** analyzed  event # %4d (%3d tracks and %3d pairs analyzed) ****", fEv, fNTracksAnalyzedOfEventAfterCut, fNPairsAnalyzedOfEventAfterCut));
+  
+  if (fMuonPairAnalysis && fMixing) {
+    for (fEvMix=fEv+1; fEvMix<=fEv+fNEventsToMix; fEvMix++) {
+      if (fEvMix<=fLastEvent) {
+
+       if (fMuonForwardTracksWithBransonMix || fMuonForwardTracksWithoutBransonMix) {
+         fMuonForwardTracksWithBransonMix    -> Delete();
+         fMuonForwardTracksWithoutBransonMix -> Delete();
+         delete fMuonForwardTracksWithBransonMix;
+         delete fMuonForwardTracksWithoutBransonMix;
+       }
+
+       fMuonForwardTracksWithBranson    -> Clear("C");
+       fMuonForwardTracksWithoutBranson -> Clear("C");
+       fMuonForwardTracksWithBranson    -> Delete();
+       fMuonForwardTracksWithoutBranson -> Delete();
+       fInputTreeWithBranson    -> GetEvent(fEvMix);
+       fInputTreeWithoutBranson -> GetEvent(fEvMix);
+
+       fMuonForwardTracksWithBransonMix    = new TClonesArray(*fMuonForwardTracksWithBranson);
+       fMuonForwardTracksWithoutBransonMix = new TClonesArray(*fMuonForwardTracksWithoutBranson);
+
+       fMuonForwardTracksWithBranson    -> Clear("C");
+       fMuonForwardTracksWithoutBranson -> Clear("C");
+       fMuonForwardTracksWithBranson    -> Delete();
+       fMuonForwardTracksWithoutBranson -> Delete();
+       fInputTreeWithBranson    -> GetEvent(fEv);
+       fInputTreeWithoutBranson -> GetEvent(fEv);
+
+       AliDebug(2,Form("**** mixing event # %4d (%3d tracks) with event # %4d (%3d tracks) ****", 
+                       fEv,    fMuonForwardTracksWithBranson->GetEntriesFast(),
+                       fEvMix, fMuonForwardTracksWithBransonMix ->GetEntriesFast()));    
+       if (fMuonForwardTrackPairsWithBranson || fMuonForwardTrackPairsWithoutBranson) {
+         fMuonForwardTrackPairsWithBranson    -> Clear("C");
+         fMuonForwardTrackPairsWithoutBranson -> Clear("C");
+         fMuonForwardTrackPairsWithBranson    -> Delete();
+         fMuonForwardTrackPairsWithoutBranson -> Delete();
+       }
+       BuildMuonPairsMix();
+       fNPairsAnalyzedOfEvent = 0;
+       fNPairsAnalyzedOfEventAfterCut = 0;
+       fNPairsOfEvent = fMuonForwardTrackPairsWithBranson->GetEntriesFast();
+       while (AnalyzeMuonPair(kMixedEvent)) continue;
+       AliDebug(2,Form("**** analyzed  mixed event pair (%4d ; %4d) : %3d pairs analyzed ****", fEv, fEvMix, fNPairsAnalyzedOfEventAfterCut));
+       
+       // delete fMuonForwardTracksWithBransonMix;
+       // delete fMuonForwardTracksWithoutBransonMix;
+
+      }
+    }
+  }
+    
   fEv++;
   
   return kTRUE;
@@ -166,61 +292,69 @@ Bool_t AliMuonForwardTrackAnalysis::AnalyzeSingleMuon() {
 
   if (fNTracksAnalyzedOfEvent>=fNTracksOfEvent) return kFALSE;
 
-  fMFTTrack = (AliMuonForwardTrack*) fMuonForwardTracks->At(fNTracksAnalyzedOfEvent);
+  fMFTTrackWithBranson    = (AliMuonForwardTrack*) fMuonForwardTracksWithBranson->At(fNTracksAnalyzedOfEvent);
+  fMFTTrackWithoutBranson = (AliMuonForwardTrack*) fMuonForwardTracksWithoutBranson->At(fNTracksAnalyzedOfEvent);
+
   fNTracksAnalyzedOfEvent++;
-  if (fMatchTrigger && !fMFTTrack->GetMatchTrigger()) return kTRUE;
-  fMCRefTrack = fMFTTrack->GetMCTrackRef();
 
-  if (!fMCRefTrack) return kTRUE;
-  if (fMFTTrack->GetNWrongClustersMC()>fMaxNWrongClustersMC) return kTRUE;
+  Bool_t passedCut = kFALSE;
+  if (fUseBransonForCut) passedCut = PassedCutSingleMuon(fMFTTrackWithBranson);
+  else passedCut = PassedCutSingleMuon(fMFTTrackWithoutBranson);
+  if (!passedCut) return kTRUE;
 
-  Double_t xOrig=gRandom->Gaus(0., fXVertResMC);
-  Double_t yOrig=gRandom->Gaus(0., fYVertResMC);
-//   Double_t xOrig = 0.;
-//   Double_t yOrig = 0.;
-  Double_t zOrig=gRandom->Gaus(0., fZVertResMC);
+  if (fUseBransonForKinematics) fMFTTrack = fMFTTrackWithBranson;
+  else fMFTTrack = fMFTTrackWithoutBranson;
+  if (fMFTTrack->GetMatchTrigger() < fTriggerLevel) return kTRUE;
+  fMCRefTrack = fMFTTrack->GetMCTrackRef();
 
-  AliMUONTrackParam *param = fMFTTrack->GetTrackParamAtMFTCluster(0);
-  AliMUONTrackExtrap::ExtrapToZCov(param, zOrig);
+  if (fMaxNWrongClustersMC<6) {
+    if (!fMCRefTrack) return kTRUE;
+    if (fMFTTrack->GetNWrongClustersMC()>fMaxNWrongClustersMC) return kTRUE;
+  }
 
-  TLorentzVector pMu;
-  Double_t mMu = TDatabasePDG::Instance()->GetParticle("mu-")->Mass();
-  Double_t energy = TMath::Sqrt(param->P()*param->P() + mMu*mMu);
-  pMu.SetPxPyPzE(param->Px(), param->Py(), param->Pz(), energy);
+  if (fMCRefTrack) {
+    fHistZOriginSingleMuonsMC  -> Fill(-1.*fMCRefTrack->Vz());
+    fHistZROriginSingleMuonsMC -> Fill(-1.*fMCRefTrack->Vz(), TMath::Sqrt(fMCRefTrack->Vx()*fMCRefTrack->Vx()+fMCRefTrack->Vy()*fMCRefTrack->Vy()));
+  }
 
-  if (fMFTTrack->Pt()<fPtMinSingleMuons) return kTRUE;
+  fMFTTrack -> EvalKinem(fPrimaryVtxZ);
 
   TMatrixD cov(5,5);
-  cov = param->GetCovariances();
+  cov = fMFTTrack->GetParamCovMatrix();
 
-  fHistErrorSingleMuonsX -> Fill(1.e4*TMath::Sqrt(cov(0,0)));
-  fHistErrorSingleMuonsY -> Fill(1.e4*TMath::Sqrt(cov(2,2)));
+  fHistXErrorSingleMuonsVsEtaVsP  -> Fill(1.e4*TMath::Sqrt(cov(0,0)), fMFTTrack->Eta(), fMFTTrack->P());
+  fHistYErrorSingleMuonsVsEtaVsP  -> Fill(1.e4*TMath::Sqrt(cov(2,2)), fMFTTrack->Eta(), fMFTTrack->P());
+  fHistXErrorSingleMuonsVsEtaVsPt -> Fill(1.e4*TMath::Sqrt(cov(0,0)), fMFTTrack->Eta(), fMFTTrack->Pt());
+  fHistYErrorSingleMuonsVsEtaVsPt -> Fill(1.e4*TMath::Sqrt(cov(2,2)), fMFTTrack->Eta(), fMFTTrack->Pt());
 
-  Double_t dX = fMFTTrack->GetOffsetX(xOrig, zOrig);
-  Double_t dY = fMFTTrack->GetOffsetY(yOrig, zOrig);
+  fHistXErrorSingleMuonsVsEtaVsP  -> Fill(1.e4*TMath::Sqrt(cov(0,0)), fMFTTrack->Eta(), fMFTTrack->P());
+  fHistYErrorSingleMuonsVsEtaVsP  -> Fill(1.e4*TMath::Sqrt(cov(2,2)), fMFTTrack->Eta(), fMFTTrack->P());
+  fHistXErrorSingleMuonsVsEtaVsPt -> Fill(1.e4*TMath::Sqrt(cov(0,0)), fMFTTrack->Eta(), fMFTTrack->Pt());
+  fHistYErrorSingleMuonsVsEtaVsPt -> Fill(1.e4*TMath::Sqrt(cov(2,2)), fMFTTrack->Eta(), fMFTTrack->Pt());
+
+  Double_t dX = fMFTTrack->GetOffsetX(fPrimaryVtxX, fPrimaryVtxZ);
+  Double_t dY = fMFTTrack->GetOffsetY(fPrimaryVtxY, fPrimaryVtxZ);
   
-  Double_t offset = fMFTTrack->GetOffset(xOrig, yOrig, zOrig);
-  Double_t weightedOffset = fMFTTrack->GetWeightedOffset(xOrig, yOrig, zOrig);
+  Double_t offset = fMFTTrack->GetOffset(fPrimaryVtxX, fPrimaryVtxY, fPrimaryVtxZ);
+  Double_t weightedOffset = fMFTTrack->GetWeightedOffset(fPrimaryVtxX, fPrimaryVtxY, fPrimaryVtxZ);
 
   //  AliDebug(2, Form("pdg code = %d\n", fMCRefTrack->GetPdgCode()));
 
-  fHistOffsetSingleMuonsX -> Fill(1.e4*dX);
-  fHistOffsetSingleMuonsY -> Fill(1.e4*dY);
-  Int_t rapBin = fHistOffsetSingleMuonsX_vsPtRapidity->GetXaxis()->FindBin(pMu.Rapidity());
-  Int_t ptBin  = fHistOffsetSingleMuonsX_vsPtRapidity->GetYaxis()->FindBin(pMu.Pt());
-  //  AliDebug(1, Form("pt = %f (%d), rap = %f (%d)\n", pMu.Pt(), pMu.Rapidity(), ptBin, rapBin));
-  if (0<rapBin && rapBin<=fNRapBinsOffsetSingleMuons && 0<ptBin && ptBin<=fNPtBinsOffsetSingleMuons) {
-    fHistOffsetSingleMuonsX_tmp[rapBin-1][ptBin-1]->Fill(1.e4*dX);
-    fHistOffsetSingleMuonsY_tmp[rapBin-1][ptBin-1]->Fill(1.e4*dY);
-  }
-  fHistSingleMuonsPtRapidity  -> Fill(pMu.Rapidity(), pMu.Pt());
-  fHistOffsetSingleMuons      -> Fill(1.e4*offset);
-  fHistWOffsetSingleMuons     -> Fill(weightedOffset);
+  fHistXOffsetSingleMuonsVsEtaVsP  -> Fill(1.e4*dX,        fMFTTrack->Eta(), fMFTTrack->P());
+  fHistYOffsetSingleMuonsVsEtaVsP  -> Fill(1.e4*dY,        fMFTTrack->Eta(), fMFTTrack->P());
+  fHistOffsetSingleMuonsVsEtaVsP   -> Fill(1.e4*offset,    fMFTTrack->Eta(), fMFTTrack->P());
+  fHistWOffsetSingleMuonsVsEtaVsP  -> Fill(weightedOffset, fMFTTrack->Eta(), fMFTTrack->P());
+  fHistXOffsetSingleMuonsVsEtaVsPt -> Fill(1.e4*dX,        fMFTTrack->Eta(), fMFTTrack->Pt());
+  fHistYOffsetSingleMuonsVsEtaVsPt -> Fill(1.e4*dY,        fMFTTrack->Eta(), fMFTTrack->Pt());
+  fHistOffsetSingleMuonsVsEtaVsPt  -> Fill(1.e4*offset,    fMFTTrack->Eta(), fMFTTrack->Pt());
+  fHistWOffsetSingleMuonsVsEtaVsPt -> Fill(weightedOffset, fMFTTrack->Eta(), fMFTTrack->Pt());
+
+  fHistSingleMuonsPtRapidity -> Fill(fMFTTrack->Rapidity(), fMFTTrack->Pt());
   Double_t chi2OverNdf = fMFTTrack->GetGlobalChi2()/Double_t(fMFTTrack->GetNMFTClusters()+fMFTTrack->GetNMUONClusters());
   fHistSingleMuonsOffsetChi2  -> Fill(1.e4*offset, chi2OverNdf);
-  fGraphSingleMuonsOffsetChi2 -> SetPoint(fGraphSingleMuonsOffsetChi2->GetN(),1.e4*offset, chi2OverNdf);
 
   fNTracksAnalyzed++;
+  fNTracksAnalyzedOfEventAfterCut++;
 
   return kTRUE;
 
@@ -228,31 +362,68 @@ Bool_t AliMuonForwardTrackAnalysis::AnalyzeSingleMuon() {
 
 //====================================================================================================================================================
 
-Bool_t AliMuonForwardTrackAnalysis::AnalyzeMuonPair() {
+Bool_t AliMuonForwardTrackAnalysis::AnalyzeMuonPair(Int_t opt) {
 
   if (fNPairsAnalyzedOfEvent>=fNPairsOfEvent) return kFALSE;
 
-  fMFTTrackPair = (AliMuonForwardTrackPair*) fMuonForwardTrackPairs->At(fNPairsAnalyzedOfEvent);
+  fMFTTrackPairWithBranson    = (AliMuonForwardTrackPair*) fMuonForwardTrackPairsWithBranson->At(fNPairsAnalyzedOfEvent);
+  fMFTTrackPairWithoutBranson = (AliMuonForwardTrackPair*) fMuonForwardTrackPairsWithoutBranson->At(fNPairsAnalyzedOfEvent);
 
-  if (fOption==kResonanceOnly && !fMFTTrackPair->IsResonance()) {
-    fNPairsAnalyzedOfEvent++;
-    return kTRUE;
-  }
+  fNPairsAnalyzedOfEvent++;
 
-  Double_t xOrig=gRandom->Gaus(0., fXVertResMC);
-  Double_t yOrig=gRandom->Gaus(0., fYVertResMC);
-  Double_t zOrig=gRandom->Gaus(0., fZVertResMC);
-  AliDebug(1, Form("origin = (%f, %f, %f)", xOrig, yOrig, zOrig));
+  if (fUseBransonForKinematics) fMFTTrackPair = fMFTTrackPairWithBranson;
+  else fMFTTrackPair = fMFTTrackPairWithoutBranson;
 
-  fHistMassMuonPairs           -> Fill(fMFTTrackPair->GetMass(zOrig));
-  fHistWOffsetMuonPairs        -> Fill(fMFTTrackPair->GetWeightedOffset(xOrig, yOrig, zOrig));
-  fHistMassMuonPairsWithoutMFT -> Fill(fMFTTrackPair->GetMassWithoutMFT(xOrig, yOrig, zOrig));
-  fHistMassMuonPairsMC         -> Fill(fMFTTrackPair->GetMassMC());
-  fHistRapidityPtMuonPairsMC   -> Fill(fMFTTrackPair->GetRapidityMC(), fMFTTrackPair->GetPtMC());
+  if ( fMFTTrackPair->GetTrack(0)->GetNWrongClustersMC()>fMaxNWrongClustersMC || 
+       fMFTTrackPair->GetTrack(1)->GetNWrongClustersMC()>fMaxNWrongClustersMC ) return kTRUE;
 
-  AliDebug(1, Form("mass = %f   MC = %f", fMFTTrackPair->GetMass(zOrig), fMFTTrackPair->GetMassMC()));
+  if (fOption==kResonanceOnly && !fMFTTrackPair->IsResonance())   return kTRUE;
+  if (fOption==kResonanceOnly && fMFTTrackPair->GetCharge() != 0) return kTRUE;
 
-  fNPairsAnalyzedOfEvent++;
+  Double_t pca[3] = {0};
+  fMFTTrackPair -> GetPointOfClosestApproach(pca);
+  Double_t distancePrimaryVtxPCA = TMath::Sqrt(TMath::Power(fPrimaryVtxX-pca[0],2)+
+                                              TMath::Power(fPrimaryVtxY-pca[1],2)+
+                                              TMath::Power(fPrimaryVtxZ-pca[2],2));
+
+  fMFTTrackPair -> SetKinem(fPrimaryVtxZ);
+
+  Bool_t passedCut = kFALSE;
+  if (fUseBransonForCut) passedCut = PassedCutMuonPair(fMFTTrackPairWithBranson);
+  else passedCut = PassedCutMuonPair(fMFTTrackPairWithoutBranson);
+
+  if (!passedCut) return kTRUE;
+
+  // --------------- Filling dimuon histograms --------------------------------
+
+  if (fMFTTrackPair->GetCharge() == 0) {
+    if (fOption==kResonanceOnly && opt==kSingleEvent) fHistMassMuonPairsMCVsPt -> Fill(fMFTTrackPair->GetMassMC(), fMFTTrackPair->GetPt());
+    fHistMassMuonPairsVsPt[opt]                -> Fill(fMFTTrackPair->GetMass(), fMFTTrackPair->GetPt());
+    fHistMassMuonPairsWithoutMFTVsPt[opt]      -> Fill(fMFTTrackPair->GetMassWithoutMFT(fPrimaryVtxX, fPrimaryVtxY, fPrimaryVtxZ), fMFTTrackPair->GetPt());
+    fHistWOffsetMuonPairsAtPrimaryVtxVsPt[opt] -> Fill(fMFTTrackPair->GetWeightedOffset(fPrimaryVtxX, fPrimaryVtxY, fPrimaryVtxZ), fMFTTrackPair->GetPt());
+    fHistWOffsetMuonPairsAtPCAVsPt[opt]        -> Fill(fMFTTrackPair->GetWeightedOffsetAtPCA(), fMFTTrackPair->GetPt());
+    fHistDistancePrimaryVtxPCAVsPt[opt]        -> Fill(distancePrimaryVtxPCA*1.e4, fMFTTrackPair->GetPt());
+    fHistPCAQualityVsPt[opt]                   -> Fill(fMFTTrackPair->GetPCAQuality(), fMFTTrackPair->GetPt());
+    fHistPseudoProperDecayLengthVsPt[opt]      -> Fill(GetPseudoProperDecayLength(fMFTTrackPair, fTrueMass), fMFTTrackPair->GetPt());
+    if (fEvalDimuonVtxResolution && opt==kSingleEvent) {
+      fHistDimuonVtxResolutionXVsPt->Fill(pca[0]*1.e4, fMFTTrackPair->GetPt());
+      fHistDimuonVtxResolutionYVsPt->Fill(pca[1]*1.e4, fMFTTrackPair->GetPt());
+      fHistDimuonVtxResolutionZVsPt->Fill(pca[2]*1.e4, fMFTTrackPair->GetPt());
+    }
+    fHistRapidityPtMuonPairs[opt] -> Fill(fMFTTrackPair->GetRapidity(), fMFTTrackPair->GetPt());
+  } 
+  else if (fMFTTrackPair->GetCharge() == -2) {
+    fHistMassMuonPairsVsPtLSm[opt]           -> Fill(fMFTTrackPair->GetMass(), fMFTTrackPair->GetPt());
+    fHistMassMuonPairsWithoutMFTVsPtLSm[opt] -> Fill(fMFTTrackPair->GetMassWithoutMFT(fPrimaryVtxX, fPrimaryVtxY, fPrimaryVtxZ), fMFTTrackPair->GetPt());
+  } 
+  else if (fMFTTrackPair->GetCharge() == 2) {
+    fHistMassMuonPairsVsPtLSp[opt]           -> Fill(fMFTTrackPair->GetMass(), fMFTTrackPair->GetPt());
+    fHistMassMuonPairsWithoutMFTVsPtLSp[opt] -> Fill(fMFTTrackPair->GetMassWithoutMFT(fPrimaryVtxX, fPrimaryVtxY, fPrimaryVtxZ), fMFTTrackPair->GetPt());
+  }
+  
+  AliDebug(1, Form("mass = %f   MC = %f", fMFTTrackPair->GetMass(), fMFTTrackPair->GetMassMC()));
+
+  fNPairsAnalyzedOfEventAfterCut++;
 
   return kTRUE;
 
@@ -262,24 +433,44 @@ Bool_t AliMuonForwardTrackAnalysis::AnalyzeMuonPair() {
 
 void AliMuonForwardTrackAnalysis::BuildMuonPairs() {
 
-  for (Int_t iTrack=0; iTrack<fMuonForwardTracks->GetEntries(); iTrack++) {
+  Int_t nMuonPairs = 0;
+
+  for (Int_t iTrack=0; iTrack<fMuonForwardTracksWithBranson->GetEntriesFast(); iTrack++) {
     for (Int_t jTrack=0; jTrack<iTrack; jTrack++) {
-    
-      AliMuonForwardTrack *track0 = (AliMuonForwardTrack*) fMuonForwardTracks->At(iTrack);
-      AliMuonForwardTrack *track1 = (AliMuonForwardTrack*) fMuonForwardTracks->At(jTrack);
-
-      if (fMatchTrigger) if (!track0->GetMatchTrigger() || !track1->GetMatchTrigger()) continue;
-      if (!track0->GetMCTrackRef() || !track1->GetMCTrackRef()) continue;
-      if (track0->GetNWrongClustersMC()>fMaxNWrongClustersMC || track1->GetNWrongClustersMC()>fMaxNWrongClustersMC) continue;
-      if (track0->Pt()<fPtMinSingleMuons || track1->Pt()<fPtMinSingleMuons) continue;
-
-      AliMuonForwardTrackPair *trackPair = new AliMuonForwardTrackPair(track0, track1);
-      if (fOption==kResonanceOnly && !trackPair->IsResonance()) {
-       delete trackPair;
-       continue;
-      }
-      new ((*fMuonForwardTrackPairs)[fMuonForwardTrackPairs->GetEntries()]) AliMuonForwardTrackPair(*trackPair);
+      
+      AliMuonForwardTrack *track0_WithBranson = (AliMuonForwardTrack*) fMuonForwardTracksWithBranson->At(iTrack);
+      AliMuonForwardTrack *track1_WithBranson = (AliMuonForwardTrack*) fMuonForwardTracksWithBranson->At(jTrack);
+      
+      AliMuonForwardTrack *track0_WithoutBranson = (AliMuonForwardTrack*) fMuonForwardTracksWithoutBranson->At(iTrack);
+      AliMuonForwardTrack *track1_WithoutBranson = (AliMuonForwardTrack*) fMuonForwardTracksWithoutBranson->At(jTrack);
+
+      if (track0_WithBranson->GetMatchTrigger()<fTriggerLevel || track1_WithBranson->GetMatchTrigger()<fTriggerLevel) continue;
+
+      // Resonance   = both mu from the same resonance
+      // Open Charm  = both mu from charmed mesons
+      // Open Beauty = both mu from beauty mesons
+      // Background1mu  = at least one mu from background
+      // Background2mu  = both mu from background
+
+      AliMuonForwardTrackPair *trackPairWithBranson = new AliMuonForwardTrackPair(track0_WithBranson, track1_WithBranson);
 
+      if (fOption==kResonanceOnly && !trackPairWithBranson->IsResonance()) { delete trackPairWithBranson; continue; }
+      if (fOption==kNoResonances  &&  trackPairWithBranson->IsResonance()) { delete trackPairWithBranson; continue; }
+      
+      if (fOption==kCharmOnly  && (!track0_WithBranson->IsFromCharm()  || !track1_WithBranson->IsFromCharm()))  { delete trackPairWithBranson; continue; }
+
+      if (fOption==kBeautyOnly && (!track0_WithBranson->IsFromBeauty() || !track1_WithBranson->IsFromBeauty())) { delete trackPairWithBranson; continue; }
+
+      if (fOption==kBackground1mu && (!track0_WithBranson->IsFromBackground() && !track1_WithBranson->IsFromBackground())) { delete trackPairWithBranson; continue; }
+      if (fOption==kBackground2mu && (!track0_WithBranson->IsFromBackground() || !track1_WithBranson->IsFromBackground())) { delete trackPairWithBranson; continue; }
+
+      delete trackPairWithBranson;
+      
+      new ((*fMuonForwardTrackPairsWithBranson)[nMuonPairs]) AliMuonForwardTrackPair(track0_WithBranson, track1_WithBranson);
+      new ((*fMuonForwardTrackPairsWithoutBranson)[nMuonPairs]) AliMuonForwardTrackPair(track0_WithoutBranson, track1_WithoutBranson);
+
+      nMuonPairs++;
+      
     }
   }
 
@@ -287,143 +478,387 @@ void AliMuonForwardTrackAnalysis::BuildMuonPairs() {
 
 //====================================================================================================================================================
 
-void AliMuonForwardTrackAnalysis::Terminate(Char_t *outputFileName) {
+void AliMuonForwardTrackAnalysis::BuildMuonPairsMix() {
 
-  for (Int_t rapBin=0; rapBin<fNRapBinsOffsetSingleMuons; rapBin++) {
-    for (Int_t ptBin=0; ptBin<fNPtBinsOffsetSingleMuons; ptBin++) {
-      Int_t binMin_x = fHistOffsetSingleMuonsX_tmp[rapBin][ptBin]->FindBin(-3*fHistOffsetSingleMuonsX_tmp[rapBin][ptBin]->GetRMS());
-      Int_t binMax_x = fHistOffsetSingleMuonsX_tmp[rapBin][ptBin]->FindBin(+3*fHistOffsetSingleMuonsX_tmp[rapBin][ptBin]->GetRMS());
-      for (Int_t bin=1; bin<=fHistOffsetSingleMuonsX_tmp[rapBin][ptBin]->GetNbinsX(); bin++) {
-       if (bin<binMin_x || bin>binMax_x) fHistOffsetSingleMuonsX_tmp[rapBin][ptBin]->SetBinContent(bin, 0);
-      }
-      Int_t binMin_y = fHistOffsetSingleMuonsY_tmp[rapBin][ptBin]->FindBin(-3*fHistOffsetSingleMuonsY_tmp[rapBin][ptBin]->GetRMS());
-      Int_t binMax_y = fHistOffsetSingleMuonsY_tmp[rapBin][ptBin]->FindBin(+3*fHistOffsetSingleMuonsY_tmp[rapBin][ptBin]->GetRMS());
-      for (Int_t bin=1; bin<=fHistOffsetSingleMuonsY_tmp[rapBin][ptBin]->GetNbinsX(); bin++) {
-       if (bin<binMin_y || bin>binMax_y) fHistOffsetSingleMuonsY_tmp[rapBin][ptBin]->SetBinContent(bin, 0);
-      }
-      fHistOffsetSingleMuonsX_vsPtRapidity->SetBinContent(rapBin+1, ptBin+1, fHistOffsetSingleMuonsX_tmp[rapBin][ptBin]->GetRMS());
-      fHistOffsetSingleMuonsX_vsPtRapidity->SetBinError(rapBin+1, ptBin+1, fHistOffsetSingleMuonsX_tmp[rapBin][ptBin]->GetRMSError());
-      fHistOffsetSingleMuonsY_vsPtRapidity->SetBinContent(rapBin+1, ptBin+1, fHistOffsetSingleMuonsY_tmp[rapBin][ptBin]->GetRMS());
-      fHistOffsetSingleMuonsY_vsPtRapidity->SetBinError(rapBin+1, ptBin+1, fHistOffsetSingleMuonsY_tmp[rapBin][ptBin]->GetRMSError());
+  Int_t nMuonPairs = 0;
+
+  for (Int_t iTrack=0; iTrack<fMuonForwardTracksWithBranson->GetEntriesFast(); iTrack++) {
+    for (Int_t jTrack=0; jTrack<fMuonForwardTracksWithBransonMix->GetEntriesFast(); jTrack++) {
+      
+      AliMuonForwardTrack *track0_WithBranson = (AliMuonForwardTrack*) fMuonForwardTracksWithBranson -> At(iTrack);
+      AliMuonForwardTrack *track1_WithBranson = (AliMuonForwardTrack*) fMuonForwardTracksWithBransonMix  -> At(jTrack);
+      
+      AliMuonForwardTrack *track0_WithoutBranson = (AliMuonForwardTrack*) fMuonForwardTracksWithoutBranson -> At(iTrack);
+      AliMuonForwardTrack *track1_WithoutBranson = (AliMuonForwardTrack*) fMuonForwardTracksWithoutBransonMix  -> At(jTrack);
+
+      if (track0_WithBranson->GetMatchTrigger()<fTriggerLevel || track1_WithBranson->GetMatchTrigger()<fTriggerLevel) continue;
+
+      // Resonance   = both mu from the same resonance
+      // Open Charm  = both mu from charmed mesons
+      // Open Beauty = both mu from beauty mesons
+      // Background1mu  = at least one mu from background
+      // Background2mu  = both mu from background
+
+      AliMuonForwardTrackPair *trackPairWithBranson = new AliMuonForwardTrackPair(track0_WithBranson, track1_WithBranson);
+
+      if (fOption==kResonanceOnly && !trackPairWithBranson->IsResonance()) { delete trackPairWithBranson; continue; }
+      if (fOption==kNoResonances  &&  trackPairWithBranson->IsResonance()) { delete trackPairWithBranson; continue; }
+      
+      if (fOption==kCharmOnly  && (!track0_WithBranson->IsFromCharm()  || !track1_WithBranson->IsFromCharm()))  { delete trackPairWithBranson; continue; }
+
+      if (fOption==kBeautyOnly && (!track0_WithBranson->IsFromBeauty() || !track1_WithBranson->IsFromBeauty())) { delete trackPairWithBranson; continue; }
+
+      if (fOption==kBackground1mu && (!track0_WithBranson->IsFromBackground() && !track1_WithBranson->IsFromBackground())) { delete trackPairWithBranson; continue; }
+      if (fOption==kBackground2mu && (!track0_WithBranson->IsFromBackground() || !track1_WithBranson->IsFromBackground())) { delete trackPairWithBranson; continue; }
+
+      delete trackPairWithBranson;
+
+      new ((*fMuonForwardTrackPairsWithBranson)[nMuonPairs])    AliMuonForwardTrackPair(track0_WithBranson,    track1_WithBranson);
+      new ((*fMuonForwardTrackPairsWithoutBranson)[nMuonPairs]) AliMuonForwardTrackPair(track0_WithoutBranson, track1_WithoutBranson);
+      
+      nMuonPairs++;
+      
     }
   }
 
-  TFile *fileOut = new TFile(Form("%s/%s",fOutputDir.Data(),outputFileName), "recreate");
+}
 
-  printf("Writing output objects to file %s\n", fileOut->GetName());
+//====================================================================================================================================================
 
-  fHistOffsetSingleMuonsX -> Write();
-  fHistOffsetSingleMuonsY -> Write();
-  fHistOffsetSingleMuons  -> Write();
-  fHistWOffsetSingleMuons -> Write();
-  fHistErrorSingleMuonsX  -> Write();
-  fHistErrorSingleMuonsY  -> Write();
+Bool_t AliMuonForwardTrackAnalysis::PassedCutSingleMuon(AliMuonForwardTrack *track) {
 
-  fHistOffsetSingleMuonsX_vsPtRapidity -> Write();
-  fHistOffsetSingleMuonsY_vsPtRapidity -> Write();
+  track -> EvalKinem(fPrimaryVtxZ);
+  if (track->Pt()<fMinPtSingleMuons) return kFALSE;
+  if (track->Eta()<fMinEtaSingleMuons || track->Eta()>fMaxEtaSingleMuons) return kFALSE;
+  
+  Double_t offset = 1.e4*track->GetOffset(fPrimaryVtxX, fPrimaryVtxY, fPrimaryVtxZ);
+  Double_t chi2OverNdf = track->GetGlobalChi2() / Double_t(track->GetNMFTClusters()+track->GetNMUONClusters()); 
+  offset /= fMaxOffsetSingleMuons;
+  chi2OverNdf /= fMaxChi2SingleMuons;
 
-//   for (Int_t rapBin=0; rapBin<fNRapBinsOffsetSingleMuons; rapBin++) {
-//     for (Int_t ptBin=0; ptBin<fNPtBinsOffsetSingleMuons; ptBin++) {
-//       fHistOffsetSingleMuonsX_tmp[rapBin][ptBin] -> Write();
-//       fHistOffsetSingleMuonsY_tmp[rapBin][ptBin] -> Write();
-//     }
-//   }
+  if (fCorrelateCutOnOffsetChi2) {
+    if (TMath::Sqrt(offset*offset + chi2OverNdf*chi2OverNdf) > 1) return kFALSE;
+  }
+  else {
+    if (offset>1 || chi2OverNdf>1) return kFALSE;
+  }
 
-  fHistSingleMuonsPtRapidity -> Write();
-  fHistSingleMuonsOffsetChi2 -> Write();
+  return kTRUE;
 
-  fGraphSingleMuonsOffsetChi2 -> Write();
+}
 
-  fHistWOffsetMuonPairs        -> Write();
-  fHistMassMuonPairs          -> Write();
-  fHistMassMuonPairsWithoutMFT -> Write();
-  fHistMassMuonPairsMC         -> Write();
-  fHistRapidityPtMuonPairsMC   -> Write();
+//====================================================================================================================================================
 
-  fileOut -> Close();
+Bool_t AliMuonForwardTrackAnalysis::PassedCutMuonPair(AliMuonForwardTrackPair *pair) {
+
+  if (!PassedCutSingleMuon(pair->GetTrack(0)) || !PassedCutSingleMuon(pair->GetTrack(1))) return kFALSE;
+
+  if (pair->GetMass()>fMassMax || pair->GetMass()<fMassMin) return kFALSE;
+
+  if (pair->GetWeightedOffset(fPrimaryVtxX, fPrimaryVtxY, fPrimaryVtxZ) > fMaxWOffsetMuonPairsAtPrimaryVtx) return kFALSE;
+  if (pair->GetWeightedOffsetAtPCA() > fMaxWOffsetMuonPairsAtPCA) return kFALSE;
+
+  Double_t pca[3] = {0};
+  pair -> GetPointOfClosestApproach(pca);
+  Double_t distancePrimaryVtxPCA = TMath::Sqrt(TMath::Power(fPrimaryVtxX-pca[0],2)+
+                                              TMath::Power(fPrimaryVtxY-pca[1],2)+
+                                              TMath::Power(fPrimaryVtxZ-pca[2],2));
+
+  if (distancePrimaryVtxPCA > fMaxDistancePrimaryVtxPCA) return kFALSE;
+  if (pair->GetPCAQuality() < fMinPCAQuality) return kFALSE;
+
+  return kTRUE;
 
 }
 
 //====================================================================================================================================================
 
-void AliMuonForwardTrackAnalysis::BookHistos() {
+Double_t AliMuonForwardTrackAnalysis::GetPseudoProperDecayLength(AliMuonForwardTrackPair *pair, Double_t trueMass) {
 
-  fHistOffsetSingleMuonsX = new TH1D("fHistOffsetSingleMuonsX", "Offset for single muons along X",  200, -1000, 1000);
-  fHistOffsetSingleMuonsY = new TH1D("fHistOffsetSingleMuonsY", "Offset for single muons along Y",  200, -1000, 1000);
-  fHistErrorSingleMuonsX  = new TH1D("fHistErrorSingleMuonsX",  "Coordinate Error for single muons along X",  200, 0, 1000);
-  fHistErrorSingleMuonsY  = new TH1D("fHistErrorSingleMuonsY",  "Coordinate Error for single muons along Y",  200, 0, 1000);
-  fHistOffsetSingleMuons  = new TH1D("fHistOffsetSingleMuons",  "Offset for single muons",          200, 0, 2000);
-  fHistWOffsetSingleMuons = new TH1D("fHistWOffsetSingleMuons", "Weighted Offset for single muons", 300, 0, 15);  
-
-  fHistOffsetSingleMuonsX_vsPtRapidity = new TH2D("fHistOffsetSingleMuonsX_vsPtRapidity", "Offset for single muons along X", 
-                                                 10, -4, -2.5, 10, 0.5, 5.5);
-  fHistOffsetSingleMuonsY_vsPtRapidity = new TH2D("fHistOffsetSingleMuonsY_vsPtRapidity", "Offset for single muons along Y", 
-                                                 10, -4, -2.5, 10, 0.5, 5.5);
-
-  for (Int_t rapBin=0; rapBin<fNRapBinsOffsetSingleMuons; rapBin++) {
-    for (Int_t ptBin=0; ptBin<fNPtBinsOffsetSingleMuons; ptBin++) {
-      fHistOffsetSingleMuonsX_tmp[rapBin][ptBin] = new TH1D(Form("fHistOffsetSingleMuonsX_tmp_%02d_%02d",rapBin,ptBin), "",  1000, -1000, 1000);
-      fHistOffsetSingleMuonsY_tmp[rapBin][ptBin] = new TH1D(Form("fHistOffsetSingleMuonsY_tmp_%02d_%02d",rapBin,ptBin), "",  1000, -1000, 1000);
-    }
+  Double_t pca[3] = {0};
+  pair -> GetPointOfClosestApproach(pca);
+
+  TVector2 vecVertexPCA(pca[0]-fPrimaryVtxX, pca[1]-fPrimaryVtxY);
+  TVector2 dimuonPt(pair->GetPx(), pair->GetPy());
+  TVector2 dimuonPtUnit = dimuonPt.Unit();
+  
+  Double_t l_xy = vecVertexPCA*dimuonPtUnit;
+  Double_t pseudoProperDecayLength = (l_xy * trueMass / pair->GetPt()) * 1.e4 ; // in micron
+
+  return pseudoProperDecayLength;
+  
+}
+
+//====================================================================================================================================================
+
+void AliMuonForwardTrackAnalysis::Terminate(Char_t *outputFileName) {
+
+  TFile *fileOut = new TFile(Form("%s/%s",fOutputDir.Data(),outputFileName), "recreate");
+
+  printf("Writing output objects to file %s\n", fileOut->GetName());
+
+  // single muons
+
+  fHistXOffsetSingleMuonsVsEtaVsP  -> Write();
+  fHistYOffsetSingleMuonsVsEtaVsP  -> Write();
+  fHistXOffsetSingleMuonsVsEtaVsPt -> Write();
+  fHistYOffsetSingleMuonsVsEtaVsPt -> Write();
+
+  fHistXErrorSingleMuonsVsEtaVsP   -> Write();
+  fHistYErrorSingleMuonsVsEtaVsP   -> Write();
+  fHistXErrorSingleMuonsVsEtaVsPt  -> Write();
+  fHistYErrorSingleMuonsVsEtaVsPt  -> Write();
+
+  fHistOffsetSingleMuonsVsEtaVsP   -> Write();
+  fHistOffsetSingleMuonsVsEtaVsPt  -> Write();
+  fHistWOffsetSingleMuonsVsEtaVsP  -> Write();
+  fHistWOffsetSingleMuonsVsEtaVsPt -> Write();
+
+  fHistSingleMuonsPtRapidity  -> Write();
+  fHistSingleMuonsOffsetChi2  -> Write();
+  fHistZOriginSingleMuonsMC   -> Write();
+  fHistZROriginSingleMuonsMC  -> Write();
+
+  // dimuons
+
+  fHistMassMuonPairsVsPt[kSingleEvent]                -> Write();
+  fHistMassMuonPairsWithoutMFTVsPt[kSingleEvent]      -> Write();
+  fHistMassMuonPairsVsPtLSp[kSingleEvent]             -> Write();
+  fHistMassMuonPairsWithoutMFTVsPtLSp[kSingleEvent]   -> Write();
+  fHistMassMuonPairsVsPtLSm[kSingleEvent]             -> Write();
+  fHistMassMuonPairsWithoutMFTVsPtLSm[kSingleEvent]   -> Write();
+  fHistWOffsetMuonPairsAtPrimaryVtxVsPt[kSingleEvent] -> Write();
+  fHistWOffsetMuonPairsAtPCAVsPt[kSingleEvent]        -> Write();
+  fHistDistancePrimaryVtxPCAVsPt[kSingleEvent]        -> Write();
+  fHistPCAQualityVsPt[kSingleEvent]                   -> Write();
+  fHistPseudoProperDecayLengthVsPt[kSingleEvent]      -> Write();
+  fHistRapidityPtMuonPairs[kSingleEvent]              -> Write();
+  fHistMassMuonPairsMCVsPt                            -> Write();
+  if (fEvalDimuonVtxResolution) {
+    fHistDimuonVtxResolutionXVsPt -> Write();
+    fHistDimuonVtxResolutionYVsPt -> Write();
+    fHistDimuonVtxResolutionZVsPt -> Write();
   }
 
-  fHistSingleMuonsPtRapidity = new TH2D("fHistSingleMuonsPtRapidity", "Phase Space for single muons", 10, -4, -2.5, 10, 0.5, 5.5);
-  fHistSingleMuonsOffsetChi2 = new TH2D("fHistSingleMuonsOffsetChi2", "Offset vs #chi^{2}/ndf for single muons", 400, 0, 4000, 100, 0, 20);
+  if (fMixing) {
+    fHistMassMuonPairsVsPt[kMixedEvent]                -> Write();
+    fHistMassMuonPairsWithoutMFTVsPt[kMixedEvent]      -> Write();
+    fHistMassMuonPairsVsPtLSp[kMixedEvent]             -> Write();
+    fHistMassMuonPairsWithoutMFTVsPtLSp[kMixedEvent]   -> Write();
+    fHistMassMuonPairsVsPtLSm[kMixedEvent]             -> Write();
+    fHistMassMuonPairsWithoutMFTVsPtLSm[kMixedEvent]   -> Write();
+    fHistWOffsetMuonPairsAtPrimaryVtxVsPt[kMixedEvent] -> Write();
+    fHistWOffsetMuonPairsAtPCAVsPt[kMixedEvent]        -> Write();
+    fHistDistancePrimaryVtxPCAVsPt[kMixedEvent]        -> Write();
+    fHistPCAQualityVsPt[kMixedEvent]                   -> Write();
+    fHistPseudoProperDecayLengthVsPt[kMixedEvent]      -> Write();
+    fHistRapidityPtMuonPairs[kMixedEvent]              -> Write();
+  }
 
-  fHistOffsetSingleMuonsX -> SetXTitle("Offset(X)  [#mum]");
-  fHistOffsetSingleMuonsY -> SetXTitle("Offset(Y)  [#mum]");
-  fHistErrorSingleMuonsX  -> SetXTitle("Err. on track position at z_{vtx} (X)  [#mum]");
-  fHistErrorSingleMuonsY  -> SetXTitle("Err. on track position at z_{vtx} (Y)  [#mum]");
-  fHistOffsetSingleMuons  -> SetXTitle("Offset  [#mum]");
-  fHistWOffsetSingleMuons -> SetXTitle("Weighted Offset");
+  fileOut -> Close();
 
-  fHistOffsetSingleMuonsX_vsPtRapidity -> SetXTitle("y^{#mu}");
-  fHistOffsetSingleMuonsY_vsPtRapidity -> SetXTitle("y^{#mu}");
-  fHistOffsetSingleMuonsX_vsPtRapidity -> SetYTitle("p_{T}^{#mu}  [GeV/c]");
-  fHistOffsetSingleMuonsY_vsPtRapidity -> SetYTitle("p_{T}^{#mu}  [GeV/c]");
+}
+
+//====================================================================================================================================================
+
+void AliMuonForwardTrackAnalysis::BookHistos() {
+
+  // -------------------- single muons
+
+  fHistXOffsetSingleMuonsVsEtaVsP  = new TH3D("fHistXOffsetSingleMuonsVsEtaVsP",  "Offset for single muons along X vs #eta vs P",      200, -1000, 1000, 100, -4.5, -2., 100, 0, 100);
+  fHistYOffsetSingleMuonsVsEtaVsP  = new TH3D("fHistYOffsetSingleMuonsVsEtaVsP",  "Offset for single muons along Y vs #eta vs P",      200, -1000, 1000, 100, -4.5, -2., 100, 0, 100);
+  fHistXOffsetSingleMuonsVsEtaVsPt = new TH3D("fHistXOffsetSingleMuonsVsEtaVsPt", "Offset for single muons along X vs #eta vs p_{T}",  200, -1000, 1000, 100, -4.5, -2., 100, 0,  10);
+  fHistYOffsetSingleMuonsVsEtaVsPt = new TH3D("fHistYOffsetSingleMuonsVsEtaVsPt", "Offset for single muons along Y vs #eta vs p_{T}",  200, -1000, 1000, 100, -4.5, -2., 100, 0,  10);
+
+  fHistXErrorSingleMuonsVsEtaVsP   = new TH3D("fHistXErrorSingleMuonsVsEtaVsP",   "Coordinate Error for single muons along X vs #eta vs P",      200, 0, 1000, 100, -4.5, -2., 100, 0, 100);
+  fHistYErrorSingleMuonsVsEtaVsP   = new TH3D("fHistYErrorSingleMuonsVsEtaVsP",   "Coordinate Error for single muons along Y vs #eta vs P",      200, 0, 1000, 100, -4.5, -2., 100, 0, 100);
+  fHistXErrorSingleMuonsVsEtaVsPt  = new TH3D("fHistXErrorSingleMuonsVsEtaVsPt",  "Coordinate Error for single muons along X vs #eta vs p_{T}",  200, 0, 1000, 100, -4.5, -2., 100, 0,  10);
+  fHistYErrorSingleMuonsVsEtaVsPt  = new TH3D("fHistYErrorSingleMuonsVsEtaVsPt",  "Coordinate Error for single muons along Y vs #eta vs p_{T}",  200, 0, 1000, 100, -4.5, -2., 100, 0,  10);
+
+  fHistOffsetSingleMuonsVsEtaVsP   = new TH3D("fHistOffsetSingleMuonsVsEtaVsP",   "Offset for single muons vs #eta vs P",              200, 0, 2000, 100, -4.5, -2., 100, 0, 100);
+  fHistOffsetSingleMuonsVsEtaVsPt  = new TH3D("fHistOffsetSingleMuonsVsEtaVsPt",  "Offset for single muons vs #eta vs p_{T}",          200, 0, 2000, 100, -4.5, -2., 100, 0,  10);
+  fHistWOffsetSingleMuonsVsEtaVsP  = new TH3D("fHistWOffsetSingleMuonsVsEtaVsP",  "Weighted Offset for single muons vs #eta vs P",     300, 0,   15, 100, -4.5, -2., 100, 0, 100);  
+  fHistWOffsetSingleMuonsVsEtaVsPt = new TH3D("fHistWOffsetSingleMuonsVsEtaVsPt", "Weighted Offset for single muons vs #eta vs p_{T}", 300, 0,   15, 100, -4.5, -2., 100, 0,  10);  
+
+  fHistSingleMuonsPtRapidity = new TH2D("fHistSingleMuonsPtRapidity", "Phase Space for single muons", 100, -4.5, -2., 100, 0., 10.);
+  fHistSingleMuonsOffsetChi2 = new TH2D("fHistSingleMuonsOffsetChi2", "Offset vs #chi^{2}/ndf for single muons", 400, 0, 4000, 200, 0, 10);
+  fHistZOriginSingleMuonsMC  = new TH1D("fHistZOriginSingleMuonsMC",  "Z origin for single muons (from MC)",   1000, 0., 500.);
+  fHistZROriginSingleMuonsMC = new TH2D("fHistZROriginSingleMuonsMC", "Z-R origin for single muons (from MC)", 1000, 0., 500., 1000, 0., 100.);
+
+  //
+
+  fHistXOffsetSingleMuonsVsEtaVsP  -> SetXTitle("Offset(X)  [#mum]");
+  fHistYOffsetSingleMuonsVsEtaVsP  -> SetXTitle("Offset(Y)  [#mum]");
+  fHistXOffsetSingleMuonsVsEtaVsPt -> SetXTitle("Offset(X)  [#mum]");
+  fHistYOffsetSingleMuonsVsEtaVsPt -> SetXTitle("Offset(Y)  [#mum]");
+
+  fHistXErrorSingleMuonsVsEtaVsP   -> SetXTitle("Err. on track position at z_{vtx} (X)  [#mum]");
+  fHistYErrorSingleMuonsVsEtaVsP   -> SetXTitle("Err. on track position at z_{vtx} (Y)  [#mum]");
+  fHistXErrorSingleMuonsVsEtaVsPt  -> SetXTitle("Err. on track position at z_{vtx} (X)  [#mum]");
+  fHistYErrorSingleMuonsVsEtaVsPt  -> SetXTitle("Err. on track position at z_{vtx} (Y)  [#mum]");
+
+  fHistOffsetSingleMuonsVsEtaVsP   -> SetXTitle("Offset  [#mum]");
+  fHistOffsetSingleMuonsVsEtaVsPt  -> SetXTitle("Offset  [#mum]");
+  fHistWOffsetSingleMuonsVsEtaVsP  -> SetXTitle("Weighted Offset");
+  fHistWOffsetSingleMuonsVsEtaVsPt -> SetXTitle("Weighted Offset");    
+
+  fHistXOffsetSingleMuonsVsEtaVsP  -> SetYTitle("#eta");
+  fHistYOffsetSingleMuonsVsEtaVsP  -> SetYTitle("#eta");
+  fHistXOffsetSingleMuonsVsEtaVsPt -> SetYTitle("#eta");
+  fHistYOffsetSingleMuonsVsEtaVsPt -> SetYTitle("#eta");
+
+  fHistXErrorSingleMuonsVsEtaVsP   -> SetYTitle("#eta");
+  fHistYErrorSingleMuonsVsEtaVsP   -> SetYTitle("#eta");
+  fHistXErrorSingleMuonsVsEtaVsPt  -> SetYTitle("#eta");
+  fHistYErrorSingleMuonsVsEtaVsPt  -> SetYTitle("#eta");
+
+  fHistOffsetSingleMuonsVsEtaVsP   -> SetYTitle("#eta");
+  fHistOffsetSingleMuonsVsEtaVsPt  -> SetYTitle("#eta");
+  fHistWOffsetSingleMuonsVsEtaVsP  -> SetYTitle("#eta");
+  fHistWOffsetSingleMuonsVsEtaVsPt -> SetYTitle("#eta");    
+
+  fHistXOffsetSingleMuonsVsEtaVsP  -> SetZTitle("P  [GeV/c]");
+  fHistYOffsetSingleMuonsVsEtaVsP  -> SetZTitle("P  [GeV/c]");
+  fHistXOffsetSingleMuonsVsEtaVsPt -> SetZTitle("p_{T}  [GeV/c]");
+  fHistYOffsetSingleMuonsVsEtaVsPt -> SetZTitle("p_{T}  [GeV/c]");
+
+  fHistXErrorSingleMuonsVsEtaVsP   -> SetZTitle("P  [GeV/c]");
+  fHistYErrorSingleMuonsVsEtaVsP   -> SetZTitle("P  [GeV/c]");
+  fHistXErrorSingleMuonsVsEtaVsPt  -> SetZTitle("p_{T}  [GeV/c]");
+  fHistYErrorSingleMuonsVsEtaVsPt  -> SetZTitle("p_{T}  [GeV/c]");
+
+  fHistOffsetSingleMuonsVsEtaVsP   -> SetZTitle("P  [GeV/c]");
+  fHistOffsetSingleMuonsVsEtaVsPt  -> SetZTitle("p_{T}  [GeV/c]");
+  fHistWOffsetSingleMuonsVsEtaVsP  -> SetZTitle("P  [GeV/c]");
+  fHistWOffsetSingleMuonsVsEtaVsPt -> SetZTitle("p_{T}  [GeV/c]");    
 
   fHistSingleMuonsPtRapidity -> SetXTitle("y^{#mu}");
   fHistSingleMuonsPtRapidity -> SetYTitle("p_{T}^{#mu}  [GeV/c]");
   fHistSingleMuonsOffsetChi2 -> SetXTitle("Offset  [#mum]");
   fHistSingleMuonsOffsetChi2 -> SetYTitle("#chi^{2}/ndf");
 
-  fHistOffsetSingleMuonsX -> Sumw2();
-  fHistOffsetSingleMuonsY -> Sumw2();
-  fHistErrorSingleMuonsX  -> Sumw2();
-  fHistErrorSingleMuonsY  -> Sumw2();
-  fHistOffsetSingleMuons  -> Sumw2();
-  fHistWOffsetSingleMuons -> Sumw2();
-
-  fHistOffsetSingleMuonsX_vsPtRapidity -> Sumw2();
-  fHistOffsetSingleMuonsY_vsPtRapidity -> Sumw2();
-  fHistSingleMuonsPtRapidity           -> Sumw2();
-  fHistSingleMuonsOffsetChi2           -> Sumw2();
+  fHistZOriginSingleMuonsMC  -> SetXTitle("Z  [cm]");
+  fHistZROriginSingleMuonsMC -> SetXTitle("Z  [cm]");
+  fHistZROriginSingleMuonsMC -> SetYTitle("R  [cm]");
+
+  //
+
+  fHistXOffsetSingleMuonsVsEtaVsP  -> Sumw2();
+  fHistYOffsetSingleMuonsVsEtaVsP  -> Sumw2();
+  fHistXOffsetSingleMuonsVsEtaVsPt -> Sumw2();
+  fHistYOffsetSingleMuonsVsEtaVsPt -> Sumw2();
+
+  fHistXErrorSingleMuonsVsEtaVsP   -> Sumw2();
+  fHistYErrorSingleMuonsVsEtaVsP   -> Sumw2();
+  fHistXErrorSingleMuonsVsEtaVsPt  -> Sumw2();
+  fHistYErrorSingleMuonsVsEtaVsPt  -> Sumw2();
+
+  fHistOffsetSingleMuonsVsEtaVsP   -> Sumw2();
+  fHistOffsetSingleMuonsVsEtaVsPt  -> Sumw2();
+  fHistWOffsetSingleMuonsVsEtaVsP  -> Sumw2();
+  fHistWOffsetSingleMuonsVsEtaVsPt -> Sumw2();
+
+  fHistSingleMuonsPtRapidity  -> Sumw2();
+  fHistSingleMuonsOffsetChi2  -> Sumw2();
+  fHistZOriginSingleMuonsMC   -> Sumw2();
+  fHistZROriginSingleMuonsMC  -> Sumw2();
+
+  // -------------------- muon pairs
+
+  Int_t nBinsPt = 20; Double_t ptMin = 0., ptMax = 10.;   // dimuon pt
+
+  fHistMassMuonPairsVsPt[kSingleEvent]               = new TH2D("fHistMassMuonPairsVsPt", "Dimuon Mass (MUON+MFT) vs p_{T}^{#mu#mu}", 1000, 0, 10, nBinsPt, ptMin, ptMax);
+  fHistMassMuonPairsWithoutMFTVsPt[kSingleEvent]      = new TH2D("fHistMassMuonPairsWithoutMFTVsPt", "Dimuon Mass (MUON only) vs p_{T}^{#mu#mu}", 1000, 0, 10, nBinsPt, ptMin, ptMax); 
+  fHistMassMuonPairsVsPtLSp[kSingleEvent]            = new TH2D("fHistMassMuonPairsVsPtLSp", "Dimuon Mass (MUON+MFT) vs p_{T}^{#mu#mu} Like-Sign ++", 1000, 0, 10, nBinsPt, ptMin, ptMax);
+  fHistMassMuonPairsWithoutMFTVsPtLSp[kSingleEvent]   = new TH2D("fHistMassMuonPairsWithoutMFTVsPtLSp", "Dimuon Mass (MUON only) vs p_{T}^{#mu#mu} Like-Sign ++", 1000, 0, 10, nBinsPt, ptMin, ptMax); 
+  fHistMassMuonPairsVsPtLSm[kSingleEvent]            = new TH2D("fHistMassMuonPairsVsPtLSm", "Dimuon Mass (MUON+MFT) vs p_{T}^{#mu#mu} Like-Sign --", 1000, 0, 10, nBinsPt, ptMin, ptMax);
+  fHistMassMuonPairsWithoutMFTVsPtLSm[kSingleEvent]   = new TH2D("fHistMassMuonPairsWithoutMFTVsPtLSm", "Dimuon Mass (MUON only) vs p_{T}^{#mu#mu} Like-Sign --", 1000, 0, 10, nBinsPt, ptMin, ptMax); 
+  fHistWOffsetMuonPairsAtPrimaryVtxVsPt[kSingleEvent] = new TH2D("fHistWOffsetMuonPairsAtPrimaryVtxVsPt", "Weighted Offset for Muon Pairs at Primary Vertex vs p_{T}^{#mu#mu}", 300, 0, 60, nBinsPt, ptMin, ptMax);
+  fHistWOffsetMuonPairsAtPCAVsPt[kSingleEvent]        = new TH2D("fHistWOffsetMuonPairsAtPCAVsPt", "Weighted Offset for Muon Pairs at PCA vs p_{T}^{#mu#mu}", 300, 0, 60, nBinsPt, ptMin, ptMax);
+  fHistDistancePrimaryVtxPCAVsPt[kSingleEvent]        = new TH2D("fHistDistancePrimaryVtxPCAVsPt", "Distance between PCA and primary vertex vs p_{T}^{#mu#mu}", 1000, 0, 50000, nBinsPt, ptMin, ptMax);
+  fHistPCAQualityVsPt[kSingleEvent]                   = new TH2D("fHistPCAQualityVsPt", "PCA Quality vs p_{T}^{#mu#mu}", 200, 0, 1, nBinsPt, ptMin, ptMax);
+  fHistPseudoProperDecayLengthVsPt[kSingleEvent]      = new TH2D("fHistPseudoProperDecayLengthVsPt", "Pseudo proper decay length vs p_{T}^{#mu#mu}", 1000, -5000, 5000, nBinsPt, ptMin, ptMax);
+
+  fHistMassMuonPairsVsPt[kMixedEvent]                = new TH2D("fHistMassMuonPairsVsPtMix", "Dimuon Mass Mix (MUON+MFT) vs p_{T}^{#mu#mu}", 1000, 0, 10, nBinsPt, ptMin, ptMax);
+  fHistMassMuonPairsWithoutMFTVsPt[kMixedEvent]       = new TH2D("fHistMassMuonPairsWithoutMFTVsPtMix", "Dimuon Mass Mix (MUON only) vs p_{T}^{#mu#mu}", 1000, 0, 10, nBinsPt, ptMin, ptMax); 
+  fHistMassMuonPairsVsPtLSp[kMixedEvent]             = new TH2D("fHistMassMuonPairsVsPtLSpMix", "Dimuon Mass Mix (MUON+MFT) vs p_{T}^{#mu#mu} Like-Sign ++", 1000, 0, 10, nBinsPt, ptMin, ptMax);
+  fHistMassMuonPairsWithoutMFTVsPtLSp[kMixedEvent]    = new TH2D("fHistMassMuonPairsWithoutMFTVsPtLSpMix", "Dimuon Mass Mix (MUON only) vs p_{T}^{#mu#mu} Like-Sign ++", 1000, 0, 10, nBinsPt, ptMin, ptMax); 
+  fHistMassMuonPairsVsPtLSm[kMixedEvent]             = new TH2D("fHistMassMuonPairsVsPtLSmMix", "Dimuon Mass Mix (MUON+MFT) vs p_{T}^{#mu#mu} Like-Sign --", 1000, 0, 10, nBinsPt, ptMin, ptMax);
+  fHistMassMuonPairsWithoutMFTVsPtLSm[kMixedEvent]    = new TH2D("fHistMassMuonPairsWithoutMFTVsPtLSmMix", "Dimuon Mass Mix (MUON only) vs p_{T}^{#mu#mu} Like-Sign --", 1000, 0, 10, nBinsPt, ptMin, ptMax); 
+  fHistWOffsetMuonPairsAtPrimaryVtxVsPt[kMixedEvent]  = new TH2D("fHistWOffsetMuonPairsAtPrimaryVtxVsPtMix", "Weighted Offset for Muon Pairs Mix at Primary Vertex vs p_{T}^{#mu#mu}", 300, 0, 60, nBinsPt, ptMin, ptMax);
+  fHistWOffsetMuonPairsAtPCAVsPt[kMixedEvent]         = new TH2D("fHistWOffsetMuonPairsAtPCAVsPtMix", "Weighted Offset for Muon Pairs Mix at PCA vs p_{T}^{#mu#mu}", 300, 0, 60, nBinsPt, ptMin, ptMax);
+  fHistDistancePrimaryVtxPCAVsPt[kMixedEvent]         = new TH2D("fHistDistancePrimaryVtxPCAVsPtMix", "Distance between PCA and primary vertex Mix vs p_{T}^{#mu#mu}", 1000, 0, 50000, nBinsPt, ptMin, ptMax);
+  fHistPCAQualityVsPt[kMixedEvent]                    = new TH2D("fHistPCAQualityVsPtMix", "PCA Quality Mix vs p_{T}^{#mu#mu}", 200, 0, 1, nBinsPt, ptMin, ptMax);
+  fHistPseudoProperDecayLengthVsPt[kMixedEvent]       = new TH2D("fHistPseudoProperDecayLengthVsPtMix", "Pseudo proper decay length Mix vs p_{T}^{#mu#mu}", 1000, -5000, 5000, nBinsPt, ptMin, ptMax);
+
+  fHistMassMuonPairsMCVsPt      = new TH2D("fHistMassMuonPairsMCVsPt", "Dimuon Mass (MC) vs p_{T}^{#mu#mu}", 1000, 0, 10, nBinsPt, ptMin, ptMax); 
+  fHistDimuonVtxResolutionXVsPt = new TH2D("fHistDimuonVtxResolutionXVsPt", "Dimuon Vtx offset along X w.r.t. true Vtx vs p_{T}^{#mu#mu}", 2000, -1000, 1000, nBinsPt, ptMin, ptMax);
+  fHistDimuonVtxResolutionYVsPt = new TH2D("fHistDimuonVtxResolutionYVsPt", "Dimuon Vtx offset along Y w.r.t. true Vtx vs p_{T}^{#mu#mu}", 2000, -1000, 1000, nBinsPt, ptMin, ptMax);
+  fHistDimuonVtxResolutionZVsPt = new TH2D("fHistDimuonVtxResolutionZVsPt", "Dimuon Vtx offset along Z w.r.t. true Vtx vs p_{T}^{#mu#mu}", 2000, -1000, 1000, nBinsPt, ptMin, ptMax);
     
-  //--------------------------------------------
-
-  fGraphSingleMuonsOffsetChi2 = new TGraph("fGraphSingleMuonsOffsetChi2");
-  fGraphSingleMuonsOffsetChi2 -> SetName("fGraphSingleMuonsOffsetChi2");
-
-  //--------------------------------------------
-
-  fHistWOffsetMuonPairs        = new TH1D("fHistWOffsetMuonPairs",        "Weighted Offset for Muon Pairs", 300, 0, 60);
-  fHistMassMuonPairs          = new TH1D("fHistMassMuonPairs",           "Dimuon Mass (MUON+MFT)",  fNMassBins, fMassMin, fMassMax);
-  fHistMassMuonPairsWithoutMFT = new TH1D("fHistMassMuonPairsWithoutMFT", "Dimuon Mass (MUON only)", fNMassBins, fMassMin, fMassMax);
-  fHistMassMuonPairsMC         = new TH1D("fHistMassMuonPairsMC",         "Dimuon Mass (MC)",        fNMassBins, fMassMin, fMassMax);
-  fHistRapidityPtMuonPairsMC   = new TH2D("fHistRapidityPtMuonPairsMC",   "Dimuon Phase Space (MC)", 100, -4.5, -2., 100, 0., 10.); 
-
-  fHistWOffsetMuonPairs        -> SetXTitle("Weighted Offset");
-  fHistMassMuonPairs          -> SetXTitle("Mass  [GeV/c^{2}]");
-  fHistMassMuonPairsWithoutMFT -> SetXTitle("Mass  [GeV/c^{2}]");
-  fHistMassMuonPairsMC         -> SetXTitle("Mass  [GeV/c^{2}]");
-  fHistRapidityPtMuonPairsMC   -> SetXTitle("Rapidity");
-  fHistRapidityPtMuonPairsMC   -> SetYTitle("p_{T}  [GeV/c]");
-
-  fHistWOffsetMuonPairs        -> Sumw2();
-  fHistMassMuonPairs          -> Sumw2();
-  fHistMassMuonPairsWithoutMFT -> Sumw2();
-  fHistMassMuonPairsMC         -> Sumw2();
-  fHistRapidityPtMuonPairsMC   -> Sumw2();
+  for (Int_t i=0; i<2; i++) {
+    fHistMassMuonPairsVsPt[i]                -> SetXTitle("Mass  [GeV/c^{2}]");
+    fHistMassMuonPairsWithoutMFTVsPt[i]      -> SetXTitle("Mass  [GeV/c^{2}]");
+    fHistMassMuonPairsVsPtLSp[i]             -> SetXTitle("Mass  [GeV/c^{2}]");
+    fHistMassMuonPairsWithoutMFTVsPtLSp[i]   -> SetXTitle("Mass  [GeV/c^{2}]");
+    fHistMassMuonPairsVsPtLSm[i]             -> SetXTitle("Mass  [GeV/c^{2}]");
+    fHistMassMuonPairsWithoutMFTVsPtLSm[i]   -> SetXTitle("Mass  [GeV/c^{2}]");
+    fHistWOffsetMuonPairsAtPrimaryVtxVsPt[i] -> SetXTitle("Weighted Offset at Primary Vtx");
+    fHistWOffsetMuonPairsAtPCAVsPt[i]        -> SetXTitle("Weighted Offset at PCA");
+    fHistDistancePrimaryVtxPCAVsPt[i]        -> SetXTitle("PCA - Primary Vtx [#mum]");
+    fHistPCAQualityVsPt[i]                   -> SetXTitle("PCA Quality");
+    fHistPseudoProperDecayLengthVsPt[i]      -> SetXTitle("l_{xy}  [#mum]");
+  }
+  fHistMassMuonPairsMCVsPt              -> SetXTitle("Mass  [GeV/c^{2}]");    
+  fHistDimuonVtxResolutionXVsPt         -> SetXTitle("Offset  [#mum]");
+  fHistDimuonVtxResolutionYVsPt         -> SetXTitle("Offset  [#mum]");
+  fHistDimuonVtxResolutionZVsPt         -> SetXTitle("Offset  [#mum]");
+
+  for (Int_t i=0; i<2; i++) {
+    fHistMassMuonPairsVsPt[i]                -> SetYTitle("p_{T}  [GeV/c]");
+    fHistMassMuonPairsWithoutMFTVsPt[i]      -> SetYTitle("p_{T}  [GeV/c]");
+    fHistMassMuonPairsVsPtLSp[i]             -> SetYTitle("p_{T}  [GeV/c]");
+    fHistMassMuonPairsWithoutMFTVsPtLSp[i]   -> SetYTitle("p_{T}  [GeV/c]");
+    fHistMassMuonPairsVsPtLSm[i]             -> SetYTitle("p_{T}  [GeV/c]");
+    fHistMassMuonPairsWithoutMFTVsPtLSm[i]   -> SetYTitle("p_{T}  [GeV/c]");
+    fHistWOffsetMuonPairsAtPrimaryVtxVsPt[i] -> SetYTitle("p_{T}  [GeV/c]");
+    fHistWOffsetMuonPairsAtPCAVsPt[i]        -> SetYTitle("p_{T}  [GeV/c]");
+    fHistDistancePrimaryVtxPCAVsPt[i]        -> SetYTitle("p_{T}  [GeV/c]");  
+    fHistPCAQualityVsPt[i]                   -> SetYTitle("p_{T}  [GeV/c]");
+    fHistPseudoProperDecayLengthVsPt[i]      -> SetYTitle("p_{T}  [GeV/c]");
+  }
+  fHistMassMuonPairsMCVsPt              -> SetYTitle("p_{T}  [GeV/c]");
+  fHistDimuonVtxResolutionXVsPt         -> SetYTitle("p_{T}  [GeV/c]");
+  fHistDimuonVtxResolutionYVsPt         -> SetYTitle("p_{T}  [GeV/c]");
+  fHistDimuonVtxResolutionZVsPt         -> SetYTitle("p_{T}  [GeV/c]");
+
+  for (Int_t i=0; i<2; i++) {
+    fHistMassMuonPairsVsPt[i]                -> Sumw2();
+    fHistMassMuonPairsWithoutMFTVsPt[i]      -> Sumw2();
+    fHistMassMuonPairsVsPtLSp[i]             -> Sumw2();
+    fHistMassMuonPairsWithoutMFTVsPtLSp[i]   -> Sumw2();
+    fHistMassMuonPairsVsPtLSm[i]             -> Sumw2();
+    fHistMassMuonPairsWithoutMFTVsPtLSm[i]   -> Sumw2();
+    fHistWOffsetMuonPairsAtPrimaryVtxVsPt[i] -> Sumw2();
+    fHistWOffsetMuonPairsAtPCAVsPt[i]        -> Sumw2();
+    fHistDistancePrimaryVtxPCAVsPt[i]        -> Sumw2();
+    fHistPCAQualityVsPt[i]                   -> Sumw2();
+    fHistPseudoProperDecayLengthVsPt[i]      -> Sumw2();
+  }
+  fHistMassMuonPairsMCVsPt              -> Sumw2();    
+  fHistDimuonVtxResolutionXVsPt         -> Sumw2();
+  fHistDimuonVtxResolutionYVsPt         -> Sumw2();
+  fHistDimuonVtxResolutionZVsPt         -> Sumw2();
+
+  fHistRapidityPtMuonPairs[kSingleEvent] = new TH2D("fHistRapidityPtMuonPairs", "Dimuon Phase Space (rec)", 20, -4.5, -2., 20, 0., 10.); 
+  fHistRapidityPtMuonPairs[kSingleEvent] -> SetXTitle("Rapidity");
+  fHistRapidityPtMuonPairs[kSingleEvent] -> SetYTitle("p_{T}  [GeV/c]");
+  fHistRapidityPtMuonPairs[kSingleEvent] -> Sumw2();
+
+  fHistRapidityPtMuonPairs[kMixedEvent] = new TH2D("fHistRapidityPtMuonPairsMix", "Dimuon Phase Space (rec)", 20, -4.5, -2., 20, 0., 10.); 
+  fHistRapidityPtMuonPairs[kMixedEvent] -> SetXTitle("Rapidity");
+  fHistRapidityPtMuonPairs[kMixedEvent] -> SetYTitle("p_{T}  [GeV/c]");
+  fHistRapidityPtMuonPairs[kMixedEvent] -> Sumw2();
 
 }