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