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