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