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