PWG3/muondep/AliAnalysisTaskMuonResolution.cxx

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