- Disentangle masks effect from trigger chamber efficiency estimation.
[u/mrichter/AliRoot.git] / MUON / MUONRecoCheck.C
1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 *                                                                        *
4 * Author: The ALICE Off-line Project.                                    *
5 * Contributors are mentioned in the code where appropriate.              *
6 *                                                                        *
7 * Permission to use, copy, modify and distribute this software and its   *
8 * documentation strictly for non-commercial purposes is hereby granted   *
9 * without fee, provided that the above copyright notice appears in all   *
10 * copies and that both the copyright notice and this permission notice   *
11 * appear in the supporting documentation. The authors make no claims     *
12 * about the suitability of this software for any purpose. It is          *
13 * provided "as is" without express or implied warranty.                  *
14 **************************************************************************/
15
16 // $Id$
17
18 /// \ingroup macros
19 /// \file MUONRecoCheck.C
20 /// \brief Utility macro to check the muon reconstruction. 
21 ///
22 /// Reconstructed tracks are compared to reference tracks. The reference tracks 
23 /// are built from AliTrackReference for the hit in chamber (0..9) and from 
24 /// kinematics (TreeK) for the vertex parameters.  
25 ///
26 /// \author Jean-Pierre Cussonneau, Philippe Pillot, Subatech  
27
28 // ROOT includes
29 #include <Riostream.h>
30 #include "TMath.h"
31 #include "TObjArray.h"
32 #include "TH1.h"
33 #include "TH2.h"
34 #include "TH3.h"
35 #include "TGraphErrors.h"
36 #include "TGraphAsymmErrors.h"
37 #include "TF1.h"
38 #include "TFile.h"
39 #include "TCanvas.h"
40 #include "TLegend.h"
41 #include "TGeoManager.h"
42
43 // STEER includes
44 #include "AliCDBManager.h"
45 #include "AliGeomManager.h"
46 #include "AliLog.h"
47
48 // MUON includes
49 #include "AliMUONCDB.h"
50 #include "AliMUONConstants.h"
51 #include "AliMUONTrack.h"
52 #include "AliMUONRecoCheck.h"
53 #include "AliMUONTrackParam.h"
54 #include "AliMUONRecoParam.h"
55 #include "AliMUONVTrackStore.h"
56 #include "AliMUONVCluster.h"
57 #include "AliMUONTrackExtrap.h"
58 #include "AliMUONESDInterface.h"
59 #include "AliMUONVTriggerTrackStore.h"
60 #include "AliMUONTriggerTrack.h"
61
62 Double_t langaufun(Double_t *x, Double_t *par);
63 void     FitGausResVsMom(TH2* h, Int_t nBins, const Double_t mean0, const Double_t sigma0, const char* fitting, TGraphAsymmErrors* gMean, TGraphAsymmErrors* gSigma);
64 void     FitPDCAVsMom(TH2* h, Int_t nBins, const char* fitting, TGraphAsymmErrors* gMean, TGraphAsymmErrors* gSigma);
65 TCanvas* DrawVsAng(const char* name, const char* title, TH1* h1, TH2* h2);
66 TCanvas* DrawVsPos(const char* name, const char* title, TH2* h1, TH2* h2, TH2* h3);
67 TCanvas* DrawResMomVsMom(const char* name, const char* title, TH2* h, Int_t nBins, TF1* f2 = 0x0, const char* fitting = "");
68
69 //------------------------------------------------------------------------------------
70 void MUONRecoCheck(Int_t nEvent = -1, const char* pathSim="./generated/", const char* esdFileName="AliESDs.root",
71                    const char* ocdbPath = "local://$ALICE_ROOT/OCDB", const Double_t pMin = 0., const Double_t pMax = 300.,
72                    const Int_t pNBins = 30, Int_t absorberRegion = -1)
73 {
74   /// Associate the reconstructed tracks with the simulated ones and check the quality of the reconstruction
75   /// (tracking/trigger efficiency; momentum, slope,... resolutions at first cluster and at vertex; cluster resolution).
76   /// - You can choose the momentum range and number of bins used to study the track resolution versus momentum.
77   /// - You can limit the calculation of track resolution at vertex to the tracks crossing the absorber in a given region
78   /// with the flag "absorberRegion": -1=all, 1=[2,3]deg, 2=[3,10]deg.
79   
80   Double_t aAbsLimits[2];
81   if (absorberRegion > -1) {
82     if (absorberRegion == 1) {
83       aAbsLimits[0] = 2.;
84       aAbsLimits[1] = 3.;
85     } else if (absorberRegion == 2) {
86       aAbsLimits[0] = 3.;
87       aAbsLimits[1] = 10.;
88     } else {
89       cout<<"Unknown absorber region. Valid choices are: -1=all, 1=[2,3]deg, 2=[3,10]deg"<<endl;
90       return;
91     }
92   } else {
93     aAbsLimits[0] = 0.;
94     aAbsLimits[1] = 90.;
95   }
96   
97   cout<<endl<<"Momentum range for track resolution studies: "<<pNBins<<" bins in ["<<pMin<<","<<pMax<<"] GeV/c"<<endl;
98   if (pMin >= pMax || pNBins <= 0) {
99     cout<<"--> incorrect momentum range"<<endl<<endl;
100     return;
101   } else cout<<endl;
102   
103   AliLog::SetClassDebugLevel("AliMCEvent",-1);
104   
105   // ###################################### define histograms ###################################### //
106   // File for histograms and histogram booking
107   TFile *histoFile = new TFile("MUONRecoCheck.root", "RECREATE");
108   
109   TH1F *hReconstructible = new TH1F("hReconstructible"," Nb of reconstructible tracks / evt",15,-0.5,14.5);
110   TH1F *hReco = new TH1F("hReco"," Nb of reconstructed tracks / evt",15,-0.5,14.5);
111   TH1F *hNClusterComp = new TH1F("hNClusterComp"," Nb of compatible clusters / track ",15,-0.5,14.5);
112   TH1F *hTrackRefID = new TH1F("hTrackRefID"," track reference ID ",100,-0.5,99.5);
113   TH1F *hTriggerable = new TH1F("hTriggerable"," Nb of triggerable tracks / evt",15,-0.5,14.5);
114   TH1F *hTriggered = new TH1F("hTriggered"," Nb of triggered tracks / evt",15,-0.5,14.5);
115   
116   // momentum resolution at vertex
117   histoFile->mkdir("momentumAtVertex","momentumAtVertex");
118   histoFile->cd("momentumAtVertex");
119   
120   const Double_t pEdges[2] = {pMin, pMax};
121   const Int_t deltaPAtVtxNBins = 250;
122   Double_t deltaPAtVtxEdges[2];
123   if (pMax < 50.) { deltaPAtVtxEdges[0] = -20.; deltaPAtVtxEdges[1] = 5.; }
124   else { deltaPAtVtxEdges[0] = -35.; deltaPAtVtxEdges[1] = 15.; }
125   
126   TH1F *hResMomVertex = new TH1F("hResMomVertex"," delta P at vertex;#Delta_{p} (GeV/c)",deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
127   
128   TH2D *hResMomVertexVsMom = new TH2D("hResMomVertexVsMom","#Delta_{p} at vertex versus p;p (GeV/c);#Delta_{p} (GeV/c)",2*pNBins,pEdges[0],pEdges[1],deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
129   TH2D *hResMomVertexVsMom_2_3_Deg = new TH2D("hResMomVertexVsMom_2_3_Deg","#Delta_{p} at vertex versus p for tracks between 2 and 3 degrees at absorber end;p (GeV/c);#Delta_{p} (GeV/c)",2*pNBins,pEdges[0],pEdges[1],deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
130   TH2D *hResMomVertexVsMom_3_10_Deg = new TH2D("hResMomVertexVsMom_3_10_Deg","#Delta_{p} at vertex versus p for tracks between 3 and 10 degrees at absorber end;p (GeV/c);#Delta_{p} (GeV/c)",2*pNBins,pEdges[0],pEdges[1],deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
131   TH2D *hResMomVertexVsMom_0_2_DegMC = new TH2D("hResMomVertexVsMom_0_2_DegMC","#Delta_{p} at vertex versus p for tracks with MC angle below 2 degrees;p (GeV/c);#Delta_{p} (GeV/c)",2*pNBins,pEdges[0],pEdges[1],deltaPAtVtxNBins/10,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
132   
133   TH2D *hResMomVertexVsPosAbsEnd_0_2_DegMC = new TH2D("hResMomVertexVsPosAbsEnd_0_2_DegMC","#Delta_{p} at vertex versus track position at absorber end for tracks with MC angle < 2 degrees;position (cm);#Delta_{p} (GeV/c)",1000,0.,100.,deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
134   TH2D *hResMomVertexVsPosAbsEnd_2_3_DegMC = new TH2D("hResMomVertexVsPosAbsEnd_2_3_DegMC","#Delta_{p} at vertex versus track position at absorber end for tracks with MC angle in [2,3[ degrees;position (cm);#Delta_{p} (GeV/c)",1000,0.,100.,deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
135   TH2D *hResMomVertexVsPosAbsEnd_3_10_DegMC = new TH2D("hResMomVertexVsPosAbsEnd_3_10_DegMC","#Delta_{p} at vertex versus track position at absorber end for tracks with MC angle in [3,10[ degrees;position (cm);#Delta_{p} (GeV/c)",1000,0.,100.,deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
136   
137   TH2D *hResMomVertexVsAngle = new TH2D("hResMomVertexVsAngle","#Delta_{p} at vertex versus track position at absorber end converted to degrees;angle (Deg);#Delta_{p} (GeV/c)",10,0.,10.,deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
138   TH2D *hResMomVertexVsMCAngle = new TH2D("hResMomVertexVsMCAngle","#Delta_{p} at vertex versus MC angle;MC angle (Deg);#Delta_{p} (GeV/c)",10,0.,10.,deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
139   TH3D *hResMomVertexVsAngleVsMom = new TH3D("hResMomVertexVsAngleVsMom","#Delta_{p} at vertex versus track position at absorber end converted to degrees versus momentum;p (GeV/c);angle (Deg);#Delta_{p} (GeV/c)",2*pNBins,pEdges[0],pEdges[1],100,0.,10.,deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
140   
141   TGraphAsymmErrors* gMeanResMomVertexVsMom = new TGraphAsymmErrors(pNBins);
142   gMeanResMomVertexVsMom->SetName("gMeanResMomVertexVsMom");
143   gMeanResMomVertexVsMom->SetTitle("<#Delta_{p}> at vertex versus p;p (GeV/c);<#Delta_{p}> (GeV/c)");
144   TGraphAsymmErrors* gMostProbResMomVertexVsMom = new TGraphAsymmErrors(pNBins);
145   gMostProbResMomVertexVsMom->SetName("gMostProbResMomVertexVsMom");
146   gMostProbResMomVertexVsMom->SetTitle("Most probable #Delta_{p} at vertex versus p;p (GeV/c);Most prob. #Delta_{p} (GeV/c)");
147   TGraphAsymmErrors* gSigmaResMomVertexVsMom = new TGraphAsymmErrors(pNBins);
148   gSigmaResMomVertexVsMom->SetName("gSigmaResMomVertexVsMom");
149   gSigmaResMomVertexVsMom->SetTitle("#sigma_{p}/p at vertex versus p;p (GeV/c);#sigma_{p}/p (%)");
150   
151   // momentum resolution at first cluster
152   histoFile->mkdir("momentumAtFirstCluster","momentumAtFirstCluster");
153   histoFile->cd("momentumAtFirstCluster");
154   
155   const Int_t deltaPAtFirstClNBins = 500;
156   Double_t deltaPAtFirstClEdges[2];
157   if (pMax < 25.) { deltaPAtFirstClEdges[0] = -5.; deltaPAtFirstClEdges[1] = 5.; }
158   else if (pMax < 50.) { deltaPAtFirstClEdges[0] = -10.; deltaPAtFirstClEdges[1] = 10.; }
159   else { deltaPAtFirstClEdges[0] = -25.; deltaPAtFirstClEdges[1] = 25.; }
160   
161   TH1F *hResMomFirstCluster = new TH1F("hResMomFirstCluster"," delta P at first cluster;#Delta_{p} (GeV/c)",deltaPAtFirstClNBins,deltaPAtFirstClEdges[0],deltaPAtFirstClEdges[1]);
162   TH2D *hResMomFirstClusterVsMom = new TH2D("hResMomFirstClusterVsMom","#Delta_{p} at first cluster versus p;p (GeV/c);#Delta_{p} (GeV/c)",2*pNBins,pEdges[0],pEdges[1],deltaPAtFirstClNBins,deltaPAtFirstClEdges[0],deltaPAtFirstClEdges[1]);
163   
164   TGraphAsymmErrors* gMeanResMomFirstClusterVsMom = new TGraphAsymmErrors(pNBins);
165   gMeanResMomFirstClusterVsMom->SetName("gMeanResMomFirstClusterVsMom");
166   gMeanResMomFirstClusterVsMom->SetTitle("<#Delta_{p}> at first cluster versus p;p (GeV/c);<#Delta_{p}> (GeV/c)");
167   TGraphAsymmErrors* gSigmaResMomFirstClusterVsMom = new TGraphAsymmErrors(pNBins);
168   gSigmaResMomFirstClusterVsMom->SetName("gSigmaResMomFirstClusterVsMom");
169   gSigmaResMomFirstClusterVsMom->SetTitle("#sigma_{p}/p at first cluster versus p;p (GeV/c);#sigma_{p}/p (%)");
170   
171   // angular resolution at vertex
172   histoFile->mkdir("slopesAtVertex","slopesAtVertex");
173   histoFile->cd("slopesAtVertex");
174   
175   const Int_t deltaSlopeAtVtxNBins = 500;
176   const Double_t deltaSlopeAtVtxEdges[2] = {-0.05, 0.05};
177   
178   TH1F *hResSlopeXVertex = new TH1F("hResSlopeXVertex","#Delta_{slope_{X}} at vertex;#Delta_{slope_{X}}", deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
179   TH1F *hResSlopeYVertex = new TH1F("hResSlopeYVertex","#Delta_{slope_{Y}} at vertex;#Delta_{slope_{Y}}", deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
180   TH2D *hResSlopeXVertexVsMom = new TH2D("hResSlopeXVertexVsMom","#Delta_{slope_{X}} at vertex versus p;p (GeV/c);#Delta_{slope_{X}}",2*pNBins,pEdges[0],pEdges[1], deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
181   TH2D *hResSlopeYVertexVsMom = new TH2D("hResSlopeYVertexVsMom","#Delta_{slope_{Y}} at vertex versus p;p (GeV/c);#Delta_{slope_{Y}}",2*pNBins,pEdges[0],pEdges[1], deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
182   
183   TH2D *hResSlopeXVertexVsPosAbsEnd_0_2_DegMC = new TH2D("hResSlopeXVertexVsPosAbsEnd_0_2_DegMC","#Delta_{slope_{X}} at vertex versus track position at absorber end for tracks with MC angle < 2 degrees;position (cm);#Delta_{slope_{X}}",1000,0.,100.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
184   TH2D *hResSlopeYVertexVsPosAbsEnd_0_2_DegMC = new TH2D("hResSlopeYVertexVsPosAbsEnd_0_2_DegMC","#Delta_{slope_{Y}} at vertex versus track position at absorber end for tracks with MC angle < 2 degrees;position (cm);#Delta_{slope_{Y}}",1000,0.,100.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
185   TH2D *hResSlopeXVertexVsPosAbsEnd_2_3_DegMC = new TH2D("hResSlopeXVertexVsPosAbsEnd_2_3_DegMC","#Delta_{slope_{X}} at vertex versus track position at absorber end for tracks with MC angle in [2,3[ degrees;position (cm);#Delta_{slope_{X}}",1000,0.,100.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
186   TH2D *hResSlopeYVertexVsPosAbsEnd_2_3_DegMC = new TH2D("hResSlopeYVertexVsPosAbsEnd_2_3_DegMC","#Delta_{slope_{Y}} at vertex versus track position at absorber end for tracks with MC angle in [2,3[ degrees;position (cm);#Delta_{slope_{Y}}",1000,0.,100.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
187   TH2D *hResSlopeXVertexVsPosAbsEnd_3_10_DegMC = new TH2D("hResSlopeXVertexVsPosAbsEnd_3_10_DegMC","#Delta_{slope_{X}} at vertex versus track position at absorber end for tracks with MC angle in [3,10[ degrees;position (cm);#Delta_{slope_{X}}",1000,0.,100.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
188   TH2D *hResSlopeYVertexVsPosAbsEnd_3_10_DegMC = new TH2D("hResSlopeYVertexVsPosAbsEnd_3_10_DegMC","#Delta_{slope_{Y}} at vertex versus track position at absorber end for tracks with MC angle in [3,10[ degrees;position (cm);#Delta_{slope_{Y}}",1000,0.,100.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
189   
190   TH2D *hResSlopeXVertexVsAngle = new TH2D("hResSlopeXVertexVsAngle","#Delta_{slope_{X}} at vertex versus track position at absorber end converted to degrees;angle (Deg);#Delta_{slope_{X}}",10,0.,10.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
191   TH2D *hResSlopeYVertexVsAngle = new TH2D("hResSlopeYVertexVsAngle","#Delta_{slope_{Y}} at vertex versus track position at absorber end converted to degrees;angle (Deg);#Delta_{slope_{Y}}",10,0.,10.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
192   TH2D *hResSlopeXVertexVsMCAngle = new TH2D("hResSlopeXVertexVsMCAngle","#Delta_{slope_{X}} at vertex versus MC angle;MC angle (Deg);#Delta_{slope_{X}}",10,0.,10.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
193   TH2D *hResSlopeYVertexVsMCAngle = new TH2D("hResSlopeYVertexVsMCAngle","#Delta_{slope_{Y}} at vertex versus MC angle;MC angle (Deg);#Delta_{slope_{Y}}",10,0.,10.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
194   
195   TGraphAsymmErrors* gMeanResSlopeXVertexVsMom = new TGraphAsymmErrors(pNBins);
196   gMeanResSlopeXVertexVsMom->SetName("gMeanResSlopeXVertexVsMom");
197   gMeanResSlopeXVertexVsMom->SetTitle("<#Delta_{slope_{X}}> at vertex versus p;p (GeV/c);<#Delta_{slope_{X}}>");
198   TGraphAsymmErrors* gSigmaResSlopeXVertexVsMom = new TGraphAsymmErrors(pNBins);
199   gSigmaResSlopeXVertexVsMom->SetName("gSigmaResSlopeXVertexVsMom");
200   gSigmaResSlopeXVertexVsMom->SetTitle("#sigma_{slope_{X}} at vertex versus p;p (GeV/c);#sigma_{slope_{X}}");
201   TGraphAsymmErrors* gMeanResSlopeYVertexVsMom = new TGraphAsymmErrors(pNBins);
202   gMeanResSlopeYVertexVsMom->SetName("gMeanResSlopeYVertexVsMom");
203   gMeanResSlopeYVertexVsMom->SetTitle("<#Delta_{slope_{Y}}> at vertex versus p;p (GeV/c);<#Delta_{slope_{Y}}>");
204   TGraphAsymmErrors* gSigmaResSlopeYVertexVsMom = new TGraphAsymmErrors(pNBins);
205   gSigmaResSlopeYVertexVsMom->SetName("gSigmaResSlopeYVertexVsMom");
206   gSigmaResSlopeYVertexVsMom->SetTitle("#sigma_{slope_{Y}} at vertex versus p;p (GeV/c);#sigma_{slope_{Y}}");
207   
208   // angular resolution at first cluster
209   histoFile->mkdir("slopesAtFirstCluster","slopesAtFirstCluster");
210   histoFile->cd("slopesAtFirstCluster");
211   
212   const Int_t deltaSlopeAtFirstClNBins = 500;
213   const Double_t deltaSlopeAtFirstClEdges[2] = {-0.01, 0.01};
214   
215   TH1F *hResSlopeXFirstCluster = new TH1F("hResSlopeXFirstCluster","#Delta_{slope_{X}} at first cluster;#Delta_{slope_{X}}", deltaSlopeAtFirstClNBins, deltaSlopeAtFirstClEdges[0], deltaSlopeAtFirstClEdges[1]);
216   TH2D *hResSlopeXFirstClusterVsMom = new TH2D("hResSlopeXFirstClusterVsMom","#Delta_{slope_{X}} at first cluster versus p;p (GeV/c);#Delta_{slope_{X}}",2*pNBins,pEdges[0],pEdges[1], deltaSlopeAtFirstClNBins, deltaSlopeAtFirstClEdges[0], deltaSlopeAtFirstClEdges[1]);
217   TH1F *hResSlopeYFirstCluster = new TH1F("hResSlopeYFirstCluster","#Delta_{slope_{Y}} at first cluster;#Delta_{slope_{Y}}", deltaSlopeAtFirstClNBins, deltaSlopeAtFirstClEdges[0], deltaSlopeAtFirstClEdges[1]);
218   TH2D *hResSlopeYFirstClusterVsMom = new TH2D("hResSlopeYFirstClusterVsMom","#Delta_{slope_{Y}} at first cluster versus p;p (GeV/c);#Delta_{slope_{Y}}",2*pNBins,pEdges[0],pEdges[1], deltaSlopeAtFirstClNBins, deltaSlopeAtFirstClEdges[0], deltaSlopeAtFirstClEdges[1]);
219   
220   TGraphAsymmErrors* gMeanResSlopeXFirstClusterVsMom = new TGraphAsymmErrors(pNBins);
221   gMeanResSlopeXFirstClusterVsMom->SetName("gMeanResSlopeXFirstClusterVsMom");
222   gMeanResSlopeXFirstClusterVsMom->SetTitle("<#Delta_{slope_{X}}> at first cluster versus p;p (GeV/c);<#Delta_{slope_{X}}>");
223   TGraphAsymmErrors* gSigmaResSlopeXFirstClusterVsMom = new TGraphAsymmErrors(pNBins);
224   gSigmaResSlopeXFirstClusterVsMom->SetName("gSigmaResSlopeXFirstClusterVsMom");
225   gSigmaResSlopeXFirstClusterVsMom->SetTitle("#sigma_{slope_{X}} at first cluster versus p;p (GeV/c);#sigma_{slope_{X}}");
226   TGraphAsymmErrors* gMeanResSlopeYFirstClusterVsMom = new TGraphAsymmErrors(pNBins);
227   gMeanResSlopeYFirstClusterVsMom->SetName("gMeanResSlopeYFirstClusterVsMom");
228   gMeanResSlopeYFirstClusterVsMom->SetTitle("<#Delta_{slope_{Y}}> at first cluster versus p;p (GeV/c);<#Delta_{slope_{Y}}>");
229   TGraphAsymmErrors* gSigmaResSlopeYFirstClusterVsMom = new TGraphAsymmErrors(pNBins);
230   gSigmaResSlopeYFirstClusterVsMom->SetName("gSigmaResSlopeYFirstClusterVsMom");
231   gSigmaResSlopeYFirstClusterVsMom->SetTitle("#sigma_{slope_{Y}} at first cluster versus p;p (GeV/c);#sigma_{slope_{Y}}");
232   
233   // DCA resolution and MCS angular dispersion
234   histoFile->mkdir("DCA","DCA");
235   histoFile->cd("DCA");
236   
237   const Int_t deltaPDCANBins = 500;
238   const Double_t deltaPDCAEdges[2] = {0., 1000.};
239   const Double_t deltaPMCSAngEdges[2] = {-1.5, 1.5};
240   
241   TH1F *hPDCA = new TH1F("hPDCA","p #times DCA at vertex;p #times DCA (GeV #times cm)", deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
242   TH2D *hPDCAVsMom_2_3_Deg = new TH2D("hPDCAVsMom_2_3_Deg","p #times DCA versus p for tracks within [2,3[ degrees at absorber end;p (GeV/c);p #times DCA (GeV #times cm)",2*pNBins,pEdges[0],pEdges[1], deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
243   TH2D *hPDCAVsMom_3_10_Deg = new TH2D("hPDCAVsMom_3_10_Deg","p #times DCA versus p for tracks within [3,10[ degrees at absorber end;p (GeV/c);p #times DCA (GeV #times cm)",2*pNBins,pEdges[0],pEdges[1], deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
244   TH2D *hPMCSAngVsMom_2_3_Deg = new TH2D("hPMCSAngVsMom_2_3_Deg","p #times #Delta#theta_{MCS} versus p for tracks within [2,3[ degrees at absorber end;p (GeV/c);p #times #Delta#theta_{MCS} (GeV)",2*pNBins,pEdges[0],pEdges[1], deltaPDCANBins, deltaPMCSAngEdges[0], deltaPMCSAngEdges[1]);
245   TH2D *hPMCSAngVsMom_3_10_Deg = new TH2D("hPMCSAngVsMom_3_10_Deg","p #times #Delta#theta_{MCS} versus p for tracks within [2,3[ degrees at absorber end;p (GeV/c);p #times #Delta#theta_{MCS} (GeV)",2*pNBins,pEdges[0],pEdges[1], deltaPDCANBins, deltaPMCSAngEdges[0], deltaPMCSAngEdges[1]);
246   
247   TH2D *hPDCAVsPosAbsEnd_0_2_DegMC = new TH2D("hPDCAVsPosAbsEnd_0_2_DegMC","p #times DCA versus track position at absorber end for tracks with MC angle < 2 degrees;position (cm);p #times DCA (GeV #times cm)",1000,0.,100.,deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
248   TH2D *hPDCAVsPosAbsEnd_2_3_DegMC = new TH2D("hPDCAVsPosAbsEnd_2_3_DegMC","p #times DCA}versus track position at absorber end for tracks with MC angle in [2,3[ degrees;position (cm);p #times DCA (GeV #times cm)",1000,0.,100.,deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
249   TH2D *hPDCAVsPosAbsEnd_3_10_DegMC = new TH2D("hPDCAVsPosAbsEnd_3_10_DegMC","p #times DCA versus track position at absorber end for tracks with MC angle in [3,10[ degrees;position (cm);p #times DCA (GeV #times cm)",1000,0.,100.,deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
250   
251   TH2D *hPDCAVsAngle = new TH2D("hPDCAVsAngle","p #times DCA versus track position at absorber end converted to degrees;angle (Deg);p #times DCA (GeV #times cm)",10,0.,10.,deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
252   TH2D *hPDCAVsMCAngle = new TH2D("hPDCAVsMCAngle","p #times DCA versus MC angle;MC angle (Deg);p #times DCA (GeV #times cm)",10,0.,10.,deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
253   
254   TGraphAsymmErrors* gMeanPDCAVsMom_2_3_Deg = new TGraphAsymmErrors(pNBins);
255   gMeanPDCAVsMom_2_3_Deg->SetName("gMeanPDCAVsMom_2_3_Deg");
256   gMeanPDCAVsMom_2_3_Deg->SetTitle("<p #times DCA> versus p for tracks within [2,3[ degrees at absorber end;p (GeV/c);<p #times DCA> (GeV #times cm)");
257   TGraphAsymmErrors* gSigmaPDCAVsMom_2_3_Deg = new TGraphAsymmErrors(pNBins);
258   gSigmaPDCAVsMom_2_3_Deg->SetName("gSigmaPDCAVsMom_2_3_Deg");
259   gSigmaPDCAVsMom_2_3_Deg->SetTitle("#sigma_{p #times DCA} versus p for tracks within [2,3[ degrees at absorber end;p (GeV/c);#sigma_{p #times DCA} (GeV #times cm)");
260   TGraphAsymmErrors* gMeanPDCAVsMom_3_10_Deg = new TGraphAsymmErrors(pNBins);
261   gMeanPDCAVsMom_3_10_Deg->SetName("gMeanPDCAVsMom_3_10_Deg");
262   gMeanPDCAVsMom_3_10_Deg->SetTitle("<p #times DCA> versus p for tracks within [3,10[ degrees at absorber end;p (GeV/c);<p #times DCA> (GeV #times cm)");
263   TGraphAsymmErrors* gSigmaPDCAVsMom_3_10_Deg = new TGraphAsymmErrors(pNBins);
264   gSigmaPDCAVsMom_3_10_Deg->SetName("gSigmaPDCAVsMom_3_10_Deg");
265   gSigmaPDCAVsMom_3_10_Deg->SetTitle("#sigma_{p #times DCA} versus p for tracks within [3,10[ degrees at absorber end;p (GeV/c);#sigma_{p #times DCA} (GeV #times cm)");
266   TGraphAsymmErrors* gMeanPMCSAngVsMom_2_3_Deg = new TGraphAsymmErrors(pNBins);
267   gMeanPMCSAngVsMom_2_3_Deg->SetName("gMeanPMCSAngVsMom_2_3_Deg");
268   gMeanPMCSAngVsMom_2_3_Deg->SetTitle("<p #times #Delta#theta_{MCS}> versus p for tracks within [2,3[ degrees at absorber end;p (GeV/c);<p #times #Delta#theta_{MCS}> (GeV)");
269   TGraphAsymmErrors* gSigmaPMCSAngVsMom_2_3_Deg = new TGraphAsymmErrors(pNBins);
270   gSigmaPMCSAngVsMom_2_3_Deg->SetName("gSigmaPMCSAngVsMom_2_3_Deg");
271   gSigmaPMCSAngVsMom_2_3_Deg->SetTitle("#sigma_{p #times #Delta#theta_{MCS}} versus p for tracks within [2,3[ degrees at absorber end;p (GeV/c);#sigma_{p #times #Delta#theta_{MCS}} (GeV)");
272   TGraphAsymmErrors* gMeanPMCSAngVsMom_3_10_Deg = new TGraphAsymmErrors(pNBins);
273   gMeanPMCSAngVsMom_3_10_Deg->SetName("gMeanPMCSAngVsMom_3_10_Deg");
274   gMeanPMCSAngVsMom_3_10_Deg->SetTitle("<p #times #Delta#theta_{MCS}> versus p for tracks within [3,10[ degrees at absorber end;p (GeV/c);<p #times #Delta#theta_{MCS}> (GeV)");
275   TGraphAsymmErrors* gSigmaPMCSAngVsMom_3_10_Deg = new TGraphAsymmErrors(pNBins);
276   gSigmaPMCSAngVsMom_3_10_Deg->SetName("gSigmaPMCSAngVsMom_3_10_Deg");
277   gSigmaPMCSAngVsMom_3_10_Deg->SetTitle("#sigma_{p #times #Delta#theta_{MCS}} versus p for tracks within [3,10[ degrees at absorber end;p (GeV/c);#sigma_{p #times #Delta#theta_{MCS}} (GeV)");
278   
279   // eta resolution at vertex
280   histoFile->mkdir("etaAtVertex","etaAtVertex");
281   histoFile->cd("etaAtVertex");
282   
283   const Int_t deltaEtaAtVtxNBins = 500;
284   const Double_t deltaEtaAtVtxEdges[2] = {-0.5, 0.5};
285   
286   TH1F *hResEtaVertex = new TH1F("hResEtaVertex","#Delta_{eta} at vertex;#Delta_{eta}", deltaEtaAtVtxNBins, deltaEtaAtVtxEdges[0], deltaEtaAtVtxEdges[1]);
287   TH2D *hResEtaVertexVsMom = new TH2D("hResEtaVertexVsMom","#Delta_{eta} at vertex versus p;p (GeV/c);#Delta_{eta}",2*pNBins,pEdges[0],pEdges[1], deltaEtaAtVtxNBins, deltaEtaAtVtxEdges[0], deltaEtaAtVtxEdges[1]);
288   
289   TH2D *hResEtaVertexVsPosAbsEnd_0_2_DegMC = new TH2D("hResEtaVertexVsPosAbsEnd_0_2_DegMC","#Delta_{eta} at vertex versus track position at absorber end for tracks with MC angle < 2 degrees;position (cm);#Delta_{eta}",1000,0.,100.,deltaEtaAtVtxNBins, deltaEtaAtVtxEdges[0], deltaEtaAtVtxEdges[1]);
290   TH2D *hResEtaVertexVsPosAbsEnd_2_3_DegMC = new TH2D("hResEtaVertexVsPosAbsEnd_2_3_DegMC","#Delta_{eta} at vertex versus track position at absorber end for tracks with MC angle in [2,3[ degrees;position (cm);#Delta_{eta}",1000,0.,100.,deltaEtaAtVtxNBins, deltaEtaAtVtxEdges[0], deltaEtaAtVtxEdges[1]);
291   TH2D *hResEtaVertexVsPosAbsEnd_3_10_DegMC = new TH2D("hResEtaVertexVsPosAbsEnd_3_10_DegMC","#Delta_{eta} at vertex versus track position at absorber end for tracks with MC angle in [3,10[ degrees;position (cm);#Delta_{eta}",1000,0.,100.,deltaEtaAtVtxNBins, deltaEtaAtVtxEdges[0], deltaEtaAtVtxEdges[1]);
292   
293   TH2D *hResEtaVertexVsAngle = new TH2D("hResEtaVertexVsAngle","#Delta_{eta} at vertex versus track position at absorber end converted to degrees;angle (Deg);#Delta_{eta}",10,0.,10.,deltaEtaAtVtxNBins, deltaEtaAtVtxEdges[0], deltaEtaAtVtxEdges[1]);
294   TH2D *hResEtaVertexVsMCAngle = new TH2D("hResEtaVertexVsMCAngle","#Delta_{eta} at vertex versus MC angle;MC angle (Deg);#Delta_{eta}",10,0.,10.,deltaEtaAtVtxNBins, deltaEtaAtVtxEdges[0], deltaEtaAtVtxEdges[1]);
295   
296   TGraphAsymmErrors* gMeanResEtaVertexVsMom = new TGraphAsymmErrors(pNBins);
297   gMeanResEtaVertexVsMom->SetName("gMeanResEtaVertexVsMom");
298   gMeanResEtaVertexVsMom->SetTitle("<#Delta_{eta}> at vertex versus p;p (GeV/c);<#Delta_{eta}>");
299   TGraphAsymmErrors* gSigmaResEtaVertexVsMom = new TGraphAsymmErrors(pNBins);
300   gSigmaResEtaVertexVsMom->SetName("gSigmaResEtaVertexVsMom");
301   gSigmaResEtaVertexVsMom->SetTitle("#sigma_{eta} at vertex versus p;p (GeV/c);#sigma_{eta}");
302   
303   // phi resolution at vertex
304   histoFile->mkdir("phiAtVertex","phiAtVertex");
305   histoFile->cd("phiAtVertex");
306   
307   const Int_t deltaPhiAtVtxNBins = 500;
308   const Double_t deltaPhiAtVtxEdges[2] = {-0.5, 0.5};
309   
310   TH1F *hResPhiVertex = new TH1F("hResPhiVertex","#Delta_{phi} at vertex;#Delta_{phi}", deltaPhiAtVtxNBins, deltaPhiAtVtxEdges[0], deltaPhiAtVtxEdges[1]);
311   TH2D *hResPhiVertexVsMom = new TH2D("hResPhiVertexVsMom","#Delta_{phi} at vertex versus p;p (GeV/c);#Delta_{phi}",2*pNBins,pEdges[0],pEdges[1], deltaPhiAtVtxNBins, deltaPhiAtVtxEdges[0], deltaPhiAtVtxEdges[1]);
312   
313   TH2D *hResPhiVertexVsPosAbsEnd_0_2_DegMC = new TH2D("hResPhiVertexVsPosAbsEnd_0_2_DegMC","#Delta_{phi} at vertex versus track position at absorber end for tracks with MC angle < 2 degrees;position (cm);#Delta_{phi}",1000,0.,100.,deltaPhiAtVtxNBins, deltaPhiAtVtxEdges[0], deltaPhiAtVtxEdges[1]);
314   TH2D *hResPhiVertexVsPosAbsEnd_2_3_DegMC = new TH2D("hResPhiVertexVsPosAbsEnd_2_3_DegMC","#Delta_{phi} at vertex versus track position at absorber end for tracks with MC angle in [2,3[ degrees;position (cm);#Delta_{phi}",1000,0.,100.,deltaPhiAtVtxNBins, deltaPhiAtVtxEdges[0], deltaPhiAtVtxEdges[1]);
315   TH2D *hResPhiVertexVsPosAbsEnd_3_10_DegMC = new TH2D("hResPhiVertexVsPosAbsEnd_3_10_DegMC","#Delta_{phi} at vertex versus track position at absorber end for tracks with MC angle in [3,10[ degrees;position (cm);#Delta_{phi}",1000,0.,100.,deltaPhiAtVtxNBins, deltaPhiAtVtxEdges[0], deltaPhiAtVtxEdges[1]);
316   
317   TH2D *hResPhiVertexVsAngle = new TH2D("hResPhiVertexVsAngle","#Delta_{phi} at vertex versus track position at absorber end converted to degrees;angle (Deg);#Delta_{phi}",10,0.,10.,deltaPhiAtVtxNBins, deltaPhiAtVtxEdges[0], deltaPhiAtVtxEdges[1]);
318   TH2D *hResPhiVertexVsMCAngle = new TH2D("hResPhiVertexVsMCAngle","#Delta_{phi} at vertex versus MC angle;MC angle (Deg);#Delta_{phi}",10,0.,10.,deltaPhiAtVtxNBins, deltaPhiAtVtxEdges[0], deltaPhiAtVtxEdges[1]);
319   
320   TGraphAsymmErrors* gMeanResPhiVertexVsMom = new TGraphAsymmErrors(pNBins);
321   gMeanResPhiVertexVsMom->SetName("gMeanResPhiVertexVsMom");
322   gMeanResPhiVertexVsMom->SetTitle("<#Delta_{phi}> at vertex versus p;p (GeV/c);<#Delta_{phi}>");
323   TGraphAsymmErrors* gSigmaResPhiVertexVsMom = new TGraphAsymmErrors(pNBins);
324   gSigmaResPhiVertexVsMom->SetName("gSigmaResPhiVertexVsMom");
325   gSigmaResPhiVertexVsMom->SetTitle("#sigma_{phi} at vertex versus p;p (GeV/c);#sigma_{phi}");
326   
327   // cluster resolution
328   histoFile->mkdir("clusters","clusters");
329   histoFile->cd("clusters");
330   
331   TH1F* hResidualXInCh[AliMUONConstants::NTrackingCh()];
332   TH1F* hResidualYInCh[AliMUONConstants::NTrackingCh()];
333   for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
334     hResidualXInCh[i] = new TH1F(Form("hResidualXInCh%d",i+1), Form("cluster-track residual-X distribution in chamber %d;#Delta_{X} (cm)",i+1), 1000, -1., 1.);
335     hResidualYInCh[i] = new TH1F(Form("hResidualYInCh%d",i+1), Form("cluster-track residual-Y distribution in chamber %d;#Delta_{Y} (cm)",i+1), 1000, -0.5, 0.5);
336   }
337   
338   TGraphErrors* gResidualXPerChMean = new TGraphErrors(AliMUONConstants::NTrackingCh());
339   gResidualXPerChMean->SetName("gResidualXPerChMean");
340   gResidualXPerChMean->SetTitle("cluster-trackRef residual-X per Ch: mean;chamber ID;<#Delta_{X}> (cm)");
341   gResidualXPerChMean->SetMarkerStyle(kFullDotLarge);
342   TGraphErrors* gResidualYPerChMean = new TGraphErrors(AliMUONConstants::NTrackingCh());
343   gResidualYPerChMean->SetName("gResidualYPerChMean");
344   gResidualYPerChMean->SetTitle("cluster-trackRef residual-Y per Ch: mean;chamber ID;<#Delta_{Y}> (cm)");
345   gResidualYPerChMean->SetMarkerStyle(kFullDotLarge);
346   TGraphErrors* gResidualXPerChSigma = new TGraphErrors(AliMUONConstants::NTrackingCh());
347   gResidualXPerChSigma->SetName("gResidualXPerChSigma");
348   gResidualXPerChSigma->SetTitle("cluster-trackRef residual-X per Ch: sigma;chamber ID;#sigma_{X} (cm)");
349   gResidualXPerChSigma->SetMarkerStyle(kFullDotLarge);
350   TGraphErrors* gResidualYPerChSigma = new TGraphErrors(AliMUONConstants::NTrackingCh());
351   gResidualYPerChSigma->SetName("gResidualYPerChSigma");
352   gResidualYPerChSigma->SetTitle("cluster-trackRef residual-Y per Ch: sigma;chamber ID;#sigma_{Y} (cm)");
353   gResidualYPerChSigma->SetMarkerStyle(kFullDotLarge);
354
355   histoFile->mkdir("trigger");
356   histoFile->cd("trigger");
357   TH1F* hResidualTrigX11 = new TH1F("hResiudalTrigX11", "Residual X11", 100, -10., 10.);
358   TH1F* hResidualTrigY11 = new TH1F("hResiudalTrigY11", "Residual Y11", 100, -10., 10.);
359   TH1F* hResidualTrigSlopeY = new TH1F("hResiudalTrigSlopeY", "Residual Y slope", 100, -0.1, 0.1);
360   TH1F* hTriggerableMatchFailed = new TH1F("hTriggerableMatchFailed", "Triggerable multiplicity for events with no match", 15, -0.5, 14.5);
361   
362   
363   // ###################################### initialize ###################################### //
364   AliMUONRecoCheck rc(esdFileName, pathSim);
365   
366   // load necessary data from OCDB
367   AliCDBManager::Instance()->SetDefaultStorage(ocdbPath);
368   AliCDBManager::Instance()->SetRun(rc.GetRunNumber());
369   if (!AliMUONCDB::LoadField()) return;
370   AliMUONTrackExtrap::SetField();
371   AliGeomManager::LoadGeometry();
372   if (!AliGeomManager::GetGeometry()) return;
373   AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam();
374   if (!recoParam) return;
375   AliMUONESDInterface::ResetTracker(recoParam);
376   
377   // get sigma cut from recoParam to associate clusters with TrackRefs in case the label are not used
378   Double_t sigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking();
379   // compute the mask of requested stations from recoParam
380   UInt_t requestedStationMask = 0;
381   for (Int_t i = 0; i < 5; i++) if (recoParam->RequestStation(i)) requestedStationMask |= ( 1 << i );
382   // get from recoParam whether a track need 2 chambers hit in the same station (4 or 5) or not to be reconstructible
383   Bool_t request2ChInSameSt45 = !recoParam->MakeMoreTrackCandidates();
384   
385   Int_t nevents = rc.NumberOfEvents();
386   
387   if (nevents < nEvent || nEvent < 0) nEvent = nevents;
388   
389   Int_t ievent;
390   Int_t nReconstructibleTracks = 0;
391   Int_t nReconstructedTracks = 0;
392   Int_t nReconstructibleTracksCheck = 0;
393   AliMUONTrackParam *trackParam;
394   Double_t x1,y1,z1,slopex1,slopey1,pX1,pY1,pZ1,p1,pT1,eta1,phi1;
395   Double_t x2,y2,z2,slopex2,slopey2,pX2,pY2,pZ2,p2,pT2,eta2,phi2;
396   Double_t dPhi;
397   Double_t xAbs,yAbs,dAbs,aAbs,aMCS,aMC;
398   Double_t xDCA,yDCA,dca,pU;
399   Double_t aMCSMoy = 0., aMCS2Moy = 0., dMCSMoy = 0., dMCS2Moy = 0., adMCSMoy = 0.;
400   Int_t nMCS = 0;
401   
402   // ###################################### fill histograms ###################################### //
403   for (ievent=0; ievent<nEvent; ievent++)
404   {
405     if ((ievent+1)%100 == 0) cout<<"\rEvent processing... "<<ievent+1<<flush;
406     
407     AliMUONVTrackStore* trackStore = rc.ReconstructedTracks(ievent, kFALSE);
408     AliMUONVTrackStore* trackRefStore = rc.ReconstructibleTracks(ievent, requestedStationMask, request2ChInSameSt45);
409     
410     hReconstructible->Fill(trackRefStore->GetSize());
411     hReco->Fill(trackStore->GetSize());
412     
413     nReconstructibleTracks += trackRefStore->GetSize();
414     nReconstructedTracks += trackStore->GetSize();
415
416     AliMUONVTriggerTrackStore* triggerTrackRefStore = rc.TriggerableTracks(ievent);
417     AliMUONVTriggerTrackStore* triggerTrackStore = rc.TriggeredTracks(ievent);
418
419     hTriggerable->Fill(triggerTrackRefStore->GetSize());
420     hTriggered->Fill(triggerTrackStore->GetSize());
421
422     // loop over trigger trackRef
423     TIter nextTrig(triggerTrackRefStore->CreateIterator());
424     AliMUONTriggerTrack* triggerTrackRef;
425     Int_t nTriggerMatches = 0;
426     while ( ( triggerTrackRef = static_cast<AliMUONTriggerTrack*>(nextTrig()) ) )
427     {
428       
429       AliMUONTriggerTrack* triggerTrackMatched = 0x0;
430       
431       // loop over trackReco and look for compatible track
432       TIter nextTrig2(triggerTrackStore->CreateIterator());
433       AliMUONTriggerTrack* triggerTrackReco;
434       while ( ( triggerTrackReco = static_cast<AliMUONTriggerTrack*>(nextTrig2()) ) )
435       {
436         
437         // check if trackReco is compatible with trackRef
438         if (triggerTrackReco->Match(*triggerTrackRef, sigmaCut)) {
439           triggerTrackMatched = triggerTrackReco;
440           nTriggerMatches++;
441           break;
442         }
443       }
444       
445       if (triggerTrackMatched) { // tracking requirements verified, track is found
446         hResidualTrigX11->Fill( triggerTrackMatched->GetX11() - triggerTrackRef->GetX11() );
447         hResidualTrigY11->Fill( triggerTrackMatched->GetY11() - triggerTrackRef->GetY11() );
448         hResidualTrigSlopeY->Fill( triggerTrackMatched->GetSlopeY() - triggerTrackRef->GetSlopeY() );
449       }
450     } // loop on trigger track ref
451     
452     if ( nTriggerMatches != triggerTrackStore->GetSize() )
453       hTriggerableMatchFailed->Fill(triggerTrackRefStore->GetSize());
454     
455     // loop over trackRef
456     TIter next(trackRefStore->CreateIterator());
457     AliMUONTrack* trackRef;
458     while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) )
459     {
460       
461       hTrackRefID->Fill(trackRef->GetUniqueID());
462       
463       AliMUONTrack* trackMatched = 0x0;
464       Int_t nMatchClusters = 0;
465       
466       // loop over trackReco and look for compatible track
467       TIter next2(trackStore->CreateIterator());
468       AliMUONTrack* trackReco;
469       while ( ( trackReco = static_cast<AliMUONTrack*>(next2()) ) )
470       {
471         
472         // check if trackReco is compatible with trackRef
473         if (trackReco->Match(*trackRef, sigmaCut, nMatchClusters)) {
474           trackMatched = trackReco;
475           break;
476         }
477         
478       }
479       
480       if (trackMatched) { // tracking requirements verified, track is found
481         nReconstructibleTracksCheck++;
482         hNClusterComp->Fill(nMatchClusters);
483         
484         // compute track position at the end of the absorber
485         AliMUONTrackParam trackParamAtAbsEnd(*((AliMUONTrackParam*)trackMatched->GetTrackParamAtCluster()->First()));
486         AliMUONTrackExtrap::ExtrapToZ(&trackParamAtAbsEnd, AliMUONConstants::AbsZEnd());
487         xAbs = trackParamAtAbsEnd.GetNonBendingCoor();
488         yAbs = trackParamAtAbsEnd.GetBendingCoor();
489         dAbs = TMath::Sqrt(xAbs*xAbs + yAbs*yAbs);
490         aAbs = TMath::ATan(-dAbs/AliMUONConstants::AbsZEnd()) * TMath::RadToDeg();
491         pX2 = trackParamAtAbsEnd.Px();
492         pY2 = trackParamAtAbsEnd.Py();
493         pZ2 = trackParamAtAbsEnd.Pz();
494         pT2 = TMath::Sqrt(pX2*pX2 + pY2*pY2);
495         aMCS = TMath::ATan(-pT2/pZ2) * TMath::RadToDeg();
496         
497         trackParam = trackRef->GetTrackParamAtVertex();
498         x1 = trackParam->GetNonBendingCoor();
499         y1 = trackParam->GetBendingCoor();
500         z1 = trackParam->GetZ();
501         slopex1 = trackParam->GetNonBendingSlope();
502         slopey1 = trackParam->GetBendingSlope();
503         pX1 = trackParam->Px();
504         pY1 = trackParam->Py();
505         pZ1 = trackParam->Pz();
506         p1  = trackParam->P();
507         pT1 = TMath::Sqrt(pX1*pX1 + pY1*pY1);
508         aMC = TMath::ATan(-pT1/pZ1) * TMath::RadToDeg();
509         eta1 = TMath::Log(TMath::Tan(0.5*TMath::ATan(-pT1/pZ1)));
510         phi1 = TMath::Pi()+TMath::ATan2(-pY1, -pX1);
511         
512         trackParam = trackMatched->GetTrackParamAtVertex();
513         x2 = trackParam->GetNonBendingCoor();
514         y2 = trackParam->GetBendingCoor();
515         z2 = trackParam->GetZ();
516         slopex2 = trackParam->GetNonBendingSlope();
517         slopey2 = trackParam->GetBendingSlope();
518         pX2 = trackParam->Px();
519         pY2 = trackParam->Py();
520         pZ2 = trackParam->Pz();
521         p2  = trackParam->P();
522         pT2 = TMath::Sqrt(pX2*pX2 + pY2*pY2);
523         eta2 = TMath::Log(TMath::Tan(0.5*TMath::ATan(-pT2/pZ2)));
524         phi2 = TMath::Pi()+TMath::ATan2(-pY2, -pX2);
525         
526         dPhi = phi2-phi1;
527         if (dPhi < -TMath::Pi()) dPhi += 2.*TMath::Pi();
528         else if (dPhi > TMath::Pi()) dPhi -= 2.*TMath::Pi();
529         
530         AliMUONTrackParam trackParamAtDCA(*((AliMUONTrackParam*) trackMatched->GetTrackParamAtCluster()->First()));
531         pU = trackParamAtDCA.P();
532         AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(&trackParamAtDCA, z2);
533         xDCA = trackParamAtDCA.GetNonBendingCoor();
534         yDCA = trackParamAtDCA.GetBendingCoor();
535         dca = TMath::Sqrt(xDCA*xDCA + yDCA*yDCA);
536         
537         hResMomVertex->Fill(p2-p1);
538         hResSlopeXVertex->Fill(slopex2-slopex1);
539         hResSlopeYVertex->Fill(slopey2-slopey1);
540         hPDCA->Fill(0.5*(p2+pU)*dca);
541         hResEtaVertex->Fill(eta2-eta1);
542         hResPhiVertex->Fill(dPhi);
543         if (aMC >= aAbsLimits[0] && aMC <= aAbsLimits[1]) {
544           hResMomVertexVsMom->Fill(p1,p2-p1);
545           hResSlopeXVertexVsMom->Fill(p1,slopex2-slopex1);
546           hResSlopeYVertexVsMom->Fill(p1,slopey2-slopey1);
547           hResEtaVertexVsMom->Fill(p1,eta2-eta1);
548           hResPhiVertexVsMom->Fill(p1,dPhi);
549         }
550         hResMomVertexVsAngleVsMom->Fill(p1,aAbs,p2-p1);
551         if (aAbs > 2. && aAbs < 3.) {
552           hResMomVertexVsMom_2_3_Deg->Fill(p1,p2-p1);
553           hPDCAVsMom_2_3_Deg->Fill(p1,0.5*(p2+pU)*dca);
554           hPMCSAngVsMom_2_3_Deg->Fill(p1,0.5*(p2+pU)*(aMCS-aMC)*TMath::DegToRad());
555         }
556         else if (aAbs >= 3. && aAbs < 10.) {
557           hResMomVertexVsMom_3_10_Deg->Fill(p1,p2-p1);
558           hPDCAVsMom_3_10_Deg->Fill(p1,0.5*(p2+pU)*dca);
559           hPMCSAngVsMom_3_10_Deg->Fill(p1,0.5*(p2+pU)*(aMCS-aMC)*TMath::DegToRad());
560           aMCSMoy += 0.5*(p2+pU)*(aMCS-aMC)*TMath::DegToRad();
561           aMCS2Moy += (0.5*(p2+pU)*(aMCS-aMC)*TMath::DegToRad()) * (0.5*(p2+pU)*(aMCS-aMC)*TMath::DegToRad());
562           dMCSMoy += 0.5*(p2+pU)*(dAbs-pT1/pZ1*AliMUONConstants::AbsZEnd());
563           dMCS2Moy += (0.5*(p2+pU)*(dAbs-pT1/pZ1*AliMUONConstants::AbsZEnd())) * (0.5*(p2+pU)*(dAbs-pT1/pZ1*AliMUONConstants::AbsZEnd()));
564           adMCSMoy += (0.5*(p2+pU)*(aMCS-aMC)*TMath::DegToRad()) * (0.5*(p2+pU)*(dAbs-pT1/pZ1*AliMUONConstants::AbsZEnd()));
565           nMCS++;
566         }
567         if (aMC < 2.) {
568           hResMomVertexVsMom_0_2_DegMC->Fill(p1,p2-p1);
569           hResMomVertexVsPosAbsEnd_0_2_DegMC->Fill(dAbs,p2-p1);
570           hResSlopeXVertexVsPosAbsEnd_0_2_DegMC->Fill(dAbs,slopex2-slopex1);
571           hResSlopeYVertexVsPosAbsEnd_0_2_DegMC->Fill(dAbs,slopey2-slopey1);
572           hPDCAVsPosAbsEnd_0_2_DegMC->Fill(dAbs,0.5*(p2+pU)*dca);
573           hResEtaVertexVsPosAbsEnd_0_2_DegMC->Fill(dAbs,eta2-eta1);
574           hResPhiVertexVsPosAbsEnd_0_2_DegMC->Fill(dAbs,dPhi);
575         }
576         else if (aMC >= 2. && aMC < 3) {
577           hResMomVertexVsPosAbsEnd_2_3_DegMC->Fill(dAbs,p2-p1);
578           hResSlopeXVertexVsPosAbsEnd_2_3_DegMC->Fill(dAbs,slopex2-slopex1);
579           hResSlopeYVertexVsPosAbsEnd_2_3_DegMC->Fill(dAbs,slopey2-slopey1);
580           hPDCAVsPosAbsEnd_2_3_DegMC->Fill(dAbs,0.5*(p2+pU)*dca);
581           hResEtaVertexVsPosAbsEnd_2_3_DegMC->Fill(dAbs,eta2-eta1);
582           hResPhiVertexVsPosAbsEnd_2_3_DegMC->Fill(dAbs,dPhi);
583         }
584         else if (aMC >= 3. && aMC < 10.) {
585           hResMomVertexVsPosAbsEnd_3_10_DegMC->Fill(dAbs,p2-p1);
586           hResSlopeXVertexVsPosAbsEnd_3_10_DegMC->Fill(dAbs,slopex2-slopex1);
587           hResSlopeYVertexVsPosAbsEnd_3_10_DegMC->Fill(dAbs,slopey2-slopey1);
588           hPDCAVsPosAbsEnd_3_10_DegMC->Fill(dAbs,0.5*(p2+pU)*dca);
589           hResEtaVertexVsPosAbsEnd_3_10_DegMC->Fill(dAbs,eta2-eta1);
590           hResPhiVertexVsPosAbsEnd_3_10_DegMC->Fill(dAbs,dPhi);
591         }
592         hResMomVertexVsAngle->Fill(aAbs,p2-p1);
593         hResSlopeXVertexVsAngle->Fill(aAbs,slopex2-slopex1);
594         hResSlopeYVertexVsAngle->Fill(aAbs,slopey2-slopey1);
595         hPDCAVsAngle->Fill(aAbs,0.5*(p2+pU)*dca);
596         hResEtaVertexVsAngle->Fill(aAbs,eta2-eta1);
597         hResPhiVertexVsAngle->Fill(aAbs,dPhi);
598         hResMomVertexVsMCAngle->Fill(aMC,p2-p1);
599         hResSlopeXVertexVsMCAngle->Fill(aMC,slopex2-slopex1);
600         hResSlopeYVertexVsMCAngle->Fill(aMC,slopey2-slopey1);
601         hPDCAVsMCAngle->Fill(aMC,0.5*(p2+pU)*dca);
602         hResEtaVertexVsMCAngle->Fill(aMC,eta2-eta1);
603         hResPhiVertexVsMCAngle->Fill(aMC,dPhi);
604         
605         trackParam = (AliMUONTrackParam*) trackRef->GetTrackParamAtCluster()->First();
606         x1 = trackParam->GetNonBendingCoor();
607         y1 = trackParam->GetBendingCoor();
608         z1 = trackParam->GetZ();
609         slopex1 = trackParam->GetNonBendingSlope();
610         slopey1 = trackParam->GetBendingSlope();
611         pX1 = trackParam->Px();
612         pY1 = trackParam->Py();
613         pZ1 = trackParam->Pz();
614         p1  = trackParam->P();
615         pT1 = TMath::Sqrt(pX1*pX1 + pY1*pY1);
616         
617         trackParam = (AliMUONTrackParam*) trackMatched->GetTrackParamAtCluster()->First();
618         x2 = trackParam->GetNonBendingCoor();
619         y2 = trackParam->GetBendingCoor();
620         z2 = trackParam->GetZ();
621         slopex2 = trackParam->GetNonBendingSlope();
622         slopey2 = trackParam->GetBendingSlope();
623         pX2 = trackParam->Px();
624         pY2 = trackParam->Py();
625         pZ2 = trackParam->Pz();
626         p2  = trackParam->P();
627         pT2 = TMath::Sqrt(pX2*pX2 + pY2*pY2);
628         
629         hResMomFirstCluster->Fill(p2-p1);
630         hResMomFirstClusterVsMom->Fill(p1,p2-p1);
631         
632         hResSlopeXFirstCluster->Fill(slopex2-slopex1);
633         hResSlopeYFirstCluster->Fill(slopey2-slopey1);
634         hResSlopeXFirstClusterVsMom->Fill(p1,slopex2-slopex1);
635         hResSlopeYFirstClusterVsMom->Fill(p1,slopey2-slopey1);
636         
637         // Fill residuals
638         // Loop over clusters of first track
639         AliMUONTrackParam* trackParamAtCluster1 = (AliMUONTrackParam*) trackMatched->GetTrackParamAtCluster()->First();
640         while (trackParamAtCluster1) {
641           AliMUONVCluster* cluster1 = trackParamAtCluster1->GetClusterPtr();
642           AliMUONTrackParam* trackParamAtCluster2 = (AliMUONTrackParam*) trackRef->GetTrackParamAtCluster()->First();
643           while (trackParamAtCluster2) {
644             AliMUONVCluster* cluster2 = trackParamAtCluster2->GetClusterPtr();
645             if (cluster1->GetDetElemId() == cluster2->GetDetElemId()) {
646               hResidualXInCh[cluster1->GetChamberId()]->Fill(cluster1->GetX() - cluster2->GetX());
647               hResidualYInCh[cluster1->GetChamberId()]->Fill(cluster1->GetY() - cluster2->GetY());
648               break;
649             }
650             trackParamAtCluster2 = (AliMUONTrackParam*) trackRef->GetTrackParamAtCluster()->After(trackParamAtCluster2);
651           }
652           trackParamAtCluster1 = (AliMUONTrackParam*) trackMatched->GetTrackParamAtCluster()->After(trackParamAtCluster1);
653         }
654         
655       }
656       
657     } // end loop track ref.
658
659   } // end loop on event  
660   cout<<"\rEvent processing... "<<nevents<<" done"<<endl;
661   
662   // ###################################### compute stuff ###################################### //
663   cout<<"\nWhen not specified, resolution at vertex is computed for ";
664   if (absorberRegion == 1) cout<<"tracks in the absorber region [2,3] deg."<<endl;
665   else if (absorberRegion == 2) cout<<"tracks in the absorber region [3,10] deg."<<endl;
666   else cout<<"all tracks"<<endl;
667   
668   // compute momentum resolution at vertex versus p
669   TF1 *f2 = new TF1("f2",langaufun,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1],4);
670   Int_t rebinFactorX = TMath::Max(hResMomVertexVsMom->GetNbinsX()/pNBins, 1);
671   for (Int_t i = rebinFactorX; i <= hResMomVertexVsMom->GetNbinsX(); i+=rebinFactorX) {
672     cout<<"\rFitting momentum residuals at vertex... "<<i/rebinFactorX<<"/"<<pNBins<<flush;
673     TH1D *tmp = hResMomVertexVsMom->ProjectionY("tmp",i-rebinFactorX+1,i,"e");
674     f2->SetParameters(0.2,0.,(Double_t)tmp->GetEntries(),1.);
675     tmp->Fit("f2","WWNQ");
676     Double_t fwhm = f2->GetParameter(0);
677     Double_t sigma = f2->GetParameter(3);
678     Double_t sigmaP = TMath::Sqrt(sigma*sigma + fwhm*fwhm/(8.*log(2.)));
679     Int_t rebin = TMath::Max(Int_t(0.5*sigmaP/tmp->GetBinWidth(1)),1);
680     while (deltaPAtVtxNBins%rebin!=0) rebin--;
681     tmp->Rebin(rebin);
682     tmp->Fit("f2","NQ");
683     fwhm = f2->GetParameter(0);
684     sigma = f2->GetParameter(3);
685     sigmaP = TMath::Sqrt(sigma*sigma + fwhm*fwhm/(8.*log(2.)));
686     Double_t fwhmErr = f2->GetParError(0);
687     Double_t sigmaErr = f2->GetParError(3);
688     Double_t sigmaPErr = TMath::Sqrt(sigma*sigma*sigmaErr*sigmaErr + fwhm*fwhm*fwhmErr*fwhmErr/(64.*log(2.)*log(2.))) / sigmaP;
689     hResMomVertexVsMom->GetXaxis()->SetRange(i-rebinFactorX+1,i);
690     Double_t p = (tmp->GetEntries() > 0) ? hResMomVertexVsMom->GetMean() : 0.5 * (hResMomVertexVsMom->GetBinLowEdge(i-rebinFactorX+1) + hResMomVertexVsMom->GetBinLowEdge(i+1));
691     hResMomVertexVsMom->GetXaxis()->SetRange();
692     Double_t pErr[2] = {p-hResMomVertexVsMom->GetBinLowEdge(i-rebinFactorX+1), hResMomVertexVsMom->GetBinLowEdge(i+1)-p};
693     gMeanResMomVertexVsMom->SetPoint(i/rebinFactorX-1, p, tmp->GetMean());
694     gMeanResMomVertexVsMom->SetPointError(i/rebinFactorX-1, pErr[0], pErr[1], tmp->GetMeanError(), tmp->GetMeanError());
695     gMostProbResMomVertexVsMom->SetPoint(i/rebinFactorX-1, p, -f2->GetParameter(1));
696     gMostProbResMomVertexVsMom->SetPointError(i/rebinFactorX-1, pErr[0], pErr[1], f2->GetParError(1), f2->GetParError(1));
697     gSigmaResMomVertexVsMom->SetPoint(i/rebinFactorX-1, p, 100.*sigmaP/p);
698     gSigmaResMomVertexVsMom->SetPointError(i/rebinFactorX-1, pErr[0], pErr[1], 100.*sigmaPErr/p, 100.*sigmaPErr/p);
699     delete tmp;
700   }
701   cout<<"\rFitting momentum residuals at vertex... "<<pNBins<<"/"<<pNBins<<endl;
702   
703   // compute momentum relative resolution at first cluster versus p
704   FitGausResVsMom(hResMomFirstClusterVsMom, pNBins, 0., 1., "momentum residuals at first cluster", gMeanResMomFirstClusterVsMom, gSigmaResMomFirstClusterVsMom);
705   rebinFactorX = TMath::Max(hResMomFirstClusterVsMom->GetNbinsX()/pNBins, 1);
706   for (Int_t i = rebinFactorX; i <= hResMomFirstClusterVsMom->GetNbinsX(); i+=rebinFactorX) {
707     Double_t x,y;
708     gSigmaResMomFirstClusterVsMom->GetPoint(i/rebinFactorX-1, x, y);
709     gSigmaResMomFirstClusterVsMom->SetPoint(i/rebinFactorX-1, x, 100.*y/x);
710     gSigmaResMomFirstClusterVsMom->SetPointEYlow(i/rebinFactorX-1, 100.*gSigmaResMomFirstClusterVsMom->GetErrorYlow(i/rebinFactorX-1)/x);
711     gSigmaResMomFirstClusterVsMom->SetPointEYhigh(i/rebinFactorX-1, 100.*gSigmaResMomFirstClusterVsMom->GetErrorYhigh(i/rebinFactorX-1)/x);
712   }
713   
714   // compute slopeX resolution at vertex versus p
715   FitGausResVsMom(hResSlopeXVertexVsMom, pNBins, 0., 2.e-3, "slopeX residuals at vertex", gMeanResSlopeXVertexVsMom, gSigmaResSlopeXVertexVsMom);
716   
717   // compute slopeY resolution at vertex versus p
718   FitGausResVsMom(hResSlopeYVertexVsMom, pNBins, 0., 2.e-3, "slopeY residuals at vertex", gMeanResSlopeYVertexVsMom, gSigmaResSlopeYVertexVsMom);
719   
720   // compute slopeX resolution at first cluster versus p
721   FitGausResVsMom(hResSlopeXFirstClusterVsMom, pNBins, 0., 3.e-4, "slopeX residuals at first cluster", gMeanResSlopeXFirstClusterVsMom, gSigmaResSlopeXFirstClusterVsMom);
722   
723   // compute slopeY resolution at first cluster versus p
724   FitGausResVsMom(hResSlopeYFirstClusterVsMom, pNBins, 0., 2.e-4, "slopeY residuals at first cluster", gMeanResSlopeYFirstClusterVsMom, gSigmaResSlopeYFirstClusterVsMom);
725   
726   // compute p*DCA resolution in the region [2,3] deg at absorber end
727   FitPDCAVsMom(hPDCAVsMom_2_3_Deg, pNBins, "p*DCA (tracks in [2,3] deg.)", gMeanPDCAVsMom_2_3_Deg, gSigmaPDCAVsMom_2_3_Deg);
728   
729   // compute p*DCA resolution in the region [3,10] deg at absorber end
730   FitPDCAVsMom(hPDCAVsMom_3_10_Deg, pNBins, "p*DCA (tracks in [3,10] deg.)", gMeanPDCAVsMom_3_10_Deg, gSigmaPDCAVsMom_3_10_Deg);
731   
732   // compute MCS angular dispersion in the region [2,3] deg at absorber end
733   FitGausResVsMom(hPMCSAngVsMom_2_3_Deg, pNBins, 0., 2.e-3, "p*MCSAngle (tracks in [2,3] deg.)", gMeanPMCSAngVsMom_2_3_Deg, gSigmaPMCSAngVsMom_2_3_Deg);
734   
735   // compute MCS angular dispersion in the region [3,10] deg at absorber end
736   FitGausResVsMom(hPMCSAngVsMom_3_10_Deg, pNBins, 0., 2.e-3, "p*MCSAngle (tracks in [3,10] deg.)", gMeanPMCSAngVsMom_3_10_Deg, gSigmaPMCSAngVsMom_3_10_Deg);
737   
738   // compute eta resolution at vertex versus p
739   FitGausResVsMom(hResEtaVertexVsMom, pNBins, 0., 0.1, "eta residuals at vertex", gMeanResEtaVertexVsMom, gSigmaResEtaVertexVsMom);
740   
741   // compute phi resolution at vertex versus p
742   FitGausResVsMom(hResPhiVertexVsMom, pNBins, 0., 0.01, "phi residuals at vertex", gMeanResPhiVertexVsMom, gSigmaResPhiVertexVsMom);
743   
744   // compute cluster-track residual mean and dispersion
745   for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
746     hResidualXInCh[i]->GetXaxis()->SetRangeUser(-3.*hResidualXInCh[i]->GetRMS(), 3.*hResidualXInCh[i]->GetRMS());
747     gResidualXPerChMean->SetPoint(i, i+1, hResidualXInCh[i]->GetMean());
748     gResidualXPerChMean->SetPointError(i, 0., hResidualXInCh[i]->GetMeanError());
749     gResidualXPerChSigma->SetPoint(i, i+1, hResidualXInCh[i]->GetRMS());
750     gResidualXPerChSigma->SetPointError(i, 0., hResidualXInCh[i]->GetRMSError());
751     hResidualXInCh[i]->GetXaxis()->SetRange(0,0);
752     hResidualYInCh[i]->GetXaxis()->SetRangeUser(-3.*hResidualYInCh[i]->GetRMS(), 3.*hResidualYInCh[i]->GetRMS());
753     gResidualYPerChMean->SetPoint(i, i+1, hResidualYInCh[i]->GetMean());
754     gResidualYPerChMean->SetPointError(i, 0., hResidualYInCh[i]->GetMeanError());
755     gResidualYPerChSigma->SetPoint(i, i+1, hResidualYInCh[i]->GetRMS());
756     gResidualYPerChSigma->SetPointError(i, 0., hResidualYInCh[i]->GetRMSError());
757     hResidualYInCh[i]->GetXaxis()->SetRange(0,0);
758   }
759   
760   // ###################################### display histograms ###################################### //
761   // diplay momentum residuals
762   TCanvas* cResMom = DrawVsAng("cResMom", "momentum residual at vertex in 3 angular regions", hResMomVertex, hResMomVertexVsAngle);
763   TCanvas* cResMomMC = DrawVsAng("cResMomMC", "momentum residual at vertex in 3 MC angular regions", hResMomVertex, hResMomVertexVsMCAngle);
764   TCanvas* cResMomVsPos = DrawVsPos("cResMomVsPos", "momentum residual at vertex versus position at absorber end in 3 MC angular regions",
765                                     hResMomVertexVsPosAbsEnd_0_2_DegMC, hResMomVertexVsPosAbsEnd_2_3_DegMC, hResMomVertexVsPosAbsEnd_3_10_DegMC);
766   TCanvas* cResMom_2_3_Deg = DrawResMomVsMom("cResMom_2_3_Deg", "momentum residual for tracks between 2 and 3 degrees",
767                                              hResMomVertexVsMom_2_3_Deg, 10, f2, "momentum residuals at vertex (tracks in [2,3] deg.)");
768   TCanvas* cResMom_3_10_Deg = DrawResMomVsMom("cResMom_3_10_Deg", "momentum residual for tracks between 3 and 10 degrees",
769                                               hResMomVertexVsMom_3_10_Deg, 10, f2, "momentum residuals at vertex (tracks in [3,10] deg.)");
770   TCanvas* cResMom_0_2_DegMC = DrawResMomVsMom("cResMom_0_2_DegMC", "momentum residuals for tracks with MC angle < 2 degrees", hResMomVertexVsMom_0_2_DegMC, 5);
771   
772   // diplay slopeX residuals
773   TCanvas* cResSlopeX = DrawVsAng("cResSlopeX", "slope_{X} residual at vertex in 3 angular regions", hResSlopeXVertex, hResSlopeXVertexVsAngle);
774   TCanvas* cResSlopeXMC = DrawVsAng("cResSlopeXMC", "slope_{X} residual at vertex in 3 MC angular regions", hResSlopeXVertex, hResSlopeXVertexVsMCAngle);
775   TCanvas* cResSlopeXVsPos = DrawVsPos("cResSlopeXVsPos", "slope_{X} residual at vertex versus position at absorber end in 3 MC angular regions",
776                                        hResSlopeXVertexVsPosAbsEnd_0_2_DegMC, hResSlopeXVertexVsPosAbsEnd_2_3_DegMC, hResSlopeXVertexVsPosAbsEnd_3_10_DegMC);
777   
778   // diplay slopeY residuals
779   TCanvas* cResSlopeY = DrawVsAng("cResSlopeY", "slope_{Y} residual at vertex in 3 angular regions", hResSlopeYVertex, hResSlopeYVertexVsAngle);
780   TCanvas* cResSlopeYMC = DrawVsAng("cResSlopeYMC", "slope_{Y} residual at vertex in 3 MC angular regions", hResSlopeYVertex, hResSlopeYVertexVsMCAngle);
781   TCanvas* cResSlopeYVsPos = DrawVsPos("cResSlopeYVsPos", "slope_{Y} residual at vertex versus position at absorber end in 3 MC angular regions",
782                                        hResSlopeYVertexVsPosAbsEnd_0_2_DegMC, hResSlopeYVertexVsPosAbsEnd_2_3_DegMC, hResSlopeYVertexVsPosAbsEnd_3_10_DegMC);
783   
784   // diplay P*DCA
785   TCanvas* cPDCA = DrawVsAng("cPDCA", "p #times DCA in 3 angular regions", hPDCA, hPDCAVsAngle);
786   TCanvas* cPDCAMC = DrawVsAng("cPDCAMC", "p #times DCA in 3 MC angular regions", hPDCA, hPDCAVsMCAngle);
787   TCanvas* cPDCAVsPos = DrawVsPos("cPDCAVsPos", "p #times DCA versus position at absorber end in 3 MC angular regions",
788                                   hPDCAVsPosAbsEnd_0_2_DegMC, hPDCAVsPosAbsEnd_2_3_DegMC, hPDCAVsPosAbsEnd_3_10_DegMC);
789   
790   // diplay eta residuals
791   TCanvas* cResEta = DrawVsAng("cResEta", "eta residual at vertex in 3 angular regions", hResEtaVertex, hResEtaVertexVsAngle);
792   TCanvas* cResEtaMC = DrawVsAng("cResEtaMC", "eta residual at vertex in 3 MC angular regions", hResEtaVertex, hResEtaVertexVsMCAngle);
793   TCanvas* cResEtaVsPos = DrawVsPos("cResEtaVsPos", "eta residual at vertex versus position at absorber end in 3 MC angular regions",
794                                     hResEtaVertexVsPosAbsEnd_0_2_DegMC, hResEtaVertexVsPosAbsEnd_2_3_DegMC, hResEtaVertexVsPosAbsEnd_3_10_DegMC);
795   
796   // diplay phi residuals
797   TCanvas* cResPhi = DrawVsAng("cResPhi", "phi residual at vertex in 3 angular regions", hResPhiVertex, hResPhiVertexVsAngle);
798   TCanvas* cResPhiMC = DrawVsAng("cResPhiMC", "phi residual at vertex in 3 MC angular regions", hResPhiVertex, hResPhiVertexVsMCAngle);
799   TCanvas* cResPhiVsPos = DrawVsPos("cResPhiVsPos", "phi residual at vertex versus position at absorber end in 3 MC angular regions",
800                                     hResPhiVertexVsPosAbsEnd_0_2_DegMC, hResPhiVertexVsPosAbsEnd_2_3_DegMC, hResPhiVertexVsPosAbsEnd_3_10_DegMC);
801   
802   // ###################################### save histogram ###################################### //
803   histoFile->Write();
804   
805   histoFile->cd("momentumAtVertex");
806   gMeanResMomVertexVsMom->Write();
807   gMostProbResMomVertexVsMom->Write();
808   gSigmaResMomVertexVsMom->Write();
809   cResMom->Write();
810   cResMomMC->Write();
811   cResMomVsPos->Write();
812   cResMom_2_3_Deg->Write();
813   cResMom_3_10_Deg->Write();
814   cResMom_0_2_DegMC->Write();
815   
816   histoFile->cd("slopesAtVertex");
817   gMeanResSlopeXVertexVsMom->Write();
818   gMeanResSlopeYVertexVsMom->Write();
819   gSigmaResSlopeXVertexVsMom->Write();
820   gSigmaResSlopeYVertexVsMom->Write();
821   cResSlopeX->Write();
822   cResSlopeY->Write();
823   cResSlopeXMC->Write();
824   cResSlopeYMC->Write();
825   cResSlopeXVsPos->Write();
826   cResSlopeYVsPos->Write();
827   
828   histoFile->cd("DCA");
829   gMeanPDCAVsMom_2_3_Deg->Write();
830   gSigmaPDCAVsMom_2_3_Deg->Write();
831   gMeanPDCAVsMom_3_10_Deg->Write();
832   gSigmaPDCAVsMom_3_10_Deg->Write();
833   gMeanPMCSAngVsMom_2_3_Deg->Write();
834   gSigmaPMCSAngVsMom_2_3_Deg->Write();
835   gMeanPMCSAngVsMom_3_10_Deg->Write();
836   gSigmaPMCSAngVsMom_3_10_Deg->Write();
837   cPDCA->Write();
838   cPDCAMC->Write();
839   cPDCAVsPos->Write();
840   
841   histoFile->cd("etaAtVertex");
842   gMeanResEtaVertexVsMom->Write();
843   gSigmaResEtaVertexVsMom->Write();
844   cResEta->Write();
845   cResEtaMC->Write();
846   cResEtaVsPos->Write();
847   
848   histoFile->cd("phiAtVertex");
849   gMeanResPhiVertexVsMom->Write();
850   gSigmaResPhiVertexVsMom->Write();
851   cResPhi->Write();
852   cResPhiMC->Write();
853   cResPhiVsPos->Write();
854   
855   histoFile->cd("momentumAtFirstCluster");
856   gMeanResMomFirstClusterVsMom->Write();
857   gSigmaResMomFirstClusterVsMom->Write();
858   
859   histoFile->cd("slopesAtFirstCluster");
860   gMeanResSlopeXFirstClusterVsMom->Write();
861   gMeanResSlopeYFirstClusterVsMom->Write();
862   gSigmaResSlopeXFirstClusterVsMom->Write();
863   gSigmaResSlopeYFirstClusterVsMom->Write();
864   
865   histoFile->cd("clusters");
866   gResidualXPerChMean->Write();
867   gResidualXPerChSigma->Write();
868   gResidualYPerChMean->Write();
869   gResidualYPerChSigma->Write();
870   
871   histoFile->Close();
872   
873   // ###################################### clean memory ###################################### //
874   delete cResMom;
875   delete cResMomMC;
876   delete cResMomVsPos;
877   delete cResMom_2_3_Deg;
878   delete cResMom_3_10_Deg;
879   delete cResMom_0_2_DegMC;
880   delete cResSlopeX;
881   delete cResSlopeY;
882   delete cResSlopeXMC;
883   delete cResSlopeYMC;
884   delete cResSlopeXVsPos;
885   delete cResSlopeYVsPos;
886   delete cPDCA;
887   delete cPDCAMC;
888   delete cPDCAVsPos;
889   delete cResEta;
890   delete cResEtaMC;
891   delete cResEtaVsPos;
892   delete cResPhi;
893   delete cResPhiMC;
894   delete cResPhiVsPos;
895   
896   // ###################################### print statistics ###################################### //
897   printf("\n");
898   printf("nb of reconstructible tracks: %d \n", nReconstructibleTracks);
899   printf("nb of reconstructed tracks: %d \n", nReconstructedTracks);
900   printf("nb of reconstructible tracks which are reconstructed: %d \n", nReconstructibleTracksCheck);
901   
902   aMCSMoy /= (Double_t) nMCS;
903   aMCS2Moy /= (Double_t) nMCS;
904   dMCSMoy /= (Double_t) nMCS;
905   dMCS2Moy /= (Double_t) nMCS;
906   adMCSMoy /= (Double_t) nMCS;
907   Double_t sigma2_ThetaMCS = aMCS2Moy - aMCSMoy*aMCSMoy;
908   Double_t sigma2_PosMCS = dMCS2Moy - dMCSMoy*dMCSMoy;
909   Double_t cov_ThetaPosMCS = - (adMCSMoy - aMCSMoy*dMCSMoy);
910   printf("\nmultiple scattering of tracks between 3 and 10 deg. at absorber end:\n");
911   printf(" sigma_ThetaMCS = %f\n", TMath::Sqrt(sigma2_ThetaMCS));
912   printf(" sigma_PosMCS = %f\n", TMath::Sqrt(sigma2_PosMCS));
913   printf(" cov_ThetaPosMCS = %f\n", cov_ThetaPosMCS);
914   printf(" --> sigma_DCA = %f\n", TMath::Sqrt(AliMUONConstants::AbsZEnd()*AliMUONConstants::AbsZEnd()*sigma2_ThetaMCS
915                                               - 2.*AliMUONConstants::AbsZEnd()*cov_ThetaPosMCS + sigma2_PosMCS));
916   printf("\n");
917 }
918
919 //------------------------------------------------------------------------------------
920 Double_t langaufun(Double_t *x, Double_t *par) {
921   
922   //Fit parameters:
923   //par[0]=Width (scale) parameter of Landau density
924   //par[1]=Most Probable (MP, location) parameter of Landau density
925   //par[2]=Total area (integral -inf to inf, normalization constant)
926   //par[3]=Width (sigma) of convoluted Gaussian function
927   //
928   //In the Landau distribution (represented by the CERNLIB approximation), 
929   //the maximum is located at x=-0.22278298 with the location parameter=0.
930   //This shift is corrected within this function, so that the actual
931   //maximum is identical to the MP parameter.
932   
933   // Numeric constants
934   Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
935   Double_t mpshift  = -0.22278298;       // Landau maximum location
936   
937   // Control constants
938   Double_t np = 100.0; // number of convolution steps
939   Double_t sc = 5.0;   // convolution extends to +-sc Gaussian sigmas
940   
941   // Variables
942   Double_t xx;
943   Double_t mpc;
944   Double_t fland;
945   Double_t sum = 0.0;
946   Double_t xlow,xupp;
947   Double_t step;
948   Double_t i;
949   
950   
951   // MP shift correction
952   mpc = par[1] - mpshift * par[0]; 
953   
954   // Range of convolution integral
955   xlow = x[0] - sc * par[3];
956   xupp = x[0] + sc * par[3];
957   
958   step = (xupp-xlow) / np;
959   
960   // Convolution integral of Landau and Gaussian by sum
961   for(i=1.0; i<=np/2; i++) {
962     xx = xlow + (i-.5) * step;
963     //change x -> -x because the tail of the Landau is at the left here...
964     fland = TMath::Landau(-xx,mpc,par[0]) / par[0];
965     sum += fland * TMath::Gaus(x[0],xx,par[3]);
966     
967     //change x -> -x because the tail of the Landau is at the left here...
968     xx = xupp - (i-.5) * step;
969     fland = TMath::Landau(-xx,mpc,par[0]) / par[0];
970     sum += fland * TMath::Gaus(x[0],xx,par[3]);
971   }
972   
973   return (par[2] * step * sum * invsq2pi / par[3]);
974 }
975
976 //------------------------------------------------------------------------------------
977 void FitGausResVsMom(TH2* h, Int_t nBins, const Double_t mean0, const Double_t sigma0,
978                      const char* fitting, TGraphAsymmErrors* gMean, TGraphAsymmErrors* gSigma)
979 {
980   /// generic function to fit residuals versus momentum with a gaussian
981   static TF1* fGaus = 0x0;
982   if (!fGaus) fGaus = new TF1("fGaus","gaus");
983   
984   Int_t rebinFactorX = TMath::Max(h->GetNbinsX()/nBins, 1);
985   for (Int_t i = rebinFactorX; i <= h->GetNbinsX(); i+=rebinFactorX) {
986     cout<<Form("\rFitting %s... %d/%d",fitting,i/rebinFactorX,nBins)<<flush;
987     TH1D *tmp = h->ProjectionY("tmp",i-rebinFactorX+1,i,"e");
988     fGaus->SetParameters(tmp->GetEntries(), mean0, sigma0);
989     tmp->Fit("fGaus","WWNQ");
990     Int_t rebin = TMath::Max(Int_t(0.5*fGaus->GetParameter(2)/tmp->GetBinWidth(1)),1);
991     while (tmp->GetNbinsX()%rebin!=0) rebin--;
992     tmp->Rebin(rebin);
993     tmp->Fit("fGaus","NQ");
994     h->GetXaxis()->SetRange(i-rebinFactorX+1,i);
995     Double_t p = (tmp->GetEntries() > 0) ? h->GetMean() : 0.5 * (h->GetBinLowEdge(i-rebinFactorX+1) + h->GetBinLowEdge(i+1));
996     h->GetXaxis()->SetRange();
997     Double_t pErr[2] = {p-h->GetBinLowEdge(i-rebinFactorX+1), h->GetBinLowEdge(i+1)-p};
998     gMean->SetPoint(i/rebinFactorX-1, p, fGaus->GetParameter(1));
999     gMean->SetPointError(i/rebinFactorX-1, pErr[0], pErr[1], fGaus->GetParError(1), fGaus->GetParError(1));
1000     gSigma->SetPoint(i/rebinFactorX-1, p, fGaus->GetParameter(2));
1001     gSigma->SetPointError(i/rebinFactorX-1, pErr[0], pErr[1], fGaus->GetParError(2), fGaus->GetParError(2));
1002     delete tmp;
1003   }
1004   cout<<Form("\rFitting %s... %d/%d",fitting,nBins,nBins)<<endl;
1005 }
1006
1007 //------------------------------------------------------------------------------------
1008 void FitPDCAVsMom(TH2* h, Int_t nBins, const char* fitting, TGraphAsymmErrors* gMean, TGraphAsymmErrors* gSigma)
1009 {
1010   /// generic function to fit p*DCA distributions
1011   static TF1* fPGaus = 0x0;
1012   if (!fPGaus) fPGaus = new TF1("fPGaus","x*gaus");
1013   
1014   Int_t rebinFactorX = TMath::Max(h->GetNbinsX()/nBins, 1);
1015   for (Int_t i = rebinFactorX; i <= h->GetNbinsX(); i+=rebinFactorX) {
1016     cout<<Form("\rFitting %s... %d/%d",fitting,i/rebinFactorX,nBins)<<flush;
1017     TH1D *tmp = h->ProjectionY("tmp",i-rebinFactorX+1,i,"e");
1018     fPGaus->SetParameters(1.,0.,80.);
1019     Int_t rebin = 25.*(tmp->GetNbinsX()/(tmp->GetBinLowEdge(tmp->GetNbinsX()+1)-tmp->GetBinLowEdge(1)));
1020     while (tmp->GetNbinsX()%rebin!=0) rebin--;
1021     tmp->Rebin(rebin);
1022     tmp->Fit("fPGaus","NQ");
1023     h->GetXaxis()->SetRange(i-rebinFactorX+1,i);
1024     Double_t p = (tmp->GetEntries() > 0) ? h->GetMean() : 0.5 * (h->GetBinLowEdge(i-rebinFactorX+1) + h->GetBinLowEdge(i+1));
1025     h->GetXaxis()->SetRange();
1026     Double_t pErr[2] = {p-h->GetBinLowEdge(i-rebinFactorX+1), h->GetBinLowEdge(i+1)-p};
1027     gMean->SetPoint(i/rebinFactorX-1, p, fPGaus->GetParameter(1));
1028     gMean->SetPointError(i/rebinFactorX-1, pErr[0], pErr[1], fPGaus->GetParError(1), fPGaus->GetParError(1));
1029     gSigma->SetPoint(i/rebinFactorX-1, p, fPGaus->GetParameter(2));
1030     gSigma->SetPointError(i/rebinFactorX-1, pErr[0], pErr[1], fPGaus->GetParError(2), fPGaus->GetParError(2));
1031     delete tmp;
1032   }
1033   cout<<Form("\rFitting %s... %d/%d",fitting,nBins,nBins)<<endl;
1034 }
1035
1036 //------------------------------------------------------------------------------------
1037 TCanvas* DrawVsAng(const char* name, const char* title, TH1* h1, TH2* h2)
1038 {
1039   /// generic function to draw histograms versus absorber angular region
1040   TCanvas* c = new TCanvas(name, title);
1041   c->cd();
1042   h1->Draw();
1043   TH1D *proj1 = h2->ProjectionY(Form("%s_proj_0_2",h2->GetName()),1,2);
1044   proj1->Draw("sames");
1045   proj1->SetLineColor(2);
1046   TH1D *proj2 = h2->ProjectionY(Form("%s_proj_2_3",h2->GetName()),3,3);
1047   proj2->Draw("sames");
1048   proj2->SetLineColor(4);
1049   TH1D *proj3 = h2->ProjectionY(Form("%s__proj_3_10",h2->GetName()),4,10);
1050   proj3->Draw("sames");
1051   proj3->SetLineColor(3);
1052   return c;
1053 }
1054
1055 //------------------------------------------------------------------------------------
1056 TCanvas* DrawVsPos(const char* name, const char* title, TH2* h1, TH2* h2, TH2* h3)
1057 {
1058   /// generic function to draw histograms versus position at absorber end
1059   TCanvas* c = new TCanvas(name, title);
1060   c->cd();
1061   h1->Draw();
1062   h1->SetMarkerColor(2);
1063   h2->Draw("sames");
1064   h2->SetMarkerColor(4);
1065   h3->Draw("sames");
1066   h3->SetMarkerColor(3);
1067   return c;
1068 }
1069
1070 //------------------------------------------------------------------------------------
1071 TCanvas* DrawResMomVsMom(const char* name, const char* title, TH2* h, Int_t nBins, TF1* f2, const char* fitting)
1072 {
1073   /// generic function to draw and eventually fit momentum residuals versus momentum
1074   TLegend* l = new TLegend(0.15,0.25,0.3,0.85);
1075   TCanvas* c = new TCanvas(name, title);
1076   c->cd();
1077   TH1D* proj = 0x0;
1078   h->Sumw2();
1079   Int_t rebinFactorX = TMath::Max(h->GetNbinsX()/nBins, 1);
1080   for (Int_t i = rebinFactorX; i <= h->GetNbinsX(); i+=rebinFactorX) {
1081     if (f2) cout<<Form("\rFitting %s... %d/%d",fitting,i/rebinFactorX,nBins)<<flush;
1082     proj = h->ProjectionY(Form("%s_%d",h->GetName(),i/rebinFactorX),i-rebinFactorX+1,i);
1083     if (proj->GetEntries() > 0) proj->Scale(1./proj->GetEntries());
1084     proj->Draw((i==rebinFactorX)?"hist":"histsames");
1085     proj->SetLineColor(i/rebinFactorX);
1086     if (f2) {
1087       f2->SetParameters(0.2,0.,1.,1.);
1088       f2->SetLineColor(i/rebinFactorX);
1089       proj->Fit("f2","WWNQ","sames");
1090       Double_t fwhm = f2->GetParameter(0);
1091       Double_t sigma = f2->GetParameter(3);
1092       Double_t sigmaP = TMath::Sqrt(sigma*sigma + fwhm*fwhm/(8.*log(2.)));
1093       Int_t rebin = TMath::Max(Int_t(0.5*sigmaP/proj->GetBinWidth(1)),1);
1094       while (proj->GetNbinsX()%rebin!=0) rebin--;
1095       proj->Rebin(rebin);
1096       proj->Scale(1./rebin);
1097       proj->Fit("f2","Q","sames");
1098     } else proj->SetLineWidth(2);
1099     Double_t p = 0.5 * (h->GetBinLowEdge(i-rebinFactorX+1) + h->GetBinLowEdge(i+1));
1100     l->AddEntry(proj,Form("%5.1f GeV",p));
1101   }
1102   if (f2) cout<<Form("\rFitting %s... %d/%d",fitting,nBins,nBins)<<endl;
1103   l->Draw("same");
1104   return c;
1105 }
1106