1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
23 #include <TMultiGraph.h>
24 #include <TGraphErrors.h>
27 #include <Riostream.h>
29 #include <TGeoManager.h>
31 #include <TObjString.h>
35 #include "AliESDEvent.h"
36 #include "AliESDMuonTrack.h"
37 #include "AliCDBManager.h"
38 #include "AliCDBStorage.h"
39 #include "AliGeomManager.h"
42 #include "AliAnalysisDataSlot.h"
43 #include "AliAnalysisManager.h"
44 #include "AliInputEventHandler.h"
45 #include "AliAnalysisTaskMuonResolution.h"
46 #include "AliCentrality.h"
49 #include "AliMUONCDB.h"
50 #include "AliMUONRecoParam.h"
51 #include "AliMUONESDInterface.h"
52 #include "AliMUONVTrackReconstructor.h"
53 #include "AliMUONTrack.h"
54 #include "AliMUONTrackParam.h"
55 #include "AliMUONTrackExtrap.h"
56 #include "AliMUONVCluster.h"
57 #include "AliMUONVDigit.h"
58 #include "AliMUONGeometryTransformer.h"
59 #include "AliMUONGeometryModuleTransformer.h"
60 #include "AliMUONGeometryDetElement.h"
61 #include "AliMpDEIterator.h"
62 #include "AliMpSegmentation.h"
63 #include "AliMpVSegmentation.h"
64 #include "AliMpConstants.h"
65 #include "AliMpDDLStore.h"
67 #include "AliMpDetElement.h"
68 #include "AliMpCathodType.h"
71 #define SafeDelete(x) if (x != NULL) { delete x; x = NULL; }
74 ClassImp(AliAnalysisTaskMuonResolution)
76 const Int_t AliAnalysisTaskMuonResolution::fgkMinEntries = 10;
78 //________________________________________________________________________
79 AliAnalysisTaskMuonResolution::AliAnalysisTaskMuonResolution() :
83 fResidualsVsCent(NULL),
90 fShowProgressBar(kFALSE),
91 fPrintClResPerCh(kFALSE),
92 fPrintClResPerDE(kFALSE),
95 fSelectPhysics(kFALSE),
98 fSelectTrigger(kFALSE),
101 fCorrectForSystematics(kTRUE),
102 fRemoveMonoCathCl(kFALSE),
103 fCheckAllPads(kFALSE),
107 fOldAlignStorage(""),
108 fNewAlignStorage(""),
109 fOldGeoTransformer(NULL),
110 fNewGeoTransformer(NULL),
111 fSelectTriggerClass(NULL)
113 /// Default constructor
116 //________________________________________________________________________
117 AliAnalysisTaskMuonResolution::AliAnalysisTaskMuonResolution(const char *name) :
118 AliAnalysisTaskSE(name),
121 fResidualsVsCent(NULL),
126 fDefaultStorage("raw://"),
128 fShowProgressBar(kFALSE),
129 fPrintClResPerCh(kFALSE),
130 fPrintClResPerDE(kFALSE),
133 fSelectPhysics(kFALSE),
135 fApplyAccCut(kFALSE),
136 fSelectTrigger(kFALSE),
139 fCorrectForSystematics(kTRUE),
140 fRemoveMonoCathCl(kFALSE),
141 fCheckAllPads(kFALSE),
145 fOldAlignStorage(""),
146 fNewAlignStorage(""),
147 fOldGeoTransformer(NULL),
148 fNewGeoTransformer(NULL),
149 fSelectTriggerClass(NULL)
153 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) SetStartingResolution(i, -1., -1.);
157 // Output slot #1 writes into a TObjArray container
158 DefineOutput(1,TObjArray::Class());
159 // Output slot #2 writes into a TObjArray container
160 DefineOutput(2,TObjArray::Class());
161 // Output slot #3 writes into a TObjArray container
162 DefineOutput(3,TObjArray::Class());
163 // Output slot #4 writes into a TObjArray container
164 DefineOutput(4,TObjArray::Class());
165 // Output slot #5 writes into a TObjArray container
166 DefineOutput(5,TObjArray::Class());
167 // Output slot #5 writes into a TObjArray container
168 DefineOutput(6,TObjArray::Class());
171 //________________________________________________________________________
172 AliAnalysisTaskMuonResolution::~AliAnalysisTaskMuonResolution()
175 if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
176 SafeDelete(fResiduals);
177 SafeDelete(fResidualsVsP);
178 SafeDelete(fResidualsVsCent);
179 SafeDelete(fTrackRes);
181 SafeDelete(fChamberRes);
182 SafeDelete(fCanvases);
184 SafeDelete(fOldGeoTransformer);
185 SafeDelete(fNewGeoTransformer);
186 SafeDelete(fSelectTriggerClass);
189 //___________________________________________________________________________
190 void AliAnalysisTaskMuonResolution::UserCreateOutputObjects()
192 /// Create histograms
194 // do it once the OCDB has been loaded (i.e. from NotifyRun())
195 if (!fOCDBLoaded) return;
197 // set the list of trigger classes that can be selected to fill histograms (in case the physics selection is not used)
198 fSelectTriggerClass = new TList();
199 fSelectTriggerClass->SetOwner();
200 fSelectTriggerClass->AddLast(new TObjString(" CINT1B-ABCE-NOPF-ALL ")); fSelectTriggerClass->Last()->SetUniqueID(AliVEvent::kMB);
201 fSelectTriggerClass->AddLast(new TObjString(" CMUS1B-ABCE-NOPF-MUON ")); fSelectTriggerClass->Last()->SetUniqueID(AliVEvent::kMUON);
202 fSelectTriggerClass->AddLast(new TObjString(" CINT1-B-")); fSelectTriggerClass->Last()->SetUniqueID(AliVEvent::kMB);
203 fSelectTriggerClass->AddLast(new TObjString(" CMUS1-B-")); fSelectTriggerClass->Last()->SetUniqueID(AliVEvent::kMUON);
204 fSelectTriggerClass->AddLast(new TObjString(" CSH1-B-")); fSelectTriggerClass->Last()->SetUniqueID(AliVEvent::kHighMult);
206 fResiduals = new TObjArray(1000);
207 fResiduals->SetOwner();
208 fResidualsVsP = new TObjArray(1000);
209 fResidualsVsP->SetOwner();
210 fResidualsVsCent = new TObjArray(1000);
211 fResidualsVsCent->SetOwner();
212 fTrackRes = new TObjArray(1000);
213 fTrackRes->SetOwner();
216 // find the highest chamber resolution and set histogram bins
217 const AliMUONRecoParam* recoParam = AliMUONESDInterface::GetTracker()->GetRecoParam();
218 Double_t maxSigma[2] = {-1., -1.};
219 for (Int_t i = 0; i < 10; i++) {
220 if (recoParam->GetDefaultNonBendingReso(i) > maxSigma[0]) maxSigma[0] = recoParam->GetDefaultNonBendingReso(i);
221 if (recoParam->GetDefaultBendingReso(i) > maxSigma[1]) maxSigma[1] = recoParam->GetDefaultBendingReso(i);
223 const char* axes[2] = {"X", "Y"};
224 const Int_t nBins = 5000;
225 const Int_t nSigma = 10;
226 const Int_t pNBins = 20;
227 const Double_t pEdges[2] = {0., 50.};
228 Int_t nCentBins = 12;
229 Double_t centRange[2] = {-10., 110.};
231 for (Int_t ia = 0; ia < 2; ia++) {
233 Double_t maxRes = nSigma*maxSigma[ia];
235 // List of residual histos
236 h2 = new TH2F(Form("hResidual%sPerCh_ClusterIn",axes[ia]), Form("cluster-track residual-%s distribution per chamber (cluster attached to the track);chamber ID;#Delta_{%s} (cm)",axes[ia],axes[ia]), 10, 0.5, 10.5, nBins, -maxRes, maxRes);
237 fResiduals->AddAtAndExpand(h2, kResidualPerChClusterIn+ia);
238 h2 = new TH2F(Form("hResidual%sPerCh_ClusterOut",axes[ia]), Form("cluster-track residual-%s distribution per chamber (cluster not attached to the track);chamber ID;#Delta_{%s} (cm)",axes[ia],axes[ia]), 10, 0.5, 10.5, nBins, -2.*maxRes, 2.*maxRes);
239 fResiduals->AddAtAndExpand(h2, kResidualPerChClusterOut+ia);
241 h2 = new TH2F(Form("hResidual%sPerHalfCh_ClusterIn",axes[ia]), Form("cluster-track residual-%s distribution per half chamber (cluster attached to the track);half chamber ID;#Delta_{%s} (cm)",axes[ia],axes[ia]), 20, 0.5, 20.5, nBins, -maxRes, maxRes);
242 for (Int_t i = 0; i < 10; i++) { h2->GetXaxis()->SetBinLabel(2*i+1, Form("%d-I",i+1)); h2->GetXaxis()->SetBinLabel(2*i+2, Form("%d-O",i+1)); }
243 fResiduals->AddAtAndExpand(h2, kResidualPerHalfChClusterIn+ia);
244 h2 = new TH2F(Form("hResidual%sPerHalfCh_ClusterOut",axes[ia]), Form("cluster-track residual-%s distribution per half chamber (cluster not attached to the track);half chamber ID;#Delta_{%s} (cm)",axes[ia],axes[ia]), 20, 0.5, 20.5, nBins, -2.*maxRes, 2.*maxRes);
245 for (Int_t i = 0; i < 10; i++) { h2->GetXaxis()->SetBinLabel(2*i+1, Form("%d-I",i+1)); h2->GetXaxis()->SetBinLabel(2*i+2, Form("%d-O",i+1)); }
246 fResiduals->AddAtAndExpand(h2, kResidualPerHalfChClusterOut+ia);
248 h2 = new TH2F(Form("hResidual%sPerDE_ClusterIn",axes[ia]), Form("cluster-track residual-%s distribution per DE (cluster attached to the track);DE ID;#Delta_{%s} (cm)",axes[ia],axes[ia]), fNDE, 0.5, fNDE+0.5, nBins, -maxRes, maxRes);
249 for (Int_t i = 1; i <= fNDE; i++) h2->GetXaxis()->SetBinLabel(i, Form("%d",fDEIds[i]));
250 fResiduals->AddAtAndExpand(h2, kResidualPerDEClusterIn+ia);
251 h2 = new TH2F(Form("hResidual%sPerDE_ClusterOut",axes[ia]), Form("cluster-track residual-%s distribution per DE (cluster not attached to the track);DE ID;#Delta_{%s} (cm)",axes[ia],axes[ia]), fNDE, 0.5, fNDE+0.5, nBins, -2.*maxRes, 2.*maxRes);
252 for (Int_t i = 1; i <= fNDE; i++) h2->GetXaxis()->SetBinLabel(i, Form("%d",fDEIds[i]));
253 fResiduals->AddAtAndExpand(h2, kResidualPerDEClusterOut+ia);
255 h2 = new TH2F(Form("hTrackRes%sPerCh",axes[ia]), Form("track #sigma_{%s} per Ch;chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]), 10, 0.5, 10.5, nBins, 0., maxRes);
256 fResiduals->AddAtAndExpand(h2, kTrackResPerCh+ia);
257 h2 = new TH2F(Form("hTrackRes%sPerHalfCh",axes[ia]), Form("track #sigma_{%s} per half Ch;half chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]), 20, 0.5, 20.5, nBins, 0., maxRes);
258 for (Int_t i = 0; i < 10; i++) { h2->GetXaxis()->SetBinLabel(2*i+1, Form("%d-I",i+1)); h2->GetXaxis()->SetBinLabel(2*i+2, Form("%d-O",i+1)); }
259 fResiduals->AddAtAndExpand(h2, kTrackResPerHalfCh+ia);
260 h2 = new TH2F(Form("hTrackRes%sPerDE",axes[ia]), Form("track #sigma_{%s} per DE;DE ID;#sigma_{%s} (cm)",axes[ia],axes[ia]), fNDE, 0.5, fNDE+0.5, nBins, 0., maxRes);
261 for (Int_t i = 1; i <= fNDE; i++) h2->GetXaxis()->SetBinLabel(i, Form("%d",fDEIds[i]));
262 fResiduals->AddAtAndExpand(h2, kTrackResPerDE+ia);
264 h2 = new TH2F(Form("hMCS%sPerCh",axes[ia]), Form("MCS %s-dispersion of extrapolated track per Ch;chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]), 10, 0.5, 10.5, nBins, 0., 0.2);
265 fResiduals->AddAtAndExpand(h2, kMCSPerCh+ia);
266 h2 = new TH2F(Form("hMCS%sPerHalfCh",axes[ia]), Form("MCS %s-dispersion of extrapolated track per half Ch;half chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]), 20, 0.5, 20.5, nBins, 0., 0.2);
267 for (Int_t i = 0; i < 10; i++) { h2->GetXaxis()->SetBinLabel(2*i+1, Form("%d-I",i+1)); h2->GetXaxis()->SetBinLabel(2*i+2, Form("%d-O",i+1)); }
268 fResiduals->AddAtAndExpand(h2, kMCSPerHalfCh+ia);
269 h2 = new TH2F(Form("hMCS%sPerDE",axes[ia]), Form("MCS %s-dispersion of extrapolated track per DE;DE ID;#sigma_{%s} (cm)",axes[ia],axes[ia]), fNDE, 0.5, fNDE+0.5, nBins, 0., 0.2);
270 for (Int_t i = 1; i <= fNDE; i++) h2->GetXaxis()->SetBinLabel(i, Form("%d",fDEIds[i]));
271 fResiduals->AddAtAndExpand(h2, kMCSPerDE+ia);
273 h2 = new TH2F(Form("hClusterRes2%sPerCh",axes[ia]), Form("cluster #sigma_{%s}^{2} per Ch;chamber ID;#sigma_{%s}^{2} (cm^{2})",axes[ia],axes[ia]), 10, 0.5, 10.5, nSigma*nBins, -0.1*maxRes*maxRes, maxRes*maxRes);
274 fResiduals->AddAtAndExpand(h2, kClusterRes2PerCh+ia);
276 // List of residual vs. p or vs. centrality histos
277 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
278 h2 = new TH2F(Form("hResidual%sInCh%dVsP_ClusterIn",axes[ia],i+1), Form("cluster-track residual-%s distribution in chamber %d versus momentum (cluster attached to the track);p (GeV/c^{2});#Delta_{%s} (cm)",axes[ia],i+1,axes[ia]), pNBins, pEdges[0], pEdges[1], nBins, -maxRes, maxRes);
279 fResidualsVsP->AddAtAndExpand(h2, kResidualInChVsPClusterIn+10*ia+i);
280 h2 = new TH2F(Form("hResidual%sInCh%dVsP_ClusterOut",axes[ia],i+1), Form("cluster-track residual-%s distribution in chamber %d versus momentum (cluster not attached to the track);p (GeV/c^{2});#Delta_{%s} (cm)",axes[ia],i+1,axes[ia]), pNBins, pEdges[0], pEdges[1], nBins, -2.*maxRes, 2.*maxRes);
281 fResidualsVsP->AddAtAndExpand(h2, kResidualInChVsPClusterOut+10*ia+i);
283 h2 = new TH2F(Form("hResidual%sInCh%dVsCent_ClusterIn",axes[ia],i+1), Form("cluster-track residual-%s distribution in chamber %d versus centrality (cluster attached to the track);centrality (%%);#Delta_{%s} (cm)",axes[ia],i+1,axes[ia]), nCentBins, centRange[0], centRange[1], nBins, -maxRes, maxRes);
284 fResidualsVsCent->AddAtAndExpand(h2, kResidualInChVsCentClusterIn+10*ia+i);
285 h2 = new TH2F(Form("hResidual%sInCh%dVsCent_ClusterOut",axes[ia],i+1), Form("cluster-track residual-%s distribution in chamber %d versus centrality (cluster not attached to the track);centrality (%%);#Delta_{%s} (cm)",axes[ia],i+1,axes[ia]), nCentBins, centRange[0], centRange[1], nBins, -2.*maxRes, 2.*maxRes);
286 fResidualsVsCent->AddAtAndExpand(h2, kResidualInChVsCentClusterOut+10*ia+i);
289 h2 = new TH2F(Form("hResidual%sVsP_ClusterIn",axes[ia]), Form("cluster-track residual-%s distribution integrated over chambers versus momentum (cluster attached to the track);p (GeV/c^{2});#Delta_{%s} (cm)",axes[ia],axes[ia]), pNBins, pEdges[0], pEdges[1], nBins, -maxRes, maxRes);
290 fResidualsVsP->AddAtAndExpand(h2, kResidualVsPClusterIn+ia);
291 h2 = new TH2F(Form("hResidual%sVsP_ClusterOut",axes[ia]), Form("cluster-track residual-%s distribution integrated over chambers versus momentum (cluster not attached to the track);p (GeV/c^{2});#Delta_{%s} (cm)",axes[ia],axes[ia]), pNBins, pEdges[0], pEdges[1], nBins, -2.*maxRes, 2.*maxRes);
292 fResidualsVsP->AddAtAndExpand(h2, kResidualVsPClusterOut+ia);
294 h2 = new TH2F(Form("hResidual%sVsCent_ClusterIn",axes[ia]), Form("cluster-track residual-%s distribution integrated over chambers versus centrality (cluster attached to the track);centrality (%%);#Delta_{%s} (cm)",axes[ia],axes[ia]), nCentBins, centRange[0], centRange[1], nBins, -maxRes, maxRes);
295 fResidualsVsCent->AddAtAndExpand(h2, kResidualVsCentClusterIn+ia);
296 h2 = new TH2F(Form("hResidual%sVsCent_ClusterOut",axes[ia]), Form("cluster-track residual-%s distribution integrated over chambers versus centrality (cluster not attached to the track);centrality (%%);#Delta_{%s} (cm)",axes[ia],axes[ia]), nCentBins, centRange[0], centRange[1], nBins, -2.*maxRes, 2.*maxRes);
297 fResidualsVsCent->AddAtAndExpand(h2, kResidualVsCentClusterOut+ia);
300 h2 = new TH2F(Form("hLocalChi2%sPerCh",axes[ia]), Form("local chi2-%s distribution per chamber;chamber ID;local #chi^{2}_{%s}", axes[ia], axes[ia]), 10, 0.5, 10.5, 1000, 0., 25.);
301 fResiduals->AddAtAndExpand(h2, kLocalChi2PerCh+ia);
302 h2 = new TH2F(Form("hLocalChi2%sPerDE",axes[ia]), Form("local chi2-%s distribution per DE;DE ID;local #chi^{2}_{%s}", axes[ia], axes[ia]), fNDE, 0.5, fNDE+0.5, 1000, 0., 25.);
303 for (Int_t i = 1; i <= fNDE; i++) h2->GetXaxis()->SetBinLabel(i, Form("%d",fDEIds[i]));
304 fResiduals->AddAtAndExpand(h2, kLocalChi2PerDE+ia);
307 h2 = new TH2F(Form("hUncorrSlope%sRes",axes[ia]), Form("muon slope_{%s} reconstructed resolution at first cluster vs p;p (GeV/c); #sigma_{slope_{%s}}", axes[ia], axes[ia]), 300, 0., 300., 1000, 0., 0.003);
308 fTrackRes->AddAtAndExpand(h2, kUncorrSlopeRes+ia);
309 h2 = new TH2F(Form("hSlope%sRes",axes[ia]), Form("muon slope_{%s} reconstructed resolution at vertex vs p;p (GeV/c); #sigma_{slope_{%s}}", axes[ia], axes[ia]), 300, 0., 300., 1000, 0., 0.02);
310 fTrackRes->AddAtAndExpand(h2, kSlopeRes+ia);
314 h2 = new TH2F("hLocalChi2PerCh", "local chi2 (~0.5*(#chi^{2}_{X}+#chi^{2}_{Y})) distribution per chamber;chamber ID;local #chi^{2}", 10, 0.5, 10.5, 1000, 0., 25.);
315 fResiduals->AddAtAndExpand(h2, kLocalChi2PerCh+2);
316 h2 = new TH2F("hLocalChi2PerDE", "local chi2 (~0.5*(#chi^{2}_{X}+#chi^{2}_{Y})) distribution per chamber;DE ID;local #chi^{2}", fNDE, 0.5, fNDE+0.5, 1000, 0., 25.);
317 for (Int_t i = 1; i <= fNDE; i++) h2->GetXaxis()->SetBinLabel(i, Form("%d",fDEIds[i]));
318 fResiduals->AddAtAndExpand(h2, kLocalChi2PerDE+2);
321 h2 = new TH2F("hUncorrPRes", "muon momentum reconstructed resolution at first cluster vs p;p (GeV/c); #sigma_{p}/p (%)", 300, 0., 300., 1000, 0., 10.);
322 fTrackRes->AddAtAndExpand(h2, kUncorrPRes);
323 h2 = new TH2F("hPRes", "muon momentum reconstructed resolution at vertex vs p;p (GeV/c); #sigma_{p}/p (%)", 300, 0., 300., 1000, 0., 10.);
324 fTrackRes->AddAtAndExpand(h2, kPRes);
325 h2 = new TH2F("hUncorrPtRes", "muon transverse momentum reconstructed resolution at first cluster vs p_{t};p_{t} (GeV/c); #sigma_{p_{t}}/p_{t} (%)", 300, 0., 30., 300, 0., 30.);
326 fTrackRes->AddAtAndExpand(h2, kUncorrPtRes);
327 h2 = new TH2F("hPtRes", "muon transverse momentum reconstructed resolution at vertex vs p_{t};p_{t} (GeV/c); #sigma_{p_{t}}/p_{t} (%)", 300, 0., 30., 300, 0., 30.);
328 fTrackRes->AddAtAndExpand(h2, kPtRes);
330 // Post data at least once per task to ensure data synchronisation (required for merging)
331 PostData(1, fResiduals);
332 PostData(2, fResidualsVsP);
333 PostData(3, fResidualsVsCent);
334 PostData(6, fTrackRes);
337 //________________________________________________________________________
338 void AliAnalysisTaskMuonResolution::UserExec(Option_t *)
342 // check if OCDB properly loaded
343 if (!fOCDBLoaded) return;
345 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
348 if (fShowProgressBar && (++fNEvents)%100 == 0) cout<<"\rEvent processing... "<<fNEvents<<"\r"<<flush;
350 // skip events that do not pass the physics selection if required
351 UInt_t triggerWord = (fInputHandler) ? fInputHandler->IsEventSelected() : 0;
352 if (fSelectPhysics && triggerWord == 0) return;
354 // skip events that do not pass the trigger selection if required
355 TString firedTriggerClasses = esd->GetFiredTriggerClasses();
356 if (!fSelectPhysics) triggerWord = BuildTriggerWord(firedTriggerClasses);
357 if (fSelectTrigger && (triggerWord & fTriggerMask) == 0) return;
359 // get the centrality
360 Float_t centrality = esd->GetCentrality()->GetCentralityPercentileUnchecked("V0M");
362 // get tracker to refit
363 AliMUONVTrackReconstructor* tracker = AliMUONESDInterface::GetTracker();
366 Int_t nTracks = (Int_t) esd->GetNumberOfMuonTracks();
367 for (Int_t iTrack = 0; iTrack < nTracks; ++iTrack) {
370 AliESDMuonTrack* esdTrack = esd->GetMuonTrack(iTrack);
373 if (!esdTrack->ContainTrackerData()) continue;
375 // skip tracks not matched with trigger if required
376 if (fMatchTrig && !esdTrack->ContainTriggerData()) continue;
378 // skip tracks that do not pass the acceptance cuts if required
379 Double_t thetaAbs = TMath::ATan(esdTrack->GetRAtAbsorberEnd()/505.) * TMath::RadToDeg();
380 Double_t eta = esdTrack->Eta();
381 if (fApplyAccCut && (thetaAbs < 2. || thetaAbs > 10. || eta < -4. || eta > -2.5)) continue;
383 // skip low momentum tracks
384 if (esdTrack->PUncorrected() < fMinMomentum) continue;
386 // get the corresponding MUON track
388 AliMUONESDInterface::ESDToMUON(*esdTrack, track, kFALSE);
390 // change the cluster resolution (and position)
391 ModifyClusters(track);
394 if (!tracker->RefitTrack(track, kFALSE)) break;
396 // save track unchanged
397 AliMUONTrack referenceTrack(track);
399 // get track param at first cluster and add MCS in first chamber
400 AliMUONTrackParam trackParamAtFirstCluster(*(static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->First())));
401 Int_t firstCh = 0; while (firstCh < 10 && !esdTrack->IsInMuonClusterMap(firstCh)) firstCh++;
402 AliMUONTrackExtrap::AddMCSEffect(&trackParamAtFirstCluster, AliMUONConstants::ChamberThicknessInX0(firstCh)/2., -1.);
404 // fill momentum error at first cluster
405 Double_t pXUncorr = trackParamAtFirstCluster.Px();
406 Double_t pYUncorr = trackParamAtFirstCluster.Py();
407 Double_t pZUncorr = trackParamAtFirstCluster.Pz();
408 Double_t pUncorr = trackParamAtFirstCluster.P();
409 TMatrixD covUncorr(5,5);
410 Cov2CovP(trackParamAtFirstCluster,covUncorr);
411 Double_t sigmaPUncorr = TMath::Sqrt(pXUncorr * (pXUncorr*covUncorr(2,2) + pYUncorr*covUncorr(2,3) + pZUncorr*covUncorr(2,4)) +
412 pYUncorr * (pXUncorr*covUncorr(3,2) + pYUncorr*covUncorr(3,3) + pZUncorr*covUncorr(3,4)) +
413 pZUncorr * (pXUncorr*covUncorr(4,2) + pYUncorr*covUncorr(4,3) + pZUncorr*covUncorr(4,4))) / pUncorr;
414 ((TH2F*)fTrackRes->UncheckedAt(kUncorrPRes))->Fill(pUncorr,100.*sigmaPUncorr/pUncorr);
416 // fill transverse momentum error at first cluster
417 Double_t ptUncorr = TMath::Sqrt(pXUncorr*pXUncorr + pYUncorr*pYUncorr);
418 Double_t sigmaPtUncorr = TMath::Sqrt(pXUncorr * (pXUncorr*covUncorr(2,2) + pYUncorr*covUncorr(2,3)) + pYUncorr * (pXUncorr*covUncorr(3,2) + pYUncorr*covUncorr(3,3))) / ptUncorr;
419 ((TH2F*)fTrackRes->UncheckedAt(kUncorrPtRes))->Fill(ptUncorr,100.*sigmaPtUncorr/ptUncorr);
421 // fill slopeX-Y error at first cluster
422 ((TH2F*)fTrackRes->UncheckedAt(kUncorrSlopeRes))->Fill(pUncorr,TMath::Sqrt(trackParamAtFirstCluster.GetCovariances()(1,1)));
423 ((TH2F*)fTrackRes->UncheckedAt(kUncorrSlopeRes+1))->Fill(pUncorr,TMath::Sqrt(trackParamAtFirstCluster.GetCovariances()(3,3)));
425 // fill momentum error at vertex
426 AliMUONTrackParam trackParamAtVtx(trackParamAtFirstCluster);
427 AliMUONTrackExtrap::ExtrapToVertex(&trackParamAtVtx, esdTrack->GetNonBendingCoor(), esdTrack->GetBendingCoor(), esdTrack->GetZ(), 0., 0.);
428 Double_t pXVtx = trackParamAtVtx.Px();
429 Double_t pYVtx = trackParamAtVtx.Py();
430 Double_t pZVtx = trackParamAtVtx.Pz();
431 Double_t pVtx = trackParamAtVtx.P();
432 TMatrixD covVtx(5,5);
433 Cov2CovP(trackParamAtVtx,covVtx);
434 Double_t sigmaPVtx = TMath::Sqrt(pXVtx * (pXVtx*covVtx(2,2) + pYVtx*covVtx(2,3) + pZVtx*covVtx(2,4)) +
435 pYVtx * (pXVtx*covVtx(3,2) + pYVtx*covVtx(3,3) + pZVtx*covVtx(3,4)) +
436 pZVtx * (pXVtx*covVtx(4,2) + pYVtx*covVtx(4,3) + pZVtx*covVtx(4,4))) / pVtx;
437 ((TH2F*)fTrackRes->UncheckedAt(kPRes))->Fill(pVtx,100.*sigmaPVtx/pVtx);
439 // fill transverse momentum error at vertex
440 Double_t ptVtx = TMath::Sqrt(pXVtx*pXVtx + pYVtx*pYVtx);
441 Double_t sigmaPtVtx = TMath::Sqrt(pXVtx * (pXVtx*covVtx(2,2) + pYVtx*covVtx(2,3)) + pYVtx * (pXVtx*covVtx(3,2) + pYVtx*covVtx(3,3))) / ptVtx;
442 ((TH2F*)fTrackRes->UncheckedAt(kPtRes))->Fill(ptVtx,100.*sigmaPtVtx/ptVtx);
444 // fill slopeX-Y error at vertex
445 ((TH2F*)fTrackRes->UncheckedAt(kSlopeRes))->Fill(pVtx,TMath::Sqrt(trackParamAtVtx.GetCovariances()(1,1)));
446 ((TH2F*)fTrackRes->UncheckedAt(kSlopeRes+1))->Fill(pVtx,TMath::Sqrt(trackParamAtVtx.GetCovariances()(3,3)));
448 // loop over clusters
449 Int_t nClusters = track.GetNClusters();
450 for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
452 // Get current, previous and next trackParams
453 AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->UncheckedAt(iCluster));
454 AliMUONTrackParam* previousTrackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->Before(trackParam));
455 AliMUONTrackParam* nextTrackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->After(trackParam));
456 if (!previousTrackParam && !nextTrackParam) {
457 AliError(Form("unable to find a cluster neither before nor after the one at position %d !?!", iCluster));
458 track.RecursiveDump();
462 // remove mono-cathod clusters if required
463 if (fRemoveMonoCathCl) {
464 Bool_t hasBending, hasNonBending;
465 if (fCheckAllPads) CheckPads(trackParam->GetClusterPtr(), hasBending, hasNonBending);
466 else CheckPadsBelow(trackParam->GetClusterPtr(), hasBending, hasNonBending);
467 if (!hasBending || !hasNonBending) continue;
470 // save current trackParam and remove it from the track
471 AliMUONTrackParam currentTrackParam(*trackParam);
472 track.RemoveTrackParamAtCluster(trackParam);
475 AliMUONVCluster* cluster = currentTrackParam.GetClusterPtr();
476 Int_t chId = cluster->GetChamberId();
477 Int_t halfChId = (cluster->GetX() > 0) ? 2*chId : 2*chId+1;
478 Int_t deId = cluster->GetDetElemId();
480 // compute residuals with cluster still attached to the track
481 AliMUONTrackParam* referenceTrackParam = static_cast<AliMUONTrackParam*>(referenceTrack.GetTrackParamAtCluster()->UncheckedAt(iCluster));
482 Double_t deltaX = cluster->GetX() - referenceTrackParam->GetNonBendingCoor();
483 Double_t deltaY = cluster->GetY() - referenceTrackParam->GetBendingCoor();
485 // compute local chi2
486 Double_t sigmaDeltaX2 = cluster->GetErrX2() - referenceTrackParam->GetCovariances()(0,0);
487 Double_t sigmaDeltaY2 = cluster->GetErrY2() - referenceTrackParam->GetCovariances()(2,2);
488 Double_t localChi2X = (sigmaDeltaX2 > 0.) ? deltaX*deltaX/sigmaDeltaX2 : 0.;
489 Double_t localChi2Y = (sigmaDeltaY2 > 0.) ? deltaY*deltaY/sigmaDeltaY2 : 0.;
490 Double_t localChi2 = 0.5 * referenceTrackParam->GetLocalChi2();
492 // fill local chi2 info at every clusters
493 ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerCh))->Fill(chId+1, localChi2X);
494 ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerCh+1))->Fill(chId+1, localChi2Y);
495 ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerCh+2))->Fill(chId+1, localChi2);
496 ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerDE))->Fill(fDEIndices[deId], localChi2X);
497 ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerDE+1))->Fill(fDEIndices[deId], localChi2Y);
498 ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerDE+2))->Fill(fDEIndices[deId], localChi2);
500 // make sure the track has another cluster in the same station and can still be refitted
501 Bool_t refit = track.IsValid( 1 << (chId/2) );
504 // refit the track and proceed if everything goes fine
505 if (tracker->RefitTrack(track, kFALSE)) {
507 // fill histograms of residuals with cluster still attached to the track
508 ((TH2F*)fResiduals->UncheckedAt(kResidualPerChClusterIn))->Fill(chId+1, deltaX);
509 ((TH2F*)fResiduals->UncheckedAt(kResidualPerChClusterIn+1))->Fill(chId+1, deltaY);
510 ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfChClusterIn))->Fill(halfChId+1, deltaX);
511 ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfChClusterIn+1))->Fill(halfChId+1, deltaY);
512 ((TH2F*)fResiduals->UncheckedAt(kResidualPerDEClusterIn))->Fill(fDEIndices[deId], deltaX);
513 ((TH2F*)fResiduals->UncheckedAt(kResidualPerDEClusterIn+1))->Fill(fDEIndices[deId], deltaY);
514 ((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsPClusterIn+chId))->Fill(pUncorr, deltaX);
515 ((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsPClusterIn+10+chId))->Fill(pUncorr, deltaY);
516 ((TH2F*)fResidualsVsP->UncheckedAt(kResidualVsPClusterIn))->Fill(pUncorr, deltaX);
517 ((TH2F*)fResidualsVsP->UncheckedAt(kResidualVsPClusterIn+1))->Fill(pUncorr, deltaY);
518 ((TH2F*)fResidualsVsCent->UncheckedAt(kResidualInChVsCentClusterIn+chId))->Fill(centrality, deltaX);
519 ((TH2F*)fResidualsVsCent->UncheckedAt(kResidualInChVsCentClusterIn+10+chId))->Fill(centrality, deltaY);
520 ((TH2F*)fResidualsVsCent->UncheckedAt(kResidualVsCentClusterIn))->Fill(centrality, deltaX);
521 ((TH2F*)fResidualsVsCent->UncheckedAt(kResidualVsCentClusterIn+1))->Fill(centrality, deltaY);
523 // find the track parameters closest to the current cluster position
524 Double_t dZWithPrevious = (previousTrackParam) ? TMath::Abs(previousTrackParam->GetClusterPtr()->GetZ() - cluster->GetZ()) : FLT_MAX;
525 Int_t previousChId = (previousTrackParam) ? previousTrackParam->GetClusterPtr()->GetChamberId() : -1;
526 Double_t dZWithNext = (nextTrackParam) ? TMath::Abs(nextTrackParam->GetClusterPtr()->GetZ() - cluster->GetZ()) : FLT_MAX;
527 AliMUONTrackParam* startingTrackParam = (nextTrackParam) ? nextTrackParam : previousTrackParam;
528 if ((fExtrapMode == 0 && previousTrackParam && dZWithPrevious < dZWithNext) ||
529 (fExtrapMode == 1 && previousTrackParam && !(chId/2 == 2 && previousChId/2 == 1) &&
530 !(chId/2 == 3 && previousChId/2 == 2))) startingTrackParam = previousTrackParam;
532 // reset current parameters
533 currentTrackParam.SetParameters(startingTrackParam->GetParameters());
534 currentTrackParam.SetZ(startingTrackParam->GetZ());
535 currentTrackParam.SetCovariances(startingTrackParam->GetCovariances());
536 currentTrackParam.ResetPropagator();
538 // extrapolate to the current cluster position and fill histograms of residuals if everything goes fine
539 if (AliMUONTrackExtrap::ExtrapToZCov(¤tTrackParam, currentTrackParam.GetClusterPtr()->GetZ(), kTRUE)) {
541 // compute MCS dispersion on the first chamber
542 TMatrixD mcsCov(5,5);
543 if (startingTrackParam == nextTrackParam && chId == 0) {
544 AliMUONTrackParam trackParamForMCS;
545 trackParamForMCS.SetParameters(nextTrackParam->GetParameters());
546 AliMUONTrackExtrap::AddMCSEffect(&trackParamForMCS,AliMUONConstants::ChamberThicknessInX0(nextTrackParam->GetClusterPtr()->GetChamberId()),-1.);
547 const TMatrixD &propagator = currentTrackParam.GetPropagator();
548 TMatrixD tmp(trackParamForMCS.GetCovariances(),TMatrixD::kMultTranspose,propagator);
549 mcsCov.Mult(propagator,tmp);
550 } else mcsCov.Zero();
553 Double_t trackResX2 = currentTrackParam.GetCovariances()(0,0) + mcsCov(0,0);
554 Double_t trackResY2 = currentTrackParam.GetCovariances()(2,2) + mcsCov(2,2);
555 deltaX = cluster->GetX() - currentTrackParam.GetNonBendingCoor();
556 deltaY = cluster->GetY() - currentTrackParam.GetBendingCoor();
558 // fill histograms with cluster not attached to the track
559 ((TH2F*)fResiduals->UncheckedAt(kResidualPerChClusterOut))->Fill(chId+1, deltaX);
560 ((TH2F*)fResiduals->UncheckedAt(kResidualPerChClusterOut+1))->Fill(chId+1, deltaY);
561 ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfChClusterOut))->Fill(halfChId+1, deltaX);
562 ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfChClusterOut+1))->Fill(halfChId+1, deltaY);
563 ((TH2F*)fResiduals->UncheckedAt(kResidualPerDEClusterOut))->Fill(fDEIndices[deId], deltaX);
564 ((TH2F*)fResiduals->UncheckedAt(kResidualPerDEClusterOut+1))->Fill(fDEIndices[deId], deltaY);
565 ((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsPClusterOut+chId))->Fill(pUncorr, deltaX);
566 ((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsPClusterOut+10+chId))->Fill(pUncorr, deltaY);
567 ((TH2F*)fResidualsVsP->UncheckedAt(kResidualVsPClusterOut))->Fill(pUncorr, deltaX);
568 ((TH2F*)fResidualsVsP->UncheckedAt(kResidualVsPClusterOut+1))->Fill(pUncorr, deltaY);
569 ((TH2F*)fResidualsVsCent->UncheckedAt(kResidualInChVsCentClusterOut+chId))->Fill(centrality, deltaX);
570 ((TH2F*)fResidualsVsCent->UncheckedAt(kResidualInChVsCentClusterOut+10+chId))->Fill(centrality, deltaY);
571 ((TH2F*)fResidualsVsCent->UncheckedAt(kResidualVsCentClusterOut))->Fill(centrality, deltaX);
572 ((TH2F*)fResidualsVsCent->UncheckedAt(kResidualVsCentClusterOut+1))->Fill(centrality, deltaY);
573 ((TH2F*)fResiduals->UncheckedAt(kTrackResPerCh))->Fill(chId+1, TMath::Sqrt(trackResX2));
574 ((TH2F*)fResiduals->UncheckedAt(kTrackResPerCh+1))->Fill(chId+1, TMath::Sqrt(trackResY2));
575 ((TH2F*)fResiduals->UncheckedAt(kTrackResPerHalfCh))->Fill(halfChId+1, TMath::Sqrt(trackResX2));
576 ((TH2F*)fResiduals->UncheckedAt(kTrackResPerHalfCh+1))->Fill(halfChId+1, TMath::Sqrt(trackResY2));
577 ((TH2F*)fResiduals->UncheckedAt(kTrackResPerDE))->Fill(fDEIndices[deId], TMath::Sqrt(trackResX2));
578 ((TH2F*)fResiduals->UncheckedAt(kTrackResPerDE+1))->Fill(fDEIndices[deId], TMath::Sqrt(trackResY2));
579 ((TH2F*)fResiduals->UncheckedAt(kMCSPerCh))->Fill(chId+1, TMath::Sqrt(mcsCov(0,0)));
580 ((TH2F*)fResiduals->UncheckedAt(kMCSPerCh+1))->Fill(chId+1, TMath::Sqrt(mcsCov(2,2)));
581 ((TH2F*)fResiduals->UncheckedAt(kMCSPerHalfCh))->Fill(halfChId+1, TMath::Sqrt(mcsCov(0,0)));
582 ((TH2F*)fResiduals->UncheckedAt(kMCSPerHalfCh+1))->Fill(halfChId+1, TMath::Sqrt(mcsCov(2,2)));
583 ((TH2F*)fResiduals->UncheckedAt(kMCSPerDE))->Fill(fDEIndices[deId], TMath::Sqrt(mcsCov(0,0)));
584 ((TH2F*)fResiduals->UncheckedAt(kMCSPerDE+1))->Fill(fDEIndices[deId], TMath::Sqrt(mcsCov(2,2)));
585 ((TH2F*)fResiduals->UncheckedAt(kClusterRes2PerCh))->Fill(chId+1, deltaX*deltaX - trackResX2);
586 ((TH2F*)fResiduals->UncheckedAt(kClusterRes2PerCh+1))->Fill(chId+1, deltaY*deltaY - trackResY2);
594 track.AddTrackParamAtCluster(currentTrackParam, *(currentTrackParam.GetClusterPtr()), kTRUE);
600 // Post final data. It will be written to a file with option "RECREATE"
601 PostData(1, fResiduals);
602 PostData(2, fResidualsVsP);
603 PostData(3, fResidualsVsCent);
604 PostData(6, fTrackRes);
607 //________________________________________________________________________
608 void AliAnalysisTaskMuonResolution::NotifyRun()
610 /// load necessary data from OCDB corresponding to the first run number and initialize analysis
612 if (fOCDBLoaded) return;
614 AliCDBManager* cdbm = AliCDBManager::Instance();
615 cdbm->SetDefaultStorage(fDefaultStorage.Data());
616 cdbm->SetRun(fCurrentRunNumber);
618 if (!AliMUONCDB::LoadField()) return;
620 if (!AliMUONCDB::LoadMapping()) return;
622 AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam();
623 if (!recoParam) return;
625 AliMUONESDInterface::ResetTracker(recoParam);
627 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
629 // set the cluster resolution to default if not already set and create output objects
630 if (fClusterResNB[i] < 0.) fClusterResNB[i] = recoParam->GetDefaultNonBendingReso(i);
631 if (fClusterResB[i] < 0.) fClusterResB[i] = recoParam->GetDefaultBendingReso(i);
633 // fill correspondence between DEId and indices for histo (starting from 1)
636 while (!it.IsDone()) {
638 fDEIndices[it.CurrentDEId()] = fNDE;
639 fDEIds[fNDE] = it.CurrentDEId();
647 // recover default storage full name (raw:// cannot be used to set specific storage)
648 TString defaultStorage(cdbm->GetDefaultStorage()->GetType());
649 if (defaultStorage == "alien") defaultStorage += Form("://folder=%s", cdbm->GetDefaultStorage()->GetBaseFolder().Data());
650 else defaultStorage += Form("://%s", cdbm->GetDefaultStorage()->GetBaseFolder().Data());
652 // reset existing geometry/alignment if any
653 if (cdbm->GetEntryCache()->Contains("GRP/Geometry/Data")) cdbm->UnloadFromCache("GRP/Geometry/Data");
654 if (cdbm->GetEntryCache()->Contains("MUON/Align/Data")) cdbm->UnloadFromCache("MUON/Align/Data");
655 if (AliGeomManager::GetGeometry()) AliGeomManager::GetGeometry()->UnlockGeometry();
657 // get original geometry transformer
658 AliGeomManager::LoadGeometry();
659 if (!AliGeomManager::GetGeometry()) return;
660 if (fOldAlignStorage != "none") {
661 if (!fOldAlignStorage.IsNull()) cdbm->SetSpecificStorage("MUON/Align/Data",fOldAlignStorage.Data());
662 else cdbm->SetSpecificStorage("MUON/Align/Data",defaultStorage.Data());
663 AliGeomManager::ApplyAlignObjsFromCDB("MUON");
665 fOldGeoTransformer = new AliMUONGeometryTransformer();
666 fOldGeoTransformer->LoadGeometryData();
668 // get new geometry transformer
669 cdbm->UnloadFromCache("GRP/Geometry/Data");
670 if (fOldAlignStorage != "none") cdbm->UnloadFromCache("MUON/Align/Data");
671 AliGeomManager::GetGeometry()->UnlockGeometry();
672 AliGeomManager::LoadGeometry();
673 if (!AliGeomManager::GetGeometry()) return;
674 if (!fNewAlignStorage.IsNull()) cdbm->SetSpecificStorage("MUON/Align/Data",fNewAlignStorage.Data());
675 else cdbm->SetSpecificStorage("MUON/Align/Data",defaultStorage.Data());
676 AliGeomManager::ApplyAlignObjsFromCDB("MUON");
677 fNewGeoTransformer = new AliMUONGeometryTransformer();
678 fNewGeoTransformer->LoadGeometryData();
682 // load geometry for track extrapolation to vertex
683 if (cdbm->GetEntryCache()->Contains("GRP/Geometry/Data")) cdbm->UnloadFromCache("GRP/Geometry/Data");
684 if (cdbm->GetEntryCache()->Contains("MUON/Align/Data")) cdbm->UnloadFromCache("MUON/Align/Data");
685 if (AliGeomManager::GetGeometry()) AliGeomManager::GetGeometry()->UnlockGeometry();
686 AliGeomManager::LoadGeometry();
687 if (!AliGeomManager::GetGeometry()) return;
688 AliGeomManager::ApplyAlignObjsFromCDB("MUON");
689 fNewGeoTransformer = new AliMUONGeometryTransformer();
690 fNewGeoTransformer->LoadGeometryData();
694 // print starting chamber resolution if required
695 if (fPrintClResPerCh) {
696 printf("\nstarting chamber resolution:\n");
697 printf(" - non-bending:");
698 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %5.3f":", %5.3f",fClusterResNB[i]);
699 printf("\n - bending:");
700 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %6.4f":", %6.4f",fClusterResB[i]);
706 UserCreateOutputObjects();
710 //________________________________________________________________________
711 void AliAnalysisTaskMuonResolution::Terminate(Option_t *)
713 /// compute final results
715 // recover output objects
716 fResiduals = static_cast<TObjArray*> (GetOutputData(1));
717 fResidualsVsP = static_cast<TObjArray*> (GetOutputData(2));
718 fResidualsVsCent = static_cast<TObjArray*> (GetOutputData(3));
719 fTrackRes = static_cast<TObjArray*> (GetOutputData(6));
720 if (!fResiduals || !fResidualsVsP || !fTrackRes) return;
723 fLocalChi2 = new TObjArray(1000);
724 fLocalChi2->SetOwner();
725 fChamberRes = new TObjArray(1000);
726 fChamberRes->SetOwner();
730 const char* axes[2] = {"X", "Y"};
731 Double_t newClusterRes[2][10], newClusterResErr[2][10];
732 fNDE = ((TH2F*)fResiduals->UncheckedAt(kResidualPerDEClusterIn))->GetXaxis()->GetNbins();
734 for (Int_t ia = 0; ia < 2; ia++) {
736 g = new TGraphErrors(AliMUONConstants::NTrackingCh());
737 g->SetName(Form("gResidual%sPerChMean_ClusterIn",axes[ia]));
738 g->SetTitle(Form("cluster-track residual-%s per Ch: mean (cluster in);chamber ID;<#Delta_{%s}> (cm)",axes[ia],axes[ia]));
739 g->SetMarkerStyle(kFullDotLarge);
740 fChamberRes->AddAtAndExpand(g, kResidualPerChMeanClusterIn+ia);
741 g = new TGraphErrors(AliMUONConstants::NTrackingCh());
742 g->SetName(Form("gResidual%sPerChMean_ClusterOut",axes[ia]));
743 g->SetTitle(Form("cluster-track residual-%s per Ch: mean (cluster out);chamber ID;<#Delta_{%s}> (cm)",axes[ia],axes[ia]));
744 g->SetMarkerStyle(kFullDotLarge);
745 fChamberRes->AddAtAndExpand(g, kResidualPerChMeanClusterOut+ia);
747 g = new TGraphErrors(2*AliMUONConstants::NTrackingCh());
748 g->SetName(Form("gResidual%sPerHalfChMean_ClusterIn",axes[ia]));
749 g->SetTitle(Form("cluster-track residual-%s per half Ch: mean (cluster in);half chamber ID;<#Delta_{%s}> (cm)",axes[ia],axes[ia]));
750 g->SetMarkerStyle(kFullDotLarge);
751 fChamberRes->AddAtAndExpand(g, kResidualPerHalfChMeanClusterIn+ia);
752 g = new TGraphErrors(2*AliMUONConstants::NTrackingCh());
753 g->SetName(Form("gResidual%sPerHalfChMean_ClusterOut",axes[ia]));
754 g->SetTitle(Form("cluster-track residual-%s per half Ch: mean (cluster out);half chamber ID;<#Delta_{%s}> (cm)",axes[ia],axes[ia]));
755 g->SetMarkerStyle(kFullDotLarge);
756 fChamberRes->AddAtAndExpand(g, kResidualPerHalfChMeanClusterOut+ia);
758 g = new TGraphErrors(fNDE);
759 g->SetName(Form("gResidual%sPerDEMean_ClusterIn",axes[ia]));
760 g->SetTitle(Form("cluster-track residual-%s per DE: mean (cluster in);DE ID;<#Delta_{%s}> (cm)",axes[ia],axes[ia]));
761 g->SetMarkerStyle(kFullDotLarge);
762 fChamberRes->AddAtAndExpand(g, kResidualPerDEMeanClusterIn+ia);
763 g = new TGraphErrors(fNDE);
764 g->SetName(Form("gResidual%sPerDEMean_ClusterOut",axes[ia]));
765 g->SetTitle(Form("cluster-track residual-%s per DE: mean (cluster out);DE ID;<#Delta_{%s}> (cm)",axes[ia],axes[ia]));
766 g->SetMarkerStyle(kFullDotLarge);
767 fChamberRes->AddAtAndExpand(g, kResidualPerDEMeanClusterOut+ia);
769 g = new TGraphErrors(AliMUONConstants::NTrackingCh());
770 g->SetName(Form("gResidual%sPerChSigma_ClusterIn",axes[ia]));
771 g->SetTitle(Form("cluster-track residual-%s per Ch: sigma (cluster in);chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
772 g->SetMarkerStyle(kFullDotLarge);
773 fChamberRes->AddAtAndExpand(g, kResidualPerChSigmaClusterIn+ia);
774 g = new TGraphErrors(AliMUONConstants::NTrackingCh());
775 g->SetName(Form("gResidual%sPerChSigma_ClusterOut",axes[ia]));
776 g->SetTitle(Form("cluster-track residual-%s per Ch: sigma (cluster out);chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
777 g->SetMarkerStyle(kFullDotLarge);
778 fChamberRes->AddAtAndExpand(g, kResidualPerChSigmaClusterOut+ia);
780 g = new TGraphErrors(AliMUONConstants::NTrackingCh());
781 g->SetName(Form("gResidual%sPerChDispersion_ClusterOut",axes[ia]));
782 g->SetTitle(Form("cluster-track residual-%s per Ch: dispersion (cluster out);chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
783 g->SetMarkerStyle(kFullDotLarge);
784 fChamberRes->AddAtAndExpand(g, kResidualPerChDispersionClusterOut+ia);
786 g = new TGraphErrors(AliMUONConstants::NTrackingCh());
787 g->SetName(Form("gCombinedResidual%sPerChSigma",axes[ia]));
788 g->SetTitle(Form("combined cluster-track residual-%s per Ch: sigma;chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
789 g->SetMarkerStyle(kFullDotLarge);
790 fChamberRes->AddAtAndExpand(g, kCombinedResidualPerChSigma+ia);
792 g = new TGraphErrors(2*AliMUONConstants::NTrackingCh());
793 g->SetName(Form("gCombinedResidual%sPerHalfChSigma",axes[ia]));
794 g->SetTitle(Form("combined cluster-track residual-%s per half Ch: sigma;half chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
795 g->SetMarkerStyle(kFullDotLarge);
796 fChamberRes->AddAtAndExpand(g, kCombinedResidualPerHalfChSigma+ia);
798 g = new TGraphErrors(fNDE);
799 g->SetName(Form("gCombinedResidual%sPerDESigma",axes[ia]));
800 g->SetTitle(Form("combined cluster-track residual-%s per DE: sigma;DE ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
801 g->SetMarkerStyle(kFullDotLarge);
802 fChamberRes->AddAtAndExpand(g, kCombinedResidualPerDESigma+ia);
804 mg = new TMultiGraph(Form("mgCombinedResidual%sSigmaVsP",axes[ia]),Form("cluster %s-resolution per chamber versus momentum;p (GeV/c^{2});#sigma_{%s} (cm)",axes[ia],axes[ia]));
805 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
806 g = new TGraphErrors(((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsPClusterIn+10*ia+i))->GetNbinsX());
807 g->SetName(Form("gRes%sVsP_ch%d",axes[ia],i+1));
808 g->SetMarkerStyle(kFullDotMedium);
809 g->SetMarkerColor(i+1+i/9);
812 fChamberRes->AddAtAndExpand(mg, kCombinedResidualSigmaVsP+ia);
814 g = new TGraphErrors(((TH2F*)fResidualsVsP->UncheckedAt(kResidualVsPClusterIn+ia))->GetNbinsX());
815 g->SetName(Form("gCombinedResidual%sSigmaVsP",axes[ia]));
816 g->SetTitle(Form("cluster %s-resolution integrated over chambers versus momentum: sigma;p (GeV/c^{2});#sigma_{%s} (cm)",axes[ia],axes[ia]));
817 g->SetMarkerStyle(kFullDotLarge);
818 fChamberRes->AddAtAndExpand(g, kCombinedResidualAllChSigmaVsP+ia);
820 mg = new TMultiGraph(Form("mgCombinedResidual%sSigmaVsCent",axes[ia]),Form("cluster %s-resolution per chamber versus centrality;centrality (%%);#sigma_{%s} (cm)",axes[ia],axes[ia]));
821 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
822 g = new TGraphErrors(((TH2F*)fResidualsVsCent->UncheckedAt(kResidualInChVsCentClusterIn+10*ia+i))->GetNbinsX());
823 g->SetName(Form("gRes%sVsCent_ch%d",axes[ia],i+1));
824 g->SetMarkerStyle(kFullDotMedium);
825 g->SetMarkerColor(i+1+i/9);
828 fChamberRes->AddAtAndExpand(mg, kCombinedResidualSigmaVsCent+ia);
830 g = new TGraphErrors(((TH2F*)fResidualsVsCent->UncheckedAt(kResidualVsCentClusterIn+ia))->GetNbinsX());
831 g->SetName(Form("gCombinedResidual%sSigmaVsCent",axes[ia]));
832 g->SetTitle(Form("cluster %s-resolution integrated over chambers versus centrality: sigma;centrality (%%);#sigma_{%s} (cm)",axes[ia],axes[ia]));
833 g->SetMarkerStyle(kFullDotLarge);
834 fChamberRes->AddAtAndExpand(g, kCombinedResidualAllChSigmaVsCent+ia);
836 g = new TGraphErrors(AliMUONConstants::NTrackingCh());
837 g->SetName(Form("gTrackRes%sPerCh",axes[ia]));
838 g->SetTitle(Form("track <#sigma_{%s}> per Ch;chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
839 g->SetMarkerStyle(kFullDotLarge);
840 fChamberRes->AddAtAndExpand(g, kTrackResPerChMean+ia);
842 g = new TGraphErrors(AliMUONConstants::NTrackingCh());
843 g->SetName(Form("gMCS%sPerCh",axes[ia]));
844 g->SetTitle(Form("MCS %s-dispersion of extrapolated track per Ch;chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
845 g->SetMarkerStyle(kFullDotLarge);
846 fChamberRes->AddAtAndExpand(g, kMCSPerChMean+ia);
848 g = new TGraphErrors(AliMUONConstants::NTrackingCh());
849 g->SetName(Form("gClusterRes%sPerCh",axes[ia]));
850 g->SetTitle(Form("cluster <#sigma_{%s}> per Ch;chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
851 g->SetMarkerStyle(kFullDotLarge);
852 fChamberRes->AddAtAndExpand(g, kClusterResPerCh+ia);
854 g = new TGraphErrors(2*AliMUONConstants::NTrackingCh());
855 g->SetName(Form("gClusterRes%sPerHalfCh",axes[ia]));
856 g->SetTitle(Form("cluster <#sigma_{%s}> per half Ch;half chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
857 g->SetMarkerStyle(kFullDotLarge);
858 fChamberRes->AddAtAndExpand(g, kClusterResPerHalfCh+ia);
860 g = new TGraphErrors(fNDE);
861 g->SetName(Form("gClusterRes%sPerDE",axes[ia]));
862 g->SetTitle(Form("cluster <#sigma_{%s}> per DE;DE ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
863 g->SetMarkerStyle(kFullDotLarge);
864 fChamberRes->AddAtAndExpand(g, kClusterResPerDE+ia);
866 g = new TGraphErrors(AliMUONConstants::NTrackingCh());
867 g->SetName(Form("gCalcClusterRes%sPerCh",axes[ia]));
868 g->SetTitle(Form("calculated cluster <#sigma_{%s}> per Ch;chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
869 g->SetMarkerStyle(kFullDotLarge);
870 fChamberRes->AddAtAndExpand(g, kCalcClusterResPerCh+ia);
872 g = new TGraphErrors(AliMUONConstants::NTrackingCh());
873 g->SetName(Form("gLocalChi2%sPerChMean",axes[ia]));
874 g->SetTitle(Form("local chi2-%s per Ch: mean;chamber ID;<local #chi^{2}_{%s}>",axes[ia],axes[ia]));
875 g->SetMarkerStyle(kFullDotLarge);
876 fLocalChi2->AddAtAndExpand(g, kLocalChi2PerChMean+ia);
878 g = new TGraphErrors(fNDE);
879 g->SetName(Form("gLocalChi2%sPerDEMean",axes[ia]));
880 g->SetTitle(Form("local chi2-%s per DE: mean;DE ID;<local #chi^{2}_{%s}>",axes[ia],axes[ia]));
881 g->SetMarkerStyle(kFullDotLarge);
882 fLocalChi2->AddAtAndExpand(g, kLocalChi2PerDEMean+ia);
884 // compute residual mean and dispersion integrated over chambers versus p
885 FillSigmaClusterVsP((TH2F*)fResidualsVsP->UncheckedAt(kResidualVsPClusterIn+ia),
886 (TH2F*)fResidualsVsP->UncheckedAt(kResidualVsPClusterOut+ia),
887 (TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualAllChSigmaVsP+ia));
889 // compute residual mean and dispersion integrated over chambers versus centrality
890 FillSigmaClusterVsCent((TH2F*)fResidualsVsCent->UncheckedAt(kResidualVsCentClusterIn+ia),
891 (TH2F*)fResidualsVsCent->UncheckedAt(kResidualVsCentClusterOut+ia),
892 (TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualAllChSigmaVsCent+ia));
894 // compute residual mean and dispersion and averaged local chi2 per chamber and half chamber
895 Double_t meanIn, meanInErr, meanOut, meanOutErr, sigma, sigmaIn, sigmaInErr, sigmaOut, sigmaOutErr;
896 Double_t sigmaTrack, sigmaTrackErr, sigmaMCS, sigmaMCSErr, clusterRes, clusterResErr, sigmaCluster, sigmaClusterErr;
897 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
900 TH1D *tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerChClusterIn+ia))->ProjectionY("tmp",i+1,i+1,"e");
901 GetMean(tmp, meanIn, meanInErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChMeanClusterIn+ia), i, i+1);
902 GetRMS(tmp, sigmaIn, sigmaInErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChSigmaClusterIn+ia), i, i+1);
905 tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerChClusterOut+ia))->ProjectionY("tmp",i+1,i+1,"e");
906 GetMean(tmp, meanOut, meanOutErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChMeanClusterOut+ia), i, i+1);
907 GetRMS(tmp, sigmaOut, sigmaOutErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChSigmaClusterOut+ia), i, i+1);
910 if (fCorrectForSystematics) {
911 sigma = TMath::Sqrt(sigmaIn*sigmaIn + meanIn*meanIn);
912 sigmaInErr = (sigma>0) ? TMath::Sqrt(sigmaIn*sigmaIn*sigmaInErr*sigmaInErr + meanIn*meanIn*meanInErr*meanInErr) / sigma : 0.;
914 sigma = TMath::Sqrt(sigmaOut*sigmaOut + meanOut*meanOut);
915 sigmaOutErr = (sigma>0) ? TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + meanOut*meanOut*meanOutErr*meanOutErr) / sigma : 0.;
918 ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChDispersionClusterOut+ia))->SetPoint(i, i+1, sigmaOut);
919 ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChDispersionClusterOut+ia))->SetPointError(i, 0., sigmaOutErr);
921 clusterRes = TMath::Sqrt(sigmaIn*sigmaOut);
922 // clusterResErr = (clusterRes > 0.) ? 0.5 * TMath::Sqrt(sigmaInErr*sigmaInErr*sigmaOut*sigmaOut + sigmaIn*sigmaIn*sigmaOutErr*sigmaOutErr) / clusterRes : 0.;
923 clusterResErr = TMath::Sqrt(sigmaInErr*sigmaOutErr);
924 ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerChSigma+ia))->SetPoint(i, i+1, clusterRes);
925 ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerChSigma+ia))->SetPointError(i, 0., clusterResErr);
926 newClusterRes[ia][i] = clusterRes;
927 newClusterResErr[ia][i] = clusterResErr;
930 tmp = ((TH2F*)fResiduals->UncheckedAt(kTrackResPerCh+ia))->ProjectionY("tmp",i+1,i+1,"e");
931 GetMean(tmp, sigmaTrack, sigmaTrackErr, (TGraphErrors*)fChamberRes->UncheckedAt(kTrackResPerChMean+ia), i, i+1, kFALSE, kFALSE);
934 tmp = ((TH2F*)fResiduals->UncheckedAt(kMCSPerCh+ia))->ProjectionY("tmp",i+1,i+1,"e");
935 GetMean(tmp, sigmaMCS, sigmaMCSErr, (TGraphErrors*)fChamberRes->UncheckedAt(kMCSPerChMean+ia), i, i+1, kFALSE, kFALSE);
938 sigmaCluster = sigmaOut*sigmaOut - sigmaTrack*sigmaTrack;
939 if (sigmaCluster > 0.) {
940 sigmaCluster = TMath::Sqrt(sigmaCluster);
941 sigmaClusterErr = TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + sigmaTrack*sigmaTrack*sigmaTrackErr*sigmaTrackErr) / sigmaCluster;
944 sigmaClusterErr = 0.;
946 ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerCh+ia))->SetPoint(i, i+1, sigmaCluster);
947 ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerCh+ia))->SetPointError(i, 0., sigmaClusterErr);
950 tmp = ((TH2F*)fResiduals->UncheckedAt(kClusterRes2PerCh+ia))->ProjectionY("tmp",i+1,i+1,"e");
952 clusterRes = tmp->GetMean();
953 if (clusterRes > 0.) {
954 ((TGraphErrors*)fChamberRes->UncheckedAt(kCalcClusterResPerCh+ia))->SetPoint(i, i+1, TMath::Sqrt(clusterRes));
955 ((TGraphErrors*)fChamberRes->UncheckedAt(kCalcClusterResPerCh+ia))->SetPointError(i, 0., 0.5 * tmp->GetMeanError() / TMath::Sqrt(clusterRes));
957 ((TGraphErrors*)fChamberRes->UncheckedAt(kCalcClusterResPerCh+ia))->SetPoint(i, i+1, 0.);
958 ((TGraphErrors*)fChamberRes->UncheckedAt(kCalcClusterResPerCh+ia))->SetPointError(i, 0., 0.);
963 FillSigmaClusterVsP((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsPClusterIn+10*ia+i),
964 (TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsPClusterOut+10*ia+i),
965 (TGraphErrors*)((TMultiGraph*)fChamberRes->UncheckedAt(kCombinedResidualSigmaVsP+ia))->GetListOfGraphs()->FindObject(Form("gRes%sVsP_ch%d",axes[ia],i+1)));
967 // method 1 versus centrality
968 FillSigmaClusterVsCent((TH2F*)fResidualsVsCent->UncheckedAt(kResidualInChVsCentClusterIn+10*ia+i),
969 (TH2F*)fResidualsVsCent->UncheckedAt(kResidualInChVsCentClusterOut+10*ia+i),
970 (TGraphErrors*)((TMultiGraph*)fChamberRes->UncheckedAt(kCombinedResidualSigmaVsCent+ia))->GetListOfGraphs()->FindObject(Form("gRes%sVsCent_ch%d",axes[ia],i+1)));
972 // compute residual mean and dispersion per half chamber
973 for (Int_t j = 0; j < 2; j++) {
977 tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfChClusterIn+ia))->ProjectionY("tmp",k+1,k+1,"e");
978 GetMean(tmp, meanIn, meanInErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMeanClusterIn+ia), k, k+1);
979 GetRMS(tmp, sigmaIn, sigmaInErr);
982 tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfChClusterOut+ia))->ProjectionY("tmp",k+1,k+1,"e");
983 GetMean(tmp, meanOut, meanOutErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMeanClusterOut+ia), k, k+1);
984 GetRMS(tmp, sigmaOut, sigmaOutErr);
987 if (fCorrectForSystematics) {
988 sigma = TMath::Sqrt(sigmaIn*sigmaIn + meanIn*meanIn);
989 sigmaInErr = (sigma>0) ? TMath::Sqrt(sigmaIn*sigmaIn*sigmaInErr*sigmaInErr + meanIn*meanIn*meanInErr*meanInErr) / sigma : 0.;
991 sigma = TMath::Sqrt(sigmaOut*sigmaOut + meanOut*meanOut);
992 sigmaOutErr = (sigma>0) ? TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + meanOut*meanOut*meanOutErr*meanOutErr) / sigma : 0.;
996 clusterRes = TMath::Sqrt(sigmaIn*sigmaOut);
997 // clusterResErr = (clusterRes > 0.) ? 0.5 * TMath::Sqrt(sigmaInErr*sigmaInErr*sigmaOut*sigmaOut + sigmaIn*sigmaIn*sigmaOutErr*sigmaOutErr) / clusterRes : 0.;
998 clusterResErr = TMath::Sqrt(sigmaInErr*sigmaOutErr);
999 ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerHalfChSigma+ia))->SetPoint(k, k+1, clusterRes);
1000 ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerHalfChSigma+ia))->SetPointError(k, 0., clusterResErr);
1003 tmp = ((TH2F*)fResiduals->UncheckedAt(kTrackResPerHalfCh+ia))->ProjectionY("tmp",k+1,k+1,"e");
1004 GetMean(tmp, sigmaTrack, sigmaTrackErr, 0x0, 0, 0, kFALSE, kFALSE);
1007 tmp = ((TH2F*)fResiduals->UncheckedAt(kMCSPerHalfCh+ia))->ProjectionY("tmp",k+1,k+1,"e");
1008 GetMean(tmp, sigmaMCS, sigmaMCSErr, 0x0, 0, 0, kFALSE, kFALSE);
1011 sigmaCluster = sigmaOut*sigmaOut - sigmaTrack*sigmaTrack;
1012 if (sigmaCluster > 0.) {
1013 sigmaCluster = TMath::Sqrt(sigmaCluster);
1014 sigmaClusterErr = TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + sigmaTrack*sigmaTrack*sigmaTrackErr*sigmaTrackErr) / sigmaCluster;
1017 sigmaClusterErr = 0.;
1019 ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerHalfCh+ia))->SetPoint(k, k+1, sigmaCluster);
1020 ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerHalfCh+ia))->SetPointError(k, 0., sigmaClusterErr);
1024 // compute averaged local chi2
1025 tmp = ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerCh+ia))->ProjectionY("tmp",i+1,i+1,"e");
1026 ((TGraphErrors*)fLocalChi2->UncheckedAt(kLocalChi2PerChMean+ia))->SetPoint(i, i+1, tmp->GetMean());
1027 ((TGraphErrors*)fLocalChi2->UncheckedAt(kLocalChi2PerChMean+ia))->SetPointError(i, 0., tmp->GetMeanError());
1032 // compute residual mean and dispersion per DE
1033 for (Int_t i = 0; i < fNDE; i++) {
1036 TH1D *tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerDEClusterIn+ia))->ProjectionY("tmp",i+1,i+1,"e");
1037 GetMean(tmp, meanIn, meanInErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMeanClusterIn+ia), i, i+1);
1038 GetRMS(tmp, sigmaIn, sigmaInErr);
1041 tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerDEClusterOut+ia))->ProjectionY("tmp",i+1,i+1,"e");
1042 GetMean(tmp, meanOut, meanOutErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMeanClusterOut+ia), i, i+1);
1043 GetRMS(tmp, sigmaOut, sigmaOutErr);
1046 if (fCorrectForSystematics) {
1047 sigma = TMath::Sqrt(sigmaIn*sigmaIn + meanIn*meanIn);
1048 sigmaInErr = (sigma>0) ? TMath::Sqrt(sigmaIn*sigmaIn*sigmaInErr*sigmaInErr + meanIn*meanIn*meanInErr*meanInErr) / sigma : 0.;
1050 sigma = TMath::Sqrt(sigmaOut*sigmaOut + meanOut*meanOut);
1051 sigmaOutErr = (sigma>0) ? TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + meanOut*meanOut*meanOutErr*meanOutErr) / sigma : 0.;
1055 clusterRes = TMath::Sqrt(sigmaIn*sigmaOut);
1056 // clusterResErr = (clusterRes > 0.) ? 0.5 * TMath::Sqrt(sigmaInErr*sigmaInErr*sigmaOut*sigmaOut + sigmaIn*sigmaIn*sigmaOutErr*sigmaOutErr) / clusterRes : 0.;
1057 clusterResErr = TMath::Sqrt(sigmaInErr*sigmaOutErr);
1058 ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerDESigma+ia))->SetPoint(i, i+1, clusterRes);
1059 ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerDESigma+ia))->SetPointError(i, 0., clusterResErr);
1062 tmp = ((TH2F*)fResiduals->UncheckedAt(kTrackResPerDE+ia))->ProjectionY("tmp",i+1,i+1,"e");
1063 GetMean(tmp, sigmaTrack, sigmaTrackErr, 0x0, 0, 0, kFALSE, kFALSE);
1066 tmp = ((TH2F*)fResiduals->UncheckedAt(kMCSPerDE+ia))->ProjectionY("tmp",i+1,i+1,"e");
1067 GetMean(tmp, sigmaMCS, sigmaMCSErr, 0x0, 0, 0, kFALSE, kFALSE);
1070 sigmaCluster = sigmaOut*sigmaOut - sigmaTrack*sigmaTrack;
1071 if (sigmaCluster > 0.) {
1072 sigmaCluster = TMath::Sqrt(sigmaCluster);
1073 sigmaClusterErr = TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + sigmaTrack*sigmaTrack*sigmaTrackErr*sigmaTrackErr) / sigmaCluster;
1076 sigmaClusterErr = 0.;
1078 ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerDE+ia))->SetPoint(i, i+1, sigmaCluster);
1079 ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerDE+ia))->SetPointError(i, 0., sigmaClusterErr);
1081 // compute averaged local chi2
1082 tmp = ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerDE+ia))->ProjectionY("tmp",i+1,i+1,"e");
1083 ((TGraphErrors*)fLocalChi2->UncheckedAt(kLocalChi2PerDEMean+ia))->SetPoint(i, i+1, tmp->GetMean());
1084 ((TGraphErrors*)fLocalChi2->UncheckedAt(kLocalChi2PerDEMean+ia))->SetPointError(i, 0., tmp->GetMeanError());
1089 // set half-chamber graph labels
1090 TAxis* xAxis = ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfChClusterIn+ia))->GetXaxis();
1091 ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMeanClusterIn+ia))->GetXaxis()->Set(20, 0.5, 20.5);
1092 ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMeanClusterOut+ia))->GetXaxis()->Set(20, 0.5, 20.5);
1093 ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerHalfChSigma+ia))->GetXaxis()->Set(20, 0.5, 20.5);
1094 ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerHalfCh+ia))->GetXaxis()->Set(20, 0.5, 20.5);
1095 for (Int_t i = 1; i <= 20; i++) {
1096 const char* label = xAxis->GetBinLabel(i);
1097 ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMeanClusterIn+ia))->GetXaxis()->SetBinLabel(i, label);
1098 ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMeanClusterOut+ia))->GetXaxis()->SetBinLabel(i, label);
1099 ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerHalfChSigma+ia))->GetXaxis()->SetBinLabel(i, label);
1100 ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerHalfCh+ia))->GetXaxis()->SetBinLabel(i, label);
1103 // set DE graph labels
1104 xAxis = ((TH2F*)fResiduals->UncheckedAt(kResidualPerDEClusterOut+ia))->GetXaxis();
1105 ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMeanClusterIn+ia))->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
1106 ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMeanClusterOut+ia))->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
1107 ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerDESigma+ia))->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
1108 ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerDE+ia))->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
1109 ((TGraphErrors*)fLocalChi2->UncheckedAt(kLocalChi2PerDEMean+ia))->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
1110 for (Int_t i = 1; i <= fNDE; i++) {
1111 const char* label = xAxis->GetBinLabel(i);
1112 ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMeanClusterIn+ia))->GetXaxis()->SetBinLabel(i, label);
1113 ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMeanClusterOut+ia))->GetXaxis()->SetBinLabel(i, label);
1114 ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerDESigma+ia))->GetXaxis()->SetBinLabel(i, label);
1115 ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerDE+ia))->GetXaxis()->SetBinLabel(i, label);
1116 ((TGraphErrors*)fLocalChi2->UncheckedAt(kLocalChi2PerDEMean+ia))->GetXaxis()->SetBinLabel(i, label);
1121 // compute averaged local chi2 per chamber (X+Y)
1122 TH2F* h2 = (TH2F*)fResiduals->UncheckedAt(kLocalChi2PerCh+2);
1123 g = new TGraphErrors(AliMUONConstants::NTrackingCh());
1124 g->SetName("gLocalChi2PerChMean");
1125 g->SetTitle("local chi2 per Ch: mean;chamber ID;<local #chi^{2}>");
1126 g->SetMarkerStyle(kFullDotLarge);
1127 fLocalChi2->AddAtAndExpand(g, kLocalChi2PerChMean+2);
1128 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
1129 TH1D* tmp = h2->ProjectionY("tmp",i+1,i+1,"e");
1130 g->SetPoint(i, i+1, tmp->GetMean());
1131 g->SetPointError(i, 0., tmp->GetMeanError());
1135 // compute averaged local chi2 per DE (X+Y)
1136 h2 = (TH2F*)fResiduals->UncheckedAt(kLocalChi2PerDE+2);
1137 g = new TGraphErrors(fNDE);
1138 g->SetName("gLocalChi2PerDEMean");
1139 g->SetTitle("local chi2 per DE: mean;DE ID;<local #chi^{2}>");
1140 g->SetMarkerStyle(kFullDotLarge);
1141 fLocalChi2->AddAtAndExpand(g, kLocalChi2PerDEMean+2);
1142 for (Int_t i = 0; i < fNDE; i++) {
1143 TH1D* tmp = h2->ProjectionY("tmp",i+1,i+1,"e");
1144 g->SetPoint(i, i+1, tmp->GetMean());
1145 g->SetPointError(i, 0., tmp->GetMeanError());
1150 g->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
1151 for (Int_t i = 1; i <= fNDE; i++) {
1152 const char* label = h2->GetXaxis()->GetBinLabel(i);
1153 g->GetXaxis()->SetBinLabel(i, label);
1157 fCanvases = new TObjArray(1000);
1158 fCanvases->SetOwner();
1160 TLegend *lResPerChMean = new TLegend(0.75,0.85,0.99,0.99);
1161 TLegend *lResPerChSigma1 = new TLegend(0.75,0.85,0.99,0.99);
1162 TLegend *lResPerChSigma2 = new TLegend(0.75,0.85,0.99,0.99);
1163 TLegend *lResPerChSigma3 = new TLegend(0.75,0.85,0.99,0.99);
1165 TCanvas* cResPerCh = new TCanvas("cResPerCh","cResPerCh",1200,500);
1166 cResPerCh->Divide(4,2);
1167 for (Int_t ia = 0; ia < 2; ia++) {
1168 cResPerCh->cd(1+4*ia);
1169 g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChMeanClusterOut+ia);
1171 g->SetMarkerColor(2);
1173 if (ia == 0) lResPerChMean->AddEntry(g,"cluster out","PL");
1174 g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChMeanClusterIn+ia);
1176 g->SetMarkerColor(4);
1178 if (ia == 0) lResPerChMean->AddEntry(g,"cluster in","PL");
1179 if (ia == 0) lResPerChMean->Draw();
1180 else lResPerChMean->DrawClone();
1181 cResPerCh->cd(2+4*ia);
1182 g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChSigmaClusterOut+ia);
1185 g->SetMarkerColor(2);
1187 if (ia == 0) lResPerChSigma1->AddEntry(g,"cluster out","PL");
1188 g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChSigmaClusterIn+ia);
1190 g->SetMarkerColor(4);
1192 if (ia == 0) lResPerChSigma1->AddEntry(g,"cluster in","PL");
1193 g = (TGraphErrors*)fChamberRes->UncheckedAt(kMCSPerChMean+ia);
1195 g->SetMarkerColor(5);
1197 if (ia == 0) lResPerChSigma1->AddEntry(g,"MCS","PL");
1198 g = (TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerChSigma+ia);
1200 g->SetMarkerColor(3);
1202 if (ia == 0) lResPerChSigma1->AddEntry(g,"combined 1","PL");
1203 if (ia == 0) lResPerChSigma1->Draw();
1204 else lResPerChSigma1->DrawClone();
1205 cResPerCh->cd(3+4*ia);
1206 g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChDispersionClusterOut+ia);
1209 g->SetMarkerColor(2);
1211 if (ia == 0) lResPerChSigma2->AddEntry(g,"cluster out","PL");
1212 g = (TGraphErrors*)fChamberRes->UncheckedAt(kMCSPerChMean+ia);
1214 if (ia == 0) lResPerChSigma2->AddEntry(g,"MCS","PL");
1215 g = (TGraphErrors*)fChamberRes->UncheckedAt(kTrackResPerChMean+ia);
1217 g->SetMarkerColor(4);
1219 if (ia == 0) lResPerChSigma2->AddEntry(g,"track res.","PL");
1220 g = (TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerCh+ia);
1222 if (ia == 0) lResPerChSigma2->AddEntry(g,"combined 2","PL");
1223 if (ia == 0) lResPerChSigma2->Draw();
1224 else lResPerChSigma2->DrawClone();
1225 cResPerCh->cd(4+4*ia);
1226 g = (TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerChSigma+ia);
1229 if (ia == 0) lResPerChSigma3->AddEntry(g,"combined 1","PL");
1230 g = (TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerCh+ia);
1232 if (ia == 0) lResPerChSigma3->AddEntry(g,"combined 2","PL");
1233 if (ia == 0) lResPerChSigma3->Draw();
1234 else lResPerChSigma3->DrawClone();
1236 fCanvases->AddAtAndExpand(cResPerCh, kResPerCh);
1238 TCanvas* cResPerHalfCh = new TCanvas("cResPerHalfCh","cResPerHalfCh",1200,500);
1239 cResPerHalfCh->Divide(2,2);
1240 for (Int_t ia = 0; ia < 2; ia++) {
1241 cResPerHalfCh->cd(1+2*ia);
1242 g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMeanClusterOut+ia);
1244 g->SetMarkerColor(2);
1246 g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMeanClusterIn+ia);
1248 g->SetMarkerColor(4);
1250 lResPerChMean->DrawClone();
1251 cResPerHalfCh->cd(2+2*ia);
1252 g = (TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerHalfChSigma+ia);
1255 g->SetMarkerColor(3);
1257 g = (TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerHalfCh+ia);
1259 lResPerChSigma3->DrawClone();
1261 fCanvases->AddAtAndExpand(cResPerHalfCh, kResPerHalfCh);
1263 TCanvas* cResPerDE = new TCanvas("cResPerDE","cResPerDE",1200,800);
1264 cResPerDE->Divide(1,4);
1265 for (Int_t ia = 0; ia < 2; ia++) {
1266 cResPerDE->cd(1+ia);
1267 g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMeanClusterOut+ia);
1269 g->SetMarkerColor(2);
1271 g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMeanClusterIn+ia);
1273 g->SetMarkerColor(4);
1275 lResPerChMean->DrawClone();
1276 cResPerDE->cd(3+ia);
1277 g = (TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerDESigma+ia);
1280 g->SetMarkerColor(3);
1282 g = (TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerDE+ia);
1284 lResPerChSigma3->DrawClone();
1286 fCanvases->AddAtAndExpand(cResPerDE, kResPerDE);
1288 TCanvas* cResPerChVsP = new TCanvas("cResPerChVsP","cResPerChVsP");
1289 cResPerChVsP->Divide(1,2);
1290 for (Int_t ia = 0; ia < 2; ia++) {
1291 cResPerChVsP->cd(1+ia);
1292 mg = (TMultiGraph*)fChamberRes->UncheckedAt(kCombinedResidualSigmaVsP+ia);
1295 fCanvases->AddAtAndExpand(cResPerChVsP, kResPerChVsP);
1297 TCanvas* cResPerChVsCent = new TCanvas("cResPerChVsCent","cResPerChVsCent");
1298 cResPerChVsCent->Divide(1,2);
1299 for (Int_t ia = 0; ia < 2; ia++) {
1300 cResPerChVsCent->cd(1+ia);
1301 mg = (TMultiGraph*)fChamberRes->UncheckedAt(kCombinedResidualSigmaVsCent+ia);
1304 fCanvases->AddAtAndExpand(cResPerChVsCent, kResPerChVsCent);
1307 if (fPrintClResPerCh) {
1308 printf("\nchamber resolution:\n");
1309 printf(" - non-bending:");
1310 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %5.3f":", %5.3f",newClusterRes[0][i]);
1311 printf("\n - bending:");
1312 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %6.4f":", %6.4f",newClusterRes[1][i]);
1316 if (fPrintClResPerDE) {
1317 Double_t iDE, clRes;
1318 printf("\nDE resolution:\n");
1319 printf(" - non-bending:");
1320 for (Int_t i = 0; i < fNDE; i++) {
1321 ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerDESigma))->GetPoint(i, iDE, clRes);
1322 printf((i==0)?" %5.3f":", %5.3f", clRes);
1324 printf("\n - bending:");
1325 for (Int_t i = 0; i < fNDE; i++) {
1326 ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerDESigma+1))->GetPoint(i, iDE, clRes);
1327 printf((i==0)?" %6.4f":", %6.4f", clRes);
1333 PostData(4, fLocalChi2);
1334 PostData(5, fChamberRes);
1337 //________________________________________________________________________
1338 void AliAnalysisTaskMuonResolution::ModifyClusters(AliMUONTrack& track)
1340 /// Reset the clusters resolution from the ones given to the task and change
1341 /// the cluster position according to the new alignment parameters if required
1343 Double_t gX,gY,gZ,lX,lY,lZ;
1345 // loop over clusters
1346 Int_t nClusters = track.GetNClusters();
1347 for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
1349 AliMUONVCluster* cl = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->UncheckedAt(iCluster))->GetClusterPtr();
1351 // change their resolution
1352 cl->SetErrXY(fClusterResNB[cl->GetChamberId()], fClusterResB[cl->GetChamberId()]);
1354 // change their position
1359 fOldGeoTransformer->Global2Local(cl->GetDetElemId(),gX,gY,gZ,lX,lY,lZ);
1360 fNewGeoTransformer->Local2Global(cl->GetDetElemId(),lX,lY,lZ,gX,gY,gZ);
1361 cl->SetXYZ(gX,gY,gZ);
1368 //________________________________________________________________________
1369 void AliAnalysisTaskMuonResolution::Zoom(TH1* h, Double_t fractionCut)
1371 /// Reduce the range of the histogram by removing a given fration of the statistic at each edge
1372 ZoomLeft(h, fractionCut);
1373 ZoomRight(h, fractionCut);
1376 //________________________________________________________________________
1377 void AliAnalysisTaskMuonResolution::ZoomLeft(TH1* h, Double_t fractionCut)
1379 /// Reduce the range of the histogram by removing a given fration of the statistic on the left side
1380 Int_t maxEventsCut = (Int_t) (fractionCut * h->GetEntries());
1381 Int_t nBins = h->GetNbinsX();
1385 Int_t eventsCut = 0;
1386 for (minBin = 1; minBin <= nBins; minBin++) {
1387 eventsCut += (Int_t) h->GetBinContent(minBin);
1388 if (eventsCut > maxEventsCut) break;
1391 // set new axis range
1392 h->GetXaxis()->SetRange(--minBin, h->GetXaxis()->GetLast());
1395 //________________________________________________________________________
1396 void AliAnalysisTaskMuonResolution::ZoomRight(TH1* h, Double_t fractionCut)
1398 /// Reduce the range of the histogram by removing a given fration of the statistic on the right side
1399 Int_t maxEventsCut = (Int_t) (fractionCut * h->GetEntries());
1400 Int_t nBins = h->GetNbinsX();
1404 Int_t eventsCut = 0;
1405 for (maxBin = nBins; maxBin >= 1; maxBin--) {
1406 eventsCut += (Int_t) h->GetBinContent(maxBin);
1407 if (eventsCut > maxEventsCut) break;
1410 // set new axis range
1411 h->GetXaxis()->SetRange(h->GetXaxis()->GetFirst(), ++maxBin);
1414 //________________________________________________________________________
1415 void AliAnalysisTaskMuonResolution::GetMean(TH1* h, Double_t& mean, Double_t& meanErr, TGraphErrors* g, Int_t i, Double_t x, Bool_t zoom, Bool_t enableFit)
1417 /// Fill graph with the mean value and the corresponding error (zooming if required)
1419 if (h->GetEntries() < fgkMinEntries) { // not enough entries
1424 } else if (enableFit && fGaus) { // take the mean of a gaussian fit
1426 fGaus->SetParameters(h->GetEntries(), 0., 0.1);
1428 if (h->GetUniqueID() != 10) {
1431 h->Fit("fGaus", "WWNQ");
1434 Int_t rebin = TMath::Max(static_cast<Int_t>(0.2*fGaus->GetParameter(2)/h->GetBinWidth(1)),1);
1435 while (h->GetNbinsX()%rebin!=0) rebin--;
1438 // use the unique ID to remember that this histogram has already been rebinned
1444 h->Fit("fGaus","NQ");
1446 mean = fGaus->GetParameter(1);
1447 meanErr = fGaus->GetParError(1);
1449 } else { // take the mean of the distribution
1453 mean = h->GetMean();
1454 meanErr = h->GetMeanError();
1456 if (zoom) h->GetXaxis()->SetRange(0,0);
1460 // fill graph if required
1462 g->SetPoint(i, x, mean);
1463 g->SetPointError(i, 0., meanErr);
1468 //________________________________________________________________________
1469 void AliAnalysisTaskMuonResolution::GetRMS(TH1* h, Double_t& rms, Double_t& rmsErr, TGraphErrors* g, Int_t i, Double_t x, Bool_t zoom)
1471 /// Return the dispersion value and the corresponding error (zooming if required) and fill graph if !=0x0
1473 if (h->GetEntries() < fgkMinEntries) { // not enough entries
1478 } else if (fGaus) { // take the sigma of a gaussian fit
1480 fGaus->SetParameters(h->GetEntries(), 0., 0.1);
1482 if (h->GetUniqueID() != 10) {
1485 h->Fit("fGaus", "WWNQ");
1488 Int_t rebin = TMath::Max(static_cast<Int_t>(0.2*fGaus->GetParameter(2)/h->GetBinWidth(1)),1);
1489 while (h->GetNbinsX()%rebin!=0) rebin--;
1492 // use the unique ID to remember that this histogram has already been rebinned
1498 h->Fit("fGaus","NQ");
1500 rms = fGaus->GetParameter(2);
1501 rmsErr = fGaus->GetParError(2);
1503 } else { // take the RMS of the distribution
1508 rmsErr = h->GetRMSError();
1510 if (zoom) h->GetXaxis()->SetRange(0,0);
1514 // fill graph if required
1516 g->SetPoint(i, x, rms);
1517 g->SetPointError(i, 0., rmsErr);
1522 //________________________________________________________________________
1523 void AliAnalysisTaskMuonResolution::FillSigmaClusterVsP(const TH2* hIn, const TH2* hOut, TGraphErrors* g, Bool_t zoom)
1525 /// Fill graph with cluster resolution from combined residuals with cluster in/out (zooming if required)
1526 Double_t sigmaIn, sigmaInErr, sigmaOut, sigmaOutErr, clusterRes, clusterResErr;
1527 for (Int_t j = 1; j <= hIn->GetNbinsX(); j++) {
1528 TH1D* tmp = hIn->ProjectionY("tmp",j,j,"e");
1529 GetRMS(tmp, sigmaIn, sigmaInErr, 0x0, 0, 0., zoom);
1531 tmp = hOut->ProjectionY("tmp",j,j,"e");
1532 GetRMS(tmp, sigmaOut, sigmaOutErr, 0x0, 0, 0., zoom);
1534 Double_t p = 0.5 * (hIn->GetBinLowEdge(j) + hIn->GetBinLowEdge(j+1));
1535 Double_t pErr = p - hIn->GetBinLowEdge(j);
1536 clusterRes = TMath::Sqrt(sigmaIn*sigmaOut);
1537 //clusterResErr = (clusterRes > 0.) ? 0.5 * TMath::Sqrt(sigmaInErr*sigmaInErr*sigmaOut*sigmaOut + sigmaIn*sigmaIn*sigmaOutErr*sigmaOutErr) / clusterRes : 0.;
1538 clusterResErr = TMath::Sqrt(sigmaInErr*sigmaOutErr);
1539 g->SetPoint(j, p, clusterRes);
1540 g->SetPointError(j, pErr, clusterResErr);
1544 //________________________________________________________________________
1545 void AliAnalysisTaskMuonResolution::FillSigmaClusterVsCent(const TH2* hIn, const TH2* hOut, TGraphErrors* g, Bool_t zoom)
1547 /// Fill graph with cluster resolution from combined residuals with cluster in/out (zooming if required)
1548 Double_t sigmaIn, sigmaInErr, sigmaOut, sigmaOutErr, clusterRes, clusterResErr;
1549 for (Int_t j = 1; j <= hIn->GetNbinsX(); j++) {
1550 TH1D* tmp = hIn->ProjectionY("tmp",j,j,"e");
1551 GetRMS(tmp, sigmaIn, sigmaInErr, 0x0, 0, 0., zoom);
1553 tmp = hOut->ProjectionY("tmp",j,j,"e");
1554 GetRMS(tmp, sigmaOut, sigmaOutErr, 0x0, 0, 0., zoom);
1556 Double_t cent = 0.5 * (hIn->GetBinLowEdge(j) + hIn->GetBinLowEdge(j+1));
1557 Double_t centErr = cent - hIn->GetBinLowEdge(j);
1558 clusterRes = TMath::Sqrt(sigmaIn*sigmaOut);
1559 //clusterResErr = (clusterRes > 0.) ? 0.5 * TMath::Sqrt(sigmaInErr*sigmaInErr*sigmaOut*sigmaOut + sigmaIn*sigmaIn*sigmaOutErr*sigmaOutErr) / clusterRes : 0.;
1560 clusterResErr = TMath::Sqrt(sigmaInErr*sigmaOutErr);
1561 g->SetPoint(j, cent, clusterRes);
1562 g->SetPointError(j, centErr, clusterResErr);
1566 //__________________________________________________________________________
1567 void AliAnalysisTaskMuonResolution::Cov2CovP(const AliMUONTrackParam ¶m, TMatrixD &covP)
1569 /// change coordinate system: (X, SlopeX, Y, SlopeY, q/Pyz) -> (X, Y, pX, pY, pZ)
1570 /// parameters (param) are given in the (X, SlopeX, Y, SlopeY, q/Pyz) coordinate system
1572 // Get useful parameters
1573 Double_t slopeX = param.GetNonBendingSlope();
1574 Double_t slopeY = param.GetBendingSlope();
1575 Double_t qOverPYZ = param.GetInverseBendingMomentum();
1576 Double_t pZ = param.Pz();
1578 // compute Jacobian to change the coordinate system from (X,SlopeX,Y,SlopeY,c/pYZ) to (X,Y,pX,pY,pZ)
1579 Double_t dpZdSlopeY = - qOverPYZ * qOverPYZ * pZ * pZ * pZ * slopeY;
1580 Double_t dpZdQOverPYZ = (qOverPYZ != 0.) ? - pZ / qOverPYZ : - FLT_MAX;
1581 TMatrixD jacob(5,5);
1586 jacob(2,3) = slopeX * dpZdSlopeY;
1587 jacob(2,4) = slopeX * dpZdQOverPYZ;
1588 jacob(3,3) = pZ + slopeY * dpZdSlopeY;
1589 jacob(3,4) = slopeY * dpZdQOverPYZ;
1590 jacob(4,3) = dpZdSlopeY;
1591 jacob(4,4) = dpZdQOverPYZ;
1593 // compute covariances in new coordinate system
1594 TMatrixD tmp(param.GetCovariances(),TMatrixD::kMultTranspose,jacob);
1595 covP.Mult(jacob,tmp);
1598 //__________________________________________________________________________
1599 UInt_t AliAnalysisTaskMuonResolution::BuildTriggerWord(const TString& firedTriggerClasses)
1601 /// build the trigger word from the fired trigger classes and the list of selectable trigger
1605 TObjString* trigClasseName = 0x0;
1606 TIter nextTrigger(fSelectTriggerClass);
1607 while ((trigClasseName = static_cast<TObjString*>(nextTrigger()))) {
1609 TRegexp GenericTriggerClasseName(trigClasseName->String());
1610 if (firedTriggerClasses.Contains(GenericTriggerClasseName)) word |= trigClasseName->GetUniqueID();
1617 //__________________________________________________________________________
1618 void AliAnalysisTaskMuonResolution::CheckPads(AliMUONVCluster *cl, Bool_t &hasBending, Bool_t &hasNonBending) const
1620 /// Check that this cluster contains pads on both cathods
1623 hasBending = kFALSE;
1624 hasNonBending = kFALSE;
1626 // loop over digits contained in the cluster
1627 for (Int_t iDigit = 0; iDigit < cl->GetNDigits(); iDigit++) {
1629 Int_t manuId = AliMUONVDigit::ManuId(cl->GetDigitId(iDigit));
1631 // check the location of the manu the digit belongs to
1633 if (manuId & AliMpConstants::ManuMask(AliMp::kNonBendingPlane)) hasNonBending = kTRUE;
1634 else hasBending = kTRUE;
1641 //________________________________________________________________________
1642 void AliAnalysisTaskMuonResolution::CheckPadsBelow(AliMUONVCluster *cl, Bool_t &hasBending, Bool_t &hasNonBending) const
1644 /// Check that this cluster contains pads on both cathods just under its position
1647 hasBending = kFALSE;
1648 hasNonBending = kFALSE;
1650 // get the cathod corresponding to the bending/non-bending plane
1651 Int_t deId = cl->GetDetElemId();
1652 AliMpDetElement* de = AliMpDDLStore::Instance()->GetDetElement(deId, kFALSE);
1654 AliMp::CathodType cath1 = de->GetCathodType(AliMp::kBendingPlane);
1655 AliMp::CathodType cath2 = de->GetCathodType(AliMp::kNonBendingPlane);
1657 // get the corresponding segmentation
1658 const AliMpVSegmentation* seg1 = AliMpSegmentation::Instance()->GetMpSegmentation(deId, cath1);
1659 const AliMpVSegmentation* seg2 = AliMpSegmentation::Instance()->GetMpSegmentation(deId, cath2);
1660 if (!seg1 || !seg2) return;
1662 // get local coordinate of the cluster
1664 Double_t gX = cl->GetX();
1665 Double_t gY = cl->GetY();
1666 Double_t gZ = cl->GetZ();
1667 fNewGeoTransformer->Global2Local(deId,gX,gY,gZ,lX,lY,lZ);
1669 // find pads below the cluster
1670 AliMpPad pad1 = seg1->PadByPosition(lX, lY, kFALSE);
1671 AliMpPad pad2 = seg2->PadByPosition(lX, lY, kFALSE);
1673 // build their ID if pads are valid
1674 UInt_t padId1 = (pad1.IsValid()) ? AliMUONVDigit::BuildUniqueID(deId, pad1.GetManuId(), pad1.GetManuChannel(), cath1) : 0;
1675 UInt_t padId2 = (pad2.IsValid()) ? AliMUONVDigit::BuildUniqueID(deId, pad2.GetManuId(), pad2.GetManuChannel(), cath2) : 0;
1677 // check if the cluster contains these pads
1678 for (Int_t iDigit = 0; iDigit < cl->GetNDigits(); iDigit++) {
1680 UInt_t digitId = cl->GetDigitId(iDigit);
1682 if (digitId == padId1) {
1685 if (hasNonBending) break;
1687 } else if (digitId == padId2) {
1689 hasNonBending = kTRUE;
1690 if (hasBending) break;