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