Coding conventions
[u/mrichter/AliRoot.git] / PWG3 / muondep / AliAnalysisTaskMuonResolution.cxx
CommitLineData
7cbb1928 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
27de2dfb 16/* $Id$ */
17
7cbb1928 18// ROOT includes
19#include <TFile.h>
20#include <TTree.h>
21#include <TH1F.h>
22#include <TH2F.h>
23#include <TMultiGraph.h>
24#include <TGraphErrors.h>
25#include <TCanvas.h>
26#include <TLegend.h>
27#include <Riostream.h>
28#include <TString.h>
29#include <TGeoManager.h>
00612ffc 30#include <TList.h>
31#include <TObjString.h>
32#include <TRegexp.h>
7cbb1928 33
34// STEER includes
35#include "AliESDEvent.h"
36#include "AliESDMuonTrack.h"
37#include "AliCDBManager.h"
38#include "AliCDBStorage.h"
39#include "AliGeomManager.h"
40
41// ANALYSIS includes
7cbb1928 42#include "AliAnalysisDataSlot.h"
43#include "AliAnalysisManager.h"
44#include "AliInputEventHandler.h"
45#include "AliAnalysisTaskMuonResolution.h"
4060971a 46#include "AliCentrality.h"
7cbb1928 47
48// MUON includes
49#include "AliMUONCDB.h"
50#include "AliMUONRecoParam.h"
51#include "AliMUONESDInterface.h"
52#include "AliMUONVTrackReconstructor.h"
53#include "AliMUONTrack.h"
54#include "AliMUONTrackParam.h"
55#include "AliMUONTrackExtrap.h"
56#include "AliMUONVCluster.h"
4060971a 57#include "AliMUONVDigit.h"
7cbb1928 58#include "AliMUONGeometryTransformer.h"
59#include "AliMUONGeometryModuleTransformer.h"
60#include "AliMUONGeometryDetElement.h"
61#include "AliMpDEIterator.h"
4060971a 62#include "AliMpSegmentation.h"
63#include "AliMpVSegmentation.h"
64#include "AliMpConstants.h"
65#include "AliMpDDLStore.h"
66#include "AliMpPad.h"
67#include "AliMpDetElement.h"
68#include "AliMpCathodType.h"
7cbb1928 69
70#ifndef SafeDelete
71#define SafeDelete(x) if (x != NULL) { delete x; x = NULL; }
72#endif
73
74ClassImp(AliAnalysisTaskMuonResolution)
75
76const Int_t AliAnalysisTaskMuonResolution::fgkMinEntries = 10;
77
78//________________________________________________________________________
79AliAnalysisTaskMuonResolution::AliAnalysisTaskMuonResolution() :
80 AliAnalysisTaskSE(),
81 fResiduals(NULL),
82 fResidualsVsP(NULL),
4060971a 83 fResidualsVsCent(NULL),
7cbb1928 84 fLocalChi2(NULL),
85 fChamberRes(NULL),
86 fTrackRes(NULL),
87 fCanvases(NULL),
00612ffc 88 fDefaultStorage(""),
7cbb1928 89 fNEvents(0),
00612ffc 90 fShowProgressBar(kFALSE),
91 fPrintClResPerCh(kFALSE),
92 fPrintClResPerDE(kFALSE),
93 fGaus(NULL),
7cbb1928 94 fMinMomentum(0.),
95 fSelectPhysics(kFALSE),
96 fMatchTrig(kFALSE),
00612ffc 97 fApplyAccCut(kFALSE),
98 fSelectTrigger(kFALSE),
99 fTriggerMask(0),
7cbb1928 100 fExtrapMode(1),
101 fCorrectForSystematics(kTRUE),
4060971a 102 fRemoveMonoCathCl(kFALSE),
103 fCheckAllPads(kFALSE),
7cbb1928 104 fOCDBLoaded(kFALSE),
105 fNDE(0),
106 fReAlign(kFALSE),
107 fOldAlignStorage(""),
108 fNewAlignStorage(""),
109 fOldGeoTransformer(NULL),
00612ffc 110 fNewGeoTransformer(NULL),
111 fSelectTriggerClass(NULL)
7cbb1928 112{
113 /// Default constructor
114}
115
116//________________________________________________________________________
117AliAnalysisTaskMuonResolution::AliAnalysisTaskMuonResolution(const char *name) :
118 AliAnalysisTaskSE(name),
119 fResiduals(NULL),
120 fResidualsVsP(NULL),
4060971a 121 fResidualsVsCent(NULL),
7cbb1928 122 fLocalChi2(NULL),
123 fChamberRes(NULL),
124 fTrackRes(NULL),
125 fCanvases(NULL),
00612ffc 126 fDefaultStorage("raw://"),
7cbb1928 127 fNEvents(0),
00612ffc 128 fShowProgressBar(kFALSE),
129 fPrintClResPerCh(kFALSE),
130 fPrintClResPerDE(kFALSE),
131 fGaus(NULL),
7cbb1928 132 fMinMomentum(0.),
133 fSelectPhysics(kFALSE),
134 fMatchTrig(kFALSE),
00612ffc 135 fApplyAccCut(kFALSE),
136 fSelectTrigger(kFALSE),
137 fTriggerMask(0),
7cbb1928 138 fExtrapMode(1),
139 fCorrectForSystematics(kTRUE),
4060971a 140 fRemoveMonoCathCl(kFALSE),
141 fCheckAllPads(kFALSE),
7cbb1928 142 fOCDBLoaded(kFALSE),
143 fNDE(0),
144 fReAlign(kFALSE),
145 fOldAlignStorage(""),
146 fNewAlignStorage(""),
147 fOldGeoTransformer(NULL),
00612ffc 148 fNewGeoTransformer(NULL),
149 fSelectTriggerClass(NULL)
7cbb1928 150{
151 /// Constructor
152
153 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) SetStartingResolution(i, -1., -1.);
154
00612ffc 155 FitResiduals();
156
7cbb1928 157 // Output slot #1 writes into a TObjArray container
158 DefineOutput(1,TObjArray::Class());
159 // Output slot #2 writes into a TObjArray container
160 DefineOutput(2,TObjArray::Class());
161 // Output slot #3 writes into a TObjArray container
162 DefineOutput(3,TObjArray::Class());
163 // Output slot #4 writes into a TObjArray container
164 DefineOutput(4,TObjArray::Class());
165 // Output slot #5 writes into a TObjArray container
166 DefineOutput(5,TObjArray::Class());
4060971a 167 // Output slot #5 writes into a TObjArray container
168 DefineOutput(6,TObjArray::Class());
7cbb1928 169}
170
171//________________________________________________________________________
172AliAnalysisTaskMuonResolution::~AliAnalysisTaskMuonResolution()
173{
174 /// Destructor
00612ffc 175 if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
176 SafeDelete(fResiduals);
177 SafeDelete(fResidualsVsP);
4060971a 178 SafeDelete(fResidualsVsCent);
00612ffc 179 SafeDelete(fTrackRes);
180 }
7cbb1928 181 SafeDelete(fChamberRes);
7cbb1928 182 SafeDelete(fCanvases);
00612ffc 183 SafeDelete(fGaus);
7cbb1928 184 SafeDelete(fOldGeoTransformer);
185 SafeDelete(fNewGeoTransformer);
00612ffc 186 SafeDelete(fSelectTriggerClass);
7cbb1928 187}
188
189//___________________________________________________________________________
190void AliAnalysisTaskMuonResolution::UserCreateOutputObjects()
191{
192 /// Create histograms
193
194 // do it once the OCDB has been loaded (i.e. from NotifyRun())
195 if (!fOCDBLoaded) return;
196
00612ffc 197 // set the list of trigger classes that can be selected to fill histograms (in case the physics selection is not used)
198 fSelectTriggerClass = new TList();
199 fSelectTriggerClass->SetOwner();
200 fSelectTriggerClass->AddLast(new TObjString(" CINT1B-ABCE-NOPF-ALL ")); fSelectTriggerClass->Last()->SetUniqueID(AliVEvent::kMB);
201 fSelectTriggerClass->AddLast(new TObjString(" CMUS1B-ABCE-NOPF-MUON ")); fSelectTriggerClass->Last()->SetUniqueID(AliVEvent::kMUON);
202 fSelectTriggerClass->AddLast(new TObjString(" CINT1-B-")); fSelectTriggerClass->Last()->SetUniqueID(AliVEvent::kMB);
203 fSelectTriggerClass->AddLast(new TObjString(" CMUS1-B-")); fSelectTriggerClass->Last()->SetUniqueID(AliVEvent::kMUON);
204 fSelectTriggerClass->AddLast(new TObjString(" CSH1-B-")); fSelectTriggerClass->Last()->SetUniqueID(AliVEvent::kHighMult);
205
7cbb1928 206 fResiduals = new TObjArray(1000);
207 fResiduals->SetOwner();
208 fResidualsVsP = new TObjArray(1000);
209 fResidualsVsP->SetOwner();
4060971a 210 fResidualsVsCent = new TObjArray(1000);
211 fResidualsVsCent->SetOwner();
7cbb1928 212 fTrackRes = new TObjArray(1000);
213 fTrackRes->SetOwner();
214 TH2F* h2;
215
216 // find the highest chamber resolution and set histogram bins
217 const AliMUONRecoParam* recoParam = AliMUONESDInterface::GetTracker()->GetRecoParam();
218 Double_t maxSigma[2] = {-1., -1.};
219 for (Int_t i = 0; i < 10; i++) {
220 if (recoParam->GetDefaultNonBendingReso(i) > maxSigma[0]) maxSigma[0] = recoParam->GetDefaultNonBendingReso(i);
221 if (recoParam->GetDefaultBendingReso(i) > maxSigma[1]) maxSigma[1] = recoParam->GetDefaultBendingReso(i);
222 }
223 const char* axes[2] = {"X", "Y"};
00612ffc 224 const Int_t nBins = 5000;
7cbb1928 225 const Int_t nSigma = 10;
226 const Int_t pNBins = 20;
227 const Double_t pEdges[2] = {0., 50.};
4060971a 228 Int_t nCentBins = 12;
229 Double_t centRange[2] = {-10., 110.};
7cbb1928 230
231 for (Int_t ia = 0; ia < 2; ia++) {
232
233 Double_t maxRes = nSigma*maxSigma[ia];
234
235 // List of residual histos
236 h2 = new TH2F(Form("hResidual%sPerCh_ClusterIn",axes[ia]), Form("cluster-track residual-%s distribution per chamber (cluster attached to the track);chamber ID;#Delta_{%s} (cm)",axes[ia],axes[ia]), 10, 0.5, 10.5, nBins, -maxRes, maxRes);
f4d9245c 237 fResiduals->AddAtAndExpand(h2, kResidualPerChClusterIn+ia);
7cbb1928 238 h2 = new TH2F(Form("hResidual%sPerCh_ClusterOut",axes[ia]), Form("cluster-track residual-%s distribution per chamber (cluster not attached to the track);chamber ID;#Delta_{%s} (cm)",axes[ia],axes[ia]), 10, 0.5, 10.5, nBins, -2.*maxRes, 2.*maxRes);
f4d9245c 239 fResiduals->AddAtAndExpand(h2, kResidualPerChClusterOut+ia);
7cbb1928 240
241 h2 = new TH2F(Form("hResidual%sPerHalfCh_ClusterIn",axes[ia]), Form("cluster-track residual-%s distribution per half chamber (cluster attached to the track);half chamber ID;#Delta_{%s} (cm)",axes[ia],axes[ia]), 20, 0.5, 20.5, nBins, -maxRes, maxRes);
12d09c5b 242 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)); }
f4d9245c 243 fResiduals->AddAtAndExpand(h2, kResidualPerHalfChClusterIn+ia);
7cbb1928 244 h2 = new TH2F(Form("hResidual%sPerHalfCh_ClusterOut",axes[ia]), 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]), 20, 0.5, 20.5, nBins, -2.*maxRes, 2.*maxRes);
12d09c5b 245 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)); }
f4d9245c 246 fResiduals->AddAtAndExpand(h2, kResidualPerHalfChClusterOut+ia);
7cbb1928 247
4060971a 248 h2 = new TH2F(Form("hResidual%sPerDE_ClusterIn",axes[ia]), Form("cluster-track residual-%s distribution per DE (cluster attached to the track);DE ID;#Delta_{%s} (cm)",axes[ia],axes[ia]), fNDE, 0.5, fNDE+0.5, nBins, -maxRes, maxRes);
7cbb1928 249 for (Int_t i = 1; i <= fNDE; i++) h2->GetXaxis()->SetBinLabel(i, Form("%d",fDEIds[i]));
f4d9245c 250 fResiduals->AddAtAndExpand(h2, kResidualPerDEClusterIn+ia);
7cbb1928 251 h2 = new TH2F(Form("hResidual%sPerDE_ClusterOut",axes[ia]), Form("cluster-track residual-%s distribution per DE (cluster not attached to the track);DE ID;#Delta_{%s} (cm)",axes[ia],axes[ia]), fNDE, 0.5, fNDE+0.5, nBins, -2.*maxRes, 2.*maxRes);
252 for (Int_t i = 1; i <= fNDE; i++) h2->GetXaxis()->SetBinLabel(i, Form("%d",fDEIds[i]));
f4d9245c 253 fResiduals->AddAtAndExpand(h2, kResidualPerDEClusterOut+ia);
7cbb1928 254
255 h2 = new TH2F(Form("hTrackRes%sPerCh",axes[ia]), Form("track #sigma_{%s} per Ch;chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]), 10, 0.5, 10.5, nBins, 0., maxRes);
256 fResiduals->AddAtAndExpand(h2, kTrackResPerCh+ia);
257 h2 = new TH2F(Form("hTrackRes%sPerHalfCh",axes[ia]), Form("track #sigma_{%s} per half Ch;half chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]), 20, 0.5, 20.5, nBins, 0., maxRes);
12d09c5b 258 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)); }
7cbb1928 259 fResiduals->AddAtAndExpand(h2, kTrackResPerHalfCh+ia);
260 h2 = new TH2F(Form("hTrackRes%sPerDE",axes[ia]), Form("track #sigma_{%s} per DE;DE ID;#sigma_{%s} (cm)",axes[ia],axes[ia]), fNDE, 0.5, fNDE+0.5, nBins, 0., maxRes);
261 for (Int_t i = 1; i <= fNDE; i++) h2->GetXaxis()->SetBinLabel(i, Form("%d",fDEIds[i]));
262 fResiduals->AddAtAndExpand(h2, kTrackResPerDE+ia);
263
264 h2 = new TH2F(Form("hMCS%sPerCh",axes[ia]), Form("MCS %s-dispersion of extrapolated track per Ch;chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]), 10, 0.5, 10.5, nBins, 0., 0.2);
265 fResiduals->AddAtAndExpand(h2, kMCSPerCh+ia);
266 h2 = new TH2F(Form("hMCS%sPerHalfCh",axes[ia]), Form("MCS %s-dispersion of extrapolated track per half Ch;half chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]), 20, 0.5, 20.5, nBins, 0., 0.2);
12d09c5b 267 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)); }
7cbb1928 268 fResiduals->AddAtAndExpand(h2, kMCSPerHalfCh+ia);
269 h2 = new TH2F(Form("hMCS%sPerDE",axes[ia]), Form("MCS %s-dispersion of extrapolated track per DE;DE ID;#sigma_{%s} (cm)",axes[ia],axes[ia]), fNDE, 0.5, fNDE+0.5, nBins, 0., 0.2);
270 for (Int_t i = 1; i <= fNDE; i++) h2->GetXaxis()->SetBinLabel(i, Form("%d",fDEIds[i]));
271 fResiduals->AddAtAndExpand(h2, kMCSPerDE+ia);
272
273 h2 = new TH2F(Form("hClusterRes2%sPerCh",axes[ia]), Form("cluster #sigma_{%s}^{2} per Ch;chamber ID;#sigma_{%s}^{2} (cm^{2})",axes[ia],axes[ia]), 10, 0.5, 10.5, nSigma*nBins, -0.1*maxRes*maxRes, maxRes*maxRes);
274 fResiduals->AddAtAndExpand(h2, kClusterRes2PerCh+ia);
275
4060971a 276 // List of residual vs. p or vs. centrality histos
7cbb1928 277 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
278 h2 = new TH2F(Form("hResidual%sInCh%dVsP_ClusterIn",axes[ia],i+1), 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]), pNBins, pEdges[0], pEdges[1], nBins, -maxRes, maxRes);
f4d9245c 279 fResidualsVsP->AddAtAndExpand(h2, kResidualInChVsPClusterIn+10*ia+i);
7cbb1928 280 h2 = new TH2F(Form("hResidual%sInCh%dVsP_ClusterOut",axes[ia],i+1), 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]), pNBins, pEdges[0], pEdges[1], nBins, -2.*maxRes, 2.*maxRes);
f4d9245c 281 fResidualsVsP->AddAtAndExpand(h2, kResidualInChVsPClusterOut+10*ia+i);
4060971a 282
283 h2 = new TH2F(Form("hResidual%sInCh%dVsCent_ClusterIn",axes[ia],i+1), 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]), nCentBins, centRange[0], centRange[1], nBins, -maxRes, maxRes);
284 fResidualsVsCent->AddAtAndExpand(h2, kResidualInChVsCentClusterIn+10*ia+i);
285 h2 = new TH2F(Form("hResidual%sInCh%dVsCent_ClusterOut",axes[ia],i+1), 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]), nCentBins, centRange[0], centRange[1], nBins, -2.*maxRes, 2.*maxRes);
286 fResidualsVsCent->AddAtAndExpand(h2, kResidualInChVsCentClusterOut+10*ia+i);
7cbb1928 287 }
288
4060971a 289 h2 = new TH2F(Form("hResidual%sVsP_ClusterIn",axes[ia]), 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]), pNBins, pEdges[0], pEdges[1], nBins, -maxRes, maxRes);
290 fResidualsVsP->AddAtAndExpand(h2, kResidualVsPClusterIn+ia);
291 h2 = new TH2F(Form("hResidual%sVsP_ClusterOut",axes[ia]), 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]), pNBins, pEdges[0], pEdges[1], nBins, -2.*maxRes, 2.*maxRes);
292 fResidualsVsP->AddAtAndExpand(h2, kResidualVsPClusterOut+ia);
293
294 h2 = new TH2F(Form("hResidual%sVsCent_ClusterIn",axes[ia]), Form("cluster-track residual-%s distribution integrated over chambers versus centrality (cluster attached to the track);centrality (%%);#Delta_{%s} (cm)",axes[ia],axes[ia]), nCentBins, centRange[0], centRange[1], nBins, -maxRes, maxRes);
295 fResidualsVsCent->AddAtAndExpand(h2, kResidualVsCentClusterIn+ia);
296 h2 = new TH2F(Form("hResidual%sVsCent_ClusterOut",axes[ia]), 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]), nCentBins, centRange[0], centRange[1], nBins, -2.*maxRes, 2.*maxRes);
297 fResidualsVsCent->AddAtAndExpand(h2, kResidualVsCentClusterOut+ia);
298
7cbb1928 299 // local chi2
300 h2 = new TH2F(Form("hLocalChi2%sPerCh",axes[ia]), Form("local chi2-%s distribution per chamber;chamber ID;local #chi^{2}_{%s}", axes[ia], axes[ia]), 10, 0.5, 10.5, 1000, 0., 25.);
301 fResiduals->AddAtAndExpand(h2, kLocalChi2PerCh+ia);
302 h2 = new TH2F(Form("hLocalChi2%sPerDE",axes[ia]), Form("local chi2-%s distribution per DE;DE ID;local #chi^{2}_{%s}", axes[ia], axes[ia]), fNDE, 0.5, fNDE+0.5, 1000, 0., 25.);
303 for (Int_t i = 1; i <= fNDE; i++) h2->GetXaxis()->SetBinLabel(i, Form("%d",fDEIds[i]));
304 fResiduals->AddAtAndExpand(h2, kLocalChi2PerDE+ia);
305
306 // track resolution
307 h2 = new TH2F(Form("hUncorrSlope%sRes",axes[ia]), Form("muon slope_{%s} reconstructed resolution at first cluster vs p;p (GeV/c); #sigma_{slope_{%s}}", axes[ia], axes[ia]), 300, 0., 300., 1000, 0., 0.003);
308 fTrackRes->AddAtAndExpand(h2, kUncorrSlopeRes+ia);
309 h2 = new TH2F(Form("hSlope%sRes",axes[ia]), Form("muon slope_{%s} reconstructed resolution at vertex vs p;p (GeV/c); #sigma_{slope_{%s}}", axes[ia], axes[ia]), 300, 0., 300., 1000, 0., 0.02);
310 fTrackRes->AddAtAndExpand(h2, kSlopeRes+ia);
311 }
312
313 // local chi2 X+Y
314 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.);
315 fResiduals->AddAtAndExpand(h2, kLocalChi2PerCh+2);
316 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.);
317 for (Int_t i = 1; i <= fNDE; i++) h2->GetXaxis()->SetBinLabel(i, Form("%d",fDEIds[i]));
318 fResiduals->AddAtAndExpand(h2, kLocalChi2PerDE+2);
319
320 // track resolution
321 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.);
322 fTrackRes->AddAtAndExpand(h2, kUncorrPRes);
323 h2 = new TH2F("hPRes", "muon momentum reconstructed resolution at vertex vs p;p (GeV/c); #sigma_{p}/p (%)", 300, 0., 300., 1000, 0., 10.);
324 fTrackRes->AddAtAndExpand(h2, kPRes);
325 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.);
326 fTrackRes->AddAtAndExpand(h2, kUncorrPtRes);
327 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.);
328 fTrackRes->AddAtAndExpand(h2, kPtRes);
329
330 // Post data at least once per task to ensure data synchronisation (required for merging)
331 PostData(1, fResiduals);
332 PostData(2, fResidualsVsP);
4060971a 333 PostData(3, fResidualsVsCent);
334 PostData(6, fTrackRes);
7cbb1928 335}
336
337//________________________________________________________________________
338void AliAnalysisTaskMuonResolution::UserExec(Option_t *)
339{
340 /// Main event loop
341
342 // check if OCDB properly loaded
343 if (!fOCDBLoaded) return;
344
345 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
346 if (!esd) return;
347
00612ffc 348 if (fShowProgressBar && (++fNEvents)%100 == 0) cout<<"\rEvent processing... "<<fNEvents<<"\r"<<flush;
349
350 // skip events that do not pass the physics selection if required
351 UInt_t triggerWord = (fInputHandler) ? fInputHandler->IsEventSelected() : 0;
352 if (fSelectPhysics && triggerWord == 0) return;
353
354 // skip events that do not pass the trigger selection if required
f4d9245c 355 TString firedTriggerClasses = esd->GetFiredTriggerClasses();
356 if (!fSelectPhysics) triggerWord = BuildTriggerWord(firedTriggerClasses);
00612ffc 357 if (fSelectTrigger && (triggerWord & fTriggerMask) == 0) return;
7cbb1928 358
4060971a 359 // get the centrality
360 Float_t centrality = esd->GetCentrality()->GetCentralityPercentileUnchecked("V0M");
361
7cbb1928 362 // get tracker to refit
363 AliMUONVTrackReconstructor* tracker = AliMUONESDInterface::GetTracker();
364
365 // loop over tracks
366 Int_t nTracks = (Int_t) esd->GetNumberOfMuonTracks();
367 for (Int_t iTrack = 0; iTrack < nTracks; ++iTrack) {
368
369 // get the ESD track
370 AliESDMuonTrack* esdTrack = esd->GetMuonTrack(iTrack);
371
372 // skip ghost tracks
373 if (!esdTrack->ContainTrackerData()) continue;
374
7cbb1928 375 // skip tracks not matched with trigger if required
376 if (fMatchTrig && !esdTrack->ContainTriggerData()) continue;
377
00612ffc 378 // skip tracks that do not pass the acceptance cuts if required
379 Double_t thetaAbs = TMath::ATan(esdTrack->GetRAtAbsorberEnd()/505.) * TMath::RadToDeg();
380 Double_t eta = esdTrack->Eta();
4060971a 381 if (fApplyAccCut && (thetaAbs < 2. || thetaAbs > 10. || eta < -4. || eta > -2.5)) continue;
00612ffc 382
7cbb1928 383 // skip low momentum tracks
384 if (esdTrack->PUncorrected() < fMinMomentum) continue;
385
386 // get the corresponding MUON track
387 AliMUONTrack track;
388 AliMUONESDInterface::ESDToMUON(*esdTrack, track, kFALSE);
389
390 // change the cluster resolution (and position)
391 ModifyClusters(track);
392
393 // refit the track
394 if (!tracker->RefitTrack(track, kFALSE)) break;
395
396 // save track unchanged
397 AliMUONTrack referenceTrack(track);
398
399 // get track param at first cluster and add MCS in first chamber
400 AliMUONTrackParam trackParamAtFirstCluster(*(static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->First())));
401 Int_t firstCh = 0; while (firstCh < 10 && !esdTrack->IsInMuonClusterMap(firstCh)) firstCh++;
402 AliMUONTrackExtrap::AddMCSEffect(&trackParamAtFirstCluster, AliMUONConstants::ChamberThicknessInX0(firstCh)/2., -1.);
403
404 // fill momentum error at first cluster
405 Double_t pXUncorr = trackParamAtFirstCluster.Px();
406 Double_t pYUncorr = trackParamAtFirstCluster.Py();
407 Double_t pZUncorr = trackParamAtFirstCluster.Pz();
408 Double_t pUncorr = trackParamAtFirstCluster.P();
409 TMatrixD covUncorr(5,5);
410 Cov2CovP(trackParamAtFirstCluster,covUncorr);
411 Double_t sigmaPUncorr = TMath::Sqrt(pXUncorr * (pXUncorr*covUncorr(2,2) + pYUncorr*covUncorr(2,3) + pZUncorr*covUncorr(2,4)) +
412 pYUncorr * (pXUncorr*covUncorr(3,2) + pYUncorr*covUncorr(3,3) + pZUncorr*covUncorr(3,4)) +
413 pZUncorr * (pXUncorr*covUncorr(4,2) + pYUncorr*covUncorr(4,3) + pZUncorr*covUncorr(4,4))) / pUncorr;
414 ((TH2F*)fTrackRes->UncheckedAt(kUncorrPRes))->Fill(pUncorr,100.*sigmaPUncorr/pUncorr);
415
416 // fill transverse momentum error at first cluster
417 Double_t ptUncorr = TMath::Sqrt(pXUncorr*pXUncorr + pYUncorr*pYUncorr);
418 Double_t sigmaPtUncorr = TMath::Sqrt(pXUncorr * (pXUncorr*covUncorr(2,2) + pYUncorr*covUncorr(2,3)) + pYUncorr * (pXUncorr*covUncorr(3,2) + pYUncorr*covUncorr(3,3))) / ptUncorr;
419 ((TH2F*)fTrackRes->UncheckedAt(kUncorrPtRes))->Fill(ptUncorr,100.*sigmaPtUncorr/ptUncorr);
420
421 // fill slopeX-Y error at first cluster
422 ((TH2F*)fTrackRes->UncheckedAt(kUncorrSlopeRes))->Fill(pUncorr,TMath::Sqrt(trackParamAtFirstCluster.GetCovariances()(1,1)));
423 ((TH2F*)fTrackRes->UncheckedAt(kUncorrSlopeRes+1))->Fill(pUncorr,TMath::Sqrt(trackParamAtFirstCluster.GetCovariances()(3,3)));
424
425 // fill momentum error at vertex
426 AliMUONTrackParam trackParamAtVtx(trackParamAtFirstCluster);
427 AliMUONTrackExtrap::ExtrapToVertex(&trackParamAtVtx, esdTrack->GetNonBendingCoor(), esdTrack->GetBendingCoor(), esdTrack->GetZ(), 0., 0.);
428 Double_t pXVtx = trackParamAtVtx.Px();
429 Double_t pYVtx = trackParamAtVtx.Py();
430 Double_t pZVtx = trackParamAtVtx.Pz();
431 Double_t pVtx = trackParamAtVtx.P();
432 TMatrixD covVtx(5,5);
433 Cov2CovP(trackParamAtVtx,covVtx);
434 Double_t sigmaPVtx = TMath::Sqrt(pXVtx * (pXVtx*covVtx(2,2) + pYVtx*covVtx(2,3) + pZVtx*covVtx(2,4)) +
435 pYVtx * (pXVtx*covVtx(3,2) + pYVtx*covVtx(3,3) + pZVtx*covVtx(3,4)) +
436 pZVtx * (pXVtx*covVtx(4,2) + pYVtx*covVtx(4,3) + pZVtx*covVtx(4,4))) / pVtx;
437 ((TH2F*)fTrackRes->UncheckedAt(kPRes))->Fill(pVtx,100.*sigmaPVtx/pVtx);
438
439 // fill transverse momentum error at vertex
440 Double_t ptVtx = TMath::Sqrt(pXVtx*pXVtx + pYVtx*pYVtx);
441 Double_t sigmaPtVtx = TMath::Sqrt(pXVtx * (pXVtx*covVtx(2,2) + pYVtx*covVtx(2,3)) + pYVtx * (pXVtx*covVtx(3,2) + pYVtx*covVtx(3,3))) / ptVtx;
442 ((TH2F*)fTrackRes->UncheckedAt(kPtRes))->Fill(ptVtx,100.*sigmaPtVtx/ptVtx);
443
444 // fill slopeX-Y error at vertex
445 ((TH2F*)fTrackRes->UncheckedAt(kSlopeRes))->Fill(pVtx,TMath::Sqrt(trackParamAtVtx.GetCovariances()(1,1)));
446 ((TH2F*)fTrackRes->UncheckedAt(kSlopeRes+1))->Fill(pVtx,TMath::Sqrt(trackParamAtVtx.GetCovariances()(3,3)));
447
448 // loop over clusters
449 Int_t nClusters = track.GetNClusters();
450 for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
451
452 // Get current, previous and next trackParams
453 AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->UncheckedAt(iCluster));
454 AliMUONTrackParam* previousTrackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->Before(trackParam));
455 AliMUONTrackParam* nextTrackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->After(trackParam));
2426f648 456 if (!previousTrackParam && !nextTrackParam) {
457 AliError(Form("unable to find a cluster neither before nor after the one at position %d !?!", iCluster));
458 track.RecursiveDump();
459 break;
460 }
7cbb1928 461
4060971a 462 // remove mono-cathod clusters if required
463 if (fRemoveMonoCathCl) {
464 Bool_t hasBending, hasNonBending;
465 if (fCheckAllPads) CheckPads(trackParam->GetClusterPtr(), hasBending, hasNonBending);
466 else CheckPadsBelow(trackParam->GetClusterPtr(), hasBending, hasNonBending);
467 if (!hasBending || !hasNonBending) continue;
468 }
469
7cbb1928 470 // save current trackParam and remove it from the track
471 AliMUONTrackParam currentTrackParam(*trackParam);
472 track.RemoveTrackParamAtCluster(trackParam);
473
474 // get cluster info
475 AliMUONVCluster* cluster = currentTrackParam.GetClusterPtr();
476 Int_t chId = cluster->GetChamberId();
477 Int_t halfChId = (cluster->GetX() > 0) ? 2*chId : 2*chId+1;
478 Int_t deId = cluster->GetDetElemId();
479
480 // compute residuals with cluster still attached to the track
481 AliMUONTrackParam* referenceTrackParam = static_cast<AliMUONTrackParam*>(referenceTrack.GetTrackParamAtCluster()->UncheckedAt(iCluster));
482 Double_t deltaX = cluster->GetX() - referenceTrackParam->GetNonBendingCoor();
483 Double_t deltaY = cluster->GetY() - referenceTrackParam->GetBendingCoor();
484
485 // compute local chi2
486 Double_t sigmaDeltaX2 = cluster->GetErrX2() - referenceTrackParam->GetCovariances()(0,0);
487 Double_t sigmaDeltaY2 = cluster->GetErrY2() - referenceTrackParam->GetCovariances()(2,2);
488 Double_t localChi2X = (sigmaDeltaX2 > 0.) ? deltaX*deltaX/sigmaDeltaX2 : 0.;
489 Double_t localChi2Y = (sigmaDeltaY2 > 0.) ? deltaY*deltaY/sigmaDeltaY2 : 0.;
490 Double_t localChi2 = 0.5 * referenceTrackParam->GetLocalChi2();
491
492 // fill local chi2 info at every clusters
493 ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerCh))->Fill(chId+1, localChi2X);
494 ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerCh+1))->Fill(chId+1, localChi2Y);
495 ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerCh+2))->Fill(chId+1, localChi2);
496 ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerDE))->Fill(fDEIndices[deId], localChi2X);
497 ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerDE+1))->Fill(fDEIndices[deId], localChi2Y);
498 ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerDE+2))->Fill(fDEIndices[deId], localChi2);
499
500 // make sure the track has another cluster in the same station and can still be refitted
501 Bool_t refit = track.IsValid( 1 << (chId/2) );
502 if (refit) {
503
504 // refit the track and proceed if everything goes fine
505 if (tracker->RefitTrack(track, kFALSE)) {
506
507 // fill histograms of residuals with cluster still attached to the track
f4d9245c 508 ((TH2F*)fResiduals->UncheckedAt(kResidualPerChClusterIn))->Fill(chId+1, deltaX);
509 ((TH2F*)fResiduals->UncheckedAt(kResidualPerChClusterIn+1))->Fill(chId+1, deltaY);
510 ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfChClusterIn))->Fill(halfChId+1, deltaX);
511 ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfChClusterIn+1))->Fill(halfChId+1, deltaY);
512 ((TH2F*)fResiduals->UncheckedAt(kResidualPerDEClusterIn))->Fill(fDEIndices[deId], deltaX);
513 ((TH2F*)fResiduals->UncheckedAt(kResidualPerDEClusterIn+1))->Fill(fDEIndices[deId], deltaY);
514 ((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsPClusterIn+chId))->Fill(pUncorr, deltaX);
515 ((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsPClusterIn+10+chId))->Fill(pUncorr, deltaY);
4060971a 516 ((TH2F*)fResidualsVsP->UncheckedAt(kResidualVsPClusterIn))->Fill(pUncorr, deltaX);
517 ((TH2F*)fResidualsVsP->UncheckedAt(kResidualVsPClusterIn+1))->Fill(pUncorr, deltaY);
518 ((TH2F*)fResidualsVsCent->UncheckedAt(kResidualInChVsCentClusterIn+chId))->Fill(centrality, deltaX);
519 ((TH2F*)fResidualsVsCent->UncheckedAt(kResidualInChVsCentClusterIn+10+chId))->Fill(centrality, deltaY);
520 ((TH2F*)fResidualsVsCent->UncheckedAt(kResidualVsCentClusterIn))->Fill(centrality, deltaX);
521 ((TH2F*)fResidualsVsCent->UncheckedAt(kResidualVsCentClusterIn+1))->Fill(centrality, deltaY);
7cbb1928 522
523 // find the track parameters closest to the current cluster position
524 Double_t dZWithPrevious = (previousTrackParam) ? TMath::Abs(previousTrackParam->GetClusterPtr()->GetZ() - cluster->GetZ()) : FLT_MAX;
525 Int_t previousChId = (previousTrackParam) ? previousTrackParam->GetClusterPtr()->GetChamberId() : -1;
526 Double_t dZWithNext = (nextTrackParam) ? TMath::Abs(nextTrackParam->GetClusterPtr()->GetZ() - cluster->GetZ()) : FLT_MAX;
12d09c5b 527 AliMUONTrackParam* startingTrackParam = (nextTrackParam) ? nextTrackParam : previousTrackParam;
528 if ((fExtrapMode == 0 && previousTrackParam && dZWithPrevious < dZWithNext) ||
7cbb1928 529 (fExtrapMode == 1 && previousTrackParam && !(chId/2 == 2 && previousChId/2 == 1) &&
530 !(chId/2 == 3 && previousChId/2 == 2))) startingTrackParam = previousTrackParam;
7cbb1928 531
532 // reset current parameters
533 currentTrackParam.SetParameters(startingTrackParam->GetParameters());
534 currentTrackParam.SetZ(startingTrackParam->GetZ());
535 currentTrackParam.SetCovariances(startingTrackParam->GetCovariances());
536 currentTrackParam.ResetPropagator();
537
538 // extrapolate to the current cluster position and fill histograms of residuals if everything goes fine
539 if (AliMUONTrackExtrap::ExtrapToZCov(&currentTrackParam, currentTrackParam.GetClusterPtr()->GetZ(), kTRUE)) {
540
541 // compute MCS dispersion on the first chamber
542 TMatrixD mcsCov(5,5);
543 if (startingTrackParam == nextTrackParam && chId == 0) {
544 AliMUONTrackParam trackParamForMCS;
545 trackParamForMCS.SetParameters(nextTrackParam->GetParameters());
546 AliMUONTrackExtrap::AddMCSEffect(&trackParamForMCS,AliMUONConstants::ChamberThicknessInX0(nextTrackParam->GetClusterPtr()->GetChamberId()),-1.);
547 const TMatrixD &propagator = currentTrackParam.GetPropagator();
548 TMatrixD tmp(trackParamForMCS.GetCovariances(),TMatrixD::kMultTranspose,propagator);
549 mcsCov.Mult(propagator,tmp);
550 } else mcsCov.Zero();
551
552 // compute residuals
553 Double_t trackResX2 = currentTrackParam.GetCovariances()(0,0) + mcsCov(0,0);
554 Double_t trackResY2 = currentTrackParam.GetCovariances()(2,2) + mcsCov(2,2);
555 deltaX = cluster->GetX() - currentTrackParam.GetNonBendingCoor();
556 deltaY = cluster->GetY() - currentTrackParam.GetBendingCoor();
557
558 // fill histograms with cluster not attached to the track
f4d9245c 559 ((TH2F*)fResiduals->UncheckedAt(kResidualPerChClusterOut))->Fill(chId+1, deltaX);
560 ((TH2F*)fResiduals->UncheckedAt(kResidualPerChClusterOut+1))->Fill(chId+1, deltaY);
561 ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfChClusterOut))->Fill(halfChId+1, deltaX);
562 ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfChClusterOut+1))->Fill(halfChId+1, deltaY);
563 ((TH2F*)fResiduals->UncheckedAt(kResidualPerDEClusterOut))->Fill(fDEIndices[deId], deltaX);
564 ((TH2F*)fResiduals->UncheckedAt(kResidualPerDEClusterOut+1))->Fill(fDEIndices[deId], deltaY);
565 ((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsPClusterOut+chId))->Fill(pUncorr, deltaX);
566 ((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsPClusterOut+10+chId))->Fill(pUncorr, deltaY);
4060971a 567 ((TH2F*)fResidualsVsP->UncheckedAt(kResidualVsPClusterOut))->Fill(pUncorr, deltaX);
568 ((TH2F*)fResidualsVsP->UncheckedAt(kResidualVsPClusterOut+1))->Fill(pUncorr, deltaY);
569 ((TH2F*)fResidualsVsCent->UncheckedAt(kResidualInChVsCentClusterOut+chId))->Fill(centrality, deltaX);
570 ((TH2F*)fResidualsVsCent->UncheckedAt(kResidualInChVsCentClusterOut+10+chId))->Fill(centrality, deltaY);
571 ((TH2F*)fResidualsVsCent->UncheckedAt(kResidualVsCentClusterOut))->Fill(centrality, deltaX);
572 ((TH2F*)fResidualsVsCent->UncheckedAt(kResidualVsCentClusterOut+1))->Fill(centrality, deltaY);
7cbb1928 573 ((TH2F*)fResiduals->UncheckedAt(kTrackResPerCh))->Fill(chId+1, TMath::Sqrt(trackResX2));
574 ((TH2F*)fResiduals->UncheckedAt(kTrackResPerCh+1))->Fill(chId+1, TMath::Sqrt(trackResY2));
575 ((TH2F*)fResiduals->UncheckedAt(kTrackResPerHalfCh))->Fill(halfChId+1, TMath::Sqrt(trackResX2));
576 ((TH2F*)fResiduals->UncheckedAt(kTrackResPerHalfCh+1))->Fill(halfChId+1, TMath::Sqrt(trackResY2));
577 ((TH2F*)fResiduals->UncheckedAt(kTrackResPerDE))->Fill(fDEIndices[deId], TMath::Sqrt(trackResX2));
578 ((TH2F*)fResiduals->UncheckedAt(kTrackResPerDE+1))->Fill(fDEIndices[deId], TMath::Sqrt(trackResY2));
579 ((TH2F*)fResiduals->UncheckedAt(kMCSPerCh))->Fill(chId+1, TMath::Sqrt(mcsCov(0,0)));
580 ((TH2F*)fResiduals->UncheckedAt(kMCSPerCh+1))->Fill(chId+1, TMath::Sqrt(mcsCov(2,2)));
581 ((TH2F*)fResiduals->UncheckedAt(kMCSPerHalfCh))->Fill(halfChId+1, TMath::Sqrt(mcsCov(0,0)));
582 ((TH2F*)fResiduals->UncheckedAt(kMCSPerHalfCh+1))->Fill(halfChId+1, TMath::Sqrt(mcsCov(2,2)));
583 ((TH2F*)fResiduals->UncheckedAt(kMCSPerDE))->Fill(fDEIndices[deId], TMath::Sqrt(mcsCov(0,0)));
584 ((TH2F*)fResiduals->UncheckedAt(kMCSPerDE+1))->Fill(fDEIndices[deId], TMath::Sqrt(mcsCov(2,2)));
585 ((TH2F*)fResiduals->UncheckedAt(kClusterRes2PerCh))->Fill(chId+1, deltaX*deltaX - trackResX2);
586 ((TH2F*)fResiduals->UncheckedAt(kClusterRes2PerCh+1))->Fill(chId+1, deltaY*deltaY - trackResY2);
587 }
588
589 }
590
591 }
592
593 // restore the track
594 track.AddTrackParamAtCluster(currentTrackParam, *(currentTrackParam.GetClusterPtr()), kTRUE);
595
596 }
597
598 }
599
600 // Post final data. It will be written to a file with option "RECREATE"
601 PostData(1, fResiduals);
602 PostData(2, fResidualsVsP);
4060971a 603 PostData(3, fResidualsVsCent);
604 PostData(6, fTrackRes);
7cbb1928 605}
606
607//________________________________________________________________________
608void AliAnalysisTaskMuonResolution::NotifyRun()
609{
610 /// load necessary data from OCDB corresponding to the first run number and initialize analysis
611
612 if (fOCDBLoaded) return;
613
614 AliCDBManager* cdbm = AliCDBManager::Instance();
00612ffc 615 cdbm->SetDefaultStorage(fDefaultStorage.Data());
7cbb1928 616 cdbm->SetRun(fCurrentRunNumber);
617
618 if (!AliMUONCDB::LoadField()) return;
619
620 if (!AliMUONCDB::LoadMapping()) return;
621
622 AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam();
623 if (!recoParam) return;
624
7cbb1928 625 AliMUONESDInterface::ResetTracker(recoParam);
626
627 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
628
629 // set the cluster resolution to default if not already set and create output objects
630 if (fClusterResNB[i] < 0.) fClusterResNB[i] = recoParam->GetDefaultNonBendingReso(i);
631 if (fClusterResB[i] < 0.) fClusterResB[i] = recoParam->GetDefaultBendingReso(i);
632
633 // fill correspondence between DEId and indices for histo (starting from 1)
634 AliMpDEIterator it;
635 it.First(i);
636 while (!it.IsDone()) {
637 fNDE++;
638 fDEIndices[it.CurrentDEId()] = fNDE;
639 fDEIds[fNDE] = it.CurrentDEId();
640 it.Next();
641 }
642
643 }
644
645 if (fReAlign) {
646
00612ffc 647 // recover default storage full name (raw:// cannot be used to set specific storage)
7cbb1928 648 TString defaultStorage(cdbm->GetDefaultStorage()->GetType());
649 if (defaultStorage == "alien") defaultStorage += Form("://folder=%s", cdbm->GetDefaultStorage()->GetBaseFolder().Data());
650 else defaultStorage += Form("://%s", cdbm->GetDefaultStorage()->GetBaseFolder().Data());
651
652 // reset existing geometry/alignment if any
653 if (cdbm->GetEntryCache()->Contains("GRP/Geometry/Data")) cdbm->UnloadFromCache("GRP/Geometry/Data");
654 if (cdbm->GetEntryCache()->Contains("MUON/Align/Data")) cdbm->UnloadFromCache("MUON/Align/Data");
655 if (AliGeomManager::GetGeometry()) AliGeomManager::GetGeometry()->UnlockGeometry();
656
657 // get original geometry transformer
658 AliGeomManager::LoadGeometry();
659 if (!AliGeomManager::GetGeometry()) return;
660 if (fOldAlignStorage != "none") {
661 if (!fOldAlignStorage.IsNull()) cdbm->SetSpecificStorage("MUON/Align/Data",fOldAlignStorage.Data());
662 else cdbm->SetSpecificStorage("MUON/Align/Data",defaultStorage.Data());
663 AliGeomManager::ApplyAlignObjsFromCDB("MUON");
664 }
665 fOldGeoTransformer = new AliMUONGeometryTransformer();
666 fOldGeoTransformer->LoadGeometryData();
667
668 // get new geometry transformer
669 cdbm->UnloadFromCache("GRP/Geometry/Data");
670 if (fOldAlignStorage != "none") cdbm->UnloadFromCache("MUON/Align/Data");
671 AliGeomManager::GetGeometry()->UnlockGeometry();
672 AliGeomManager::LoadGeometry();
673 if (!AliGeomManager::GetGeometry()) return;
674 if (!fNewAlignStorage.IsNull()) cdbm->SetSpecificStorage("MUON/Align/Data",fNewAlignStorage.Data());
675 else cdbm->SetSpecificStorage("MUON/Align/Data",defaultStorage.Data());
676 AliGeomManager::ApplyAlignObjsFromCDB("MUON");
677 fNewGeoTransformer = new AliMUONGeometryTransformer();
678 fNewGeoTransformer->LoadGeometryData();
679
680 } else {
681
682 // load geometry for track extrapolation to vertex
683 if (cdbm->GetEntryCache()->Contains("GRP/Geometry/Data")) cdbm->UnloadFromCache("GRP/Geometry/Data");
4060971a 684 if (cdbm->GetEntryCache()->Contains("MUON/Align/Data")) cdbm->UnloadFromCache("MUON/Align/Data");
7cbb1928 685 if (AliGeomManager::GetGeometry()) AliGeomManager::GetGeometry()->UnlockGeometry();
686 AliGeomManager::LoadGeometry();
687 if (!AliGeomManager::GetGeometry()) return;
4060971a 688 AliGeomManager::ApplyAlignObjsFromCDB("MUON");
689 fNewGeoTransformer = new AliMUONGeometryTransformer();
690 fNewGeoTransformer->LoadGeometryData();
7cbb1928 691
692 }
693
00612ffc 694 // print starting chamber resolution if required
695 if (fPrintClResPerCh) {
696 printf("\nstarting chamber resolution:\n");
697 printf(" - non-bending:");
698 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %5.3f":", %5.3f",fClusterResNB[i]);
699 printf("\n - bending:");
700 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %6.4f":", %6.4f",fClusterResB[i]);
701 printf("\n\n");
702 }
703
704 fOCDBLoaded = kTRUE;
7cbb1928 705
706 UserCreateOutputObjects();
707
708}
709
710//________________________________________________________________________
711void AliAnalysisTaskMuonResolution::Terminate(Option_t *)
712{
713 /// compute final results
714
715 // recover output objects
716 fResiduals = static_cast<TObjArray*> (GetOutputData(1));
717 fResidualsVsP = static_cast<TObjArray*> (GetOutputData(2));
4060971a 718 fResidualsVsCent = static_cast<TObjArray*> (GetOutputData(3));
719 fTrackRes = static_cast<TObjArray*> (GetOutputData(6));
7cbb1928 720 if (!fResiduals || !fResidualsVsP || !fTrackRes) return;
721
722 // summary graphs
723 fLocalChi2 = new TObjArray(1000);
724 fLocalChi2->SetOwner();
725 fChamberRes = new TObjArray(1000);
726 fChamberRes->SetOwner();
727 TGraphErrors* g;
728 TMultiGraph* mg;
729
730 const char* axes[2] = {"X", "Y"};
731 Double_t newClusterRes[2][10], newClusterResErr[2][10];
f4d9245c 732 fNDE = ((TH2F*)fResiduals->UncheckedAt(kResidualPerDEClusterIn))->GetXaxis()->GetNbins();
7cbb1928 733
734 for (Int_t ia = 0; ia < 2; ia++) {
735
736 g = new TGraphErrors(AliMUONConstants::NTrackingCh());
737 g->SetName(Form("gResidual%sPerChMean_ClusterIn",axes[ia]));
738 g->SetTitle(Form("cluster-track residual-%s per Ch: mean (cluster in);chamber ID;<#Delta_{%s}> (cm)",axes[ia],axes[ia]));
739 g->SetMarkerStyle(kFullDotLarge);
f4d9245c 740 fChamberRes->AddAtAndExpand(g, kResidualPerChMeanClusterIn+ia);
7cbb1928 741 g = new TGraphErrors(AliMUONConstants::NTrackingCh());
742 g->SetName(Form("gResidual%sPerChMean_ClusterOut",axes[ia]));
743 g->SetTitle(Form("cluster-track residual-%s per Ch: mean (cluster out);chamber ID;<#Delta_{%s}> (cm)",axes[ia],axes[ia]));
744 g->SetMarkerStyle(kFullDotLarge);
f4d9245c 745 fChamberRes->AddAtAndExpand(g, kResidualPerChMeanClusterOut+ia);
7cbb1928 746
747 g = new TGraphErrors(2*AliMUONConstants::NTrackingCh());
748 g->SetName(Form("gResidual%sPerHalfChMean_ClusterIn",axes[ia]));
749 g->SetTitle(Form("cluster-track residual-%s per half Ch: mean (cluster in);half chamber ID;<#Delta_{%s}> (cm)",axes[ia],axes[ia]));
750 g->SetMarkerStyle(kFullDotLarge);
f4d9245c 751 fChamberRes->AddAtAndExpand(g, kResidualPerHalfChMeanClusterIn+ia);
7cbb1928 752 g = new TGraphErrors(2*AliMUONConstants::NTrackingCh());
753 g->SetName(Form("gResidual%sPerHalfChMean_ClusterOut",axes[ia]));
754 g->SetTitle(Form("cluster-track residual-%s per half Ch: mean (cluster out);half chamber ID;<#Delta_{%s}> (cm)",axes[ia],axes[ia]));
755 g->SetMarkerStyle(kFullDotLarge);
f4d9245c 756 fChamberRes->AddAtAndExpand(g, kResidualPerHalfChMeanClusterOut+ia);
7cbb1928 757
758 g = new TGraphErrors(fNDE);
759 g->SetName(Form("gResidual%sPerDEMean_ClusterIn",axes[ia]));
760 g->SetTitle(Form("cluster-track residual-%s per DE: mean (cluster in);DE ID;<#Delta_{%s}> (cm)",axes[ia],axes[ia]));
761 g->SetMarkerStyle(kFullDotLarge);
f4d9245c 762 fChamberRes->AddAtAndExpand(g, kResidualPerDEMeanClusterIn+ia);
7cbb1928 763 g = new TGraphErrors(fNDE);
764 g->SetName(Form("gResidual%sPerDEMean_ClusterOut",axes[ia]));
765 g->SetTitle(Form("cluster-track residual-%s per DE: mean (cluster out);DE ID;<#Delta_{%s}> (cm)",axes[ia],axes[ia]));
766 g->SetMarkerStyle(kFullDotLarge);
f4d9245c 767 fChamberRes->AddAtAndExpand(g, kResidualPerDEMeanClusterOut+ia);
7cbb1928 768
769 g = new TGraphErrors(AliMUONConstants::NTrackingCh());
770 g->SetName(Form("gResidual%sPerChSigma_ClusterIn",axes[ia]));
771 g->SetTitle(Form("cluster-track residual-%s per Ch: sigma (cluster in);chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
772 g->SetMarkerStyle(kFullDotLarge);
f4d9245c 773 fChamberRes->AddAtAndExpand(g, kResidualPerChSigmaClusterIn+ia);
7cbb1928 774 g = new TGraphErrors(AliMUONConstants::NTrackingCh());
775 g->SetName(Form("gResidual%sPerChSigma_ClusterOut",axes[ia]));
776 g->SetTitle(Form("cluster-track residual-%s per Ch: sigma (cluster out);chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
777 g->SetMarkerStyle(kFullDotLarge);
f4d9245c 778 fChamberRes->AddAtAndExpand(g, kResidualPerChSigmaClusterOut+ia);
7cbb1928 779
780 g = new TGraphErrors(AliMUONConstants::NTrackingCh());
781 g->SetName(Form("gResidual%sPerChDispersion_ClusterOut",axes[ia]));
782 g->SetTitle(Form("cluster-track residual-%s per Ch: dispersion (cluster out);chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
783 g->SetMarkerStyle(kFullDotLarge);
f4d9245c 784 fChamberRes->AddAtAndExpand(g, kResidualPerChDispersionClusterOut+ia);
7cbb1928 785
786 g = new TGraphErrors(AliMUONConstants::NTrackingCh());
787 g->SetName(Form("gCombinedResidual%sPerChSigma",axes[ia]));
788 g->SetTitle(Form("combined cluster-track residual-%s per Ch: sigma;chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
789 g->SetMarkerStyle(kFullDotLarge);
790 fChamberRes->AddAtAndExpand(g, kCombinedResidualPerChSigma+ia);
791
792 g = new TGraphErrors(2*AliMUONConstants::NTrackingCh());
793 g->SetName(Form("gCombinedResidual%sPerHalfChSigma",axes[ia]));
794 g->SetTitle(Form("combined cluster-track residual-%s per half Ch: sigma;half chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
795 g->SetMarkerStyle(kFullDotLarge);
796 fChamberRes->AddAtAndExpand(g, kCombinedResidualPerHalfChSigma+ia);
797
798 g = new TGraphErrors(fNDE);
799 g->SetName(Form("gCombinedResidual%sPerDESigma",axes[ia]));
800 g->SetTitle(Form("combined cluster-track residual-%s per DE: sigma;DE ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
801 g->SetMarkerStyle(kFullDotLarge);
802 fChamberRes->AddAtAndExpand(g, kCombinedResidualPerDESigma+ia);
803
804 mg = new TMultiGraph(Form("mgCombinedResidual%sSigmaVsP",axes[ia]),Form("cluster %s-resolution per chamber versus momentum;p (GeV/c^{2});#sigma_{%s} (cm)",axes[ia],axes[ia]));
805 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
f4d9245c 806 g = new TGraphErrors(((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsPClusterIn+10*ia+i))->GetNbinsX());
7cbb1928 807 g->SetName(Form("gRes%sVsP_ch%d",axes[ia],i+1));
808 g->SetMarkerStyle(kFullDotMedium);
809 g->SetMarkerColor(i+1+i/9);
810 mg->Add(g,"p");
811 }
812 fChamberRes->AddAtAndExpand(mg, kCombinedResidualSigmaVsP+ia);
813
4060971a 814 g = new TGraphErrors(((TH2F*)fResidualsVsP->UncheckedAt(kResidualVsPClusterIn+ia))->GetNbinsX());
815 g->SetName(Form("gCombinedResidual%sSigmaVsP",axes[ia]));
816 g->SetTitle(Form("cluster %s-resolution integrated over chambers versus momentum: sigma;p (GeV/c^{2});#sigma_{%s} (cm)",axes[ia],axes[ia]));
817 g->SetMarkerStyle(kFullDotLarge);
818 fChamberRes->AddAtAndExpand(g, kCombinedResidualAllChSigmaVsP+ia);
819
820 mg = new TMultiGraph(Form("mgCombinedResidual%sSigmaVsCent",axes[ia]),Form("cluster %s-resolution per chamber versus centrality;centrality (%%);#sigma_{%s} (cm)",axes[ia],axes[ia]));
821 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
822 g = new TGraphErrors(((TH2F*)fResidualsVsCent->UncheckedAt(kResidualInChVsCentClusterIn+10*ia+i))->GetNbinsX());
823 g->SetName(Form("gRes%sVsCent_ch%d",axes[ia],i+1));
824 g->SetMarkerStyle(kFullDotMedium);
825 g->SetMarkerColor(i+1+i/9);
826 mg->Add(g,"p");
827 }
828 fChamberRes->AddAtAndExpand(mg, kCombinedResidualSigmaVsCent+ia);
829
830 g = new TGraphErrors(((TH2F*)fResidualsVsCent->UncheckedAt(kResidualVsCentClusterIn+ia))->GetNbinsX());
831 g->SetName(Form("gCombinedResidual%sSigmaVsCent",axes[ia]));
832 g->SetTitle(Form("cluster %s-resolution integrated over chambers versus centrality: sigma;centrality (%%);#sigma_{%s} (cm)",axes[ia],axes[ia]));
833 g->SetMarkerStyle(kFullDotLarge);
834 fChamberRes->AddAtAndExpand(g, kCombinedResidualAllChSigmaVsCent+ia);
835
7cbb1928 836 g = new TGraphErrors(AliMUONConstants::NTrackingCh());
837 g->SetName(Form("gTrackRes%sPerCh",axes[ia]));
838 g->SetTitle(Form("track <#sigma_{%s}> per Ch;chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
839 g->SetMarkerStyle(kFullDotLarge);
840 fChamberRes->AddAtAndExpand(g, kTrackResPerChMean+ia);
841
842 g = new TGraphErrors(AliMUONConstants::NTrackingCh());
843 g->SetName(Form("gMCS%sPerCh",axes[ia]));
844 g->SetTitle(Form("MCS %s-dispersion of extrapolated track per Ch;chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
845 g->SetMarkerStyle(kFullDotLarge);
846 fChamberRes->AddAtAndExpand(g, kMCSPerChMean+ia);
847
848 g = new TGraphErrors(AliMUONConstants::NTrackingCh());
849 g->SetName(Form("gClusterRes%sPerCh",axes[ia]));
850 g->SetTitle(Form("cluster <#sigma_{%s}> per Ch;chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
851 g->SetMarkerStyle(kFullDotLarge);
852 fChamberRes->AddAtAndExpand(g, kClusterResPerCh+ia);
853
854 g = new TGraphErrors(2*AliMUONConstants::NTrackingCh());
855 g->SetName(Form("gClusterRes%sPerHalfCh",axes[ia]));
856 g->SetTitle(Form("cluster <#sigma_{%s}> per half Ch;half chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
857 g->SetMarkerStyle(kFullDotLarge);
858 fChamberRes->AddAtAndExpand(g, kClusterResPerHalfCh+ia);
859
860 g = new TGraphErrors(fNDE);
861 g->SetName(Form("gClusterRes%sPerDE",axes[ia]));
862 g->SetTitle(Form("cluster <#sigma_{%s}> per DE;DE ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
863 g->SetMarkerStyle(kFullDotLarge);
864 fChamberRes->AddAtAndExpand(g, kClusterResPerDE+ia);
865
866 g = new TGraphErrors(AliMUONConstants::NTrackingCh());
867 g->SetName(Form("gCalcClusterRes%sPerCh",axes[ia]));
868 g->SetTitle(Form("calculated cluster <#sigma_{%s}> per Ch;chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
869 g->SetMarkerStyle(kFullDotLarge);
870 fChamberRes->AddAtAndExpand(g, kCalcClusterResPerCh+ia);
871
872 g = new TGraphErrors(AliMUONConstants::NTrackingCh());
873 g->SetName(Form("gLocalChi2%sPerChMean",axes[ia]));
874 g->SetTitle(Form("local chi2-%s per Ch: mean;chamber ID;<local #chi^{2}_{%s}>",axes[ia],axes[ia]));
875 g->SetMarkerStyle(kFullDotLarge);
876 fLocalChi2->AddAtAndExpand(g, kLocalChi2PerChMean+ia);
877
878 g = new TGraphErrors(fNDE);
879 g->SetName(Form("gLocalChi2%sPerDEMean",axes[ia]));
880 g->SetTitle(Form("local chi2-%s per DE: mean;DE ID;<local #chi^{2}_{%s}>",axes[ia],axes[ia]));
881 g->SetMarkerStyle(kFullDotLarge);
882 fLocalChi2->AddAtAndExpand(g, kLocalChi2PerDEMean+ia);
883
4060971a 884 // compute residual mean and dispersion integrated over chambers versus p
885 FillSigmaClusterVsP((TH2F*)fResidualsVsP->UncheckedAt(kResidualVsPClusterIn+ia),
886 (TH2F*)fResidualsVsP->UncheckedAt(kResidualVsPClusterOut+ia),
887 (TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualAllChSigmaVsP+ia));
888
889 // compute residual mean and dispersion integrated over chambers versus centrality
890 FillSigmaClusterVsCent((TH2F*)fResidualsVsCent->UncheckedAt(kResidualVsCentClusterIn+ia),
891 (TH2F*)fResidualsVsCent->UncheckedAt(kResidualVsCentClusterOut+ia),
892 (TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualAllChSigmaVsCent+ia));
893
7cbb1928 894 // compute residual mean and dispersion and averaged local chi2 per chamber and half chamber
00612ffc 895 Double_t meanIn, meanInErr, meanOut, meanOutErr, sigma, sigmaIn, sigmaInErr, sigmaOut, sigmaOutErr;
7cbb1928 896 Double_t sigmaTrack, sigmaTrackErr, sigmaMCS, sigmaMCSErr, clusterRes, clusterResErr, sigmaCluster, sigmaClusterErr;
897 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
898
899 // method 1
f4d9245c 900 TH1D *tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerChClusterIn+ia))->ProjectionY("tmp",i+1,i+1,"e");
901 GetMean(tmp, meanIn, meanInErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChMeanClusterIn+ia), i, i+1);
902 GetRMS(tmp, sigmaIn, sigmaInErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChSigmaClusterIn+ia), i, i+1);
7cbb1928 903 delete tmp;
904
f4d9245c 905 tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerChClusterOut+ia))->ProjectionY("tmp",i+1,i+1,"e");
906 GetMean(tmp, meanOut, meanOutErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChMeanClusterOut+ia), i, i+1);
907 GetRMS(tmp, sigmaOut, sigmaOutErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChSigmaClusterOut+ia), i, i+1);
7cbb1928 908 delete tmp;
909
910 if (fCorrectForSystematics) {
00612ffc 911 sigma = TMath::Sqrt(sigmaIn*sigmaIn + meanIn*meanIn);
912 sigmaInErr = (sigma>0) ? TMath::Sqrt(sigmaIn*sigmaIn*sigmaInErr*sigmaInErr + meanIn*meanIn*meanInErr*meanInErr) / sigma : 0.;
913 sigmaIn = sigma;
914 sigma = TMath::Sqrt(sigmaOut*sigmaOut + meanOut*meanOut);
915 sigmaOutErr = (sigma>0) ? TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + meanOut*meanOut*meanOutErr*meanOutErr) / sigma : 0.;
916 sigmaOut = sigma;
7cbb1928 917 }
f4d9245c 918 ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChDispersionClusterOut+ia))->SetPoint(i, i+1, sigmaOut);
919 ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChDispersionClusterOut+ia))->SetPointError(i, 0., sigmaOutErr);
7cbb1928 920
921 clusterRes = TMath::Sqrt(sigmaIn*sigmaOut);
922 // clusterResErr = (clusterRes > 0.) ? 0.5 * TMath::Sqrt(sigmaInErr*sigmaInErr*sigmaOut*sigmaOut + sigmaIn*sigmaIn*sigmaOutErr*sigmaOutErr) / clusterRes : 0.;
923 clusterResErr = TMath::Sqrt(sigmaInErr*sigmaOutErr);
924 ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerChSigma+ia))->SetPoint(i, i+1, clusterRes);
925 ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerChSigma+ia))->SetPointError(i, 0., clusterResErr);
926 newClusterRes[ia][i] = clusterRes;
927 newClusterResErr[ia][i] = clusterResErr;
928
929 // method 2
930 tmp = ((TH2F*)fResiduals->UncheckedAt(kTrackResPerCh+ia))->ProjectionY("tmp",i+1,i+1,"e");
00612ffc 931 GetMean(tmp, sigmaTrack, sigmaTrackErr, (TGraphErrors*)fChamberRes->UncheckedAt(kTrackResPerChMean+ia), i, i+1, kFALSE, kFALSE);
7cbb1928 932 delete tmp;
933
934 tmp = ((TH2F*)fResiduals->UncheckedAt(kMCSPerCh+ia))->ProjectionY("tmp",i+1,i+1,"e");
00612ffc 935 GetMean(tmp, sigmaMCS, sigmaMCSErr, (TGraphErrors*)fChamberRes->UncheckedAt(kMCSPerChMean+ia), i, i+1, kFALSE, kFALSE);
7cbb1928 936 delete tmp;
937
938 sigmaCluster = sigmaOut*sigmaOut - sigmaTrack*sigmaTrack;
939 if (sigmaCluster > 0.) {
940 sigmaCluster = TMath::Sqrt(sigmaCluster);
941 sigmaClusterErr = TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + sigmaTrack*sigmaTrack*sigmaTrackErr*sigmaTrackErr) / sigmaCluster;
942 } else {
943 sigmaCluster = 0.;
944 sigmaClusterErr = 0.;
945 }
946 ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerCh+ia))->SetPoint(i, i+1, sigmaCluster);
947 ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerCh+ia))->SetPointError(i, 0., sigmaClusterErr);
948
949 // method 3
950 tmp = ((TH2F*)fResiduals->UncheckedAt(kClusterRes2PerCh+ia))->ProjectionY("tmp",i+1,i+1,"e");
951 ZoomRight(tmp);
952 clusterRes = tmp->GetMean();
953 if (clusterRes > 0.) {
954 ((TGraphErrors*)fChamberRes->UncheckedAt(kCalcClusterResPerCh+ia))->SetPoint(i, i+1, TMath::Sqrt(clusterRes));
955 ((TGraphErrors*)fChamberRes->UncheckedAt(kCalcClusterResPerCh+ia))->SetPointError(i, 0., 0.5 * tmp->GetMeanError() / TMath::Sqrt(clusterRes));
956 } else {
957 ((TGraphErrors*)fChamberRes->UncheckedAt(kCalcClusterResPerCh+ia))->SetPoint(i, i+1, 0.);
958 ((TGraphErrors*)fChamberRes->UncheckedAt(kCalcClusterResPerCh+ia))->SetPointError(i, 0., 0.);
959 }
960 delete tmp;
961
962 // method 1 versus p
f4d9245c 963 FillSigmaClusterVsP((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsPClusterIn+10*ia+i),
964 (TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsPClusterOut+10*ia+i),
7cbb1928 965 (TGraphErrors*)((TMultiGraph*)fChamberRes->UncheckedAt(kCombinedResidualSigmaVsP+ia))->GetListOfGraphs()->FindObject(Form("gRes%sVsP_ch%d",axes[ia],i+1)));
966
4060971a 967 // method 1 versus centrality
968 FillSigmaClusterVsCent((TH2F*)fResidualsVsCent->UncheckedAt(kResidualInChVsCentClusterIn+10*ia+i),
969 (TH2F*)fResidualsVsCent->UncheckedAt(kResidualInChVsCentClusterOut+10*ia+i),
970 (TGraphErrors*)((TMultiGraph*)fChamberRes->UncheckedAt(kCombinedResidualSigmaVsCent+ia))->GetListOfGraphs()->FindObject(Form("gRes%sVsCent_ch%d",axes[ia],i+1)));
971
7cbb1928 972 // compute residual mean and dispersion per half chamber
973 for (Int_t j = 0; j < 2; j++) {
974 Int_t k = 2*i+j;
975
976 // method 1
f4d9245c 977 tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfChClusterIn+ia))->ProjectionY("tmp",k+1,k+1,"e");
978 GetMean(tmp, meanIn, meanInErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMeanClusterIn+ia), k, k+1);
7cbb1928 979 GetRMS(tmp, sigmaIn, sigmaInErr);
980 delete tmp;
981
f4d9245c 982 tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfChClusterOut+ia))->ProjectionY("tmp",k+1,k+1,"e");
983 GetMean(tmp, meanOut, meanOutErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMeanClusterOut+ia), k, k+1);
7cbb1928 984 GetRMS(tmp, sigmaOut, sigmaOutErr);
985 delete tmp;
986
987 if (fCorrectForSystematics) {
00612ffc 988 sigma = TMath::Sqrt(sigmaIn*sigmaIn + meanIn*meanIn);
989 sigmaInErr = (sigma>0) ? TMath::Sqrt(sigmaIn*sigmaIn*sigmaInErr*sigmaInErr + meanIn*meanIn*meanInErr*meanInErr) / sigma : 0.;
990 sigmaIn = sigma;
991 sigma = TMath::Sqrt(sigmaOut*sigmaOut + meanOut*meanOut);
992 sigmaOutErr = (sigma>0) ? TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + meanOut*meanOut*meanOutErr*meanOutErr) / sigma : 0.;
993 sigmaOut = sigma;
7cbb1928 994 }
995
996 clusterRes = TMath::Sqrt(sigmaIn*sigmaOut);
997 // clusterResErr = (clusterRes > 0.) ? 0.5 * TMath::Sqrt(sigmaInErr*sigmaInErr*sigmaOut*sigmaOut + sigmaIn*sigmaIn*sigmaOutErr*sigmaOutErr) / clusterRes : 0.;
998 clusterResErr = TMath::Sqrt(sigmaInErr*sigmaOutErr);
999 ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerHalfChSigma+ia))->SetPoint(k, k+1, clusterRes);
1000 ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerHalfChSigma+ia))->SetPointError(k, 0., clusterResErr);
1001
1002 // method 2
1003 tmp = ((TH2F*)fResiduals->UncheckedAt(kTrackResPerHalfCh+ia))->ProjectionY("tmp",k+1,k+1,"e");
00612ffc 1004 GetMean(tmp, sigmaTrack, sigmaTrackErr, 0x0, 0, 0, kFALSE, kFALSE);
7cbb1928 1005 delete tmp;
1006
1007 tmp = ((TH2F*)fResiduals->UncheckedAt(kMCSPerHalfCh+ia))->ProjectionY("tmp",k+1,k+1,"e");
00612ffc 1008 GetMean(tmp, sigmaMCS, sigmaMCSErr, 0x0, 0, 0, kFALSE, kFALSE);
7cbb1928 1009 delete tmp;
1010
1011 sigmaCluster = sigmaOut*sigmaOut - sigmaTrack*sigmaTrack;
1012 if (sigmaCluster > 0.) {
1013 sigmaCluster = TMath::Sqrt(sigmaCluster);
1014 sigmaClusterErr = TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + sigmaTrack*sigmaTrack*sigmaTrackErr*sigmaTrackErr) / sigmaCluster;
1015 } else {
1016 sigmaCluster = 0.;
1017 sigmaClusterErr = 0.;
1018 }
1019 ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerHalfCh+ia))->SetPoint(k, k+1, sigmaCluster);
1020 ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerHalfCh+ia))->SetPointError(k, 0., sigmaClusterErr);
1021
1022 }
1023
1024 // compute averaged local chi2
1025 tmp = ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerCh+ia))->ProjectionY("tmp",i+1,i+1,"e");
1026 ((TGraphErrors*)fLocalChi2->UncheckedAt(kLocalChi2PerChMean+ia))->SetPoint(i, i+1, tmp->GetMean());
1027 ((TGraphErrors*)fLocalChi2->UncheckedAt(kLocalChi2PerChMean+ia))->SetPointError(i, 0., tmp->GetMeanError());
1028 delete tmp;
1029
1030 }
1031
1032 // compute residual mean and dispersion per DE
1033 for (Int_t i = 0; i < fNDE; i++) {
1034
1035 // method 1
f4d9245c 1036 TH1D *tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerDEClusterIn+ia))->ProjectionY("tmp",i+1,i+1,"e");
1037 GetMean(tmp, meanIn, meanInErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMeanClusterIn+ia), i, i+1);
7cbb1928 1038 GetRMS(tmp, sigmaIn, sigmaInErr);
1039 delete tmp;
1040
f4d9245c 1041 tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerDEClusterOut+ia))->ProjectionY("tmp",i+1,i+1,"e");
1042 GetMean(tmp, meanOut, meanOutErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMeanClusterOut+ia), i, i+1);
7cbb1928 1043 GetRMS(tmp, sigmaOut, sigmaOutErr);
1044 delete tmp;
1045
1046 if (fCorrectForSystematics) {
00612ffc 1047 sigma = TMath::Sqrt(sigmaIn*sigmaIn + meanIn*meanIn);
1048 sigmaInErr = (sigma>0) ? TMath::Sqrt(sigmaIn*sigmaIn*sigmaInErr*sigmaInErr + meanIn*meanIn*meanInErr*meanInErr) / sigma : 0.;
1049 sigmaIn = sigma;
1050 sigma = TMath::Sqrt(sigmaOut*sigmaOut + meanOut*meanOut);
1051 sigmaOutErr = (sigma>0) ? TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + meanOut*meanOut*meanOutErr*meanOutErr) / sigma : 0.;
1052 sigmaOut = sigma;
7cbb1928 1053 }
1054
1055 clusterRes = TMath::Sqrt(sigmaIn*sigmaOut);
1056 // clusterResErr = (clusterRes > 0.) ? 0.5 * TMath::Sqrt(sigmaInErr*sigmaInErr*sigmaOut*sigmaOut + sigmaIn*sigmaIn*sigmaOutErr*sigmaOutErr) / clusterRes : 0.;
1057 clusterResErr = TMath::Sqrt(sigmaInErr*sigmaOutErr);
1058 ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerDESigma+ia))->SetPoint(i, i+1, clusterRes);
1059 ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerDESigma+ia))->SetPointError(i, 0., clusterResErr);
1060
1061 // method 2
1062 tmp = ((TH2F*)fResiduals->UncheckedAt(kTrackResPerDE+ia))->ProjectionY("tmp",i+1,i+1,"e");
00612ffc 1063 GetMean(tmp, sigmaTrack, sigmaTrackErr, 0x0, 0, 0, kFALSE, kFALSE);
7cbb1928 1064 delete tmp;
1065
1066 tmp = ((TH2F*)fResiduals->UncheckedAt(kMCSPerDE+ia))->ProjectionY("tmp",i+1,i+1,"e");
00612ffc 1067 GetMean(tmp, sigmaMCS, sigmaMCSErr, 0x0, 0, 0, kFALSE, kFALSE);
7cbb1928 1068 delete tmp;
1069
1070 sigmaCluster = sigmaOut*sigmaOut - sigmaTrack*sigmaTrack;
1071 if (sigmaCluster > 0.) {
1072 sigmaCluster = TMath::Sqrt(sigmaCluster);
1073 sigmaClusterErr = TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + sigmaTrack*sigmaTrack*sigmaTrackErr*sigmaTrackErr) / sigmaCluster;
1074 } else {
1075 sigmaCluster = 0.;
1076 sigmaClusterErr = 0.;
1077 }
1078 ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerDE+ia))->SetPoint(i, i+1, sigmaCluster);
1079 ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerDE+ia))->SetPointError(i, 0., sigmaClusterErr);
1080
1081 // compute averaged local chi2
1082 tmp = ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerDE+ia))->ProjectionY("tmp",i+1,i+1,"e");
1083 ((TGraphErrors*)fLocalChi2->UncheckedAt(kLocalChi2PerDEMean+ia))->SetPoint(i, i+1, tmp->GetMean());
1084 ((TGraphErrors*)fLocalChi2->UncheckedAt(kLocalChi2PerDEMean+ia))->SetPointError(i, 0., tmp->GetMeanError());
1085 delete tmp;
1086
1087 }
1088
12d09c5b 1089 // set half-chamber graph labels
1090 TAxis* xAxis = ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfChClusterIn+ia))->GetXaxis();
1091 ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMeanClusterIn+ia))->GetXaxis()->Set(20, 0.5, 20.5);
1092 ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMeanClusterOut+ia))->GetXaxis()->Set(20, 0.5, 20.5);
1093 ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerHalfChSigma+ia))->GetXaxis()->Set(20, 0.5, 20.5);
1094 ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerHalfCh+ia))->GetXaxis()->Set(20, 0.5, 20.5);
1095 for (Int_t i = 1; i <= 20; i++) {
1096 const char* label = xAxis->GetBinLabel(i);
1097 ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMeanClusterIn+ia))->GetXaxis()->SetBinLabel(i, label);
1098 ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMeanClusterOut+ia))->GetXaxis()->SetBinLabel(i, label);
1099 ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerHalfChSigma+ia))->GetXaxis()->SetBinLabel(i, label);
1100 ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerHalfCh+ia))->GetXaxis()->SetBinLabel(i, label);
1101 }
1102
1103 // set DE graph labels
1104 xAxis = ((TH2F*)fResiduals->UncheckedAt(kResidualPerDEClusterOut+ia))->GetXaxis();
f4d9245c 1105 ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMeanClusterIn+ia))->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
1106 ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMeanClusterOut+ia))->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
7cbb1928 1107 ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerDESigma+ia))->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
1108 ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerDE+ia))->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
1109 ((TGraphErrors*)fLocalChi2->UncheckedAt(kLocalChi2PerDEMean+ia))->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
1110 for (Int_t i = 1; i <= fNDE; i++) {
1111 const char* label = xAxis->GetBinLabel(i);
f4d9245c 1112 ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMeanClusterIn+ia))->GetXaxis()->SetBinLabel(i, label);
1113 ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMeanClusterOut+ia))->GetXaxis()->SetBinLabel(i, label);
7cbb1928 1114 ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerDESigma+ia))->GetXaxis()->SetBinLabel(i, label);
1115 ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerDE+ia))->GetXaxis()->SetBinLabel(i, label);
1116 ((TGraphErrors*)fLocalChi2->UncheckedAt(kLocalChi2PerDEMean+ia))->GetXaxis()->SetBinLabel(i, label);
1117 }
1118
1119 }
1120
1121 // compute averaged local chi2 per chamber (X+Y)
1122 TH2F* h2 = (TH2F*)fResiduals->UncheckedAt(kLocalChi2PerCh+2);
1123 g = new TGraphErrors(AliMUONConstants::NTrackingCh());
1124 g->SetName("gLocalChi2PerChMean");
1125 g->SetTitle("local chi2 per Ch: mean;chamber ID;<local #chi^{2}>");
1126 g->SetMarkerStyle(kFullDotLarge);
1127 fLocalChi2->AddAtAndExpand(g, kLocalChi2PerChMean+2);
1128 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
1129 TH1D* tmp = h2->ProjectionY("tmp",i+1,i+1,"e");
1130 g->SetPoint(i, i+1, tmp->GetMean());
1131 g->SetPointError(i, 0., tmp->GetMeanError());
1132 delete tmp;
1133 }
1134
1135 // compute averaged local chi2 per DE (X+Y)
1136 h2 = (TH2F*)fResiduals->UncheckedAt(kLocalChi2PerDE+2);
1137 g = new TGraphErrors(fNDE);
1138 g->SetName("gLocalChi2PerDEMean");
1139 g->SetTitle("local chi2 per DE: mean;DE ID;<local #chi^{2}>");
1140 g->SetMarkerStyle(kFullDotLarge);
1141 fLocalChi2->AddAtAndExpand(g, kLocalChi2PerDEMean+2);
1142 for (Int_t i = 0; i < fNDE; i++) {
1143 TH1D* tmp = h2->ProjectionY("tmp",i+1,i+1,"e");
1144 g->SetPoint(i, i+1, tmp->GetMean());
1145 g->SetPointError(i, 0., tmp->GetMeanError());
1146 delete tmp;
1147 }
1148
1149 // set graph labels
1150 g->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
1151 for (Int_t i = 1; i <= fNDE; i++) {
1152 const char* label = h2->GetXaxis()->GetBinLabel(i);
1153 g->GetXaxis()->SetBinLabel(i, label);
1154 }
1155
1156 // display
1157 fCanvases = new TObjArray(1000);
1158 fCanvases->SetOwner();
1159
1160 TLegend *lResPerChMean = new TLegend(0.75,0.85,0.99,0.99);
1161 TLegend *lResPerChSigma1 = new TLegend(0.75,0.85,0.99,0.99);
1162 TLegend *lResPerChSigma2 = new TLegend(0.75,0.85,0.99,0.99);
1163 TLegend *lResPerChSigma3 = new TLegend(0.75,0.85,0.99,0.99);
1164
1165 TCanvas* cResPerCh = new TCanvas("cResPerCh","cResPerCh",1200,500);
1166 cResPerCh->Divide(4,2);
1167 for (Int_t ia = 0; ia < 2; ia++) {
1168 cResPerCh->cd(1+4*ia);
f4d9245c 1169 g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChMeanClusterOut+ia);
7cbb1928 1170 g->Draw("ap");
1171 g->SetMarkerColor(2);
1172 g->SetLineColor(2);
1173 if (ia == 0) lResPerChMean->AddEntry(g,"cluster out","PL");
f4d9245c 1174 g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChMeanClusterIn+ia);
7cbb1928 1175 g->Draw("p");
1176 g->SetMarkerColor(4);
1177 g->SetLineColor(4);
1178 if (ia == 0) lResPerChMean->AddEntry(g,"cluster in","PL");
1179 if (ia == 0) lResPerChMean->Draw();
1180 else lResPerChMean->DrawClone();
1181 cResPerCh->cd(2+4*ia);
f4d9245c 1182 g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChSigmaClusterOut+ia);
7cbb1928 1183 g->Draw("ap");
1184 g->SetMinimum(0.);
1185 g->SetMarkerColor(2);
1186 g->SetLineColor(2);
1187 if (ia == 0) lResPerChSigma1->AddEntry(g,"cluster out","PL");
f4d9245c 1188 g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChSigmaClusterIn+ia);
7cbb1928 1189 g->Draw("p");
1190 g->SetMarkerColor(4);
1191 g->SetLineColor(4);
1192 if (ia == 0) lResPerChSigma1->AddEntry(g,"cluster in","PL");
1193 g = (TGraphErrors*)fChamberRes->UncheckedAt(kMCSPerChMean+ia);
1194 g->Draw("p");
1195 g->SetMarkerColor(5);
1196 g->SetLineColor(5);
1197 if (ia == 0) lResPerChSigma1->AddEntry(g,"MCS","PL");
1198 g = (TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerChSigma+ia);
1199 g->Draw("p");
1200 g->SetMarkerColor(3);
1201 g->SetLineColor(3);
1202 if (ia == 0) lResPerChSigma1->AddEntry(g,"combined 1","PL");
1203 if (ia == 0) lResPerChSigma1->Draw();
1204 else lResPerChSigma1->DrawClone();
1205 cResPerCh->cd(3+4*ia);
f4d9245c 1206 g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChDispersionClusterOut+ia);
7cbb1928 1207 g->Draw("ap");
1208 g->SetMinimum(0.);
1209 g->SetMarkerColor(2);
1210 g->SetLineColor(2);
1211 if (ia == 0) lResPerChSigma2->AddEntry(g,"cluster out","PL");
1212 g = (TGraphErrors*)fChamberRes->UncheckedAt(kMCSPerChMean+ia);
1213 g->Draw("p");
1214 if (ia == 0) lResPerChSigma2->AddEntry(g,"MCS","PL");
1215 g = (TGraphErrors*)fChamberRes->UncheckedAt(kTrackResPerChMean+ia);
1216 g->Draw("p");
1217 g->SetMarkerColor(4);
1218 g->SetLineColor(4);
1219 if (ia == 0) lResPerChSigma2->AddEntry(g,"track res.","PL");
1220 g = (TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerCh+ia);
1221 g->Draw("p");
1222 if (ia == 0) lResPerChSigma2->AddEntry(g,"combined 2","PL");
1223 if (ia == 0) lResPerChSigma2->Draw();
1224 else lResPerChSigma2->DrawClone();
1225 cResPerCh->cd(4+4*ia);
1226 g = (TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerChSigma+ia);
1227 g->Draw("ap");
1228 g->SetMinimum(0.);
1229 if (ia == 0) lResPerChSigma3->AddEntry(g,"combined 1","PL");
1230 g = (TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerCh+ia);
1231 g->Draw("p");
1232 if (ia == 0) lResPerChSigma3->AddEntry(g,"combined 2","PL");
1233 if (ia == 0) lResPerChSigma3->Draw();
1234 else lResPerChSigma3->DrawClone();
1235 }
1236 fCanvases->AddAtAndExpand(cResPerCh, kResPerCh);
1237
1238 TCanvas* cResPerHalfCh = new TCanvas("cResPerHalfCh","cResPerHalfCh",1200,500);
1239 cResPerHalfCh->Divide(2,2);
1240 for (Int_t ia = 0; ia < 2; ia++) {
1241 cResPerHalfCh->cd(1+2*ia);
f4d9245c 1242 g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMeanClusterOut+ia);
7cbb1928 1243 g->Draw("ap");
1244 g->SetMarkerColor(2);
1245 g->SetLineColor(2);
f4d9245c 1246 g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMeanClusterIn+ia);
7cbb1928 1247 g->Draw("p");
1248 g->SetMarkerColor(4);
1249 g->SetLineColor(4);
1250 lResPerChMean->DrawClone();
1251 cResPerHalfCh->cd(2+2*ia);
1252 g = (TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerHalfChSigma+ia);
1253 g->Draw("ap");
1254 g->SetMinimum(0.);
1255 g->SetMarkerColor(3);
1256 g->SetLineColor(3);
1257 g = (TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerHalfCh+ia);
1258 g->Draw("p");
1259 lResPerChSigma3->DrawClone();
1260 }
1261 fCanvases->AddAtAndExpand(cResPerHalfCh, kResPerHalfCh);
1262
1263 TCanvas* cResPerDE = new TCanvas("cResPerDE","cResPerDE",1200,800);
1264 cResPerDE->Divide(1,4);
1265 for (Int_t ia = 0; ia < 2; ia++) {
1266 cResPerDE->cd(1+ia);
f4d9245c 1267 g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMeanClusterOut+ia);
7cbb1928 1268 g->Draw("ap");
1269 g->SetMarkerColor(2);
1270 g->SetLineColor(2);
f4d9245c 1271 g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMeanClusterIn+ia);
7cbb1928 1272 g->Draw("p");
1273 g->SetMarkerColor(4);
1274 g->SetLineColor(4);
1275 lResPerChMean->DrawClone();
1276 cResPerDE->cd(3+ia);
1277 g = (TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerDESigma+ia);
1278 g->Draw("ap");
1279 g->SetMinimum(0.);
1280 g->SetMarkerColor(3);
1281 g->SetLineColor(3);
1282 g = (TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerDE+ia);
1283 g->Draw("p");
1284 lResPerChSigma3->DrawClone();
1285 }
1286 fCanvases->AddAtAndExpand(cResPerDE, kResPerDE);
1287
1288 TCanvas* cResPerChVsP = new TCanvas("cResPerChVsP","cResPerChVsP");
1289 cResPerChVsP->Divide(1,2);
1290 for (Int_t ia = 0; ia < 2; ia++) {
1291 cResPerChVsP->cd(1+ia);
1292 mg = (TMultiGraph*)fChamberRes->UncheckedAt(kCombinedResidualSigmaVsP+ia);
1293 mg->Draw("ap");
1294 }
1295 fCanvases->AddAtAndExpand(cResPerChVsP, kResPerChVsP);
1296
4060971a 1297 TCanvas* cResPerChVsCent = new TCanvas("cResPerChVsCent","cResPerChVsCent");
1298 cResPerChVsCent->Divide(1,2);
1299 for (Int_t ia = 0; ia < 2; ia++) {
1300 cResPerChVsCent->cd(1+ia);
1301 mg = (TMultiGraph*)fChamberRes->UncheckedAt(kCombinedResidualSigmaVsCent+ia);
1302 mg->Draw("ap");
1303 }
1304 fCanvases->AddAtAndExpand(cResPerChVsCent, kResPerChVsCent);
1305
7cbb1928 1306 // print results
00612ffc 1307 if (fPrintClResPerCh) {
1308 printf("\nchamber resolution:\n");
1309 printf(" - non-bending:");
1310 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %5.3f":", %5.3f",newClusterRes[0][i]);
1311 printf("\n - bending:");
1312 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %6.4f":", %6.4f",newClusterRes[1][i]);
1313 printf("\n\n");
1314 }
1315
1316 if (fPrintClResPerDE) {
1317 Double_t iDE, clRes;
1318 printf("\nDE resolution:\n");
1319 printf(" - non-bending:");
1320 for (Int_t i = 0; i < fNDE; i++) {
1321 ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerDESigma))->GetPoint(i, iDE, clRes);
1322 printf((i==0)?" %5.3f":", %5.3f", clRes);
1323 }
1324 printf("\n - bending:");
1325 for (Int_t i = 0; i < fNDE; i++) {
1326 ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerDESigma+1))->GetPoint(i, iDE, clRes);
1327 printf((i==0)?" %6.4f":", %6.4f", clRes);
1328 }
1329 printf("\n\n");
1330 }
7cbb1928 1331
1332 // Post final data.
4060971a 1333 PostData(4, fLocalChi2);
1334 PostData(5, fChamberRes);
7cbb1928 1335}
1336
1337//________________________________________________________________________
1338void AliAnalysisTaskMuonResolution::ModifyClusters(AliMUONTrack& track)
1339{
1340 /// Reset the clusters resolution from the ones given to the task and change
1341 /// the cluster position according to the new alignment parameters if required
1342
1343 Double_t gX,gY,gZ,lX,lY,lZ;
1344
1345 // loop over clusters
1346 Int_t nClusters = track.GetNClusters();
1347 for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
1348
1349 AliMUONVCluster* cl = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->UncheckedAt(iCluster))->GetClusterPtr();
1350
1351 // change their resolution
1352 cl->SetErrXY(fClusterResNB[cl->GetChamberId()], fClusterResB[cl->GetChamberId()]);
1353
1354 // change their position
1355 if (fReAlign) {
1356 gX = cl->GetX();
1357 gY = cl->GetY();
1358 gZ = cl->GetZ();
1359 fOldGeoTransformer->Global2Local(cl->GetDetElemId(),gX,gY,gZ,lX,lY,lZ);
1360 fNewGeoTransformer->Local2Global(cl->GetDetElemId(),lX,lY,lZ,gX,gY,gZ);
1361 cl->SetXYZ(gX,gY,gZ);
1362 }
1363
1364 }
1365
1366}
1367
1368//________________________________________________________________________
1369void AliAnalysisTaskMuonResolution::Zoom(TH1* h, Double_t fractionCut)
1370{
1371 /// Reduce the range of the histogram by removing a given fration of the statistic at each edge
1372 ZoomLeft(h, fractionCut);
1373 ZoomRight(h, fractionCut);
1374}
1375
1376//________________________________________________________________________
1377void AliAnalysisTaskMuonResolution::ZoomLeft(TH1* h, Double_t fractionCut)
1378{
1379 /// Reduce the range of the histogram by removing a given fration of the statistic on the left side
1380 Int_t maxEventsCut = (Int_t) (fractionCut * h->GetEntries());
1381 Int_t nBins = h->GetNbinsX();
1382
1383 // set low edge
1384 Int_t minBin;
1385 Int_t eventsCut = 0;
1386 for (minBin = 1; minBin <= nBins; minBin++) {
1387 eventsCut += (Int_t) h->GetBinContent(minBin);
1388 if (eventsCut > maxEventsCut) break;
1389 }
1390
1391 // set new axis range
1392 h->GetXaxis()->SetRange(--minBin, h->GetXaxis()->GetLast());
1393}
1394
1395//________________________________________________________________________
1396void AliAnalysisTaskMuonResolution::ZoomRight(TH1* h, Double_t fractionCut)
1397{
1398 /// Reduce the range of the histogram by removing a given fration of the statistic on the right side
1399 Int_t maxEventsCut = (Int_t) (fractionCut * h->GetEntries());
1400 Int_t nBins = h->GetNbinsX();
1401
1402 // set high edge
1403 Int_t maxBin;
1404 Int_t eventsCut = 0;
1405 for (maxBin = nBins; maxBin >= 1; maxBin--) {
1406 eventsCut += (Int_t) h->GetBinContent(maxBin);
1407 if (eventsCut > maxEventsCut) break;
1408 }
1409
1410 // set new axis range
1411 h->GetXaxis()->SetRange(h->GetXaxis()->GetFirst(), ++maxBin);
1412}
1413
1414//________________________________________________________________________
00612ffc 1415void AliAnalysisTaskMuonResolution::GetMean(TH1* h, Double_t& mean, Double_t& meanErr, TGraphErrors* g, Int_t i, Double_t x, Bool_t zoom, Bool_t enableFit)
7cbb1928 1416{
00612ffc 1417 /// Fill graph with the mean value and the corresponding error (zooming if required)
1418
1419 if (h->GetEntries() < fgkMinEntries) { // not enough entries
1420
1421 mean = 0.;
1422 meanErr = 0.;
1423
1424 } else if (enableFit && fGaus) { // take the mean of a gaussian fit
1425
1426 fGaus->SetParameters(h->GetEntries(), 0., 0.1);
1427
4060971a 1428 if (h->GetUniqueID() != 10) {
1429
1430 // first fit
1431 h->Fit("fGaus", "WWNQ");
1432
1433 // rebin histo
1434 Int_t rebin = TMath::Max(static_cast<Int_t>(0.2*fGaus->GetParameter(2)/h->GetBinWidth(1)),1);
1435 while (h->GetNbinsX()%rebin!=0) rebin--;
1436 h->Rebin(rebin);
1437
1438 // use the unique ID to remember that this histogram has already been rebinned
1439 h->SetUniqueID(10);
1440
1441 }
1442
1443 // second fit
1444 h->Fit("fGaus","NQ");
00612ffc 1445
1446 mean = fGaus->GetParameter(1);
1447 meanErr = fGaus->GetParError(1);
1448
1449 } else { // take the mean of the distribution
1450
00612ffc 1451 if (zoom) Zoom(h);
1452
1453 mean = h->GetMean();
1454 meanErr = h->GetMeanError();
1455
4060971a 1456 if (zoom) h->GetXaxis()->SetRange(0,0);
00612ffc 1457
1458 }
1459
1460 // fill graph if required
7cbb1928 1461 if (g) {
1462 g->SetPoint(i, x, mean);
1463 g->SetPointError(i, 0., meanErr);
1464 }
00612ffc 1465
7cbb1928 1466}
1467
1468//________________________________________________________________________
1469void AliAnalysisTaskMuonResolution::GetRMS(TH1* h, Double_t& rms, Double_t& rmsErr, TGraphErrors* g, Int_t i, Double_t x, Bool_t zoom)
1470{
00612ffc 1471 /// Return the dispersion value and the corresponding error (zooming if required) and fill graph if !=0x0
1472
1473 if (h->GetEntries() < fgkMinEntries) { // not enough entries
1474
1475 rms = 0.;
1476 rmsErr = 0.;
1477
1478 } else if (fGaus) { // take the sigma of a gaussian fit
1479
1480 fGaus->SetParameters(h->GetEntries(), 0., 0.1);
1481
4060971a 1482 if (h->GetUniqueID() != 10) {
1483
1484 // first fit
1485 h->Fit("fGaus", "WWNQ");
1486
1487 // rebin histo
1488 Int_t rebin = TMath::Max(static_cast<Int_t>(0.2*fGaus->GetParameter(2)/h->GetBinWidth(1)),1);
1489 while (h->GetNbinsX()%rebin!=0) rebin--;
1490 h->Rebin(rebin);
1491
1492 // use the unique ID to remember that this histogram has already been rebinned
1493 h->SetUniqueID(10);
1494
1495 }
1496
1497 // second fit
1498 h->Fit("fGaus","NQ");
00612ffc 1499
1500 rms = fGaus->GetParameter(2);
1501 rmsErr = fGaus->GetParError(2);
1502
1503 } else { // take the RMS of the distribution
1504
00612ffc 1505 if (zoom) Zoom(h);
1506
1507 rms = h->GetRMS();
1508 rmsErr = h->GetRMSError();
1509
4060971a 1510 if (zoom) h->GetXaxis()->SetRange(0,0);
00612ffc 1511
1512 }
1513
1514 // fill graph if required
7cbb1928 1515 if (g) {
1516 g->SetPoint(i, x, rms);
1517 g->SetPointError(i, 0., rmsErr);
1518 }
00612ffc 1519
7cbb1928 1520}
1521
1522//________________________________________________________________________
f4d9245c 1523void AliAnalysisTaskMuonResolution::FillSigmaClusterVsP(const TH2* hIn, const TH2* hOut, TGraphErrors* g, Bool_t zoom)
7cbb1928 1524{
1525 /// Fill graph with cluster resolution from combined residuals with cluster in/out (zooming if required)
1526 Double_t sigmaIn, sigmaInErr, sigmaOut, sigmaOutErr, clusterRes, clusterResErr;
1527 for (Int_t j = 1; j <= hIn->GetNbinsX(); j++) {
1528 TH1D* tmp = hIn->ProjectionY("tmp",j,j,"e");
1529 GetRMS(tmp, sigmaIn, sigmaInErr, 0x0, 0, 0., zoom);
1530 delete tmp;
1531 tmp = hOut->ProjectionY("tmp",j,j,"e");
1532 GetRMS(tmp, sigmaOut, sigmaOutErr, 0x0, 0, 0., zoom);
1533 delete tmp;
1534 Double_t p = 0.5 * (hIn->GetBinLowEdge(j) + hIn->GetBinLowEdge(j+1));
1535 Double_t pErr = p - hIn->GetBinLowEdge(j);
1536 clusterRes = TMath::Sqrt(sigmaIn*sigmaOut);
00612ffc 1537 //clusterResErr = (clusterRes > 0.) ? 0.5 * TMath::Sqrt(sigmaInErr*sigmaInErr*sigmaOut*sigmaOut + sigmaIn*sigmaIn*sigmaOutErr*sigmaOutErr) / clusterRes : 0.;
1538 clusterResErr = TMath::Sqrt(sigmaInErr*sigmaOutErr);
7cbb1928 1539 g->SetPoint(j, p, clusterRes);
1540 g->SetPointError(j, pErr, clusterResErr);
1541 }
1542}
1543
4060971a 1544//________________________________________________________________________
1545void AliAnalysisTaskMuonResolution::FillSigmaClusterVsCent(const TH2* hIn, const TH2* hOut, TGraphErrors* g, Bool_t zoom)
1546{
1547 /// Fill graph with cluster resolution from combined residuals with cluster in/out (zooming if required)
1548 Double_t sigmaIn, sigmaInErr, sigmaOut, sigmaOutErr, clusterRes, clusterResErr;
1549 for (Int_t j = 1; j <= hIn->GetNbinsX(); j++) {
1550 TH1D* tmp = hIn->ProjectionY("tmp",j,j,"e");
1551 GetRMS(tmp, sigmaIn, sigmaInErr, 0x0, 0, 0., zoom);
1552 delete tmp;
1553 tmp = hOut->ProjectionY("tmp",j,j,"e");
1554 GetRMS(tmp, sigmaOut, sigmaOutErr, 0x0, 0, 0., zoom);
1555 delete tmp;
1556 Double_t cent = 0.5 * (hIn->GetBinLowEdge(j) + hIn->GetBinLowEdge(j+1));
1557 Double_t centErr = cent - hIn->GetBinLowEdge(j);
1558 clusterRes = TMath::Sqrt(sigmaIn*sigmaOut);
1559 //clusterResErr = (clusterRes > 0.) ? 0.5 * TMath::Sqrt(sigmaInErr*sigmaInErr*sigmaOut*sigmaOut + sigmaIn*sigmaIn*sigmaOutErr*sigmaOutErr) / clusterRes : 0.;
1560 clusterResErr = TMath::Sqrt(sigmaInErr*sigmaOutErr);
1561 g->SetPoint(j, cent, clusterRes);
1562 g->SetPointError(j, centErr, clusterResErr);
1563 }
1564}
1565
7cbb1928 1566//__________________________________________________________________________
1567void AliAnalysisTaskMuonResolution::Cov2CovP(const AliMUONTrackParam &param, TMatrixD &covP)
1568{
1569 /// change coordinate system: (X, SlopeX, Y, SlopeY, q/Pyz) -> (X, Y, pX, pY, pZ)
1570 /// parameters (param) are given in the (X, SlopeX, Y, SlopeY, q/Pyz) coordinate system
1571
1572 // Get useful parameters
1573 Double_t slopeX = param.GetNonBendingSlope();
1574 Double_t slopeY = param.GetBendingSlope();
1575 Double_t qOverPYZ = param.GetInverseBendingMomentum();
1576 Double_t pZ = param.Pz();
1577
1578 // compute Jacobian to change the coordinate system from (X,SlopeX,Y,SlopeY,c/pYZ) to (X,Y,pX,pY,pZ)
1579 Double_t dpZdSlopeY = - qOverPYZ * qOverPYZ * pZ * pZ * pZ * slopeY;
1580 Double_t dpZdQOverPYZ = (qOverPYZ != 0.) ? - pZ / qOverPYZ : - FLT_MAX;
1581 TMatrixD jacob(5,5);
1582 jacob.Zero();
1583 jacob(0,0) = 1.;
1584 jacob(1,2) = 1.;
1585 jacob(2,1) = pZ;
1586 jacob(2,3) = slopeX * dpZdSlopeY;
1587 jacob(2,4) = slopeX * dpZdQOverPYZ;
1588 jacob(3,3) = pZ + slopeY * dpZdSlopeY;
1589 jacob(3,4) = slopeY * dpZdQOverPYZ;
1590 jacob(4,3) = dpZdSlopeY;
1591 jacob(4,4) = dpZdQOverPYZ;
1592
1593 // compute covariances in new coordinate system
1594 TMatrixD tmp(param.GetCovariances(),TMatrixD::kMultTranspose,jacob);
1595 covP.Mult(jacob,tmp);
1596}
1597
00612ffc 1598//__________________________________________________________________________
f4d9245c 1599UInt_t AliAnalysisTaskMuonResolution::BuildTriggerWord(const TString& firedTriggerClasses)
00612ffc 1600{
1601 /// build the trigger word from the fired trigger classes and the list of selectable trigger
1602
1603 UInt_t word = 0;
1604
1605 TObjString* trigClasseName = 0x0;
1606 TIter nextTrigger(fSelectTriggerClass);
1607 while ((trigClasseName = static_cast<TObjString*>(nextTrigger()))) {
1608
1609 TRegexp GenericTriggerClasseName(trigClasseName->String());
f4d9245c 1610 if (firedTriggerClasses.Contains(GenericTriggerClasseName)) word |= trigClasseName->GetUniqueID();
00612ffc 1611
1612 }
1613
1614 return word;
1615}
1616
4060971a 1617//__________________________________________________________________________
1618void AliAnalysisTaskMuonResolution::CheckPads(AliMUONVCluster *cl, Bool_t &hasBending, Bool_t &hasNonBending) const
1619{
1620 /// Check that this cluster contains pads on both cathods
1621
1622 // reset
1623 hasBending = kFALSE;
1624 hasNonBending = kFALSE;
1625
1626 // loop over digits contained in the cluster
1627 for (Int_t iDigit = 0; iDigit < cl->GetNDigits(); iDigit++) {
1628
1629 Int_t manuId = AliMUONVDigit::ManuId(cl->GetDigitId(iDigit));
1630
1631 // check the location of the manu the digit belongs to
1632 if (manuId > 0) {
1633 if (manuId & AliMpConstants::ManuMask(AliMp::kNonBendingPlane)) hasNonBending = kTRUE;
1634 else hasBending = kTRUE;
1635 }
1636
1637 }
1638
1639}
1640
1641//________________________________________________________________________
1642void AliAnalysisTaskMuonResolution::CheckPadsBelow(AliMUONVCluster *cl, Bool_t &hasBending, Bool_t &hasNonBending) const
1643{
1644 /// Check that this cluster contains pads on both cathods just under its position
1645
1646 // reset
1647 hasBending = kFALSE;
1648 hasNonBending = kFALSE;
1649
1650 // get the cathod corresponding to the bending/non-bending plane
1651 Int_t deId = cl->GetDetElemId();
1652 AliMpDetElement* de = AliMpDDLStore::Instance()->GetDetElement(deId, kFALSE);
1653 if (!de) return;
1654 AliMp::CathodType cath1 = de->GetCathodType(AliMp::kBendingPlane);
1655 AliMp::CathodType cath2 = de->GetCathodType(AliMp::kNonBendingPlane);
1656
1657 // get the corresponding segmentation
1658 const AliMpVSegmentation* seg1 = AliMpSegmentation::Instance()->GetMpSegmentation(deId, cath1);
1659 const AliMpVSegmentation* seg2 = AliMpSegmentation::Instance()->GetMpSegmentation(deId, cath2);
1660 if (!seg1 || !seg2) return;
1661
1662 // get local coordinate of the cluster
1663 Double_t lX,lY,lZ;
1664 Double_t gX = cl->GetX();
1665 Double_t gY = cl->GetY();
1666 Double_t gZ = cl->GetZ();
1667 fNewGeoTransformer->Global2Local(deId,gX,gY,gZ,lX,lY,lZ);
1668
1669 // find pads below the cluster
1670 AliMpPad pad1 = seg1->PadByPosition(lX, lY, kFALSE);
1671 AliMpPad pad2 = seg2->PadByPosition(lX, lY, kFALSE);
1672
1673 // build their ID if pads are valid
1674 UInt_t padId1 = (pad1.IsValid()) ? AliMUONVDigit::BuildUniqueID(deId, pad1.GetManuId(), pad1.GetManuChannel(), cath1) : 0;
1675 UInt_t padId2 = (pad2.IsValid()) ? AliMUONVDigit::BuildUniqueID(deId, pad2.GetManuId(), pad2.GetManuChannel(), cath2) : 0;
1676
1677 // check if the cluster contains these pads
1678 for (Int_t iDigit = 0; iDigit < cl->GetNDigits(); iDigit++) {
1679
1680 UInt_t digitId = cl->GetDigitId(iDigit);
1681
1682 if (digitId == padId1) {
1683
1684 hasBending = kTRUE;
1685 if (hasNonBending) break;
1686
1687 } else if (digitId == padId2) {
1688
1689 hasNonBending = kTRUE;
1690 if (hasBending) break;
1691
1692 }
1693
1694 }
1695
1696}
1697