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