--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+// ROOT includes
+#include <TFile.h>
+#include <TTree.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TMultiGraph.h>
+#include <TGraphErrors.h>
+#include <TCanvas.h>
+#include <TLegend.h>
+#include <Riostream.h>
+#include <TString.h>
+#include <TGeoManager.h>
+
+// STEER includes
+#include "AliESDEvent.h"
+#include "AliESDMuonTrack.h"
+#include "AliCDBManager.h"
+#include "AliCDBStorage.h"
+#include "AliGeomManager.h"
+
+// ANALYSIS includes
+#include "AliAnalysisTaskSE.h"
+#include "AliAnalysisDataSlot.h"
+#include "AliAnalysisManager.h"
+#include "AliInputEventHandler.h"
+#include "AliAnalysisTaskMuonResolution.h"
+
+// MUON includes
+#include "AliMUONCDB.h"
+#include "AliMUONRecoParam.h"
+#include "AliMUONESDInterface.h"
+#include "AliMUONVTrackReconstructor.h"
+#include "AliMUONTrack.h"
+#include "AliMUONTrackParam.h"
+#include "AliMUONTrackExtrap.h"
+#include "AliMUONVCluster.h"
+#include "AliMUONGeometryTransformer.h"
+#include "AliMUONGeometryModuleTransformer.h"
+#include "AliMUONGeometryDetElement.h"
+#include "AliMpDEIterator.h"
+
+#ifndef SafeDelete
+#define SafeDelete(x) if (x != NULL) { delete x; x = NULL; }
+#endif
+
+ClassImp(AliAnalysisTaskMuonResolution)
+
+const Int_t AliAnalysisTaskMuonResolution::fgkMinEntries = 10;
+
+//________________________________________________________________________
+AliAnalysisTaskMuonResolution::AliAnalysisTaskMuonResolution() :
+ AliAnalysisTaskSE(),
+ fResiduals(NULL),
+ fResidualsVsP(NULL),
+ fLocalChi2(NULL),
+ fChamberRes(NULL),
+ fTrackRes(NULL),
+ fCanvases(NULL),
+ fNEvents(0),
+ fMinMomentum(0.),
+ fSelectPhysics(kFALSE),
+ fMatchTrig(kFALSE),
+ fExtrapMode(1),
+ fCorrectForSystematics(kTRUE),
+ fOCDBLoaded(kFALSE),
+ fNDE(0),
+ fReAlign(kFALSE),
+ fOldAlignStorage(""),
+ fNewAlignStorage(""),
+ fOldGeoTransformer(NULL),
+ fNewGeoTransformer(NULL)
+{
+ /// Default constructor
+}
+
+//________________________________________________________________________
+AliAnalysisTaskMuonResolution::AliAnalysisTaskMuonResolution(const char *name) :
+ AliAnalysisTaskSE(name),
+ fResiduals(NULL),
+ fResidualsVsP(NULL),
+ fLocalChi2(NULL),
+ fChamberRes(NULL),
+ fTrackRes(NULL),
+ fCanvases(NULL),
+ fNEvents(0),
+ fMinMomentum(0.),
+ fSelectPhysics(kFALSE),
+ fMatchTrig(kFALSE),
+ fExtrapMode(1),
+ fCorrectForSystematics(kTRUE),
+ fOCDBLoaded(kFALSE),
+ fNDE(0),
+ fReAlign(kFALSE),
+ fOldAlignStorage(""),
+ fNewAlignStorage(""),
+ fOldGeoTransformer(NULL),
+ fNewGeoTransformer(NULL)
+{
+ /// Constructor
+
+ for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) SetStartingResolution(i, -1., -1.);
+
+ // Output slot #1 writes into a TObjArray container
+ DefineOutput(1,TObjArray::Class());
+ // Output slot #2 writes into a TObjArray container
+ DefineOutput(2,TObjArray::Class());
+ // Output slot #3 writes into a TObjArray container
+ DefineOutput(3,TObjArray::Class());
+ // Output slot #4 writes into a TObjArray container
+ DefineOutput(4,TObjArray::Class());
+ // Output slot #5 writes into a TObjArray container
+ DefineOutput(5,TObjArray::Class());
+}
+
+//________________________________________________________________________
+AliAnalysisTaskMuonResolution::~AliAnalysisTaskMuonResolution()
+{
+ /// Destructor
+ SafeDelete(fResiduals);
+ SafeDelete(fResidualsVsP);
+ SafeDelete(fChamberRes);
+ SafeDelete(fTrackRes);
+ SafeDelete(fCanvases);
+ SafeDelete(fOldGeoTransformer);
+ SafeDelete(fNewGeoTransformer);
+}
+
+//___________________________________________________________________________
+void AliAnalysisTaskMuonResolution::UserCreateOutputObjects()
+{
+ /// Create histograms
+
+ // do it once the OCDB has been loaded (i.e. from NotifyRun())
+ if (!fOCDBLoaded) return;
+
+ fResiduals = new TObjArray(1000);
+ fResiduals->SetOwner();
+ fResidualsVsP = new TObjArray(1000);
+ fResidualsVsP->SetOwner();
+ fTrackRes = new TObjArray(1000);
+ fTrackRes->SetOwner();
+ TH2F* h2;
+
+ // find the highest chamber resolution and set histogram bins
+ const AliMUONRecoParam* recoParam = AliMUONESDInterface::GetTracker()->GetRecoParam();
+ Double_t maxSigma[2] = {-1., -1.};
+ for (Int_t i = 0; i < 10; i++) {
+ if (recoParam->GetDefaultNonBendingReso(i) > maxSigma[0]) maxSigma[0] = recoParam->GetDefaultNonBendingReso(i);
+ if (recoParam->GetDefaultBendingReso(i) > maxSigma[1]) maxSigma[1] = recoParam->GetDefaultBendingReso(i);
+ }
+ const char* axes[2] = {"X", "Y"};
+ const Int_t nBins = 2000;
+ const Int_t nSigma = 10;
+ const Int_t pNBins = 20;
+ const Double_t pEdges[2] = {0., 50.};
+
+ for (Int_t ia = 0; ia < 2; ia++) {
+
+ Double_t maxRes = nSigma*maxSigma[ia];
+
+ // List of residual histos
+ 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);
+ fResiduals->AddAtAndExpand(h2, kResidualPerCh_ClusterIn+ia);
+ 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);
+ fResiduals->AddAtAndExpand(h2, kResidualPerCh_ClusterOut+ia);
+
+ 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);
+ fResiduals->AddAtAndExpand(h2, kResidualPerHalfCh_ClusterIn+ia);
+ 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);
+ fResiduals->AddAtAndExpand(h2, kResidualPerHalfCh_ClusterOut+ia);
+
+ 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);
+ for (Int_t i = 1; i <= fNDE; i++) h2->GetXaxis()->SetBinLabel(i, Form("%d",fDEIds[i]));
+ fResiduals->AddAtAndExpand(h2, kResidualPerDE_ClusterIn+ia);
+ 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);
+ for (Int_t i = 1; i <= fNDE; i++) h2->GetXaxis()->SetBinLabel(i, Form("%d",fDEIds[i]));
+ fResiduals->AddAtAndExpand(h2, kResidualPerDE_ClusterOut+ia);
+
+ 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);
+ fResiduals->AddAtAndExpand(h2, kTrackResPerCh+ia);
+ 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);
+ fResiduals->AddAtAndExpand(h2, kTrackResPerHalfCh+ia);
+ 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);
+ for (Int_t i = 1; i <= fNDE; i++) h2->GetXaxis()->SetBinLabel(i, Form("%d",fDEIds[i]));
+ fResiduals->AddAtAndExpand(h2, kTrackResPerDE+ia);
+
+ 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);
+ fResiduals->AddAtAndExpand(h2, kMCSPerCh+ia);
+ 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);
+ fResiduals->AddAtAndExpand(h2, kMCSPerHalfCh+ia);
+ 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);
+ for (Int_t i = 1; i <= fNDE; i++) h2->GetXaxis()->SetBinLabel(i, Form("%d",fDEIds[i]));
+ fResiduals->AddAtAndExpand(h2, kMCSPerDE+ia);
+
+ 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);
+ fResiduals->AddAtAndExpand(h2, kClusterRes2PerCh+ia);
+
+ // List of residual vs. p histos
+ for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
+ 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);
+ fResidualsVsP->AddAtAndExpand(h2, kResidualInChVsP_ClusterIn+10*ia+i);
+ 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);
+ fResidualsVsP->AddAtAndExpand(h2, kResidualInChVsP_ClusterOut+10*ia+i);
+ }
+
+ // local chi2
+ 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.);
+ fResiduals->AddAtAndExpand(h2, kLocalChi2PerCh+ia);
+ 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.);
+ for (Int_t i = 1; i <= fNDE; i++) h2->GetXaxis()->SetBinLabel(i, Form("%d",fDEIds[i]));
+ fResiduals->AddAtAndExpand(h2, kLocalChi2PerDE+ia);
+
+ // track resolution
+ 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);
+ fTrackRes->AddAtAndExpand(h2, kUncorrSlopeRes+ia);
+ 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);
+ fTrackRes->AddAtAndExpand(h2, kSlopeRes+ia);
+ }
+
+ // local chi2 X+Y
+ 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.);
+ fResiduals->AddAtAndExpand(h2, kLocalChi2PerCh+2);
+ 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.);
+ for (Int_t i = 1; i <= fNDE; i++) h2->GetXaxis()->SetBinLabel(i, Form("%d",fDEIds[i]));
+ fResiduals->AddAtAndExpand(h2, kLocalChi2PerDE+2);
+
+ // track resolution
+ 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.);
+ fTrackRes->AddAtAndExpand(h2, kUncorrPRes);
+ h2 = new TH2F("hPRes", "muon momentum reconstructed resolution at vertex vs p;p (GeV/c); #sigma_{p}/p (%)", 300, 0., 300., 1000, 0., 10.);
+ fTrackRes->AddAtAndExpand(h2, kPRes);
+ 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.);
+ fTrackRes->AddAtAndExpand(h2, kUncorrPtRes);
+ 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.);
+ fTrackRes->AddAtAndExpand(h2, kPtRes);
+
+ // Post data at least once per task to ensure data synchronisation (required for merging)
+ PostData(1, fResiduals);
+ PostData(2, fResidualsVsP);
+ PostData(5, fTrackRes);
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskMuonResolution::UserExec(Option_t *)
+{
+ /// Main event loop
+
+ // check if OCDB properly loaded
+ if (!fOCDBLoaded) return;
+
+ AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
+ if (!esd) return;
+
+ if ((++fNEvents)%100 == 0) cout<<"\rEvent processing... "<<fNEvents<<"\r"<<flush;
+
+ // get tracker to refit
+ AliMUONVTrackReconstructor* tracker = AliMUONESDInterface::GetTracker();
+
+ // loop over tracks
+ Int_t nTracks = (Int_t) esd->GetNumberOfMuonTracks();
+ for (Int_t iTrack = 0; iTrack < nTracks; ++iTrack) {
+
+ // get the ESD track
+ AliESDMuonTrack* esdTrack = esd->GetMuonTrack(iTrack);
+
+ // skip ghost tracks
+ if (!esdTrack->ContainTrackerData()) continue;
+
+ // skip tracks that do not pass the physics selection if required
+ if (fSelectPhysics && fInputHandler && !fInputHandler->IsEventSelected()) continue;
+
+ // skip tracks not matched with trigger if required
+ if (fMatchTrig && !esdTrack->ContainTriggerData()) continue;
+
+ // skip low momentum tracks
+ if (esdTrack->PUncorrected() < fMinMomentum) continue;
+
+ // get the corresponding MUON track
+ AliMUONTrack track;
+ AliMUONESDInterface::ESDToMUON(*esdTrack, track, kFALSE);
+
+ // change the cluster resolution (and position)
+ ModifyClusters(track);
+
+ // refit the track
+ if (!tracker->RefitTrack(track, kFALSE)) break;
+
+ // save track unchanged
+ AliMUONTrack referenceTrack(track);
+
+ // get track param at first cluster and add MCS in first chamber
+ AliMUONTrackParam trackParamAtFirstCluster(*(static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->First())));
+ Int_t firstCh = 0; while (firstCh < 10 && !esdTrack->IsInMuonClusterMap(firstCh)) firstCh++;
+ AliMUONTrackExtrap::AddMCSEffect(&trackParamAtFirstCluster, AliMUONConstants::ChamberThicknessInX0(firstCh)/2., -1.);
+
+ // fill momentum error at first cluster
+ Double_t pXUncorr = trackParamAtFirstCluster.Px();
+ Double_t pYUncorr = trackParamAtFirstCluster.Py();
+ Double_t pZUncorr = trackParamAtFirstCluster.Pz();
+ Double_t pUncorr = trackParamAtFirstCluster.P();
+ TMatrixD covUncorr(5,5);
+ Cov2CovP(trackParamAtFirstCluster,covUncorr);
+ Double_t sigmaPUncorr = TMath::Sqrt(pXUncorr * (pXUncorr*covUncorr(2,2) + pYUncorr*covUncorr(2,3) + pZUncorr*covUncorr(2,4)) +
+ pYUncorr * (pXUncorr*covUncorr(3,2) + pYUncorr*covUncorr(3,3) + pZUncorr*covUncorr(3,4)) +
+ pZUncorr * (pXUncorr*covUncorr(4,2) + pYUncorr*covUncorr(4,3) + pZUncorr*covUncorr(4,4))) / pUncorr;
+ ((TH2F*)fTrackRes->UncheckedAt(kUncorrPRes))->Fill(pUncorr,100.*sigmaPUncorr/pUncorr);
+
+ // fill transverse momentum error at first cluster
+ Double_t ptUncorr = TMath::Sqrt(pXUncorr*pXUncorr + pYUncorr*pYUncorr);
+ Double_t sigmaPtUncorr = TMath::Sqrt(pXUncorr * (pXUncorr*covUncorr(2,2) + pYUncorr*covUncorr(2,3)) + pYUncorr * (pXUncorr*covUncorr(3,2) + pYUncorr*covUncorr(3,3))) / ptUncorr;
+ ((TH2F*)fTrackRes->UncheckedAt(kUncorrPtRes))->Fill(ptUncorr,100.*sigmaPtUncorr/ptUncorr);
+
+ // fill slopeX-Y error at first cluster
+ ((TH2F*)fTrackRes->UncheckedAt(kUncorrSlopeRes))->Fill(pUncorr,TMath::Sqrt(trackParamAtFirstCluster.GetCovariances()(1,1)));
+ ((TH2F*)fTrackRes->UncheckedAt(kUncorrSlopeRes+1))->Fill(pUncorr,TMath::Sqrt(trackParamAtFirstCluster.GetCovariances()(3,3)));
+
+ // fill momentum error at vertex
+ AliMUONTrackParam trackParamAtVtx(trackParamAtFirstCluster);
+ AliMUONTrackExtrap::ExtrapToVertex(&trackParamAtVtx, esdTrack->GetNonBendingCoor(), esdTrack->GetBendingCoor(), esdTrack->GetZ(), 0., 0.);
+ Double_t pXVtx = trackParamAtVtx.Px();
+ Double_t pYVtx = trackParamAtVtx.Py();
+ Double_t pZVtx = trackParamAtVtx.Pz();
+ Double_t pVtx = trackParamAtVtx.P();
+ TMatrixD covVtx(5,5);
+ Cov2CovP(trackParamAtVtx,covVtx);
+ Double_t sigmaPVtx = TMath::Sqrt(pXVtx * (pXVtx*covVtx(2,2) + pYVtx*covVtx(2,3) + pZVtx*covVtx(2,4)) +
+ pYVtx * (pXVtx*covVtx(3,2) + pYVtx*covVtx(3,3) + pZVtx*covVtx(3,4)) +
+ pZVtx * (pXVtx*covVtx(4,2) + pYVtx*covVtx(4,3) + pZVtx*covVtx(4,4))) / pVtx;
+ ((TH2F*)fTrackRes->UncheckedAt(kPRes))->Fill(pVtx,100.*sigmaPVtx/pVtx);
+
+ // fill transverse momentum error at vertex
+ Double_t ptVtx = TMath::Sqrt(pXVtx*pXVtx + pYVtx*pYVtx);
+ Double_t sigmaPtVtx = TMath::Sqrt(pXVtx * (pXVtx*covVtx(2,2) + pYVtx*covVtx(2,3)) + pYVtx * (pXVtx*covVtx(3,2) + pYVtx*covVtx(3,3))) / ptVtx;
+ ((TH2F*)fTrackRes->UncheckedAt(kPtRes))->Fill(ptVtx,100.*sigmaPtVtx/ptVtx);
+
+ // fill slopeX-Y error at vertex
+ ((TH2F*)fTrackRes->UncheckedAt(kSlopeRes))->Fill(pVtx,TMath::Sqrt(trackParamAtVtx.GetCovariances()(1,1)));
+ ((TH2F*)fTrackRes->UncheckedAt(kSlopeRes+1))->Fill(pVtx,TMath::Sqrt(trackParamAtVtx.GetCovariances()(3,3)));
+
+ // loop over clusters
+ Int_t nClusters = track.GetNClusters();
+ for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
+
+ // Get current, previous and next trackParams
+ AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->UncheckedAt(iCluster));
+ AliMUONTrackParam* previousTrackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->Before(trackParam));
+ AliMUONTrackParam* nextTrackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->After(trackParam));
+
+ // save current trackParam and remove it from the track
+ AliMUONTrackParam currentTrackParam(*trackParam);
+ track.RemoveTrackParamAtCluster(trackParam);
+
+ // get cluster info
+ AliMUONVCluster* cluster = currentTrackParam.GetClusterPtr();
+ Int_t chId = cluster->GetChamberId();
+ Int_t halfChId = (cluster->GetX() > 0) ? 2*chId : 2*chId+1;
+ Int_t deId = cluster->GetDetElemId();
+
+ // compute residuals with cluster still attached to the track
+ AliMUONTrackParam* referenceTrackParam = static_cast<AliMUONTrackParam*>(referenceTrack.GetTrackParamAtCluster()->UncheckedAt(iCluster));
+ Double_t deltaX = cluster->GetX() - referenceTrackParam->GetNonBendingCoor();
+ Double_t deltaY = cluster->GetY() - referenceTrackParam->GetBendingCoor();
+
+ // compute local chi2
+ Double_t sigmaDeltaX2 = cluster->GetErrX2() - referenceTrackParam->GetCovariances()(0,0);
+ Double_t sigmaDeltaY2 = cluster->GetErrY2() - referenceTrackParam->GetCovariances()(2,2);
+ Double_t localChi2X = (sigmaDeltaX2 > 0.) ? deltaX*deltaX/sigmaDeltaX2 : 0.;
+ Double_t localChi2Y = (sigmaDeltaY2 > 0.) ? deltaY*deltaY/sigmaDeltaY2 : 0.;
+ Double_t localChi2 = 0.5 * referenceTrackParam->GetLocalChi2();
+
+ // fill local chi2 info at every clusters
+ ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerCh))->Fill(chId+1, localChi2X);
+ ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerCh+1))->Fill(chId+1, localChi2Y);
+ ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerCh+2))->Fill(chId+1, localChi2);
+ ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerDE))->Fill(fDEIndices[deId], localChi2X);
+ ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerDE+1))->Fill(fDEIndices[deId], localChi2Y);
+ ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerDE+2))->Fill(fDEIndices[deId], localChi2);
+
+ // make sure the track has another cluster in the same station and can still be refitted
+ Bool_t refit = track.IsValid( 1 << (chId/2) );
+ if (refit) {
+
+ // refit the track and proceed if everything goes fine
+ if (tracker->RefitTrack(track, kFALSE)) {
+
+ // fill histograms of residuals with cluster still attached to the track
+ ((TH2F*)fResiduals->UncheckedAt(kResidualPerCh_ClusterIn))->Fill(chId+1, deltaX);
+ ((TH2F*)fResiduals->UncheckedAt(kResidualPerCh_ClusterIn+1))->Fill(chId+1, deltaY);
+ ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfCh_ClusterIn))->Fill(halfChId+1, deltaX);
+ ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfCh_ClusterIn+1))->Fill(halfChId+1, deltaY);
+ ((TH2F*)fResiduals->UncheckedAt(kResidualPerDE_ClusterIn))->Fill(fDEIndices[deId], deltaX);
+ ((TH2F*)fResiduals->UncheckedAt(kResidualPerDE_ClusterIn+1))->Fill(fDEIndices[deId], deltaY);
+ ((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsP_ClusterIn+chId))->Fill(pUncorr, deltaX);
+ ((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsP_ClusterIn+10+chId))->Fill(pUncorr, deltaY);
+
+ // find the track parameters closest to the current cluster position
+ Double_t dZWithPrevious = (previousTrackParam) ? TMath::Abs(previousTrackParam->GetClusterPtr()->GetZ() - cluster->GetZ()) : FLT_MAX;
+ Int_t previousChId = (previousTrackParam) ? previousTrackParam->GetClusterPtr()->GetChamberId() : -1;
+ Double_t dZWithNext = (nextTrackParam) ? TMath::Abs(nextTrackParam->GetClusterPtr()->GetZ() - cluster->GetZ()) : FLT_MAX;
+ AliMUONTrackParam* startingTrackParam = 0x0;
+ if ((fExtrapMode == 0 && dZWithPrevious < dZWithNext) ||
+ (fExtrapMode == 1 && previousTrackParam && !(chId/2 == 2 && previousChId/2 == 1) &&
+ !(chId/2 == 3 && previousChId/2 == 2))) startingTrackParam = previousTrackParam;
+ else startingTrackParam = nextTrackParam;
+
+ // reset current parameters
+ currentTrackParam.SetParameters(startingTrackParam->GetParameters());
+ currentTrackParam.SetZ(startingTrackParam->GetZ());
+ currentTrackParam.SetCovariances(startingTrackParam->GetCovariances());
+ currentTrackParam.ResetPropagator();
+
+ // extrapolate to the current cluster position and fill histograms of residuals if everything goes fine
+ if (AliMUONTrackExtrap::ExtrapToZCov(¤tTrackParam, currentTrackParam.GetClusterPtr()->GetZ(), kTRUE)) {
+
+ // compute MCS dispersion on the first chamber
+ TMatrixD mcsCov(5,5);
+ if (startingTrackParam == nextTrackParam && chId == 0) {
+ AliMUONTrackParam trackParamForMCS;
+ trackParamForMCS.SetParameters(nextTrackParam->GetParameters());
+ AliMUONTrackExtrap::AddMCSEffect(&trackParamForMCS,AliMUONConstants::ChamberThicknessInX0(nextTrackParam->GetClusterPtr()->GetChamberId()),-1.);
+ const TMatrixD &propagator = currentTrackParam.GetPropagator();
+ TMatrixD tmp(trackParamForMCS.GetCovariances(),TMatrixD::kMultTranspose,propagator);
+ mcsCov.Mult(propagator,tmp);
+ } else mcsCov.Zero();
+
+ // compute residuals
+ Double_t trackResX2 = currentTrackParam.GetCovariances()(0,0) + mcsCov(0,0);
+ Double_t trackResY2 = currentTrackParam.GetCovariances()(2,2) + mcsCov(2,2);
+ deltaX = cluster->GetX() - currentTrackParam.GetNonBendingCoor();
+ deltaY = cluster->GetY() - currentTrackParam.GetBendingCoor();
+
+ // fill histograms with cluster not attached to the track
+ ((TH2F*)fResiduals->UncheckedAt(kResidualPerCh_ClusterOut))->Fill(chId+1, deltaX);
+ ((TH2F*)fResiduals->UncheckedAt(kResidualPerCh_ClusterOut+1))->Fill(chId+1, deltaY);
+ ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfCh_ClusterOut))->Fill(halfChId+1, deltaX);
+ ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfCh_ClusterOut+1))->Fill(halfChId+1, deltaY);
+ ((TH2F*)fResiduals->UncheckedAt(kResidualPerDE_ClusterOut))->Fill(fDEIndices[deId], deltaX);
+ ((TH2F*)fResiduals->UncheckedAt(kResidualPerDE_ClusterOut+1))->Fill(fDEIndices[deId], deltaY);
+ ((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsP_ClusterOut+chId))->Fill(pUncorr, deltaX);
+ ((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsP_ClusterOut+10+chId))->Fill(pUncorr, deltaY);
+ ((TH2F*)fResiduals->UncheckedAt(kTrackResPerCh))->Fill(chId+1, TMath::Sqrt(trackResX2));
+ ((TH2F*)fResiduals->UncheckedAt(kTrackResPerCh+1))->Fill(chId+1, TMath::Sqrt(trackResY2));
+ ((TH2F*)fResiduals->UncheckedAt(kTrackResPerHalfCh))->Fill(halfChId+1, TMath::Sqrt(trackResX2));
+ ((TH2F*)fResiduals->UncheckedAt(kTrackResPerHalfCh+1))->Fill(halfChId+1, TMath::Sqrt(trackResY2));
+ ((TH2F*)fResiduals->UncheckedAt(kTrackResPerDE))->Fill(fDEIndices[deId], TMath::Sqrt(trackResX2));
+ ((TH2F*)fResiduals->UncheckedAt(kTrackResPerDE+1))->Fill(fDEIndices[deId], TMath::Sqrt(trackResY2));
+ ((TH2F*)fResiduals->UncheckedAt(kMCSPerCh))->Fill(chId+1, TMath::Sqrt(mcsCov(0,0)));
+ ((TH2F*)fResiduals->UncheckedAt(kMCSPerCh+1))->Fill(chId+1, TMath::Sqrt(mcsCov(2,2)));
+ ((TH2F*)fResiduals->UncheckedAt(kMCSPerHalfCh))->Fill(halfChId+1, TMath::Sqrt(mcsCov(0,0)));
+ ((TH2F*)fResiduals->UncheckedAt(kMCSPerHalfCh+1))->Fill(halfChId+1, TMath::Sqrt(mcsCov(2,2)));
+ ((TH2F*)fResiduals->UncheckedAt(kMCSPerDE))->Fill(fDEIndices[deId], TMath::Sqrt(mcsCov(0,0)));
+ ((TH2F*)fResiduals->UncheckedAt(kMCSPerDE+1))->Fill(fDEIndices[deId], TMath::Sqrt(mcsCov(2,2)));
+ ((TH2F*)fResiduals->UncheckedAt(kClusterRes2PerCh))->Fill(chId+1, deltaX*deltaX - trackResX2);
+ ((TH2F*)fResiduals->UncheckedAt(kClusterRes2PerCh+1))->Fill(chId+1, deltaY*deltaY - trackResY2);
+ }
+
+ }
+
+ }
+
+ // restore the track
+ track.AddTrackParamAtCluster(currentTrackParam, *(currentTrackParam.GetClusterPtr()), kTRUE);
+
+ }
+
+ }
+
+ // Post final data. It will be written to a file with option "RECREATE"
+ PostData(1, fResiduals);
+ PostData(2, fResidualsVsP);
+ PostData(5, fTrackRes);
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskMuonResolution::NotifyRun()
+{
+ /// load necessary data from OCDB corresponding to the first run number and initialize analysis
+
+ if (fOCDBLoaded) return;
+
+ AliCDBManager* cdbm = AliCDBManager::Instance();
+ cdbm->SetRun(fCurrentRunNumber);
+
+ if (!AliMUONCDB::LoadField()) return;
+
+ if (!AliMUONCDB::LoadMapping()) return;
+
+ AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam();
+ if (!recoParam) return;
+
+ fOCDBLoaded = kTRUE;
+
+ AliMUONESDInterface::ResetTracker(recoParam);
+
+ for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
+
+ // set the cluster resolution to default if not already set and create output objects
+ if (fClusterResNB[i] < 0.) fClusterResNB[i] = recoParam->GetDefaultNonBendingReso(i);
+ if (fClusterResB[i] < 0.) fClusterResB[i] = recoParam->GetDefaultBendingReso(i);
+
+ // fill correspondence between DEId and indices for histo (starting from 1)
+ AliMpDEIterator it;
+ it.First(i);
+ while (!it.IsDone()) {
+ fNDE++;
+ fDEIndices[it.CurrentDEId()] = fNDE;
+ fDEIds[fNDE] = it.CurrentDEId();
+ it.Next();
+ }
+
+ }
+
+ if (fReAlign) {
+
+ // recover default storage name
+ TString defaultStorage(cdbm->GetDefaultStorage()->GetType());
+ if (defaultStorage == "alien") defaultStorage += Form("://folder=%s", cdbm->GetDefaultStorage()->GetBaseFolder().Data());
+ else defaultStorage += Form("://%s", cdbm->GetDefaultStorage()->GetBaseFolder().Data());
+
+ // reset existing geometry/alignment if any
+ if (cdbm->GetEntryCache()->Contains("GRP/Geometry/Data")) cdbm->UnloadFromCache("GRP/Geometry/Data");
+ if (cdbm->GetEntryCache()->Contains("MUON/Align/Data")) cdbm->UnloadFromCache("MUON/Align/Data");
+ if (AliGeomManager::GetGeometry()) AliGeomManager::GetGeometry()->UnlockGeometry();
+
+ // get original geometry transformer
+ AliGeomManager::LoadGeometry();
+ if (!AliGeomManager::GetGeometry()) return;
+ if (fOldAlignStorage != "none") {
+ if (!fOldAlignStorage.IsNull()) cdbm->SetSpecificStorage("MUON/Align/Data",fOldAlignStorage.Data());
+ else cdbm->SetSpecificStorage("MUON/Align/Data",defaultStorage.Data());
+ AliGeomManager::ApplyAlignObjsFromCDB("MUON");
+ }
+ fOldGeoTransformer = new AliMUONGeometryTransformer();
+ fOldGeoTransformer->LoadGeometryData();
+
+ // get new geometry transformer
+ cdbm->UnloadFromCache("GRP/Geometry/Data");
+ if (fOldAlignStorage != "none") cdbm->UnloadFromCache("MUON/Align/Data");
+ AliGeomManager::GetGeometry()->UnlockGeometry();
+ AliGeomManager::LoadGeometry();
+ if (!AliGeomManager::GetGeometry()) return;
+ if (!fNewAlignStorage.IsNull()) cdbm->SetSpecificStorage("MUON/Align/Data",fNewAlignStorage.Data());
+ else cdbm->SetSpecificStorage("MUON/Align/Data",defaultStorage.Data());
+ AliGeomManager::ApplyAlignObjsFromCDB("MUON");
+ fNewGeoTransformer = new AliMUONGeometryTransformer();
+ fNewGeoTransformer->LoadGeometryData();
+
+ } else {
+
+ // load geometry for track extrapolation to vertex
+ if (cdbm->GetEntryCache()->Contains("GRP/Geometry/Data")) cdbm->UnloadFromCache("GRP/Geometry/Data");
+ if (AliGeomManager::GetGeometry()) AliGeomManager::GetGeometry()->UnlockGeometry();
+ AliGeomManager::LoadGeometry();
+ if (!AliGeomManager::GetGeometry()) return;
+
+ }
+
+ // print starting chamber resolution
+ printf("\nstarting chamber resolution:\n");
+ printf(" - non-bending:");
+ for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %5.3f":", %5.3f",fClusterResNB[i]);
+ printf("\n - bending:");
+ for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %6.4f":", %6.4f",fClusterResB[i]);
+ printf("\n\n");
+
+ UserCreateOutputObjects();
+
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskMuonResolution::Terminate(Option_t *)
+{
+ /// compute final results
+
+ // recover output objects
+ fResiduals = static_cast<TObjArray*> (GetOutputData(1));
+ fResidualsVsP = static_cast<TObjArray*> (GetOutputData(2));
+ fTrackRes = static_cast<TObjArray*> (GetOutputData(5));
+ if (!fResiduals || !fResidualsVsP || !fTrackRes) return;
+
+ // summary graphs
+ fLocalChi2 = new TObjArray(1000);
+ fLocalChi2->SetOwner();
+ fChamberRes = new TObjArray(1000);
+ fChamberRes->SetOwner();
+ TGraphErrors* g;
+ TMultiGraph* mg;
+
+ const char* axes[2] = {"X", "Y"};
+ Double_t newClusterRes[2][10], newClusterResErr[2][10];
+ fNDE = ((TH2F*)fResiduals->UncheckedAt(kResidualPerDE_ClusterIn))->GetXaxis()->GetNbins();
+
+ for (Int_t ia = 0; ia < 2; ia++) {
+
+ g = new TGraphErrors(AliMUONConstants::NTrackingCh());
+ g->SetName(Form("gResidual%sPerChMean_ClusterIn",axes[ia]));
+ g->SetTitle(Form("cluster-track residual-%s per Ch: mean (cluster in);chamber ID;<#Delta_{%s}> (cm)",axes[ia],axes[ia]));
+ g->SetMarkerStyle(kFullDotLarge);
+ fChamberRes->AddAtAndExpand(g, kResidualPerChMean_ClusterIn+ia);
+ g = new TGraphErrors(AliMUONConstants::NTrackingCh());
+ g->SetName(Form("gResidual%sPerChMean_ClusterOut",axes[ia]));
+ g->SetTitle(Form("cluster-track residual-%s per Ch: mean (cluster out);chamber ID;<#Delta_{%s}> (cm)",axes[ia],axes[ia]));
+ g->SetMarkerStyle(kFullDotLarge);
+ fChamberRes->AddAtAndExpand(g, kResidualPerChMean_ClusterOut+ia);
+
+ g = new TGraphErrors(2*AliMUONConstants::NTrackingCh());
+ g->SetName(Form("gResidual%sPerHalfChMean_ClusterIn",axes[ia]));
+ g->SetTitle(Form("cluster-track residual-%s per half Ch: mean (cluster in);half chamber ID;<#Delta_{%s}> (cm)",axes[ia],axes[ia]));
+ g->SetMarkerStyle(kFullDotLarge);
+ fChamberRes->AddAtAndExpand(g, kResidualPerHalfChMean_ClusterIn+ia);
+ g = new TGraphErrors(2*AliMUONConstants::NTrackingCh());
+ g->SetName(Form("gResidual%sPerHalfChMean_ClusterOut",axes[ia]));
+ g->SetTitle(Form("cluster-track residual-%s per half Ch: mean (cluster out);half chamber ID;<#Delta_{%s}> (cm)",axes[ia],axes[ia]));
+ g->SetMarkerStyle(kFullDotLarge);
+ fChamberRes->AddAtAndExpand(g, kResidualPerHalfChMean_ClusterOut+ia);
+
+ g = new TGraphErrors(fNDE);
+ g->SetName(Form("gResidual%sPerDEMean_ClusterIn",axes[ia]));
+ g->SetTitle(Form("cluster-track residual-%s per DE: mean (cluster in);DE ID;<#Delta_{%s}> (cm)",axes[ia],axes[ia]));
+ g->SetMarkerStyle(kFullDotLarge);
+ fChamberRes->AddAtAndExpand(g, kResidualPerDEMean_ClusterIn+ia);
+ g = new TGraphErrors(fNDE);
+ g->SetName(Form("gResidual%sPerDEMean_ClusterOut",axes[ia]));
+ g->SetTitle(Form("cluster-track residual-%s per DE: mean (cluster out);DE ID;<#Delta_{%s}> (cm)",axes[ia],axes[ia]));
+ g->SetMarkerStyle(kFullDotLarge);
+ fChamberRes->AddAtAndExpand(g, kResidualPerDEMean_ClusterOut+ia);
+
+ g = new TGraphErrors(AliMUONConstants::NTrackingCh());
+ g->SetName(Form("gResidual%sPerChSigma_ClusterIn",axes[ia]));
+ g->SetTitle(Form("cluster-track residual-%s per Ch: sigma (cluster in);chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
+ g->SetMarkerStyle(kFullDotLarge);
+ fChamberRes->AddAtAndExpand(g, kResidualPerChSigma_ClusterIn+ia);
+ g = new TGraphErrors(AliMUONConstants::NTrackingCh());
+ g->SetName(Form("gResidual%sPerChSigma_ClusterOut",axes[ia]));
+ g->SetTitle(Form("cluster-track residual-%s per Ch: sigma (cluster out);chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
+ g->SetMarkerStyle(kFullDotLarge);
+ fChamberRes->AddAtAndExpand(g, kResidualPerChSigma_ClusterOut+ia);
+
+ g = new TGraphErrors(AliMUONConstants::NTrackingCh());
+ g->SetName(Form("gResidual%sPerChDispersion_ClusterOut",axes[ia]));
+ g->SetTitle(Form("cluster-track residual-%s per Ch: dispersion (cluster out);chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
+ g->SetMarkerStyle(kFullDotLarge);
+ fChamberRes->AddAtAndExpand(g, kResidualPerChDispersion_ClusterOut+ia);
+
+ g = new TGraphErrors(AliMUONConstants::NTrackingCh());
+ g->SetName(Form("gCombinedResidual%sPerChSigma",axes[ia]));
+ g->SetTitle(Form("combined cluster-track residual-%s per Ch: sigma;chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
+ g->SetMarkerStyle(kFullDotLarge);
+ fChamberRes->AddAtAndExpand(g, kCombinedResidualPerChSigma+ia);
+
+ g = new TGraphErrors(2*AliMUONConstants::NTrackingCh());
+ g->SetName(Form("gCombinedResidual%sPerHalfChSigma",axes[ia]));
+ g->SetTitle(Form("combined cluster-track residual-%s per half Ch: sigma;half chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
+ g->SetMarkerStyle(kFullDotLarge);
+ fChamberRes->AddAtAndExpand(g, kCombinedResidualPerHalfChSigma+ia);
+
+ g = new TGraphErrors(fNDE);
+ g->SetName(Form("gCombinedResidual%sPerDESigma",axes[ia]));
+ g->SetTitle(Form("combined cluster-track residual-%s per DE: sigma;DE ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
+ g->SetMarkerStyle(kFullDotLarge);
+ fChamberRes->AddAtAndExpand(g, kCombinedResidualPerDESigma+ia);
+
+ 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]));
+ for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
+ g = new TGraphErrors(((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsP_ClusterIn+10*ia+i))->GetNbinsX());
+ g->SetName(Form("gRes%sVsP_ch%d",axes[ia],i+1));
+ g->SetMarkerStyle(kFullDotMedium);
+ g->SetMarkerColor(i+1+i/9);
+ mg->Add(g,"p");
+ }
+ fChamberRes->AddAtAndExpand(mg, kCombinedResidualSigmaVsP+ia);
+
+ g = new TGraphErrors(AliMUONConstants::NTrackingCh());
+ g->SetName(Form("gTrackRes%sPerCh",axes[ia]));
+ g->SetTitle(Form("track <#sigma_{%s}> per Ch;chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
+ g->SetMarkerStyle(kFullDotLarge);
+ fChamberRes->AddAtAndExpand(g, kTrackResPerChMean+ia);
+
+ g = new TGraphErrors(AliMUONConstants::NTrackingCh());
+ g->SetName(Form("gMCS%sPerCh",axes[ia]));
+ g->SetTitle(Form("MCS %s-dispersion of extrapolated track per Ch;chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
+ g->SetMarkerStyle(kFullDotLarge);
+ fChamberRes->AddAtAndExpand(g, kMCSPerChMean+ia);
+
+ g = new TGraphErrors(AliMUONConstants::NTrackingCh());
+ g->SetName(Form("gClusterRes%sPerCh",axes[ia]));
+ g->SetTitle(Form("cluster <#sigma_{%s}> per Ch;chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
+ g->SetMarkerStyle(kFullDotLarge);
+ fChamberRes->AddAtAndExpand(g, kClusterResPerCh+ia);
+
+ g = new TGraphErrors(2*AliMUONConstants::NTrackingCh());
+ g->SetName(Form("gClusterRes%sPerHalfCh",axes[ia]));
+ g->SetTitle(Form("cluster <#sigma_{%s}> per half Ch;half chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
+ g->SetMarkerStyle(kFullDotLarge);
+ fChamberRes->AddAtAndExpand(g, kClusterResPerHalfCh+ia);
+
+ g = new TGraphErrors(fNDE);
+ g->SetName(Form("gClusterRes%sPerDE",axes[ia]));
+ g->SetTitle(Form("cluster <#sigma_{%s}> per DE;DE ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
+ g->SetMarkerStyle(kFullDotLarge);
+ fChamberRes->AddAtAndExpand(g, kClusterResPerDE+ia);
+
+ g = new TGraphErrors(AliMUONConstants::NTrackingCh());
+ g->SetName(Form("gCalcClusterRes%sPerCh",axes[ia]));
+ g->SetTitle(Form("calculated cluster <#sigma_{%s}> per Ch;chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
+ g->SetMarkerStyle(kFullDotLarge);
+ fChamberRes->AddAtAndExpand(g, kCalcClusterResPerCh+ia);
+
+ g = new TGraphErrors(AliMUONConstants::NTrackingCh());
+ g->SetName(Form("gLocalChi2%sPerChMean",axes[ia]));
+ g->SetTitle(Form("local chi2-%s per Ch: mean;chamber ID;<local #chi^{2}_{%s}>",axes[ia],axes[ia]));
+ g->SetMarkerStyle(kFullDotLarge);
+ fLocalChi2->AddAtAndExpand(g, kLocalChi2PerChMean+ia);
+
+ g = new TGraphErrors(fNDE);
+ g->SetName(Form("gLocalChi2%sPerDEMean",axes[ia]));
+ g->SetTitle(Form("local chi2-%s per DE: mean;DE ID;<local #chi^{2}_{%s}>",axes[ia],axes[ia]));
+ g->SetMarkerStyle(kFullDotLarge);
+ fLocalChi2->AddAtAndExpand(g, kLocalChi2PerDEMean+ia);
+
+ // compute residual mean and dispersion and averaged local chi2 per chamber and half chamber
+ Double_t meanIn, meanInErr, meanOut, meanOutErr, sigmaIn, sigmaInErr, sigmaOut, sigmaOutErr;
+ Double_t sigmaTrack, sigmaTrackErr, sigmaMCS, sigmaMCSErr, clusterRes, clusterResErr, sigmaCluster, sigmaClusterErr;
+ for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
+
+ // method 1
+ TH1D *tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerCh_ClusterIn+ia))->ProjectionY("tmp",i+1,i+1,"e");
+ GetMean(tmp, meanIn, meanInErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChMean_ClusterIn+ia), i, i+1);
+ GetRMS(tmp, sigmaIn, sigmaInErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChSigma_ClusterIn+ia), i, i+1);
+ delete tmp;
+
+ tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerCh_ClusterOut+ia))->ProjectionY("tmp",i+1,i+1,"e");
+ GetMean(tmp, meanOut, meanOutErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChMean_ClusterOut+ia), i, i+1);
+ GetRMS(tmp, sigmaOut, sigmaOutErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChSigma_ClusterOut+ia), i, i+1);
+ delete tmp;
+
+ if (fCorrectForSystematics) {
+ sigmaIn = TMath::Sqrt(sigmaIn*sigmaIn + meanIn*meanIn);
+ sigmaInErr = (sigmaIn>0) ? TMath::Sqrt(sigmaIn*sigmaIn*sigmaInErr*sigmaInErr + meanIn*meanIn*meanInErr*meanInErr) / sigmaIn : 0.;
+ sigmaOut = TMath::Sqrt(sigmaOut*sigmaOut + meanOut*meanOut);
+ sigmaOutErr = (sigmaOut>0) ? TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + meanOut*meanOut*meanOutErr*meanOutErr) / sigmaOut : 0.;
+ }
+ ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChDispersion_ClusterOut+ia))->SetPoint(i, i+1, sigmaOut);
+ ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChDispersion_ClusterOut+ia))->SetPointError(i, 0., sigmaOutErr);
+
+ clusterRes = TMath::Sqrt(sigmaIn*sigmaOut);
+ // clusterResErr = (clusterRes > 0.) ? 0.5 * TMath::Sqrt(sigmaInErr*sigmaInErr*sigmaOut*sigmaOut + sigmaIn*sigmaIn*sigmaOutErr*sigmaOutErr) / clusterRes : 0.;
+ clusterResErr = TMath::Sqrt(sigmaInErr*sigmaOutErr);
+ ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerChSigma+ia))->SetPoint(i, i+1, clusterRes);
+ ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerChSigma+ia))->SetPointError(i, 0., clusterResErr);
+ newClusterRes[ia][i] = clusterRes;
+ newClusterResErr[ia][i] = clusterResErr;
+
+ // method 2
+ tmp = ((TH2F*)fResiduals->UncheckedAt(kTrackResPerCh+ia))->ProjectionY("tmp",i+1,i+1,"e");
+ GetMean(tmp, sigmaTrack, sigmaTrackErr, (TGraphErrors*)fChamberRes->UncheckedAt(kTrackResPerChMean+ia), i, i+1, kFALSE);
+ delete tmp;
+
+ tmp = ((TH2F*)fResiduals->UncheckedAt(kMCSPerCh+ia))->ProjectionY("tmp",i+1,i+1,"e");
+ GetMean(tmp, sigmaMCS, sigmaMCSErr, (TGraphErrors*)fChamberRes->UncheckedAt(kMCSPerChMean+ia), i, i+1, kFALSE);
+ delete tmp;
+
+ sigmaCluster = sigmaOut*sigmaOut - sigmaTrack*sigmaTrack;
+ if (sigmaCluster > 0.) {
+ sigmaCluster = TMath::Sqrt(sigmaCluster);
+ sigmaClusterErr = TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + sigmaTrack*sigmaTrack*sigmaTrackErr*sigmaTrackErr) / sigmaCluster;
+ } else {
+ sigmaCluster = 0.;
+ sigmaClusterErr = 0.;
+ }
+ ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerCh+ia))->SetPoint(i, i+1, sigmaCluster);
+ ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerCh+ia))->SetPointError(i, 0., sigmaClusterErr);
+
+ // method 3
+ tmp = ((TH2F*)fResiduals->UncheckedAt(kClusterRes2PerCh+ia))->ProjectionY("tmp",i+1,i+1,"e");
+ ZoomRight(tmp);
+ clusterRes = tmp->GetMean();
+ if (clusterRes > 0.) {
+ ((TGraphErrors*)fChamberRes->UncheckedAt(kCalcClusterResPerCh+ia))->SetPoint(i, i+1, TMath::Sqrt(clusterRes));
+ ((TGraphErrors*)fChamberRes->UncheckedAt(kCalcClusterResPerCh+ia))->SetPointError(i, 0., 0.5 * tmp->GetMeanError() / TMath::Sqrt(clusterRes));
+ } else {
+ ((TGraphErrors*)fChamberRes->UncheckedAt(kCalcClusterResPerCh+ia))->SetPoint(i, i+1, 0.);
+ ((TGraphErrors*)fChamberRes->UncheckedAt(kCalcClusterResPerCh+ia))->SetPointError(i, 0., 0.);
+ }
+ delete tmp;
+
+ // method 1 versus p
+ FillSigmaClusterVsP((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsP_ClusterIn+10*ia+i),
+ (TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsP_ClusterOut+10*ia+i),
+ (TGraphErrors*)((TMultiGraph*)fChamberRes->UncheckedAt(kCombinedResidualSigmaVsP+ia))->GetListOfGraphs()->FindObject(Form("gRes%sVsP_ch%d",axes[ia],i+1)));
+
+ // compute residual mean and dispersion per half chamber
+ for (Int_t j = 0; j < 2; j++) {
+ Int_t k = 2*i+j;
+
+ // method 1
+ tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfCh_ClusterIn+ia))->ProjectionY("tmp",k+1,k+1,"e");
+ GetMean(tmp, meanIn, meanInErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMean_ClusterIn+ia), k, k+1);
+ GetRMS(tmp, sigmaIn, sigmaInErr);
+ delete tmp;
+
+ tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfCh_ClusterOut+ia))->ProjectionY("tmp",k+1,k+1,"e");
+ GetMean(tmp, meanOut, meanOutErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMean_ClusterOut+ia), k, k+1);
+ GetRMS(tmp, sigmaOut, sigmaOutErr);
+ delete tmp;
+
+ if (fCorrectForSystematics) {
+ sigmaIn = TMath::Sqrt(sigmaIn*sigmaIn + meanIn*meanIn);
+ sigmaInErr = (sigmaIn>0) ? TMath::Sqrt(sigmaIn*sigmaIn*sigmaInErr*sigmaInErr + meanIn*meanIn*meanInErr*meanInErr) / sigmaIn : 0.;
+ sigmaOut = TMath::Sqrt(sigmaOut*sigmaOut + meanOut*meanOut);
+ sigmaOutErr = (sigmaOut>0) ? TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + meanOut*meanOut*meanOutErr*meanOutErr) / sigmaOut : 0.;
+ }
+
+ clusterRes = TMath::Sqrt(sigmaIn*sigmaOut);
+ // clusterResErr = (clusterRes > 0.) ? 0.5 * TMath::Sqrt(sigmaInErr*sigmaInErr*sigmaOut*sigmaOut + sigmaIn*sigmaIn*sigmaOutErr*sigmaOutErr) / clusterRes : 0.;
+ clusterResErr = TMath::Sqrt(sigmaInErr*sigmaOutErr);
+ ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerHalfChSigma+ia))->SetPoint(k, k+1, clusterRes);
+ ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerHalfChSigma+ia))->SetPointError(k, 0., clusterResErr);
+
+ // method 2
+ tmp = ((TH2F*)fResiduals->UncheckedAt(kTrackResPerHalfCh+ia))->ProjectionY("tmp",k+1,k+1,"e");
+ GetMean(tmp, sigmaTrack, sigmaTrackErr, 0x0, 0, 0, kFALSE);
+ delete tmp;
+
+ tmp = ((TH2F*)fResiduals->UncheckedAt(kMCSPerHalfCh+ia))->ProjectionY("tmp",k+1,k+1,"e");
+ GetMean(tmp, sigmaMCS, sigmaMCSErr, 0x0, 0, 0, kFALSE);
+ delete tmp;
+
+ sigmaCluster = sigmaOut*sigmaOut - sigmaTrack*sigmaTrack;
+ if (sigmaCluster > 0.) {
+ sigmaCluster = TMath::Sqrt(sigmaCluster);
+ sigmaClusterErr = TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + sigmaTrack*sigmaTrack*sigmaTrackErr*sigmaTrackErr) / sigmaCluster;
+ } else {
+ sigmaCluster = 0.;
+ sigmaClusterErr = 0.;
+ }
+ ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerHalfCh+ia))->SetPoint(k, k+1, sigmaCluster);
+ ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerHalfCh+ia))->SetPointError(k, 0., sigmaClusterErr);
+
+ }
+
+ // compute averaged local chi2
+ tmp = ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerCh+ia))->ProjectionY("tmp",i+1,i+1,"e");
+ ((TGraphErrors*)fLocalChi2->UncheckedAt(kLocalChi2PerChMean+ia))->SetPoint(i, i+1, tmp->GetMean());
+ ((TGraphErrors*)fLocalChi2->UncheckedAt(kLocalChi2PerChMean+ia))->SetPointError(i, 0., tmp->GetMeanError());
+ delete tmp;
+
+ }
+
+ // compute residual mean and dispersion per DE
+ for (Int_t i = 0; i < fNDE; i++) {
+
+ // method 1
+ TH1D *tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerDE_ClusterIn+ia))->ProjectionY("tmp",i+1,i+1,"e");
+ GetMean(tmp, meanIn, meanInErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMean_ClusterIn+ia), i, i+1);
+ GetRMS(tmp, sigmaIn, sigmaInErr);
+ delete tmp;
+
+ tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerDE_ClusterOut+ia))->ProjectionY("tmp",i+1,i+1,"e");
+ GetMean(tmp, meanOut, meanOutErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMean_ClusterOut+ia), i, i+1);
+ GetRMS(tmp, sigmaOut, sigmaOutErr);
+ delete tmp;
+
+ if (fCorrectForSystematics) {
+ sigmaIn = TMath::Sqrt(sigmaIn*sigmaIn + meanIn*meanIn);
+ sigmaInErr = (sigmaIn>0) ? TMath::Sqrt(sigmaIn*sigmaIn*sigmaInErr*sigmaInErr + meanIn*meanIn*meanInErr*meanInErr) / sigmaIn : 0.;
+ sigmaOut = TMath::Sqrt(sigmaOut*sigmaOut + meanOut*meanOut);
+ sigmaOutErr = (sigmaOut>0) ? TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + meanOut*meanOut*meanOutErr*meanOutErr) / sigmaOut : 0.;
+ }
+
+ clusterRes = TMath::Sqrt(sigmaIn*sigmaOut);
+ // clusterResErr = (clusterRes > 0.) ? 0.5 * TMath::Sqrt(sigmaInErr*sigmaInErr*sigmaOut*sigmaOut + sigmaIn*sigmaIn*sigmaOutErr*sigmaOutErr) / clusterRes : 0.;
+ clusterResErr = TMath::Sqrt(sigmaInErr*sigmaOutErr);
+ ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerDESigma+ia))->SetPoint(i, i+1, clusterRes);
+ ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerDESigma+ia))->SetPointError(i, 0., clusterResErr);
+
+ // method 2
+ tmp = ((TH2F*)fResiduals->UncheckedAt(kTrackResPerDE+ia))->ProjectionY("tmp",i+1,i+1,"e");
+ GetMean(tmp, sigmaTrack, sigmaTrackErr, 0x0, 0, 0, kFALSE);
+ delete tmp;
+
+ tmp = ((TH2F*)fResiduals->UncheckedAt(kMCSPerDE+ia))->ProjectionY("tmp",i+1,i+1,"e");
+ GetMean(tmp, sigmaMCS, sigmaMCSErr, 0x0, 0, 0, kFALSE);
+ delete tmp;
+
+ sigmaCluster = sigmaOut*sigmaOut - sigmaTrack*sigmaTrack;
+ if (sigmaCluster > 0.) {
+ sigmaCluster = TMath::Sqrt(sigmaCluster);
+ sigmaClusterErr = TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + sigmaTrack*sigmaTrack*sigmaTrackErr*sigmaTrackErr) / sigmaCluster;
+ } else {
+ sigmaCluster = 0.;
+ sigmaClusterErr = 0.;
+ }
+ ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerDE+ia))->SetPoint(i, i+1, sigmaCluster);
+ ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerDE+ia))->SetPointError(i, 0., sigmaClusterErr);
+
+ // compute averaged local chi2
+ tmp = ((TH2F*)fResiduals->UncheckedAt(kLocalChi2PerDE+ia))->ProjectionY("tmp",i+1,i+1,"e");
+ ((TGraphErrors*)fLocalChi2->UncheckedAt(kLocalChi2PerDEMean+ia))->SetPoint(i, i+1, tmp->GetMean());
+ ((TGraphErrors*)fLocalChi2->UncheckedAt(kLocalChi2PerDEMean+ia))->SetPointError(i, 0., tmp->GetMeanError());
+ delete tmp;
+
+ }
+
+ // set graph labels
+ TAxis* xAxis = ((TH2F*)fResiduals->UncheckedAt(kResidualPerDE_ClusterOut+ia))->GetXaxis();
+ ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMean_ClusterIn+ia))->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
+ ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMean_ClusterOut+ia))->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
+ ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerDESigma+ia))->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
+ ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerDE+ia))->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
+ ((TGraphErrors*)fLocalChi2->UncheckedAt(kLocalChi2PerDEMean+ia))->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
+ for (Int_t i = 1; i <= fNDE; i++) {
+ const char* label = xAxis->GetBinLabel(i);
+ ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMean_ClusterIn+ia))->GetXaxis()->SetBinLabel(i, label);
+ ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMean_ClusterOut+ia))->GetXaxis()->SetBinLabel(i, label);
+ ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerDESigma+ia))->GetXaxis()->SetBinLabel(i, label);
+ ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerDE+ia))->GetXaxis()->SetBinLabel(i, label);
+ ((TGraphErrors*)fLocalChi2->UncheckedAt(kLocalChi2PerDEMean+ia))->GetXaxis()->SetBinLabel(i, label);
+ }
+
+ }
+
+ // compute averaged local chi2 per chamber (X+Y)
+ TH2F* h2 = (TH2F*)fResiduals->UncheckedAt(kLocalChi2PerCh+2);
+ g = new TGraphErrors(AliMUONConstants::NTrackingCh());
+ g->SetName("gLocalChi2PerChMean");
+ g->SetTitle("local chi2 per Ch: mean;chamber ID;<local #chi^{2}>");
+ g->SetMarkerStyle(kFullDotLarge);
+ fLocalChi2->AddAtAndExpand(g, kLocalChi2PerChMean+2);
+ for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
+ TH1D* tmp = h2->ProjectionY("tmp",i+1,i+1,"e");
+ g->SetPoint(i, i+1, tmp->GetMean());
+ g->SetPointError(i, 0., tmp->GetMeanError());
+ delete tmp;
+ }
+
+ // compute averaged local chi2 per DE (X+Y)
+ h2 = (TH2F*)fResiduals->UncheckedAt(kLocalChi2PerDE+2);
+ g = new TGraphErrors(fNDE);
+ g->SetName("gLocalChi2PerDEMean");
+ g->SetTitle("local chi2 per DE: mean;DE ID;<local #chi^{2}>");
+ g->SetMarkerStyle(kFullDotLarge);
+ fLocalChi2->AddAtAndExpand(g, kLocalChi2PerDEMean+2);
+ for (Int_t i = 0; i < fNDE; i++) {
+ TH1D* tmp = h2->ProjectionY("tmp",i+1,i+1,"e");
+ g->SetPoint(i, i+1, tmp->GetMean());
+ g->SetPointError(i, 0., tmp->GetMeanError());
+ delete tmp;
+ }
+
+ // set graph labels
+ g->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
+ for (Int_t i = 1; i <= fNDE; i++) {
+ const char* label = h2->GetXaxis()->GetBinLabel(i);
+ g->GetXaxis()->SetBinLabel(i, label);
+ }
+
+ // display
+ fCanvases = new TObjArray(1000);
+ fCanvases->SetOwner();
+
+ TLegend *lResPerChMean = new TLegend(0.75,0.85,0.99,0.99);
+ TLegend *lResPerChSigma1 = new TLegend(0.75,0.85,0.99,0.99);
+ TLegend *lResPerChSigma2 = new TLegend(0.75,0.85,0.99,0.99);
+ TLegend *lResPerChSigma3 = new TLegend(0.75,0.85,0.99,0.99);
+
+ TCanvas* cResPerCh = new TCanvas("cResPerCh","cResPerCh",1200,500);
+ cResPerCh->Divide(4,2);
+ for (Int_t ia = 0; ia < 2; ia++) {
+ cResPerCh->cd(1+4*ia);
+ g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChMean_ClusterOut+ia);
+ g->Draw("ap");
+ g->SetMarkerColor(2);
+ g->SetLineColor(2);
+ if (ia == 0) lResPerChMean->AddEntry(g,"cluster out","PL");
+ g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChMean_ClusterIn+ia);
+ g->Draw("p");
+ g->SetMarkerColor(4);
+ g->SetLineColor(4);
+ if (ia == 0) lResPerChMean->AddEntry(g,"cluster in","PL");
+ if (ia == 0) lResPerChMean->Draw();
+ else lResPerChMean->DrawClone();
+ cResPerCh->cd(2+4*ia);
+ g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChSigma_ClusterOut+ia);
+ g->Draw("ap");
+ g->SetMinimum(0.);
+ g->SetMarkerColor(2);
+ g->SetLineColor(2);
+ if (ia == 0) lResPerChSigma1->AddEntry(g,"cluster out","PL");
+ g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChSigma_ClusterIn+ia);
+ g->Draw("p");
+ g->SetMarkerColor(4);
+ g->SetLineColor(4);
+ if (ia == 0) lResPerChSigma1->AddEntry(g,"cluster in","PL");
+ g = (TGraphErrors*)fChamberRes->UncheckedAt(kMCSPerChMean+ia);
+ g->Draw("p");
+ g->SetMarkerColor(5);
+ g->SetLineColor(5);
+ if (ia == 0) lResPerChSigma1->AddEntry(g,"MCS","PL");
+ g = (TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerChSigma+ia);
+ g->Draw("p");
+ g->SetMarkerColor(3);
+ g->SetLineColor(3);
+ if (ia == 0) lResPerChSigma1->AddEntry(g,"combined 1","PL");
+ if (ia == 0) lResPerChSigma1->Draw();
+ else lResPerChSigma1->DrawClone();
+ cResPerCh->cd(3+4*ia);
+ g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChDispersion_ClusterOut+ia);
+ g->Draw("ap");
+ g->SetMinimum(0.);
+ g->SetMarkerColor(2);
+ g->SetLineColor(2);
+ if (ia == 0) lResPerChSigma2->AddEntry(g,"cluster out","PL");
+ g = (TGraphErrors*)fChamberRes->UncheckedAt(kMCSPerChMean+ia);
+ g->Draw("p");
+ if (ia == 0) lResPerChSigma2->AddEntry(g,"MCS","PL");
+ g = (TGraphErrors*)fChamberRes->UncheckedAt(kTrackResPerChMean+ia);
+ g->Draw("p");
+ g->SetMarkerColor(4);
+ g->SetLineColor(4);
+ if (ia == 0) lResPerChSigma2->AddEntry(g,"track res.","PL");
+ g = (TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerCh+ia);
+ g->Draw("p");
+ if (ia == 0) lResPerChSigma2->AddEntry(g,"combined 2","PL");
+ if (ia == 0) lResPerChSigma2->Draw();
+ else lResPerChSigma2->DrawClone();
+ cResPerCh->cd(4+4*ia);
+ g = (TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerChSigma+ia);
+ g->Draw("ap");
+ g->SetMinimum(0.);
+ if (ia == 0) lResPerChSigma3->AddEntry(g,"combined 1","PL");
+ g = (TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerCh+ia);
+ g->Draw("p");
+ if (ia == 0) lResPerChSigma3->AddEntry(g,"combined 2","PL");
+ if (ia == 0) lResPerChSigma3->Draw();
+ else lResPerChSigma3->DrawClone();
+ }
+ fCanvases->AddAtAndExpand(cResPerCh, kResPerCh);
+
+ TCanvas* cResPerHalfCh = new TCanvas("cResPerHalfCh","cResPerHalfCh",1200,500);
+ cResPerHalfCh->Divide(2,2);
+ for (Int_t ia = 0; ia < 2; ia++) {
+ cResPerHalfCh->cd(1+2*ia);
+ g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMean_ClusterOut+ia);
+ g->Draw("ap");
+ g->SetMarkerColor(2);
+ g->SetLineColor(2);
+ g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMean_ClusterIn+ia);
+ g->Draw("p");
+ g->SetMarkerColor(4);
+ g->SetLineColor(4);
+ lResPerChMean->DrawClone();
+ cResPerHalfCh->cd(2+2*ia);
+ g = (TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerHalfChSigma+ia);
+ g->Draw("ap");
+ g->SetMinimum(0.);
+ g->SetMarkerColor(3);
+ g->SetLineColor(3);
+ g = (TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerHalfCh+ia);
+ g->Draw("p");
+ lResPerChSigma3->DrawClone();
+ }
+ fCanvases->AddAtAndExpand(cResPerHalfCh, kResPerHalfCh);
+
+ TCanvas* cResPerDE = new TCanvas("cResPerDE","cResPerDE",1200,800);
+ cResPerDE->Divide(1,4);
+ for (Int_t ia = 0; ia < 2; ia++) {
+ cResPerDE->cd(1+ia);
+ g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMean_ClusterOut+ia);
+ g->Draw("ap");
+ g->SetMarkerColor(2);
+ g->SetLineColor(2);
+ g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMean_ClusterIn+ia);
+ g->Draw("p");
+ g->SetMarkerColor(4);
+ g->SetLineColor(4);
+ lResPerChMean->DrawClone();
+ cResPerDE->cd(3+ia);
+ g = (TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerDESigma+ia);
+ g->Draw("ap");
+ g->SetMinimum(0.);
+ g->SetMarkerColor(3);
+ g->SetLineColor(3);
+ g = (TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerDE+ia);
+ g->Draw("p");
+ lResPerChSigma3->DrawClone();
+ }
+ fCanvases->AddAtAndExpand(cResPerDE, kResPerDE);
+
+ TCanvas* cResPerChVsP = new TCanvas("cResPerChVsP","cResPerChVsP");
+ cResPerChVsP->Divide(1,2);
+ for (Int_t ia = 0; ia < 2; ia++) {
+ cResPerChVsP->cd(1+ia);
+ mg = (TMultiGraph*)fChamberRes->UncheckedAt(kCombinedResidualSigmaVsP+ia);
+ mg->Draw("ap");
+ }
+ fCanvases->AddAtAndExpand(cResPerChVsP, kResPerChVsP);
+
+ // print results
+ printf("\nchamber resolution:\n");
+ printf(" - non-bending:");
+ for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %5.3f":", %5.3f",newClusterRes[0][i]);
+ printf("\n - bending:");
+ for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %6.4f":", %6.4f",newClusterRes[1][i]);
+ printf("\n\n");
+
+ // Post final data.
+ PostData(3, fLocalChi2);
+ PostData(4, fChamberRes);
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskMuonResolution::ModifyClusters(AliMUONTrack& track)
+{
+ /// Reset the clusters resolution from the ones given to the task and change
+ /// the cluster position according to the new alignment parameters if required
+
+ Double_t gX,gY,gZ,lX,lY,lZ;
+
+ // loop over clusters
+ Int_t nClusters = track.GetNClusters();
+ for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
+
+ AliMUONVCluster* cl = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->UncheckedAt(iCluster))->GetClusterPtr();
+
+ // change their resolution
+ cl->SetErrXY(fClusterResNB[cl->GetChamberId()], fClusterResB[cl->GetChamberId()]);
+
+ // change their position
+ if (fReAlign) {
+ gX = cl->GetX();
+ gY = cl->GetY();
+ gZ = cl->GetZ();
+ fOldGeoTransformer->Global2Local(cl->GetDetElemId(),gX,gY,gZ,lX,lY,lZ);
+ fNewGeoTransformer->Local2Global(cl->GetDetElemId(),lX,lY,lZ,gX,gY,gZ);
+ cl->SetXYZ(gX,gY,gZ);
+ }
+
+ }
+
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskMuonResolution::Zoom(TH1* h, Double_t fractionCut)
+{
+ /// Reduce the range of the histogram by removing a given fration of the statistic at each edge
+ ZoomLeft(h, fractionCut);
+ ZoomRight(h, fractionCut);
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskMuonResolution::ZoomLeft(TH1* h, Double_t fractionCut)
+{
+ /// Reduce the range of the histogram by removing a given fration of the statistic on the left side
+ Int_t maxEventsCut = (Int_t) (fractionCut * h->GetEntries());
+ Int_t nBins = h->GetNbinsX();
+
+ // set low edge
+ Int_t minBin;
+ Int_t eventsCut = 0;
+ for (minBin = 1; minBin <= nBins; minBin++) {
+ eventsCut += (Int_t) h->GetBinContent(minBin);
+ if (eventsCut > maxEventsCut) break;
+ }
+
+ // set new axis range
+ h->GetXaxis()->SetRange(--minBin, h->GetXaxis()->GetLast());
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskMuonResolution::ZoomRight(TH1* h, Double_t fractionCut)
+{
+ /// Reduce the range of the histogram by removing a given fration of the statistic on the right side
+ Int_t maxEventsCut = (Int_t) (fractionCut * h->GetEntries());
+ Int_t nBins = h->GetNbinsX();
+
+ // set high edge
+ Int_t maxBin;
+ Int_t eventsCut = 0;
+ for (maxBin = nBins; maxBin >= 1; maxBin--) {
+ eventsCut += (Int_t) h->GetBinContent(maxBin);
+ if (eventsCut > maxEventsCut) break;
+ }
+
+ // set new axis range
+ h->GetXaxis()->SetRange(h->GetXaxis()->GetFirst(), ++maxBin);
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskMuonResolution::GetMean(TH1* h, Double_t& mean, Double_t& meanErr, TGraphErrors* g, Int_t i, Double_t x, Bool_t zoom)
+{
+ /// Fill graph with the mean value of the histogram and the corresponding error (zooming if required)
+ Int_t firstBin = h->GetXaxis()->GetFirst();
+ Int_t lastBin = h->GetXaxis()->GetLast();
+ if (zoom) Zoom(h);
+ mean = (h->GetEntries() > fgkMinEntries) ? h->GetMean() : 0.;
+ meanErr = (h->GetEntries() > fgkMinEntries) ? h->GetMeanError() : 0.;
+ if (g) {
+ g->SetPoint(i, x, mean);
+ g->SetPointError(i, 0., meanErr);
+ }
+ if (zoom) h->GetXaxis()->SetRange(firstBin,lastBin);
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskMuonResolution::GetRMS(TH1* h, Double_t& rms, Double_t& rmsErr, TGraphErrors* g, Int_t i, Double_t x, Bool_t zoom)
+{
+ /// Return the RMS of the histogram and the corresponding error (zooming if required) and fill graph if !=0x0
+ Int_t firstBin = h->GetXaxis()->GetFirst();
+ Int_t lastBin = h->GetXaxis()->GetLast();
+ if (zoom) Zoom(h);
+ rms = (h->GetEntries() > fgkMinEntries) ? h->GetRMS() : 0.;
+ rmsErr = (h->GetEntries() > fgkMinEntries) ? h->GetRMSError() : 0.;
+ if (g) {
+ g->SetPoint(i, x, rms);
+ g->SetPointError(i, 0., rmsErr);
+ }
+ if (zoom) h->GetXaxis()->SetRange(firstBin,lastBin);
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskMuonResolution::FillSigmaClusterVsP(TH2* hIn, TH2* hOut, TGraphErrors* g, Bool_t zoom)
+{
+ /// Fill graph with cluster resolution from combined residuals with cluster in/out (zooming if required)
+ Double_t sigmaIn, sigmaInErr, sigmaOut, sigmaOutErr, clusterRes, clusterResErr;
+ for (Int_t j = 1; j <= hIn->GetNbinsX(); j++) {
+ TH1D* tmp = hIn->ProjectionY("tmp",j,j,"e");
+ GetRMS(tmp, sigmaIn, sigmaInErr, 0x0, 0, 0., zoom);
+ delete tmp;
+ tmp = hOut->ProjectionY("tmp",j,j,"e");
+ GetRMS(tmp, sigmaOut, sigmaOutErr, 0x0, 0, 0., zoom);
+ delete tmp;
+ Double_t p = 0.5 * (hIn->GetBinLowEdge(j) + hIn->GetBinLowEdge(j+1));
+ Double_t pErr = p - hIn->GetBinLowEdge(j);
+ clusterRes = TMath::Sqrt(sigmaIn*sigmaOut);
+ clusterResErr = (clusterRes > 0.) ? 0.5 * TMath::Sqrt(sigmaInErr*sigmaInErr*sigmaOut*sigmaOut + sigmaIn*sigmaIn*sigmaOutErr*sigmaOutErr) / clusterRes : 0.;
+ g->SetPoint(j, p, clusterRes);
+ g->SetPointError(j, pErr, clusterResErr);
+ }
+}
+
+//__________________________________________________________________________
+void AliAnalysisTaskMuonResolution::Cov2CovP(const AliMUONTrackParam ¶m, TMatrixD &covP)
+{
+ /// change coordinate system: (X, SlopeX, Y, SlopeY, q/Pyz) -> (X, Y, pX, pY, pZ)
+ /// parameters (param) are given in the (X, SlopeX, Y, SlopeY, q/Pyz) coordinate system
+
+ // Get useful parameters
+ Double_t slopeX = param.GetNonBendingSlope();
+ Double_t slopeY = param.GetBendingSlope();
+ Double_t qOverPYZ = param.GetInverseBendingMomentum();
+ Double_t pZ = param.Pz();
+
+ // compute Jacobian to change the coordinate system from (X,SlopeX,Y,SlopeY,c/pYZ) to (X,Y,pX,pY,pZ)
+ Double_t dpZdSlopeY = - qOverPYZ * qOverPYZ * pZ * pZ * pZ * slopeY;
+ Double_t dpZdQOverPYZ = (qOverPYZ != 0.) ? - pZ / qOverPYZ : - FLT_MAX;
+ TMatrixD jacob(5,5);
+ jacob.Zero();
+ jacob(0,0) = 1.;
+ jacob(1,2) = 1.;
+ jacob(2,1) = pZ;
+ jacob(2,3) = slopeX * dpZdSlopeY;
+ jacob(2,4) = slopeX * dpZdQOverPYZ;
+ jacob(3,3) = pZ + slopeY * dpZdSlopeY;
+ jacob(3,4) = slopeY * dpZdQOverPYZ;
+ jacob(4,3) = dpZdSlopeY;
+ jacob(4,4) = dpZdQOverPYZ;
+
+ // compute covariances in new coordinate system
+ TMatrixD tmp(param.GetCovariances(),TMatrixD::kMultTranspose,jacob);
+ covP.Mult(jacob,tmp);
+}
+
--- /dev/null
+//--------------------------------------------------------------------------
+// Base macro for submitting muon Resolution analysis.
+//
+// In case it is not run with full aliroot, it needs the following libraries:
+// - libSTEERBase.so
+// - libESD.so
+// - libAOD.so
+// - libANALYSIS.so
+// - libANALYSISalice.so
+// - libGui.so
+// - libMinuit.so
+// - libProofPlayer.so
+// - libXMLParser.so
+// - libRAWDatabase.so
+// - libCDB.so
+// - libSTEER.so
+// - libMUONcore.so
+// - libMUONmapping.so
+// - libMUONcalib.so
+// - libMUONgeometry.so
+// - libMUONtrigger.so
+// - libMUONraw.so
+// - libMUONbase.so
+// - libMUONrec.so
+// - libCORRFW.so
+// - libPWG3muondep.so
+//
+// The macro reads ESDs and store outputs in standard output file (AnalysisResults.root)
+// It also needs to load magnetic field, mapping, geometry (+alignment), and recoonstruction parameters from the OCDB
+//
+// Author: Philippe Pillot - SUBATECH Nantes
+//--------------------------------------------------------------------------
+
+TString aliroot="VO_ALICE@AliRoot::v4-19-15-AN";
+
+enum {kLocal, kInteractif_xml, kInteractif_ESDList, kProof};
+
+// global declaration needed in Proof mode, certainly because of the interpretor
+TStopwatch* localTimer;
+TMultiGraph* mgClusterResXVsStep;
+TMultiGraph* mgClusterResYVsStep;
+Double_t clusterResNB[10];
+Double_t clusterResB[10];
+Double_t clusterResNBErr[10];
+Double_t clusterResBErr[10];
+
+void RunMuonResolution(TString smode = "local", TString inputFileName = "AliESDs.root", Int_t nSteps = 10,
+ Bool_t selectPhysics = kFALSE, Bool_t matchTrig = kFALSE, Double_t minMomentum = 0.,
+ Bool_t correctForSystematics = kTRUE, Int_t extrapMode = 1)
+{
+ /// Compute the cluster resolution by studying cluster-track residual, deconvoluting from track resolution
+ /// if extrapMode == 0: extrapolate from the closest cluster
+ /// if extrapMode == 1: extrapolate from the previous cluster except between stations 2-3-4
+ /// if correctForSystematics == kTRUE: the systematic shifts of the residuals is included in the resolution
+
+ // timer start...
+ localTimer = new TStopwatch;
+
+ // check parameters
+ nSteps = TMath::Max(nSteps,1);
+ if (extrapMode != 0 && extrapMode != 1) {
+ Error("RunMuonResolution","incorrect extrapolation mode!");
+ return;
+ }
+
+ // Check runing mode
+ Int_t mode = GetMode(smode, inputFileName);
+ if(mode < 0){
+ Error("RunMuonResolution","Please provide either an ESD root file a collection of ESDs or a dataset.");
+ return;
+ }
+
+ gSystem->Load("libVMC");
+ gSystem->Load("libSTEERBase");
+ gSystem->Load("libESD");
+ gSystem->Load("libAOD");
+ gSystem->Load("libANALYSIS");
+ gSystem->Load("libANALYSISalice");
+
+ // Load additional libraries
+ gSystem->Load("libProofPlayer");
+ TString extraLibs="Physics:Minuit:XMLParser:Gui:RAWDatabase:CDB:STEER:MUONcore:MUONmapping:MUONcalib:MUONgeometry:MUONtrigger:MUONraw:MUONbase:MUONrec:CORRFW:PWG3muondep";
+ TObjArray* libs = extraLibs.Tokenize(":");
+ for (Int_t i = 0; i < libs->GetEntriesFast(); i++)
+ gSystem->Load(Form("lib%s",static_cast<TObjString*>(libs->UncheckedAt(i))->GetName()));
+
+ // check for old output file to removed
+ char remove = '';
+ if (!gSystem->Exec("ls chamberResolution_step*[0-9].root")) {
+ cout<<"above files must be removed from the current directory. Delete? [y=yes, n=no] "<<flush;
+ while (remove != 'y' && remove != 'n') cin>>remove;
+ if (remove == 'y') gSystem->Exec("rm -f chamberResolution_step*[0-9].root");
+ else {
+ Error("RunMuonResolution","cannot proceed with these files there otherwise results will be mixed up!");
+ return;
+ }
+ }
+
+ // OCDB access
+ if(mode == kLocal) AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
+ else if (mode != kProof) {
+ if (!TGrid::Connect("alien://")) return;
+ if(mode == kInteractif_ESDList || !gSystem->AccessPathName("ConfigureCuts.C"))
+ AliCDBManager::Instance()->SetDefaultStorage("raw://");
+ }
+
+ // Create input object
+ TObject* inputObj = 0x0;
+ if (mode == kProof) inputObj = new TObjString(inputFileName);
+ else inputObj = CreateChain(mode, inputFileName);
+ if (!inputObj) return;
+
+ // set starting chamber resolution (if -1 they will be loaded from recoParam in the task)
+ for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
+ clusterResNB[i] = -1.;
+ clusterResB[i] = -1.;
+ clusterResNBErr[i] = 0.;
+ clusterResBErr[i] = 0.;
+ }
+
+ // output graphs
+ mgClusterResXVsStep = new TMultiGraph("mgClusterResXVsStep","cluster X-resolution versus step;step;#sigma_{X} (cm)");
+ mgClusterResYVsStep = new TMultiGraph("mgClusterResYVsStep","cluster Y-resolution versus step;step;#sigma_{Y} (cm)");
+ TGraphErrors* gClusterResXVsStep[10];
+ TGraphErrors* gClusterResYVsStep[10];
+ for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
+ gClusterResXVsStep[i] = new TGraphErrors(nSteps+1);
+ gClusterResXVsStep[i]->SetName(Form("gResX_ch%d",i+1));
+ gClusterResXVsStep[i]->SetMarkerStyle(kFullDotMedium);
+ gClusterResXVsStep[i]->SetMarkerColor(i+1+i/9);
+ mgClusterResXVsStep->Add(gClusterResXVsStep[i],"lp");
+
+ gClusterResYVsStep[i] = new TGraphErrors(nSteps+1);
+ gClusterResYVsStep[i]->SetName(Form("gResY_ch%d",i+1));
+ gClusterResYVsStep[i]->SetMarkerStyle(kFullDotMedium);
+ gClusterResYVsStep[i]->SetMarkerColor(i+1+i/9);
+ mgClusterResYVsStep->Add(gClusterResYVsStep[i],"lp");
+ }
+
+ // loop over step
+ for (Int_t iStep = 0; iStep < nSteps; iStep++) {
+ cout<<"step "<<iStep+1<<"/"<<nSteps<<endl;
+
+ // Connect to proof if needed and prepare environment
+ if (mode == kProof) {
+
+ // set general environment and close previous session
+ if (iStep == 0) {
+ gEnv->SetValue("XSec.GSI.DelegProxy","2");
+ TProof::AddEnvVar("ALIROOT_EXTRA_LIBS", extraLibs.Data());
+ } else gProof->Close("s");
+
+ // connect
+ if (gSystem->Getenv("alien_API_USER") == NULL) TProof::Open("alice-caf.cern.ch","workers=20");
+ else TProof::Open(Form("%s@alice-caf.cern.ch",gSystem->Getenv("alien_API_USER")),"workers=20");
+ if (!gProof) return;
+
+ // set environment and compile task on workers
+ gProof->EnablePackage(aliroot.Data());
+
+ // prepare OCDB access on workers
+ gProof->Exec("AliCDBManager::Instance()->SetDefaultStorage(\"raw://\")");
+
+ }
+
+ // run one step
+ run(mode, iStep, inputObj, extrapMode, correctForSystematics, minMomentum, matchTrig, selectPhysics, clusterResNB, clusterResB);
+
+ // fill graph with starting resolutions from the task at first step
+ if (iStep == 0) for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
+ gClusterResXVsStep[i]->SetPoint(0, 0, clusterResNB[i]);
+ gClusterResXVsStep[i]->SetPointError(0, 0., clusterResNBErr[i]);
+ gClusterResYVsStep[i]->SetPoint(0, 0, clusterResB[i]);
+ gClusterResYVsStep[i]->SetPointError(0, 0., clusterResBErr[i]);
+ }
+
+ // read the chamber resolution from the output file
+ if (!GetChamberResolution(iStep, clusterResNB, clusterResB, clusterResNBErr, clusterResBErr)) return;
+
+ // fill graphs with computed resolutions
+ for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
+ gClusterResXVsStep[i]->SetPoint(iStep+1, iStep+1, clusterResNB[i]);
+ gClusterResXVsStep[i]->SetPointError(iStep+1, 0., clusterResNBErr[i]);
+ gClusterResYVsStep[i]->SetPoint(iStep+1, iStep+1, clusterResB[i]);
+ gClusterResYVsStep[i]->SetPointError(iStep+1, 0., clusterResBErr[i]);
+ }
+
+ }
+
+ // copy final results in results.root file
+ gSystem->Exec(Form("cp chamberResolution_step%d.root results.root", nSteps-1));
+
+ // display convergence
+ TCanvas* convergence = new TCanvas("convergence","convergence");
+ convergence->Divide(1,2);
+ convergence->cd(1);
+ mgClusterResXVsStep->Draw("ap");
+ convergence->cd(2);
+ mgClusterResYVsStep->Draw("ap");
+
+ // save convergence plots
+ TFile* outFile = TFile::Open("results.root","UPDATE");
+ if (!outFile || !outFile->IsOpen()) return;
+ outFile->cd();
+ mgClusterResXVsStep->Write();
+ mgClusterResYVsStep->Write();
+ convergence->Write();
+ outFile->Close();
+
+ // print results
+ printf("\nchamber resolution:\n");
+ printf(" - non-bending:");
+ for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %5.3f":", %5.3f",clusterResNB[i]);
+ printf("\n - bending:");
+ for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %6.4f":", %6.4f",clusterResB[i]);
+ printf("\n\n");
+
+ // ...timer stop
+ localTimer->Stop();
+ localTimer->Print();
+
+}
+
+//______________________________________________________________________________
+void run(Int_t mode, Int_t iStep, TObject* input, Int_t extrapMode, Bool_t correctForSystematics,
+ Double_t minMomentum, Bool_t matchTrig, Bool_t selectPhysics, Double_t clusterResNB[10], Double_t clusterResB[10])
+{
+ /// launch the analysis with these parameters
+
+ // Create the analysis manager
+ AliAnalysisManager *mgr = new AliAnalysisManager("MuonResolutionAnalysis");
+
+ // ESD input handler
+ AliESDInputHandler* esdH = new AliESDInputHandler();
+ esdH->SetReadFriends(kFALSE);
+ mgr->SetInputEventHandler(esdH);
+
+ // event selection
+ if (selectPhysics) {
+ gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
+ AliPhysicsSelectionTask* physicsSelection = AddTaskPhysicsSelection();
+ if(!physicsSelection) {
+ Error("RunMuonResolution","AliPhysicsSelectionTask not created!");
+ return;
+ }
+ physicsSelection->GetPhysicsSelection()->SetUseMuonTriggers();
+ }
+
+ // Muon Resolution analysis
+ gROOT->LoadMacro("$ALICE_ROOT/PWG3/muondep/AddTaskMuonResolution.C");
+ AliAnalysisManager::SetCommonFileName(Form("chamberResolution_step%d.root", iStep));
+ AliAnalysisTaskMuonResolution* muonResolution = AddTaskMuonResolution(selectPhysics, matchTrig, minMomentum, correctForSystematics, extrapMode);
+ if(!muonResolution) {
+ Error("RunMuonResolution","AliAnalysisTaskMuonResolution not created!");
+ return;
+ }
+ muonResolution->SetStartingResolution(clusterResNB, clusterResB);
+
+ // Enable debug printouts
+ //mgr->SetDebugLevel(2);
+
+ // start local analysis
+ if (mgr->InitAnalysis()) {
+ mgr->PrintStatus();
+ if (mode == kProof) mgr->StartAnalysis("proof", Form("%s",static_cast<TObjString*>(input)->GetName()));
+ else mgr->StartAnalysis("local", static_cast<TChain*>(input));
+ }
+
+ // save the summary canvases
+ if (muonResolution->GetCanvases()) {
+ TFile* outFile = TFile::Open(Form("chamberResolution_step%d.root", iStep),"UPDATE");
+ if (outFile && outFile->IsOpen()) {
+ muonResolution->GetCanvases()->Write();
+ outFile->Close();
+ }
+ }
+
+ // return starting chamber resolution from the task
+ muonResolution->GetStartingResolution(clusterResNB, clusterResB);
+
+ // clean memory
+ delete mgr;
+ TObject::SetObjectStat(kFALSE);
+}
+
+//______________________________________________________________________________
+Bool_t GetChamberResolution(Int_t iStep, Double_t clusterResNB[10], Double_t clusterResB[10], Double_t clusterResNBErr[10], Double_t clusterResBErr[10])
+{
+ /// read the chamber resolution from the output file
+
+ TFile* outFile = TFile::Open(Form("chamberResolution_step%d.root", iStep),"READ");
+
+ if (!outFile || !outFile->IsOpen()) {
+ Error("GetChamberResolution","output file does not exist!");
+ return kFALSE;
+ }
+
+ TObjArray* summary = static_cast<TObjArray*>(outFile->FindObjectAny("ChamberRes"));
+ TGraphErrors* gCombinedResidualXPerChSigma = (summary) ? static_cast<TGraphErrors*>(summary->FindObject("gCombinedResidualXPerChSigma")) : 0x0;
+ TGraphErrors* gCombinedResidualYPerChSigma = (summary) ? static_cast<TGraphErrors*>(summary->FindObject("gCombinedResidualYPerChSigma")) : 0x0;
+
+ if (!gCombinedResidualXPerChSigma || !gCombinedResidualYPerChSigma) {
+ Error("GetChamberResolution","resolution graphs do not exist!");
+ return kFALSE;
+ }
+
+ Double_t dummy;
+ for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
+ gCombinedResidualXPerChSigma->GetPoint(i, dummy, clusterResNB[i]);
+ gCombinedResidualYPerChSigma->GetPoint(i, dummy, clusterResB[i]);
+ clusterResNBErr[i] = gCombinedResidualXPerChSigma->GetErrorY(i);
+ clusterResBErr[i] = gCombinedResidualYPerChSigma->GetErrorY(i);
+ }
+
+ outFile->Close();
+
+ return kTRUE;
+}
+
+//______________________________________________________________________________
+Int_t GetMode(TString smode, TString input)
+{
+ if (smode == "local") {
+ if ( input.EndsWith(".xml") ) return kInteractif_xml;
+ else if ( input.EndsWith(".txt") ) return kInteractif_ESDList;
+ else if ( input.EndsWith(".root") ) return kLocal;
+ } else if (smode == "proof") return kProof;
+ return -1;
+}
+
+//______________________________________________________________________________
+TChain* CreateChainFromCollection(const char *xmlfile)
+{
+ // Create a chain from the collection of tags.
+ TAlienCollection* coll = TAlienCollection::Open(xmlfile);
+ if (!coll) {
+ Error("CreateChainFromCollection", Form("Cannot create an AliEn collection from %s!", xmlfile));
+ return NULL;
+ }
+
+ TGridResult* tagResult = coll->GetGridResult("",kFALSE,kFALSE);
+ AliTagAnalysis *tagAna = new AliTagAnalysis("ESD");
+ tagAna->ChainGridTags(tagResult);
+
+ AliRunTagCuts *runCuts = new AliRunTagCuts();
+ AliLHCTagCuts *lhcCuts = new AliLHCTagCuts();
+ AliDetectorTagCuts *detCuts = new AliDetectorTagCuts();
+ AliEventTagCuts *evCuts = new AliEventTagCuts();
+
+ // Check if the cuts configuration file was provided
+ if (!gSystem->AccessPathName("ConfigureCuts.C")) {
+ gROOT->LoadMacro("ConfigureCuts.C");
+ ConfigureCuts(runCuts, lhcCuts, detCuts, evCuts);
+ }
+
+ TChain *chain = tagAna->QueryTags(runCuts, lhcCuts, detCuts, evCuts);
+ if (!chain || !chain->GetNtrees()) return NULL;
+ chain->ls();
+ return chain;
+}
+
+//______________________________________________________________________________
+TChain* CreateChainFromFile(const char *rootfile)
+{
+ // Create a chain using the root file.
+ TChain* chain = new TChain("esdTree");
+ chain->Add(rootfile);
+ if (!chain->GetNtrees()) return NULL;
+ chain->ls();
+ return chain;
+}
+
+//______________________________________________________________________________
+TChain* CreateChainFromESDList(const char *esdList)
+{
+ // Create a chain using tags from the run list.
+ TChain* chain = new TChain("esdTree");
+ ifstream inFile(esdList);
+ TString inFileName;
+ if (inFile.is_open()) {
+ while (! inFile.eof() ) {
+ inFileName.ReadLine(inFile,kFALSE);
+ if(!inFileName.EndsWith(".root")) continue;
+ chain->Add(inFileName.Data());
+ }
+ }
+ inFile.close();
+ if (!chain->GetNtrees()) return NULL;
+ chain->ls();
+ return chain;
+}
+
+//______________________________________________________________________________
+TChain* CreateChain(Int_t mode, TString input)
+{
+ printf("*******************************\n");
+ printf("*** Getting the Chain ***\n");
+ printf("*******************************\n");
+ if(mode == kInteractif_xml) return CreateChainFromCollection(input.Data());
+ else if (mode == kInteractif_ESDList) return CreateChainFromESDList(input.Data());
+ else if (mode == kLocal) return CreateChainFromFile(input.Data());
+ else return NULL;
+}
+
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff