AliRoot Core  v5-06-15 (45dab64)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MUONRecoCheck.C
Go to the documentation of this file.
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 
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 {
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 {
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 {
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 {
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 {
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 {
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 {
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 {
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 
AliMUONConstants::AbsZEnd
static Double_t AbsZEnd()
Return z-position of absorber end.
Definition: AliMUONConstants.h:82
AliMUONVTrackStore
Base class of a track container.
Definition: AliMUONVTrackStore.h:21
MUONRecoCheck
void MUONRecoCheck(Int_t nEvent=-1, const char *pathSim="./generated/", const char *esdFileName="AliESDs.root", const char *ocdbPath="local://$ALICE_ROOT/OCDB", const Double_t pMin=0., const Double_t pMax=300., const Int_t pNBins=30, Int_t absorberRegion=-1, Bool_t fitResiduals=kTRUE)
Definition: MUONRecoCheck.C:73
printf
printf("Chi2/npoints = %f\n", TMath::Sqrt(chi2/npoints))
DrawVsAng
TCanvas * DrawVsAng(const char *name, const char *title, TH1 *h1, TH2 *h2)
Definition: MUONRecoCheck.C:1155
FitPDCAVsMom
void FitPDCAVsMom(TH2 *h, Int_t nBins, const char *fitting, TGraphAsymmErrors *gMean, TGraphAsymmErrors *gSigma)
Definition: MUONRecoCheck.C:1087
AliMpDEIterator
The iterator over detection elements.
Definition: AliMpDEIterator.h:29
AliMUONRecoCheck::ReconstructibleTracks
AliMUONVTrackStore * ReconstructibleTracks(Int_t event, UInt_t requestedStationMask=0x1F, Bool_t request2ChInSameSt45=kTRUE, Bool_t hitInFrontOfPad=kFALSE)
Return reconstructible reference tracks.
Definition: AliMUONRecoCheck.cxx:343
FitGausResVsMom
void FitGausResVsMom(TH2 *h, Int_t nBins, const Double_t mean0, const Double_t sigma0, const char *fitting, TGraphAsymmErrors *gMean, TGraphAsymmErrors *gSigma)
Definition: MUONRecoCheck.C:1056
AliMUONRecoParam::ImproveTracks
void ImproveTracks(Bool_t flag)
switch on/off the track improvement and keep the default cut in sigma to apply on cluster (local chi2...
Definition: AliMUONRecoParam.h:120
Zoom
void Zoom(TH1 *h, Double_t fractionCut=0.01)
Definition: MUONRecoCheck.C:1226
AliMUONRecoParam::MakeMoreTrackCandidates
void MakeMoreTrackCandidates(Bool_t flag)
switch on/off the building of track candidates starting from 1 cluster in each of the stations 4 and ...
Definition: AliMUONRecoParam.h:160
AliMUONTrackParam.h
AliMUONVTriggerTrackStore::CreateIterator
virtual TIterator * CreateIterator() const =0
Iterator to loop over tracks.
AliMUONRecoParam::GetSigmaCutForTracking
Double_t GetSigmaCutForTracking() const
return the cut in sigma to apply on cluster (local chi2) and track (global chi2) during tracking ...
Definition: AliMUONRecoParam.h:117
AliMUONRecoCheck.h
AliMUONRecoCheck::NumberOfEvents
Int_t NumberOfEvents() const
Return the total number of events.
Definition: AliMUONRecoCheck.cxx:218
AliMUONConstants::NTrackingCh
static Int_t NTrackingCh()
Return number of tracking chambers.
Definition: AliMUONConstants.cxx:134
DrawVsPos
TCanvas * DrawVsPos(const char *name, const char *title, TH2 *h1, TH2 *h2, TH2 *h3)
Definition: MUONRecoCheck.C:1174
AliMUONTrackParam::GetBendingCoor
Double_t GetBendingCoor() const
return bending coordinate (cm)
Definition: AliMUONTrackParam.h:45
AliMUONVTrackStore::CreateIterator
virtual TIterator * CreateIterator() const =0
Create an iterator to loop over tracks.
AliMUONESDInterface.h
AliMUONRecoParam::GetDefaultNonBendingReso
Double_t GetDefaultNonBendingReso(Int_t iCh) const
Get the default non bending resolution of chamber iCh.
Definition: AliMUONRecoParam.h:257
AliMUONESDInterface::ResetTracker
static void ResetTracker(const AliMUONRecoParam *recoParam=0x0, Bool_t info=kTRUE)
Definition: AliMUONESDInterface.cxx:436
AliMpDEIterator::IsDone
Bool_t IsDone() const
Definition: AliMpDEIterator.cxx:123
AliMUONRecoCheck::GetRunNumber
Int_t GetRunNumber()
Return the run number of the current ESD event.
Definition: AliMUONRecoCheck.cxx:203
AliMUONVStore::GetSize
virtual Int_t GetSize() const =0
The number of objects stored.
AliMUONTrackParam::GetZ
Double_t GetZ() const
return Z coordinate (cm)
Definition: AliMUONTrackParam.h:33
AliMUONTrackParam
Track parameters in ALICE dimuon spectrometer.
Definition: AliMUONTrackParam.h:22
AliMUONRecoParam.h
AliMpDEIterator::CurrentDEId
Int_t CurrentDEId() const
Definition: AliMpDEIterator.cxx:140
AliMUONRecoParam
Class with MUON reconstruction parameters.
Definition: AliMUONRecoParam.h:18
AliMUONTrackExtrap::ExtrapToZ
static Bool_t ExtrapToZ(AliMUONTrackParam *trackParam, Double_t zEnd)
Definition: AliMUONTrackExtrap.cxx:159
AliMUONTriggerTrack
Reconstructed trigger track in ALICE dimuon spectrometer.
Definition: AliMUONTriggerTrack.h:20
AliMUONRecoCheck::TriggeredTracks
AliMUONVTriggerTrackStore * TriggeredTracks(Int_t event)
Return the list of reconstructed trigger tracks.
Definition: AliMUONRecoCheck.cxx:257
request2ChInSameSt45
Bool_t request2ChInSameSt45
Definition: MUONFakes.C:38
AliMUONVTrackStore.h
AliMUONTriggerTrack::GetX11
Float_t GetX11() const
Return x position of fired Y strip in MC11.
Definition: AliMUONTriggerTrack.h:36
AliMUONRecoCheck
Utility class to check reconstruction.
Definition: AliMUONRecoCheck.h:31
AliMUONVCluster.h
AliMUONCDB::LoadField
Bool_t LoadField()
Definition: AliMUONCDB.cxx:233
DrawResMomVsMom
TCanvas * DrawResMomVsMom(const char *name, const char *title, TH2 *h, Int_t nBins, TF1 *f2=0x0, const char *fitting="")
Definition: MUONRecoCheck.C:1189
sum
void sum()
Definition: AliTPCclusterAnalysis.C:78
AliMUONTriggerTrack::GetSlopeY
Float_t GetSlopeY() const
Return track slope in Y.
Definition: AliMUONTriggerTrack.h:50
AliMUONTrackParam::GetBendingSlope
Double_t GetBendingSlope() const
return bending slope (cm ** -1)
Definition: AliMUONTrackParam.h:49
sigmaCut
Double_t sigmaCut
Definition: DIMUONFakes.C:36
AliMUONTrack.h
AliMUONVCluster
abstract base class for clusters
Definition: AliMUONVCluster.h:18
AliMUONTrackParam::Py
Double_t Py() const
Definition: AliMUONTrackParam.cxx:210
AliMpDEIterator.h
AliMUONTriggerTrack::GetY11
Float_t GetY11() const
Return y position of fired X strip in MC11.
Definition: AliMUONTriggerTrack.h:38
AliMUONCDB.h
AliMUONConstants.h
AliMUONRecoParam::GetDefaultBendingReso
Double_t GetDefaultBendingReso(Int_t iCh) const
Get the default bending resolution of chamber iCh.
Definition: AliMUONRecoParam.h:261
AliMUONTrackParam::P
Double_t P() const
Definition: AliMUONTrackParam.cxx:234
AliMUONVTriggerTrackStore
Base class of a trigger track store.
Definition: AliMUONVTriggerTrackStore.h:21
AliMUONTrackExtrap.h
AliMUONTrackParam::Px
Double_t Px() const
Definition: AliMUONTrackParam.cxx:196
AliMUONTriggerTrack.h
AliMUONRecoCheck::ReconstructedTracks
AliMUONVTrackStore * ReconstructedTracks(Int_t event, Bool_t refit=kTRUE)
Return the list of reconstructed tracks.
Definition: AliMUONRecoCheck.cxx:226
AliMUONVCluster::GetChamberId
static Int_t GetChamberId(UInt_t uniqueID)
Return chamber id (0..), part of the uniqueID.
Definition: AliMUONVCluster.h:56
AliMUONVCluster::GetY
virtual Double_t GetY() const =0
Return coordinate Y (cm)
requestedStationMask
UInt_t requestedStationMask
Definition: MUONFakes.C:37
AliMUONTrackExtrap::ExtrapToVertexWithoutBranson
static void ExtrapToVertexWithoutBranson(AliMUONTrackParam *trackParam, Double_t zVtx)
Definition: AliMUONTrackExtrap.cxx:902
AliMUONRecoParam::RequestStation
void RequestStation(Int_t iSt, Bool_t flag)
request or not at least one cluster in the station to validate the track
Definition: AliMUONRecoParam.h:185
AliMUONTrackParam::GetNonBendingCoor
Double_t GetNonBendingCoor() const
return non bending coordinate (cm)
Definition: AliMUONTrackParam.h:37
AliMUONRecoCheck::TriggerableTracks
AliMUONVTriggerTrackStore * TriggerableTracks(Int_t event)
Return triggerable reference tracks.
Definition: AliMUONRecoCheck.cxx:314
AliMUONTrackParam::GetClusterPtr
AliMUONVCluster * GetClusterPtr() const
get pointeur to associated cluster
Definition: AliMUONTrackParam.h:99
AliMUONVCluster::GetDetElemId
static Int_t GetDetElemId(UInt_t uniqueID)
Return detection element id, part of the uniqueID.
Definition: AliMUONVCluster.h:58
AliMUONVCluster::GetX
virtual Double_t GetX() const =0
Return coordinate X (cm)
AliMUONTrackParam::Pz
Double_t Pz() const
Definition: AliMUONTrackParam.cxx:224
AliMUONTrack
Reconstructed track in ALICE dimuon spectrometer.
Definition: AliMUONTrack.h:24
AliMUONRecoParam::GetSigmaCutForImprovement
Double_t GetSigmaCutForImprovement() const
return the cut in sigma to apply on cluster (local chi2) during track improvement ...
Definition: AliMUONRecoParam.h:126
proj
TEveProjectionManager * proj
Definition: tpc_tracks.C:9
dPhi
void dPhi()
Definition: CalibAlign.C:462
AliMUONCDB::LoadMapping
Bool_t LoadMapping(Bool_t segmentationOnly=kFALSE)
Definition: AliMUONCDB.cxx:257
FillResidual
void FillResidual(TH1 *h, Int_t i, Double_t &sigma, TGraphErrors *gMean, TGraphErrors *gSigma, Bool_t correctForSystematics, Bool_t fitResiduals)
Definition: MUONRecoCheck.C:1116
AliMUONVTriggerTrackStore.h
AliMUONCDB::LoadRecoParam
AliMUONRecoParam * LoadRecoParam()
Definition: AliMUONCDB.cxx:290
AliMUONTrack::GetTrackParamAtCluster
TObjArray * GetTrackParamAtCluster() const
Definition: AliMUONTrack.cxx:339
langaufun
Double_t langaufun(Double_t *x, Double_t *par)
Definition: MUONRecoCheck.C:999
AliMUONTrack::GetTrackParamAtVertex
AliMUONTrackParam * GetTrackParamAtVertex() const
return pointer to track parameters at vertex (can be 0x0)
Definition: AliMUONTrack.h:98
AliMUONTrackParam::GetNonBendingSlope
Double_t GetNonBendingSlope() const
return non bending slope (cm ** -1)
Definition: AliMUONTrackParam.h:41