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