PWG3/muondep/AliAnalysisTaskMuonResolution.cxx

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