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