Analysis code updated
[u/mrichter/AliRoot.git] / MFT / AliMuonForwardTrackAnalysis.cxx
1 //====================================================================================================================================================
2 //
3 //      Class for the analysis of the ALICE muon forward tracks (MUON + MFT)
4 //
5 //      Contact author: antonio.uras@cern.ch
6 //
7 //====================================================================================================================================================
8
9 #include "TObject.h"
10 #include "TClonesArray.h"
11 #include "AliMuonForwardTrack.h"
12 #include "AliMuonForwardTrackPair.h"
13 #include "TMatrixD.h"
14 #include "TTree.h"
15 #include "TH1D.h"
16 #include "AliLog.h"
17 #include "TFile.h"
18 #include "TParticle.h"
19 #include "AliMUONTrackParam.h"
20 #include "AliMUONTrackExtrap.h"
21 #include "TGeoManager.h"
22 #include "TRandom.h"
23 #include "TLorentzVector.h"
24 #include "TDatabasePDG.h"
25 #include "TGraph.h"
26 #include "AliMuonForwardTrackAnalysis.h"
27 #include "TSystem.h"
28
29 ClassImp(AliMuonForwardTrackAnalysis)
30
31 //====================================================================================================================================================
32
33 AliMuonForwardTrackAnalysis::AliMuonForwardTrackAnalysis():
34   TObject(),
35   fInputDir(0),
36   fOutputDir(0),
37   fInputTreeWithBranson(0x0),
38   fInputTreeWithoutBranson(0x0),
39   fMuonForwardTracksWithBranson(0),
40   fMuonForwardTracksWithBransonMix(0),
41   fMuonForwardTrackPairsWithBranson(0),
42   fMuonForwardTracksWithoutBranson(0),
43   fMuonForwardTracksWithoutBransonMix(0),
44   fMuonForwardTrackPairsWithoutBranson(0),
45   fMFTTrackWithBranson(0),
46   fMFTTrackWithoutBranson(0),
47   fMFTTrack(0),
48   fMFTTrackPairWithBranson(0),
49   fMFTTrackPairWithoutBranson(0),
50   fMFTTrackPair(0),
51   fMCRefTrack(0),
52   fEv(0),
53   fFirstEvent(-1),
54   fLastEvent(-1),
55   fNTracksOfEvent(0),
56   fNTracksAnalyzedOfEvent(0),
57   fNTracksAnalyzed(0),
58   fNPairsOfEvent(0),
59   fNPairsAnalyzedOfEvent(0),
60   fNTracksAnalyzedOfEventAfterCut(0),
61   fNPairsAnalyzedOfEventAfterCut(0),
62   fHistXOffsetSingleMuonsVsP(0x0),
63   fHistYOffsetSingleMuonsVsP(0x0),
64   fHistOffsetSingleMuonsVsP(0x0),
65   fHistWOffsetSingleMuonsVsP(0x0),
66   fHistXOffsetSingleMuonsVsPt(0x0),
67   fHistYOffsetSingleMuonsVsPt(0x0),
68   fHistOffsetSingleMuonsVsPt(0x0),
69   fHistWOffsetSingleMuonsVsPt(0x0),
70   fHistXErrorSingleMuonsVsP(0x0),
71   fHistYErrorSingleMuonsVsP(0x0),
72   fHistXErrorSingleMuonsVsPt(0x0),
73   fHistYErrorSingleMuonsVsPt(0x0),
74   fHistZOriginSingleMuonsMC(0x0),
75   fHistZROriginSingleMuonsMC(0x0), 
76   fHistSingleMuonsPtRapidity(0x0), 
77   fHistSingleMuonsOffsetChi2(0x0),
78   fHistMassMuonPairsMCVsPt(0x0),
79   fHistDimuonVtxResolutionXVsPt(0x0),
80   fHistDimuonVtxResolutionYVsPt(0x0),
81   fHistDimuonVtxResolutionZVsPt(0x0),
82   fTrueMass(0.),
83   fMassMin(0),
84   fMassMax(9999),
85   fSingleMuonAnalysis(1),
86   fMuonPairAnalysis(1),
87   fOption(0),
88   fTriggerLevel(0),
89   fXVertResMC(50.e-4),
90   fYVertResMC(50.e-4),
91   fZVertResMC(50.e-4),
92   fPrimaryVtxX(0.),
93   fPrimaryVtxY(0.),
94   fPrimaryVtxZ(0.),
95   fMaxNWrongClustersMC(999),
96   fMinPtSingleMuons(0),
97   fUseBransonForCut(kFALSE),
98   fUseBransonForKinematics(kFALSE),
99   fCorrelateCutOnOffsetChi2(kFALSE),
100   fMaxChi2SingleMuons(1.e9), 
101   fMaxOffsetSingleMuons(1.e9),
102   fMaxWOffsetMuonPairsAtPrimaryVtx(1.e9), 
103   fMaxWOffsetMuonPairsAtPCA(1.e9), 
104   fMaxDistancePrimaryVtxPCA(1.e9), 
105   fMinPCAQuality(0.),
106   fMixing(kFALSE),
107   fNEventsToMix(10)
108 {
109
110   // default constructor
111
112   for (Int_t i=0; i<2; i++) {
113     fHistRapidityPtMuonPairs[i]              = 0;
114     fHistMassMuonPairsVsPt[i]                = 0;
115     fHistMassMuonPairsWithoutMFTVsPt[i]      = 0;
116     fHistMassMuonPairsVsPtLSp[i]             = 0;
117     fHistMassMuonPairsWithoutMFTVsPtLSp[i]   = 0;
118     fHistMassMuonPairsVsPtLSm[i]             = 0;
119     fHistMassMuonPairsWithoutMFTVsPtLSm[i]   = 0;
120     fHistWOffsetMuonPairsAtPrimaryVtxVsPt[i] = 0;
121     fHistWOffsetMuonPairsAtPCAVsPt[i]        = 0;
122     fHistDistancePrimaryVtxPCAVsPt[i]        = 0;
123     fHistPCAQualityVsPt[i]                   = 0;
124     fHistPseudoProperDecayLengthVsPt[i]      = 0;
125   }
126
127 }
128
129 //====================================================================================================================================================
130
131 Bool_t AliMuonForwardTrackAnalysis::Init(Char_t *inputFileName) {
132
133   BookHistos();
134   
135   TFile *inputFile = new TFile(Form("%s/%s",fInputDir.Data(),inputFileName));
136   if (!inputFile || !inputFile->IsOpen()) {
137     AliError(Form("Error opening file %s", inputFileName));
138     return kFALSE;
139   }
140   
141   fInputTreeWithBranson    = (TTree*) inputFile->Get("AliMuonForwardTracksWithBranson");
142   if (!fInputTreeWithBranson) {
143     AliError("Error reading input tree");
144     return kFALSE;
145   }
146   
147   fInputTreeWithoutBranson    = (TTree*) inputFile->Get("AliMuonForwardTracksWithoutBranson");
148   if (!fInputTreeWithoutBranson) {
149     AliError("Error reading input tree");
150     return kFALSE;
151   }
152   
153   if (fFirstEvent>=fInputTreeWithBranson->GetEntriesFast()) return kFALSE;
154   else if (fFirstEvent<0 || fLastEvent<0 || fFirstEvent>fLastEvent || fFirstEvent>=fInputTreeWithBranson->GetEntriesFast()) {
155     fFirstEvent = 0;
156     fLastEvent  = fInputTreeWithBranson->GetEntriesFast()-1;
157   }
158   else {
159     fLastEvent = TMath::Min(fLastEvent, Int_t(fInputTreeWithBranson->GetEntriesFast()-1));
160   }
161   
162   AliInfo(Form("Analysing events %d to %d", fFirstEvent, fLastEvent));
163   
164   fMuonForwardTracksWithBranson = new TClonesArray("AliMuonForwardTrack",30);
165   fInputTreeWithBranson->SetBranchAddress("tracks", &fMuonForwardTracksWithBranson);  
166   
167   fMuonForwardTracksWithoutBranson = new TClonesArray("AliMuonForwardTrack",30);
168   fInputTreeWithoutBranson->SetBranchAddress("tracks", &fMuonForwardTracksWithoutBranson);  
169   
170   TGeoManager::Import(Form("%s/geometry.root",fInputDir.Data()));
171   
172   AliMUONTrackExtrap::SetField();
173   
174   fMuonForwardTrackPairsWithBranson    = new TClonesArray("AliMuonForwardTrackPair",10);
175   fMuonForwardTrackPairsWithoutBranson = new TClonesArray("AliMuonForwardTrackPair",10);
176   fMuonForwardTrackPairsWithBranson    -> SetOwner(kTRUE);
177   fMuonForwardTrackPairsWithoutBranson -> SetOwner(kTRUE);
178   
179   return kTRUE;
180   
181 }
182   
183 //====================================================================================================================================================
184
185 Bool_t AliMuonForwardTrackAnalysis::LoadNextEvent() {
186
187   if (fEv>fLastEvent) return kFALSE;
188   if (fEv<fFirstEvent) { fEv++; return kTRUE; }
189
190   fMuonForwardTracksWithBranson    -> Clear("C");
191   fMuonForwardTracksWithoutBranson -> Clear("C");
192   fMuonForwardTracksWithBranson    -> Delete();
193   fMuonForwardTracksWithoutBranson -> Delete();
194   fInputTreeWithBranson    -> GetEvent(fEv);
195   fInputTreeWithoutBranson -> GetEvent(fEv);
196
197   AliDebug(2,Form("**** analyzing event # %4d (%3d tracks) ****", fEv, fMuonForwardTracksWithBranson->GetEntriesFast()));
198
199   AliInfo(Form("**** analyzing event # %6d of %6d ****", fEv, fLastEvent));
200
201   fPrimaryVtxX = gRandom->Gaus(0., fXVertResMC);
202   fPrimaryVtxY = gRandom->Gaus(0., fYVertResMC);
203   fPrimaryVtxZ = gRandom->Gaus(0., fZVertResMC);
204
205   if (fSingleMuonAnalysis) {
206     fNTracksAnalyzedOfEvent = 0;
207     fNTracksAnalyzedOfEventAfterCut = 0;
208     fNTracksOfEvent = fMuonForwardTracksWithBranson->GetEntriesFast();
209     while (AnalyzeSingleMuon()) continue;
210   }
211   
212   if (fMuonPairAnalysis) {
213     if (fMuonForwardTrackPairsWithBranson || fMuonForwardTrackPairsWithoutBranson) {
214       fMuonForwardTrackPairsWithBranson    -> Clear("C");
215       fMuonForwardTrackPairsWithoutBranson -> Clear("C");
216       fMuonForwardTrackPairsWithBranson    -> Delete();
217       fMuonForwardTrackPairsWithoutBranson -> Delete();
218     }
219     BuildMuonPairs();
220     fNPairsAnalyzedOfEvent = 0;
221     fNPairsAnalyzedOfEventAfterCut = 0;
222     fNPairsOfEvent = fMuonForwardTrackPairsWithBranson->GetEntriesFast();
223     while (AnalyzeMuonPair(kSingleEvent)) continue;
224   }
225
226   AliDebug(2,Form("**** analyzed  event # %4d (%3d tracks and %3d pairs analyzed) ****", fEv, fNTracksAnalyzedOfEventAfterCut, fNPairsAnalyzedOfEventAfterCut));
227   
228   if (fMuonPairAnalysis && fMixing) {
229     for (fEvMix=fEv+1; fEvMix<=fEv+fNEventsToMix; fEvMix++) {
230       if (fEvMix<=fLastEvent) {
231
232         if (fMuonForwardTracksWithBransonMix || fMuonForwardTracksWithoutBransonMix) {
233           fMuonForwardTracksWithBransonMix    -> Delete();
234           fMuonForwardTracksWithoutBransonMix -> Delete();
235           delete fMuonForwardTracksWithBransonMix;
236           delete fMuonForwardTracksWithoutBransonMix;
237         }
238
239         fMuonForwardTracksWithBranson    -> Clear("C");
240         fMuonForwardTracksWithoutBranson -> Clear("C");
241         fMuonForwardTracksWithBranson    -> Delete();
242         fMuonForwardTracksWithoutBranson -> Delete();
243         fInputTreeWithBranson    -> GetEvent(fEvMix);
244         fInputTreeWithoutBranson -> GetEvent(fEvMix);
245
246         fMuonForwardTracksWithBransonMix    = new TClonesArray(*fMuonForwardTracksWithBranson);
247         fMuonForwardTracksWithoutBransonMix = new TClonesArray(*fMuonForwardTracksWithoutBranson);
248
249         fMuonForwardTracksWithBranson    -> Clear("C");
250         fMuonForwardTracksWithoutBranson -> Clear("C");
251         fMuonForwardTracksWithBranson    -> Delete();
252         fMuonForwardTracksWithoutBranson -> Delete();
253         fInputTreeWithBranson    -> GetEvent(fEv);
254         fInputTreeWithoutBranson -> GetEvent(fEv);
255
256         AliDebug(2,Form("**** mixing event # %4d (%3d tracks) with event # %4d (%3d tracks) ****", 
257                         fEv,    fMuonForwardTracksWithBranson->GetEntriesFast(),
258                         fEvMix, fMuonForwardTracksWithBransonMix ->GetEntriesFast()));    
259         if (fMuonForwardTrackPairsWithBranson || fMuonForwardTrackPairsWithoutBranson) {
260           fMuonForwardTrackPairsWithBranson    -> Clear("C");
261           fMuonForwardTrackPairsWithoutBranson -> Clear("C");
262           fMuonForwardTrackPairsWithBranson    -> Delete();
263           fMuonForwardTrackPairsWithoutBranson -> Delete();
264         }
265         BuildMuonPairsMix();
266         fNPairsAnalyzedOfEvent = 0;
267         fNPairsAnalyzedOfEventAfterCut = 0;
268         fNPairsOfEvent = fMuonForwardTrackPairsWithBranson->GetEntriesFast();
269         while (AnalyzeMuonPair(kMixedEvent)) continue;
270         AliDebug(2,Form("**** analyzed  mixed event pair (%4d ; %4d) : %3d pairs analyzed ****", fEv, fEvMix, fNPairsAnalyzedOfEventAfterCut));
271         
272         // delete fMuonForwardTracksWithBransonMix;
273         // delete fMuonForwardTracksWithoutBransonMix;
274
275       }
276     }
277   }
278     
279   fEv++;
280   
281   return kTRUE;
282
283 }
284
285 //====================================================================================================================================================
286
287 Bool_t AliMuonForwardTrackAnalysis::AnalyzeSingleMuon() {
288
289   if (fNTracksAnalyzedOfEvent>=fNTracksOfEvent) return kFALSE;
290
291   fMFTTrackWithBranson    = (AliMuonForwardTrack*) fMuonForwardTracksWithBranson->At(fNTracksAnalyzedOfEvent);
292   fMFTTrackWithoutBranson = (AliMuonForwardTrack*) fMuonForwardTracksWithoutBranson->At(fNTracksAnalyzedOfEvent);
293
294   fNTracksAnalyzedOfEvent++;
295
296   Bool_t passedCut = kFALSE;
297   if (fUseBransonForCut) passedCut = PassedCutSingleMuon(fMFTTrackWithBranson);
298   else passedCut = PassedCutSingleMuon(fMFTTrackWithoutBranson);
299   if (!passedCut) return kTRUE;
300
301   if (fUseBransonForKinematics) fMFTTrack = fMFTTrackWithBranson;
302   else fMFTTrack = fMFTTrackWithoutBranson;
303   if (fMFTTrack->GetMatchTrigger() < fTriggerLevel) return kTRUE;
304   fMCRefTrack = fMFTTrack->GetMCTrackRef();
305
306   if (fMaxNWrongClustersMC<6) {
307     if (!fMCRefTrack) return kTRUE;
308     if (fMFTTrack->GetNWrongClustersMC()>fMaxNWrongClustersMC) return kTRUE;
309   }
310
311   if (fMCRefTrack) {
312     fHistZOriginSingleMuonsMC  -> Fill(-1.*fMCRefTrack->Vz());
313     fHistZROriginSingleMuonsMC -> Fill(-1.*fMCRefTrack->Vz(), TMath::Sqrt(fMCRefTrack->Vx()*fMCRefTrack->Vx()+fMCRefTrack->Vy()*fMCRefTrack->Vy()));
314   }
315
316   AliMUONTrackParam *param = fMFTTrack->GetTrackParamAtMFTCluster(0);
317   AliMUONTrackExtrap::ExtrapToZCov(param, fPrimaryVtxZ);
318
319   TLorentzVector pMu;
320   Double_t mMu = TDatabasePDG::Instance()->GetParticle("mu-")->Mass();
321   Double_t energy = TMath::Sqrt(param->P()*param->P() + mMu*mMu);
322   pMu.SetPxPyPzE(param->Px(), param->Py(), param->Pz(), energy);
323
324   TMatrixD cov(5,5);
325   cov = param->GetCovariances();
326
327   fHistXErrorSingleMuonsVsP  -> Fill(1.e4*TMath::Sqrt(cov(0,0)), pMu.P());
328   fHistYErrorSingleMuonsVsP  -> Fill(1.e4*TMath::Sqrt(cov(2,2)), pMu.P());
329   fHistXErrorSingleMuonsVsPt -> Fill(1.e4*TMath::Sqrt(cov(0,0)), pMu.Pt());
330   fHistYErrorSingleMuonsVsPt -> Fill(1.e4*TMath::Sqrt(cov(2,2)), pMu.Pt());
331
332   Double_t dX = fMFTTrack->GetOffsetX(fPrimaryVtxX, fPrimaryVtxZ);
333   Double_t dY = fMFTTrack->GetOffsetY(fPrimaryVtxY, fPrimaryVtxZ);
334   
335   Double_t offset = fMFTTrack->GetOffset(fPrimaryVtxX, fPrimaryVtxY, fPrimaryVtxZ);
336   Double_t weightedOffset = fMFTTrack->GetWeightedOffset(fPrimaryVtxX, fPrimaryVtxY, fPrimaryVtxZ);
337
338   //  AliDebug(2, Form("pdg code = %d\n", fMCRefTrack->GetPdgCode()));
339
340   fHistXOffsetSingleMuonsVsP  -> Fill(1.e4*dX,        pMu.P());
341   fHistYOffsetSingleMuonsVsP  -> Fill(1.e4*dY,        pMu.P());
342   fHistOffsetSingleMuonsVsP   -> Fill(1.e4*offset,    pMu.P());
343   fHistWOffsetSingleMuonsVsP  -> Fill(weightedOffset, pMu.P());
344   fHistXOffsetSingleMuonsVsPt -> Fill(1.e4*dX,        pMu.Pt());
345   fHistYOffsetSingleMuonsVsPt -> Fill(1.e4*dY,        pMu.Pt());
346   fHistOffsetSingleMuonsVsPt  -> Fill(1.e4*offset,    pMu.Pt());
347   fHistWOffsetSingleMuonsVsPt -> Fill(weightedOffset, pMu.Pt());
348
349   fHistSingleMuonsPtRapidity -> Fill(pMu.Rapidity(), pMu.Pt());
350   Double_t chi2OverNdf = fMFTTrack->GetGlobalChi2()/Double_t(fMFTTrack->GetNMFTClusters()+fMFTTrack->GetNMUONClusters());
351   fHistSingleMuonsOffsetChi2  -> Fill(1.e4*offset, chi2OverNdf);
352
353   fNTracksAnalyzed++;
354   fNTracksAnalyzedOfEventAfterCut++;
355
356   return kTRUE;
357
358 }
359
360 //====================================================================================================================================================
361
362 Bool_t AliMuonForwardTrackAnalysis::AnalyzeMuonPair(Int_t opt) {
363
364   if (fNPairsAnalyzedOfEvent>=fNPairsOfEvent) return kFALSE;
365
366   fMFTTrackPairWithBranson    = (AliMuonForwardTrackPair*) fMuonForwardTrackPairsWithBranson->At(fNPairsAnalyzedOfEvent);
367   fMFTTrackPairWithoutBranson = (AliMuonForwardTrackPair*) fMuonForwardTrackPairsWithoutBranson->At(fNPairsAnalyzedOfEvent);
368
369   fNPairsAnalyzedOfEvent++;
370
371   if (fUseBransonForKinematics) fMFTTrackPair = fMFTTrackPairWithBranson;
372   else fMFTTrackPair = fMFTTrackPairWithoutBranson;
373
374   if ( fMFTTrackPair->GetTrack(0)->GetNWrongClustersMC()>fMaxNWrongClustersMC || 
375        fMFTTrackPair->GetTrack(1)->GetNWrongClustersMC()>fMaxNWrongClustersMC ) return kTRUE;
376
377   if (fOption==kResonanceOnly && !fMFTTrackPair->IsResonance())   return kTRUE;
378   if (fOption==kResonanceOnly && fMFTTrackPair->GetCharge() != 0) return kTRUE;
379
380   Double_t pca[3] = {0};
381   fMFTTrackPair -> GetPointOfClosestApproach(pca);
382   Double_t distancePrimaryVtxPCA = TMath::Sqrt(TMath::Power(fPrimaryVtxX-pca[0],2)+
383                                                TMath::Power(fPrimaryVtxY-pca[1],2)+
384                                                TMath::Power(fPrimaryVtxZ-pca[2],2));
385
386   fMFTTrackPair -> SetKinem(fPrimaryVtxZ);
387
388   Bool_t passedCut = kFALSE;
389   if (fUseBransonForCut) passedCut = PassedCutMuonPair(fMFTTrackPairWithBranson);
390   else passedCut = PassedCutMuonPair(fMFTTrackPairWithoutBranson);
391
392   if (!passedCut) return kTRUE;
393
394   // --------------- Filling dimuon histograms --------------------------------
395
396   if (fMFTTrackPair->GetCharge() == 0) {
397     if (fOption==kResonanceOnly && opt==kSingleEvent) fHistMassMuonPairsMCVsPt -> Fill(fMFTTrackPair->GetMassMC(), fMFTTrackPair->GetPt());
398     fHistMassMuonPairsVsPt[opt]                -> Fill(fMFTTrackPair->GetMass(), fMFTTrackPair->GetPt());
399     fHistMassMuonPairsWithoutMFTVsPt[opt]      -> Fill(fMFTTrackPair->GetMassWithoutMFT(fPrimaryVtxX, fPrimaryVtxY, fPrimaryVtxZ), fMFTTrackPair->GetPt());
400     fHistWOffsetMuonPairsAtPrimaryVtxVsPt[opt] -> Fill(fMFTTrackPair->GetWeightedOffset(fPrimaryVtxX, fPrimaryVtxY, fPrimaryVtxZ), fMFTTrackPair->GetPt());
401     fHistWOffsetMuonPairsAtPCAVsPt[opt]        -> Fill(fMFTTrackPair->GetWeightedOffsetAtPCA(), fMFTTrackPair->GetPt());
402     fHistDistancePrimaryVtxPCAVsPt[opt]        -> Fill(distancePrimaryVtxPCA*1.e4, fMFTTrackPair->GetPt());
403     fHistPCAQualityVsPt[opt]                   -> Fill(fMFTTrackPair->GetPCAQuality(), fMFTTrackPair->GetPt());
404     fHistPseudoProperDecayLengthVsPt[opt]      -> Fill(GetPseudoProperDecayLength(fMFTTrackPair, fTrueMass), fMFTTrackPair->GetPt());
405     if (fEvalDimuonVtxResolution && opt==kSingleEvent) {
406       fHistDimuonVtxResolutionXVsPt->Fill(pca[0]*1.e4, fMFTTrackPair->GetPt());
407       fHistDimuonVtxResolutionYVsPt->Fill(pca[1]*1.e4, fMFTTrackPair->GetPt());
408       fHistDimuonVtxResolutionZVsPt->Fill(pca[2]*1.e4, fMFTTrackPair->GetPt());
409     }
410     fHistRapidityPtMuonPairs[opt] -> Fill(fMFTTrackPair->GetRapidity(), fMFTTrackPair->GetPt());
411   } 
412   else if (fMFTTrackPair->GetCharge() == -2) {
413     fHistMassMuonPairsVsPtLSm[opt]           -> Fill(fMFTTrackPair->GetMass(), fMFTTrackPair->GetPt());
414     fHistMassMuonPairsWithoutMFTVsPtLSm[opt] -> Fill(fMFTTrackPair->GetMassWithoutMFT(fPrimaryVtxX, fPrimaryVtxY, fPrimaryVtxZ), fMFTTrackPair->GetPt());
415   } 
416   else if (fMFTTrackPair->GetCharge() == 2) {
417     fHistMassMuonPairsVsPtLSp[opt]           -> Fill(fMFTTrackPair->GetMass(), fMFTTrackPair->GetPt());
418     fHistMassMuonPairsWithoutMFTVsPtLSp[opt] -> Fill(fMFTTrackPair->GetMassWithoutMFT(fPrimaryVtxX, fPrimaryVtxY, fPrimaryVtxZ), fMFTTrackPair->GetPt());
419   }
420   
421   AliDebug(1, Form("mass = %f   MC = %f", fMFTTrackPair->GetMass(), fMFTTrackPair->GetMassMC()));
422
423   fNPairsAnalyzedOfEventAfterCut++;
424
425   return kTRUE;
426
427 }
428
429 //====================================================================================================================================================
430
431 void AliMuonForwardTrackAnalysis::BuildMuonPairs() {
432
433   Int_t nMuonPairs = 0;
434
435   for (Int_t iTrack=0; iTrack<fMuonForwardTracksWithBranson->GetEntriesFast(); iTrack++) {
436     for (Int_t jTrack=0; jTrack<iTrack; jTrack++) {
437       
438       AliMuonForwardTrack *track0_WithBranson = (AliMuonForwardTrack*) fMuonForwardTracksWithBranson->At(iTrack);
439       AliMuonForwardTrack *track1_WithBranson = (AliMuonForwardTrack*) fMuonForwardTracksWithBranson->At(jTrack);
440       
441       AliMuonForwardTrack *track0_WithoutBranson = (AliMuonForwardTrack*) fMuonForwardTracksWithoutBranson->At(iTrack);
442       AliMuonForwardTrack *track1_WithoutBranson = (AliMuonForwardTrack*) fMuonForwardTracksWithoutBranson->At(jTrack);
443
444       if (track0_WithBranson->GetMatchTrigger()<fTriggerLevel || track1_WithBranson->GetMatchTrigger()<fTriggerLevel) continue;
445
446       // Resonance   = both mu from the same resonance
447       // Open Charm  = both mu from charmed mesons
448       // Open Beauty = both mu from beauty mesons
449       // Background1mu  = at least one mu from background
450       // Background2mu  = both mu from background
451
452       AliMuonForwardTrackPair *trackPairWithBranson = new AliMuonForwardTrackPair(track0_WithBranson, track1_WithBranson);
453
454       if (fOption==kResonanceOnly && !trackPairWithBranson->IsResonance()) { delete trackPairWithBranson; continue; }
455       if (fOption==kNoResonances  &&  trackPairWithBranson->IsResonance()) { delete trackPairWithBranson; continue; }
456       
457       if (fOption==kCharmOnly  && (!track0_WithBranson->IsFromCharm()  || !track1_WithBranson->IsFromCharm()))  { delete trackPairWithBranson; continue; }
458
459       if (fOption==kBeautyOnly && (!track0_WithBranson->IsFromBeauty() || !track1_WithBranson->IsFromBeauty())) { delete trackPairWithBranson; continue; }
460
461       if (fOption==kBackground1mu && (!track0_WithBranson->IsFromBackground() && !track1_WithBranson->IsFromBackground())) { delete trackPairWithBranson; continue; }
462       if (fOption==kBackground2mu && (!track0_WithBranson->IsFromBackground() || !track1_WithBranson->IsFromBackground())) { delete trackPairWithBranson; continue; }
463
464       delete trackPairWithBranson;
465       
466       new ((*fMuonForwardTrackPairsWithBranson)[nMuonPairs]) AliMuonForwardTrackPair(track0_WithBranson, track1_WithBranson);
467       new ((*fMuonForwardTrackPairsWithoutBranson)[nMuonPairs]) AliMuonForwardTrackPair(track0_WithoutBranson, track1_WithoutBranson);
468
469       nMuonPairs++;
470       
471     }
472   }
473
474 }
475
476 //====================================================================================================================================================
477
478 void AliMuonForwardTrackAnalysis::BuildMuonPairsMix() {
479
480   Int_t nMuonPairs = 0;
481
482   for (Int_t iTrack=0; iTrack<fMuonForwardTracksWithBranson->GetEntriesFast(); iTrack++) {
483     for (Int_t jTrack=0; jTrack<fMuonForwardTracksWithBransonMix->GetEntriesFast(); jTrack++) {
484       
485       AliMuonForwardTrack *track0_WithBranson = (AliMuonForwardTrack*) fMuonForwardTracksWithBranson -> At(iTrack);
486       AliMuonForwardTrack *track1_WithBranson = (AliMuonForwardTrack*) fMuonForwardTracksWithBransonMix  -> At(jTrack);
487       
488       AliMuonForwardTrack *track0_WithoutBranson = (AliMuonForwardTrack*) fMuonForwardTracksWithoutBranson -> At(iTrack);
489       AliMuonForwardTrack *track1_WithoutBranson = (AliMuonForwardTrack*) fMuonForwardTracksWithoutBransonMix  -> At(jTrack);
490
491       if (track0_WithBranson->GetMatchTrigger()<fTriggerLevel || track1_WithBranson->GetMatchTrigger()<fTriggerLevel) continue;
492
493       // Resonance   = both mu from the same resonance
494       // Open Charm  = both mu from charmed mesons
495       // Open Beauty = both mu from beauty mesons
496       // Background1mu  = at least one mu from background
497       // Background2mu  = both mu from background
498
499       AliMuonForwardTrackPair *trackPairWithBranson = new AliMuonForwardTrackPair(track0_WithBranson, track1_WithBranson);
500
501       if (fOption==kResonanceOnly && !trackPairWithBranson->IsResonance()) { delete trackPairWithBranson; continue; }
502       if (fOption==kNoResonances  &&  trackPairWithBranson->IsResonance()) { delete trackPairWithBranson; continue; }
503       
504       if (fOption==kCharmOnly  && (!track0_WithBranson->IsFromCharm()  || !track1_WithBranson->IsFromCharm()))  { delete trackPairWithBranson; continue; }
505
506       if (fOption==kBeautyOnly && (!track0_WithBranson->IsFromBeauty() || !track1_WithBranson->IsFromBeauty())) { delete trackPairWithBranson; continue; }
507
508       if (fOption==kBackground1mu && (!track0_WithBranson->IsFromBackground() && !track1_WithBranson->IsFromBackground())) { delete trackPairWithBranson; continue; }
509       if (fOption==kBackground2mu && (!track0_WithBranson->IsFromBackground() || !track1_WithBranson->IsFromBackground())) { delete trackPairWithBranson; continue; }
510
511       delete trackPairWithBranson;
512
513       new ((*fMuonForwardTrackPairsWithBranson)[nMuonPairs])    AliMuonForwardTrackPair(track0_WithBranson,    track1_WithBranson);
514       new ((*fMuonForwardTrackPairsWithoutBranson)[nMuonPairs]) AliMuonForwardTrackPair(track0_WithoutBranson, track1_WithoutBranson);
515       
516       nMuonPairs++;
517       
518     }
519   }
520
521 }
522
523 //====================================================================================================================================================
524
525 Bool_t AliMuonForwardTrackAnalysis::PassedCutSingleMuon(AliMuonForwardTrack *track) {
526
527   AliMUONTrackParam *param = track->GetTrackParamAtMFTCluster(0);
528   AliMUONTrackExtrap::ExtrapToZCov(param, fPrimaryVtxZ);
529
530   if (track->Pt()<fMinPtSingleMuons) return kFALSE;
531   
532   Double_t offset = 1.e4*track->GetOffset(fPrimaryVtxX, fPrimaryVtxY, fPrimaryVtxZ);
533   Double_t chi2OverNdf = track->GetGlobalChi2() / Double_t(track->GetNMFTClusters()+track->GetNMUONClusters()); 
534   offset /= fMaxOffsetSingleMuons;
535   chi2OverNdf /= fMaxChi2SingleMuons;
536
537   if (fCorrelateCutOnOffsetChi2) {
538     if (TMath::Sqrt(offset*offset + chi2OverNdf*chi2OverNdf) > 1) return kFALSE;
539   }
540   else {
541     if (offset>1 || chi2OverNdf>1) return kFALSE;
542   }
543
544   return kTRUE;
545
546 }
547
548 //====================================================================================================================================================
549
550 Bool_t AliMuonForwardTrackAnalysis::PassedCutMuonPair(AliMuonForwardTrackPair *pair) {
551
552   if (!PassedCutSingleMuon(pair->GetTrack(0)) || !PassedCutSingleMuon(pair->GetTrack(1))) return kFALSE;
553
554   if (pair->GetMass()>fMassMax || pair->GetMass()<fMassMin) return kFALSE;
555
556   if (pair->GetWeightedOffset(fPrimaryVtxX, fPrimaryVtxY, fPrimaryVtxZ) > fMaxWOffsetMuonPairsAtPrimaryVtx) return kFALSE;
557   if (pair->GetWeightedOffsetAtPCA() > fMaxWOffsetMuonPairsAtPCA) return kFALSE;
558
559   Double_t pca[3] = {0};
560   pair -> GetPointOfClosestApproach(pca);
561   Double_t distancePrimaryVtxPCA = TMath::Sqrt(TMath::Power(fPrimaryVtxX-pca[0],2)+
562                                                TMath::Power(fPrimaryVtxY-pca[1],2)+
563                                                TMath::Power(fPrimaryVtxZ-pca[2],2));
564
565   if (distancePrimaryVtxPCA > fMaxDistancePrimaryVtxPCA) return kFALSE;
566   if (pair->GetPCAQuality() < fMinPCAQuality) return kFALSE;
567
568   return kTRUE;
569
570 }
571
572 //====================================================================================================================================================
573
574 Double_t AliMuonForwardTrackAnalysis::GetPseudoProperDecayLength(AliMuonForwardTrackPair *pair, Double_t trueMass) {
575
576   Double_t pca[3] = {0};
577   pair -> GetPointOfClosestApproach(pca);
578
579   TVector2 vecVertexPCA(pca[0]-fPrimaryVtxX, pca[1]-fPrimaryVtxY);
580   TVector2 dimuonPt(pair->GetPx(), pair->GetPy());
581   TVector2 dimuonPtUnit = dimuonPt.Unit();
582   
583   Double_t l_xy = vecVertexPCA*dimuonPtUnit;
584   Double_t pseudoProperDecayLength = (l_xy * trueMass / pair->GetPt()) * 1.e4 ; // in micron
585
586   return pseudoProperDecayLength;
587   
588 }
589
590 //====================================================================================================================================================
591
592 void AliMuonForwardTrackAnalysis::Terminate(Char_t *outputFileName) {
593
594   TFile *fileOut = new TFile(Form("%s/%s",fOutputDir.Data(),outputFileName), "recreate");
595
596   printf("Writing output objects to file %s\n", fileOut->GetName());
597
598   // single muons
599
600   fHistXOffsetSingleMuonsVsP  -> Write();
601   fHistYOffsetSingleMuonsVsP  -> Write();
602   fHistXOffsetSingleMuonsVsPt -> Write();
603   fHistYOffsetSingleMuonsVsPt -> Write();
604
605   fHistXErrorSingleMuonsVsP   -> Write();
606   fHistYErrorSingleMuonsVsP   -> Write();
607   fHistXErrorSingleMuonsVsPt  -> Write();
608   fHistYErrorSingleMuonsVsPt  -> Write();
609
610   fHistOffsetSingleMuonsVsP   -> Write();
611   fHistOffsetSingleMuonsVsPt  -> Write();
612   fHistWOffsetSingleMuonsVsP  -> Write();
613   fHistWOffsetSingleMuonsVsPt -> Write();
614
615   fHistSingleMuonsPtRapidity  -> Write();
616   fHistSingleMuonsOffsetChi2  -> Write();
617   fHistZOriginSingleMuonsMC   -> Write();
618   fHistZROriginSingleMuonsMC  -> Write();
619
620   // dimuons
621
622   fHistMassMuonPairsVsPt[kSingleEvent]                -> Write();
623   fHistMassMuonPairsWithoutMFTVsPt[kSingleEvent]      -> Write();
624   fHistMassMuonPairsVsPtLSp[kSingleEvent]             -> Write();
625   fHistMassMuonPairsWithoutMFTVsPtLSp[kSingleEvent]   -> Write();
626   fHistMassMuonPairsVsPtLSm[kSingleEvent]             -> Write();
627   fHistMassMuonPairsWithoutMFTVsPtLSm[kSingleEvent]   -> Write();
628   fHistWOffsetMuonPairsAtPrimaryVtxVsPt[kSingleEvent] -> Write();
629   fHistWOffsetMuonPairsAtPCAVsPt[kSingleEvent]        -> Write();
630   fHistDistancePrimaryVtxPCAVsPt[kSingleEvent]        -> Write();
631   fHistPCAQualityVsPt[kSingleEvent]                   -> Write();
632   fHistPseudoProperDecayLengthVsPt[kSingleEvent]      -> Write();
633   fHistRapidityPtMuonPairs[kSingleEvent]              -> Write();
634   fHistMassMuonPairsMCVsPt                            -> Write();
635   if (fEvalDimuonVtxResolution) {
636     fHistDimuonVtxResolutionXVsPt -> Write();
637     fHistDimuonVtxResolutionYVsPt -> Write();
638     fHistDimuonVtxResolutionZVsPt -> Write();
639   }
640
641   if (fMixing) {
642     fHistMassMuonPairsVsPt[kMixedEvent]                -> Write();
643     fHistMassMuonPairsWithoutMFTVsPt[kMixedEvent]      -> Write();
644     fHistMassMuonPairsVsPtLSp[kMixedEvent]             -> Write();
645     fHistMassMuonPairsWithoutMFTVsPtLSp[kMixedEvent]   -> Write();
646     fHistMassMuonPairsVsPtLSm[kMixedEvent]             -> Write();
647     fHistMassMuonPairsWithoutMFTVsPtLSm[kMixedEvent]   -> Write();
648     fHistWOffsetMuonPairsAtPrimaryVtxVsPt[kMixedEvent] -> Write();
649     fHistWOffsetMuonPairsAtPCAVsPt[kMixedEvent]        -> Write();
650     fHistDistancePrimaryVtxPCAVsPt[kMixedEvent]        -> Write();
651     fHistPCAQualityVsPt[kMixedEvent]                   -> Write();
652     fHistPseudoProperDecayLengthVsPt[kMixedEvent]      -> Write();
653     fHistRapidityPtMuonPairs[kMixedEvent]              -> Write();
654   }
655
656   fileOut -> Close();
657
658 }
659
660 //====================================================================================================================================================
661
662 void AliMuonForwardTrackAnalysis::BookHistos() {
663
664   // -------------------- single muons
665
666   fHistXOffsetSingleMuonsVsP  = new TH2D("fHistXOffsetSingleMuonsVsP",  "Offset for single muons along X vs P",      200, -1000, 1000, 100, 0, 100);
667   fHistYOffsetSingleMuonsVsP  = new TH2D("fHistYOffsetSingleMuonsVsP",  "Offset for single muons along Y vs P",      200, -1000, 1000, 100, 0, 100);
668   fHistXOffsetSingleMuonsVsPt = new TH2D("fHistXOffsetSingleMuonsVsPt", "Offset for single muons along X vs p_{T}",  200, -1000, 1000, 100, 0,  10);
669   fHistYOffsetSingleMuonsVsPt = new TH2D("fHistYOffsetSingleMuonsVsPt", "Offset for single muons along Y vs p_{T}",  200, -1000, 1000, 100, 0,  10);
670
671   fHistXErrorSingleMuonsVsP   = new TH2D("fHistXErrorSingleMuonsVsP",   "Coordinate Error for single muons along X vs P",      200, 0, 1000, 100, 0, 100);
672   fHistYErrorSingleMuonsVsP   = new TH2D("fHistYErrorSingleMuonsVsP",   "Coordinate Error for single muons along Y vs P",      200, 0, 1000, 100, 0, 100);
673   fHistXErrorSingleMuonsVsPt  = new TH2D("fHistXErrorSingleMuonsVsPt",  "Coordinate Error for single muons along X vs p_{T}",  200, 0, 1000, 100, 0,  10);
674   fHistYErrorSingleMuonsVsPt  = new TH2D("fHistYErrorSingleMuonsVsPt",  "Coordinate Error for single muons along Y vs p_{T}",  200, 0, 1000, 100, 0,  10);
675
676   fHistOffsetSingleMuonsVsP   = new TH2D("fHistOffsetSingleMuonsVsP",   "Offset for single muons vs P",              200, 0, 2000, 100, 0, 100);
677   fHistOffsetSingleMuonsVsPt  = new TH2D("fHistOffsetSingleMuonsVsPt",  "Offset for single muons vs p_{T}",          200, 0, 2000, 100, 0,  10);
678   fHistWOffsetSingleMuonsVsP  = new TH2D("fHistWOffsetSingleMuonsVsP",  "Weighted Offset for single muons vs P",     300, 0,   15, 100, 0, 100);  
679   fHistWOffsetSingleMuonsVsPt = new TH2D("fHistWOffsetSingleMuonsVsPt", "Weighted Offset for single muons vs p_{T}", 300, 0,   15, 100, 0,  10);  
680
681   fHistSingleMuonsPtRapidity = new TH2D("fHistSingleMuonsPtRapidity", "Phase Space for single muons", 100, -4.5, -2., 100, 0., 10.);
682   fHistSingleMuonsOffsetChi2 = new TH2D("fHistSingleMuonsOffsetChi2", "Offset vs #chi^{2}/ndf for single muons", 400, 0, 4000, 200, 0, 10);
683   fHistZOriginSingleMuonsMC  = new TH1D("fHistZOriginSingleMuonsMC",  "Z origin for single muons (from MC)",   1000, 0., 500.);
684   fHistZROriginSingleMuonsMC = new TH2D("fHistZROriginSingleMuonsMC", "Z-R origin for single muons (from MC)", 1000, 0., 500., 1000, 0., 100.);
685
686   //
687
688   fHistXOffsetSingleMuonsVsP  -> SetXTitle("Offset(X)  [#mum]");
689   fHistYOffsetSingleMuonsVsP  -> SetXTitle("Offset(Y)  [#mum]");
690   fHistXOffsetSingleMuonsVsPt -> SetXTitle("Offset(X)  [#mum]");
691   fHistYOffsetSingleMuonsVsPt -> SetXTitle("Offset(Y)  [#mum]");
692
693   fHistXErrorSingleMuonsVsP   -> SetXTitle("Err. on track position at z_{vtx} (X)  [#mum]");
694   fHistYErrorSingleMuonsVsP   -> SetXTitle("Err. on track position at z_{vtx} (Y)  [#mum]");
695   fHistXErrorSingleMuonsVsPt  -> SetXTitle("Err. on track position at z_{vtx} (X)  [#mum]");
696   fHistYErrorSingleMuonsVsPt  -> SetXTitle("Err. on track position at z_{vtx} (Y)  [#mum]");
697
698   fHistOffsetSingleMuonsVsP   -> SetXTitle("Offset  [#mum]");
699   fHistOffsetSingleMuonsVsPt  -> SetXTitle("Offset  [#mum]");
700   fHistWOffsetSingleMuonsVsP  -> SetXTitle("Weighted Offset");
701   fHistWOffsetSingleMuonsVsPt -> SetXTitle("Weighted Offset");    
702
703   fHistXOffsetSingleMuonsVsP  -> SetYTitle("P  [GeV/c]");
704   fHistYOffsetSingleMuonsVsP  -> SetYTitle("P  [GeV/c]");
705   fHistXOffsetSingleMuonsVsPt -> SetYTitle("p_{T}  [GeV/c]");
706   fHistYOffsetSingleMuonsVsPt -> SetYTitle("p_{T}  [GeV/c]");
707
708   fHistXErrorSingleMuonsVsP   -> SetYTitle("P  [GeV/c]");
709   fHistYErrorSingleMuonsVsP   -> SetYTitle("P  [GeV/c]");
710   fHistXErrorSingleMuonsVsPt  -> SetYTitle("p_{T}  [GeV/c]");
711   fHistYErrorSingleMuonsVsPt  -> SetYTitle("p_{T}  [GeV/c]");
712
713   fHistOffsetSingleMuonsVsP   -> SetYTitle("P  [GeV/c]");
714   fHistOffsetSingleMuonsVsPt  -> SetYTitle("p_{T}  [GeV/c]");
715   fHistWOffsetSingleMuonsVsP  -> SetYTitle("P  [GeV/c]");
716   fHistWOffsetSingleMuonsVsPt -> SetYTitle("p_{T}  [GeV/c]");    
717
718   fHistSingleMuonsPtRapidity -> SetXTitle("y^{#mu}");
719   fHistSingleMuonsPtRapidity -> SetYTitle("p_{T}^{#mu}  [GeV/c]");
720   fHistSingleMuonsOffsetChi2 -> SetXTitle("Offset  [#mum]");
721   fHistSingleMuonsOffsetChi2 -> SetYTitle("#chi^{2}/ndf");
722
723   fHistZOriginSingleMuonsMC  -> SetXTitle("Z  [cm]");
724   fHistZROriginSingleMuonsMC -> SetXTitle("Z  [cm]");
725   fHistZROriginSingleMuonsMC -> SetXTitle("R  [cm]");
726
727   //
728
729   fHistXOffsetSingleMuonsVsP  -> Sumw2();
730   fHistYOffsetSingleMuonsVsP  -> Sumw2();
731   fHistXOffsetSingleMuonsVsPt -> Sumw2();
732   fHistYOffsetSingleMuonsVsPt -> Sumw2();
733
734   fHistXErrorSingleMuonsVsP   -> Sumw2();
735   fHistYErrorSingleMuonsVsP   -> Sumw2();
736   fHistXErrorSingleMuonsVsPt  -> Sumw2();
737   fHistYErrorSingleMuonsVsPt  -> Sumw2();
738
739   fHistOffsetSingleMuonsVsP   -> Sumw2();
740   fHistOffsetSingleMuonsVsPt  -> Sumw2();
741   fHistWOffsetSingleMuonsVsP  -> Sumw2();
742   fHistWOffsetSingleMuonsVsPt -> Sumw2();
743
744   fHistSingleMuonsPtRapidity  -> Sumw2();
745   fHistSingleMuonsOffsetChi2  -> Sumw2();
746   fHistZOriginSingleMuonsMC   -> Sumw2();
747   fHistZROriginSingleMuonsMC  -> Sumw2();
748
749   // -------------------- muon pairs
750
751   Int_t nBinsPt = 20; Double_t ptMin = 0., ptMax = 10.;   // dimuon pt
752
753   fHistMassMuonPairsVsPt[kSingleEvent]                = new TH2D("fHistMassMuonPairsVsPt", "Dimuon Mass (MUON+MFT) vs p_{T}^{#mu#mu}", 1000, 0, 10, nBinsPt, ptMin, ptMax);
754   fHistMassMuonPairsWithoutMFTVsPt[kSingleEvent]      = new TH2D("fHistMassMuonPairsWithoutMFTVsPt", "Dimuon Mass (MUON only) vs p_{T}^{#mu#mu}", 1000, 0, 10, nBinsPt, ptMin, ptMax); 
755   fHistMassMuonPairsVsPtLSp[kSingleEvent]             = new TH2D("fHistMassMuonPairsVsPtLSp", "Dimuon Mass (MUON+MFT) vs p_{T}^{#mu#mu} Like-Sign ++", 1000, 0, 10, nBinsPt, ptMin, ptMax);
756   fHistMassMuonPairsWithoutMFTVsPtLSp[kSingleEvent]   = new TH2D("fHistMassMuonPairsWithoutMFTVsPtLSp", "Dimuon Mass (MUON only) vs p_{T}^{#mu#mu} Like-Sign ++", 1000, 0, 10, nBinsPt, ptMin, ptMax); 
757   fHistMassMuonPairsVsPtLSm[kSingleEvent]             = new TH2D("fHistMassMuonPairsVsPtLSm", "Dimuon Mass (MUON+MFT) vs p_{T}^{#mu#mu} Like-Sign --", 1000, 0, 10, nBinsPt, ptMin, ptMax);
758   fHistMassMuonPairsWithoutMFTVsPtLSm[kSingleEvent]   = new TH2D("fHistMassMuonPairsWithoutMFTVsPtLSm", "Dimuon Mass (MUON only) vs p_{T}^{#mu#mu} Like-Sign --", 1000, 0, 10, nBinsPt, ptMin, ptMax); 
759   fHistWOffsetMuonPairsAtPrimaryVtxVsPt[kSingleEvent] = new TH2D("fHistWOffsetMuonPairsAtPrimaryVtxVsPt", "Weighted Offset for Muon Pairs at Primary Vertex vs p_{T}^{#mu#mu}", 300, 0, 60, nBinsPt, ptMin, ptMax);
760   fHistWOffsetMuonPairsAtPCAVsPt[kSingleEvent]        = new TH2D("fHistWOffsetMuonPairsAtPCAVsPt", "Weighted Offset for Muon Pairs at PCA vs p_{T}^{#mu#mu}", 300, 0, 60, nBinsPt, ptMin, ptMax);
761   fHistDistancePrimaryVtxPCAVsPt[kSingleEvent]        = new TH2D("fHistDistancePrimaryVtxPCAVsPt", "Distance between PCA and primary vertex vs p_{T}^{#mu#mu}", 1000, 0, 50000, nBinsPt, ptMin, ptMax);
762   fHistPCAQualityVsPt[kSingleEvent]                   = new TH2D("fHistPCAQualityVsPt", "PCA Quality vs p_{T}^{#mu#mu}", 200, 0, 1, nBinsPt, ptMin, ptMax);
763   fHistPseudoProperDecayLengthVsPt[kSingleEvent]      = new TH2D("fHistPseudoProperDecayLengthVsPt", "Pseudo proper decay length vs p_{T}^{#mu#mu}", 1000, -5000, 5000, nBinsPt, ptMin, ptMax);
764
765   fHistMassMuonPairsVsPt[kMixedEvent]                 = new TH2D("fHistMassMuonPairsVsPtMix", "Dimuon Mass Mix (MUON+MFT) vs p_{T}^{#mu#mu}", 1000, 0, 10, nBinsPt, ptMin, ptMax);
766   fHistMassMuonPairsWithoutMFTVsPt[kMixedEvent]       = new TH2D("fHistMassMuonPairsWithoutMFTVsPtMix", "Dimuon Mass Mix (MUON only) vs p_{T}^{#mu#mu}", 1000, 0, 10, nBinsPt, ptMin, ptMax); 
767   fHistMassMuonPairsVsPtLSp[kMixedEvent]              = new TH2D("fHistMassMuonPairsVsPtLSpMix", "Dimuon Mass Mix (MUON+MFT) vs p_{T}^{#mu#mu} Like-Sign ++", 1000, 0, 10, nBinsPt, ptMin, ptMax);
768   fHistMassMuonPairsWithoutMFTVsPtLSp[kMixedEvent]    = new TH2D("fHistMassMuonPairsWithoutMFTVsPtLSpMix", "Dimuon Mass Mix (MUON only) vs p_{T}^{#mu#mu} Like-Sign ++", 1000, 0, 10, nBinsPt, ptMin, ptMax); 
769   fHistMassMuonPairsVsPtLSm[kMixedEvent]              = new TH2D("fHistMassMuonPairsVsPtLSmMix", "Dimuon Mass Mix (MUON+MFT) vs p_{T}^{#mu#mu} Like-Sign --", 1000, 0, 10, nBinsPt, ptMin, ptMax);
770   fHistMassMuonPairsWithoutMFTVsPtLSm[kMixedEvent]    = new TH2D("fHistMassMuonPairsWithoutMFTVsPtLSmMix", "Dimuon Mass Mix (MUON only) vs p_{T}^{#mu#mu} Like-Sign --", 1000, 0, 10, nBinsPt, ptMin, ptMax); 
771   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);
772   fHistWOffsetMuonPairsAtPCAVsPt[kMixedEvent]         = new TH2D("fHistWOffsetMuonPairsAtPCAVsPtMix", "Weighted Offset for Muon Pairs Mix at PCA vs p_{T}^{#mu#mu}", 300, 0, 60, nBinsPt, ptMin, ptMax);
773   fHistDistancePrimaryVtxPCAVsPt[kMixedEvent]         = new TH2D("fHistDistancePrimaryVtxPCAVsPtMix", "Distance between PCA and primary vertex Mix vs p_{T}^{#mu#mu}", 1000, 0, 50000, nBinsPt, ptMin, ptMax);
774   fHistPCAQualityVsPt[kMixedEvent]                    = new TH2D("fHistPCAQualityVsPtMix", "PCA Quality Mix vs p_{T}^{#mu#mu}", 200, 0, 1, nBinsPt, ptMin, ptMax);
775   fHistPseudoProperDecayLengthVsPt[kMixedEvent]       = new TH2D("fHistPseudoProperDecayLengthVsPtMix", "Pseudo proper decay length Mix vs p_{T}^{#mu#mu}", 1000, -5000, 5000, nBinsPt, ptMin, ptMax);
776
777   fHistMassMuonPairsMCVsPt      = new TH2D("fHistMassMuonPairsMCVsPt", "Dimuon Mass (MC) vs p_{T}^{#mu#mu}", 1000, 0, 10, nBinsPt, ptMin, ptMax); 
778   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);
779   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);
780   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);
781     
782   for (Int_t i=0; i<2; i++) {
783     fHistMassMuonPairsVsPt[i]                -> SetXTitle("Mass  [GeV/c^{2}]");
784     fHistMassMuonPairsWithoutMFTVsPt[i]      -> SetXTitle("Mass  [GeV/c^{2}]");
785     fHistMassMuonPairsVsPtLSp[i]             -> SetXTitle("Mass  [GeV/c^{2}]");
786     fHistMassMuonPairsWithoutMFTVsPtLSp[i]   -> SetXTitle("Mass  [GeV/c^{2}]");
787     fHistMassMuonPairsVsPtLSm[i]             -> SetXTitle("Mass  [GeV/c^{2}]");
788     fHistMassMuonPairsWithoutMFTVsPtLSm[i]   -> SetXTitle("Mass  [GeV/c^{2}]");
789     fHistWOffsetMuonPairsAtPrimaryVtxVsPt[i] -> SetXTitle("Weighted Offset at Primary Vtx");
790     fHistWOffsetMuonPairsAtPCAVsPt[i]        -> SetXTitle("Weighted Offset at PCA");
791     fHistDistancePrimaryVtxPCAVsPt[i]        -> SetXTitle("PCA - Primary Vtx [#mum]");
792     fHistPCAQualityVsPt[i]                   -> SetXTitle("PCA Quality");
793     fHistPseudoProperDecayLengthVsPt[i]      -> SetXTitle("l_{xy}  [#mum]");
794   }
795   fHistMassMuonPairsMCVsPt              -> SetXTitle("Mass  [GeV/c^{2}]");    
796   fHistDimuonVtxResolutionXVsPt         -> SetXTitle("Offset  [#mum]");
797   fHistDimuonVtxResolutionYVsPt         -> SetXTitle("Offset  [#mum]");
798   fHistDimuonVtxResolutionZVsPt         -> SetXTitle("Offset  [#mum]");
799
800   for (Int_t i=0; i<2; i++) {
801     fHistMassMuonPairsVsPt[i]                -> SetYTitle("p_{T}  [GeV/c]");
802     fHistMassMuonPairsWithoutMFTVsPt[i]      -> SetYTitle("p_{T}  [GeV/c]");
803     fHistMassMuonPairsVsPtLSp[i]             -> SetYTitle("p_{T}  [GeV/c]");
804     fHistMassMuonPairsWithoutMFTVsPtLSp[i]   -> SetYTitle("p_{T}  [GeV/c]");
805     fHistMassMuonPairsVsPtLSm[i]             -> SetYTitle("p_{T}  [GeV/c]");
806     fHistMassMuonPairsWithoutMFTVsPtLSm[i]   -> SetYTitle("p_{T}  [GeV/c]");
807     fHistWOffsetMuonPairsAtPrimaryVtxVsPt[i] -> SetYTitle("p_{T}  [GeV/c]");
808     fHistWOffsetMuonPairsAtPCAVsPt[i]        -> SetYTitle("p_{T}  [GeV/c]");
809     fHistDistancePrimaryVtxPCAVsPt[i]        -> SetYTitle("p_{T}  [GeV/c]");  
810     fHistPCAQualityVsPt[i]                   -> SetYTitle("p_{T}  [GeV/c]");
811     fHistPseudoProperDecayLengthVsPt[i]      -> SetYTitle("p_{T}  [GeV/c]");
812   }
813   fHistMassMuonPairsMCVsPt              -> SetYTitle("p_{T}  [GeV/c]");
814   fHistDimuonVtxResolutionXVsPt         -> SetYTitle("p_{T}  [GeV/c]");
815   fHistDimuonVtxResolutionYVsPt         -> SetYTitle("p_{T}  [GeV/c]");
816   fHistDimuonVtxResolutionZVsPt         -> SetYTitle("p_{T}  [GeV/c]");
817
818   for (Int_t i=0; i<2; i++) {
819     fHistMassMuonPairsVsPt[i]                -> Sumw2();
820     fHistMassMuonPairsWithoutMFTVsPt[i]      -> Sumw2();
821     fHistMassMuonPairsVsPtLSp[i]             -> Sumw2();
822     fHistMassMuonPairsWithoutMFTVsPtLSp[i]   -> Sumw2();
823     fHistMassMuonPairsVsPtLSm[i]             -> Sumw2();
824     fHistMassMuonPairsWithoutMFTVsPtLSm[i]   -> Sumw2();
825     fHistWOffsetMuonPairsAtPrimaryVtxVsPt[i] -> Sumw2();
826     fHistWOffsetMuonPairsAtPCAVsPt[i]        -> Sumw2();
827     fHistDistancePrimaryVtxPCAVsPt[i]        -> Sumw2();
828     fHistPCAQualityVsPt[i]                   -> Sumw2();
829     fHistPseudoProperDecayLengthVsPt[i]      -> Sumw2();
830   }
831   fHistMassMuonPairsMCVsPt              -> Sumw2();    
832   fHistDimuonVtxResolutionXVsPt         -> Sumw2();
833   fHistDimuonVtxResolutionYVsPt         -> Sumw2();
834   fHistDimuonVtxResolutionZVsPt         -> Sumw2();
835
836   fHistRapidityPtMuonPairs[kSingleEvent] = new TH2D("fHistRapidityPtMuonPairs", "Dimuon Phase Space (rec)", 20, -4.5, -2., 20, 0., 10.); 
837   fHistRapidityPtMuonPairs[kSingleEvent] -> SetXTitle("Rapidity");
838   fHistRapidityPtMuonPairs[kSingleEvent] -> SetYTitle("p_{T}  [GeV/c]");
839   fHistRapidityPtMuonPairs[kSingleEvent] -> Sumw2();
840
841   fHistRapidityPtMuonPairs[kMixedEvent] = new TH2D("fHistRapidityPtMuonPairsMix", "Dimuon Phase Space (rec)", 20, -4.5, -2., 20, 0., 10.); 
842   fHistRapidityPtMuonPairs[kMixedEvent] -> SetXTitle("Rapidity");
843   fHistRapidityPtMuonPairs[kMixedEvent] -> SetYTitle("p_{T}  [GeV/c]");
844   fHistRapidityPtMuonPairs[kMixedEvent] -> Sumw2();
845
846 }
847
848 //====================================================================================================================================================
849