AliPhysics  3f43add (3f43add)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskMuonResolution.cxx
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 
18 // ROOT includes
19 #include <TROOT.h>
20 #include <TFile.h>
21 #include <TTree.h>
22 #include <TH1F.h>
23 #include <TH2F.h>
24 #include <TMultiGraph.h>
25 #include <TGraphErrors.h>
26 #include <TCanvas.h>
27 #include <TLegend.h>
28 #include <Riostream.h>
29 #include <TString.h>
30 #include <TGeoManager.h>
31 #include <TList.h>
32 #include <TObjString.h>
33 #include <TRegexp.h>
34 
35 // STEER includes
36 #include "AliESDEvent.h"
37 #include "AliESDMuonTrack.h"
38 #include "AliCDBManager.h"
39 #include "AliCDBStorage.h"
40 #include "AliGeomManager.h"
41 #include "AliVVertex.h"
42 
43 // ANALYSIS includes
44 #include "AliAnalysisDataSlot.h"
45 #include "AliAnalysisManager.h"
46 #include "AliInputEventHandler.h"
48 #include "AliCentrality.h"
49 
50 // MUON includes
51 #include "AliMUONCDB.h"
52 #include "AliMUONConstants.h"
53 #include "AliMUONRecoParam.h"
54 #include "AliMUONESDInterface.h"
55 #include "AliMUONVTrackReconstructor.h"
56 #include "AliMUONTrack.h"
57 #include "AliMUONTrackParam.h"
58 #include "AliMUONTrackExtrap.h"
59 #include "AliMUONVCluster.h"
60 #include "AliMUONVDigit.h"
61 #include "AliMUONGeometryTransformer.h"
62 #include "AliMUONGeometryModuleTransformer.h"
63 #include "AliMUONGeometryDetElement.h"
64 #include "AliMpDEIterator.h"
65 #include "AliMpSegmentation.h"
66 #include "AliMpVSegmentation.h"
67 #include "AliMpConstants.h"
68 #include "AliMpDDLStore.h"
69 #include "AliMpPad.h"
70 #include "AliMpDetElement.h"
71 #include "AliMpCathodType.h"
72 
73 #ifndef SafeDelete
74 #define SafeDelete(x) if (x != NULL) { delete x; x = NULL; }
75 #endif
76 
77 using std::cout;
78 using std::endl;
79 using std::flush;
80 
82 
83 const Int_t AliAnalysisTaskMuonResolution::fgkMinEntries = 10;
84 
85 //________________________________________________________________________
87  AliAnalysisTaskSE(),
88  fResiduals(NULL),
89  fResidualsVsP(NULL),
90  fResidualsVsCent(NULL),
91  fResidualsVsAngle(NULL),
92  fLocalChi2(NULL),
93  fChamberRes(NULL),
94  fTrackRes(NULL),
95  fCanvases(NULL),
96  fTmpHists(NULL),
97  fDefaultStorage(""),
98  fNEvents(0),
99  fShowProgressBar(kFALSE),
100  fPrintClResPerCh(kFALSE),
101  fPrintClResPerDE(kFALSE),
102  fGaus(NULL),
103  fMinMomentum(0.),
104  fMinPt(0.),
105  fExtrapMode(1),
106  fCorrectForSystematics(kTRUE),
107  fRemoveMonoCathCl(kFALSE),
108  fCheckAllPads(kFALSE),
109  fImproveTracks(kFALSE),
110  fShiftHalfCh(kFALSE),
111  fPrintHalfChShift(kFALSE),
112  fShiftDE(kFALSE),
113  fPrintDEShift(kFALSE),
114  fOCDBLoaded(kFALSE),
115  fNDE(0),
116  fReAlign(kFALSE),
117  fOldAlignStorage(""),
118  fOldAlignVersion(-1),
119  fOldAlignSubVersion(-1),
120  fNewAlignStorage(""),
121  fNewAlignVersion(-1),
122  fNewAlignSubVersion(-1),
123  fOldGeoTransformer(NULL),
124  fNewGeoTransformer(NULL),
125  fMuonEventCuts(0x0),
126  fMuonTrackCuts(0x0)
127 {
129 
130  for (Int_t i = 0; i < 10; i++) SetStartingResolution(i, -1., -1.);
131  for (Int_t i = 0; i < 1100; i++) fDEIndices[i] = 0;
132  for (Int_t i = 0; i < 200; i++) fDEIds[i] = 0;
133  for (Int_t i = 0; i < 20; i++) SetHalfChShift(i, 0., 0.);
134  for (Int_t i = 0; i < 200; i++) SetDEShift(i, 0., 0.);
135 
136 }
137 
138 //________________________________________________________________________
140  AliAnalysisTaskSE(name),
141  fResiduals(NULL),
142  fResidualsVsP(NULL),
143  fResidualsVsCent(NULL),
144  fResidualsVsAngle(NULL),
145  fLocalChi2(NULL),
146  fChamberRes(NULL),
147  fTrackRes(NULL),
148  fCanvases(NULL),
149  fTmpHists(NULL),
150  fDefaultStorage("raw://"),
151  fNEvents(0),
152  fShowProgressBar(kFALSE),
153  fPrintClResPerCh(kFALSE),
154  fPrintClResPerDE(kFALSE),
155  fGaus(NULL),
156  fMinMomentum(0.),
157  fMinPt(0.),
158  fExtrapMode(1),
159  fCorrectForSystematics(kTRUE),
160  fRemoveMonoCathCl(kFALSE),
161  fCheckAllPads(kFALSE),
162  fImproveTracks(kFALSE),
163  fShiftHalfCh(kFALSE),
164  fPrintHalfChShift(kFALSE),
165  fShiftDE(kFALSE),
166  fPrintDEShift(kFALSE),
167  fOCDBLoaded(kFALSE),
168  fNDE(0),
169  fReAlign(kFALSE),
170  fOldAlignStorage(""),
171  fOldAlignVersion(-1),
172  fOldAlignSubVersion(-1),
173  fNewAlignStorage(""),
174  fNewAlignVersion(-1),
175  fNewAlignSubVersion(-1),
176  fOldGeoTransformer(NULL),
177  fNewGeoTransformer(NULL),
178  fMuonEventCuts(0x0),
179  fMuonTrackCuts(0x0)
180 {
182 
183  for (Int_t i = 0; i < 10; i++) SetStartingResolution(i, -1., -1.);
184  for (Int_t i = 0; i < 1100; i++) fDEIndices[i] = 0;
185  for (Int_t i = 0; i < 200; i++) fDEIds[i] = 0;
186  for (Int_t i = 0; i < 20; i++) SetHalfChShift(i, 0., 0.);
187  for (Int_t i = 0; i < 200; i++) SetDEShift(i, 0., 0.);
188 
189  FitResiduals();
190 
191  // Output slot #1 writes into a TObjArray container
192  DefineOutput(1,TObjArray::Class());
193  // Output slot #2 writes into a TObjArray container
194  DefineOutput(2,TObjArray::Class());
195  // Output slot #3 writes into a TObjArray container
196  DefineOutput(3,TObjArray::Class());
197  // Output slot #4 writes into a TObjArray container
198  DefineOutput(4,TObjArray::Class());
199  // Output slot #5 writes into a TObjArray container
200  DefineOutput(5,TObjArray::Class());
201  // Output slot #6 writes into a TObjArray container
202  DefineOutput(6,TObjArray::Class());
203  // Output slot #7 writes into a TObjArray container
204  DefineOutput(7,TObjArray::Class());
205 }
206 
207 //________________________________________________________________________
209 {
211  if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
217  }
222  SafeDelete(fGaus);
227 }
228 
229 //___________________________________________________________________________
231 {
233 
234  // do it once the OCDB has been loaded (i.e. from NotifyRun())
235  if (!fOCDBLoaded) return;
236 
237  fResiduals = new TObjArray(1000);
238  fResiduals->SetOwner();
239  fResidualsVsP = new TObjArray(1000);
240  fResidualsVsP->SetOwner();
241  fResidualsVsCent = new TObjArray(1000);
242  fResidualsVsCent->SetOwner();
243  fResidualsVsAngle = new TObjArray(1000);
244  fResidualsVsAngle->SetOwner();
245  fTrackRes = new TObjArray(1000);
246  fTrackRes->SetOwner();
247  TH2F* h2;
248 
249  // find the highest chamber resolution and set histogram bins
250  const AliMUONRecoParam* recoParam = AliMUONESDInterface::GetTracker()->GetRecoParam();
251  Double_t maxSigma[2] = {-1., -1.};
252  for (Int_t i = 0; i < 10; i++) {
253  if (recoParam->GetDefaultNonBendingReso(i) > maxSigma[0]) maxSigma[0] = recoParam->GetDefaultNonBendingReso(i);
254  if (recoParam->GetDefaultBendingReso(i) > maxSigma[1]) maxSigma[1] = recoParam->GetDefaultBendingReso(i);
255  }
256  const char* axes[2] = {"X", "Y"};
257  const char* side[2] = {"I", "O"};
258  const Int_t nBins = 5000;
259  const Int_t nSigma = 10;
260  const Int_t pNBins = 20;
261  const Double_t pEdges[2] = {0., 50.};
262  Int_t nCentBins = 12;
263  Double_t centRange[2] = {-10., 110.};
264  Int_t nAngleBins = 20;
265  Double_t angleRange[2][2] = {{-15., 15.}, {-40., 40.}};
266  TString name, title;
267 
268  for (Int_t ia = 0; ia < 2; ia++) {
269 
270  Double_t maxRes = nSigma*maxSigma[ia];
271 
272  // List of residual histos
273  name = Form("hResidual%sPerCh_ClusterIn",axes[ia]); title = Form("cluster-track residual-%s distribution per chamber (cluster attached to the track);chamber ID;#Delta_{%s} (cm)",axes[ia],axes[ia]);
274  h2 = new TH2F(name.Data(), title.Data(), 10, 0.5, 10.5, nBins, -maxRes, maxRes);
275  fResiduals->AddAtAndExpand(h2, kResidualPerChClusterIn+ia);
276  name = Form("hResidual%sPerCh_ClusterOut",axes[ia]); title = Form("cluster-track residual-%s distribution per chamber (cluster not attached to the track);chamber ID;#Delta_{%s} (cm)",axes[ia],axes[ia]);
277  h2 = new TH2F(name.Data(), title.Data(), 10, 0.5, 10.5, nBins, -2.*maxRes, 2.*maxRes);
278  fResiduals->AddAtAndExpand(h2, kResidualPerChClusterOut+ia);
279 
280  name = Form("hResidual%sPerHalfCh_ClusterIn",axes[ia]); title = Form("cluster-track residual-%s distribution per half chamber (cluster attached to the track);half chamber ID;#Delta_{%s} (cm)",axes[ia],axes[ia]);
281  h2 = new TH2F(name.Data(), title.Data(), 20, 0.5, 20.5, nBins, -maxRes, maxRes);
282  for (Int_t i = 0; i < 10; i++) { h2->GetXaxis()->SetBinLabel(2*i+1, Form("%d-I",i+1)); h2->GetXaxis()->SetBinLabel(2*i+2, Form("%d-O",i+1)); }
283  fResiduals->AddAtAndExpand(h2, kResidualPerHalfChClusterIn+ia);
284  name = Form("hResidual%sPerHalfCh_ClusterOut",axes[ia]); title = Form("cluster-track residual-%s distribution per half chamber (cluster not attached to the track);half chamber ID;#Delta_{%s} (cm)",axes[ia],axes[ia]);
285  h2 = new TH2F(name.Data(), title.Data(), 20, 0.5, 20.5, nBins, -2.*maxRes, 2.*maxRes);
286  for (Int_t i = 0; i < 10; i++) { h2->GetXaxis()->SetBinLabel(2*i+1, Form("%d-I",i+1)); h2->GetXaxis()->SetBinLabel(2*i+2, Form("%d-O",i+1)); }
287  fResiduals->AddAtAndExpand(h2, kResidualPerHalfChClusterOut+ia);
288 
289  name = Form("hResidual%sPerDE_ClusterIn",axes[ia]); title = Form("cluster-track residual-%s distribution per DE (cluster attached to the track);DE ID;#Delta_{%s} (cm)",axes[ia],axes[ia]);
290  h2 = new TH2F(name.Data(), title.Data(), fNDE, 0.5, fNDE+0.5, nBins, -maxRes, maxRes);
291  for (Int_t i = 1; i <= fNDE; i++) h2->GetXaxis()->SetBinLabel(i, Form("%d",fDEIds[i]));
292  fResiduals->AddAtAndExpand(h2, kResidualPerDEClusterIn+ia);
293  name = Form("hResidual%sPerDE_ClusterOut",axes[ia]); title = Form("cluster-track residual-%s distribution per DE (cluster not attached to the track);DE ID;#Delta_{%s} (cm)",axes[ia],axes[ia]);
294  h2 = new TH2F(name.Data(), title.Data(), fNDE, 0.5, fNDE+0.5, nBins, -2.*maxRes, 2.*maxRes);
295  for (Int_t i = 1; i <= fNDE; i++) h2->GetXaxis()->SetBinLabel(i, Form("%d",fDEIds[i]));
296  fResiduals->AddAtAndExpand(h2, kResidualPerDEClusterOut+ia);
297 
298  name = Form("hTrackRes%sPerCh",axes[ia]); title = Form("track #sigma_{%s} per Ch;chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]);
299  h2 = new TH2F(name.Data(), title.Data(), 10, 0.5, 10.5, nBins, 0., maxRes);
300  fResiduals->AddAtAndExpand(h2, kTrackResPerCh+ia);
301  name = Form("hTrackRes%sPerHalfCh",axes[ia]); title = Form("track #sigma_{%s} per half Ch;half chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]);
302  h2 = new TH2F(name.Data(), title.Data(), 20, 0.5, 20.5, nBins, 0., maxRes);
303  for (Int_t i = 0; i < 10; i++) { h2->GetXaxis()->SetBinLabel(2*i+1, Form("%d-I",i+1)); h2->GetXaxis()->SetBinLabel(2*i+2, Form("%d-O",i+1)); }
304  fResiduals->AddAtAndExpand(h2, kTrackResPerHalfCh+ia);
305  name = Form("hTrackRes%sPerDE",axes[ia]); title = Form("track #sigma_{%s} per DE;DE ID;#sigma_{%s} (cm)",axes[ia],axes[ia]);
306  h2 = new TH2F(name.Data(), title.Data(), fNDE, 0.5, fNDE+0.5, nBins, 0., maxRes);
307  for (Int_t i = 1; i <= fNDE; i++) h2->GetXaxis()->SetBinLabel(i, Form("%d",fDEIds[i]));
308  fResiduals->AddAtAndExpand(h2, kTrackResPerDE+ia);
309 
310  name = Form("hMCS%sPerCh",axes[ia]); title = Form("MCS %s-dispersion of extrapolated track per Ch;chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]);
311  h2 = new TH2F(name.Data(), title.Data(), 10, 0.5, 10.5, nBins, 0., 0.2);
312  fResiduals->AddAtAndExpand(h2, kMCSPerCh+ia);
313  name = Form("hMCS%sPerHalfCh",axes[ia]); title = Form("MCS %s-dispersion of extrapolated track per half Ch;half chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]);
314  h2 = new TH2F(name.Data(), title.Data(), 20, 0.5, 20.5, nBins, 0., 0.2);
315  for (Int_t i = 0; i < 10; i++) { h2->GetXaxis()->SetBinLabel(2*i+1, Form("%d-I",i+1)); h2->GetXaxis()->SetBinLabel(2*i+2, Form("%d-O",i+1)); }
316  fResiduals->AddAtAndExpand(h2, kMCSPerHalfCh+ia);
317  name = Form("hMCS%sPerDE",axes[ia]); title = Form("MCS %s-dispersion of extrapolated track per DE;DE ID;#sigma_{%s} (cm)",axes[ia],axes[ia]);
318  h2 = new TH2F(name.Data(), title.Data(), fNDE, 0.5, fNDE+0.5, nBins, 0., 0.2);
319  for (Int_t i = 1; i <= fNDE; i++) h2->GetXaxis()->SetBinLabel(i, Form("%d",fDEIds[i]));
320  fResiduals->AddAtAndExpand(h2, kMCSPerDE+ia);
321 
322  name = Form("hClusterRes2%sPerCh",axes[ia]); title = Form("cluster #sigma_{%s}^{2} per Ch;chamber ID;#sigma_{%s}^{2} (cm^{2})",axes[ia],axes[ia]);
323  h2 = new TH2F(name.Data(), title.Data(), 10, 0.5, 10.5, nSigma*nBins, -0.1*maxRes*maxRes, maxRes*maxRes);
324  fResiduals->AddAtAndExpand(h2, kClusterRes2PerCh+ia);
325 
326  // List of residual histos vs. p or vs. centrality or vs. angle
327  for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
328  name = Form("hResidual%sInCh%dVsP_ClusterIn",axes[ia],i+1); title = Form("cluster-track residual-%s distribution in chamber %d versus momentum (cluster attached to the track);p (GeV/c^{2});#Delta_{%s} (cm)",axes[ia],i+1,axes[ia]);
329  h2 = new TH2F(name.Data(), title.Data(), pNBins, pEdges[0], pEdges[1], nBins, -maxRes, maxRes);
330  fResidualsVsP->AddAtAndExpand(h2, kResidualInChVsPClusterIn+10*ia+i);
331  name = Form("hResidual%sInCh%dVsP_ClusterOut",axes[ia],i+1); title = Form("cluster-track residual-%s distribution in chamber %d versus momentum (cluster not attached to the track);p (GeV/c^{2});#Delta_{%s} (cm)",axes[ia],i+1,axes[ia]);
332  h2 = new TH2F(name.Data(), title.Data(), pNBins, pEdges[0], pEdges[1], nBins, -2.*maxRes, 2.*maxRes);
333  fResidualsVsP->AddAtAndExpand(h2, kResidualInChVsPClusterOut+10*ia+i);
334 
335  name = Form("hResidual%sInCh%dVsCent_ClusterIn",axes[ia],i+1); title = Form("cluster-track residual-%s distribution in chamber %d versus centrality (cluster attached to the track);centrality (%%);#Delta_{%s} (cm)",axes[ia],i+1,axes[ia]);
336  h2 = new TH2F(name.Data(), title.Data(), nCentBins, centRange[0], centRange[1], nBins, -maxRes, maxRes);
337  fResidualsVsCent->AddAtAndExpand(h2, kResidualInChVsCentClusterIn+10*ia+i);
338  name = Form("hResidual%sInCh%dVsCent_ClusterOut",axes[ia],i+1); title = Form("cluster-track residual-%s distribution in chamber %d versus centrality (cluster not attached to the track);centrality (%%);#Delta_{%s} (cm)",axes[ia],i+1,axes[ia]);
339  h2 = new TH2F(name.Data(), title.Data(), nCentBins, centRange[0], centRange[1], nBins, -2.*maxRes, 2.*maxRes);
340  fResidualsVsCent->AddAtAndExpand(h2, kResidualInChVsCentClusterOut+10*ia+i);
341 
342  name = Form("hResidual%sInCh%dVsAngle_ClusterIn",axes[ia],i+1); title = Form("cluster-track residual-%s distribution in chamber %d versus track angle (cluster attached to the track);%s-angle (deg);#Delta_{%s} (cm)",axes[ia],i+1,axes[ia],axes[ia]);
343  h2 = new TH2F(name.Data(), title.Data(), nAngleBins, angleRange[ia][0], angleRange[ia][1], nBins, -maxRes, maxRes);
344  fResidualsVsAngle->AddAtAndExpand(h2, kResidualInChVsAngleClusterIn+10*ia+i);
345  name = Form("hResidual%sInCh%dVsAngle_ClusterOut",axes[ia],i+1); title = Form("cluster-track residual-%s distribution in chamber %d versus track angle (cluster not attached to the track);%s-angle (deg);#Delta_{%s} (cm)",axes[ia],i+1,axes[ia],axes[ia]);
346  h2 = new TH2F(name.Data(), title.Data(), nAngleBins, angleRange[ia][0], angleRange[ia][1], nBins, -2.*maxRes, 2.*maxRes);
347  fResidualsVsAngle->AddAtAndExpand(h2, kResidualInChVsAngleClusterOut+10*ia+i);
348 
349  for (Int_t j = 0; j < 2; j++) {
350  name = Form("hResidual%sInHalfCh%d%sVsAngle_ClusterIn",axes[ia],i+1,side[j]); title = Form("cluster-track residual-%s distribution in half-chamber %d-%s versus track angle (cluster attached to the track);%s-angle (deg);#Delta_{%s} (cm)",axes[ia],i+1,side[j],axes[ia],axes[ia]);
351  h2 = new TH2F(name.Data(), title.Data(), nAngleBins, angleRange[ia][0], angleRange[ia][1], nBins, -maxRes, maxRes);
352  fResidualsVsAngle->AddAtAndExpand(h2, kResidualInHalfChVsAngleClusterIn+20*ia+2*i+j);
353  }
354  }
355 
356  name = Form("hResidual%sVsP_ClusterIn",axes[ia]); title = Form("cluster-track residual-%s distribution integrated over chambers versus momentum (cluster attached to the track);p (GeV/c^{2});#Delta_{%s} (cm)",axes[ia],axes[ia]);
357  h2 = new TH2F(name.Data(), title.Data(), pNBins, pEdges[0], pEdges[1], nBins, -maxRes, maxRes);
358  fResidualsVsP->AddAtAndExpand(h2, kResidualVsPClusterIn+ia);
359  name = Form("hResidual%sVsP_ClusterOut",axes[ia]); title = Form("cluster-track residual-%s distribution integrated over chambers versus momentum (cluster not attached to the track);p (GeV/c^{2});#Delta_{%s} (cm)",axes[ia],axes[ia]);
360  h2 = new TH2F(name.Data(), title.Data(), pNBins, pEdges[0], pEdges[1], nBins, -2.*maxRes, 2.*maxRes);
361  fResidualsVsP->AddAtAndExpand(h2, kResidualVsPClusterOut+ia);
362 
363  name = Form("hResidual%sVsCent_ClusterIn",axes[ia]); title = Form("cluster-track residual-%s distribution integrated over chambers versus centrality (cluster attached to the track);centrality (%%);#Delta_{%s} (cm)",axes[ia],axes[ia]);
364  h2 = new TH2F(name.Data(), title.Data(), nCentBins, centRange[0], centRange[1], nBins, -maxRes, maxRes);
365  fResidualsVsCent->AddAtAndExpand(h2, kResidualVsCentClusterIn+ia);
366  name = Form("hResidual%sVsCent_ClusterOut",axes[ia]); title = Form("cluster-track residual-%s distribution integrated over chambers versus centrality (cluster not attached to the track);centrality (%%);#Delta_{%s} (cm)",axes[ia],axes[ia]);
367  h2 = new TH2F(name.Data(), title.Data(), nCentBins, centRange[0], centRange[1], nBins, -2.*maxRes, 2.*maxRes);
368  fResidualsVsCent->AddAtAndExpand(h2, kResidualVsCentClusterOut+ia);
369 
370  name = Form("hResidual%sVsAngle_ClusterIn",axes[ia]); title = Form("cluster-track residual-%s distribution integrated over chambers versus track angle (cluster attached to the track);%s-angle (deg);#Delta_{%s} (cm)",axes[ia],axes[ia],axes[ia]);
371  h2 = new TH2F(name.Data(), title.Data(), nAngleBins, angleRange[ia][0], angleRange[ia][1], nBins, -maxRes, maxRes);
372  fResidualsVsAngle->AddAtAndExpand(h2, kResidualVsAngleClusterIn+ia);
373  name = Form("hResidual%sVsAngle_ClusterOut",axes[ia]); title = Form("cluster-track residual-%s distribution integrated over chambers versus track angle (cluster not attached to the track);%s-angle (deg);#Delta_{%s} (cm)",axes[ia],axes[ia],axes[ia]);
374  h2 = new TH2F(name.Data(), title.Data(), nAngleBins, angleRange[ia][0], angleRange[ia][1], nBins, -2.*maxRes, 2.*maxRes);
375  fResidualsVsAngle->AddAtAndExpand(h2, kResidualVsAngleClusterOut+ia);
376 
377  // local chi2
378  name = Form("hLocalChi2%sPerCh",axes[ia]); title = Form("local chi2-%s distribution per chamber;chamber ID;local #chi^{2}_{%s}", axes[ia], axes[ia]);
379  h2 = new TH2F(name.Data(), title.Data(), 10, 0.5, 10.5, 1000, 0., 25.);
380  fResiduals->AddAtAndExpand(h2, kLocalChi2PerCh+ia);
381  name = Form("hLocalChi2%sPerDE",axes[ia]); title = Form("local chi2-%s distribution per DE;DE ID;local #chi^{2}_{%s}", axes[ia], axes[ia]);
382  h2 = new TH2F(name.Data(), title.Data(), fNDE, 0.5, fNDE+0.5, 1000, 0., 25.);
383  for (Int_t i = 1; i <= fNDE; i++) h2->GetXaxis()->SetBinLabel(i, Form("%d",fDEIds[i]));
384  fResiduals->AddAtAndExpand(h2, kLocalChi2PerDE+ia);
385 
386  // track resolution
387  name = Form("hUncorrSlope%sRes",axes[ia]); title = Form("muon slope_{%s} reconstructed resolution at first cluster vs p;p (GeV/c); #sigma_{slope_{%s}}", axes[ia], axes[ia]);
388  h2 = new TH2F(name.Data(), title.Data(), 300, 0., 300., 1000, 0., 0.003);
389  fTrackRes->AddAtAndExpand(h2, kUncorrSlopeRes+ia);
390  name = Form("hSlope%sRes",axes[ia]); title = Form("muon slope_{%s} reconstructed resolution at vertex vs p;p (GeV/c); #sigma_{slope_{%s}}", axes[ia], axes[ia]);
391  h2 = new TH2F(name.Data(), title.Data(), 300, 0., 300., 1000, 0., 0.02);
392  fTrackRes->AddAtAndExpand(h2, kSlopeRes+ia);
393  }
394 
395  // local chi2 X+Y
396  h2 = new TH2F("hLocalChi2PerCh", "local chi2 (~0.5*(#chi^{2}_{X}+#chi^{2}_{Y})) distribution per chamber;chamber ID;local #chi^{2}", 10, 0.5, 10.5, 1000, 0., 25.);
397  fResiduals->AddAtAndExpand(h2, kLocalChi2PerCh+2);
398  h2 = new TH2F("hLocalChi2PerDE", "local chi2 (~0.5*(#chi^{2}_{X}+#chi^{2}_{Y})) distribution per chamber;DE ID;local #chi^{2}", fNDE, 0.5, fNDE+0.5, 1000, 0., 25.);
399  for (Int_t i = 1; i <= fNDE; i++) h2->GetXaxis()->SetBinLabel(i, Form("%d",fDEIds[i]));
400  fResiduals->AddAtAndExpand(h2, kLocalChi2PerDE+2);
401 
402  // track resolution
403  h2 = new TH2F("hUncorrPRes", "muon momentum reconstructed resolution at first cluster vs p;p (GeV/c); #sigma_{p}/p (%)", 300, 0., 300., 1000, 0., 10.);
404  fTrackRes->AddAtAndExpand(h2, kUncorrPRes);
405  h2 = new TH2F("hPRes", "muon momentum reconstructed resolution at vertex vs p;p (GeV/c); #sigma_{p}/p (%)", 300, 0., 300., 1000, 0., 10.);
406  fTrackRes->AddAtAndExpand(h2, kPRes);
407  h2 = new TH2F("hUncorrPtRes", "muon transverse momentum reconstructed resolution at first cluster vs p_{t};p_{t} (GeV/c); #sigma_{p_{t}}/p_{t} (%)", 300, 0., 30., 300, 0., 30.);
408  fTrackRes->AddAtAndExpand(h2, kUncorrPtRes);
409  h2 = new TH2F("hPtRes", "muon transverse momentum reconstructed resolution at vertex vs p_{t};p_{t} (GeV/c); #sigma_{p_{t}}/p_{t} (%)", 300, 0., 30., 300, 0., 30.);
410  fTrackRes->AddAtAndExpand(h2, kPtRes);
411 
412  // Post data at least once per task to ensure data synchronisation (required for merging)
413  PostData(1, fResiduals);
414  PostData(2, fResidualsVsP);
415  PostData(3, fResidualsVsCent);
416  PostData(4, fResidualsVsAngle);
417  PostData(5, fTrackRes);
418 }
419 
420 //________________________________________________________________________
422 {
424 
425  // check if OCDB properly loaded
426  if (!fOCDBLoaded) AliFatal("Problem occur while loading OCDB objects");
427 
428  AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
429  if (!esd) return;
430 
431  if (fShowProgressBar && (++fNEvents)%100 == 0) cout<<"\rEvent processing... "<<fNEvents<<"\r"<<flush;
432 
433  // skip events that do not pass required selections
434  if (fMuonEventCuts && !fMuonEventCuts->IsSelected(fInputHandler)) return;
435 
436  // get the centrality
437  Float_t centrality = esd->GetCentrality()->GetCentralityPercentileUnchecked("V0M");
438 
439  // get tracker to refit
440  AliMUONVTrackReconstructor* tracker = AliMUONESDInterface::GetTracker();
441 
442  // loop over tracks
443  Int_t nTracks = (Int_t) esd->GetNumberOfMuonTracks();
444  for (Int_t iTrack = 0; iTrack < nTracks; ++iTrack) {
445 
446  // get the ESD track
447  AliESDMuonTrack* esdTrack = esd->GetMuonTrack(iTrack);
448 
449  // skip ghost tracks
450  if (!esdTrack->ContainTrackerData()) continue;
451 
452  // apply standard track cuts if any
453  if (fMuonTrackCuts && !fMuonTrackCuts->IsSelected(esdTrack)) continue;
454 
455  // select positive or negative tracks
456  //if (esdTrack->Charge() < 0) continue;
457 
458  // skip tracks with not enough clusters
459  //if (esdTrack->GetNClusters() < 8) continue;
460 
461  // get the corresponding MUON track
462  AliMUONTrack track;
463  AliMUONESDInterface::ESDToMUON(*esdTrack, track, kFALSE);
464 
465  // change the cluster resolution (and position)
466  ModifyClusters(track);
467 
468  // refit the track
469  if (!tracker->RefitTrack(track, fImproveTracks) || (fImproveTracks && !track.IsImproved())) continue;
470 
471  // save track unchanged
472  AliMUONTrack referenceTrack(track);
473 
474  // get track param at first cluster and add MCS in first chamber
475  AliMUONTrackParam trackParamAtFirstCluster(*(static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->First())));
476  Int_t firstCh = trackParamAtFirstCluster.GetClusterPtr()->GetChamberId();
477  AliMUONTrackExtrap::AddMCSEffect(&trackParamAtFirstCluster, AliMUONConstants::ChamberThicknessInX0(firstCh)/2., -1.);
478 
479  // fill momentum error at first cluster
480  Double_t pXUncorr = trackParamAtFirstCluster.Px();
481  Double_t pYUncorr = trackParamAtFirstCluster.Py();
482  Double_t pZUncorr = trackParamAtFirstCluster.Pz();
483  Double_t pUncorr = trackParamAtFirstCluster.P();
484  TMatrixD covUncorr(5,5);
485  Cov2CovP(trackParamAtFirstCluster,covUncorr);
486  Double_t sigmaPUncorr = TMath::Sqrt(pXUncorr * (pXUncorr*covUncorr(2,2) + pYUncorr*covUncorr(2,3) + pZUncorr*covUncorr(2,4)) +
487  pYUncorr * (pXUncorr*covUncorr(3,2) + pYUncorr*covUncorr(3,3) + pZUncorr*covUncorr(3,4)) +
488  pZUncorr * (pXUncorr*covUncorr(4,2) + pYUncorr*covUncorr(4,3) + pZUncorr*covUncorr(4,4))) / pUncorr;
489  ((TH2F*)fTrackRes->UncheckedAt(kUncorrPRes))->Fill(pUncorr,100.*sigmaPUncorr/pUncorr);
490 
491  // fill transverse momentum error at first cluster
492  Double_t ptUncorr = TMath::Sqrt(pXUncorr*pXUncorr + pYUncorr*pYUncorr);
493  Double_t sigmaPtUncorr = TMath::Sqrt(pXUncorr * (pXUncorr*covUncorr(2,2) + pYUncorr*covUncorr(2,3)) + pYUncorr * (pXUncorr*covUncorr(3,2) + pYUncorr*covUncorr(3,3))) / ptUncorr;
494  ((TH2F*)fTrackRes->UncheckedAt(kUncorrPtRes))->Fill(ptUncorr,100.*sigmaPtUncorr/ptUncorr);
495 
496  // fill slopeX-Y error at first cluster
497  ((TH2F*)fTrackRes->UncheckedAt(kUncorrSlopeRes))->Fill(pUncorr,TMath::Sqrt(trackParamAtFirstCluster.GetCovariances()(1,1)));
498  ((TH2F*)fTrackRes->UncheckedAt(kUncorrSlopeRes+1))->Fill(pUncorr,TMath::Sqrt(trackParamAtFirstCluster.GetCovariances()(3,3)));
499 
500  // fill momentum error at vertex
501  AliMUONTrackParam trackParamAtVtx(trackParamAtFirstCluster);
502  AliMUONTrackExtrap::ExtrapToVertex(&trackParamAtVtx, esdTrack->GetNonBendingCoor(), esdTrack->GetBendingCoor(), esdTrack->GetZ(), 0., 0.);
503  Double_t pXVtx = trackParamAtVtx.Px();
504  Double_t pYVtx = trackParamAtVtx.Py();
505  Double_t pZVtx = trackParamAtVtx.Pz();
506  Double_t pVtx = trackParamAtVtx.P();
507  TMatrixD covVtx(5,5);
508  Cov2CovP(trackParamAtVtx,covVtx);
509  Double_t sigmaPVtx = TMath::Sqrt(pXVtx * (pXVtx*covVtx(2,2) + pYVtx*covVtx(2,3) + pZVtx*covVtx(2,4)) +
510  pYVtx * (pXVtx*covVtx(3,2) + pYVtx*covVtx(3,3) + pZVtx*covVtx(3,4)) +
511  pZVtx * (pXVtx*covVtx(4,2) + pYVtx*covVtx(4,3) + pZVtx*covVtx(4,4))) / pVtx;
512  ((TH2F*)fTrackRes->UncheckedAt(kPRes))->Fill(pVtx,100.*sigmaPVtx/pVtx);
513 
514  // fill transverse momentum error at vertex
515  Double_t ptVtx = TMath::Sqrt(pXVtx*pXVtx + pYVtx*pYVtx);
516  Double_t sigmaPtVtx = TMath::Sqrt(pXVtx * (pXVtx*covVtx(2,2) + pYVtx*covVtx(2,3)) + pYVtx * (pXVtx*covVtx(3,2) + pYVtx*covVtx(3,3))) / ptVtx;
517  ((TH2F*)fTrackRes->UncheckedAt(kPtRes))->Fill(ptVtx,100.*sigmaPtVtx/ptVtx);
518 
519  // fill slopeX-Y error at vertex
520  ((TH2F*)fTrackRes->UncheckedAt(kSlopeRes))->Fill(pVtx,TMath::Sqrt(trackParamAtVtx.GetCovariances()(1,1)));
521  ((TH2F*)fTrackRes->UncheckedAt(kSlopeRes+1))->Fill(pVtx,TMath::Sqrt(trackParamAtVtx.GetCovariances()(3,3)));
522 
523  // skip low momentum tracks in the computation of resolution (cut on initial momentum to avoid correlations)
524  if (esdTrack->PUncorrected() < fMinMomentum) continue;
525 
526  // skip low pT tracks in the computation of resolution (cut on initial pT to avoid correlations)
527  Double_t pXUncorr0 = esdTrack->PxUncorrected();
528  Double_t pYUncorr0 = esdTrack->PyUncorrected();
529  if (TMath::Sqrt(pXUncorr0*pXUncorr0+pYUncorr0*pYUncorr0) < fMinPt) continue;
530 
531  // loop over clusters
532  Int_t nClusters = track.GetNClusters();
533  for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
534 
535  // Get current, previous and next trackParams
536  AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->UncheckedAt(iCluster));
537  AliMUONTrackParam* previousTrackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->Before(trackParam));
538  AliMUONTrackParam* nextTrackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->After(trackParam));
539  if (!previousTrackParam && !nextTrackParam) {
540  AliError(Form("unable to find a cluster neither before nor after the one at position %d !?!", iCluster));
541  track.RecursiveDump();
542  break;
543  }
544 
545  // remove mono-cathod clusters on stations 3-4-5 if required
546  if (fRemoveMonoCathCl && trackParam->GetClusterPtr()->GetChamberId() > 3) {
547  Bool_t hasBending, hasNonBending;
548  if (fCheckAllPads) CheckPads(trackParam->GetClusterPtr(), hasBending, hasNonBending);
549  else CheckPadsBelow(trackParam->GetClusterPtr(), hasBending, hasNonBending);
550  if (!hasBending || !hasNonBending) continue;
551  }
552 
553  // save current trackParam and remove it from the track
554  AliMUONTrackParam currentTrackParam(*trackParam);
555  track.RemoveTrackParamAtCluster(trackParam);
556 
557  // get cluster info
558  AliMUONVCluster* cluster = currentTrackParam.GetClusterPtr();
559  Int_t chId = cluster->GetChamberId();
560  Int_t halfChId = (cluster->GetX() > 0) ? 2*chId : 2*chId+1;
561  Int_t deId = cluster->GetDetElemId();
562 
563  // compute residuals with cluster still attached to the track
564  AliMUONTrackParam* referenceTrackParam = static_cast<AliMUONTrackParam*>(referenceTrack.GetTrackParamAtCluster()->UncheckedAt(iCluster));
565  Double_t deltaX = cluster->GetX() - referenceTrackParam->GetNonBendingCoor();
566  Double_t deltaY = cluster->GetY() - referenceTrackParam->GetBendingCoor();
567 
568  // compute local track angles
569  Double_t angleX = TMath::ATan(referenceTrackParam->GetNonBendingSlope())*TMath::RadToDeg();
570  Double_t angleY = TMath::ATan(referenceTrackParam->GetBendingSlope())*TMath::RadToDeg();
571 
572  // compute local chi2
573  Double_t sigmaDeltaX2 = cluster->GetErrX2() - referenceTrackParam->GetCovariances()(0,0);
574  Double_t sigmaDeltaY2 = cluster->GetErrY2() - referenceTrackParam->GetCovariances()(2,2);
575  Double_t localChi2X = (sigmaDeltaX2 > 0.) ? deltaX*deltaX/sigmaDeltaX2 : 0.;
576  Double_t localChi2Y = (sigmaDeltaY2 > 0.) ? deltaY*deltaY/sigmaDeltaY2 : 0.;
577  Double_t localChi2 = 0.5 * referenceTrackParam->GetLocalChi2();
578 
579  // fill local chi2 info at every clusters
580  ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerCh))->Fill(chId+1, localChi2X);
581  ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerCh+1))->Fill(chId+1, localChi2Y);
582  ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerCh+2))->Fill(chId+1, localChi2);
583  ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerDE))->Fill(fDEIndices[deId], localChi2X);
584  ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerDE+1))->Fill(fDEIndices[deId], localChi2Y);
585  ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerDE+2))->Fill(fDEIndices[deId], localChi2);
586 
587  // make sure the track has another cluster in the same station and can still be refitted
588  Bool_t refit = track.IsValid( 1 << (chId/2) );
589  if (refit) {
590 
591  // refit the track and proceed if everything goes fine
592  if (tracker->RefitTrack(track, kFALSE)) {
593 
594  // fill histograms of residuals with cluster still attached to the track
595  ((TH2F*)fResiduals->UncheckedAt(kResidualPerChClusterIn))->Fill(chId+1, deltaX);
596  ((TH2F*)fResiduals->UncheckedAt(kResidualPerChClusterIn+1))->Fill(chId+1, deltaY);
597  ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfChClusterIn))->Fill(halfChId+1, deltaX);
598  ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfChClusterIn+1))->Fill(halfChId+1, deltaY);
599  ((TH2F*)fResiduals->UncheckedAt(kResidualPerDEClusterIn))->Fill(fDEIndices[deId], deltaX);
600  ((TH2F*)fResiduals->UncheckedAt(kResidualPerDEClusterIn+1))->Fill(fDEIndices[deId], deltaY);
601  ((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsPClusterIn+chId))->Fill(pUncorr, deltaX);
602  ((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsPClusterIn+10+chId))->Fill(pUncorr, deltaY);
603  ((TH2F*)fResidualsVsP->UncheckedAt(kResidualVsPClusterIn))->Fill(pUncorr, deltaX);
604  ((TH2F*)fResidualsVsP->UncheckedAt(kResidualVsPClusterIn+1))->Fill(pUncorr, deltaY);
605  ((TH2F*)fResidualsVsCent->UncheckedAt(kResidualInChVsCentClusterIn+chId))->Fill(centrality, deltaX);
606  ((TH2F*)fResidualsVsCent->UncheckedAt(kResidualInChVsCentClusterIn+10+chId))->Fill(centrality, deltaY);
607  ((TH2F*)fResidualsVsCent->UncheckedAt(kResidualVsCentClusterIn))->Fill(centrality, deltaX);
608  ((TH2F*)fResidualsVsCent->UncheckedAt(kResidualVsCentClusterIn+1))->Fill(centrality, deltaY);
609  ((TH2F*)fResidualsVsAngle->UncheckedAt(kResidualInChVsAngleClusterIn+chId))->Fill(angleX, deltaX);
610  ((TH2F*)fResidualsVsAngle->UncheckedAt(kResidualInChVsAngleClusterIn+10+chId))->Fill(angleY, deltaY);
611  ((TH2F*)fResidualsVsAngle->UncheckedAt(kResidualVsAngleClusterIn))->Fill(angleX, deltaX);
612  ((TH2F*)fResidualsVsAngle->UncheckedAt(kResidualVsAngleClusterIn+1))->Fill(angleY, deltaY);
613  ((TH2F*)fResidualsVsAngle->UncheckedAt(kResidualInHalfChVsAngleClusterIn+halfChId))->Fill(angleX, deltaX);
614  ((TH2F*)fResidualsVsAngle->UncheckedAt(kResidualInHalfChVsAngleClusterIn+20+halfChId))->Fill(angleY, deltaY);
615 
616  // find the track parameters closest to the current cluster position
617  Double_t dZWithPrevious = (previousTrackParam) ? TMath::Abs(previousTrackParam->GetClusterPtr()->GetZ() - cluster->GetZ()) : FLT_MAX;
618  Int_t previousChId = (previousTrackParam) ? previousTrackParam->GetClusterPtr()->GetChamberId() : -1;
619  Double_t dZWithNext = (nextTrackParam) ? TMath::Abs(nextTrackParam->GetClusterPtr()->GetZ() - cluster->GetZ()) : FLT_MAX;
620  AliMUONTrackParam* startingTrackParam = (nextTrackParam) ? nextTrackParam : previousTrackParam;
621  if ((fExtrapMode == 0 && previousTrackParam && dZWithPrevious < dZWithNext) ||
622  (fExtrapMode == 1 && previousTrackParam && !(chId/2 == 2 && previousChId/2 == 1) &&
623  !(chId/2 == 3 && previousChId/2 == 2))) startingTrackParam = previousTrackParam;
624 
625  // reset current parameters
626  currentTrackParam.SetParameters(startingTrackParam->GetParameters());
627  currentTrackParam.SetZ(startingTrackParam->GetZ());
628  currentTrackParam.SetCovariances(startingTrackParam->GetCovariances());
629  currentTrackParam.ResetPropagator();
630 
631  // extrapolate to the current cluster position and fill histograms of residuals if everything goes fine
632  if (AliMUONTrackExtrap::ExtrapToZCov(&currentTrackParam, currentTrackParam.GetClusterPtr()->GetZ(), kTRUE)) {
633 
634  // compute MCS dispersion on the first chamber
635  TMatrixD mcsCov(5,5);
636  if (startingTrackParam == nextTrackParam && chId == 0) {
637  AliMUONTrackParam trackParamForMCS;
638  trackParamForMCS.SetParameters(nextTrackParam->GetParameters());
639  AliMUONTrackExtrap::AddMCSEffect(&trackParamForMCS,AliMUONConstants::ChamberThicknessInX0(nextTrackParam->GetClusterPtr()->GetChamberId()),-1.);
640  const TMatrixD &propagator = currentTrackParam.GetPropagator();
641  TMatrixD tmp(trackParamForMCS.GetCovariances(),TMatrixD::kMultTranspose,propagator);
642  mcsCov.Mult(propagator,tmp);
643  } else mcsCov.Zero();
644 
645  // compute residuals
646  Double_t trackResX2 = currentTrackParam.GetCovariances()(0,0) + mcsCov(0,0);
647  Double_t trackResY2 = currentTrackParam.GetCovariances()(2,2) + mcsCov(2,2);
648  deltaX = cluster->GetX() - currentTrackParam.GetNonBendingCoor();
649  deltaY = cluster->GetY() - currentTrackParam.GetBendingCoor();
650 
651  // fill histograms with cluster not attached to the track
652  ((TH2F*)fResiduals->UncheckedAt(kResidualPerChClusterOut))->Fill(chId+1, deltaX);
653  ((TH2F*)fResiduals->UncheckedAt(kResidualPerChClusterOut+1))->Fill(chId+1, deltaY);
654  ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfChClusterOut))->Fill(halfChId+1, deltaX);
655  ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfChClusterOut+1))->Fill(halfChId+1, deltaY);
656  ((TH2F*)fResiduals->UncheckedAt(kResidualPerDEClusterOut))->Fill(fDEIndices[deId], deltaX);
657  ((TH2F*)fResiduals->UncheckedAt(kResidualPerDEClusterOut+1))->Fill(fDEIndices[deId], deltaY);
658  ((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsPClusterOut+chId))->Fill(pUncorr, deltaX);
659  ((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsPClusterOut+10+chId))->Fill(pUncorr, deltaY);
660  ((TH2F*)fResidualsVsP->UncheckedAt(kResidualVsPClusterOut))->Fill(pUncorr, deltaX);
661  ((TH2F*)fResidualsVsP->UncheckedAt(kResidualVsPClusterOut+1))->Fill(pUncorr, deltaY);
662  ((TH2F*)fResidualsVsCent->UncheckedAt(kResidualInChVsCentClusterOut+chId))->Fill(centrality, deltaX);
663  ((TH2F*)fResidualsVsCent->UncheckedAt(kResidualInChVsCentClusterOut+10+chId))->Fill(centrality, deltaY);
664  ((TH2F*)fResidualsVsCent->UncheckedAt(kResidualVsCentClusterOut))->Fill(centrality, deltaX);
665  ((TH2F*)fResidualsVsCent->UncheckedAt(kResidualVsCentClusterOut+1))->Fill(centrality, deltaY);
666  ((TH2F*)fResidualsVsAngle->UncheckedAt(kResidualInChVsAngleClusterOut+chId))->Fill(angleX, deltaX);
667  ((TH2F*)fResidualsVsAngle->UncheckedAt(kResidualInChVsAngleClusterOut+10+chId))->Fill(angleY, deltaY);
668  ((TH2F*)fResidualsVsAngle->UncheckedAt(kResidualVsAngleClusterOut))->Fill(angleX, deltaX);
669  ((TH2F*)fResidualsVsAngle->UncheckedAt(kResidualVsAngleClusterOut+1))->Fill(angleY, deltaY);
670  ((TH2F*)fResiduals->UncheckedAt(kTrackResPerCh))->Fill(chId+1, TMath::Sqrt(trackResX2));
671  ((TH2F*)fResiduals->UncheckedAt(kTrackResPerCh+1))->Fill(chId+1, TMath::Sqrt(trackResY2));
672  ((TH2F*)fResiduals->UncheckedAt(kTrackResPerHalfCh))->Fill(halfChId+1, TMath::Sqrt(trackResX2));
673  ((TH2F*)fResiduals->UncheckedAt(kTrackResPerHalfCh+1))->Fill(halfChId+1, TMath::Sqrt(trackResY2));
674  ((TH2F*)fResiduals->UncheckedAt(kTrackResPerDE))->Fill(fDEIndices[deId], TMath::Sqrt(trackResX2));
675  ((TH2F*)fResiduals->UncheckedAt(kTrackResPerDE+1))->Fill(fDEIndices[deId], TMath::Sqrt(trackResY2));
676  ((TH2F*)fResiduals->UncheckedAt(kMCSPerCh))->Fill(chId+1, TMath::Sqrt(mcsCov(0,0)));
677  ((TH2F*)fResiduals->UncheckedAt(kMCSPerCh+1))->Fill(chId+1, TMath::Sqrt(mcsCov(2,2)));
678  ((TH2F*)fResiduals->UncheckedAt(kMCSPerHalfCh))->Fill(halfChId+1, TMath::Sqrt(mcsCov(0,0)));
679  ((TH2F*)fResiduals->UncheckedAt(kMCSPerHalfCh+1))->Fill(halfChId+1, TMath::Sqrt(mcsCov(2,2)));
680  ((TH2F*)fResiduals->UncheckedAt(kMCSPerDE))->Fill(fDEIndices[deId], TMath::Sqrt(mcsCov(0,0)));
681  ((TH2F*)fResiduals->UncheckedAt(kMCSPerDE+1))->Fill(fDEIndices[deId], TMath::Sqrt(mcsCov(2,2)));
682  ((TH2F*)fResiduals->UncheckedAt(kClusterRes2PerCh))->Fill(chId+1, deltaX*deltaX - trackResX2);
683  ((TH2F*)fResiduals->UncheckedAt(kClusterRes2PerCh+1))->Fill(chId+1, deltaY*deltaY - trackResY2);
684  }
685 
686  }
687 
688  }
689 
690  // restore the track
691  track.AddTrackParamAtCluster(currentTrackParam, *(currentTrackParam.GetClusterPtr()), kTRUE);
692 
693  }
694 
695  }
696 
697  // Post final data. It will be written to a file with option "RECREATE"
698  PostData(1, fResiduals);
699  PostData(2, fResidualsVsP);
700  PostData(3, fResidualsVsCent);
701  PostData(4, fResidualsVsAngle);
702  PostData(5, fTrackRes);
703 }
704 
705 //________________________________________________________________________
707 {
709 
710  if (fOCDBLoaded) return;
711 
712  AliCDBManager* cdbm = AliCDBManager::Instance();
713  cdbm->SetDefaultStorage(fDefaultStorage.Data());
714  cdbm->SetRun(fCurrentRunNumber);
715 
716  if (!AliMUONCDB::LoadField()) return;
717 
718  if (!AliMUONCDB::LoadMapping()) return;
719 
720  AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam();
721  if (!recoParam) return;
722 
723  if (fImproveTracks) recoParam->ImproveTracks(kTRUE);
724  else recoParam->ImproveTracks(kFALSE);
725 
726  AliMUONESDInterface::ResetTracker(recoParam);
727 
728  for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
729 
730  // set the cluster resolution to default if not already set and create output objects
731  if (fClusterResNB[i] < 0.) fClusterResNB[i] = recoParam->GetDefaultNonBendingReso(i);
732  if (fClusterResB[i] < 0.) fClusterResB[i] = recoParam->GetDefaultBendingReso(i);
733 
734  // fill correspondence between DEId and indices for histo (starting from 1)
735  AliMpDEIterator it;
736  it.First(i);
737  while (!it.IsDone()) {
738  fNDE++;
739  fDEIndices[it.CurrentDEId()] = fNDE;
740  fDEIds[fNDE] = it.CurrentDEId();
741  it.Next();
742  }
743 
744  }
745 
746  // recover default storage full name (raw:// cannot be used to set specific storage)
747  TString defaultStorage(cdbm->GetDefaultStorage()->GetType());
748  if (defaultStorage == "alien") defaultStorage += Form("://folder=%s", cdbm->GetDefaultStorage()->GetBaseFolder().Data());
749  else defaultStorage += Form("://%s", cdbm->GetDefaultStorage()->GetBaseFolder().Data());
750 
751  if (fReAlign) {
752 
753  // reset existing geometry/alignment if any
754  if (cdbm->GetEntryCache()->Contains("GRP/Geometry/Data")) cdbm->UnloadFromCache("GRP/Geometry/Data");
755  if (cdbm->GetEntryCache()->Contains("MUON/Align/Data")) cdbm->UnloadFromCache("MUON/Align/Data");
756  if (AliGeomManager::GetGeometry()) AliGeomManager::GetGeometry()->UnlockGeometry();
757 
758  // get original geometry transformer
759  AliGeomManager::LoadGeometry();
760  if (!AliGeomManager::GetGeometry()) return;
761  if (fOldAlignStorage != "none") {
762  cdbm->SetSpecificStorage("MUON/Align/Data", fOldAlignStorage.IsNull() ? defaultStorage.Data() : fOldAlignStorage.Data(), fOldAlignVersion, fOldAlignSubVersion);
763  AliGeomManager::ApplyAlignObjsFromCDB("MUON");
764  }
765  fOldGeoTransformer = new AliMUONGeometryTransformer();
766  fOldGeoTransformer->LoadGeometryData();
767 
768  }
769 
770  // reset existing geometry/alignment if any
771  if (cdbm->GetEntryCache()->Contains("GRP/Geometry/Data")) cdbm->UnloadFromCache("GRP/Geometry/Data");
772  if (cdbm->GetEntryCache()->Contains("MUON/Align/Data")) cdbm->UnloadFromCache("MUON/Align/Data");
773  if (AliGeomManager::GetGeometry()) AliGeomManager::GetGeometry()->UnlockGeometry();
774 
775  // load original or realigned geometry (also used to check pads below clusters and for track extrapolation to vertex)
776  AliGeomManager::LoadGeometry();
777  if (!AliGeomManager::GetGeometry()) return;
778  if (fNewAlignStorage != "none") {
779  cdbm->SetSpecificStorage("MUON/Align/Data", fNewAlignStorage.IsNull() ? defaultStorage.Data() : fNewAlignStorage.Data(), fNewAlignVersion, fNewAlignSubVersion);
780  AliGeomManager::ApplyAlignObjsFromCDB("MUON");
781  }
782  fNewGeoTransformer = new AliMUONGeometryTransformer();
783  fNewGeoTransformer->LoadGeometryData();
784 
785  // print starting chamber resolution if required
786  if (fPrintClResPerCh) {
787  printf("\nstarting chamber resolution:\n");
788  printf(" - non-bending:");
789  for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %5.3f":", %5.3f",fClusterResNB[i]);
790  printf("\n - bending:");
791  for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %6.4f":", %6.4f",fClusterResB[i]);
792  printf("\n\n");
793  }
794 
795  if (fMuonEventCuts) fMuonEventCuts->Print();
796 
797  // get the trackCuts for this run
798  if (fMuonTrackCuts) {
799  fMuonTrackCuts->SetRun(fInputHandler);
800  fMuonTrackCuts->Print();
801  }
802 
803  fOCDBLoaded = kTRUE;
804 
806 
807 }
808 
809 //________________________________________________________________________
811 {
813 
814  // recover output objects
815  fResiduals = static_cast<TObjArray*> (GetOutputData(1));
816  fResidualsVsP = static_cast<TObjArray*> (GetOutputData(2));
817  fResidualsVsCent = static_cast<TObjArray*> (GetOutputData(3));
818  fResidualsVsAngle = static_cast<TObjArray*> (GetOutputData(4));
819  fTrackRes = static_cast<TObjArray*> (GetOutputData(5));
821 
822  // summary graphs
823  fLocalChi2 = new TObjArray(1000);
824  fLocalChi2->SetOwner();
825  fChamberRes = new TObjArray(1000);
826  fChamberRes->SetOwner();
827  fCanvases = new TObjArray(1000);
828  fCanvases->SetOwner();
829  fTmpHists = new TObjArray(1000);
830  fTmpHists->SetOwner();
831  TGraphErrors* g;
832  TMultiGraph* mg;
833  TCanvas* c;
834 
835  const char* axes[2] = {"X", "Y"};
836  const char* side[2] = {"I", "O"};
837  Double_t newClusterRes[2][10], newClusterResErr[2][10];
838  fNDE = ((TH2F*)fResiduals->UncheckedAt(kResidualPerDEClusterIn))->GetXaxis()->GetNbins();
839  TString name, title;
840 
841  for (Int_t ia = 0; ia < 2; ia++) {
842 
843  // shifts per chamber
844  g = new TGraphErrors(AliMUONConstants::NTrackingCh());
845  g->SetName(Form("gResidual%sPerChMean_ClusterIn",axes[ia]));
846  g->SetTitle(Form("cluster-track residual-%s per Ch: mean (cluster in);chamber ID;<#Delta_{%s}> (cm)",axes[ia],axes[ia]));
847  g->SetMarkerStyle(kFullDotLarge);
848  fChamberRes->AddAtAndExpand(g, kResidualPerChMeanClusterIn+ia);
849  g = new TGraphErrors(AliMUONConstants::NTrackingCh());
850  g->SetName(Form("gResidual%sPerChMean_ClusterOut",axes[ia]));
851  g->SetTitle(Form("cluster-track residual-%s per Ch: mean (cluster out);chamber ID;<#Delta_{%s}> (cm)",axes[ia],axes[ia]));
852  g->SetMarkerStyle(kFullDotLarge);
853  fChamberRes->AddAtAndExpand(g, kResidualPerChMeanClusterOut+ia);
854 
855  // shifts per half-chamber
856  g = new TGraphErrors(2*AliMUONConstants::NTrackingCh());
857  g->SetName(Form("gResidual%sPerHalfChMean_ClusterIn",axes[ia]));
858  g->SetTitle(Form("cluster-track residual-%s per half Ch: mean (cluster in);half chamber ID;<#Delta_{%s}> (cm)",axes[ia],axes[ia]));
859  g->SetMarkerStyle(kFullDotLarge);
860  fChamberRes->AddAtAndExpand(g, kResidualPerHalfChMeanClusterIn+ia);
861  g = new TGraphErrors(2*AliMUONConstants::NTrackingCh());
862  g->SetName(Form("gResidual%sPerHalfChMean_ClusterOut",axes[ia]));
863  g->SetTitle(Form("cluster-track residual-%s per half Ch: mean (cluster out);half chamber ID;<#Delta_{%s}> (cm)",axes[ia],axes[ia]));
864  g->SetMarkerStyle(kFullDotLarge);
865  fChamberRes->AddAtAndExpand(g, kResidualPerHalfChMeanClusterOut+ia);
866 
867  // shifts per DE
868  g = new TGraphErrors(fNDE);
869  g->SetName(Form("gResidual%sPerDEMean_ClusterIn",axes[ia]));
870  g->SetTitle(Form("cluster-track residual-%s per DE: mean (cluster in);DE ID;<#Delta_{%s}> (cm)",axes[ia],axes[ia]));
871  g->SetMarkerStyle(kFullDotLarge);
872  fChamberRes->AddAtAndExpand(g, kResidualPerDEMeanClusterIn+ia);
873  g = new TGraphErrors(fNDE);
874  g->SetName(Form("gResidual%sPerDEMean_ClusterOut",axes[ia]));
875  g->SetTitle(Form("cluster-track residual-%s per DE: mean (cluster out);DE ID;<#Delta_{%s}> (cm)",axes[ia],axes[ia]));
876  g->SetMarkerStyle(kFullDotLarge);
877  fChamberRes->AddAtAndExpand(g, kResidualPerDEMeanClusterOut+ia);
878 
879  // resolution per chamber
880  g = new TGraphErrors(AliMUONConstants::NTrackingCh());
881  g->SetName(Form("gResidual%sPerChSigma_ClusterIn",axes[ia]));
882  g->SetTitle(Form("cluster-track residual-%s per Ch: sigma (cluster in);chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
883  g->SetMarkerStyle(kFullDotLarge);
884  fChamberRes->AddAtAndExpand(g, kResidualPerChSigmaClusterIn+ia);
885  g = new TGraphErrors(AliMUONConstants::NTrackingCh());
886  g->SetName(Form("gResidual%sPerChSigma_ClusterOut",axes[ia]));
887  g->SetTitle(Form("cluster-track residual-%s per Ch: sigma (cluster out);chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
888  g->SetMarkerStyle(kFullDotLarge);
889  fChamberRes->AddAtAndExpand(g, kResidualPerChSigmaClusterOut+ia);
890 
891  g = new TGraphErrors(AliMUONConstants::NTrackingCh());
892  g->SetName(Form("gResidual%sPerChDispersion_ClusterOut",axes[ia]));
893  g->SetTitle(Form("cluster-track residual-%s per Ch: dispersion (cluster out);chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
894  g->SetMarkerStyle(kFullDotLarge);
895  fChamberRes->AddAtAndExpand(g, kResidualPerChDispersionClusterOut+ia);
896 
897  g = new TGraphErrors(AliMUONConstants::NTrackingCh());
898  g->SetName(Form("gCombinedResidual%sPerChSigma",axes[ia]));
899  g->SetTitle(Form("combined cluster-track residual-%s per Ch: sigma;chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
900  g->SetMarkerStyle(kFullDotLarge);
901  fChamberRes->AddAtAndExpand(g, kCombinedResidualPerChSigma+ia);
902 
903  // resolution per half-chamber
904  g = new TGraphErrors(2*AliMUONConstants::NTrackingCh());
905  g->SetName(Form("gCombinedResidual%sPerHalfChSigma",axes[ia]));
906  g->SetTitle(Form("combined cluster-track residual-%s per half Ch: sigma;half chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
907  g->SetMarkerStyle(kFullDotLarge);
908  fChamberRes->AddAtAndExpand(g, kCombinedResidualPerHalfChSigma+ia);
909 
910  // resolution per DE
911  g = new TGraphErrors(fNDE);
912  g->SetName(Form("gCombinedResidual%sPerDESigma",axes[ia]));
913  g->SetTitle(Form("combined cluster-track residual-%s per DE: sigma;DE ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
914  g->SetMarkerStyle(kFullDotLarge);
915  fChamberRes->AddAtAndExpand(g, kCombinedResidualPerDESigma+ia);
916 
917  // shifts versus p
918  name = Form("mgResidual%sMeanVsP_ClusterIn",axes[ia]); title = Form("cluster-track residual-%s per chamber versus momentum: mean (cluster in);p (GeV/c^{2});<#Delta_{%s}> (cm)",axes[ia],axes[ia]);
919  mg = new TMultiGraph(name.Data(), title.Data());
920  for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
921  g = new TGraphErrors(((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsPClusterIn+10*ia+i))->GetNbinsX());
922  g->SetName(Form("gShift%sVsP_ch%d",axes[ia],i+1));
923  g->SetMarkerStyle(kFullDotMedium);
924  g->SetMarkerColor(i+1+i/9);
925  mg->Add(g,"p");
926  }
927  fChamberRes->AddAtAndExpand(mg, kResidualMeanClusterInVsP+ia);
928 
929  // resolutions versus p
930  name = Form("mgCombinedResidual%sSigmaVsP",axes[ia]); title = Form("cluster %s-resolution per chamber versus momentum;p (GeV/c^{2});#sigma_{%s} (cm)",axes[ia],axes[ia]);
931  mg = new TMultiGraph(name.Data(), title.Data());
932  for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
933  g = new TGraphErrors(((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsPClusterIn+10*ia+i))->GetNbinsX());
934  g->SetName(Form("gRes%sVsP_ch%d",axes[ia],i+1));
935  g->SetMarkerStyle(kFullDotMedium);
936  g->SetMarkerColor(i+1+i/9);
937  mg->Add(g,"p");
938  }
939  fChamberRes->AddAtAndExpand(mg, kCombinedResidualSigmaVsP+ia);
940 
941  g = new TGraphErrors(((TH2F*)fResidualsVsP->UncheckedAt(kResidualVsPClusterIn+ia))->GetNbinsX());
942  g->SetName(Form("gCombinedResidual%sSigmaVsP",axes[ia]));
943  g->SetTitle(Form("cluster %s-resolution integrated over chambers versus momentum: sigma;p (GeV/c^{2});#sigma_{%s} (cm)",axes[ia],axes[ia]));
944  g->SetMarkerStyle(kFullDotLarge);
945  fChamberRes->AddAtAndExpand(g, kCombinedResidualAllChSigmaVsP+ia);
946 
947  // shifts versus centrality
948  name = Form("mgResidual%sMeanVsCent_ClusterIn",axes[ia]); title = Form("cluster-track residual-%s per chamber versus centrality: mean (cluster in);centrality (%%);<#Delta_{%s}> (cm)",axes[ia],axes[ia]);
949  mg = new TMultiGraph(name.Data(), title.Data());
950  for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
951  g = new TGraphErrors(((TH2F*)fResidualsVsCent->UncheckedAt(kResidualInChVsCentClusterIn+10*ia+i))->GetNbinsX());
952  g->SetName(Form("gShift%sVsCent_ch%d",axes[ia],i+1));
953  g->SetMarkerStyle(kFullDotMedium);
954  g->SetMarkerColor(i+1+i/9);
955  mg->Add(g,"p");
956  }
957  fChamberRes->AddAtAndExpand(mg, kResidualMeanClusterInVsCent+ia);
958 
959  // resolutions versus centrality
960  name = Form("mgCombinedResidual%sSigmaVsCent",axes[ia]); title = Form("cluster %s-resolution per chamber versus centrality;centrality (%%);#sigma_{%s} (cm)",axes[ia],axes[ia]);
961  mg = new TMultiGraph(name.Data(), title.Data());
962  for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
963  g = new TGraphErrors(((TH2F*)fResidualsVsCent->UncheckedAt(kResidualInChVsCentClusterIn+10*ia+i))->GetNbinsX());
964  g->SetName(Form("gRes%sVsCent_ch%d",axes[ia],i+1));
965  g->SetMarkerStyle(kFullDotMedium);
966  g->SetMarkerColor(i+1+i/9);
967  mg->Add(g,"p");
968  }
969  fChamberRes->AddAtAndExpand(mg, kCombinedResidualSigmaVsCent+ia);
970 
971  g = new TGraphErrors(((TH2F*)fResidualsVsCent->UncheckedAt(kResidualVsCentClusterIn+ia))->GetNbinsX());
972  g->SetName(Form("gCombinedResidual%sSigmaVsCent",axes[ia]));
973  g->SetTitle(Form("cluster %s-resolution integrated over chambers versus centrality: sigma;centrality (%%);#sigma_{%s} (cm)",axes[ia],axes[ia]));
974  g->SetMarkerStyle(kFullDotLarge);
975  fChamberRes->AddAtAndExpand(g, kCombinedResidualAllChSigmaVsCent+ia);
976 
977  // shifts versus track angle
978  name = Form("mgResidual%sMeanVsAngle_ClusterIn",axes[ia]); title = Form("cluster-track residual-%s per chamber versus track angle: mean (cluster in);%s-angle (deg);<#Delta_{%s}> (cm)",axes[ia],axes[ia],axes[ia]);
979  mg = new TMultiGraph(name.Data(), title.Data());
980  for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
981  g = new TGraphErrors(((TH2F*)fResidualsVsAngle->UncheckedAt(kResidualInChVsAngleClusterIn+10*ia+i))->GetNbinsX());
982  g->SetName(Form("gShift%sVsAngle_ch%d",axes[ia],i+1));
983  g->SetMarkerStyle(kFullDotMedium);
984  g->SetMarkerColor(i+1+i/9);
985  mg->Add(g,"p");
986  }
987  fChamberRes->AddAtAndExpand(mg, kResidualMeanClusterInVsAngle+ia);
988 
989  name = Form("mgHChResidual%sMeanVsAngle_ClusterIn",axes[ia]); title = Form("cluster-track residual-%s per half-chamber versus track angle: mean (cluster in);%s-angle (deg);<#Delta_{%s}> (cm)",axes[ia],axes[ia],axes[ia]);
990  mg = new TMultiGraph(name.Data(), title.Data());
991  for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
992  for (Int_t j = 0; j < 2; j++) {
993  g = new TGraphErrors(((TH2F*)fResidualsVsAngle->UncheckedAt(kResidualInHalfChVsAngleClusterIn+20*ia+2*i+j))->GetNbinsX());
994  g->SetName(Form("gShift%sVsAngle_halfCh%d%s",axes[ia],i+1,side[j]));
995  g->SetMarkerStyle(kFullDotMedium);
996  g->SetMarkerColor(2*i+j+1+(2*i+j)/9);
997  mg->Add(g,"p");
998  }
999  }
1000  fChamberRes->AddAtAndExpand(mg, kHChResidualMeanClusterInVsAngle+ia);
1001 
1002  // resolutions versus track angle
1003  name = Form("mgCombinedResidual%sSigmaVsAngle",axes[ia]); title = Form("cluster %s-resolution per chamber versus track angle;%s-angle (deg);#sigma_{%s} (cm)",axes[ia],axes[ia],axes[ia]);
1004  mg = new TMultiGraph(name.Data(), title.Data());
1005  for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
1006  g = new TGraphErrors(((TH2F*)fResidualsVsAngle->UncheckedAt(kResidualInChVsAngleClusterIn+10*ia+i))->GetNbinsX());
1007  g->SetName(Form("gRes%sVsAngle_ch%d",axes[ia],i+1));
1008  g->SetMarkerStyle(kFullDotMedium);
1009  g->SetMarkerColor(i+1+i/9);
1010  mg->Add(g,"p");
1011  }
1012  fChamberRes->AddAtAndExpand(mg, kCombinedResidualSigmaVsAngle+ia);
1013 
1014  g = new TGraphErrors(((TH2F*)fResidualsVsAngle->UncheckedAt(kResidualVsAngleClusterIn+ia))->GetNbinsX());
1015  g->SetName(Form("gCombinedResidual%sSigmaVsAngle",axes[ia]));
1016  g->SetTitle(Form("cluster %s-resolution integrated over chambers versus track angle: sigma;%s-angle (deg);#sigma_{%s} (cm)",axes[ia],axes[ia],axes[ia]));
1017  g->SetMarkerStyle(kFullDotLarge);
1018  fChamberRes->AddAtAndExpand(g, kCombinedResidualAllChSigmaVsAngle+ia);
1019 
1020  g = new TGraphErrors(AliMUONConstants::NTrackingCh());
1021  g->SetName(Form("gTrackRes%sPerCh",axes[ia]));
1022  g->SetTitle(Form("track <#sigma_{%s}> per Ch;chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
1023  g->SetMarkerStyle(kFullDotLarge);
1024  fChamberRes->AddAtAndExpand(g, kTrackResPerChMean+ia);
1025 
1026  g = new TGraphErrors(AliMUONConstants::NTrackingCh());
1027  g->SetName(Form("gMCS%sPerCh",axes[ia]));
1028  g->SetTitle(Form("MCS %s-dispersion of extrapolated track per Ch;chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
1029  g->SetMarkerStyle(kFullDotLarge);
1030  fChamberRes->AddAtAndExpand(g, kMCSPerChMean+ia);
1031 
1032  g = new TGraphErrors(AliMUONConstants::NTrackingCh());
1033  g->SetName(Form("gClusterRes%sPerCh",axes[ia]));
1034  g->SetTitle(Form("cluster <#sigma_{%s}> per Ch;chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
1035  g->SetMarkerStyle(kFullDotLarge);
1036  fChamberRes->AddAtAndExpand(g, kClusterResPerCh+ia);
1037 
1038  g = new TGraphErrors(2*AliMUONConstants::NTrackingCh());
1039  g->SetName(Form("gClusterRes%sPerHalfCh",axes[ia]));
1040  g->SetTitle(Form("cluster <#sigma_{%s}> per half Ch;half chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
1041  g->SetMarkerStyle(kFullDotLarge);
1042  fChamberRes->AddAtAndExpand(g, kClusterResPerHalfCh+ia);
1043 
1044  g = new TGraphErrors(fNDE);
1045  g->SetName(Form("gClusterRes%sPerDE",axes[ia]));
1046  g->SetTitle(Form("cluster <#sigma_{%s}> per DE;DE ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
1047  g->SetMarkerStyle(kFullDotLarge);
1048  fChamberRes->AddAtAndExpand(g, kClusterResPerDE+ia);
1049 
1050  g = new TGraphErrors(AliMUONConstants::NTrackingCh());
1051  g->SetName(Form("gCalcClusterRes%sPerCh",axes[ia]));
1052  g->SetTitle(Form("calculated cluster <#sigma_{%s}> per Ch;chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
1053  g->SetMarkerStyle(kFullDotLarge);
1054  fChamberRes->AddAtAndExpand(g, kCalcClusterResPerCh+ia);
1055 
1056  g = new TGraphErrors(AliMUONConstants::NTrackingCh());
1057  g->SetName(Form("gLocalChi2%sPerChMean",axes[ia]));
1058  g->SetTitle(Form("local chi2-%s per Ch: mean;chamber ID;<local #chi^{2}_{%s}>",axes[ia],axes[ia]));
1059  g->SetMarkerStyle(kFullDotLarge);
1060  fLocalChi2->AddAtAndExpand(g, kLocalChi2PerChMean+ia);
1061 
1062  g = new TGraphErrors(fNDE);
1063  g->SetName(Form("gLocalChi2%sPerDEMean",axes[ia]));
1064  g->SetTitle(Form("local chi2-%s per DE: mean;DE ID;<local #chi^{2}_{%s}>",axes[ia],axes[ia]));
1065  g->SetMarkerStyle(kFullDotLarge);
1066  fLocalChi2->AddAtAndExpand(g, kLocalChi2PerDEMean+ia);
1067 
1068  // canvases
1069  name = Form("cDetailRes%sPerChClIn",axes[ia]); title = Form("cDetailRes%sPerChClIn",axes[ia]);
1070  c = new TCanvas(name.Data(), title.Data(), 1200, 500);
1071  c->Divide(5,2);
1072  fCanvases->AddAtAndExpand(c, kDetailResPerCh+2*ia);
1073 
1074  name = Form("cDetailRes%sPerChClOut",axes[ia]); title = Form("cDetailRes%sPerChClOut",axes[ia]);
1075  c = new TCanvas(name.Data(), title.Data(), 1200, 500);
1076  c->Divide(5,2);
1077  fCanvases->AddAtAndExpand(c, kDetailResPerCh+1+2*ia);
1078 
1079  name = Form("cDetailRes%sPerHalfChClIn",axes[ia]); title = Form("cDetailRes%sPerHalfChClIn",axes[ia]);
1080  c = new TCanvas(name.Data(), title.Data(), 1200, 800);
1081  c->Divide(5,4);
1082  fCanvases->AddAtAndExpand(c, kDetailResPerHalfCh+2*ia);
1083 
1084  name = Form("cDetailRes%sPerHalfChClOut",axes[ia]); title = Form("cDetailRes%sPerHalfChClOut",axes[ia]);
1085  c = new TCanvas(name.Data(), title.Data(), 1200, 800);
1086  c->Divide(5,4);
1087  fCanvases->AddAtAndExpand(c, kDetailResPerHalfCh+1+2*ia);
1088 
1089  name = Form("cDetailRes%sPerDEClIn",axes[ia]); title = Form("cDetailRes%sPerDEClIn",axes[ia]);
1090  c = new TCanvas(name.Data(), title.Data(), 1200, 800);
1091  c->Divide(13,12);
1092  fCanvases->AddAtAndExpand(c, kDetailResPerDE+2*ia);
1093 
1094  name = Form("cDetailRes%sPerDEClOut",axes[ia]); title = Form("cDetailRes%sPerDEClOut",axes[ia]);
1095  c = new TCanvas(name.Data(), title.Data(), 1200, 800);
1096  c->Divide(13,12);
1097  fCanvases->AddAtAndExpand(c, kDetailResPerDE+1+2*ia);
1098 
1099  // compute residual mean and dispersion integrated over chambers versus p
1101  (TH2F*)fResidualsVsP->UncheckedAt(kResidualVsPClusterOut+ia),
1102  0x0, (TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualAllChSigmaVsP+ia));
1103 
1104  // compute residual mean and dispersion integrated over chambers versus centrality
1106  (TH2F*)fResidualsVsCent->UncheckedAt(kResidualVsCentClusterOut+ia),
1107  0x0, (TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualAllChSigmaVsCent+ia));
1108 
1109  // compute residual mean and dispersion integrated over chambers versus track angle
1111  (TH2F*)fResidualsVsAngle->UncheckedAt(kResidualVsAngleClusterOut+ia),
1112  0x0, (TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualAllChSigmaVsAngle+ia));
1113 
1114  // compute residual mean and dispersion and averaged local chi2 per chamber and half chamber
1115  Double_t meanIn, meanInErr, meanOut, meanOutErr, sigma, sigmaIn, sigmaInErr, sigmaOut, sigmaOutErr;
1116  Double_t sigmaTrack, sigmaTrackErr, sigmaMCS, sigmaMCSErr, clusterRes, clusterResErr, sigmaCluster, sigmaClusterErr;
1117  Double_t dumy1, dumy2;
1118  for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
1119 
1120  // method 1
1121  TH1D *tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerChClusterIn+ia))->ProjectionY(Form("hRes%sCh%dClIn",axes[ia],i+1),i+1,i+1,"e");
1122  tmp->SetTitle(Form("chamber %d",i+1));
1123  GetMeanRMS(tmp, meanIn, meanInErr, sigmaIn, sigmaInErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChMeanClusterIn+ia),
1124  (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChSigmaClusterIn+ia), i, i+1);
1125  fTmpHists->AddLast(tmp);
1126 
1127  tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerChClusterOut+ia))->ProjectionY(Form("hRes%sCh%dClOut",axes[ia],i+1),i+1,i+1,"e");
1128  tmp->SetTitle(Form("chamber %d",i+1));
1129  GetMeanRMS(tmp, meanOut, meanOutErr, sigmaOut, sigmaOutErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChMeanClusterOut+ia),
1130  (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChSigmaClusterOut+ia), i, i+1);
1131  fTmpHists->AddLast(tmp);
1132 
1133  if (fCorrectForSystematics) {
1134  sigma = TMath::Sqrt(sigmaIn*sigmaIn + meanIn*meanIn);
1135  sigmaInErr = (sigma>0) ? TMath::Sqrt(sigmaIn*sigmaIn*sigmaInErr*sigmaInErr + meanIn*meanIn*meanInErr*meanInErr) / sigma : 0.;
1136  sigmaIn = sigma;
1137  sigma = TMath::Sqrt(sigmaOut*sigmaOut + meanOut*meanOut);
1138  sigmaOutErr = (sigma>0) ? TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + meanOut*meanOut*meanOutErr*meanOutErr) / sigma : 0.;
1139  sigmaOut = sigma;
1140  }
1141  ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChDispersionClusterOut+ia))->SetPoint(i, i+1, sigmaOut);
1142  ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChDispersionClusterOut+ia))->SetPointError(i, 0., sigmaOutErr);
1143 
1144  clusterRes = TMath::Sqrt(sigmaIn*sigmaOut);
1145  // clusterResErr = (clusterRes > 0.) ? 0.5 * TMath::Sqrt(sigmaInErr*sigmaInErr*sigmaOut*sigmaOut + sigmaIn*sigmaIn*sigmaOutErr*sigmaOutErr) / clusterRes : 0.;
1146  clusterResErr = TMath::Sqrt(sigmaInErr*sigmaOutErr);
1147  ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerChSigma+ia))->SetPoint(i, i+1, clusterRes);
1148  ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerChSigma+ia))->SetPointError(i, 0., clusterResErr);
1149  newClusterRes[ia][i] = clusterRes;
1150  newClusterResErr[ia][i] = clusterResErr;
1151 
1152  // method 2
1153  tmp = ((TH2F*)fResiduals->UncheckedAt(kTrackResPerCh+ia))->ProjectionY("tmp",i+1,i+1,"e");
1154  GetMeanRMS(tmp, sigmaTrack, sigmaTrackErr, dumy1, dumy2, (TGraphErrors*)fChamberRes->UncheckedAt(kTrackResPerChMean+ia), 0x0, i, i+1, kFALSE, kFALSE);
1155  delete tmp;
1156 
1157  tmp = ((TH2F*)fResiduals->UncheckedAt(kMCSPerCh+ia))->ProjectionY("tmp",i+1,i+1,"e");
1158  GetMeanRMS(tmp, sigmaMCS, sigmaMCSErr, dumy1, dumy2, (TGraphErrors*)fChamberRes->UncheckedAt(kMCSPerChMean+ia), 0x0, i, i+1, kFALSE, kFALSE);
1159  delete tmp;
1160 
1161  sigmaCluster = sigmaOut*sigmaOut - sigmaTrack*sigmaTrack;
1162  if (sigmaCluster > 0.) {
1163  sigmaCluster = TMath::Sqrt(sigmaCluster);
1164  sigmaClusterErr = TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + sigmaTrack*sigmaTrack*sigmaTrackErr*sigmaTrackErr) / sigmaCluster;
1165  } else {
1166  sigmaCluster = 0.;
1167  sigmaClusterErr = 0.;
1168  }
1169  ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerCh+ia))->SetPoint(i, i+1, sigmaCluster);
1170  ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerCh+ia))->SetPointError(i, 0., sigmaClusterErr);
1171 
1172  // method 3
1173  tmp = ((TH2F*)fResiduals->UncheckedAt(kClusterRes2PerCh+ia))->ProjectionY("tmp",i+1,i+1,"e");
1174  ZoomRight(tmp);
1175  clusterRes = tmp->GetMean();
1176  if (clusterRes > 0.) {
1177  ((TGraphErrors*)fChamberRes->UncheckedAt(kCalcClusterResPerCh+ia))->SetPoint(i, i+1, TMath::Sqrt(clusterRes));
1178  ((TGraphErrors*)fChamberRes->UncheckedAt(kCalcClusterResPerCh+ia))->SetPointError(i, 0., 0.5 * tmp->GetMeanError() / TMath::Sqrt(clusterRes));
1179  } else {
1180  ((TGraphErrors*)fChamberRes->UncheckedAt(kCalcClusterResPerCh+ia))->SetPoint(i, i+1, 0.);
1181  ((TGraphErrors*)fChamberRes->UncheckedAt(kCalcClusterResPerCh+ia))->SetPointError(i, 0., 0.);
1182  }
1183  delete tmp;
1184 
1185  // method 1 versus p
1187  (TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsPClusterOut+10*ia+i),
1188  (TGraphErrors*)((TMultiGraph*)fChamberRes->UncheckedAt(kResidualMeanClusterInVsP+ia))->GetListOfGraphs()->FindObject(Form("gShift%sVsP_ch%d",axes[ia],i+1)),
1189  (TGraphErrors*)((TMultiGraph*)fChamberRes->UncheckedAt(kCombinedResidualSigmaVsP+ia))->GetListOfGraphs()->FindObject(Form("gRes%sVsP_ch%d",axes[ia],i+1)));
1190 
1191  // method 1 versus centrality
1193  (TH2F*)fResidualsVsCent->UncheckedAt(kResidualInChVsCentClusterOut+10*ia+i),
1194  (TGraphErrors*)((TMultiGraph*)fChamberRes->UncheckedAt(kResidualMeanClusterInVsCent+ia))->GetListOfGraphs()->FindObject(Form("gShift%sVsCent_ch%d",axes[ia],i+1)),
1195  (TGraphErrors*)((TMultiGraph*)fChamberRes->UncheckedAt(kCombinedResidualSigmaVsCent+ia))->GetListOfGraphs()->FindObject(Form("gRes%sVsCent_ch%d",axes[ia],i+1)));
1196 
1197  // method 1 versus track angle
1199  (TH2F*)fResidualsVsAngle->UncheckedAt(kResidualInChVsAngleClusterOut+10*ia+i),
1200  (TGraphErrors*)((TMultiGraph*)fChamberRes->UncheckedAt(kResidualMeanClusterInVsAngle+ia))->GetListOfGraphs()->FindObject(Form("gShift%sVsAngle_ch%d",axes[ia],i+1)),
1201  (TGraphErrors*)((TMultiGraph*)fChamberRes->UncheckedAt(kCombinedResidualSigmaVsAngle+ia))->GetListOfGraphs()->FindObject(Form("gRes%sVsAngle_ch%d",axes[ia],i+1)));
1202 
1203  // compute residual mean and dispersion per half chamber
1204  for (Int_t j = 0; j < 2; j++) {
1205  Int_t k = 2*i+j;
1206 
1207  // method 1
1208  tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfChClusterIn+ia))->ProjectionY(Form("hRes%sHalfCh%dClIn",axes[ia],k+1),k+1,k+1,"e");
1209  tmp->SetTitle(Form("half chamber %s",((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfChClusterIn+ia))->GetXaxis()->GetBinLabel(k+1)));
1210  GetMeanRMS(tmp, meanIn, meanInErr, sigmaIn, sigmaInErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMeanClusterIn+ia), 0x0, k, k+1);
1211  fTmpHists->AddLast(tmp);
1212 
1213  tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfChClusterOut+ia))->ProjectionY(Form("hRes%sHalfCh%dClOut",axes[ia],k+1),k+1,k+1,"e");
1214  tmp->SetTitle(Form("half chamber %s",((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfChClusterOut+ia))->GetXaxis()->GetBinLabel(k+1)));
1215  GetMeanRMS(tmp, meanOut, meanOutErr, sigmaOut, sigmaOutErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMeanClusterOut+ia), 0x0, k, k+1);
1216  fTmpHists->AddLast(tmp);
1217 
1218  if (fCorrectForSystematics) {
1219  sigma = TMath::Sqrt(sigmaIn*sigmaIn + meanIn*meanIn);
1220  sigmaInErr = (sigma>0) ? TMath::Sqrt(sigmaIn*sigmaIn*sigmaInErr*sigmaInErr + meanIn*meanIn*meanInErr*meanInErr) / sigma : 0.;
1221  sigmaIn = sigma;
1222  sigma = TMath::Sqrt(sigmaOut*sigmaOut + meanOut*meanOut);
1223  sigmaOutErr = (sigma>0) ? TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + meanOut*meanOut*meanOutErr*meanOutErr) / sigma : 0.;
1224  sigmaOut = sigma;
1225  }
1226 
1227  clusterRes = TMath::Sqrt(sigmaIn*sigmaOut);
1228  // clusterResErr = (clusterRes > 0.) ? 0.5 * TMath::Sqrt(sigmaInErr*sigmaInErr*sigmaOut*sigmaOut + sigmaIn*sigmaIn*sigmaOutErr*sigmaOutErr) / clusterRes : 0.;
1229  clusterResErr = TMath::Sqrt(sigmaInErr*sigmaOutErr);
1230  ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerHalfChSigma+ia))->SetPoint(k, k+1, clusterRes);
1231  ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerHalfChSigma+ia))->SetPointError(k, 0., clusterResErr);
1232 
1233  // method 2
1234  tmp = ((TH2F*)fResiduals->UncheckedAt(kTrackResPerHalfCh+ia))->ProjectionY("tmp",k+1,k+1,"e");
1235  GetMeanRMS(tmp, sigmaTrack, sigmaTrackErr, dumy1, dumy2, 0x0, 0x0, 0, 0, kFALSE, kFALSE);
1236  delete tmp;
1237 
1238  tmp = ((TH2F*)fResiduals->UncheckedAt(kMCSPerHalfCh+ia))->ProjectionY("tmp",k+1,k+1,"e");
1239  GetMeanRMS(tmp, sigmaMCS, sigmaMCSErr, dumy1, dumy2, 0x0, 0x0, 0, 0, kFALSE, kFALSE);
1240  delete tmp;
1241 
1242  sigmaCluster = sigmaOut*sigmaOut - sigmaTrack*sigmaTrack;
1243  if (sigmaCluster > 0.) {
1244  sigmaCluster = TMath::Sqrt(sigmaCluster);
1245  sigmaClusterErr = TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + sigmaTrack*sigmaTrack*sigmaTrackErr*sigmaTrackErr) / sigmaCluster;
1246  } else {
1247  sigmaCluster = 0.;
1248  sigmaClusterErr = 0.;
1249  }
1250  ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerHalfCh+ia))->SetPoint(k, k+1, sigmaCluster);
1251  ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerHalfCh+ia))->SetPointError(k, 0., sigmaClusterErr);
1252 
1253  // method 1 versus track angle
1254  FillMeanSigmaClusterVsX((TH2F*)fResidualsVsAngle->UncheckedAt(kResidualInHalfChVsAngleClusterIn+20*ia+2*i+j), 0x0,
1255  (TGraphErrors*)((TMultiGraph*)fChamberRes->UncheckedAt(kHChResidualMeanClusterInVsAngle+ia))->GetListOfGraphs()->FindObject(Form("gShift%sVsAngle_halfCh%d%s",axes[ia],i+1,side[j])), 0x0);
1256 
1257  }
1258 
1259  // compute averaged local chi2
1260  tmp = ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerCh+ia))->ProjectionY("tmp",i+1,i+1,"e");
1261  ((TGraphErrors*)fLocalChi2->UncheckedAt(kLocalChi2PerChMean+ia))->SetPoint(i, i+1, tmp->GetMean());
1262  ((TGraphErrors*)fLocalChi2->UncheckedAt(kLocalChi2PerChMean+ia))->SetPointError(i, 0., tmp->GetMeanError());
1263  delete tmp;
1264 
1265  }
1266 
1267  // compute residual mean and dispersion per DE
1268  for (Int_t i = 0; i < fNDE; i++) {
1269 
1270  // method 1
1271  TH1D *tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerDEClusterIn+ia))->ProjectionY(Form("hRes%sDE%dClIn",axes[ia],i+1),i+1,i+1,"e");
1272  tmp->SetTitle(Form("DE %s",((TH2F*)fResiduals->UncheckedAt(kResidualPerDEClusterIn+ia))->GetXaxis()->GetBinLabel(i+1)));
1273  GetMeanRMS(tmp, meanIn, meanInErr, sigmaIn, sigmaInErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMeanClusterIn+ia), 0x0, i, i+1);
1274  fTmpHists->AddLast(tmp);
1275 
1276  tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerDEClusterOut+ia))->ProjectionY(Form("hRes%sDE%dClOut",axes[ia],i+1),i+1,i+1,"e");
1277  tmp->SetTitle(Form("DE %s",((TH2F*)fResiduals->UncheckedAt(kResidualPerDEClusterOut+ia))->GetXaxis()->GetBinLabel(i+1)));
1278  GetMeanRMS(tmp, meanOut, meanOutErr, sigmaOut, sigmaOutErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMeanClusterOut+ia), 0x0, i, i+1);
1279  fTmpHists->AddLast(tmp);
1280 
1281  if (fCorrectForSystematics) {
1282  sigma = TMath::Sqrt(sigmaIn*sigmaIn + meanIn*meanIn);
1283  sigmaInErr = (sigma>0) ? TMath::Sqrt(sigmaIn*sigmaIn*sigmaInErr*sigmaInErr + meanIn*meanIn*meanInErr*meanInErr) / sigma : 0.;
1284  sigmaIn = sigma;
1285  sigma = TMath::Sqrt(sigmaOut*sigmaOut + meanOut*meanOut);
1286  sigmaOutErr = (sigma>0) ? TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + meanOut*meanOut*meanOutErr*meanOutErr) / sigma : 0.;
1287  sigmaOut = sigma;
1288  }
1289 
1290  clusterRes = TMath::Sqrt(sigmaIn*sigmaOut);
1291  // clusterResErr = (clusterRes > 0.) ? 0.5 * TMath::Sqrt(sigmaInErr*sigmaInErr*sigmaOut*sigmaOut + sigmaIn*sigmaIn*sigmaOutErr*sigmaOutErr) / clusterRes : 0.;
1292  clusterResErr = TMath::Sqrt(sigmaInErr*sigmaOutErr);
1293  ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerDESigma+ia))->SetPoint(i, i+1, clusterRes);
1294  ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerDESigma+ia))->SetPointError(i, 0., clusterResErr);
1295 
1296  // method 2
1297  tmp = ((TH2F*)fResiduals->UncheckedAt(kTrackResPerDE+ia))->ProjectionY("tmp",i+1,i+1,"e");
1298  GetMeanRMS(tmp, sigmaTrack, sigmaTrackErr, dumy1, dumy2, 0x0, 0x0, 0, 0, kFALSE, kFALSE);
1299  delete tmp;
1300 
1301  tmp = ((TH2F*)fResiduals->UncheckedAt(kMCSPerDE+ia))->ProjectionY("tmp",i+1,i+1,"e");
1302  GetMeanRMS(tmp, sigmaMCS, sigmaMCSErr, dumy1, dumy2, 0x0, 0x0, 0, 0, kFALSE, kFALSE);
1303  delete tmp;
1304 
1305  sigmaCluster = sigmaOut*sigmaOut - sigmaTrack*sigmaTrack;
1306  if (sigmaCluster > 0.) {
1307  sigmaCluster = TMath::Sqrt(sigmaCluster);
1308  sigmaClusterErr = TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + sigmaTrack*sigmaTrack*sigmaTrackErr*sigmaTrackErr) / sigmaCluster;
1309  } else {
1310  sigmaCluster = 0.;
1311  sigmaClusterErr = 0.;
1312  }
1313  ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerDE+ia))->SetPoint(i, i+1, sigmaCluster);
1314  ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerDE+ia))->SetPointError(i, 0., sigmaClusterErr);
1315 
1316  // compute averaged local chi2
1317  tmp = ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerDE+ia))->ProjectionY("tmp",i+1,i+1,"e");
1318  ((TGraphErrors*)fLocalChi2->UncheckedAt(kLocalChi2PerDEMean+ia))->SetPoint(i, i+1, tmp->GetMean());
1319  ((TGraphErrors*)fLocalChi2->UncheckedAt(kLocalChi2PerDEMean+ia))->SetPointError(i, 0., tmp->GetMeanError());
1320  delete tmp;
1321 
1322  }
1323 
1324  // set half-chamber graph labels
1325  TAxis* xAxis = ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfChClusterIn+ia))->GetXaxis();
1326  ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMeanClusterIn+ia))->GetXaxis()->Set(20, 0.5, 20.5);
1327  ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMeanClusterOut+ia))->GetXaxis()->Set(20, 0.5, 20.5);
1328  ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerHalfChSigma+ia))->GetXaxis()->Set(20, 0.5, 20.5);
1329  ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerHalfCh+ia))->GetXaxis()->Set(20, 0.5, 20.5);
1330  for (Int_t i = 1; i <= 20; i++) {
1331  const char* label = xAxis->GetBinLabel(i);
1332  ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMeanClusterIn+ia))->GetXaxis()->SetBinLabel(i, label);
1333  ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMeanClusterOut+ia))->GetXaxis()->SetBinLabel(i, label);
1334  ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerHalfChSigma+ia))->GetXaxis()->SetBinLabel(i, label);
1335  ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerHalfCh+ia))->GetXaxis()->SetBinLabel(i, label);
1336  }
1337 
1338  // set DE graph labels
1339  xAxis = ((TH2F*)fResiduals->UncheckedAt(kResidualPerDEClusterOut+ia))->GetXaxis();
1340  ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMeanClusterIn+ia))->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
1341  ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMeanClusterOut+ia))->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
1342  ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerDESigma+ia))->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
1343  ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerDE+ia))->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
1344  ((TGraphErrors*)fLocalChi2->UncheckedAt(kLocalChi2PerDEMean+ia))->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
1345  for (Int_t i = 1; i <= fNDE; i++) {
1346  const char* label = xAxis->GetBinLabel(i);
1347  ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMeanClusterIn+ia))->GetXaxis()->SetBinLabel(i, label);
1348  ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMeanClusterOut+ia))->GetXaxis()->SetBinLabel(i, label);
1349  ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerDESigma+ia))->GetXaxis()->SetBinLabel(i, label);
1350  ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerDE+ia))->GetXaxis()->SetBinLabel(i, label);
1351  ((TGraphErrors*)fLocalChi2->UncheckedAt(kLocalChi2PerDEMean+ia))->GetXaxis()->SetBinLabel(i, label);
1352  }
1353 
1354  }
1355 
1356  // compute averaged local chi2 per chamber (X+Y)
1357  TH2F* h2 = (TH2F*)fResiduals->UncheckedAt(kLocalChi2PerCh+2);
1358  g = new TGraphErrors(AliMUONConstants::NTrackingCh());
1359  g->SetName("gLocalChi2PerChMean");
1360  g->SetTitle("local chi2 per Ch: mean;chamber ID;<local #chi^{2}>");
1361  g->SetMarkerStyle(kFullDotLarge);
1362  fLocalChi2->AddAtAndExpand(g, kLocalChi2PerChMean+2);
1363  for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
1364  TH1D* tmp = h2->ProjectionY("tmp",i+1,i+1,"e");
1365  g->SetPoint(i, i+1, tmp->GetMean());
1366  g->SetPointError(i, 0., tmp->GetMeanError());
1367  delete tmp;
1368  }
1369 
1370  // compute averaged local chi2 per DE (X+Y)
1371  h2 = (TH2F*)fResiduals->UncheckedAt(kLocalChi2PerDE+2);
1372  g = new TGraphErrors(fNDE);
1373  g->SetName("gLocalChi2PerDEMean");
1374  g->SetTitle("local chi2 per DE: mean;DE ID;<local #chi^{2}>");
1375  g->SetMarkerStyle(kFullDotLarge);
1376  fLocalChi2->AddAtAndExpand(g, kLocalChi2PerDEMean+2);
1377  for (Int_t i = 0; i < fNDE; i++) {
1378  TH1D* tmp = h2->ProjectionY("tmp",i+1,i+1,"e");
1379  g->SetPoint(i, i+1, tmp->GetMean());
1380  g->SetPointError(i, 0., tmp->GetMeanError());
1381  delete tmp;
1382  }
1383 
1384  // set graph labels
1385  g->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
1386  for (Int_t i = 1; i <= fNDE; i++) {
1387  const char* label = h2->GetXaxis()->GetBinLabel(i);
1388  g->GetXaxis()->SetBinLabel(i, label);
1389  }
1390 
1391  // display
1392  TLegend *lResPerChMean = new TLegend(0.75,0.85,0.99,0.99);
1393  TLegend *lResPerChSigma1 = new TLegend(0.75,0.85,0.99,0.99);
1394  TLegend *lResPerChSigma2 = new TLegend(0.75,0.85,0.99,0.99);
1395  TLegend *lResPerChSigma3 = new TLegend(0.75,0.85,0.99,0.99);
1396 
1397  TCanvas* cResPerCh = new TCanvas("cResPerCh","cResPerCh",1200,500);
1398  cResPerCh->Divide(4,2);
1399  for (Int_t ia = 0; ia < 2; ia++) {
1400  gROOT->SetSelectedPad(cResPerCh->cd(1+4*ia));
1401  g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChMeanClusterOut+ia);
1402  g->Draw("ap");
1403  g->SetMarkerColor(2);
1404  g->SetLineColor(2);
1405  if (ia == 0) lResPerChMean->AddEntry(g,"cluster out","PL");
1406  g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChMeanClusterIn+ia);
1407  g->Draw("p");
1408  g->SetMarkerColor(4);
1409  g->SetLineColor(4);
1410  if (ia == 0) lResPerChMean->AddEntry(g,"cluster in","PL");
1411  if (ia == 0) lResPerChMean->Draw();
1412  else lResPerChMean->DrawClone();
1413  gROOT->SetSelectedPad(cResPerCh->cd(2+4*ia));
1414  g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChSigmaClusterOut+ia);
1415  g->Draw("ap");
1416  g->SetMinimum(0.);
1417  g->SetMarkerColor(2);
1418  g->SetLineColor(2);
1419  if (ia == 0) lResPerChSigma1->AddEntry(g,"cluster out","PL");
1420  g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChSigmaClusterIn+ia);
1421  g->Draw("p");
1422  g->SetMarkerColor(4);
1423  g->SetLineColor(4);
1424  if (ia == 0) lResPerChSigma1->AddEntry(g,"cluster in","PL");
1425  g = (TGraphErrors*)fChamberRes->UncheckedAt(kMCSPerChMean+ia);
1426  g->Draw("p");
1427  g->SetMarkerColor(5);
1428  g->SetLineColor(5);
1429  if (ia == 0) lResPerChSigma1->AddEntry(g,"MCS","PL");
1430  g = (TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerChSigma+ia);
1431  g->Draw("p");
1432  g->SetMarkerColor(3);
1433  g->SetLineColor(3);
1434  if (ia == 0) lResPerChSigma1->AddEntry(g,"combined 1","PL");
1435  if (ia == 0) lResPerChSigma1->Draw();
1436  else lResPerChSigma1->DrawClone();
1437  gROOT->SetSelectedPad(cResPerCh->cd(3+4*ia));
1438  g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChDispersionClusterOut+ia);
1439  g->Draw("ap");
1440  g->SetMinimum(0.);
1441  g->SetMarkerColor(2);
1442  g->SetLineColor(2);
1443  if (ia == 0) lResPerChSigma2->AddEntry(g,"cluster out","PL");
1444  g = (TGraphErrors*)fChamberRes->UncheckedAt(kMCSPerChMean+ia);
1445  g->Draw("p");
1446  if (ia == 0) lResPerChSigma2->AddEntry(g,"MCS","PL");
1447  g = (TGraphErrors*)fChamberRes->UncheckedAt(kTrackResPerChMean+ia);
1448  g->Draw("p");
1449  g->SetMarkerColor(4);
1450  g->SetLineColor(4);
1451  if (ia == 0) lResPerChSigma2->AddEntry(g,"track res.","PL");
1452  g = (TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerCh+ia);
1453  g->Draw("p");
1454  if (ia == 0) lResPerChSigma2->AddEntry(g,"combined 2","PL");
1455  if (ia == 0) lResPerChSigma2->Draw();
1456  else lResPerChSigma2->DrawClone();
1457  gROOT->SetSelectedPad(cResPerCh->cd(4+4*ia));
1458  g = (TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerChSigma+ia);
1459  g->Draw("ap");
1460  g->SetMinimum(0.);
1461  if (ia == 0) lResPerChSigma3->AddEntry(g,"combined 1","PL");
1462  g = (TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerCh+ia);
1463  g->Draw("p");
1464  if (ia == 0) lResPerChSigma3->AddEntry(g,"combined 2","PL");
1465  if (ia == 0) lResPerChSigma3->Draw();
1466  else lResPerChSigma3->DrawClone();
1467  }
1468  fCanvases->AddAtAndExpand(cResPerCh, kResPerCh);
1469 
1470  TCanvas* cResPerHalfCh = new TCanvas("cResPerHalfCh","cResPerHalfCh",1200,500);
1471  cResPerHalfCh->Divide(2,2);
1472  for (Int_t ia = 0; ia < 2; ia++) {
1473  gROOT->SetSelectedPad(cResPerHalfCh->cd(1+2*ia));
1474  g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMeanClusterOut+ia);
1475  g->Draw("ap");
1476  g->SetMarkerColor(2);
1477  g->SetLineColor(2);
1478  g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMeanClusterIn+ia);
1479  g->Draw("p");
1480  g->SetMarkerColor(4);
1481  g->SetLineColor(4);
1482  lResPerChMean->DrawClone();
1483  gROOT->SetSelectedPad(cResPerHalfCh->cd(2+2*ia));
1484  g = (TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerHalfChSigma+ia);
1485  g->Draw("ap");
1486  g->SetMinimum(0.);
1487  g->SetMarkerColor(3);
1488  g->SetLineColor(3);
1489  g = (TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerHalfCh+ia);
1490  g->Draw("p");
1491  lResPerChSigma3->DrawClone();
1492  }
1493  fCanvases->AddAtAndExpand(cResPerHalfCh, kResPerHalfCh);
1494 
1495  TCanvas* cResPerDE = new TCanvas("cResPerDE","cResPerDE",1200,800);
1496  cResPerDE->Divide(1,4);
1497  for (Int_t ia = 0; ia < 2; ia++) {
1498  gROOT->SetSelectedPad(cResPerDE->cd(1+ia));
1499  g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMeanClusterOut+ia);
1500  g->Draw("ap");
1501  g->SetMarkerColor(2);
1502  g->SetLineColor(2);
1503  g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMeanClusterIn+ia);
1504  g->Draw("p");
1505  g->SetMarkerColor(4);
1506  g->SetLineColor(4);
1507  lResPerChMean->DrawClone();
1508  gROOT->SetSelectedPad(cResPerDE->cd(3+ia));
1509  g = (TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerDESigma+ia);
1510  g->Draw("ap");
1511  g->SetMinimum(0.);
1512  g->SetMarkerColor(3);
1513  g->SetLineColor(3);
1514  g = (TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerDE+ia);
1515  g->Draw("p");
1516  lResPerChSigma3->DrawClone();
1517  }
1518  fCanvases->AddAtAndExpand(cResPerDE, kResPerDE);
1519 
1520  TCanvas* cResPerChVsP = new TCanvas("cResPerChVsP","cResPerChVsP");
1521  cResPerChVsP->Divide(1,2);
1522  for (Int_t ia = 0; ia < 2; ia++) {
1523  gROOT->SetSelectedPad(cResPerChVsP->cd(1+ia));
1524  mg = (TMultiGraph*)fChamberRes->UncheckedAt(kCombinedResidualSigmaVsP+ia);
1525  mg->Draw("ap");
1526  }
1527  fCanvases->AddAtAndExpand(cResPerChVsP, kResPerChVsP);
1528 
1529  TCanvas* cResPerChVsCent = new TCanvas("cResPerChVsCent","cResPerChVsCent");
1530  cResPerChVsCent->Divide(1,2);
1531  for (Int_t ia = 0; ia < 2; ia++) {
1532  gROOT->SetSelectedPad(cResPerChVsCent->cd(1+ia));
1533  mg = (TMultiGraph*)fChamberRes->UncheckedAt(kCombinedResidualSigmaVsCent+ia);
1534  mg->Draw("ap");
1535  }
1536  fCanvases->AddAtAndExpand(cResPerChVsCent, kResPerChVsCent);
1537 
1538  TCanvas* cResPerChVsAngle = new TCanvas("cResPerChVsAngle","cResPerChVsAngle");
1539  cResPerChVsAngle->Divide(1,2);
1540  for (Int_t ia = 0; ia < 2; ia++) {
1541  gROOT->SetSelectedPad(cResPerChVsAngle->cd(1+ia));
1542  mg = (TMultiGraph*)fChamberRes->UncheckedAt(kCombinedResidualSigmaVsAngle+ia);
1543  mg->Draw("ap");
1544  }
1545  fCanvases->AddAtAndExpand(cResPerChVsAngle, kResPerChVsAngle);
1546 
1547  TCanvas* cShiftPerChVsP = new TCanvas("cShiftPerChVsP","cShiftPerChVsP");
1548  cShiftPerChVsP->Divide(1,2);
1549  for (Int_t ia = 0; ia < 2; ia++) {
1550  gROOT->SetSelectedPad(cShiftPerChVsP->cd(1+ia));
1551  mg = (TMultiGraph*)fChamberRes->UncheckedAt(kResidualMeanClusterInVsP+ia);
1552  mg->Draw("ap");
1553  }
1554  fCanvases->AddAtAndExpand(cShiftPerChVsP, kShiftPerChVsP);
1555 
1556  TCanvas* cShiftPerChVsCent = new TCanvas("cShiftPerChVsCent","cShiftPerChVsCent");
1557  cShiftPerChVsCent->Divide(1,2);
1558  for (Int_t ia = 0; ia < 2; ia++) {
1559  gROOT->SetSelectedPad(cShiftPerChVsCent->cd(1+ia));
1560  mg = (TMultiGraph*)fChamberRes->UncheckedAt(kResidualMeanClusterInVsCent+ia);
1561  mg->Draw("ap");
1562  }
1563  fCanvases->AddAtAndExpand(cShiftPerChVsCent, kShiftPerChVsCent);
1564 
1565  TCanvas* cShiftPerChVsAngle = new TCanvas("cShiftPerChVsAngle","cShiftPerChVsAngle");
1566  cShiftPerChVsAngle->Divide(1,2);
1567  for (Int_t ia = 0; ia < 2; ia++) {
1568  gROOT->SetSelectedPad(cShiftPerChVsAngle->cd(1+ia));
1569  mg = (TMultiGraph*)fChamberRes->UncheckedAt(kResidualMeanClusterInVsAngle+ia);
1570  mg->Draw("ap");
1571  }
1572  fCanvases->AddAtAndExpand(cShiftPerChVsAngle, kShiftPerChVsAngle);
1573 
1574  // print results
1575  if (fPrintClResPerCh) {
1576  printf("\nchamber resolution:\n");
1577  printf(" - non-bending:");
1578  for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %5.3f":", %5.3f",newClusterRes[0][i]);
1579  printf("\n - bending:");
1580  for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %6.4f":", %6.4f",newClusterRes[1][i]);
1581  printf("\n\n");
1582  }
1583 
1584  if (fPrintClResPerDE) {
1585  Double_t iDE, clRes;
1586  printf("\nDE resolution:\n");
1587  printf(" - non-bending:");
1588  for (Int_t i = 0; i < fNDE; i++) {
1589  ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerDESigma))->GetPoint(i, iDE, clRes);
1590  printf((i==0)?" %5.3f":", %5.3f", clRes);
1591  }
1592  printf("\n - bending:");
1593  for (Int_t i = 0; i < fNDE; i++) {
1594  ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerDESigma+1))->GetPoint(i, iDE, clRes);
1595  printf((i==0)?" %6.4f":", %6.4f", clRes);
1596  }
1597  printf("\n\n");
1598  }
1599 
1600  if (fPrintHalfChShift) {
1601  Double_t iHCh, hChShift;
1602  printf("\nhalf-chamber residual displacements:\n");
1603  printf(" - non-bending:");
1604  for (Int_t i = 0; i < 2*AliMUONConstants::NTrackingCh(); i++) {
1605  ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMeanClusterIn))->GetPoint(i, iHCh, hChShift);
1606  printf((i==0)?" %6.4f":", %6.4f", hChShift);
1607  }
1608  printf("\n - bending:");
1609  for (Int_t i = 0; i < 2*AliMUONConstants::NTrackingCh(); i++) {
1610  ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMeanClusterIn+1))->GetPoint(i, iHCh, hChShift);
1611  printf((i==0)?" %6.4f":", %6.4f", hChShift);
1612  }
1613  printf("\n\n");
1614  }
1615 
1616  if (fPrintDEShift) {
1617  Double_t iDE, deShift;
1618  printf("\nDE residual displacements:\n");
1619  printf(" - non-bending:");
1620  for (Int_t i = 0; i < fNDE; i++) {
1621  ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMeanClusterIn))->GetPoint(i, iDE, deShift);
1622  printf((i==0)?" %6.4f":", %6.4f", deShift);
1623  }
1624  printf("\n - bending:");
1625  for (Int_t i = 0; i < fNDE; i++) {
1626  ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMeanClusterIn+1))->GetPoint(i, iDE, deShift);
1627  printf((i==0)?" %6.4f":", %6.4f", deShift);
1628  }
1629  printf("\n\n");
1630  }
1631 
1632  // Post final data.
1633  PostData(6, fLocalChi2);
1634  PostData(7, fChamberRes);
1635 }
1636 
1637 //________________________________________________________________________
1639 {
1642 
1643  Double_t gX,gY,gZ,lX,lY,lZ;
1644 
1645  // loop over clusters
1646  Int_t nClusters = track.GetNClusters();
1647  for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
1648 
1649  AliMUONVCluster* cl = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->UncheckedAt(iCluster))->GetClusterPtr();
1650  Int_t chId = cl->GetChamberId();
1651  Int_t halfChId = (cl->GetX() > 0) ? 2*chId : 2*chId+1;
1652  Int_t deId = cl->GetDetElemId();
1653 
1654  // change their resolution
1655  cl->SetErrXY(fClusterResNB[chId], fClusterResB[chId]);
1656 
1657  // change their position
1658  gX = cl->GetX();
1659  gY = cl->GetY();
1660  gZ = cl->GetZ();
1661  if (fReAlign) { // change the alignement
1662  fOldGeoTransformer->Global2Local(deId,gX,gY,gZ,lX,lY,lZ);
1663  fNewGeoTransformer->Local2Global(deId,lX,lY,lZ,gX,gY,gZ);
1664  }
1665  if (fShiftHalfCh) { // correct for half-chamber displacement
1666  gX -= fHalfChShiftNB[halfChId];
1667  gY -= fHalfChShiftB[halfChId];
1668  }
1669  if (fShiftDE) { // correct for DE displacement
1670  gX -= fDEShiftNB[fDEIndices[deId]-1];
1671  gY -= fDEShiftB[fDEIndices[deId]-1];
1672  }
1673  cl->SetXYZ(gX,gY,gZ);
1674 
1675  // "remove" mono-cathod clusters on stations 3-4-5 if required
1676  // (to be done after moving clusters to the new position)
1677  if (fRemoveMonoCathCl && chId > 3) {
1678  Bool_t hasBending, hasNonBending;
1679  if (fCheckAllPads) CheckPads(cl, hasBending, hasNonBending);
1680  else CheckPadsBelow(cl, hasBending, hasNonBending);
1681  if (!hasNonBending) cl->SetErrXY(10., cl->GetErrY());
1682  if (!hasBending) cl->SetErrXY(cl->GetErrX(), 10.);
1683  }
1684 
1685  }
1686 
1687 }
1688 
1689 //________________________________________________________________________
1690 void AliAnalysisTaskMuonResolution::Zoom(TH1* h, Double_t fractionCut)
1691 {
1693  ZoomLeft(h, fractionCut);
1694  ZoomRight(h, fractionCut);
1695 }
1696 
1697 //________________________________________________________________________
1698 void AliAnalysisTaskMuonResolution::ZoomLeft(TH1* h, Double_t fractionCut)
1699 {
1701  Int_t maxEventsCut = (Int_t) (fractionCut * h->GetEntries());
1702  Int_t nBins = h->GetNbinsX();
1703 
1704  // set low edge
1705  Int_t minBin;
1706  Int_t eventsCut = 0;
1707  for (minBin = 1; minBin <= nBins; minBin++) {
1708  eventsCut += (Int_t) h->GetBinContent(minBin);
1709  if (eventsCut > maxEventsCut) break;
1710  }
1711 
1712  // set new axis range
1713  h->GetXaxis()->SetRange(--minBin, h->GetXaxis()->GetLast());
1714 }
1715 
1716 //________________________________________________________________________
1717 void AliAnalysisTaskMuonResolution::ZoomRight(TH1* h, Double_t fractionCut)
1718 {
1720  Int_t maxEventsCut = (Int_t) (fractionCut * h->GetEntries());
1721  Int_t nBins = h->GetNbinsX();
1722 
1723  // set high edge
1724  Int_t maxBin;
1725  Int_t eventsCut = 0;
1726  for (maxBin = nBins; maxBin >= 1; maxBin--) {
1727  eventsCut += (Int_t) h->GetBinContent(maxBin);
1728  if (eventsCut > maxEventsCut) break;
1729  }
1730 
1731  // set new axis range
1732  h->GetXaxis()->SetRange(h->GetXaxis()->GetFirst(), ++maxBin);
1733 }
1734 
1735 //________________________________________________________________________
1736 void AliAnalysisTaskMuonResolution::GetMeanRMS(TH1* h, Double_t& mean, Double_t& meanErr,
1737  Double_t& rms, Double_t& rmsErr,
1738  TGraphErrors* gMean, TGraphErrors* gRMS,
1739  Int_t i, Double_t x, Bool_t zoom, Bool_t enableFit)
1740 {
1742 
1743  if (h->GetEntries() < fgkMinEntries) { // not enough entries
1744 
1745  mean = 0.;
1746  meanErr = 0.;
1747  rms = 0.;
1748  rmsErr = 0.;
1749 
1750  } else if (enableFit && fGaus) { // take the mean of a gaussian fit
1751 
1752  // first fit
1753  Double_t xMin = h->GetXaxis()->GetXmin();
1754  Double_t xMax = h->GetXaxis()->GetXmax();
1755  fGaus->SetRange(xMin, xMax);
1756  fGaus->SetParameters(h->GetEntries(), 0., 0.1);
1757  fGaus->SetParLimits(1, xMin, xMax);
1758  h->Fit("fGaus", "WWNQ");
1759 
1760  // rebin histo
1761  Int_t rebin = static_cast<Int_t>(TMath::Min(0.1*h->GetNbinsX(),TMath::Max(0.3*fGaus->GetParameter(2)/h->GetBinWidth(1),1.)));
1762  while (h->GetNbinsX()%rebin!=0) rebin--;
1763  h->Rebin(rebin);
1764 
1765  // second fit
1766  xMin = TMath::Max(fGaus->GetParameter(1)-10.*fGaus->GetParameter(2), h->GetXaxis()->GetXmin());
1767  xMax = TMath::Min(fGaus->GetParameter(1)+10.*fGaus->GetParameter(2), h->GetXaxis()->GetXmax());
1768  fGaus->SetRange(xMin, xMax);
1769  fGaus->SetParLimits(1, xMin, xMax);
1770  h->Fit("fGaus","NQR");
1771 
1772  mean = fGaus->GetParameter(1);
1773  meanErr = fGaus->GetParError(1);
1774  rms = fGaus->GetParameter(2);
1775  rmsErr = fGaus->GetParError(2);
1776 
1777  // display the detail of the fit
1778  if (!strstr(h->GetName(),"tmp")) {
1779  Int_t ia = (strstr(h->GetName(),"ResX")) ? 0 : 1;
1780  Int_t ib = (strstr(h->GetName(),"ClIn")) ? 0 : 1;
1781  if (strstr(h->GetName(),"Half"))
1782  gROOT->SetSelectedPad(((TCanvas*)fCanvases->UncheckedAt(kDetailResPerHalfCh+ib+2*ia))->cd(i+1));
1783  else if (strstr(h->GetName(),"Ch"))
1784  gROOT->SetSelectedPad(((TCanvas*)fCanvases->UncheckedAt(kDetailResPerCh+ib+2*ia))->cd(i+1));
1785  else
1786  gROOT->SetSelectedPad(((TCanvas*)fCanvases->UncheckedAt(kDetailResPerDE+ib+2*ia))->cd(i+1));
1787  gPad->SetLogy();
1788  h->Draw("hist");
1789  TF1* f = (TF1*)fGaus->DrawClone("same");
1790  f->SetLineWidth(1);
1791  f->SetLineColor(2);
1792  fTmpHists->AddLast(f);
1793  }
1794 
1795  } else { // take the mean of the distribution
1796 
1797  if (zoom) Zoom(h);
1798 
1799  mean = h->GetMean();
1800  meanErr = h->GetMeanError();
1801  rms = h->GetRMS();
1802  rmsErr = h->GetRMSError();
1803 
1804  if (zoom) h->GetXaxis()->SetRange(0,0);
1805 
1806  }
1807 
1808  // fill graph if required
1809  if (gMean) {
1810  gMean->SetPoint(i, x, mean);
1811  gMean->SetPointError(i, 0., meanErr);
1812  }
1813  if (gRMS) {
1814  gRMS->SetPoint(i, x, rms);
1815  gRMS->SetPointError(i, 0., rmsErr);
1816  }
1817 
1818 }
1819 
1820 //________________________________________________________________________
1821 void AliAnalysisTaskMuonResolution::FillMeanSigmaClusterVsX(const TH2* hIn, const TH2* hOut,
1822  TGraphErrors* gMean, TGraphErrors* gSigma)
1823 {
1826  Double_t meanIn, meanInErr, sigmaIn, sigmaInErr, sigmaOut, sigmaOutErr, clusterRes, clusterResErr, dumy1, dumy2;
1827  for (Int_t j = 1; j <= hIn->GetNbinsX(); j++) {
1828  TH1D* tmp = hIn->ProjectionY("tmp",j,j,"e");
1829  GetMeanRMS(tmp, meanIn, meanInErr, sigmaIn, sigmaInErr, 0x0, 0x0, 0, 0.);
1830  delete tmp;
1831  if (hOut) {
1832  tmp = hOut->ProjectionY("tmp",j,j,"e");
1833  GetMeanRMS(tmp, dumy1, dumy2, sigmaOut, sigmaOutErr, 0x0, 0x0, 0, 0.);
1834  delete tmp;
1835  }
1836  Double_t x = 0.5 * (hIn->GetXaxis()->GetBinLowEdge(j) + hIn->GetXaxis()->GetBinLowEdge(j+1));
1837  Double_t xErr = x - hIn->GetXaxis()->GetBinLowEdge(j);
1838  if (gMean) {
1839  gMean->SetPoint(j-1, x, meanIn);
1840  gMean->SetPointError(j-1, xErr, meanInErr);
1841  }
1842  if (hOut && gSigma) {
1843  clusterRes = TMath::Sqrt(sigmaIn*sigmaOut);
1844  //clusterResErr = (clusterRes > 0.) ? 0.5 * TMath::Sqrt(sigmaInErr*sigmaInErr*sigmaOut*sigmaOut + sigmaIn*sigmaIn*sigmaOutErr*sigmaOutErr) / clusterRes : 0.;
1845  clusterResErr = TMath::Sqrt(sigmaInErr*sigmaOutErr);
1846  gSigma->SetPoint(j-1, x, clusterRes);
1847  gSigma->SetPointError(j-1, xErr, clusterResErr);
1848  }
1849  }
1850 }
1851 
1852 //__________________________________________________________________________
1853 void AliAnalysisTaskMuonResolution::Cov2CovP(const AliMUONTrackParam &param, TMatrixD &covP)
1854 {
1857 
1858  // Get useful parameters
1859  Double_t slopeX = param.GetNonBendingSlope();
1860  Double_t slopeY = param.GetBendingSlope();
1861  Double_t qOverPYZ = param.GetInverseBendingMomentum();
1862  Double_t pZ = param.Pz();
1863 
1864  // compute Jacobian to change the coordinate system from (X,SlopeX,Y,SlopeY,c/pYZ) to (X,Y,pX,pY,pZ)
1865  Double_t dpZdSlopeY = - qOverPYZ * qOverPYZ * pZ * pZ * pZ * slopeY;
1866  Double_t dpZdQOverPYZ = (qOverPYZ != 0.) ? - pZ / qOverPYZ : - FLT_MAX;
1867  TMatrixD jacob(5,5);
1868  jacob.Zero();
1869  jacob(0,0) = 1.;
1870  jacob(1,2) = 1.;
1871  jacob(2,1) = pZ;
1872  jacob(2,3) = slopeX * dpZdSlopeY;
1873  jacob(2,4) = slopeX * dpZdQOverPYZ;
1874  jacob(3,3) = pZ + slopeY * dpZdSlopeY;
1875  jacob(3,4) = slopeY * dpZdQOverPYZ;
1876  jacob(4,3) = dpZdSlopeY;
1877  jacob(4,4) = dpZdQOverPYZ;
1878 
1879  // compute covariances in new coordinate system
1880  TMatrixD tmp(param.GetCovariances(),TMatrixD::kMultTranspose,jacob);
1881  covP.Mult(jacob,tmp);
1882 }
1883 
1884 //__________________________________________________________________________
1885 void AliAnalysisTaskMuonResolution::CheckPads(AliMUONVCluster *cl, Bool_t &hasBending, Bool_t &hasNonBending) const
1886 {
1888 
1889  // reset
1890  hasBending = kFALSE;
1891  hasNonBending = kFALSE;
1892 
1893  // loop over digits contained in the cluster
1894  for (Int_t iDigit = 0; iDigit < cl->GetNDigits(); iDigit++) {
1895 
1896  Int_t manuId = AliMUONVDigit::ManuId(cl->GetDigitId(iDigit));
1897 
1898  // check the location of the manu the digit belongs to
1899  if (manuId > 0) {
1900  if (manuId & AliMpConstants::ManuMask(AliMp::kNonBendingPlane)) hasNonBending = kTRUE;
1901  else hasBending = kTRUE;
1902  }
1903 
1904  }
1905 
1906 }
1907 
1908 //________________________________________________________________________
1909 void AliAnalysisTaskMuonResolution::CheckPadsBelow(AliMUONVCluster *cl, Bool_t &hasBending, Bool_t &hasNonBending) const
1910 {
1912 
1913  // reset
1914  hasBending = kFALSE;
1915  hasNonBending = kFALSE;
1916 
1917  // get the cathod corresponding to the bending/non-bending plane
1918  Int_t deId = cl->GetDetElemId();
1919  AliMpDetElement* de = AliMpDDLStore::Instance()->GetDetElement(deId, kFALSE);
1920  if (!de) return;
1921  AliMp::CathodType cath1 = de->GetCathodType(AliMp::kBendingPlane);
1922  AliMp::CathodType cath2 = de->GetCathodType(AliMp::kNonBendingPlane);
1923 
1924  // get the corresponding segmentation
1925  const AliMpVSegmentation* seg1 = AliMpSegmentation::Instance()->GetMpSegmentation(deId, cath1);
1926  const AliMpVSegmentation* seg2 = AliMpSegmentation::Instance()->GetMpSegmentation(deId, cath2);
1927  if (!seg1 || !seg2) return;
1928 
1929  // get global coordinate of the cluster
1930  Double_t gX = cl->GetX();
1931  Double_t gY = cl->GetY();
1932  Double_t gZ = cl->GetZ();
1933 
1934  // revert half-chamber or DE displacement if any
1935  Int_t chId = cl->GetChamberId();
1936  Int_t halfChId = (cl->GetX() > 0) ? 2*chId : 2*chId+1;
1937  if (fShiftHalfCh) {
1938  gX += fHalfChShiftNB[halfChId];
1939  gY += fHalfChShiftB[halfChId];
1940  }
1941  if (fShiftDE) {
1942  gX += fDEShiftNB[fDEIndices[deId]-1];
1943  gY += fDEShiftB[fDEIndices[deId]-1];
1944  }
1945 
1946  // get local coordinate of the cluster
1947  Double_t lX,lY,lZ;
1948  fNewGeoTransformer->Global2Local(deId,gX,gY,gZ,lX,lY,lZ);
1949 
1950  // find pads below the cluster
1951  AliMpPad pad1 = seg1->PadByPosition(lX, lY, kFALSE);
1952  AliMpPad pad2 = seg2->PadByPosition(lX, lY, kFALSE);
1953 
1954  // build their ID if pads are valid
1955  UInt_t padId1 = (pad1.IsValid()) ? AliMUONVDigit::BuildUniqueID(deId, pad1.GetManuId(), pad1.GetManuChannel(), cath1) : 0;
1956  UInt_t padId2 = (pad2.IsValid()) ? AliMUONVDigit::BuildUniqueID(deId, pad2.GetManuId(), pad2.GetManuChannel(), cath2) : 0;
1957 
1958  // check if the cluster contains these pads
1959  for (Int_t iDigit = 0; iDigit < cl->GetNDigits(); iDigit++) {
1960 
1961  UInt_t digitId = cl->GetDigitId(iDigit);
1962 
1963  if (digitId == padId1) {
1964 
1965  hasBending = kTRUE;
1966  if (hasNonBending) break;
1967 
1968  } else if (digitId == padId2) {
1969 
1970  hasNonBending = kTRUE;
1971  if (hasBending) break;
1972 
1973  }
1974 
1975  }
1976 
1977 }
1978 
Double_t fHalfChShiftB[20]
half-chamber deplacements in bending direction
muon slope-X/Y reconstructed resolution at first cluster vs p
cluster-track residual-X/Y distribution in chamber i versus momentum (cluster attached to the track) ...
muon momentum reconstructed resolution at first cluster vs p
cluster X/Y-resolution per chamber versus momentum
cluster-track residual-X/Y distribution in half-chamber i versus track angle in X/Y direction (cluste...
cluster-track residual-X/Y distribution integrated over chambers versus centrality (cluster not attac...
static const Int_t fgkMinEntries
minimum number of entries needed to compute resolution
AliMUONGeometryTransformer * fOldGeoTransformer
geometry transformer used to recontruct the present data
Bool_t fShiftDE
flag telling wether to displace DEs by fDEShift(N)B[i] or not
cluster-track residual-X/Y per DE: mean (cluster in)
const char * title
Definition: MakeQAPdf.C:26
cluster-track residual-X/Y distribution in chamber i versus centrality (cluster not attached to the t...
cluster-track residual-X/Y distribution per DE (cluster not attached to the track) ...
Bool_t fPrintClResPerCh
print the cluster resolution per chamber
cluster X/Y-resolution integrated over chambers versus centrality
cluster-track residual-X/Y per chamber: sigma (cluster out)
Bool_t fShowProgressBar
show the progression bar
TString fOldAlignStorage
location of the OCDB storage where to find old MUON/Align/Data (use the default one if empty) ...
TString fNewAlignStorage
location of the OCDB storage where to find new MUON/Align/Data (use the default one if empty) ...
Bool_t fShiftHalfCh
flag telling wether to displace half-chambers by fHalfChShift(N)B[i] or not
cluster X/Y-resolution integrated over chambers versus track angle in X/Y direction ...
Bool_t fImproveTracks
flag telling whether to improve or not the track before measuring the resolution
centrality
Bool_t fPrintHalfChShift
print the half-chamber displacements
cluster X/Y-resolution integrated over chambers versus momentum
AliMuonEventCuts * fMuonEventCuts
cuts to select events to be considered
cluster-track residual-X/Y distribution per chamber (cluster not attached to the track) ...
AliMuonTrackCuts * fMuonTrackCuts
cuts to select tracks to be considered
Bool_t fPrintDEShift
print the DE displacements
cluster-track residual-X/Y per chamber: mean (cluster out)
void ZoomLeft(TH1 *h, Double_t fractionCut=0.02)
Int_t fDEIds[200]
ID of DE refered by index in histograms.
cluster-track residual-X/Y distribution integrated over chambers versus track angle in X/Y direction ...
Double_t fClusterResNB[10]
cluster resolution in non-bending direction
void Zoom(TH1 *h, Double_t fractionCut=0.01)
Bool_t fRemoveMonoCathCl
remove or not the mono-cathod clusters
void FillMeanSigmaClusterVsX(const TH2 *hIn, const TH2 *hOut, TGraphErrors *gMean, TGraphErrors *gSigma)
cluster-track residual-X/Y distribution in chamber i versus track angle in X/Y direction (cluster att...
Double_t * sigma
TObjArray * fCanvases
List of canvases summarizing the results.
muon momentum reconstructed resolution at vertex vs p
Double_t fHalfChShiftNB[20]
half-chamber deplacements in non-bending direction
cluster-track residual-X/Y distribution in chamber i versus centrality (cluster attached to the track...
AliMUONGeometryTransformer * fNewGeoTransformer
new geometry transformer containing the new alignment to be applied
muon slope-X/Y reconstructed resolution at vertex vs p
Double_t fMinMomentum
use only tracks with momentum higher than this value
cluster-track residual-X/Y per chamber: dispersion (cluster out)
Bool_t fOCDBLoaded
flag telling if the OCDB has been properly loaded or not
void SetHalfChShift(Int_t hchId, Double_t valNB, Double_t valB)
TObjArray * fResidualsVsP
List of residual vs. p histos.
calculated cluster X/Y-resolution per chamber
void CheckPadsBelow(AliMUONVCluster *cl, Bool_t &hasBending, Bool_t &hasNonBending) const
cluster-track residual-X/Y distribution per half chamber (cluster not attached to the track) ...
cluster-track residual-X/Y distribution in chamber i versus momentum (cluster not attached to the tra...
combined cluster-track residual-X/Y per half chamber
void GetMeanRMS(TH1 *h, Double_t &mean, Double_t &meanErr, Double_t &rms, Double_t &rmsErr, TGraphErrors *gMean=0x0, TGraphErrors *gRMS=0x0, Int_t i=0, Double_t x=0, Bool_t zoom=kTRUE, Bool_t enableFit=kTRUE)
cluster-track residual-X/Y distribution integrated over chambers versus centrality (cluster attached ...
cluster-track residual-X/Y per half-chamber versus track angle in X/Y direction: mean (cluster in) ...
local chi2-X/Y/total distribution per chamber
cluster-track residual-X/Y per chamber versus centrality: mean (cluster in)
Int_t fNEvents
number of processed events
TObjArray * fResiduals
List of residual histos.
TObjArray * fResidualsVsAngle
List of residual vs. track angle histos.
muon transverse momentum reconstructed resolution at vertex vs p
cluster-track residual-X/Y distribution integrated over chambers versus momentum (cluster not attache...
TObjArray * fTrackRes
List of plots related to track resolution (p, pT, ...)
Int_t fOldAlignVersion
specific version of the old MUON/Align/Data/object to load
Double_t fMinPt
use only tracks with pT higher than this value
MCS X/Y-dispersion of extrapolated track per DE.
TObjArray * fLocalChi2
List of plots related to local chi2 per chamber/DE.
cluster-track residual-X/Y distribution per chamber (cluster attached to the track) ...
cluster-track residual-X/Y distribution in chamber i versus track angle in X/Y direction (cluster not...
TObjArray * fResidualsVsCent
List of residual vs. centrality histos.
Int_t fNewAlignSubVersion
specific subversion of the new MUON/Align/Data/object to load
cluster-track residual-X/Y per DE: mean (cluster out)
MCS X/Y-dispersion of extrapolated track per chamber.
Double_t fDEShiftB[200]
DE deplacements in bending direction.
Double_t fClusterResB[10]
cluster resolution in bending direction
cluster-track residual-X/Y per half chamber: mean (cluster in)
cluster-track residual-X/Y per chamber versus momentum: mean (cluster in)
cluster X/Y-resolution per chamber versus centrality
Muon spectrometer resolution.
Double_t fDEShiftNB[200]
DE deplacements in non-bending direction.
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
cluster-track residual-X/Y per chamber: sigma (cluster in)
MCS X/Y-dispersion of extrapolated track per chamber.
cluster-track residual-X/Y per chamber: mean (cluster in)
cluster-track residual-X/Y distribution integrated over chambers versus momentum (cluster attached to...
cluster-track residual-X/Y per half chamber: mean (cluster out)
Bool_t fCheckAllPads
use all pads or only the ones directly below the cluster to look for mono-cathods ...
Int_t rebin
void SetDEShift(Int_t iDE, Double_t valNB, Double_t valB)
TObjArray * fChamberRes
List of plots related to chamber/DE resolution.
cluster-track residual-X/Y per chamber versus track angle in X/Y direction: mean (cluster in) ...
cluster-track residual-X/Y distribution per half chamber (cluster attached to the track) ...
void ZoomRight(TH1 *h, Double_t fractionCut=0.02)
TString fDefaultStorage
location of the default OCDB storage
#define SafeDelete(x)
MCS X/Y-dispersion of extrapolated track per half chamber.
Bool_t fPrintClResPerDE
print the cluster resolution per DE
Int_t fOldAlignSubVersion
specific subversion of the old MUON/Align/Data/object to load
Bool_t fCorrectForSystematics
add or not the systematic shifts of the residuals to the resolution
Int_t fDEIndices[1100]
index of DE in histograms refered by ID
void Cov2CovP(const AliMUONTrackParam &param, TMatrixD &covP)
cluster X/Y-resolution per chamber versus track angle in X/Y direction
TObjArray * fTmpHists
List of temporary histograms.
TF1 * fGaus
gaussian function to fit the residuals
Int_t fNewAlignVersion
specific version of the new MUON/Align/Data/object to load
void SetStartingResolution(Int_t chId, Double_t valNB, Double_t valB)
muon transverse momentum reconstructed resolution at first cluster vs p
Bool_t fReAlign
flag telling whether to re-align the spectrometer or not before computing resolution ...
void CheckPads(AliMUONVCluster *cl, Bool_t &hasBending, Bool_t &hasNonBending) const
cluster-track residual-X/Y distribution per DE (cluster attached to the track)
Int_t fExtrapMode
extrapolation mode to get the track parameters and covariances at a given cluster ...
cluster-track residual-X/Y distribution integrated over chambers versus track angle in X/Y direction ...