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