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