]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/muondep/AliAnalysisTaskMuonResolution.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
small fix to remove TF1Helper which sneaked in again, inadvertently...
[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
16// ROOT includes
17#include <TFile.h>
18#include <TTree.h>
19#include <TH1F.h>
20#include <TH2F.h>
21#include <TMultiGraph.h>
22#include <TGraphErrors.h>
23#include <TCanvas.h>
24#include <TLegend.h>
25#include <Riostream.h>
26#include <TString.h>
27#include <TGeoManager.h>
28
29// STEER includes
30#include "AliESDEvent.h"
31#include "AliESDMuonTrack.h"
32#include "AliCDBManager.h"
33#include "AliCDBStorage.h"
34#include "AliGeomManager.h"
35
36// ANALYSIS includes
37#include "AliAnalysisTaskSE.h"
38#include "AliAnalysisDataSlot.h"
39#include "AliAnalysisManager.h"
40#include "AliInputEventHandler.h"
41#include "AliAnalysisTaskMuonResolution.h"
42
43// MUON includes
44#include "AliMUONCDB.h"
45#include "AliMUONRecoParam.h"
46#include "AliMUONESDInterface.h"
47#include "AliMUONVTrackReconstructor.h"
48#include "AliMUONTrack.h"
49#include "AliMUONTrackParam.h"
50#include "AliMUONTrackExtrap.h"
51#include "AliMUONVCluster.h"
52#include "AliMUONGeometryTransformer.h"
53#include "AliMUONGeometryModuleTransformer.h"
54#include "AliMUONGeometryDetElement.h"
55#include "AliMpDEIterator.h"
56
57#ifndef SafeDelete
58#define SafeDelete(x) if (x != NULL) { delete x; x = NULL; }
59#endif
60
61ClassImp(AliAnalysisTaskMuonResolution)
62
63const Int_t AliAnalysisTaskMuonResolution::fgkMinEntries = 10;
64
65//________________________________________________________________________
66AliAnalysisTaskMuonResolution::AliAnalysisTaskMuonResolution() :
67 AliAnalysisTaskSE(),
68 fResiduals(NULL),
69 fResidualsVsP(NULL),
70 fLocalChi2(NULL),
71 fChamberRes(NULL),
72 fTrackRes(NULL),
73 fCanvases(NULL),
74 fNEvents(0),
75 fMinMomentum(0.),
76 fSelectPhysics(kFALSE),
77 fMatchTrig(kFALSE),
78 fExtrapMode(1),
79 fCorrectForSystematics(kTRUE),
80 fOCDBLoaded(kFALSE),
81 fNDE(0),
82 fReAlign(kFALSE),
83 fOldAlignStorage(""),
84 fNewAlignStorage(""),
85 fOldGeoTransformer(NULL),
86 fNewGeoTransformer(NULL)
87{
88 /// Default constructor
89}
90
91//________________________________________________________________________
92AliAnalysisTaskMuonResolution::AliAnalysisTaskMuonResolution(const char *name) :
93 AliAnalysisTaskSE(name),
94 fResiduals(NULL),
95 fResidualsVsP(NULL),
96 fLocalChi2(NULL),
97 fChamberRes(NULL),
98 fTrackRes(NULL),
99 fCanvases(NULL),
100 fNEvents(0),
101 fMinMomentum(0.),
102 fSelectPhysics(kFALSE),
103 fMatchTrig(kFALSE),
104 fExtrapMode(1),
105 fCorrectForSystematics(kTRUE),
106 fOCDBLoaded(kFALSE),
107 fNDE(0),
108 fReAlign(kFALSE),
109 fOldAlignStorage(""),
110 fNewAlignStorage(""),
111 fOldGeoTransformer(NULL),
112 fNewGeoTransformer(NULL)
113{
114 /// Constructor
115
116 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) SetStartingResolution(i, -1., -1.);
117
118 // Output slot #1 writes into a TObjArray container
119 DefineOutput(1,TObjArray::Class());
120 // Output slot #2 writes into a TObjArray container
121 DefineOutput(2,TObjArray::Class());
122 // Output slot #3 writes into a TObjArray container
123 DefineOutput(3,TObjArray::Class());
124 // Output slot #4 writes into a TObjArray container
125 DefineOutput(4,TObjArray::Class());
126 // Output slot #5 writes into a TObjArray container
127 DefineOutput(5,TObjArray::Class());
128}
129
130//________________________________________________________________________
131AliAnalysisTaskMuonResolution::~AliAnalysisTaskMuonResolution()
132{
133 /// Destructor
134 SafeDelete(fResiduals);
135 SafeDelete(fResidualsVsP);
136 SafeDelete(fChamberRes);
137 SafeDelete(fTrackRes);
138 SafeDelete(fCanvases);
139 SafeDelete(fOldGeoTransformer);
140 SafeDelete(fNewGeoTransformer);
141}
142
143//___________________________________________________________________________
144void AliAnalysisTaskMuonResolution::UserCreateOutputObjects()
145{
146 /// Create histograms
147
148 // do it once the OCDB has been loaded (i.e. from NotifyRun())
149 if (!fOCDBLoaded) return;
150
151 fResiduals = new TObjArray(1000);
152 fResiduals->SetOwner();
153 fResidualsVsP = new TObjArray(1000);
154 fResidualsVsP->SetOwner();
155 fTrackRes = new TObjArray(1000);
156 fTrackRes->SetOwner();
157 TH2F* h2;
158
159 // find the highest chamber resolution and set histogram bins
160 const AliMUONRecoParam* recoParam = AliMUONESDInterface::GetTracker()->GetRecoParam();
161 Double_t maxSigma[2] = {-1., -1.};
162 for (Int_t i = 0; i < 10; i++) {
163 if (recoParam->GetDefaultNonBendingReso(i) > maxSigma[0]) maxSigma[0] = recoParam->GetDefaultNonBendingReso(i);
164 if (recoParam->GetDefaultBendingReso(i) > maxSigma[1]) maxSigma[1] = recoParam->GetDefaultBendingReso(i);
165 }
166 const char* axes[2] = {"X", "Y"};
167 const Int_t nBins = 2000;
168 const Int_t nSigma = 10;
169 const Int_t pNBins = 20;
170 const Double_t pEdges[2] = {0., 50.};
171
172 for (Int_t ia = 0; ia < 2; ia++) {
173
174 Double_t maxRes = nSigma*maxSigma[ia];
175
176 // List of residual histos
177 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);
178 fResiduals->AddAtAndExpand(h2, kResidualPerCh_ClusterIn+ia);
179 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);
180 fResiduals->AddAtAndExpand(h2, kResidualPerCh_ClusterOut+ia);
181
182 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);
183 fResiduals->AddAtAndExpand(h2, kResidualPerHalfCh_ClusterIn+ia);
184 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);
185 fResiduals->AddAtAndExpand(h2, kResidualPerHalfCh_ClusterOut+ia);
186
187 h2 = new TH2F(Form("hResidual%sPerDE_ClusterIn",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, -maxRes, maxRes);
188 for (Int_t i = 1; i <= fNDE; i++) h2->GetXaxis()->SetBinLabel(i, Form("%d",fDEIds[i]));
189 fResiduals->AddAtAndExpand(h2, kResidualPerDE_ClusterIn+ia);
190 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);
191 for (Int_t i = 1; i <= fNDE; i++) h2->GetXaxis()->SetBinLabel(i, Form("%d",fDEIds[i]));
192 fResiduals->AddAtAndExpand(h2, kResidualPerDE_ClusterOut+ia);
193
194 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);
195 fResiduals->AddAtAndExpand(h2, kTrackResPerCh+ia);
196 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);
197 fResiduals->AddAtAndExpand(h2, kTrackResPerHalfCh+ia);
198 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);
199 for (Int_t i = 1; i <= fNDE; i++) h2->GetXaxis()->SetBinLabel(i, Form("%d",fDEIds[i]));
200 fResiduals->AddAtAndExpand(h2, kTrackResPerDE+ia);
201
202 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);
203 fResiduals->AddAtAndExpand(h2, kMCSPerCh+ia);
204 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);
205 fResiduals->AddAtAndExpand(h2, kMCSPerHalfCh+ia);
206 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);
207 for (Int_t i = 1; i <= fNDE; i++) h2->GetXaxis()->SetBinLabel(i, Form("%d",fDEIds[i]));
208 fResiduals->AddAtAndExpand(h2, kMCSPerDE+ia);
209
210 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);
211 fResiduals->AddAtAndExpand(h2, kClusterRes2PerCh+ia);
212
213 // List of residual vs. p histos
214 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
215 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);
216 fResidualsVsP->AddAtAndExpand(h2, kResidualInChVsP_ClusterIn+10*ia+i);
217 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);
218 fResidualsVsP->AddAtAndExpand(h2, kResidualInChVsP_ClusterOut+10*ia+i);
219 }
220
221 // local chi2
222 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.);
223 fResiduals->AddAtAndExpand(h2, kLocalChi2PerCh+ia);
224 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.);
225 for (Int_t i = 1; i <= fNDE; i++) h2->GetXaxis()->SetBinLabel(i, Form("%d",fDEIds[i]));
226 fResiduals->AddAtAndExpand(h2, kLocalChi2PerDE+ia);
227
228 // track resolution
229 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);
230 fTrackRes->AddAtAndExpand(h2, kUncorrSlopeRes+ia);
231 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);
232 fTrackRes->AddAtAndExpand(h2, kSlopeRes+ia);
233 }
234
235 // local chi2 X+Y
236 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.);
237 fResiduals->AddAtAndExpand(h2, kLocalChi2PerCh+2);
238 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.);
239 for (Int_t i = 1; i <= fNDE; i++) h2->GetXaxis()->SetBinLabel(i, Form("%d",fDEIds[i]));
240 fResiduals->AddAtAndExpand(h2, kLocalChi2PerDE+2);
241
242 // track resolution
243 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.);
244 fTrackRes->AddAtAndExpand(h2, kUncorrPRes);
245 h2 = new TH2F("hPRes", "muon momentum reconstructed resolution at vertex vs p;p (GeV/c); #sigma_{p}/p (%)", 300, 0., 300., 1000, 0., 10.);
246 fTrackRes->AddAtAndExpand(h2, kPRes);
247 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.);
248 fTrackRes->AddAtAndExpand(h2, kUncorrPtRes);
249 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.);
250 fTrackRes->AddAtAndExpand(h2, kPtRes);
251
252 // Post data at least once per task to ensure data synchronisation (required for merging)
253 PostData(1, fResiduals);
254 PostData(2, fResidualsVsP);
255 PostData(5, fTrackRes);
256}
257
258//________________________________________________________________________
259void AliAnalysisTaskMuonResolution::UserExec(Option_t *)
260{
261 /// Main event loop
262
263 // check if OCDB properly loaded
264 if (!fOCDBLoaded) return;
265
266 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
267 if (!esd) return;
268
269 if ((++fNEvents)%100 == 0) cout<<"\rEvent processing... "<<fNEvents<<"\r"<<flush;
270
271 // get tracker to refit
272 AliMUONVTrackReconstructor* tracker = AliMUONESDInterface::GetTracker();
273
274 // loop over tracks
275 Int_t nTracks = (Int_t) esd->GetNumberOfMuonTracks();
276 for (Int_t iTrack = 0; iTrack < nTracks; ++iTrack) {
277
278 // get the ESD track
279 AliESDMuonTrack* esdTrack = esd->GetMuonTrack(iTrack);
280
281 // skip ghost tracks
282 if (!esdTrack->ContainTrackerData()) continue;
283
284 // skip tracks that do not pass the physics selection if required
285 if (fSelectPhysics && fInputHandler && !fInputHandler->IsEventSelected()) continue;
286
287 // skip tracks not matched with trigger if required
288 if (fMatchTrig && !esdTrack->ContainTriggerData()) continue;
289
290 // skip low momentum tracks
291 if (esdTrack->PUncorrected() < fMinMomentum) continue;
292
293 // get the corresponding MUON track
294 AliMUONTrack track;
295 AliMUONESDInterface::ESDToMUON(*esdTrack, track, kFALSE);
296
297 // change the cluster resolution (and position)
298 ModifyClusters(track);
299
300 // refit the track
301 if (!tracker->RefitTrack(track, kFALSE)) break;
302
303 // save track unchanged
304 AliMUONTrack referenceTrack(track);
305
306 // get track param at first cluster and add MCS in first chamber
307 AliMUONTrackParam trackParamAtFirstCluster(*(static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->First())));
308 Int_t firstCh = 0; while (firstCh < 10 && !esdTrack->IsInMuonClusterMap(firstCh)) firstCh++;
309 AliMUONTrackExtrap::AddMCSEffect(&trackParamAtFirstCluster, AliMUONConstants::ChamberThicknessInX0(firstCh)/2., -1.);
310
311 // fill momentum error at first cluster
312 Double_t pXUncorr = trackParamAtFirstCluster.Px();
313 Double_t pYUncorr = trackParamAtFirstCluster.Py();
314 Double_t pZUncorr = trackParamAtFirstCluster.Pz();
315 Double_t pUncorr = trackParamAtFirstCluster.P();
316 TMatrixD covUncorr(5,5);
317 Cov2CovP(trackParamAtFirstCluster,covUncorr);
318 Double_t sigmaPUncorr = TMath::Sqrt(pXUncorr * (pXUncorr*covUncorr(2,2) + pYUncorr*covUncorr(2,3) + pZUncorr*covUncorr(2,4)) +
319 pYUncorr * (pXUncorr*covUncorr(3,2) + pYUncorr*covUncorr(3,3) + pZUncorr*covUncorr(3,4)) +
320 pZUncorr * (pXUncorr*covUncorr(4,2) + pYUncorr*covUncorr(4,3) + pZUncorr*covUncorr(4,4))) / pUncorr;
321 ((TH2F*)fTrackRes->UncheckedAt(kUncorrPRes))->Fill(pUncorr,100.*sigmaPUncorr/pUncorr);
322
323 // fill transverse momentum error at first cluster
324 Double_t ptUncorr = TMath::Sqrt(pXUncorr*pXUncorr + pYUncorr*pYUncorr);
325 Double_t sigmaPtUncorr = TMath::Sqrt(pXUncorr * (pXUncorr*covUncorr(2,2) + pYUncorr*covUncorr(2,3)) + pYUncorr * (pXUncorr*covUncorr(3,2) + pYUncorr*covUncorr(3,3))) / ptUncorr;
326 ((TH2F*)fTrackRes->UncheckedAt(kUncorrPtRes))->Fill(ptUncorr,100.*sigmaPtUncorr/ptUncorr);
327
328 // fill slopeX-Y error at first cluster
329 ((TH2F*)fTrackRes->UncheckedAt(kUncorrSlopeRes))->Fill(pUncorr,TMath::Sqrt(trackParamAtFirstCluster.GetCovariances()(1,1)));
330 ((TH2F*)fTrackRes->UncheckedAt(kUncorrSlopeRes+1))->Fill(pUncorr,TMath::Sqrt(trackParamAtFirstCluster.GetCovariances()(3,3)));
331
332 // fill momentum error at vertex
333 AliMUONTrackParam trackParamAtVtx(trackParamAtFirstCluster);
334 AliMUONTrackExtrap::ExtrapToVertex(&trackParamAtVtx, esdTrack->GetNonBendingCoor(), esdTrack->GetBendingCoor(), esdTrack->GetZ(), 0., 0.);
335 Double_t pXVtx = trackParamAtVtx.Px();
336 Double_t pYVtx = trackParamAtVtx.Py();
337 Double_t pZVtx = trackParamAtVtx.Pz();
338 Double_t pVtx = trackParamAtVtx.P();
339 TMatrixD covVtx(5,5);
340 Cov2CovP(trackParamAtVtx,covVtx);
341 Double_t sigmaPVtx = TMath::Sqrt(pXVtx * (pXVtx*covVtx(2,2) + pYVtx*covVtx(2,3) + pZVtx*covVtx(2,4)) +
342 pYVtx * (pXVtx*covVtx(3,2) + pYVtx*covVtx(3,3) + pZVtx*covVtx(3,4)) +
343 pZVtx * (pXVtx*covVtx(4,2) + pYVtx*covVtx(4,3) + pZVtx*covVtx(4,4))) / pVtx;
344 ((TH2F*)fTrackRes->UncheckedAt(kPRes))->Fill(pVtx,100.*sigmaPVtx/pVtx);
345
346 // fill transverse momentum error at vertex
347 Double_t ptVtx = TMath::Sqrt(pXVtx*pXVtx + pYVtx*pYVtx);
348 Double_t sigmaPtVtx = TMath::Sqrt(pXVtx * (pXVtx*covVtx(2,2) + pYVtx*covVtx(2,3)) + pYVtx * (pXVtx*covVtx(3,2) + pYVtx*covVtx(3,3))) / ptVtx;
349 ((TH2F*)fTrackRes->UncheckedAt(kPtRes))->Fill(ptVtx,100.*sigmaPtVtx/ptVtx);
350
351 // fill slopeX-Y error at vertex
352 ((TH2F*)fTrackRes->UncheckedAt(kSlopeRes))->Fill(pVtx,TMath::Sqrt(trackParamAtVtx.GetCovariances()(1,1)));
353 ((TH2F*)fTrackRes->UncheckedAt(kSlopeRes+1))->Fill(pVtx,TMath::Sqrt(trackParamAtVtx.GetCovariances()(3,3)));
354
355 // loop over clusters
356 Int_t nClusters = track.GetNClusters();
357 for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
358
359 // Get current, previous and next trackParams
360 AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->UncheckedAt(iCluster));
361 AliMUONTrackParam* previousTrackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->Before(trackParam));
362 AliMUONTrackParam* nextTrackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->After(trackParam));
363
364 // save current trackParam and remove it from the track
365 AliMUONTrackParam currentTrackParam(*trackParam);
366 track.RemoveTrackParamAtCluster(trackParam);
367
368 // get cluster info
369 AliMUONVCluster* cluster = currentTrackParam.GetClusterPtr();
370 Int_t chId = cluster->GetChamberId();
371 Int_t halfChId = (cluster->GetX() > 0) ? 2*chId : 2*chId+1;
372 Int_t deId = cluster->GetDetElemId();
373
374 // compute residuals with cluster still attached to the track
375 AliMUONTrackParam* referenceTrackParam = static_cast<AliMUONTrackParam*>(referenceTrack.GetTrackParamAtCluster()->UncheckedAt(iCluster));
376 Double_t deltaX = cluster->GetX() - referenceTrackParam->GetNonBendingCoor();
377 Double_t deltaY = cluster->GetY() - referenceTrackParam->GetBendingCoor();
378
379 // compute local chi2
380 Double_t sigmaDeltaX2 = cluster->GetErrX2() - referenceTrackParam->GetCovariances()(0,0);
381 Double_t sigmaDeltaY2 = cluster->GetErrY2() - referenceTrackParam->GetCovariances()(2,2);
382 Double_t localChi2X = (sigmaDeltaX2 > 0.) ? deltaX*deltaX/sigmaDeltaX2 : 0.;
383 Double_t localChi2Y = (sigmaDeltaY2 > 0.) ? deltaY*deltaY/sigmaDeltaY2 : 0.;
384 Double_t localChi2 = 0.5 * referenceTrackParam->GetLocalChi2();
385
386 // fill local chi2 info at every clusters
387 ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerCh))->Fill(chId+1, localChi2X);
388 ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerCh+1))->Fill(chId+1, localChi2Y);
389 ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerCh+2))->Fill(chId+1, localChi2);
390 ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerDE))->Fill(fDEIndices[deId], localChi2X);
391 ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerDE+1))->Fill(fDEIndices[deId], localChi2Y);
392 ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerDE+2))->Fill(fDEIndices[deId], localChi2);
393
394 // make sure the track has another cluster in the same station and can still be refitted
395 Bool_t refit = track.IsValid( 1 << (chId/2) );
396 if (refit) {
397
398 // refit the track and proceed if everything goes fine
399 if (tracker->RefitTrack(track, kFALSE)) {
400
401 // fill histograms of residuals with cluster still attached to the track
402 ((TH2F*)fResiduals->UncheckedAt(kResidualPerCh_ClusterIn))->Fill(chId+1, deltaX);
403 ((TH2F*)fResiduals->UncheckedAt(kResidualPerCh_ClusterIn+1))->Fill(chId+1, deltaY);
404 ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfCh_ClusterIn))->Fill(halfChId+1, deltaX);
405 ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfCh_ClusterIn+1))->Fill(halfChId+1, deltaY);
406 ((TH2F*)fResiduals->UncheckedAt(kResidualPerDE_ClusterIn))->Fill(fDEIndices[deId], deltaX);
407 ((TH2F*)fResiduals->UncheckedAt(kResidualPerDE_ClusterIn+1))->Fill(fDEIndices[deId], deltaY);
408 ((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsP_ClusterIn+chId))->Fill(pUncorr, deltaX);
409 ((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsP_ClusterIn+10+chId))->Fill(pUncorr, deltaY);
410
411 // find the track parameters closest to the current cluster position
412 Double_t dZWithPrevious = (previousTrackParam) ? TMath::Abs(previousTrackParam->GetClusterPtr()->GetZ() - cluster->GetZ()) : FLT_MAX;
413 Int_t previousChId = (previousTrackParam) ? previousTrackParam->GetClusterPtr()->GetChamberId() : -1;
414 Double_t dZWithNext = (nextTrackParam) ? TMath::Abs(nextTrackParam->GetClusterPtr()->GetZ() - cluster->GetZ()) : FLT_MAX;
415 AliMUONTrackParam* startingTrackParam = 0x0;
416 if ((fExtrapMode == 0 && dZWithPrevious < dZWithNext) ||
417 (fExtrapMode == 1 && previousTrackParam && !(chId/2 == 2 && previousChId/2 == 1) &&
418 !(chId/2 == 3 && previousChId/2 == 2))) startingTrackParam = previousTrackParam;
419 else startingTrackParam = nextTrackParam;
420
421 // reset current parameters
422 currentTrackParam.SetParameters(startingTrackParam->GetParameters());
423 currentTrackParam.SetZ(startingTrackParam->GetZ());
424 currentTrackParam.SetCovariances(startingTrackParam->GetCovariances());
425 currentTrackParam.ResetPropagator();
426
427 // extrapolate to the current cluster position and fill histograms of residuals if everything goes fine
428 if (AliMUONTrackExtrap::ExtrapToZCov(&currentTrackParam, currentTrackParam.GetClusterPtr()->GetZ(), kTRUE)) {
429
430 // compute MCS dispersion on the first chamber
431 TMatrixD mcsCov(5,5);
432 if (startingTrackParam == nextTrackParam && chId == 0) {
433 AliMUONTrackParam trackParamForMCS;
434 trackParamForMCS.SetParameters(nextTrackParam->GetParameters());
435 AliMUONTrackExtrap::AddMCSEffect(&trackParamForMCS,AliMUONConstants::ChamberThicknessInX0(nextTrackParam->GetClusterPtr()->GetChamberId()),-1.);
436 const TMatrixD &propagator = currentTrackParam.GetPropagator();
437 TMatrixD tmp(trackParamForMCS.GetCovariances(),TMatrixD::kMultTranspose,propagator);
438 mcsCov.Mult(propagator,tmp);
439 } else mcsCov.Zero();
440
441 // compute residuals
442 Double_t trackResX2 = currentTrackParam.GetCovariances()(0,0) + mcsCov(0,0);
443 Double_t trackResY2 = currentTrackParam.GetCovariances()(2,2) + mcsCov(2,2);
444 deltaX = cluster->GetX() - currentTrackParam.GetNonBendingCoor();
445 deltaY = cluster->GetY() - currentTrackParam.GetBendingCoor();
446
447 // fill histograms with cluster not attached to the track
448 ((TH2F*)fResiduals->UncheckedAt(kResidualPerCh_ClusterOut))->Fill(chId+1, deltaX);
449 ((TH2F*)fResiduals->UncheckedAt(kResidualPerCh_ClusterOut+1))->Fill(chId+1, deltaY);
450 ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfCh_ClusterOut))->Fill(halfChId+1, deltaX);
451 ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfCh_ClusterOut+1))->Fill(halfChId+1, deltaY);
452 ((TH2F*)fResiduals->UncheckedAt(kResidualPerDE_ClusterOut))->Fill(fDEIndices[deId], deltaX);
453 ((TH2F*)fResiduals->UncheckedAt(kResidualPerDE_ClusterOut+1))->Fill(fDEIndices[deId], deltaY);
454 ((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsP_ClusterOut+chId))->Fill(pUncorr, deltaX);
455 ((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsP_ClusterOut+10+chId))->Fill(pUncorr, deltaY);
456 ((TH2F*)fResiduals->UncheckedAt(kTrackResPerCh))->Fill(chId+1, TMath::Sqrt(trackResX2));
457 ((TH2F*)fResiduals->UncheckedAt(kTrackResPerCh+1))->Fill(chId+1, TMath::Sqrt(trackResY2));
458 ((TH2F*)fResiduals->UncheckedAt(kTrackResPerHalfCh))->Fill(halfChId+1, TMath::Sqrt(trackResX2));
459 ((TH2F*)fResiduals->UncheckedAt(kTrackResPerHalfCh+1))->Fill(halfChId+1, TMath::Sqrt(trackResY2));
460 ((TH2F*)fResiduals->UncheckedAt(kTrackResPerDE))->Fill(fDEIndices[deId], TMath::Sqrt(trackResX2));
461 ((TH2F*)fResiduals->UncheckedAt(kTrackResPerDE+1))->Fill(fDEIndices[deId], TMath::Sqrt(trackResY2));
462 ((TH2F*)fResiduals->UncheckedAt(kMCSPerCh))->Fill(chId+1, TMath::Sqrt(mcsCov(0,0)));
463 ((TH2F*)fResiduals->UncheckedAt(kMCSPerCh+1))->Fill(chId+1, TMath::Sqrt(mcsCov(2,2)));
464 ((TH2F*)fResiduals->UncheckedAt(kMCSPerHalfCh))->Fill(halfChId+1, TMath::Sqrt(mcsCov(0,0)));
465 ((TH2F*)fResiduals->UncheckedAt(kMCSPerHalfCh+1))->Fill(halfChId+1, TMath::Sqrt(mcsCov(2,2)));
466 ((TH2F*)fResiduals->UncheckedAt(kMCSPerDE))->Fill(fDEIndices[deId], TMath::Sqrt(mcsCov(0,0)));
467 ((TH2F*)fResiduals->UncheckedAt(kMCSPerDE+1))->Fill(fDEIndices[deId], TMath::Sqrt(mcsCov(2,2)));
468 ((TH2F*)fResiduals->UncheckedAt(kClusterRes2PerCh))->Fill(chId+1, deltaX*deltaX - trackResX2);
469 ((TH2F*)fResiduals->UncheckedAt(kClusterRes2PerCh+1))->Fill(chId+1, deltaY*deltaY - trackResY2);
470 }
471
472 }
473
474 }
475
476 // restore the track
477 track.AddTrackParamAtCluster(currentTrackParam, *(currentTrackParam.GetClusterPtr()), kTRUE);
478
479 }
480
481 }
482
483 // Post final data. It will be written to a file with option "RECREATE"
484 PostData(1, fResiduals);
485 PostData(2, fResidualsVsP);
486 PostData(5, fTrackRes);
487}
488
489//________________________________________________________________________
490void AliAnalysisTaskMuonResolution::NotifyRun()
491{
492 /// load necessary data from OCDB corresponding to the first run number and initialize analysis
493
494 if (fOCDBLoaded) return;
495
496 AliCDBManager* cdbm = AliCDBManager::Instance();
497 cdbm->SetRun(fCurrentRunNumber);
498
499 if (!AliMUONCDB::LoadField()) return;
500
501 if (!AliMUONCDB::LoadMapping()) return;
502
503 AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam();
504 if (!recoParam) return;
505
506 fOCDBLoaded = kTRUE;
507
508 AliMUONESDInterface::ResetTracker(recoParam);
509
510 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
511
512 // set the cluster resolution to default if not already set and create output objects
513 if (fClusterResNB[i] < 0.) fClusterResNB[i] = recoParam->GetDefaultNonBendingReso(i);
514 if (fClusterResB[i] < 0.) fClusterResB[i] = recoParam->GetDefaultBendingReso(i);
515
516 // fill correspondence between DEId and indices for histo (starting from 1)
517 AliMpDEIterator it;
518 it.First(i);
519 while (!it.IsDone()) {
520 fNDE++;
521 fDEIndices[it.CurrentDEId()] = fNDE;
522 fDEIds[fNDE] = it.CurrentDEId();
523 it.Next();
524 }
525
526 }
527
528 if (fReAlign) {
529
530 // recover default storage name
531 TString defaultStorage(cdbm->GetDefaultStorage()->GetType());
532 if (defaultStorage == "alien") defaultStorage += Form("://folder=%s", cdbm->GetDefaultStorage()->GetBaseFolder().Data());
533 else defaultStorage += Form("://%s", cdbm->GetDefaultStorage()->GetBaseFolder().Data());
534
535 // reset existing geometry/alignment if any
536 if (cdbm->GetEntryCache()->Contains("GRP/Geometry/Data")) cdbm->UnloadFromCache("GRP/Geometry/Data");
537 if (cdbm->GetEntryCache()->Contains("MUON/Align/Data")) cdbm->UnloadFromCache("MUON/Align/Data");
538 if (AliGeomManager::GetGeometry()) AliGeomManager::GetGeometry()->UnlockGeometry();
539
540 // get original geometry transformer
541 AliGeomManager::LoadGeometry();
542 if (!AliGeomManager::GetGeometry()) return;
543 if (fOldAlignStorage != "none") {
544 if (!fOldAlignStorage.IsNull()) cdbm->SetSpecificStorage("MUON/Align/Data",fOldAlignStorage.Data());
545 else cdbm->SetSpecificStorage("MUON/Align/Data",defaultStorage.Data());
546 AliGeomManager::ApplyAlignObjsFromCDB("MUON");
547 }
548 fOldGeoTransformer = new AliMUONGeometryTransformer();
549 fOldGeoTransformer->LoadGeometryData();
550
551 // get new geometry transformer
552 cdbm->UnloadFromCache("GRP/Geometry/Data");
553 if (fOldAlignStorage != "none") cdbm->UnloadFromCache("MUON/Align/Data");
554 AliGeomManager::GetGeometry()->UnlockGeometry();
555 AliGeomManager::LoadGeometry();
556 if (!AliGeomManager::GetGeometry()) return;
557 if (!fNewAlignStorage.IsNull()) cdbm->SetSpecificStorage("MUON/Align/Data",fNewAlignStorage.Data());
558 else cdbm->SetSpecificStorage("MUON/Align/Data",defaultStorage.Data());
559 AliGeomManager::ApplyAlignObjsFromCDB("MUON");
560 fNewGeoTransformer = new AliMUONGeometryTransformer();
561 fNewGeoTransformer->LoadGeometryData();
562
563 } else {
564
565 // load geometry for track extrapolation to vertex
566 if (cdbm->GetEntryCache()->Contains("GRP/Geometry/Data")) cdbm->UnloadFromCache("GRP/Geometry/Data");
567 if (AliGeomManager::GetGeometry()) AliGeomManager::GetGeometry()->UnlockGeometry();
568 AliGeomManager::LoadGeometry();
569 if (!AliGeomManager::GetGeometry()) return;
570
571 }
572
573 // print starting chamber resolution
574 printf("\nstarting chamber resolution:\n");
575 printf(" - non-bending:");
576 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %5.3f":", %5.3f",fClusterResNB[i]);
577 printf("\n - bending:");
578 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %6.4f":", %6.4f",fClusterResB[i]);
579 printf("\n\n");
580
581 UserCreateOutputObjects();
582
583}
584
585//________________________________________________________________________
586void AliAnalysisTaskMuonResolution::Terminate(Option_t *)
587{
588 /// compute final results
589
590 // recover output objects
591 fResiduals = static_cast<TObjArray*> (GetOutputData(1));
592 fResidualsVsP = static_cast<TObjArray*> (GetOutputData(2));
593 fTrackRes = static_cast<TObjArray*> (GetOutputData(5));
594 if (!fResiduals || !fResidualsVsP || !fTrackRes) return;
595
596 // summary graphs
597 fLocalChi2 = new TObjArray(1000);
598 fLocalChi2->SetOwner();
599 fChamberRes = new TObjArray(1000);
600 fChamberRes->SetOwner();
601 TGraphErrors* g;
602 TMultiGraph* mg;
603
604 const char* axes[2] = {"X", "Y"};
605 Double_t newClusterRes[2][10], newClusterResErr[2][10];
606 fNDE = ((TH2F*)fResiduals->UncheckedAt(kResidualPerDE_ClusterIn))->GetXaxis()->GetNbins();
607
608 for (Int_t ia = 0; ia < 2; ia++) {
609
610 g = new TGraphErrors(AliMUONConstants::NTrackingCh());
611 g->SetName(Form("gResidual%sPerChMean_ClusterIn",axes[ia]));
612 g->SetTitle(Form("cluster-track residual-%s per Ch: mean (cluster in);chamber ID;<#Delta_{%s}> (cm)",axes[ia],axes[ia]));
613 g->SetMarkerStyle(kFullDotLarge);
614 fChamberRes->AddAtAndExpand(g, kResidualPerChMean_ClusterIn+ia);
615 g = new TGraphErrors(AliMUONConstants::NTrackingCh());
616 g->SetName(Form("gResidual%sPerChMean_ClusterOut",axes[ia]));
617 g->SetTitle(Form("cluster-track residual-%s per Ch: mean (cluster out);chamber ID;<#Delta_{%s}> (cm)",axes[ia],axes[ia]));
618 g->SetMarkerStyle(kFullDotLarge);
619 fChamberRes->AddAtAndExpand(g, kResidualPerChMean_ClusterOut+ia);
620
621 g = new TGraphErrors(2*AliMUONConstants::NTrackingCh());
622 g->SetName(Form("gResidual%sPerHalfChMean_ClusterIn",axes[ia]));
623 g->SetTitle(Form("cluster-track residual-%s per half Ch: mean (cluster in);half chamber ID;<#Delta_{%s}> (cm)",axes[ia],axes[ia]));
624 g->SetMarkerStyle(kFullDotLarge);
625 fChamberRes->AddAtAndExpand(g, kResidualPerHalfChMean_ClusterIn+ia);
626 g = new TGraphErrors(2*AliMUONConstants::NTrackingCh());
627 g->SetName(Form("gResidual%sPerHalfChMean_ClusterOut",axes[ia]));
628 g->SetTitle(Form("cluster-track residual-%s per half Ch: mean (cluster out);half chamber ID;<#Delta_{%s}> (cm)",axes[ia],axes[ia]));
629 g->SetMarkerStyle(kFullDotLarge);
630 fChamberRes->AddAtAndExpand(g, kResidualPerHalfChMean_ClusterOut+ia);
631
632 g = new TGraphErrors(fNDE);
633 g->SetName(Form("gResidual%sPerDEMean_ClusterIn",axes[ia]));
634 g->SetTitle(Form("cluster-track residual-%s per DE: mean (cluster in);DE ID;<#Delta_{%s}> (cm)",axes[ia],axes[ia]));
635 g->SetMarkerStyle(kFullDotLarge);
636 fChamberRes->AddAtAndExpand(g, kResidualPerDEMean_ClusterIn+ia);
637 g = new TGraphErrors(fNDE);
638 g->SetName(Form("gResidual%sPerDEMean_ClusterOut",axes[ia]));
639 g->SetTitle(Form("cluster-track residual-%s per DE: mean (cluster out);DE ID;<#Delta_{%s}> (cm)",axes[ia],axes[ia]));
640 g->SetMarkerStyle(kFullDotLarge);
641 fChamberRes->AddAtAndExpand(g, kResidualPerDEMean_ClusterOut+ia);
642
643 g = new TGraphErrors(AliMUONConstants::NTrackingCh());
644 g->SetName(Form("gResidual%sPerChSigma_ClusterIn",axes[ia]));
645 g->SetTitle(Form("cluster-track residual-%s per Ch: sigma (cluster in);chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
646 g->SetMarkerStyle(kFullDotLarge);
647 fChamberRes->AddAtAndExpand(g, kResidualPerChSigma_ClusterIn+ia);
648 g = new TGraphErrors(AliMUONConstants::NTrackingCh());
649 g->SetName(Form("gResidual%sPerChSigma_ClusterOut",axes[ia]));
650 g->SetTitle(Form("cluster-track residual-%s per Ch: sigma (cluster out);chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
651 g->SetMarkerStyle(kFullDotLarge);
652 fChamberRes->AddAtAndExpand(g, kResidualPerChSigma_ClusterOut+ia);
653
654 g = new TGraphErrors(AliMUONConstants::NTrackingCh());
655 g->SetName(Form("gResidual%sPerChDispersion_ClusterOut",axes[ia]));
656 g->SetTitle(Form("cluster-track residual-%s per Ch: dispersion (cluster out);chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
657 g->SetMarkerStyle(kFullDotLarge);
658 fChamberRes->AddAtAndExpand(g, kResidualPerChDispersion_ClusterOut+ia);
659
660 g = new TGraphErrors(AliMUONConstants::NTrackingCh());
661 g->SetName(Form("gCombinedResidual%sPerChSigma",axes[ia]));
662 g->SetTitle(Form("combined cluster-track residual-%s per Ch: sigma;chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
663 g->SetMarkerStyle(kFullDotLarge);
664 fChamberRes->AddAtAndExpand(g, kCombinedResidualPerChSigma+ia);
665
666 g = new TGraphErrors(2*AliMUONConstants::NTrackingCh());
667 g->SetName(Form("gCombinedResidual%sPerHalfChSigma",axes[ia]));
668 g->SetTitle(Form("combined cluster-track residual-%s per half Ch: sigma;half chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
669 g->SetMarkerStyle(kFullDotLarge);
670 fChamberRes->AddAtAndExpand(g, kCombinedResidualPerHalfChSigma+ia);
671
672 g = new TGraphErrors(fNDE);
673 g->SetName(Form("gCombinedResidual%sPerDESigma",axes[ia]));
674 g->SetTitle(Form("combined cluster-track residual-%s per DE: sigma;DE ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
675 g->SetMarkerStyle(kFullDotLarge);
676 fChamberRes->AddAtAndExpand(g, kCombinedResidualPerDESigma+ia);
677
678 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]));
679 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
680 g = new TGraphErrors(((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsP_ClusterIn+10*ia+i))->GetNbinsX());
681 g->SetName(Form("gRes%sVsP_ch%d",axes[ia],i+1));
682 g->SetMarkerStyle(kFullDotMedium);
683 g->SetMarkerColor(i+1+i/9);
684 mg->Add(g,"p");
685 }
686 fChamberRes->AddAtAndExpand(mg, kCombinedResidualSigmaVsP+ia);
687
688 g = new TGraphErrors(AliMUONConstants::NTrackingCh());
689 g->SetName(Form("gTrackRes%sPerCh",axes[ia]));
690 g->SetTitle(Form("track <#sigma_{%s}> per Ch;chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
691 g->SetMarkerStyle(kFullDotLarge);
692 fChamberRes->AddAtAndExpand(g, kTrackResPerChMean+ia);
693
694 g = new TGraphErrors(AliMUONConstants::NTrackingCh());
695 g->SetName(Form("gMCS%sPerCh",axes[ia]));
696 g->SetTitle(Form("MCS %s-dispersion of extrapolated track per Ch;chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
697 g->SetMarkerStyle(kFullDotLarge);
698 fChamberRes->AddAtAndExpand(g, kMCSPerChMean+ia);
699
700 g = new TGraphErrors(AliMUONConstants::NTrackingCh());
701 g->SetName(Form("gClusterRes%sPerCh",axes[ia]));
702 g->SetTitle(Form("cluster <#sigma_{%s}> per Ch;chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
703 g->SetMarkerStyle(kFullDotLarge);
704 fChamberRes->AddAtAndExpand(g, kClusterResPerCh+ia);
705
706 g = new TGraphErrors(2*AliMUONConstants::NTrackingCh());
707 g->SetName(Form("gClusterRes%sPerHalfCh",axes[ia]));
708 g->SetTitle(Form("cluster <#sigma_{%s}> per half Ch;half chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
709 g->SetMarkerStyle(kFullDotLarge);
710 fChamberRes->AddAtAndExpand(g, kClusterResPerHalfCh+ia);
711
712 g = new TGraphErrors(fNDE);
713 g->SetName(Form("gClusterRes%sPerDE",axes[ia]));
714 g->SetTitle(Form("cluster <#sigma_{%s}> per DE;DE ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
715 g->SetMarkerStyle(kFullDotLarge);
716 fChamberRes->AddAtAndExpand(g, kClusterResPerDE+ia);
717
718 g = new TGraphErrors(AliMUONConstants::NTrackingCh());
719 g->SetName(Form("gCalcClusterRes%sPerCh",axes[ia]));
720 g->SetTitle(Form("calculated cluster <#sigma_{%s}> per Ch;chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
721 g->SetMarkerStyle(kFullDotLarge);
722 fChamberRes->AddAtAndExpand(g, kCalcClusterResPerCh+ia);
723
724 g = new TGraphErrors(AliMUONConstants::NTrackingCh());
725 g->SetName(Form("gLocalChi2%sPerChMean",axes[ia]));
726 g->SetTitle(Form("local chi2-%s per Ch: mean;chamber ID;<local #chi^{2}_{%s}>",axes[ia],axes[ia]));
727 g->SetMarkerStyle(kFullDotLarge);
728 fLocalChi2->AddAtAndExpand(g, kLocalChi2PerChMean+ia);
729
730 g = new TGraphErrors(fNDE);
731 g->SetName(Form("gLocalChi2%sPerDEMean",axes[ia]));
732 g->SetTitle(Form("local chi2-%s per DE: mean;DE ID;<local #chi^{2}_{%s}>",axes[ia],axes[ia]));
733 g->SetMarkerStyle(kFullDotLarge);
734 fLocalChi2->AddAtAndExpand(g, kLocalChi2PerDEMean+ia);
735
736 // compute residual mean and dispersion and averaged local chi2 per chamber and half chamber
737 Double_t meanIn, meanInErr, meanOut, meanOutErr, sigmaIn, sigmaInErr, sigmaOut, sigmaOutErr;
738 Double_t sigmaTrack, sigmaTrackErr, sigmaMCS, sigmaMCSErr, clusterRes, clusterResErr, sigmaCluster, sigmaClusterErr;
739 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
740
741 // method 1
742 TH1D *tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerCh_ClusterIn+ia))->ProjectionY("tmp",i+1,i+1,"e");
743 GetMean(tmp, meanIn, meanInErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChMean_ClusterIn+ia), i, i+1);
744 GetRMS(tmp, sigmaIn, sigmaInErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChSigma_ClusterIn+ia), i, i+1);
745 delete tmp;
746
747 tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerCh_ClusterOut+ia))->ProjectionY("tmp",i+1,i+1,"e");
748 GetMean(tmp, meanOut, meanOutErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChMean_ClusterOut+ia), i, i+1);
749 GetRMS(tmp, sigmaOut, sigmaOutErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChSigma_ClusterOut+ia), i, i+1);
750 delete tmp;
751
752 if (fCorrectForSystematics) {
753 sigmaIn = TMath::Sqrt(sigmaIn*sigmaIn + meanIn*meanIn);
754 sigmaInErr = (sigmaIn>0) ? TMath::Sqrt(sigmaIn*sigmaIn*sigmaInErr*sigmaInErr + meanIn*meanIn*meanInErr*meanInErr) / sigmaIn : 0.;
755 sigmaOut = TMath::Sqrt(sigmaOut*sigmaOut + meanOut*meanOut);
756 sigmaOutErr = (sigmaOut>0) ? TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + meanOut*meanOut*meanOutErr*meanOutErr) / sigmaOut : 0.;
757 }
758 ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChDispersion_ClusterOut+ia))->SetPoint(i, i+1, sigmaOut);
759 ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChDispersion_ClusterOut+ia))->SetPointError(i, 0., sigmaOutErr);
760
761 clusterRes = TMath::Sqrt(sigmaIn*sigmaOut);
762 // clusterResErr = (clusterRes > 0.) ? 0.5 * TMath::Sqrt(sigmaInErr*sigmaInErr*sigmaOut*sigmaOut + sigmaIn*sigmaIn*sigmaOutErr*sigmaOutErr) / clusterRes : 0.;
763 clusterResErr = TMath::Sqrt(sigmaInErr*sigmaOutErr);
764 ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerChSigma+ia))->SetPoint(i, i+1, clusterRes);
765 ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerChSigma+ia))->SetPointError(i, 0., clusterResErr);
766 newClusterRes[ia][i] = clusterRes;
767 newClusterResErr[ia][i] = clusterResErr;
768
769 // method 2
770 tmp = ((TH2F*)fResiduals->UncheckedAt(kTrackResPerCh+ia))->ProjectionY("tmp",i+1,i+1,"e");
771 GetMean(tmp, sigmaTrack, sigmaTrackErr, (TGraphErrors*)fChamberRes->UncheckedAt(kTrackResPerChMean+ia), i, i+1, kFALSE);
772 delete tmp;
773
774 tmp = ((TH2F*)fResiduals->UncheckedAt(kMCSPerCh+ia))->ProjectionY("tmp",i+1,i+1,"e");
775 GetMean(tmp, sigmaMCS, sigmaMCSErr, (TGraphErrors*)fChamberRes->UncheckedAt(kMCSPerChMean+ia), i, i+1, kFALSE);
776 delete tmp;
777
778 sigmaCluster = sigmaOut*sigmaOut - sigmaTrack*sigmaTrack;
779 if (sigmaCluster > 0.) {
780 sigmaCluster = TMath::Sqrt(sigmaCluster);
781 sigmaClusterErr = TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + sigmaTrack*sigmaTrack*sigmaTrackErr*sigmaTrackErr) / sigmaCluster;
782 } else {
783 sigmaCluster = 0.;
784 sigmaClusterErr = 0.;
785 }
786 ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerCh+ia))->SetPoint(i, i+1, sigmaCluster);
787 ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerCh+ia))->SetPointError(i, 0., sigmaClusterErr);
788
789 // method 3
790 tmp = ((TH2F*)fResiduals->UncheckedAt(kClusterRes2PerCh+ia))->ProjectionY("tmp",i+1,i+1,"e");
791 ZoomRight(tmp);
792 clusterRes = tmp->GetMean();
793 if (clusterRes > 0.) {
794 ((TGraphErrors*)fChamberRes->UncheckedAt(kCalcClusterResPerCh+ia))->SetPoint(i, i+1, TMath::Sqrt(clusterRes));
795 ((TGraphErrors*)fChamberRes->UncheckedAt(kCalcClusterResPerCh+ia))->SetPointError(i, 0., 0.5 * tmp->GetMeanError() / TMath::Sqrt(clusterRes));
796 } else {
797 ((TGraphErrors*)fChamberRes->UncheckedAt(kCalcClusterResPerCh+ia))->SetPoint(i, i+1, 0.);
798 ((TGraphErrors*)fChamberRes->UncheckedAt(kCalcClusterResPerCh+ia))->SetPointError(i, 0., 0.);
799 }
800 delete tmp;
801
802 // method 1 versus p
803 FillSigmaClusterVsP((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsP_ClusterIn+10*ia+i),
804 (TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsP_ClusterOut+10*ia+i),
805 (TGraphErrors*)((TMultiGraph*)fChamberRes->UncheckedAt(kCombinedResidualSigmaVsP+ia))->GetListOfGraphs()->FindObject(Form("gRes%sVsP_ch%d",axes[ia],i+1)));
806
807 // compute residual mean and dispersion per half chamber
808 for (Int_t j = 0; j < 2; j++) {
809 Int_t k = 2*i+j;
810
811 // method 1
812 tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfCh_ClusterIn+ia))->ProjectionY("tmp",k+1,k+1,"e");
813 GetMean(tmp, meanIn, meanInErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMean_ClusterIn+ia), k, k+1);
814 GetRMS(tmp, sigmaIn, sigmaInErr);
815 delete tmp;
816
817 tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfCh_ClusterOut+ia))->ProjectionY("tmp",k+1,k+1,"e");
818 GetMean(tmp, meanOut, meanOutErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMean_ClusterOut+ia), k, k+1);
819 GetRMS(tmp, sigmaOut, sigmaOutErr);
820 delete tmp;
821
822 if (fCorrectForSystematics) {
823 sigmaIn = TMath::Sqrt(sigmaIn*sigmaIn + meanIn*meanIn);
824 sigmaInErr = (sigmaIn>0) ? TMath::Sqrt(sigmaIn*sigmaIn*sigmaInErr*sigmaInErr + meanIn*meanIn*meanInErr*meanInErr) / sigmaIn : 0.;
825 sigmaOut = TMath::Sqrt(sigmaOut*sigmaOut + meanOut*meanOut);
826 sigmaOutErr = (sigmaOut>0) ? TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + meanOut*meanOut*meanOutErr*meanOutErr) / sigmaOut : 0.;
827 }
828
829 clusterRes = TMath::Sqrt(sigmaIn*sigmaOut);
830 // clusterResErr = (clusterRes > 0.) ? 0.5 * TMath::Sqrt(sigmaInErr*sigmaInErr*sigmaOut*sigmaOut + sigmaIn*sigmaIn*sigmaOutErr*sigmaOutErr) / clusterRes : 0.;
831 clusterResErr = TMath::Sqrt(sigmaInErr*sigmaOutErr);
832 ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerHalfChSigma+ia))->SetPoint(k, k+1, clusterRes);
833 ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerHalfChSigma+ia))->SetPointError(k, 0., clusterResErr);
834
835 // method 2
836 tmp = ((TH2F*)fResiduals->UncheckedAt(kTrackResPerHalfCh+ia))->ProjectionY("tmp",k+1,k+1,"e");
837 GetMean(tmp, sigmaTrack, sigmaTrackErr, 0x0, 0, 0, kFALSE);
838 delete tmp;
839
840 tmp = ((TH2F*)fResiduals->UncheckedAt(kMCSPerHalfCh+ia))->ProjectionY("tmp",k+1,k+1,"e");
841 GetMean(tmp, sigmaMCS, sigmaMCSErr, 0x0, 0, 0, kFALSE);
842 delete tmp;
843
844 sigmaCluster = sigmaOut*sigmaOut - sigmaTrack*sigmaTrack;
845 if (sigmaCluster > 0.) {
846 sigmaCluster = TMath::Sqrt(sigmaCluster);
847 sigmaClusterErr = TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + sigmaTrack*sigmaTrack*sigmaTrackErr*sigmaTrackErr) / sigmaCluster;
848 } else {
849 sigmaCluster = 0.;
850 sigmaClusterErr = 0.;
851 }
852 ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerHalfCh+ia))->SetPoint(k, k+1, sigmaCluster);
853 ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerHalfCh+ia))->SetPointError(k, 0., sigmaClusterErr);
854
855 }
856
857 // compute averaged local chi2
858 tmp = ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerCh+ia))->ProjectionY("tmp",i+1,i+1,"e");
859 ((TGraphErrors*)fLocalChi2->UncheckedAt(kLocalChi2PerChMean+ia))->SetPoint(i, i+1, tmp->GetMean());
860 ((TGraphErrors*)fLocalChi2->UncheckedAt(kLocalChi2PerChMean+ia))->SetPointError(i, 0., tmp->GetMeanError());
861 delete tmp;
862
863 }
864
865 // compute residual mean and dispersion per DE
866 for (Int_t i = 0; i < fNDE; i++) {
867
868 // method 1
869 TH1D *tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerDE_ClusterIn+ia))->ProjectionY("tmp",i+1,i+1,"e");
870 GetMean(tmp, meanIn, meanInErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMean_ClusterIn+ia), i, i+1);
871 GetRMS(tmp, sigmaIn, sigmaInErr);
872 delete tmp;
873
874 tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerDE_ClusterOut+ia))->ProjectionY("tmp",i+1,i+1,"e");
875 GetMean(tmp, meanOut, meanOutErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMean_ClusterOut+ia), i, i+1);
876 GetRMS(tmp, sigmaOut, sigmaOutErr);
877 delete tmp;
878
879 if (fCorrectForSystematics) {
880 sigmaIn = TMath::Sqrt(sigmaIn*sigmaIn + meanIn*meanIn);
881 sigmaInErr = (sigmaIn>0) ? TMath::Sqrt(sigmaIn*sigmaIn*sigmaInErr*sigmaInErr + meanIn*meanIn*meanInErr*meanInErr) / sigmaIn : 0.;
882 sigmaOut = TMath::Sqrt(sigmaOut*sigmaOut + meanOut*meanOut);
883 sigmaOutErr = (sigmaOut>0) ? TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + meanOut*meanOut*meanOutErr*meanOutErr) / sigmaOut : 0.;
884 }
885
886 clusterRes = TMath::Sqrt(sigmaIn*sigmaOut);
887 // clusterResErr = (clusterRes > 0.) ? 0.5 * TMath::Sqrt(sigmaInErr*sigmaInErr*sigmaOut*sigmaOut + sigmaIn*sigmaIn*sigmaOutErr*sigmaOutErr) / clusterRes : 0.;
888 clusterResErr = TMath::Sqrt(sigmaInErr*sigmaOutErr);
889 ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerDESigma+ia))->SetPoint(i, i+1, clusterRes);
890 ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerDESigma+ia))->SetPointError(i, 0., clusterResErr);
891
892 // method 2
893 tmp = ((TH2F*)fResiduals->UncheckedAt(kTrackResPerDE+ia))->ProjectionY("tmp",i+1,i+1,"e");
894 GetMean(tmp, sigmaTrack, sigmaTrackErr, 0x0, 0, 0, kFALSE);
895 delete tmp;
896
897 tmp = ((TH2F*)fResiduals->UncheckedAt(kMCSPerDE+ia))->ProjectionY("tmp",i+1,i+1,"e");
898 GetMean(tmp, sigmaMCS, sigmaMCSErr, 0x0, 0, 0, kFALSE);
899 delete tmp;
900
901 sigmaCluster = sigmaOut*sigmaOut - sigmaTrack*sigmaTrack;
902 if (sigmaCluster > 0.) {
903 sigmaCluster = TMath::Sqrt(sigmaCluster);
904 sigmaClusterErr = TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + sigmaTrack*sigmaTrack*sigmaTrackErr*sigmaTrackErr) / sigmaCluster;
905 } else {
906 sigmaCluster = 0.;
907 sigmaClusterErr = 0.;
908 }
909 ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerDE+ia))->SetPoint(i, i+1, sigmaCluster);
910 ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerDE+ia))->SetPointError(i, 0., sigmaClusterErr);
911
912 // compute averaged local chi2
913 tmp = ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerDE+ia))->ProjectionY("tmp",i+1,i+1,"e");
914 ((TGraphErrors*)fLocalChi2->UncheckedAt(kLocalChi2PerDEMean+ia))->SetPoint(i, i+1, tmp->GetMean());
915 ((TGraphErrors*)fLocalChi2->UncheckedAt(kLocalChi2PerDEMean+ia))->SetPointError(i, 0., tmp->GetMeanError());
916 delete tmp;
917
918 }
919
920 // set graph labels
921 TAxis* xAxis = ((TH2F*)fResiduals->UncheckedAt(kResidualPerDE_ClusterOut+ia))->GetXaxis();
922 ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMean_ClusterIn+ia))->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
923 ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMean_ClusterOut+ia))->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
924 ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerDESigma+ia))->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
925 ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerDE+ia))->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
926 ((TGraphErrors*)fLocalChi2->UncheckedAt(kLocalChi2PerDEMean+ia))->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
927 for (Int_t i = 1; i <= fNDE; i++) {
928 const char* label = xAxis->GetBinLabel(i);
929 ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMean_ClusterIn+ia))->GetXaxis()->SetBinLabel(i, label);
930 ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMean_ClusterOut+ia))->GetXaxis()->SetBinLabel(i, label);
931 ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerDESigma+ia))->GetXaxis()->SetBinLabel(i, label);
932 ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerDE+ia))->GetXaxis()->SetBinLabel(i, label);
933 ((TGraphErrors*)fLocalChi2->UncheckedAt(kLocalChi2PerDEMean+ia))->GetXaxis()->SetBinLabel(i, label);
934 }
935
936 }
937
938 // compute averaged local chi2 per chamber (X+Y)
939 TH2F* h2 = (TH2F*)fResiduals->UncheckedAt(kLocalChi2PerCh+2);
940 g = new TGraphErrors(AliMUONConstants::NTrackingCh());
941 g->SetName("gLocalChi2PerChMean");
942 g->SetTitle("local chi2 per Ch: mean;chamber ID;<local #chi^{2}>");
943 g->SetMarkerStyle(kFullDotLarge);
944 fLocalChi2->AddAtAndExpand(g, kLocalChi2PerChMean+2);
945 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
946 TH1D* tmp = h2->ProjectionY("tmp",i+1,i+1,"e");
947 g->SetPoint(i, i+1, tmp->GetMean());
948 g->SetPointError(i, 0., tmp->GetMeanError());
949 delete tmp;
950 }
951
952 // compute averaged local chi2 per DE (X+Y)
953 h2 = (TH2F*)fResiduals->UncheckedAt(kLocalChi2PerDE+2);
954 g = new TGraphErrors(fNDE);
955 g->SetName("gLocalChi2PerDEMean");
956 g->SetTitle("local chi2 per DE: mean;DE ID;<local #chi^{2}>");
957 g->SetMarkerStyle(kFullDotLarge);
958 fLocalChi2->AddAtAndExpand(g, kLocalChi2PerDEMean+2);
959 for (Int_t i = 0; i < fNDE; i++) {
960 TH1D* tmp = h2->ProjectionY("tmp",i+1,i+1,"e");
961 g->SetPoint(i, i+1, tmp->GetMean());
962 g->SetPointError(i, 0., tmp->GetMeanError());
963 delete tmp;
964 }
965
966 // set graph labels
967 g->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
968 for (Int_t i = 1; i <= fNDE; i++) {
969 const char* label = h2->GetXaxis()->GetBinLabel(i);
970 g->GetXaxis()->SetBinLabel(i, label);
971 }
972
973 // display
974 fCanvases = new TObjArray(1000);
975 fCanvases->SetOwner();
976
977 TLegend *lResPerChMean = new TLegend(0.75,0.85,0.99,0.99);
978 TLegend *lResPerChSigma1 = new TLegend(0.75,0.85,0.99,0.99);
979 TLegend *lResPerChSigma2 = new TLegend(0.75,0.85,0.99,0.99);
980 TLegend *lResPerChSigma3 = new TLegend(0.75,0.85,0.99,0.99);
981
982 TCanvas* cResPerCh = new TCanvas("cResPerCh","cResPerCh",1200,500);
983 cResPerCh->Divide(4,2);
984 for (Int_t ia = 0; ia < 2; ia++) {
985 cResPerCh->cd(1+4*ia);
986 g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChMean_ClusterOut+ia);
987 g->Draw("ap");
988 g->SetMarkerColor(2);
989 g->SetLineColor(2);
990 if (ia == 0) lResPerChMean->AddEntry(g,"cluster out","PL");
991 g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChMean_ClusterIn+ia);
992 g->Draw("p");
993 g->SetMarkerColor(4);
994 g->SetLineColor(4);
995 if (ia == 0) lResPerChMean->AddEntry(g,"cluster in","PL");
996 if (ia == 0) lResPerChMean->Draw();
997 else lResPerChMean->DrawClone();
998 cResPerCh->cd(2+4*ia);
999 g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChSigma_ClusterOut+ia);
1000 g->Draw("ap");
1001 g->SetMinimum(0.);
1002 g->SetMarkerColor(2);
1003 g->SetLineColor(2);
1004 if (ia == 0) lResPerChSigma1->AddEntry(g,"cluster out","PL");
1005 g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChSigma_ClusterIn+ia);
1006 g->Draw("p");
1007 g->SetMarkerColor(4);
1008 g->SetLineColor(4);
1009 if (ia == 0) lResPerChSigma1->AddEntry(g,"cluster in","PL");
1010 g = (TGraphErrors*)fChamberRes->UncheckedAt(kMCSPerChMean+ia);
1011 g->Draw("p");
1012 g->SetMarkerColor(5);
1013 g->SetLineColor(5);
1014 if (ia == 0) lResPerChSigma1->AddEntry(g,"MCS","PL");
1015 g = (TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerChSigma+ia);
1016 g->Draw("p");
1017 g->SetMarkerColor(3);
1018 g->SetLineColor(3);
1019 if (ia == 0) lResPerChSigma1->AddEntry(g,"combined 1","PL");
1020 if (ia == 0) lResPerChSigma1->Draw();
1021 else lResPerChSigma1->DrawClone();
1022 cResPerCh->cd(3+4*ia);
1023 g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChDispersion_ClusterOut+ia);
1024 g->Draw("ap");
1025 g->SetMinimum(0.);
1026 g->SetMarkerColor(2);
1027 g->SetLineColor(2);
1028 if (ia == 0) lResPerChSigma2->AddEntry(g,"cluster out","PL");
1029 g = (TGraphErrors*)fChamberRes->UncheckedAt(kMCSPerChMean+ia);
1030 g->Draw("p");
1031 if (ia == 0) lResPerChSigma2->AddEntry(g,"MCS","PL");
1032 g = (TGraphErrors*)fChamberRes->UncheckedAt(kTrackResPerChMean+ia);
1033 g->Draw("p");
1034 g->SetMarkerColor(4);
1035 g->SetLineColor(4);
1036 if (ia == 0) lResPerChSigma2->AddEntry(g,"track res.","PL");
1037 g = (TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerCh+ia);
1038 g->Draw("p");
1039 if (ia == 0) lResPerChSigma2->AddEntry(g,"combined 2","PL");
1040 if (ia == 0) lResPerChSigma2->Draw();
1041 else lResPerChSigma2->DrawClone();
1042 cResPerCh->cd(4+4*ia);
1043 g = (TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerChSigma+ia);
1044 g->Draw("ap");
1045 g->SetMinimum(0.);
1046 if (ia == 0) lResPerChSigma3->AddEntry(g,"combined 1","PL");
1047 g = (TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerCh+ia);
1048 g->Draw("p");
1049 if (ia == 0) lResPerChSigma3->AddEntry(g,"combined 2","PL");
1050 if (ia == 0) lResPerChSigma3->Draw();
1051 else lResPerChSigma3->DrawClone();
1052 }
1053 fCanvases->AddAtAndExpand(cResPerCh, kResPerCh);
1054
1055 TCanvas* cResPerHalfCh = new TCanvas("cResPerHalfCh","cResPerHalfCh",1200,500);
1056 cResPerHalfCh->Divide(2,2);
1057 for (Int_t ia = 0; ia < 2; ia++) {
1058 cResPerHalfCh->cd(1+2*ia);
1059 g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMean_ClusterOut+ia);
1060 g->Draw("ap");
1061 g->SetMarkerColor(2);
1062 g->SetLineColor(2);
1063 g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMean_ClusterIn+ia);
1064 g->Draw("p");
1065 g->SetMarkerColor(4);
1066 g->SetLineColor(4);
1067 lResPerChMean->DrawClone();
1068 cResPerHalfCh->cd(2+2*ia);
1069 g = (TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerHalfChSigma+ia);
1070 g->Draw("ap");
1071 g->SetMinimum(0.);
1072 g->SetMarkerColor(3);
1073 g->SetLineColor(3);
1074 g = (TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerHalfCh+ia);
1075 g->Draw("p");
1076 lResPerChSigma3->DrawClone();
1077 }
1078 fCanvases->AddAtAndExpand(cResPerHalfCh, kResPerHalfCh);
1079
1080 TCanvas* cResPerDE = new TCanvas("cResPerDE","cResPerDE",1200,800);
1081 cResPerDE->Divide(1,4);
1082 for (Int_t ia = 0; ia < 2; ia++) {
1083 cResPerDE->cd(1+ia);
1084 g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMean_ClusterOut+ia);
1085 g->Draw("ap");
1086 g->SetMarkerColor(2);
1087 g->SetLineColor(2);
1088 g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMean_ClusterIn+ia);
1089 g->Draw("p");
1090 g->SetMarkerColor(4);
1091 g->SetLineColor(4);
1092 lResPerChMean->DrawClone();
1093 cResPerDE->cd(3+ia);
1094 g = (TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerDESigma+ia);
1095 g->Draw("ap");
1096 g->SetMinimum(0.);
1097 g->SetMarkerColor(3);
1098 g->SetLineColor(3);
1099 g = (TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerDE+ia);
1100 g->Draw("p");
1101 lResPerChSigma3->DrawClone();
1102 }
1103 fCanvases->AddAtAndExpand(cResPerDE, kResPerDE);
1104
1105 TCanvas* cResPerChVsP = new TCanvas("cResPerChVsP","cResPerChVsP");
1106 cResPerChVsP->Divide(1,2);
1107 for (Int_t ia = 0; ia < 2; ia++) {
1108 cResPerChVsP->cd(1+ia);
1109 mg = (TMultiGraph*)fChamberRes->UncheckedAt(kCombinedResidualSigmaVsP+ia);
1110 mg->Draw("ap");
1111 }
1112 fCanvases->AddAtAndExpand(cResPerChVsP, kResPerChVsP);
1113
1114 // print results
1115 printf("\nchamber resolution:\n");
1116 printf(" - non-bending:");
1117 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %5.3f":", %5.3f",newClusterRes[0][i]);
1118 printf("\n - bending:");
1119 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %6.4f":", %6.4f",newClusterRes[1][i]);
1120 printf("\n\n");
1121
1122 // Post final data.
1123 PostData(3, fLocalChi2);
1124 PostData(4, fChamberRes);
1125}
1126
1127//________________________________________________________________________
1128void AliAnalysisTaskMuonResolution::ModifyClusters(AliMUONTrack& track)
1129{
1130 /// Reset the clusters resolution from the ones given to the task and change
1131 /// the cluster position according to the new alignment parameters if required
1132
1133 Double_t gX,gY,gZ,lX,lY,lZ;
1134
1135 // loop over clusters
1136 Int_t nClusters = track.GetNClusters();
1137 for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
1138
1139 AliMUONVCluster* cl = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->UncheckedAt(iCluster))->GetClusterPtr();
1140
1141 // change their resolution
1142 cl->SetErrXY(fClusterResNB[cl->GetChamberId()], fClusterResB[cl->GetChamberId()]);
1143
1144 // change their position
1145 if (fReAlign) {
1146 gX = cl->GetX();
1147 gY = cl->GetY();
1148 gZ = cl->GetZ();
1149 fOldGeoTransformer->Global2Local(cl->GetDetElemId(),gX,gY,gZ,lX,lY,lZ);
1150 fNewGeoTransformer->Local2Global(cl->GetDetElemId(),lX,lY,lZ,gX,gY,gZ);
1151 cl->SetXYZ(gX,gY,gZ);
1152 }
1153
1154 }
1155
1156}
1157
1158//________________________________________________________________________
1159void AliAnalysisTaskMuonResolution::Zoom(TH1* h, Double_t fractionCut)
1160{
1161 /// Reduce the range of the histogram by removing a given fration of the statistic at each edge
1162 ZoomLeft(h, fractionCut);
1163 ZoomRight(h, fractionCut);
1164}
1165
1166//________________________________________________________________________
1167void AliAnalysisTaskMuonResolution::ZoomLeft(TH1* h, Double_t fractionCut)
1168{
1169 /// Reduce the range of the histogram by removing a given fration of the statistic on the left side
1170 Int_t maxEventsCut = (Int_t) (fractionCut * h->GetEntries());
1171 Int_t nBins = h->GetNbinsX();
1172
1173 // set low edge
1174 Int_t minBin;
1175 Int_t eventsCut = 0;
1176 for (minBin = 1; minBin <= nBins; minBin++) {
1177 eventsCut += (Int_t) h->GetBinContent(minBin);
1178 if (eventsCut > maxEventsCut) break;
1179 }
1180
1181 // set new axis range
1182 h->GetXaxis()->SetRange(--minBin, h->GetXaxis()->GetLast());
1183}
1184
1185//________________________________________________________________________
1186void AliAnalysisTaskMuonResolution::ZoomRight(TH1* h, Double_t fractionCut)
1187{
1188 /// Reduce the range of the histogram by removing a given fration of the statistic on the right side
1189 Int_t maxEventsCut = (Int_t) (fractionCut * h->GetEntries());
1190 Int_t nBins = h->GetNbinsX();
1191
1192 // set high edge
1193 Int_t maxBin;
1194 Int_t eventsCut = 0;
1195 for (maxBin = nBins; maxBin >= 1; maxBin--) {
1196 eventsCut += (Int_t) h->GetBinContent(maxBin);
1197 if (eventsCut > maxEventsCut) break;
1198 }
1199
1200 // set new axis range
1201 h->GetXaxis()->SetRange(h->GetXaxis()->GetFirst(), ++maxBin);
1202}
1203
1204//________________________________________________________________________
1205void AliAnalysisTaskMuonResolution::GetMean(TH1* h, Double_t& mean, Double_t& meanErr, TGraphErrors* g, Int_t i, Double_t x, Bool_t zoom)
1206{
1207 /// Fill graph with the mean value of the histogram and the corresponding error (zooming if required)
1208 Int_t firstBin = h->GetXaxis()->GetFirst();
1209 Int_t lastBin = h->GetXaxis()->GetLast();
1210 if (zoom) Zoom(h);
1211 mean = (h->GetEntries() > fgkMinEntries) ? h->GetMean() : 0.;
1212 meanErr = (h->GetEntries() > fgkMinEntries) ? h->GetMeanError() : 0.;
1213 if (g) {
1214 g->SetPoint(i, x, mean);
1215 g->SetPointError(i, 0., meanErr);
1216 }
1217 if (zoom) h->GetXaxis()->SetRange(firstBin,lastBin);
1218}
1219
1220//________________________________________________________________________
1221void AliAnalysisTaskMuonResolution::GetRMS(TH1* h, Double_t& rms, Double_t& rmsErr, TGraphErrors* g, Int_t i, Double_t x, Bool_t zoom)
1222{
1223 /// Return the RMS of the histogram and the corresponding error (zooming if required) and fill graph if !=0x0
1224 Int_t firstBin = h->GetXaxis()->GetFirst();
1225 Int_t lastBin = h->GetXaxis()->GetLast();
1226 if (zoom) Zoom(h);
1227 rms = (h->GetEntries() > fgkMinEntries) ? h->GetRMS() : 0.;
1228 rmsErr = (h->GetEntries() > fgkMinEntries) ? h->GetRMSError() : 0.;
1229 if (g) {
1230 g->SetPoint(i, x, rms);
1231 g->SetPointError(i, 0., rmsErr);
1232 }
1233 if (zoom) h->GetXaxis()->SetRange(firstBin,lastBin);
1234}
1235
1236//________________________________________________________________________
1237void AliAnalysisTaskMuonResolution::FillSigmaClusterVsP(TH2* hIn, TH2* hOut, TGraphErrors* g, Bool_t zoom)
1238{
1239 /// Fill graph with cluster resolution from combined residuals with cluster in/out (zooming if required)
1240 Double_t sigmaIn, sigmaInErr, sigmaOut, sigmaOutErr, clusterRes, clusterResErr;
1241 for (Int_t j = 1; j <= hIn->GetNbinsX(); j++) {
1242 TH1D* tmp = hIn->ProjectionY("tmp",j,j,"e");
1243 GetRMS(tmp, sigmaIn, sigmaInErr, 0x0, 0, 0., zoom);
1244 delete tmp;
1245 tmp = hOut->ProjectionY("tmp",j,j,"e");
1246 GetRMS(tmp, sigmaOut, sigmaOutErr, 0x0, 0, 0., zoom);
1247 delete tmp;
1248 Double_t p = 0.5 * (hIn->GetBinLowEdge(j) + hIn->GetBinLowEdge(j+1));
1249 Double_t pErr = p - hIn->GetBinLowEdge(j);
1250 clusterRes = TMath::Sqrt(sigmaIn*sigmaOut);
1251 clusterResErr = (clusterRes > 0.) ? 0.5 * TMath::Sqrt(sigmaInErr*sigmaInErr*sigmaOut*sigmaOut + sigmaIn*sigmaIn*sigmaOutErr*sigmaOutErr) / clusterRes : 0.;
1252 g->SetPoint(j, p, clusterRes);
1253 g->SetPointError(j, pErr, clusterResErr);
1254 }
1255}
1256
1257//__________________________________________________________________________
1258void AliAnalysisTaskMuonResolution::Cov2CovP(const AliMUONTrackParam &param, TMatrixD &covP)
1259{
1260 /// change coordinate system: (X, SlopeX, Y, SlopeY, q/Pyz) -> (X, Y, pX, pY, pZ)
1261 /// parameters (param) are given in the (X, SlopeX, Y, SlopeY, q/Pyz) coordinate system
1262
1263 // Get useful parameters
1264 Double_t slopeX = param.GetNonBendingSlope();
1265 Double_t slopeY = param.GetBendingSlope();
1266 Double_t qOverPYZ = param.GetInverseBendingMomentum();
1267 Double_t pZ = param.Pz();
1268
1269 // compute Jacobian to change the coordinate system from (X,SlopeX,Y,SlopeY,c/pYZ) to (X,Y,pX,pY,pZ)
1270 Double_t dpZdSlopeY = - qOverPYZ * qOverPYZ * pZ * pZ * pZ * slopeY;
1271 Double_t dpZdQOverPYZ = (qOverPYZ != 0.) ? - pZ / qOverPYZ : - FLT_MAX;
1272 TMatrixD jacob(5,5);
1273 jacob.Zero();
1274 jacob(0,0) = 1.;
1275 jacob(1,2) = 1.;
1276 jacob(2,1) = pZ;
1277 jacob(2,3) = slopeX * dpZdSlopeY;
1278 jacob(2,4) = slopeX * dpZdQOverPYZ;
1279 jacob(3,3) = pZ + slopeY * dpZdSlopeY;
1280 jacob(3,4) = slopeY * dpZdQOverPYZ;
1281 jacob(4,3) = dpZdSlopeY;
1282 jacob(4,4) = dpZdQOverPYZ;
1283
1284 // compute covariances in new coordinate system
1285 TMatrixD tmp(param.GetCovariances(),TMatrixD::kMultTranspose,jacob);
1286 covP.Mult(jacob,tmp);
1287}
1288