PWG3/muondep/AliAnalysisTaskMuonPerformance.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2007, 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 //-----------------------------------------------------------------------------
  17 /// \class AliAnalysisTaskMuonPerformance
  18 /// Analysis task to chek the tracker and trigger reconstruction
  19 /// The output is a list of histograms.
  20 ///
  21 /// \author Diego Stocco and Philippe Pillot
  22 //-----------------------------------------------------------------------------
  23
  24 // ROOT includes
  25 #include "TROOT.h"
  26 #include "TH1.h"
  27 #include "TH2.h"
  28 #include "TH3.h"
  29 #include "TF1.h"
  30 #include "TLegend.h"
  31 #include "TGraphErrors.h"
  32 #include "TGraphAsymmErrors.h"
  33 #include "TAxis.h"
  34 #include "TObjArray.h"
  35 #include "TCanvas.h"
  36 #include "TMath.h"
  37 #include "TMCProcess.h"
  38 #include "TGeoManager.h"
  39
  40 // STEER includes
  41 #include "AliMCEvent.h"
  42 #include "AliMCParticle.h"
  43 #include "AliMCEventHandler.h"
  44 #include "AliESDEvent.h"
  45 #include "AliESDMuonTrack.h"
  46 #include "AliCDBManager.h"
  47 #include "AliGeomManager.h"
  48
  49 // ANALYSIS includes
  50 #include "AliAnalysisManager.h"
  51 #include "AliAnalysisDataSlot.h"
  52 #include "AliAnalysisDataContainer.h"
  53
  54 // CORRFW includes
  55 #include "AliCFContainer.h"
  56 #include "AliCFGridSparse.h"
  57 #include "AliCFEffGrid.h"
  58
  59 // MUON includes
  60 #include "AliMUONConstants.h"
  61 #include "AliMUONTrack.h"
  62 #include "AliMUONTriggerTrack.h"
  63 #include "AliMUONLocalTrigger.h"
  64 #include "AliMUONVTrackStore.h"
  65 #include "AliMUONVTriggerTrackStore.h"
  66 #include "AliMUONESDInterface.h"
  67 #include "AliMUONRecoParam.h"
  68 #include "AliMUONCDB.h"
  69 #include "AliMUONTrackExtrap.h"
  70 #include "AliMUONTrackParam.h"
  71 #include "AliMUONRecoCheck.h"
  72 #include "AliMUONVCluster.h"
  73 #include "AliMpDEIterator.h"
  74
  75 #include "AliAnalysisTaskMuonPerformance.h"
  76
  77
  78 /// \cond CLASSIMP
  79 ClassImp(AliAnalysisTaskMuonPerformance) // Class implementation in ROOT context
  80 /// \endcond
  81
  82
  83 //________________________________________________________________________
  84 AliAnalysisTaskMuonPerformance::AliAnalysisTaskMuonPerformance() :
  85 AliAnalysisTaskSE(),
  86 fDefaultStorage(""),
  87 fNPBins(30),
  88 fCorrectForSystematics(kFALSE),
  89 fFitResiduals(kFALSE),
  90 fRequestedStationMask(0),
  91 fRequest2ChInSameSt45(0),
  92 fSigmaCut(-1.),
  93 fSigmaCutTrig(-1.),
  94 fNDE(0),
  95 fCFContainer(0x0),
  96 fEfficiencyList(0x0),
  97 fTriggerList(0x0),
  98 fTrackerList(0x0),
  99 fPAtVtxList(0x0),
 100 fSlopeAtVtxList(0x0),
 101 fEtaAtVtxList(0x0),
 102 fPhiAtVtxList(0x0),
 103 fPAt1stClList(0x0),
 104 fSlopeAt1stClList(0x0),
 105 fDCAList(0x0),
 106 fClusterList(0x0)
 107 {
 108   //
 109   /// Default Constructor.
 110   //
 111   fPRange[0] = 0.;
 112   fPRange[1] = 0.;
 113   fClusterMaxRes[0] = 0.;
 114   fClusterMaxRes[1] = 0.;
 115   for (Int_t i = 0; i < 1100; i++) fDEIndices[i] = 0;
 116   for (Int_t i = 0; i < 200; i++) fDEIds[i] = 0;
 117 }
 118
 119
 120 //________________________________________________________________________
 121 AliAnalysisTaskMuonPerformance::AliAnalysisTaskMuonPerformance(const char *name) :
 122 AliAnalysisTaskSE(name),
 123 fDefaultStorage("raw://"),
 124 fNPBins(30),
 125 fCorrectForSystematics(kFALSE),
 126 fFitResiduals(kFALSE),
 127 fRequestedStationMask(0),
 128 fRequest2ChInSameSt45(0),
 129 fSigmaCut(-1.),
 130 fSigmaCutTrig(-1.),
 131 fNDE(0),
 132 fCFContainer(0x0),
 133 fEfficiencyList(0x0),
 134 fTriggerList(0x0),
 135 fTrackerList(0x0),
 136 fPAtVtxList(0x0),
 137 fSlopeAtVtxList(0x0),
 138 fEtaAtVtxList(0x0),
 139 fPhiAtVtxList(0x0),
 140 fPAt1stClList(0x0),
 141 fSlopeAt1stClList(0x0),
 142 fDCAList(0x0),
 143 fClusterList(0x0)
 144 {
 145   //
 146   /// Constructor.
 147   //
 148   fPRange[0] = 0.;
 149   fPRange[1] = 300.;
 150   fClusterMaxRes[0] = 0.;
 151   fClusterMaxRes[1] = 0.;
 152   for (Int_t i = 0; i < 1100; i++) fDEIndices[i] = 0;
 153   for (Int_t i = 0; i < 200; i++) fDEIds[i] = 0;
 154
 155   DefineOutput(1, AliCFContainer::Class());
 156   DefineOutput(2, TObjArray::Class());
 157   DefineOutput(3, TObjArray::Class());
 158   DefineOutput(4, TObjArray::Class());
 159   DefineOutput(5, TObjArray::Class());
 160   DefineOutput(6, TObjArray::Class());
 161   DefineOutput(7, TObjArray::Class());
 162   DefineOutput(8, TObjArray::Class());
 163   DefineOutput(9, TObjArray::Class());
 164   DefineOutput(10, TObjArray::Class());
 165   DefineOutput(11, TObjArray::Class());
 166   DefineOutput(12, TObjArray::Class());
 167 }
 168
 169
 170 //________________________________________________________________________
 171 AliAnalysisTaskMuonPerformance::~AliAnalysisTaskMuonPerformance()
 172 {
 173   //
 174   /// Destructor
 175   //
 176   if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
 177     delete fCFContainer;
 178     delete fEfficiencyList;
 179     delete fTriggerList;
 180     delete fTrackerList;
 181     delete fPAtVtxList;
 182     delete fSlopeAtVtxList;
 183     delete fEtaAtVtxList;
 184     delete fPhiAtVtxList;
 185     delete fPAt1stClList;
 186     delete fSlopeAt1stClList;
 187     delete fDCAList;
 188     delete fClusterList;
 189   }
 190 }
 191
 192
 193 //___________________________________________________________________________
 194 void AliAnalysisTaskMuonPerformance::NotifyRun()
 195 {
 196   //
 197   /// Use the event handler information to correctly fill the analysis flags:
 198   /// - check if Monte Carlo information is present
 199   //
 200
 201   // load OCDB objects only once
 202   if (fSigmaCut > 0) return;
 203
 204   // set OCDB location
 205   AliCDBManager* cdbm = AliCDBManager::Instance();
 206   cdbm->SetDefaultStorage(fDefaultStorage.Data());
 207   cdbm->SetRun(fCurrentRunNumber);
 208
 209   // load magnetic field for track extrapolation
 210   if (!AliMUONCDB::LoadField()) return;
 211
 212   // load mapping
 213   if (!AliMUONCDB::LoadMapping()) return;
 214
 215   // load geometry for track extrapolation to vertex
 216   if (!AliGeomManager::GetGeometry()) AliGeomManager::LoadGeometry();
 217   if (!AliGeomManager::GetGeometry()) return;
 218
 219   // load recoParam
 220   AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam();
 221   if (!recoParam) {
 222     fRequestedStationMask = 0;
 223     fRequest2ChInSameSt45 = kFALSE;
 224     fSigmaCut = -1.;
 225     fSigmaCutTrig = -1.;
 226     AliError("--> skip this run");
 227     return;
 228   }
 229
 230   // compute the mask of requested stations from recoParam
 231   fRequestedStationMask = 0;
 232   for (Int_t i = 0; i < 5; i++) if (recoParam->RequestStation(i)) fRequestedStationMask |= ( 1 << i );
 233
 234   // get from recoParam whether a track need 2 chambers hit in the same station (4 or 5) or not to be reconstructible
 235   fRequest2ChInSameSt45 = !recoParam->MakeMoreTrackCandidates();
 236
 237   // get sigma cut from recoParam to associate clusters with TrackRefs in case the labels are not used
 238   fSigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking();
 239
 240   // get sigma cut from recoParam to associate trigger track to triggerable track
 241   fSigmaCutTrig = recoParam->GetSigmaCutForTrigger();
 242
 243   AliMUONESDInterface::ResetTracker(recoParam);
 244
 245   for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
 246
 247     // find the highest chamber resolution and set histogram bins
 248     if (recoParam->GetDefaultNonBendingReso(i) > fClusterMaxRes[0]) fClusterMaxRes[0] = recoParam->GetDefaultNonBendingReso(i);
 249     if (recoParam->GetDefaultBendingReso(i) > fClusterMaxRes[1]) fClusterMaxRes[1] = recoParam->GetDefaultBendingReso(i);
 250
 251     // fill correspondence between DEId and indices for histo (starting from 1)
 252     AliMpDEIterator it;
 253     it.First(i);
 254     while (!it.IsDone()) {
 255       fNDE++;
 256       fDEIndices[it.CurrentDEId()] = fNDE;
 257       fDEIds[fNDE] = it.CurrentDEId();
 258       it.Next();
 259     }
 260
 261   }
 262
 263   UserCreateOutputObjects();
 264 }
 265
 266
 267 //___________________________________________________________________________
 268 void AliAnalysisTaskMuonPerformance::UserCreateOutputObjects()
 269 {
 270   //
 271   /// Create output objects
 272   //
 273
 274   // do it once the OCDB has been loaded (i.e. from NotifyRun())
 275   if (fSigmaCut < 0) return;
 276
 277   // ------ CFContainer ------
 278
 279   // define axes of particle container
 280   Int_t nPtBins = 60;
 281   Double_t ptMin = 0., ptMax = 30.;
 282   TString ptName("Pt"), ptTitle("p_{t}"), ptUnits("GeV/c");
 283
 284   Int_t nEtaBins = 25;
 285   Double_t etaMin = -4.5, etaMax = -2.;
 286   TString etaName("Eta"), etaTitle("#eta"), etaUnits("a.u.");
 287
 288   Int_t nPhiBins = 36;
 289   Double_t phiMin = 0.; Double_t phiMax = 2*TMath::Pi();
 290   TString phiName("Phi"), phiTitle("#phi"), phiUnits("rad");
 291
 292   Int_t nThetaAbsEndBins = 4;
 293   Double_t thetaAbsEndMin = -0.5, thetaAbsEndMax = 3.5;
 294   TString thetaAbsEndName("ThetaAbsEnd"), thetaAbsEndTitle("#theta_{abs}"), thetaAbsEndUnits("a.u.");
 295
 296   Int_t nChargeBins = 2;
 297   Double_t chargeMin = -2, chargeMax = 2.;
 298   TString chargeName("Charge"), chargeTitle("charge"), chargeUnits("e");
 299
 300   Int_t nHasTrackerBins = 2;
 301   Double_t hasTrackerMin = -0.5, hasTrackerMax = (Double_t)nHasTrackerBins - 0.5;
 302   TString hasTrackerName("HasTracker"), hasTrackerTitle("Has tracker"), hasTrackerUnits("");
 303
 304   Int_t nTriggerBins = kNtrigCuts;
 305   Double_t triggerMin = -0.5, triggerMax = (Double_t)nTriggerBins - 0.5;
 306   TString triggerName("MatchTrig"), triggerTitle("Trigger match"), triggerUnits("");
 307
 308   Int_t nMotherTypeBins = kNtrackSources;
 309   Double_t motherTypeMin = -0.5, motherTypeMax = (Double_t)kNtrackSources - 0.5;
 310   TString motherTypeName("MotherType"), motherTypeTitle("motherType"), motherTypeUnits("");
 311
 312   Int_t nMatchMCBins = kNMatchMC;
 313   Double_t matchMCMin = -0.5, matchMCMax = (Double_t)kNMatchMC - 0.5;
 314   TString matchMCName("MatchMC"), matchMCTitle("MatchMC"), matchMCUnits("");
 315
 316   Int_t nbins[kNvars] = {nPtBins, nEtaBins, nPhiBins, nThetaAbsEndBins, nChargeBins, nHasTrackerBins, nTriggerBins, nMotherTypeBins, nMatchMCBins};
 317   Double_t xmin[kNvars] = {ptMin, etaMin, phiMin, thetaAbsEndMin, chargeMin, hasTrackerMin, triggerMin, motherTypeMin, matchMCMin};
 318   Double_t xmax[kNvars] = {ptMax, etaMax, phiMax, thetaAbsEndMax, chargeMax, hasTrackerMax, triggerMax, motherTypeMax, matchMCMax};
 319   TString axisTitle[kNvars] = {ptTitle, etaTitle, phiTitle, thetaAbsEndTitle, chargeTitle, hasTrackerTitle, triggerTitle, motherTypeTitle, matchMCTitle};
 320   TString axisUnits[kNvars] = {ptUnits, etaUnits, phiUnits, thetaAbsEndUnits, chargeUnits, hasTrackerUnits, triggerUnits, motherTypeUnits, matchMCUnits};
 321
 322   // create particle container
 323   fCFContainer = new AliCFContainer(GetOutputSlot(1)->GetContainer()->GetName(),"container for tracks",kNsteps,kNvars,nbins);
 324
 325   // set axes title and limits
 326   for ( Int_t idim = 0; idim<kNvars; idim++) {
 327     TString histoTitle = Form("%s (%s)", axisTitle[idim].Data(), axisUnits[idim].Data());
 328     histoTitle.ReplaceAll("()","");
 329     fCFContainer->SetVarTitle(idim, histoTitle.Data());
 330     fCFContainer->SetBinLimits(idim, xmin[idim], xmax[idim]);
 331   }
 332
 333   // define histo names
 334   TString stepTitle[kNsteps] = {"reconstructed", "generated"};
 335
 336   // define axes labels if any
 337   TString trigName[kNtrigCuts];
 338   trigName[kNoMatchTrig] = "NoMatch";
 339   trigName[kAllPtTrig]   = "AllPt";
 340   trigName[kLowPtTrig]   = "LowPt";
 341   trigName[kHighPtTrig]  = "HighPt";
 342
 343   TString srcName[kNtrackSources];
 344   srcName[kCharmMu]     = "Charm";
 345   srcName[kBeautyMu]    = "Beauty";
 346   srcName[kPrimaryMu]   = "Decay";
 347   srcName[kSecondaryMu] = "Secondary";
 348   srcName[kRecoHadron]  = "Hadrons";
 349   srcName[kUnknownPart] = "Fakes";
 350
 351   TString mMCName[kNMatchMC];
 352   mMCName[kNoMatch]     = "NoMatch";
 353   mMCName[kTrackerOnly] = "TrackerOnly";
 354   mMCName[kMatchedSame] = "MatchedSame";
 355   mMCName[kMatchedDiff] = "MatchedDiff";
 356   mMCName[kTriggerOnly] = "TriggerOnly";
 357
 358   // set histo name and axis labels if any
 359   for (Int_t istep=0; istep<kNsteps; istep++) {
 360
 361     // histo name
 362     fCFContainer->SetStepTitle(istep, stepTitle[istep].Data());
 363     AliCFGridSparse* gridSparse = fCFContainer->GetGrid(istep);
 364
 365     // trigger info
 366     TAxis* triggerAxis =  gridSparse->GetAxis(kVarTrigger);
 367     for ( Int_t ibin=0; ibin<kNtrigCuts; ibin++ ) {
 368       triggerAxis->SetBinLabel(ibin+1,trigName[ibin]);
 369     }
 370
 371     // mother type
 372     TAxis* motherTypeAxis = gridSparse->GetAxis(kVarMotherType);
 373     for ( Int_t ibin=0; ibin<kNtrackSources; ibin++ ) {
 374       motherTypeAxis->SetBinLabel(ibin+1,srcName[ibin]);
 375     }
 376
 377     // MC matching flag
 378     TAxis* matchMCAxis = gridSparse->GetAxis(kVarMatchMC);
 379     for ( Int_t ibin=0; ibin<kNMatchMC; ibin++ ) {
 380       matchMCAxis->SetBinLabel(ibin+1,mMCName[ibin]);
 381     }
 382
 383   }
 384
 385   // ------ trigger histograms ------
 386
 387   fTriggerList = new TObjArray(100);
 388   fTriggerList->SetOwner();
 389
 390   TH1F* h1 = new TH1F("hResTrigX11", "Residual X11;X11_{reco} - X11_{MC} (cm)", 100, -10., 10.);
 391   fTriggerList->AddAt(h1, kResTrigX11);
 392   h1 = new TH1F("hResTrigY11", "Residual Y11;Y11_{reco} - Y11_{MC} (cm)", 100, -10., 10.);
 393   fTriggerList->AddAt(h1, kResTrigY11);
 394   h1 = new TH1F("hResTrigSlopeY", "Residual slope y;ySlope_{reco} - ySlope_{MC} (rad)", 100, -0.1, 0.1);
 395   fTriggerList->AddAt(h1, kResTrigSlopeY);
 396
 397   // ------ tracker histograms ------
 398
 399   fTrackerList = new TObjArray(100);
 400   fTrackerList->SetOwner();
 401
 402   // momentum resolution at vertex
 403   const Int_t deltaPAtVtxNBins = 250;
 404   Double_t deltaPAtVtxEdges[2];
 405   deltaPAtVtxEdges[0] = -20. - 0.05 * fPRange[1];
 406   deltaPAtVtxEdges[1] =   5. + 0.05 * fPRange[1];
 407
 408   h1 = new TH1F("hResPAtVtx"," delta P at vertex;#Delta_{p} (GeV/c)",deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
 409   fTrackerList->AddAt(h1, kResPAtVtx);
 410   TH2F *h2 = new TH2F("hResPAtVtxVsP","#Delta_{p} at vertex versus p;p (GeV/c);#Delta_{p} (GeV/c)",2*fNPBins,fPRange[0],fPRange[1],deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
 411   fTrackerList->AddAt(h2, kResPAtVtxVsP);
 412   h2 = new TH2F("hResPAtVtxVsPIn23deg","#Delta_{p} at vertex versus p for tracks between 2 and 3 degrees at absorber end;p (GeV/c);#Delta_{p} (GeV/c)",2*fNPBins,fPRange[0],fPRange[1],deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
 413   fTrackerList->AddAt(h2, kResPAtVtxVsPIn23deg);
 414   h2 = new TH2F("hResPAtVtxVsPIn310deg","#Delta_{p} at vertex versus p for tracks between 3 and 10 degrees at absorber end;p (GeV/c);#Delta_{p} (GeV/c)",2*fNPBins,fPRange[0],fPRange[1],deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
 415   fTrackerList->AddAt(h2, kResPAtVtxVsPIn310deg);
 416   h2 = new TH2F("hResPAtVtxVsPIn02degMC","#Delta_{p} at vertex versus p for tracks with MC angle below 2 degrees;p (GeV/c);#Delta_{p} (GeV/c)",2*fNPBins,fPRange[0],fPRange[1],deltaPAtVtxNBins/10,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
 417   fTrackerList->AddAt(h2, kResPAtVtxVsPIn02degMC);
 418   h2 = new TH2F("hResPAtVtxVsPosAbsEndIn02degMC","#Delta_{p} at vertex versus track position at absorber end for tracks with MC angle < 2 degrees;position (cm);#Delta_{p} (GeV/c)",1000,0.,100.,deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
 419   fTrackerList->AddAt(h2, kResPAtVtxVsPosAbsEndIn02degMC);
 420   h2 = new TH2F("hResPAtVtxVsPosAbsEndIn23degMC","#Delta_{p} at vertex versus track position at absorber end for tracks with MC angle in [2,3[ degrees;position (cm);#Delta_{p} (GeV/c)",1000,0.,100.,deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
 421   fTrackerList->AddAt(h2, kResPAtVtxVsPosAbsEndIn23degMC);
 422   h2 = new TH2F("hResPAtVtxVsPosAbsEndIn310degMC","#Delta_{p} at vertex versus track position at absorber end for tracks with MC angle in [3,10[ degrees;position (cm);#Delta_{p} (GeV/c)",1000,0.,100.,deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
 423   fTrackerList->AddAt(h2, kResPAtVtxVsPosAbsEndIn310degMC);
 424   h2 = new TH2F("hResPAtVtxVsAngleAtAbsEnd","#Delta_{p} at vertex versus track position at absorber end converted to degrees;angle (Deg);#Delta_{p} (GeV/c)",10,0.,10.,deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
 425   fTrackerList->AddAt(h2, kResPAtVtxVsAngleAtAbsEnd);
 426   h2 = new TH2F("hResPAtVtxVsMCAngle","#Delta_{p} at vertex versus MC angle;MC angle (Deg);#Delta_{p} (GeV/c)",10,0.,10.,deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
 427   fTrackerList->AddAt(h2, kResPAtVtxVsMCAngle);
 428   TH3F *h3 = new TH3F("hResPAtVtxVsAngleAtAbsEndVsP","#Delta_{p} at vertex versus track position at absorber end converted to degrees versus momentum;p (GeV/c);angle (Deg);#Delta_{p} (GeV/c)",2*fNPBins,fPRange[0],fPRange[1],100,0.,10.,deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
 429   fTrackerList->AddAt(h3, kResPAtVtxVsAngleAtAbsEndVsP);
 430
 431   // transverse momentum resolution at vertex
 432   h2 = new TH2F("hResPtAtVtxVsPt","#Delta_{p_{t}} at vertex versus p_{t};p_{t} (GeV/c);#Delta_{p_{t}} (GeV/c)",2*fNPBins,fPRange[0]/10.,fPRange[1]/10.,deltaPAtVtxNBins,deltaPAtVtxEdges[0]/10.,deltaPAtVtxEdges[1]/10.);
 433   fTrackerList->AddAt(h2, kResPtAtVtxVsPt);
 434
 435   // momentum resolution at first cluster
 436   const Int_t deltaPAtFirstClNBins = 500;
 437   Double_t deltaPAtFirstClEdges[2];
 438   deltaPAtFirstClEdges[0] = -5. - 0.05 * fPRange[1];
 439   deltaPAtFirstClEdges[1] =  5. + 0.05 * fPRange[1];
 440
 441   h1 = new TH1F("hResPAt1stCl"," delta P at first cluster;#Delta_{p} (GeV/c)",deltaPAtFirstClNBins,deltaPAtFirstClEdges[0],deltaPAtFirstClEdges[1]);
 442   fTrackerList->AddAt(h1, kResPAt1stCl);
 443   h2 = new TH2F("hResPAt1stClVsP","#Delta_{p} at first cluster versus p;p (GeV/c);#Delta_{p} (GeV/c)",2*fNPBins,fPRange[0],fPRange[1],deltaPAtFirstClNBins,deltaPAtFirstClEdges[0],deltaPAtFirstClEdges[1]);
 444   fTrackerList->AddAt(h2, kResPAt1stClVsP);
 445
 446   // transverse momentum resolution at vertex
 447   h2 = new TH2F("hResPtAt1stClVsPt","#Delta_{p_{t}} at first cluster versus p_{t};p_{t} (GeV/c);#Delta_{p_{t}} (GeV/c)",2*fNPBins,fPRange[0]/10.,fPRange[1]/10.,deltaPAtFirstClNBins,deltaPAtFirstClEdges[0]/10.,deltaPAtFirstClEdges[1]/10.);
 448   fTrackerList->AddAt(h2, kResPtAt1stClVsPt);
 449
 450   // angular resolution at vertex
 451   const Int_t deltaSlopeAtVtxNBins = 500;
 452   const Double_t deltaSlopeAtVtxEdges[2] = {-0.05, 0.05};
 453
 454   h1 = new TH1F("hResSlopeXAtVtx","#Delta_{slope_{X}} at vertex;#Delta_{slope_{X}}", deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
 455   fTrackerList->AddAt(h1, kResSlopeXAtVtx);
 456   h1 = new TH1F("hResSlopeYAtVtx","#Delta_{slope_{Y}} at vertex;#Delta_{slope_{Y}}", deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
 457   fTrackerList->AddAt(h1, kResSlopeYAtVtx);
 458   h2 = new TH2F("hResSlopeXAtVtxVsP","#Delta_{slope_{X}} at vertex versus p;p (GeV/c);#Delta_{slope_{X}}",2*fNPBins,fPRange[0],fPRange[1], deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
 459   fTrackerList->AddAt(h2, kResSlopeXAtVtxVsP);
 460   h2 = new TH2F("hResSlopeYAtVtxVsP","#Delta_{slope_{Y}} at vertex versus p;p (GeV/c);#Delta_{slope_{Y}}",2*fNPBins,fPRange[0],fPRange[1], deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
 461   fTrackerList->AddAt(h2, kResSlopeYAtVtxVsP);
 462   h2 = new TH2F("hResSlopeXAtVtxVsPosAbsEndIn02degMC","#Delta_{slope_{X}} at vertex versus track position at absorber end for tracks with MC angle < 2 degrees;position (cm);#Delta_{slope_{X}}",1000,0.,100.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
 463   fTrackerList->AddAt(h2, kResSlopeXAtVtxVsPosAbsEndIn02degMC);
 464   h2 = new TH2F("hResSlopeYAtVtxVsPosAbsEndIn02degMC","#Delta_{slope_{Y}} at vertex versus track position at absorber end for tracks with MC angle < 2 degrees;position (cm);#Delta_{slope_{Y}}",1000,0.,100.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
 465   fTrackerList->AddAt(h2, kResSlopeYAtVtxVsPosAbsEndIn02degMC);
 466   h2 = new TH2F("hResSlopeXAtVtxVsPosAbsEndIn23degMC","#Delta_{slope_{X}} at vertex versus track position at absorber end for tracks with MC angle in [2,3[ degrees;position (cm);#Delta_{slope_{X}}",1000,0.,100.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
 467   fTrackerList->AddAt(h2, kResSlopeXAtVtxVsPosAbsEndIn23degMC);
 468   h2 = new TH2F("hResSlopeYAtVtxVsPosAbsEndIn23degMC","#Delta_{slope_{Y}} at vertex versus track position at absorber end for tracks with MC angle in [2,3[ degrees;position (cm);#Delta_{slope_{Y}}",1000,0.,100.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
 469   fTrackerList->AddAt(h2, kResSlopeYAtVtxVsPosAbsEndIn23degMC);
 470   h2 = new TH2F("hResSlopeXAtVtxVsPosAbsEndIn310degMC","#Delta_{slope_{X}} at vertex versus track position at absorber end for tracks with MC angle in [3,10[ degrees;position (cm);#Delta_{slope_{X}}",1000,0.,100.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
 471   fTrackerList->AddAt(h2, kResSlopeXAtVtxVsPosAbsEndIn310degMC);
 472   h2 = new TH2F("hResSlopeYAtVtxVsPosAbsEndIn310degMC","#Delta_{slope_{Y}} at vertex versus track position at absorber end for tracks with MC angle in [3,10[ degrees;position (cm);#Delta_{slope_{Y}}",1000,0.,100.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
 473   fTrackerList->AddAt(h2, kResSlopeYAtVtxVsPosAbsEndIn310degMC);
 474   h2 = new TH2F("hResSlopeXAtVtxVsAngleAtAbsEnd","#Delta_{slope_{X}} at vertex versus track position at absorber end converted to degrees;angle (Deg);#Delta_{slope_{X}}",10,0.,10.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
 475   fTrackerList->AddAt(h2, kResSlopeXAtVtxVsAngleAtAbsEnd);
 476   h2 = new TH2F("hResSlopeYAtVtxVsAngleAtAbsEnd","#Delta_{slope_{Y}} at vertex versus track position at absorber end converted to degrees;angle (Deg);#Delta_{slope_{Y}}",10,0.,10.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
 477   fTrackerList->AddAt(h2, kResSlopeYAtVtxVsAngleAtAbsEnd);
 478   h2 = new TH2F("hResSlopeXAtVtxVsMCAngle","#Delta_{slope_{X}} at vertex versus MC angle;MC angle (Deg);#Delta_{slope_{X}}",10,0.,10.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
 479   fTrackerList->AddAt(h2, kResSlopeXAtVtxVsMCAngle);
 480   h2 = new TH2F("hResSlopeYAtVtxVsMCAngle","#Delta_{slope_{Y}} at vertex versus MC angle;MC angle (Deg);#Delta_{slope_{Y}}",10,0.,10.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
 481   fTrackerList->AddAt(h2, kResSlopeYAtVtxVsMCAngle);
 482
 483   // angular resolution at first cluster
 484   const Int_t deltaSlopeAtFirstClNBins = 500;
 485   const Double_t deltaSlopeAtFirstClEdges[2] = {-0.01, 0.01};
 486
 487   h1 = new TH1F("hResSlopeXAt1stCl","#Delta_{slope_{X}} at first cluster;#Delta_{slope_{X}}", deltaSlopeAtFirstClNBins, deltaSlopeAtFirstClEdges[0], deltaSlopeAtFirstClEdges[1]);
 488   fTrackerList->AddAt(h1, kResSlopeXAt1stCl);
 489   h1 = new TH1F("hResSlopeYAt1stCl","#Delta_{slope_{Y}} at first cluster;#Delta_{slope_{Y}}", deltaSlopeAtFirstClNBins, deltaSlopeAtFirstClEdges[0], deltaSlopeAtFirstClEdges[1]);
 490   fTrackerList->AddAt(h1, kResSlopeYAt1stCl);
 491   h2 = new TH2F("hResSlopeXAt1stClVsP","#Delta_{slope_{X}} at first cluster versus p;p (GeV/c);#Delta_{slope_{X}}",2*fNPBins,fPRange[0],fPRange[1], deltaSlopeAtFirstClNBins, deltaSlopeAtFirstClEdges[0], deltaSlopeAtFirstClEdges[1]);
 492   fTrackerList->AddAt(h2, kResSlopeXAt1stClVsP);
 493   h2 = new TH2F("hResSlopeYAt1stClVsP","#Delta_{slope_{Y}} at first cluster versus p;p (GeV/c);#Delta_{slope_{Y}}",2*fNPBins,fPRange[0],fPRange[1], deltaSlopeAtFirstClNBins, deltaSlopeAtFirstClEdges[0], deltaSlopeAtFirstClEdges[1]);
 494   fTrackerList->AddAt(h2, kResSlopeYAt1stClVsP);
 495
 496   // eta resolution at vertex
 497   const Int_t deltaEtaAtVtxNBins = 500;
 498   const Double_t deltaEtaAtVtxEdges[2] = {-0.5, 0.5};
 499
 500   h1 = new TH1F("hResEtaAtVtx","#Delta_{eta} at vertex;#Delta_{eta}", deltaEtaAtVtxNBins, deltaEtaAtVtxEdges[0], deltaEtaAtVtxEdges[1]);
 501   fTrackerList->AddAt(h1, kResEtaAtVtx);
 502   h2 = new TH2F("hResEtaAtVtxVsP","#Delta_{eta} at vertex versus p;p (GeV/c);#Delta_{eta}",2*fNPBins,fPRange[0],fPRange[1], deltaEtaAtVtxNBins, deltaEtaAtVtxEdges[0], deltaEtaAtVtxEdges[1]);
 503   fTrackerList->AddAt(h2, kResEtaAtVtxVsP);
 504   h2 = new TH2F("hResEtaAtVtxVsPosAbsEndIn02degMC","#Delta_{eta} at vertex versus track position at absorber end for tracks with MC angle < 2 degrees;position (cm);#Delta_{eta}",1000,0.,100.,deltaEtaAtVtxNBins, deltaEtaAtVtxEdges[0], deltaEtaAtVtxEdges[1]);
 505   fTrackerList->AddAt(h2, kResEtaAtVtxVsPosAbsEndIn02degMC);
 506   h2 = new TH2F("hResEtaAtVtxVsPosAbsEndIn23degMC","#Delta_{eta} at vertex versus track position at absorber end for tracks with MC angle in [2,3[ degrees;position (cm);#Delta_{eta}",1000,0.,100.,deltaEtaAtVtxNBins, deltaEtaAtVtxEdges[0], deltaEtaAtVtxEdges[1]);
 507   fTrackerList->AddAt(h2, kResEtaAtVtxVsPosAbsEndIn23degMC);
 508   h2 = new TH2F("hResEtaAtVtxVsPosAbsEndIn310degMC","#Delta_{eta} at vertex versus track position at absorber end for tracks with MC angle in [3,10[ degrees;position (cm);#Delta_{eta}",1000,0.,100.,deltaEtaAtVtxNBins, deltaEtaAtVtxEdges[0], deltaEtaAtVtxEdges[1]);
 509   fTrackerList->AddAt(h2, kResEtaAtVtxVsPosAbsEndIn310degMC);
 510   h2 = new TH2F("hResEtaAtVtxVsAngleAtAbsEnd","#Delta_{eta} at vertex versus track position at absorber end converted to degrees;angle (Deg);#Delta_{eta}",10,0.,10.,deltaEtaAtVtxNBins, deltaEtaAtVtxEdges[0], deltaEtaAtVtxEdges[1]);
 511   fTrackerList->AddAt(h2, kResEtaAtVtxVsAngleAtAbsEnd);
 512   h2 = new TH2F("hResEtaAtVtxVsMCAngle","#Delta_{eta} at vertex versus MC angle;MC angle (Deg);#Delta_{eta}",10,0.,10.,deltaEtaAtVtxNBins, deltaEtaAtVtxEdges[0], deltaEtaAtVtxEdges[1]);
 513   fTrackerList->AddAt(h2, kResEtaAtVtxVsMCAngle);
 514
 515   // phi resolution at vertex
 516   const Int_t deltaPhiAtVtxNBins = 500;
 517   const Double_t deltaPhiAtVtxEdges[2] = {-0.5, 0.5};
 518
 519   h1 = new TH1F("hResPhiAtVtx","#Delta_{phi} at vertex;#Delta_{phi}", deltaPhiAtVtxNBins, deltaPhiAtVtxEdges[0], deltaPhiAtVtxEdges[1]);
 520   fTrackerList->AddAt(h1, kResPhiAtVtx);
 521   h2 = new TH2F("hResPhiAtVtxVsP","#Delta_{phi} at vertex versus p;p (GeV/c);#Delta_{phi}",2*fNPBins,fPRange[0],fPRange[1], deltaPhiAtVtxNBins, deltaPhiAtVtxEdges[0], deltaPhiAtVtxEdges[1]);
 522   fTrackerList->AddAt(h2, kResPhiAtVtxVsP);
 523   h2 = new TH2F("hResPhiAtVtxVsPosAbsEndIn02degMC","#Delta_{phi} at vertex versus track position at absorber end for tracks with MC angle < 2 degrees;position (cm);#Delta_{phi}",1000,0.,100.,deltaPhiAtVtxNBins, deltaPhiAtVtxEdges[0], deltaPhiAtVtxEdges[1]);
 524   fTrackerList->AddAt(h2, kResPhiAtVtxVsPosAbsEndIn02degMC);
 525   h2 = new TH2F("hResPhiAtVtxVsPosAbsEndIn23degMC","#Delta_{phi} at vertex versus track position at absorber end for tracks with MC angle in [2,3[ degrees;position (cm);#Delta_{phi}",1000,0.,100.,deltaPhiAtVtxNBins, deltaPhiAtVtxEdges[0], deltaPhiAtVtxEdges[1]);
 526   fTrackerList->AddAt(h2, kResPhiAtVtxVsPosAbsEndIn23degMC);
 527   h2 = new TH2F("hResPhiAtVtxVsPosAbsEndIn310degMC","#Delta_{phi} at vertex versus track position at absorber end for tracks with MC angle in [3,10[ degrees;position (cm);#Delta_{phi}",1000,0.,100.,deltaPhiAtVtxNBins, deltaPhiAtVtxEdges[0], deltaPhiAtVtxEdges[1]);
 528   fTrackerList->AddAt(h2, kResPhiAtVtxVsPosAbsEndIn310degMC);
 529   h2 = new TH2F("hResPhiAtVtxVsAngleAtAbsEnd","#Delta_{phi} at vertex versus track position at absorber end converted to degrees;angle (Deg);#Delta_{phi}",10,0.,10.,deltaPhiAtVtxNBins, deltaPhiAtVtxEdges[0], deltaPhiAtVtxEdges[1]);
 530   fTrackerList->AddAt(h2, kResPhiAtVtxVsAngleAtAbsEnd);
 531   h2 = new TH2F("hResPhiAtVtxVsMCAngle","#Delta_{phi} at vertex versus MC angle;MC angle (Deg);#Delta_{phi}",10,0.,10.,deltaPhiAtVtxNBins, deltaPhiAtVtxEdges[0], deltaPhiAtVtxEdges[1]);
 532   fTrackerList->AddAt(h2, kResPhiAtVtxVsMCAngle);
 533
 534   // DCA resolution
 535   const Int_t deltaPDCANBins = 500;
 536   const Double_t deltaPDCAEdges[2] = {0., 1000.};
 537   const Double_t deltaPMCSAngEdges[2] = {-1.5, 1.5};
 538
 539   h1 = new TH1F("hPDCA","p #times DCA at vertex;p #times DCA (GeV #times cm)", deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
 540   fTrackerList->AddAt(h1, kPDCA);
 541   h2 = new TH2F("hPDCAVsPIn23deg","p #times DCA versus p for tracks within [2,3[ degrees at absorber end;p (GeV/c);p #times DCA (GeV #times cm)",2*fNPBins,fPRange[0],fPRange[1], deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
 542   fTrackerList->AddAt(h2, kPDCAVsPIn23deg);
 543   h2 = new TH2F("hPDCAVsPIn310deg","p #times DCA versus p for tracks within [3,10[ degrees at absorber end;p (GeV/c);p #times DCA (GeV #times cm)",2*fNPBins,fPRange[0],fPRange[1], deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
 544   fTrackerList->AddAt(h2, kPDCAVsPIn310deg);
 545   h2 = new TH2F("hPDCAVsPosAbsEndIn02degMC","p #times DCA versus track position at absorber end for tracks with MC angle < 2 degrees;position (cm);p #times DCA (GeV #times cm)",1000,0.,100.,deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
 546   fTrackerList->AddAt(h2, kPDCAVsPosAbsEndIn02degMC);
 547   h2 = new TH2F("hPDCAVsPosAbsEndIn23degMC","p #times DCA}versus track position at absorber end for tracks with MC angle in [2,3[ degrees;position (cm);p #times DCA (GeV #times cm)",1000,0.,100.,deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
 548   fTrackerList->AddAt(h2, kPDCAVsPosAbsEndIn23degMC);
 549   h2 = new TH2F("hPDCAVsPosAbsEndIn310degMC","p #times DCA versus track position at absorber end for tracks with MC angle in [3,10[ degrees;position (cm);p #times DCA (GeV #times cm)",1000,0.,100.,deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
 550   fTrackerList->AddAt(h2, kPDCAVsPosAbsEndIn310degMC);
 551   h2 = new TH2F("hPDCAVsAngleAtAbsEnd","p #times DCA versus track position at absorber end converted to degrees;angle (Deg);p #times DCA (GeV #times cm)",10,0.,10.,deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
 552   fTrackerList->AddAt(h2, kPDCAVsAngleAtAbsEnd);
 553   h2 = new TH2F("hPDCAVsMCAngle","p #times DCA versus MC angle;MC angle (Deg);p #times DCA (GeV #times cm)",10,0.,10.,deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
 554   fTrackerList->AddAt(h2, kPDCAVsMCAngle);
 555
 556   // MCS angular dispersion
 557   h2 = new TH2F("hPMCSAngVsPIn23deg","p #times #Delta#theta_{MCS} versus p for tracks within [2,3[ degrees at absorber end;p (GeV/c);p #times #Delta#theta_{MCS} (GeV)",2*fNPBins,fPRange[0],fPRange[1], deltaPDCANBins, deltaPMCSAngEdges[0], deltaPMCSAngEdges[1]);
 558   fTrackerList->AddAt(h2, kPMCSAngVsPIn23deg);
 559   h2 = new TH2F("hPMCSAngVsPIn310deg","p #times #Delta#theta_{MCS} versus p for tracks within [2,3[ degrees at absorber end;p (GeV/c);p #times #Delta#theta_{MCS} (GeV)",2*fNPBins,fPRange[0],fPRange[1], deltaPDCANBins, deltaPMCSAngEdges[0], deltaPMCSAngEdges[1]);
 560   fTrackerList->AddAt(h2, kPMCSAngVsPIn310deg);
 561
 562   // cluster resolution
 563   Int_t nCh = AliMUONConstants::NTrackingCh();
 564   const Int_t clusterResNBins = 5000;
 565   Double_t clusterResMaxX = 10.*fClusterMaxRes[0];
 566   Double_t clusterResMaxY = 10.*fClusterMaxRes[1];
 567
 568   h2 = new TH2F("hResClXVsCh", "cluster-track residual-X distribution per chamber;chamber ID;#Delta_{X} (cm)", nCh, 0.5, nCh+0.5, clusterResNBins, -clusterResMaxX, clusterResMaxX);
 569   fTrackerList->AddAt(h2, kResClXVsCh);
 570   h2 = new TH2F("hResClYVsCh", "cluster-track residual-Y distribution per chamber;chamber ID;#Delta_{Y} (cm)", nCh, 0.5, nCh+0.5, clusterResNBins, -clusterResMaxY, clusterResMaxY);
 571   fTrackerList->AddAt(h2, kResClYVsCh);
 572   h2 = new TH2F("hResClXVsDE", "cluster-track residual-X distribution per DE;DE ID;#Delta_{X} (cm)", fNDE, 0.5, fNDE+0.5, clusterResNBins, -clusterResMaxX, clusterResMaxX);
 573   for (Int_t i = 1; i <= fNDE; i++) h2->GetXaxis()->SetBinLabel(i, Form("%d",fDEIds[i]));
 574   fTrackerList->AddAt(h2, kResClXVsDE);
 575   h2 = new TH2F("hResClYVsDE", "cluster-track residual-Y distribution per DE;DE ID;#Delta_{Y} (cm)", fNDE, 0.5, fNDE+0.5, clusterResNBins, -clusterResMaxY, clusterResMaxY);
 576   for (Int_t i = 1; i <= fNDE; i++) h2->GetXaxis()->SetBinLabel(i, Form("%d",fDEIds[i]));
 577   fTrackerList->AddAt(h2, kResClYVsDE);
 578
 579   // Disable printout of AliMCEvent
 580   AliLog::SetClassDebugLevel("AliMCEvent",-1);
 581
 582   PostData(1, fCFContainer);
 583   PostData(3, fTriggerList);
 584   PostData(4, fTrackerList);
 585 }
 586
 587 //________________________________________________________________________
 588 void AliAnalysisTaskMuonPerformance::UserExec(Option_t * /*option*/)
 589 {
 590   //
 591   /// Main loop
 592   /// Called for each event
 593   //
 594
 595   // check that OCDB objects have been properly set
 596   if (fSigmaCut < 0) return;
 597
 598   // Load ESD event
 599   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 600
 601   if ( ! esd ) {
 602     AliError ("ESD event not found. Nothing done!");
 603     return;
 604   }
 605
 606   // Load MC event
 607   AliMCEventHandler* mcH = static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 608   if ( ! mcH ) {
 609     AliError ("MCH event handler not found. Nothing done!");
 610     return;
 611   }
 612
 613   // Get reference tracks
 614   AliMUONRecoCheck rc(esd,mcH);
 615   AliMUONVTriggerTrackStore* triggerTrackRefStore = rc.TriggerableTracks(-1);
 616   AliMUONVTrackStore* reconstructibleStore = rc.ReconstructibleTracks(-1, fRequestedStationMask, fRequest2ChInSameSt45);
 617
 618   Double_t containerInput[kNvars];
 619   AliMUONTrackParam *trackParam;
 620   Double_t x1,y1,z1,slopex1,slopey1,pX1,pY1,pZ1,p1,pT1,eta1,phi1;
 621   Double_t x2,y2,z2,slopex2,slopey2,pX2,pY2,pZ2,p2,pT2,eta2,phi2;
 622   Double_t dPhi;
 623   Double_t xAbs,yAbs,dAbs,aAbs,aMCS,aMC;
 624   Double_t xDCA,yDCA,dca,pU;
 625
 626   // ------ Loop over reconstructed tracks ------
 627   AliESDMuonTrack *esdTrack = 0x0;
 628   Int_t nMuTracks = esd->GetNumberOfMuonTracks();
 629   for (Int_t iMuTrack = 0; iMuTrack < nMuTracks; ++iMuTrack) {
 630
 631     esdTrack = esd->GetMuonTrack(iMuTrack);
 632
 633     AliMUONTrack* matchedTrackRef = 0x0;
 634     AliMUONTriggerTrack* matchedTrigTrackRef = 0x0;
 635     containerInput[kVarMatchMC] = static_cast<Double_t>(kNoMatch);
 636
 637     // Tracker
 638     if (esdTrack->ContainTrackerData()) {
 639
 640       // convert ESD track to MUON track (without recomputing track parameters at each clusters)
 641       AliMUONTrack muonTrack;
 642       AliMUONESDInterface::ESDToMUON(*esdTrack, muonTrack, kFALSE);
 643
 644       // try to match the reconstructed track with a simulated one
 645       Int_t nMatchClusters = 0;
 646       matchedTrackRef = rc.FindCompatibleTrack(muonTrack, *reconstructibleStore, nMatchClusters, kFALSE, fSigmaCut);
 647       if (matchedTrackRef) {
 648
 649         containerInput[kVarMatchMC] = static_cast<Double_t>(kTrackerOnly);
 650
 651         // simulated track parameters at vertex
 652         trackParam = matchedTrackRef->GetTrackParamAtVertex();
 653         x1 = trackParam->GetNonBendingCoor();
 654         y1 = trackParam->GetBendingCoor();
 655         z1 = trackParam->GetZ();
 656         slopex1 = trackParam->GetNonBendingSlope();
 657         slopey1 = trackParam->GetBendingSlope();
 658         pX1 = trackParam->Px();
 659         pY1 = trackParam->Py();
 660         pZ1 = trackParam->Pz();
 661         p1  = trackParam->P();
 662         pT1 = TMath::Sqrt(pX1*pX1 + pY1*pY1);
 663         aMC = TMath::ATan(-pT1/pZ1) * TMath::RadToDeg();
 664         eta1 = TMath::Log(TMath::Tan(0.5*TMath::ATan(-pT1/pZ1)));
 665         phi1 = TMath::Pi()+TMath::ATan2(-pY1, -pX1);
 666
 667         // reconstructed track parameters at the end of the absorber
 668         AliMUONTrackParam trackParamAtAbsEnd(*((AliMUONTrackParam*)muonTrack.GetTrackParamAtCluster()->First()));
 669         AliMUONTrackExtrap::ExtrapToZ(&trackParamAtAbsEnd, AliMUONConstants::AbsZEnd());
 670         xAbs = trackParamAtAbsEnd.GetNonBendingCoor();
 671         yAbs = trackParamAtAbsEnd.GetBendingCoor();
 672         dAbs = TMath::Sqrt(xAbs*xAbs + yAbs*yAbs);
 673         aAbs = TMath::ATan(-dAbs/AliMUONConstants::AbsZEnd()) * TMath::RadToDeg();
 674         pX2 = trackParamAtAbsEnd.Px();
 675         pY2 = trackParamAtAbsEnd.Py();
 676         pZ2 = trackParamAtAbsEnd.Pz();
 677         pT2 = TMath::Sqrt(pX2*pX2 + pY2*pY2);
 678         aMCS = TMath::ATan(-pT2/pZ2) * TMath::RadToDeg();
 679
 680         // reconstructed track parameters at vertex
 681         trackParam = muonTrack.GetTrackParamAtVertex();
 682         x2 = trackParam->GetNonBendingCoor();
 683         y2 = trackParam->GetBendingCoor();
 684         z2 = trackParam->GetZ();
 685         slopex2 = trackParam->GetNonBendingSlope();
 686         slopey2 = trackParam->GetBendingSlope();
 687         pX2 = trackParam->Px();
 688         pY2 = trackParam->Py();
 689         pZ2 = trackParam->Pz();
 690         p2  = trackParam->P();
 691         pT2 = TMath::Sqrt(pX2*pX2 + pY2*pY2);
 692         eta2 = TMath::Log(TMath::Tan(0.5*TMath::ATan(-pT2/pZ2)));
 693         phi2 = TMath::Pi()+TMath::ATan2(-pY2, -pX2);
 694
 695         // reconstructed track parameters at DCA
 696         AliMUONTrackParam trackParamAtDCA(*((AliMUONTrackParam*) muonTrack.GetTrackParamAtCluster()->First()));
 697         pU = trackParamAtDCA.P();
 698         AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(&trackParamAtDCA, z2);
 699         xDCA = trackParamAtDCA.GetNonBendingCoor();
 700         yDCA = trackParamAtDCA.GetBendingCoor();
 701         dca = TMath::Sqrt(xDCA*xDCA + yDCA*yDCA);
 702
 703         dPhi = phi2-phi1;
 704         if (dPhi < -TMath::Pi()) dPhi += 2.*TMath::Pi();
 705         else if (dPhi > TMath::Pi()) dPhi -= 2.*TMath::Pi();
 706
 707         // fill histograms
 708         static_cast<TH1*>(fTrackerList->At(kResPAtVtx))->Fill(p2-p1);
 709         static_cast<TH2*>(fTrackerList->At(kResPAtVtxVsP))->Fill(p1,p2-p1);
 710         static_cast<TH2*>(fTrackerList->At(kResPAtVtxVsAngleAtAbsEnd))->Fill(aAbs,p2-p1);
 711         static_cast<TH2*>(fTrackerList->At(kResPAtVtxVsMCAngle))->Fill(aMC,p2-p1);
 712         static_cast<TH3*>(fTrackerList->At(kResPAtVtxVsAngleAtAbsEndVsP))->Fill(p1,aAbs,p2-p1);
 713         static_cast<TH2*>(fTrackerList->At(kResPtAtVtxVsPt))->Fill(pT1,pT2-pT1);
 714
 715         static_cast<TH1*>(fTrackerList->At(kResSlopeXAtVtx))->Fill(slopex2-slopex1);
 716         static_cast<TH1*>(fTrackerList->At(kResSlopeYAtVtx))->Fill(slopey2-slopey1);
 717         static_cast<TH2*>(fTrackerList->At(kResSlopeXAtVtxVsP))->Fill(p1,slopex2-slopex1);
 718         static_cast<TH2*>(fTrackerList->At(kResSlopeYAtVtxVsP))->Fill(p1,slopey2-slopey1);
 719         static_cast<TH2*>(fTrackerList->At(kResSlopeXAtVtxVsAngleAtAbsEnd))->Fill(aAbs,slopex2-slopex1);
 720         static_cast<TH2*>(fTrackerList->At(kResSlopeYAtVtxVsAngleAtAbsEnd))->Fill(aAbs,slopey2-slopey1);
 721         static_cast<TH2*>(fTrackerList->At(kResSlopeXAtVtxVsMCAngle))->Fill(aMC,slopex2-slopex1);
 722         static_cast<TH2*>(fTrackerList->At(kResSlopeYAtVtxVsMCAngle))->Fill(aMC,slopey2-slopey1);
 723
 724         static_cast<TH1*>(fTrackerList->At(kResEtaAtVtx))->Fill(eta2-eta1);
 725         static_cast<TH2*>(fTrackerList->At(kResEtaAtVtxVsP))->Fill(p1,eta2-eta1);
 726         static_cast<TH2*>(fTrackerList->At(kResEtaAtVtxVsAngleAtAbsEnd))->Fill(aAbs,eta2-eta1);
 727         static_cast<TH2*>(fTrackerList->At(kResEtaAtVtxVsMCAngle))->Fill(aMC,eta2-eta1);
 728
 729         static_cast<TH1*>(fTrackerList->At(kResPhiAtVtx))->Fill(dPhi);
 730         static_cast<TH2*>(fTrackerList->At(kResPhiAtVtxVsP))->Fill(p1,dPhi);
 731         static_cast<TH2*>(fTrackerList->At(kResPhiAtVtxVsAngleAtAbsEnd))->Fill(aAbs,dPhi);
 732         static_cast<TH2*>(fTrackerList->At(kResPhiAtVtxVsMCAngle))->Fill(aMC,dPhi);
 733
 734         static_cast<TH1*>(fTrackerList->At(kPDCA))->Fill(0.5*(p2+pU)*dca);
 735         static_cast<TH2*>(fTrackerList->At(kPDCAVsAngleAtAbsEnd))->Fill(aAbs,0.5*(p2+pU)*dca);
 736         static_cast<TH2*>(fTrackerList->At(kPDCAVsMCAngle))->Fill(aMC,0.5*(p2+pU)*dca);
 737
 738         if (aAbs > 2. && aAbs < 3.) {
 739
 740           static_cast<TH2*>(fTrackerList->At(kResPAtVtxVsPIn23deg))->Fill(p1,p2-p1);
 741           static_cast<TH2*>(fTrackerList->At(kPDCAVsPIn23deg))->Fill(p1,0.5*(p2+pU)*dca);
 742           static_cast<TH2*>(fTrackerList->At(kPMCSAngVsPIn23deg))->Fill(p1,0.5*(p2+pU)*(aMCS-aMC)*TMath::DegToRad());
 743
 744         } else if (aAbs >= 3. && aAbs < 10.) {
 745
 746           static_cast<TH2*>(fTrackerList->At(kResPAtVtxVsPIn310deg))->Fill(p1,p2-p1);
 747           static_cast<TH2*>(fTrackerList->At(kPDCAVsPIn310deg))->Fill(p1,0.5*(p2+pU)*dca);
 748           static_cast<TH2*>(fTrackerList->At(kPMCSAngVsPIn310deg))->Fill(p1,0.5*(p2+pU)*(aMCS-aMC)*TMath::DegToRad());
 749
 750         }
 751
 752         if (aMC < 2.) {
 753
 754           static_cast<TH2*>(fTrackerList->At(kResPAtVtxVsPIn02degMC))->Fill(p1,p2-p1);
 755           static_cast<TH2*>(fTrackerList->At(kResPAtVtxVsPosAbsEndIn02degMC))->Fill(dAbs,p2-p1);
 756           static_cast<TH2*>(fTrackerList->At(kResSlopeXAtVtxVsPosAbsEndIn02degMC))->Fill(dAbs,slopex2-slopex1);
 757           static_cast<TH2*>(fTrackerList->At(kResSlopeYAtVtxVsPosAbsEndIn02degMC))->Fill(dAbs,slopey2-slopey1);
 758           static_cast<TH2*>(fTrackerList->At(kResEtaAtVtxVsPosAbsEndIn02degMC))->Fill(dAbs,eta2-eta1);
 759           static_cast<TH2*>(fTrackerList->At(kResPhiAtVtxVsPosAbsEndIn02degMC))->Fill(dAbs,dPhi);
 760           static_cast<TH2*>(fTrackerList->At(kPDCAVsPosAbsEndIn02degMC))->Fill(dAbs,0.5*(p2+pU)*dca);
 761
 762         } else if (aMC >= 2. && aMC < 3) {
 763
 764           static_cast<TH2*>(fTrackerList->At(kResPAtVtxVsPosAbsEndIn23degMC))->Fill(dAbs,p2-p1);
 765           static_cast<TH2*>(fTrackerList->At(kResSlopeXAtVtxVsPosAbsEndIn23degMC))->Fill(dAbs,slopex2-slopex1);
 766           static_cast<TH2*>(fTrackerList->At(kResSlopeYAtVtxVsPosAbsEndIn23degMC))->Fill(dAbs,slopey2-slopey1);
 767           static_cast<TH2*>(fTrackerList->At(kResEtaAtVtxVsPosAbsEndIn23degMC))->Fill(dAbs,eta2-eta1);
 768           static_cast<TH2*>(fTrackerList->At(kResPhiAtVtxVsPosAbsEndIn23degMC))->Fill(dAbs,dPhi);
 769           static_cast<TH2*>(fTrackerList->At(kPDCAVsPosAbsEndIn23degMC))->Fill(dAbs,0.5*(p2+pU)*dca);
 770
 771         } else if (aMC >= 3. && aMC < 10.) {
 772
 773           static_cast<TH2*>(fTrackerList->At(kResPAtVtxVsPosAbsEndIn310degMC))->Fill(dAbs,p2-p1);
 774           static_cast<TH2*>(fTrackerList->At(kResSlopeXAtVtxVsPosAbsEndIn310degMC))->Fill(dAbs,slopex2-slopex1);
 775           static_cast<TH2*>(fTrackerList->At(kResSlopeYAtVtxVsPosAbsEndIn310degMC))->Fill(dAbs,slopey2-slopey1);
 776           static_cast<TH2*>(fTrackerList->At(kResEtaAtVtxVsPosAbsEndIn310degMC))->Fill(dAbs,eta2-eta1);
 777           static_cast<TH2*>(fTrackerList->At(kResPhiAtVtxVsPosAbsEndIn310degMC))->Fill(dAbs,dPhi);
 778           static_cast<TH2*>(fTrackerList->At(kPDCAVsPosAbsEndIn310degMC))->Fill(dAbs,0.5*(p2+pU)*dca);
 779
 780         }
 781
 782         // simulated track parameters at vertex
 783         trackParam = (AliMUONTrackParam*) matchedTrackRef->GetTrackParamAtCluster()->First();
 784         x1 = trackParam->GetNonBendingCoor();
 785         y1 = trackParam->GetBendingCoor();
 786         z1 = trackParam->GetZ();
 787         slopex1 = trackParam->GetNonBendingSlope();
 788         slopey1 = trackParam->GetBendingSlope();
 789         pX1 = trackParam->Px();
 790         pY1 = trackParam->Py();
 791         pZ1 = trackParam->Pz();
 792         p1  = trackParam->P();
 793         pT1 = TMath::Sqrt(pX1*pX1 + pY1*pY1);
 794
 795         // reconstructed track parameters at vertex
 796         trackParam = (AliMUONTrackParam*) muonTrack.GetTrackParamAtCluster()->First();
 797         x2 = trackParam->GetNonBendingCoor();
 798         y2 = trackParam->GetBendingCoor();
 799         z2 = trackParam->GetZ();
 800         slopex2 = trackParam->GetNonBendingSlope();
 801         slopey2 = trackParam->GetBendingSlope();
 802         pX2 = trackParam->Px();
 803         pY2 = trackParam->Py();
 804         pZ2 = trackParam->Pz();
 805         p2  = trackParam->P();
 806         pT2 = TMath::Sqrt(pX2*pX2 + pY2*pY2);
 807
 808         // fill histograms
 809         static_cast<TH1*>(fTrackerList->At(kResPAt1stCl))->Fill(p2-p1);
 810         static_cast<TH2*>(fTrackerList->At(kResPAt1stClVsP))->Fill(p1,p2-p1);
 811         static_cast<TH2*>(fTrackerList->At(kResPtAt1stClVsPt))->Fill(pT1,pT2-pT1);
 812         static_cast<TH1*>(fTrackerList->At(kResSlopeXAt1stCl))->Fill(slopex2-slopex1);
 813         static_cast<TH1*>(fTrackerList->At(kResSlopeYAt1stCl))->Fill(slopey2-slopey1);
 814         static_cast<TH2*>(fTrackerList->At(kResSlopeXAt1stClVsP))->Fill(p1,slopex2-slopex1);
 815         static_cast<TH2*>(fTrackerList->At(kResSlopeYAt1stClVsP))->Fill(p1,slopey2-slopey1);
 816
 817         // clusters residuals
 818         for (Int_t iCl1 = 0; iCl1 < muonTrack.GetNClusters(); iCl1++) {
 819
 820           AliMUONVCluster* cluster1 = static_cast<AliMUONTrackParam*>(muonTrack.GetTrackParamAtCluster()->UncheckedAt(iCl1))->GetClusterPtr();
 821           Int_t chId = cluster1->GetChamberId();
 822           Int_t deId = cluster1->GetDetElemId();
 823
 824           for (Int_t iCl2 = 0; iCl2 < matchedTrackRef->GetNClusters(); iCl2++) {
 825
 826             AliMUONVCluster* cluster2 = static_cast<AliMUONTrackParam*>(matchedTrackRef->GetTrackParamAtCluster()->UncheckedAt(iCl2))->GetClusterPtr();
 827
 828             if (cluster2->GetDetElemId() == deId) {
 829
 830               // fill histograms
 831               static_cast<TH2*>(fTrackerList->At(kResClXVsCh))->Fill(chId+1, cluster1->GetX()-cluster2->GetX());
 832               static_cast<TH2*>(fTrackerList->At(kResClYVsCh))->Fill(chId+1, cluster1->GetY()-cluster2->GetY());
 833               static_cast<TH2*>(fTrackerList->At(kResClXVsDE))->Fill(fDEIndices[deId], cluster1->GetX()-cluster2->GetX());
 834               static_cast<TH2*>(fTrackerList->At(kResClYVsDE))->Fill(fDEIndices[deId], cluster1->GetY()-cluster2->GetY());
 835
 836               break;
 837
 838             }
 839
 840           }
 841
 842         }
 843
 844       } // end if (matchedTrackRef)
 845
 846     } // end if (esdTrack->ContainTrackerData())
 847
 848     // Trigger
 849     if (esdTrack->ContainTriggerData()) {
 850
 851       // Convert ESD track to trigger track
 852       AliMUONLocalTrigger locTrg;
 853       AliMUONESDInterface::ESDToMUON(*esdTrack, locTrg);
 854       AliMUONTriggerTrack trigTrack;
 855       rc.TriggerToTrack(locTrg, trigTrack);
 856
 857       // try to match the trigger track with the same MC track if any
 858       if (matchedTrackRef) {
 859         AliMUONTriggerTrack* triggerTrackRef = static_cast<AliMUONTriggerTrack*>(triggerTrackRefStore->FindObject(matchedTrackRef->GetUniqueID()));
 860         if (triggerTrackRef && trigTrack.Match(*triggerTrackRef, fSigmaCutTrig)) matchedTrigTrackRef = triggerTrackRef;
 861         if (matchedTrigTrackRef) containerInput[kVarMatchMC] = static_cast<Double_t>(kMatchedSame);
 862       }
 863
 864       // or try to match with any triggerable track
 865       if (!matchedTrigTrackRef) {
 866         matchedTrigTrackRef = rc.FindCompatibleTrack(trigTrack, *triggerTrackRefStore, fSigmaCutTrig);
 867         if (matchedTrigTrackRef) containerInput[kVarMatchMC] = (matchedTrackRef) ? static_cast<Double_t>(kMatchedDiff) : static_cast<Double_t>(kTriggerOnly);
 868       }
 869
 870       // fill histograms
 871       if (matchedTrigTrackRef) {
 872         static_cast<TH1*>(fTriggerList->At(kResTrigX11))->Fill(trigTrack.GetX11() - matchedTrigTrackRef->GetX11());
 873         static_cast<TH1*>(fTriggerList->At(kResTrigY11))->Fill(trigTrack.GetY11() - matchedTrigTrackRef->GetY11());
 874         static_cast<TH1*>(fTriggerList->At(kResTrigSlopeY))->Fill(trigTrack.GetSlopeY() - matchedTrigTrackRef->GetSlopeY());
 875       }
 876
 877     }
 878
 879     // get MC particle ID
 880     Int_t mcID = -1;
 881     if (matchedTrackRef) mcID = static_cast<Int_t>(matchedTrackRef->GetUniqueID());
 882     else if (matchedTrigTrackRef) mcID = static_cast<Int_t>(matchedTrigTrackRef->GetUniqueID());
 883
 884     // fill particle container
 885     FillContainerInfo(containerInput, esdTrack, mcID);
 886     fCFContainer->Fill(containerInput, kStepReconstructed);
 887
 888   }
 889
 890   // ------ Loop over reconstructible tracks ------
 891   AliMUONTrack* trackRef = 0x0;
 892   TIter next(reconstructibleStore->CreateIterator());
 893   while ((trackRef = static_cast<AliMUONTrack*>(next()))) {
 894
 895     // find the corresponding triggerable track if any
 896     AliMUONTriggerTrack* trigTrackRef = static_cast<AliMUONTriggerTrack*>(triggerTrackRefStore->FindObject(trackRef->GetUniqueID()));
 897
 898     // fill particle container
 899     FillContainerInfo(containerInput, 0x0, static_cast<Int_t>(trackRef->GetUniqueID()));
 900     containerInput[kVarHasTracker] = 1.;
 901     if (trigTrackRef) containerInput[kVarTrigger] = static_cast<Double_t>(kAllPtTrig);
 902     containerInput[kVarMatchMC] = static_cast<Double_t>(kNoMatch);
 903     fCFContainer->Fill(containerInput, kStepGeneratedMC);
 904
 905   }
 906
 907   // ------ Loop over triggerable ghosts ------
 908   AliMUONTriggerTrack* trigTrackRef = 0x0;
 909   TIter nextTrig(triggerTrackRefStore->CreateIterator());
 910   while ((trigTrackRef = static_cast<AliMUONTriggerTrack*>(nextTrig()))) {
 911
 912     // discard tracks also reconstructible
 913     if (reconstructibleStore->FindObject(trigTrackRef->GetUniqueID())) continue;
 914
 915     // fill particle container
 916     FillContainerInfo(containerInput, 0x0, static_cast<Int_t>(trigTrackRef->GetUniqueID()));
 917     containerInput[kVarTrigger] = static_cast<Double_t>(kAllPtTrig);
 918     containerInput[kVarMatchMC] = static_cast<Double_t>(kNoMatch);
 919     fCFContainer->Fill(containerInput, kStepGeneratedMC);
 920
 921   }
 922
 923   // Post final data
 924   PostData(1, fCFContainer);
 925   PostData(3, fTriggerList);
 926   PostData(4, fTrackerList);
 927 }
 928
 929 //________________________________________________________________________
 930 void AliAnalysisTaskMuonPerformance::Terminate(Option_t *)
 931 {
 932   //
 933   /// Draw some histogram at the end.
 934   //
 935
 936   // do it only locally
 937   //if (gROOT->IsBatch()) return;
 938
 939   // get output containers
 940   fCFContainer = dynamic_cast<AliCFContainer*>(GetOutputData(1));
 941   fTriggerList = dynamic_cast<TObjArray*>(GetOutputData(3));
 942   fTrackerList = dynamic_cast<TObjArray*>(GetOutputData(4));
 943   if (!fCFContainer || !fTriggerList || !fTrackerList) {
 944     AliWarning("Output containers not found: summary histograms are not created");
 945     return;
 946   }
 947
 948   // --- compute efficiencies ---
 949   fEfficiencyList = new TObjArray(100);
 950   fEfficiencyList->SetOwner();
 951
 952   AliCFEffGrid* efficiency = new AliCFEffGrid("eff","",*fCFContainer);
 953   efficiency->CalculateEfficiency(kStepReconstructed, kStepGeneratedMC);
 954   Double_t totalEff = 0., totalEffErr = 0.;
 955   Int_t sumEffBin = 0;
 956   TH1* auxHisto = 0x0;
 957
 958   // add histogram summarizing global efficiencies
 959   TH1D* effSummary = new TH1D("effSummary", "Efficiency summary", 1, 0., 0.);
 960   effSummary->GetYaxis()->SetTitle("Efficiency");
 961   fEfficiencyList->AddLast(effSummary);
 962
 963   // Tracker efficiency using all reconstructed tracks
 964   efficiency->GetNum()->SetRangeUser(kVarHasTracker, 1., 1.);
 965   efficiency->GetDen()->SetRangeUser(kVarHasTracker, 1., 1.);
 966   auxHisto = efficiency->Project(kVarPt);
 967   auxHisto->SetName("effVsPt_trackerTracks");
 968   fEfficiencyList->AddLast(auxHisto);
 969   auxHisto = efficiency->Project(kVarEta);
 970   auxHisto->SetName("effVsEta_trackerTracks");
 971   fEfficiencyList->AddLast(auxHisto);
 972   GetEfficiency(efficiency, totalEff, totalEffErr);
 973   sumEffBin++;
 974   effSummary->Fill("Tracker_all", totalEff);
 975   effSummary->SetBinError(sumEffBin, totalEffErr);
 976   printf("Tracker efficiency using all reconstructed tracks = %f +- %f\n", totalEff, totalEffErr);
 977
 978   // Tracker efficiency using tracks matched with reconstructible ones
 979   efficiency->GetNum()->SetRangeUser(kVarMatchMC, kTrackerOnly, kMatchedDiff);
 980   auxHisto = efficiency->Project(kVarPt);
 981   auxHisto->SetName("effVsPt_trackerTracksMatchMC");
 982   fEfficiencyList->AddLast(auxHisto);
 983   auxHisto = efficiency->Project(kVarEta);
 984   auxHisto->SetName("effVsEta_trackerTracksMatchMC");
 985   fEfficiencyList->AddLast(auxHisto);
 986   GetEfficiency(efficiency, totalEff, totalEffErr);
 987   sumEffBin++;
 988   effSummary->Fill("Tracker_MCId", totalEff);
 989   effSummary->SetBinError(sumEffBin, totalEffErr);
 990   printf("Tracker efficiency using reconstructed tracks matching MC = %f +- %f\n", totalEff, totalEffErr);
 991
 992   // Matched efficiency using all reconstructed tracks
 993   efficiency->GetNum()->SetRangeUser(kVarTrigger, kAllPtTrig, kHighPtTrig);
 994   efficiency->GetDen()->SetRangeUser(kVarTrigger, kAllPtTrig, kAllPtTrig);
 995   efficiency->GetNum()->GetAxis(kVarMatchMC)->SetRange();
 996   auxHisto = efficiency->Project(kVarPt);
 997   auxHisto->SetName("effVsPt_matchedTracks");
 998   fEfficiencyList->AddLast(auxHisto);
 999   auxHisto = efficiency->Project(kVarEta);
1000   auxHisto->SetName("effVsEta_matchedTracks");
1001   fEfficiencyList->AddLast(auxHisto);
1002   GetEfficiency(efficiency, totalEff, totalEffErr);
1003   sumEffBin++;
1004   effSummary->Fill("Matched_all", totalEff);
1005   effSummary->SetBinError(sumEffBin, totalEffErr);
1006   printf("Matched efficiency using all reconstructed tracks = %f +- %f\n", totalEff, totalEffErr);
1007
1008   // Matched efficiency using tracks matched with reconstructible & triggerable ones
1009   efficiency->GetNum()->SetRangeUser(kVarMatchMC, kMatchedSame, kMatchedSame);
1010   auxHisto = efficiency->Project(kVarPt);
1011   auxHisto->SetName("effVsPt_matchedTracksMatchMC");
1012   fEfficiencyList->AddLast(auxHisto);
1013   auxHisto = efficiency->Project(kVarEta);
1014   auxHisto->SetName("effVsEta_matchedTracksMatchMC");
1015   fEfficiencyList->AddLast(auxHisto);
1016   GetEfficiency(efficiency, totalEff, totalEffErr);
1017   sumEffBin++;
1018   effSummary->Fill("Matched_MCId", totalEff);
1019   effSummary->SetBinError(sumEffBin, totalEffErr);
1020   printf("Matched efficiency using reconstructed tracks matching MC = %f +- %f\n", totalEff, totalEffErr);
1021
1022   // Trigger efficiency using all reconstructed tracks
1023   efficiency->GetNum()->GetAxis(kVarHasTracker)->SetRange();
1024   efficiency->GetDen()->GetAxis(kVarHasTracker)->SetRange();
1025   efficiency->GetNum()->GetAxis(kVarMatchMC)->SetRange();
1026   auxHisto = efficiency->Project(kVarPt);
1027   auxHisto->SetName("effVsPt_triggerTracks");
1028   fEfficiencyList->AddLast(auxHisto);
1029   auxHisto = efficiency->Project(kVarEta);
1030   auxHisto->SetName("effVsEta_triggerTracks");
1031   fEfficiencyList->AddLast(auxHisto);
1032   GetEfficiency(efficiency, totalEff, totalEffErr);
1033   sumEffBin++;
1034   effSummary->Fill("Trigger_all", totalEff);
1035   effSummary->SetBinError(sumEffBin, totalEffErr);
1036   printf("Trigger efficiency using all reconstructed tracks = %f +- %f\n", totalEff, totalEffErr);
1037
1038   // Trigger efficiency using tracks matched with triggerable ones
1039   efficiency->GetNum()->SetRangeUser(kVarMatchMC, kMatchedSame, kTriggerOnly);
1040   auxHisto = efficiency->Project(kVarPt);
1041   auxHisto->SetName("effVsPt_triggerTracksMatchMC");
1042   fEfficiencyList->AddLast(auxHisto);
1043   auxHisto = efficiency->Project(kVarEta);
1044   auxHisto->SetName("effVsEta_triggerTracksMatchMC");
1045   fEfficiencyList->AddLast(auxHisto);
1046   GetEfficiency(efficiency, totalEff, totalEffErr);
1047   sumEffBin++;
1048   effSummary->Fill("Trigger_MCId", totalEff);
1049   effSummary->SetBinError(sumEffBin, totalEffErr);
1050   printf("Trigger efficiency using reconstructed tracks matching MC = %f +- %f\n", totalEff, totalEffErr);
1051
1052   // plot histogram summarizing global efficiencies
1053   TCanvas* cEffSummary = new TCanvas("cEffSummary","Efficiency summary",20,20,310,310);
1054   cEffSummary->SetFillColor(10); cEffSummary->SetHighLightColor(10);
1055   cEffSummary->SetLeftMargin(0.15); cEffSummary->SetBottomMargin(0.15);
1056   effSummary->DrawCopy("etext");
1057
1058   // --- plot trigger resolution ---
1059   TCanvas* cTriggerResolution = new TCanvas("cTriggerResolution","Trigger resolution",10,10,310,310);
1060   cTriggerResolution->SetFillColor(10); cTriggerResolution->SetHighLightColor(10);
1061   cTriggerResolution->SetLeftMargin(0.15); cTriggerResolution->SetBottomMargin(0.15);
1062   cTriggerResolution->Divide(2,2);
1063   cTriggerResolution->cd(1);
1064   static_cast<TH1*>(fTriggerList->UncheckedAt(kResTrigX11))->DrawCopy();
1065   cTriggerResolution->cd(2);
1066   static_cast<TH1*>(fTriggerList->UncheckedAt(kResTrigY11))->DrawCopy();
1067   cTriggerResolution->cd(3);
1068   static_cast<TH1*>(fTriggerList->UncheckedAt(kResTrigSlopeY))->DrawCopy();
1069
1070   // --- compute momentum resolution at vertex ---
1071   fPAtVtxList = new TObjArray(100);
1072   fPAtVtxList->SetOwner();
1073
1074   // define graphs
1075   TGraphAsymmErrors* gMeanResPAtVtxVsP = new TGraphAsymmErrors(fNPBins);
1076   gMeanResPAtVtxVsP->SetName("gMeanResPAtVtxVsP");
1077   gMeanResPAtVtxVsP->SetTitle("<#Delta_{p}> at vertex versus p;p (GeV/c);<#Delta_{p}> (GeV/c)");
1078   fPAtVtxList->AddAt(gMeanResPAtVtxVsP, kMeanResPAtVtxVsP);
1079   TGraphAsymmErrors* gMostProbResPAtVtxVsP = new TGraphAsymmErrors(fNPBins);
1080   gMostProbResPAtVtxVsP->SetName("gMostProbResPAtVtxVsP");
1081   gMostProbResPAtVtxVsP->SetTitle("Most probable #Delta_{p} at vertex versus p;p (GeV/c);Most prob. #Delta_{p} (GeV/c)");
1082   fPAtVtxList->AddAt(gMostProbResPAtVtxVsP, kMostProbResPAtVtxVsP);
1083   TGraphAsymmErrors* gSigmaResPAtVtxVsP = new TGraphAsymmErrors(fNPBins);
1084   gSigmaResPAtVtxVsP->SetName("gSigmaResPAtVtxVsP");
1085   gSigmaResPAtVtxVsP->SetTitle("#sigma_{p}/p at vertex versus p;p (GeV/c);#sigma_{p}/p (%)");
1086   fPAtVtxList->AddAt(gSigmaResPAtVtxVsP, kSigmaResPAtVtxVsP);
1087
1088   // fit histo and fill graphs
1089   TH2* h = static_cast<TH2*>(fTrackerList->UncheckedAt(kResPAtVtxVsP));
1090   FitLandauGausResVsP(h, "momentum residuals at vertex", gMeanResPAtVtxVsP, gMostProbResPAtVtxVsP, gSigmaResPAtVtxVsP);
1091
1092   // convert resolution into relative resolution
1093   Int_t rebinFactorX = TMath::Max(h->GetNbinsX()/fNPBins, 1);
1094   for (Int_t i = rebinFactorX; i <= h->GetNbinsX(); i+=rebinFactorX) {
1095     Double_t x,y;
1096     gSigmaResPAtVtxVsP->GetPoint(i/rebinFactorX-1, x, y);
1097     gSigmaResPAtVtxVsP->SetPoint(i/rebinFactorX-1, x, 100.*y/x);
1098     gSigmaResPAtVtxVsP->SetPointEYlow(i/rebinFactorX-1, 100.*gSigmaResPAtVtxVsP->GetErrorYlow(i/rebinFactorX-1)/x);
1099     gSigmaResPAtVtxVsP->SetPointEYhigh(i/rebinFactorX-1, 100.*gSigmaResPAtVtxVsP->GetErrorYhigh(i/rebinFactorX-1)/x);
1100   }
1101
1102   // --- compute momentum resolution at first cluster ---
1103   fPAt1stClList = new TObjArray(100);
1104   fPAt1stClList->SetOwner();
1105
1106   // define graphs
1107   TGraphAsymmErrors* gMeanResPAt1stClVsP = new TGraphAsymmErrors(fNPBins);
1108   gMeanResPAt1stClVsP->SetName("gMeanResPAt1stClVsP");
1109   gMeanResPAt1stClVsP->SetTitle("<#Delta_{p}> at first cluster versus p;p (GeV/c);<#Delta_{p}> (GeV/c)");
1110   fPAt1stClList->AddAt(gMeanResPAt1stClVsP, kMeanResPAt1stClVsP);
1111   TGraphAsymmErrors* gSigmaResPAt1stClVsP = new TGraphAsymmErrors(fNPBins);
1112   gSigmaResPAt1stClVsP->SetName("gSigmaResPAt1stClVsP");
1113   gSigmaResPAt1stClVsP->SetTitle("#sigma_{p}/p at first cluster versus p;p (GeV/c);#sigma_{p}/p (%)");
1114   fPAt1stClList->AddAt(gSigmaResPAt1stClVsP, kSigmaResPAt1stClVsP);
1115
1116   // fit histo and fill graphs
1117   h = static_cast<TH2*>(fTrackerList->UncheckedAt(kResPAt1stClVsP));
1118   FitGausResVsMom(h, 0., 1., "momentum residuals at first cluster", gMeanResPAt1stClVsP, gSigmaResPAt1stClVsP);
1119
1120   // convert resolution into relative resolution
1121   rebinFactorX = TMath::Max(h->GetNbinsX()/fNPBins, 1);
1122   for (Int_t i = rebinFactorX; i <= h->GetNbinsX(); i+=rebinFactorX) {
1123     Double_t x,y;
1124     gSigmaResPAt1stClVsP->GetPoint(i/rebinFactorX-1, x, y);
1125     gSigmaResPAt1stClVsP->SetPoint(i/rebinFactorX-1, x, 100.*y/x);
1126     gSigmaResPAt1stClVsP->SetPointEYlow(i/rebinFactorX-1, 100.*gSigmaResPAt1stClVsP->GetErrorYlow(i/rebinFactorX-1)/x);
1127     gSigmaResPAt1stClVsP->SetPointEYhigh(i/rebinFactorX-1, 100.*gSigmaResPAt1stClVsP->GetErrorYhigh(i/rebinFactorX-1)/x);
1128   }
1129
1130   // --- compute slope resolution at vertex ---
1131   fSlopeAtVtxList = new TObjArray(100);
1132   fSlopeAtVtxList->SetOwner();
1133
1134   // define graphs
1135   TGraphAsymmErrors* gMeanResSlopeXAtVtxVsP = new TGraphAsymmErrors(fNPBins);
1136   gMeanResSlopeXAtVtxVsP->SetName("gMeanResSlopeXAtVtxVsP");
1137   gMeanResSlopeXAtVtxVsP->SetTitle("<#Delta_{slope_{X}}> at vertex versus p;p (GeV/c);<#Delta_{slope_{X}}>");
1138   fSlopeAtVtxList->AddAt(gMeanResSlopeXAtVtxVsP, kMeanResSlopeXAtVtxVsP);
1139   TGraphAsymmErrors* gMeanResSlopeYAtVtxVsP = new TGraphAsymmErrors(fNPBins);
1140   gMeanResSlopeYAtVtxVsP->SetName("gMeanResSlopeYAtVtxVsP");
1141   gMeanResSlopeYAtVtxVsP->SetTitle("<#Delta_{slope_{Y}}> at vertex versus p;p (GeV/c);<#Delta_{slope_{Y}}>");
1142   fSlopeAtVtxList->AddAt(gMeanResSlopeYAtVtxVsP, kMeanResSlopeYAtVtxVsP);
1143   TGraphAsymmErrors* gSigmaResSlopeXAtVtxVsP = new TGraphAsymmErrors(fNPBins);
1144   gSigmaResSlopeXAtVtxVsP->SetName("gSigmaResSlopeXAtVtxVsP");
1145   gSigmaResSlopeXAtVtxVsP->SetTitle("#sigma_{slope_{X}} at vertex versus p;p (GeV/c);#sigma_{slope_{X}}");
1146   fSlopeAtVtxList->AddAt(gSigmaResSlopeXAtVtxVsP, kSigmaResSlopeXAtVtxVsP);
1147   TGraphAsymmErrors* gSigmaResSlopeYAtVtxVsP = new TGraphAsymmErrors(fNPBins);
1148   gSigmaResSlopeYAtVtxVsP->SetName("gSigmaResSlopeYAtVtxVsP");
1149   gSigmaResSlopeYAtVtxVsP->SetTitle("#sigma_{slope_{Y}} at vertex versus p;p (GeV/c);#sigma_{slope_{Y}}");
1150   fSlopeAtVtxList->AddAt(gSigmaResSlopeYAtVtxVsP, kSigmaResSlopeYAtVtxVsP);
1151
1152   // fit histo and fill graphs
1153   FitGausResVsMom(static_cast<TH2*>(fTrackerList->UncheckedAt(kResSlopeXAtVtxVsP)), 0., 2.e-3,
1154                   "slopeX residuals at vertex", gMeanResSlopeXAtVtxVsP, gSigmaResSlopeXAtVtxVsP);
1155   FitGausResVsMom(static_cast<TH2*>(fTrackerList->UncheckedAt(kResSlopeYAtVtxVsP)), 0., 2.e-3,
1156                   "slopeY residuals at vertex", gMeanResSlopeYAtVtxVsP, gSigmaResSlopeYAtVtxVsP);
1157
1158   // --- compute slope resolution at first cluster ---
1159   fSlopeAt1stClList = new TObjArray(100);
1160   fSlopeAt1stClList->SetOwner();
1161
1162   // define graphs
1163   TGraphAsymmErrors* gMeanResSlopeXAt1stClVsP = new TGraphAsymmErrors(fNPBins);
1164   gMeanResSlopeXAt1stClVsP->SetName("gMeanResSlopeXAt1stClVsP");
1165   gMeanResSlopeXAt1stClVsP->SetTitle("<#Delta_{slope_{X}}> at first cluster versus p;p (GeV/c);<#Delta_{slope_{X}}>");
1166   fSlopeAt1stClList->AddAt(gMeanResSlopeXAt1stClVsP, kMeanResSlopeXAt1stClVsP);
1167   TGraphAsymmErrors* gMeanResSlopeYAt1stClVsP = new TGraphAsymmErrors(fNPBins);
1168   gMeanResSlopeYAt1stClVsP->SetName("gMeanResSlopeYAt1stClVsP");
1169   gMeanResSlopeYAt1stClVsP->SetTitle("<#Delta_{slope_{Y}}> at first cluster versus p;p (GeV/c);<#Delta_{slope_{Y}}>");
1170   fSlopeAt1stClList->AddAt(gMeanResSlopeYAt1stClVsP, kMeanResSlopeYAt1stClVsP);
1171   TGraphAsymmErrors* gSigmaResSlopeXAt1stClVsP = new TGraphAsymmErrors(fNPBins);
1172   gSigmaResSlopeXAt1stClVsP->SetName("gSigmaResSlopeXAt1stClVsP");
1173   gSigmaResSlopeXAt1stClVsP->SetTitle("#sigma_{slope_{X}} at first cluster versus p;p (GeV/c);#sigma_{slope_{X}}");
1174   fSlopeAt1stClList->AddAt(gSigmaResSlopeXAt1stClVsP, kSigmaResSlopeXAt1stClVsP);
1175   TGraphAsymmErrors* gSigmaResSlopeYAt1stClVsP = new TGraphAsymmErrors(fNPBins);
1176   gSigmaResSlopeYAt1stClVsP->SetName("gSigmaResSlopeYAt1stClVsP");
1177   gSigmaResSlopeYAt1stClVsP->SetTitle("#sigma_{slope_{Y}} at first cluster versus p;p (GeV/c);#sigma_{slope_{Y}}");
1178   fSlopeAt1stClList->AddAt(gSigmaResSlopeYAt1stClVsP, kSigmaResSlopeYAt1stClVsP);
1179
1180   // fit histo and fill graphs
1181   FitGausResVsMom(static_cast<TH2*>(fTrackerList->UncheckedAt(kResSlopeXAt1stClVsP)), 0., 3.e-4,
1182                   "slopeX residuals at first cluster", gMeanResSlopeXAt1stClVsP, gSigmaResSlopeXAt1stClVsP);
1183   FitGausResVsMom(static_cast<TH2*>(fTrackerList->UncheckedAt(kResSlopeYAt1stClVsP)), 0., 3.e-4,
1184                   "slopeY residuals at first cluster", gMeanResSlopeYAt1stClVsP, gSigmaResSlopeYAt1stClVsP);
1185
1186   // --- compute eta resolution at vertex ---
1187   fEtaAtVtxList = new TObjArray(100);
1188   fEtaAtVtxList->SetOwner();
1189
1190   // define graphs
1191   TGraphAsymmErrors* gMeanResEtaAtVtxVsP = new TGraphAsymmErrors(fNPBins);
1192   gMeanResEtaAtVtxVsP->SetName("gMeanResEtaAtVtxVsP");
1193   gMeanResEtaAtVtxVsP->SetTitle("<#Delta_{eta}> at vertex versus p;p (GeV/c);<#Delta_{eta}>");
1194   fEtaAtVtxList->AddAt(gMeanResEtaAtVtxVsP, kMeanResEtaAtVtxVsP);
1195   TGraphAsymmErrors* gSigmaResEtaAtVtxVsP = new TGraphAsymmErrors(fNPBins);
1196   gSigmaResEtaAtVtxVsP->SetName("gSigmaResEtaAtVtxVsP");
1197   gSigmaResEtaAtVtxVsP->SetTitle("#sigma_{eta} at vertex versus p;p (GeV/c);#sigma_{eta}");
1198   fEtaAtVtxList->AddAt(gSigmaResEtaAtVtxVsP, kSigmaResEtaAtVtxVsP);
1199
1200   // fit histo and fill graphs
1201   FitGausResVsMom(static_cast<TH2*>(fTrackerList->UncheckedAt(kResEtaAtVtxVsP)), 0., 0.1,
1202                   "eta residuals at vertex", gMeanResEtaAtVtxVsP, gSigmaResEtaAtVtxVsP);
1203
1204   // --- compute phi resolution at vertex ---
1205   fPhiAtVtxList = new TObjArray(100);
1206   fPhiAtVtxList->SetOwner();
1207
1208   // define graphs
1209   TGraphAsymmErrors* gMeanResPhiAtVtxVsP = new TGraphAsymmErrors(fNPBins);
1210   gMeanResPhiAtVtxVsP->SetName("gMeanResPhiAtVtxVsP");
1211   gMeanResPhiAtVtxVsP->SetTitle("<#Delta_{phi}> at vertex versus p;p (GeV/c);<#Delta_{phi}>");
1212   fPhiAtVtxList->AddAt(gMeanResPhiAtVtxVsP, kMeanResPhiAtVtxVsP);
1213   TGraphAsymmErrors* gSigmaResPhiAtVtxVsP = new TGraphAsymmErrors(fNPBins);
1214   gSigmaResPhiAtVtxVsP->SetName("gSigmaResPhiAtVtxVsP");
1215   gSigmaResPhiAtVtxVsP->SetTitle("#sigma_{phi} at vertex versus p;p (GeV/c);#sigma_{phi}");
1216   fPhiAtVtxList->AddAt(gSigmaResPhiAtVtxVsP, kSigmaResPhiAtVtxVsP);
1217
1218   // fit histo and fill graphs
1219   FitGausResVsMom(static_cast<TH2*>(fTrackerList->UncheckedAt(kResPhiAtVtxVsP)), 0., 0.01,
1220                   "phi residuals at vertex", gMeanResPhiAtVtxVsP, gSigmaResPhiAtVtxVsP);
1221
1222   // --- compute DCA resolution and MCS dispersion ---
1223   fDCAList = new TObjArray(100);
1224   fDCAList->SetOwner();
1225
1226   // define graphs
1227   TGraphAsymmErrors* gMeanPDCAVsPIn23deg = new TGraphAsymmErrors(fNPBins);
1228   gMeanPDCAVsPIn23deg->SetName("gMeanPDCAVsPIn23deg");
1229   gMeanPDCAVsPIn23deg->SetTitle("<p #times DCA> versus p for tracks within [2,3[ degrees at absorber end;p (GeV/c);<p #times DCA> (GeV #times cm)");
1230   fDCAList->AddAt(gMeanPDCAVsPIn23deg, kMeanPDCAVsPIn23deg);
1231   TGraphAsymmErrors* gSigmaPDCAVsPIn23deg = new TGraphAsymmErrors(fNPBins);
1232   gSigmaPDCAVsPIn23deg->SetName("gSigmaPDCAVsPIn23deg");
1233   gSigmaPDCAVsPIn23deg->SetTitle("#sigma_{p #times DCA} versus p for tracks within [2,3[ degrees at absorber end;p (GeV/c);#sigma_{p #times DCA} (GeV #times cm)");
1234   fDCAList->AddAt(gSigmaPDCAVsPIn23deg, kSigmaPDCAVsPIn23deg);
1235   TGraphAsymmErrors* gMeanPDCAVsPIn310deg = new TGraphAsymmErrors(fNPBins);
1236   gMeanPDCAVsPIn310deg->SetName("gMeanPDCAVsPIn310deg");
1237   gMeanPDCAVsPIn310deg->SetTitle("<p #times DCA> versus p for tracks within [3,10[ degrees at absorber end;p (GeV/c);<p #times DCA> (GeV #times cm)");
1238   fDCAList->AddAt(gMeanPDCAVsPIn310deg, kMeanPDCAVsPIn310deg);
1239   TGraphAsymmErrors* gSigmaPDCAVsPIn310deg = new TGraphAsymmErrors(fNPBins);
1240   gSigmaPDCAVsPIn310deg->SetName("gSigmaPDCAVsPIn310deg");
1241   gSigmaPDCAVsPIn310deg->SetTitle("#sigma_{p #times DCA} versus p for tracks within [3,10[ degrees at absorber end;p (GeV/c);#sigma_{p #times DCA} (GeV #times cm)");
1242   fDCAList->AddAt(gSigmaPDCAVsPIn310deg, kSigmaPDCAVsPIn310deg);
1243
1244   TGraphAsymmErrors* gMeanPMCSAngVsPIn23deg = new TGraphAsymmErrors(fNPBins);
1245   gMeanPMCSAngVsPIn23deg->SetName("gMeanPMCSAngVsPIn23deg");
1246   gMeanPMCSAngVsPIn23deg->SetTitle("<p #times #Delta#theta_{MCS}> versus p for tracks within [2,3[ degrees at absorber end;p (GeV/c);<p #times #Delta#theta_{MCS}> (GeV)");
1247   fDCAList->AddAt(gMeanPMCSAngVsPIn23deg, kMeanPMCSAngVsPIn23deg);
1248   TGraphAsymmErrors* gSigmaPMCSAngVsPIn23deg = new TGraphAsymmErrors(fNPBins);
1249   gSigmaPMCSAngVsPIn23deg->SetName("gSigmaPMCSAngVsPIn23deg");
1250   gSigmaPMCSAngVsPIn23deg->SetTitle("#sigma_{p #times #Delta#theta_{MCS}} versus p for tracks within [2,3[ degrees at absorber end;p (GeV/c);#sigma_{p #times #Delta#theta_{MCS}} (GeV)");
1251   fDCAList->AddAt(gSigmaPMCSAngVsPIn23deg, kSigmaPMCSAngVsPIn23deg);
1252   TGraphAsymmErrors* gMeanPMCSAngVsPIn310deg = new TGraphAsymmErrors(fNPBins);
1253   gMeanPMCSAngVsPIn310deg->SetName("gMeanPMCSAngVsPIn310deg");
1254   gMeanPMCSAngVsPIn310deg->SetTitle("<p #times #Delta#theta_{MCS}> versus p for tracks within [3,10[ degrees at absorber end;p (GeV/c);<p #times #Delta#theta_{MCS}> (GeV)");
1255   fDCAList->AddAt(gMeanPMCSAngVsPIn310deg, kMeanPMCSAngVsPIn310deg);
1256   TGraphAsymmErrors* gSigmaPMCSAngVsPIn310deg = new TGraphAsymmErrors(fNPBins);
1257   gSigmaPMCSAngVsPIn310deg->SetName("gSigmaPMCSAngVsPIn310deg");
1258   gSigmaPMCSAngVsPIn310deg->SetTitle("#sigma_{p #times #Delta#theta_{MCS}} versus p for tracks within [3,10[ degrees at absorber end;p (GeV/c);#sigma_{p #times #Delta#theta_{MCS}} (GeV)");
1259   fDCAList->AddAt(gSigmaPMCSAngVsPIn310deg, kSigmaPMCSAngVsPIn310deg);
1260
1261   // fit histo and fill graphs
1262   FitPDCAVsMom(static_cast<TH2*>(fTrackerList->UncheckedAt(kPDCAVsPIn23deg)),
1263                "p*DCA (tracks in [2,3] deg.)", gMeanPDCAVsPIn23deg, gSigmaPDCAVsPIn23deg);
1264   FitPDCAVsMom(static_cast<TH2*>(fTrackerList->UncheckedAt(kPDCAVsPIn310deg)),
1265                "p*DCA (tracks in [3,10] deg.)", gMeanPDCAVsPIn310deg, gSigmaPDCAVsPIn310deg);
1266   FitGausResVsMom(static_cast<TH2*>(fTrackerList->UncheckedAt(kPMCSAngVsPIn23deg)), 0., 2.e-3,
1267                   "p*MCSAngle (tracks in [2,3] deg.)", gMeanPMCSAngVsPIn23deg, gSigmaPMCSAngVsPIn23deg);
1268   FitGausResVsMom(static_cast<TH2*>(fTrackerList->UncheckedAt(kPMCSAngVsPIn310deg)), 0., 2.e-3,
1269                   "p*MCSAngle (tracks in [3,10] deg.)", gMeanPMCSAngVsPIn310deg, gSigmaPMCSAngVsPIn310deg);
1270
1271   // --- compute cluster resolution ---
1272   fClusterList = new TObjArray(100);
1273   fClusterList->SetOwner();
1274
1275   // define graphs per chamber
1276   TGraphErrors* gMeanResClXVsCh = new TGraphErrors(AliMUONConstants::NTrackingCh());
1277   gMeanResClXVsCh->SetName("gMeanResClXVsCh");
1278   gMeanResClXVsCh->SetTitle("cluster-trackRef residual-X per Ch: mean;chamber ID;<#Delta_{X}> (cm)");
1279   gMeanResClXVsCh->SetMarkerStyle(kFullDotLarge);
1280   fClusterList->AddAt(gMeanResClXVsCh, kMeanResClXVsCh);
1281   TGraphErrors* gMeanResClYVsCh = new TGraphErrors(AliMUONConstants::NTrackingCh());
1282   gMeanResClYVsCh->SetName("gMeanResClYVsCh");
1283   gMeanResClYVsCh->SetTitle("cluster-trackRef residual-Y per Ch: mean;chamber ID;<#Delta_{Y}> (cm)");
1284   gMeanResClYVsCh->SetMarkerStyle(kFullDotLarge);
1285   fClusterList->AddAt(gMeanResClYVsCh, kMeanResClYVsCh);
1286   TGraphErrors* gSigmaResClXVsCh = new TGraphErrors(AliMUONConstants::NTrackingCh());
1287   gSigmaResClXVsCh->SetName("gSigmaResClXVsCh");
1288   gSigmaResClXVsCh->SetTitle("cluster-trackRef residual-X per Ch: sigma;chamber ID;#sigma_{X} (cm)");
1289   gSigmaResClXVsCh->SetMarkerStyle(kFullDotLarge);
1290   fClusterList->AddAt(gSigmaResClXVsCh, kSigmaResClXVsCh);
1291   TGraphErrors* gSigmaResClYVsCh = new TGraphErrors(AliMUONConstants::NTrackingCh());
1292   gSigmaResClYVsCh->SetName("gSigmaResClYVsCh");
1293   gSigmaResClYVsCh->SetTitle("cluster-trackRef residual-Y per Ch: sigma;chamber ID;#sigma_{Y} (cm)");
1294   gSigmaResClYVsCh->SetMarkerStyle(kFullDotLarge);
1295   fClusterList->AddAt(gSigmaResClYVsCh, kSigmaResClYVsCh);
1296
1297   // define graphs per DE
1298   TGraphErrors* gMeanResClXVsDE = new TGraphErrors(fNDE);
1299   gMeanResClXVsDE->SetName("gMeanResClXVsDE");
1300   gMeanResClXVsDE->SetTitle("cluster-trackRef residual-X per DE: mean;DE ID;<#Delta_{X}> (cm)");
1301   gMeanResClXVsDE->SetMarkerStyle(kFullDotLarge);
1302   fClusterList->AddAt(gMeanResClXVsDE, kMeanResClXVsDE);
1303   TGraphErrors* gMeanResClYVsDE = new TGraphErrors(fNDE);
1304   gMeanResClYVsDE->SetName("gMeanResClYVsDE");
1305   gMeanResClYVsDE->SetTitle("cluster-trackRef residual-Y per dE: mean;DE ID;<#Delta_{Y}> (cm)");
1306   gMeanResClYVsDE->SetMarkerStyle(kFullDotLarge);
1307   fClusterList->AddAt(gMeanResClYVsDE, kMeanResClYVsDE);
1308   TGraphErrors* gSigmaResClXVsDE = new TGraphErrors(fNDE);
1309   gSigmaResClXVsDE->SetName("gSigmaResClXVsDE");
1310   gSigmaResClXVsDE->SetTitle("cluster-trackRef residual-X per DE: sigma;DE ID;#sigma_{X} (cm)");
1311   gSigmaResClXVsDE->SetMarkerStyle(kFullDotLarge);
1312   fClusterList->AddAt(gSigmaResClXVsDE, kSigmaResClXVsDE);
1313   TGraphErrors* gSigmaResClYVsDE = new TGraphErrors(fNDE);
1314   gSigmaResClYVsDE->SetName("gSigmaResClYVsDE");
1315   gSigmaResClYVsDE->SetTitle("cluster-trackRef residual-Y per DE: sigma;DE ID;#sigma_{Y} (cm)");
1316   gSigmaResClYVsDE->SetMarkerStyle(kFullDotLarge);
1317   fClusterList->AddAt(gSigmaResClYVsDE, kSigmaResClYVsDE);
1318
1319   // fit histo and fill graphs per chamber
1320   Double_t clusterResPerCh[10][2];
1321   for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
1322     TH1D *tmp = static_cast<TH2*>(fTrackerList->UncheckedAt(kResClXVsCh))->ProjectionY("tmp",i+1,i+1,"e");
1323     FitClusterResidual(tmp, i, clusterResPerCh[i][0], gMeanResClXVsCh, gSigmaResClXVsCh);
1324     delete tmp;
1325     tmp = static_cast<TH2*>(fTrackerList->UncheckedAt(kResClYVsCh))->ProjectionY("tmp",i+1,i+1,"e");
1326     FitClusterResidual(tmp, i, clusterResPerCh[i][1], gMeanResClYVsCh, gSigmaResClYVsCh);
1327     delete tmp;
1328   }
1329
1330   // fit histo and fill graphs per DE
1331   Double_t clusterResPerDE[200][2];
1332   for (Int_t i = 0; i < fNDE; i++) {
1333     TH1D *tmp = static_cast<TH2*>(fTrackerList->UncheckedAt(kResClXVsDE))->ProjectionY("tmp",i+1,i+1,"e");
1334     FitClusterResidual(tmp, i, clusterResPerDE[i][0], gMeanResClXVsDE, gSigmaResClXVsDE);
1335     delete tmp;
1336     tmp = static_cast<TH2*>(fTrackerList->UncheckedAt(kResClYVsDE))->ProjectionY("tmp",i+1,i+1,"e");
1337     FitClusterResidual(tmp, i, clusterResPerDE[i][1], gMeanResClYVsDE, gSigmaResClYVsDE);
1338     delete tmp;
1339   }
1340
1341   // set DE graph labels
1342   TAxis* xAxis = static_cast<TH2*>(fTrackerList->UncheckedAt(kResClXVsDE))->GetXaxis();
1343   gMeanResClXVsDE->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
1344   gMeanResClYVsDE->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
1345   gSigmaResClXVsDE->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
1346   gSigmaResClYVsDE->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
1347   for (Int_t i = 1; i <= fNDE; i++) {
1348     const char* label = xAxis->GetBinLabel(i);
1349     gMeanResClXVsDE->GetXaxis()->SetBinLabel(i, label);
1350     gMeanResClYVsDE->GetXaxis()->SetBinLabel(i, label);
1351     gSigmaResClXVsDE->GetXaxis()->SetBinLabel(i, label);
1352     gSigmaResClYVsDE->GetXaxis()->SetBinLabel(i, label);
1353   }
1354
1355   // --- diplay momentum residuals at vertex ---
1356   TCanvas* cResPAtVtx = DrawVsAng("cResPAtVtx", "momentum residual at vertex in 3 angular regions",
1357                                   static_cast<TH1*>(fTrackerList->UncheckedAt(kResPAtVtx)),
1358                                   static_cast<TH2*>(fTrackerList->UncheckedAt(kResPAtVtxVsAngleAtAbsEnd)));
1359   fPAtVtxList->AddAt(cResPAtVtx, kcResPAtVtx);
1360   TCanvas* cResPAtVtxMC = DrawVsAng("cResPAtVtxMC", "momentum residual at vertex in 3 MC angular regions",
1361                                     static_cast<TH1*>(fTrackerList->UncheckedAt(kResPAtVtx)),
1362                                     static_cast<TH2*>(fTrackerList->UncheckedAt(kResPAtVtxVsMCAngle)));
1363   fPAtVtxList->AddAt(cResPAtVtxMC, kcResPAtVtxMC);
1364   TCanvas* cResPAtVtxVsPosAbsEndMC = DrawVsPos("cResPAtVtxVsPosAbsEndMC", "momentum residual at vertex versus position at absorber end in 3 MC angular regions",
1365                                                static_cast<TH2*>(fTrackerList->UncheckedAt(kResPAtVtxVsPosAbsEndIn02degMC)),
1366                                                static_cast<TH2*>(fTrackerList->UncheckedAt(kResPAtVtxVsPosAbsEndIn23degMC)),
1367                                                static_cast<TH2*>(fTrackerList->UncheckedAt(kResPAtVtxVsPosAbsEndIn310degMC)));
1368   fPAtVtxList->AddAt(cResPAtVtxVsPosAbsEndMC, kcResPAtVtxVsPosAbsEndMC);
1369   TCanvas* cResPAtVtxVsPIn23deg = DrawFitLandauGausResPVsP("cResPAtVtxVsPIn23deg", "momentum residual for tracks between 2 and 3 degrees",
1370                                                       static_cast<TH2*>(fTrackerList->UncheckedAt(kResPAtVtxVsPIn23deg)),
1371                                                       10, "momentum residuals at vertex (tracks in [2,3] deg.)");
1372   fPAtVtxList->AddAt(cResPAtVtxVsPIn23deg, kcResPAtVtxVsPIn23deg);
1373   TCanvas* cResPAtVtxVsPIn310deg = DrawFitLandauGausResPVsP("cResPAtVtxVsPIn310deg", "momentum residual for tracks between 3 and 10 degrees",
1374                                                        static_cast<TH2*>(fTrackerList->UncheckedAt(kResPAtVtxVsPIn310deg)),
1375                                                        10, "momentum residuals at vertex (tracks in [3,10] deg.)");
1376   fPAtVtxList->AddAt(cResPAtVtxVsPIn310deg, kcResPAtVtxVsPIn310deg);
1377   TCanvas* cResPAtVtxVsPIn02degMC = DrawResPVsP("cResPAtVtxVsPIn02degMC", "momentum residuals for tracks with MC angle < 2 degrees",
1378                                            static_cast<TH2*>(fTrackerList->UncheckedAt(kResPAtVtxVsPIn02degMC)), 5);
1379   fPAtVtxList->AddAt(cResPAtVtxVsPIn02degMC, kcResPAtVtxVsPIn02degMC);
1380
1381   // --- diplay slopeX residuals at vertex ---
1382   TCanvas* cResSlopeXAtVtx = DrawVsAng("cResSlopeXAtVtx", "slope_{X} residual at vertex in 3 angular regions",
1383                                        static_cast<TH1*>(fTrackerList->UncheckedAt(kResSlopeXAtVtx)),
1384                                        static_cast<TH2*>(fTrackerList->UncheckedAt(kResSlopeXAtVtxVsAngleAtAbsEnd)));
1385   fSlopeAtVtxList->AddAt(cResSlopeXAtVtx, kcResSlopeXAtVtx);
1386   TCanvas* cResSlopeYAtVtx = DrawVsAng("cResSlopeYAtVtx", "slope_{Y} residual at vertex in 3 angular regions",
1387                                        static_cast<TH1*>(fTrackerList->UncheckedAt(kResSlopeYAtVtx)),
1388                                        static_cast<TH2*>(fTrackerList->UncheckedAt(kResSlopeYAtVtxVsAngleAtAbsEnd)));
1389   fSlopeAtVtxList->AddAt(cResSlopeYAtVtx, kcResSlopeYAtVtx);
1390   TCanvas* cResSlopeXAtVtxMC = DrawVsAng("cResSlopeXAtVtxMC", "slope_{X} residual at vertex in 3 MC angular regions",
1391                                          static_cast<TH1*>(fTrackerList->UncheckedAt(kResSlopeXAtVtx)),
1392                                          static_cast<TH2*>(fTrackerList->UncheckedAt(kResSlopeXAtVtxVsMCAngle)));
1393   fSlopeAtVtxList->AddAt(cResSlopeXAtVtxMC, kcResSlopeXAtVtxMC);
1394   TCanvas* cResSlopeYAtVtxMC = DrawVsAng("cResSlopeYAtVtxMC", "slope_{Y} residual at vertex in 3 MC angular regions",
1395                                          static_cast<TH1*>(fTrackerList->UncheckedAt(kResSlopeYAtVtx)),
1396                                          static_cast<TH2*>(fTrackerList->UncheckedAt(kResSlopeYAtVtxVsMCAngle)));
1397   fSlopeAtVtxList->AddAt(cResSlopeYAtVtxMC, kcResSlopeYAtVtxMC);
1398   TCanvas* cResSlopeXAtVtxVsPosAbsEndMC = DrawVsPos("cResSlopeXAtVtxVsPosAbsEndMC", "slope_{X} residual at vertex versus position at absorber end in 3 MC angular regions",
1399                                                     static_cast<TH2*>(fTrackerList->UncheckedAt(kResSlopeXAtVtxVsPosAbsEndIn02degMC)),
1400                                                     static_cast<TH2*>(fTrackerList->UncheckedAt(kResSlopeXAtVtxVsPosAbsEndIn23degMC)),
1401                                                     static_cast<TH2*>(fTrackerList->UncheckedAt(kResSlopeXAtVtxVsPosAbsEndIn310degMC)));
1402   fSlopeAtVtxList->AddAt(cResSlopeXAtVtxVsPosAbsEndMC, kcResSlopeXAtVtxVsPosAbsEndMC);
1403   TCanvas* cResSlopeYAtVtxVsPosAbsEndMC = DrawVsPos("cResSlopeYAtVtxVsPosAbsEndMC", "slope_{Y} residual at vertex versus position at absorber end in 3 MC angular regions",
1404                                                     static_cast<TH2*>(fTrackerList->UncheckedAt(kResSlopeYAtVtxVsPosAbsEndIn02degMC)),
1405                                                     static_cast<TH2*>(fTrackerList->UncheckedAt(kResSlopeYAtVtxVsPosAbsEndIn23degMC)),
1406                                                     static_cast<TH2*>(fTrackerList->UncheckedAt(kResSlopeYAtVtxVsPosAbsEndIn310degMC)));
1407   fSlopeAtVtxList->AddAt(cResSlopeYAtVtxVsPosAbsEndMC, kcResSlopeYAtVtxVsPosAbsEndMC);
1408
1409   // --- diplay eta residuals at vertex ---
1410   TCanvas* cResEtaAtVtx = DrawVsAng("cResEtaAtVtx", "eta residual at vertex in 3 angular regions",
1411                                     static_cast<TH1*>(fTrackerList->UncheckedAt(kResEtaAtVtx)),
1412                                     static_cast<TH2*>(fTrackerList->UncheckedAt(kResEtaAtVtxVsAngleAtAbsEnd)));
1413   fEtaAtVtxList->AddAt(cResEtaAtVtx, kcResEtaAtVtx);
1414   TCanvas* cResEtaAtVtxMC = DrawVsAng("cResEtaAtVtxMC", "eta residual at vertex in 3 MC angular regions",
1415                                       static_cast<TH1*>(fTrackerList->UncheckedAt(kResEtaAtVtx)),
1416                                       static_cast<TH2*>(fTrackerList->UncheckedAt(kResEtaAtVtxVsMCAngle)));
1417   fEtaAtVtxList->AddAt(cResEtaAtVtxMC, kcResEtaAtVtxMC);
1418   TCanvas* cResEtaAtVtxVsPosAbsEndMC = DrawVsPos("cResEtaAtVtxVsPosAbsEndMC", "eta residual at vertex versus position at absorber end in 3 MC angular regions",
1419                                                  static_cast<TH2*>(fTrackerList->UncheckedAt(kResEtaAtVtxVsPosAbsEndIn02degMC)),
1420                                                  static_cast<TH2*>(fTrackerList->UncheckedAt(kResEtaAtVtxVsPosAbsEndIn23degMC)),
1421                                                  static_cast<TH2*>(fTrackerList->UncheckedAt(kResEtaAtVtxVsPosAbsEndIn310degMC)));
1422   fEtaAtVtxList->AddAt(cResEtaAtVtxVsPosAbsEndMC, kcResEtaAtVtxVsPosAbsEndMC);
1423
1424   // --- diplay phi residuals at vertex ---
1425   TCanvas* cResPhiAtVtx = DrawVsAng("cResPhiAtVtx", "phi residual at vertex in 3 angular regions",
1426                                     static_cast<TH1*>(fTrackerList->UncheckedAt(kResPhiAtVtx)),
1427                                     static_cast<TH2*>(fTrackerList->UncheckedAt(kResPhiAtVtxVsAngleAtAbsEnd)));
1428   fPhiAtVtxList->AddAt(cResPhiAtVtx, kcResPhiAtVtx);
1429   TCanvas* cResPhiAtVtxMC = DrawVsAng("cResPhiAtVtxMC", "phi residual at vertex in 3 MC angular regions",
1430                                       static_cast<TH1*>(fTrackerList->UncheckedAt(kResPhiAtVtx)),
1431                                       static_cast<TH2*>(fTrackerList->UncheckedAt(kResPhiAtVtxVsMCAngle)));
1432   fPhiAtVtxList->AddAt(cResPhiAtVtxMC, kcResPhiAtVtxMC);
1433   TCanvas* cResPhiAtVtxVsPosAbsEndMC = DrawVsPos("cResPhiAtVtxVsPosAbsEndMC", "phi residual at vertex versus position at absorber end in 3 MC angular regions",
1434                                                  static_cast<TH2*>(fTrackerList->UncheckedAt(kResPhiAtVtxVsPosAbsEndIn02degMC)),
1435                                                  static_cast<TH2*>(fTrackerList->UncheckedAt(kResPhiAtVtxVsPosAbsEndIn23degMC)),
1436                                                  static_cast<TH2*>(fTrackerList->UncheckedAt(kResPhiAtVtxVsPosAbsEndIn310degMC)));
1437   fPhiAtVtxList->AddAt(cResPhiAtVtxVsPosAbsEndMC, kcResPhiAtVtxVsPosAbsEndMC);
1438
1439   // --- diplay P*DCA ---
1440   TCanvas* cPDCA = DrawVsAng("cPDCA", "p #times DCA in 3 angular regions",
1441                              static_cast<TH1*>(fTrackerList->UncheckedAt(kPDCA)),
1442                              static_cast<TH2*>(fTrackerList->UncheckedAt(kPDCAVsAngleAtAbsEnd)));
1443   fDCAList->AddAt(cPDCA, kcPDCA);
1444   TCanvas* cPDCAMC = DrawVsAng("cPDCAMC", "p #times DCA in 3 MC angular regions",
1445                                static_cast<TH1*>(fTrackerList->UncheckedAt(kPDCA)),
1446                                static_cast<TH2*>(fTrackerList->UncheckedAt(kPDCAVsMCAngle)));
1447   fDCAList->AddAt(cPDCAMC, kcPDCAMC);
1448   TCanvas* cPDCAVsPosAbsEndMC = DrawVsPos("cPDCAVsPosAbsEndMC", "p #times DCA versus position at absorber end in 3 MC angular regions",
1449                                   static_cast<TH2*>(fTrackerList->UncheckedAt(kPDCAVsPosAbsEndIn02degMC)),
1450                                   static_cast<TH2*>(fTrackerList->UncheckedAt(kPDCAVsPosAbsEndIn23degMC)),
1451                                   static_cast<TH2*>(fTrackerList->UncheckedAt(kPDCAVsPosAbsEndIn310degMC)));
1452   fDCAList->AddAt(cPDCAVsPosAbsEndMC, kcPDCAVsPosAbsEndMC);
1453
1454   // post param containers
1455   PostData(2, fEfficiencyList);
1456   PostData(5, fPAtVtxList);
1457   PostData(6, fSlopeAtVtxList);
1458   PostData(7, fEtaAtVtxList);
1459   PostData(8, fPhiAtVtxList);
1460   PostData(9, fPAt1stClList);
1461   PostData(10, fSlopeAt1stClList);
1462   PostData(11, fDCAList);
1463   PostData(12, fClusterList);
1464 }
1465
1466
1467 //________________________________________________________________________
1468 Bool_t AliAnalysisTaskMuonPerformance::GetEfficiency(AliCFEffGrid* efficiency, Double_t& calcEff, Double_t& calcEffErr)
1469 {
1470   //
1471   /// Calculate the efficiency when cuts on the THnSparse are applied
1472   //
1473
1474   Bool_t isGood = kTRUE;
1475
1476   TH1* histo = 0x0;
1477   Double_t sum[2] = {0., 0.};
1478   for ( Int_t ihisto=0; ihisto<2; ihisto++ ) {
1479     histo = ( ihisto==0 ) ? efficiency->GetNum()->Project(kVarCharge) : efficiency->GetDen()->Project(kVarCharge);
1480     sum[ihisto] = histo->Integral();
1481     delete histo;
1482   }
1483
1484   if ( sum[1] == 0. ) isGood = kFALSE;
1485
1486   calcEff = ( isGood ) ? sum[0]/sum[1] : 0.;
1487   if ( calcEff > 1. ) isGood = kFALSE;
1488
1489   calcEffErr = ( isGood ) ? TMath::Sqrt(calcEff*(1-calcEff)/sum[1]) : 0.;
1490
1491   return isGood;
1492 }
1493
1494 //________________________________________________________________________
1495 Int_t AliAnalysisTaskMuonPerformance::RecoTrackMother(AliMCParticle* mcParticle)
1496 {
1497   //
1498   /// Find track mother from kinematics
1499   //
1500
1501   if ( ! mcParticle )
1502     return kUnknownPart;
1503
1504   Int_t recoPdg = mcParticle->PdgCode();
1505
1506   // Track is not a muon
1507   if ( TMath::Abs(recoPdg) != 13 ) return kRecoHadron;
1508
1509   Int_t imother = mcParticle->GetMother();
1510
1511   Bool_t isFirstMotherHF = kFALSE;
1512   Int_t step = 0;
1513
1514   //printf("\n"); // REMEMBER TO CUT
1515   while ( imother >= 0 ) {
1516     TParticle* part = static_cast<AliMCParticle*>(fMCEvent->GetTrack(imother))->Particle();
1517     //printf("Particle %s  -> %s\n", part->GetName(), TMCProcessName[part->GetUniqueID()]); // REMEMBER TO CUT
1518
1519     if ( part->GetUniqueID() == kPHadronic )
1520       return kSecondaryMu;
1521
1522     Int_t absPdg = TMath::Abs(part->GetPdgCode());
1523
1524     step++;
1525     if ( step == 1 ) // only for first mother
1526       isFirstMotherHF = ( ( absPdg >= 400  && absPdg < 600 ) ||
1527                           ( absPdg >= 4000 && absPdg < 6000 ) );
1528
1529     if ( absPdg < 4 )
1530       return kPrimaryMu;
1531
1532     if ( isFirstMotherHF) {
1533       if ( absPdg == 4 )
1534         return kCharmMu;
1535       else if ( absPdg == 5 )
1536         return kBeautyMu;
1537     }
1538
1539     imother = part->GetFirstMother();
1540   }
1541
1542   return kPrimaryMu;
1543
1544 }
1545
1546 //________________________________________________________________________
1547 Float_t AliAnalysisTaskMuonPerformance::GetBinThetaAbsEnd(Float_t RAtAbsEnd, Bool_t isTheta)
1548 {
1549   //
1550   /// Get bin of theta at absorber end region
1551   //
1552   Float_t thetaDeg = ( isTheta ) ? RAtAbsEnd : TMath::ATan( RAtAbsEnd / 505. );
1553   thetaDeg *= TMath::RadToDeg();
1554   if ( thetaDeg < 2. )
1555     return 0.;
1556   else if ( thetaDeg < 3. )
1557     return 1.;
1558   else if ( thetaDeg < 9. )
1559     return 2.;
1560
1561   return 3.;
1562 }
1563
1564 //________________________________________________________________________
1565 void AliAnalysisTaskMuonPerformance::FillContainerInfo(Double_t* containerInput, AliESDMuonTrack* esdTrack, Int_t mcID)
1566 {
1567   //
1568   /// Fill container info (except kVarMatchMC)
1569   //
1570
1571   AliMCParticle* mcPart = 0x0;
1572   if (mcID >= 0) mcPart = static_cast<AliMCParticle*>(fMCEvent->GetTrack(mcID));
1573
1574   AliVParticle* track = esdTrack;
1575   if (!track) track = mcPart;
1576
1577   if (track) {
1578     containerInput[kVarPt] = track->Pt();
1579     containerInput[kVarEta] = track->Eta();
1580     containerInput[kVarPhi] = track->Phi();
1581     containerInput[kVarThetaZones] = (esdTrack) ? GetBinThetaAbsEnd(esdTrack->GetRAtAbsorberEnd()) : GetBinThetaAbsEnd(TMath::Pi()-track->Theta(),kTRUE);
1582     containerInput[kVarCharge] = (esdTrack) ? static_cast<Double_t>(track->Charge()) : static_cast<Double_t>(track->Charge())/3.;
1583     containerInput[kVarHasTracker] = (esdTrack) ? static_cast<Double_t>(esdTrack->ContainTrackerData()) : 0.;
1584     containerInput[kVarTrigger] = (esdTrack) ? static_cast<Double_t>(esdTrack->GetMatchTrigger()) : 0.;
1585     containerInput[kVarMotherType] = static_cast<Double_t>(RecoTrackMother(mcPart));
1586   }
1587
1588   if (esdTrack) esdTrack->SetLabel(mcID);
1589 }
1590
1591 //________________________________________________________________________
1592 Double_t langaufun(Double_t *x, Double_t *par) {
1593
1594   //Fit parameters:
1595   //par[0]=Width (scale) parameter of Landau density
1596   //par[1]=Most Probable (MP, location) parameter of Landau density
1597   //par[2]=Total area (integral -inf to inf, normalization constant)
1598   //par[3]=Width (sigma) of convoluted Gaussian function
1599   //
1600   //In the Landau distribution (represented by the CERNLIB approximation),
1601   //the maximum is located at x=-0.22278298 with the location parameter=0.
1602   //This shift is corrected within this function, so that the actual
1603   //maximum is identical to the MP parameter.
1604
1605   // Numeric constants
1606   Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
1607   Double_t mpshift  = -0.22278298;       // Landau maximum location
1608
1609   // Control constants
1610   Double_t np = 100.0; // number of convolution steps
1611   Double_t sc = 5.0;   // convolution extends to +-sc Gaussian sigmas
1612
1613   // Variables
1614   Double_t xx;
1615   Double_t mpc;
1616   Double_t fland;
1617   Double_t sum = 0.0;
1618   Double_t xlow,xupp;
1619   Double_t step;
1620   Double_t i;
1621
1622
1623   // MP shift correction
1624   mpc = par[1] - mpshift * par[0];
1625
1626   // Range of convolution integral
1627   xlow = x[0] - sc * par[3];
1628   xupp = x[0] + sc * par[3];
1629
1630   step = (xupp-xlow) / np;
1631
1632   // Convolution integral of Landau and Gaussian by sum
1633   for(i=1.0; i<=np/2; i++) {
1634     xx = xlow + (i-.5) * step;
1635     //change x -> -x because the tail of the Landau is at the left here...
1636     fland = TMath::Landau(-xx,mpc,par[0]) / par[0];
1637     sum += fland * TMath::Gaus(x[0],xx,par[3]);
1638
1639     //change x -> -x because the tail of the Landau is at the left here...
1640     xx = xupp - (i-.5) * step;
1641     fland = TMath::Landau(-xx,mpc,par[0]) / par[0];
1642     sum += fland * TMath::Gaus(x[0],xx,par[3]);
1643   }
1644
1645   return (par[2] * step * sum * invsq2pi / par[3]);
1646 }
1647
1648 //________________________________________________________________________
1649 void AliAnalysisTaskMuonPerformance::FitLandauGausResVsP(TH2* h, const char* fitting, TGraphAsymmErrors* gMean,
1650                                                          TGraphAsymmErrors* gMostProb, TGraphAsymmErrors* gSigma)
1651 {
1652   /// generic function to fit residuals versus momentum with a landau convoluted with a gaussian
1653
1654   static TF1* fLandauGaus = 0x0;
1655   if (!fLandauGaus) fLandauGaus = new TF1("fLandauGaus",langaufun,h->GetYaxis()->GetBinLowEdge(1),h->GetYaxis()->GetBinLowEdge(h->GetNbinsY()+1),4);
1656
1657   Int_t rebinFactorX = TMath::Max(h->GetNbinsX()/fNPBins, 1);
1658   for (Int_t i = rebinFactorX; i <= h->GetNbinsX(); i+=rebinFactorX) {
1659
1660     cout<<Form("\rFitting %s... %d/%d",fitting,i/rebinFactorX,fNPBins)<<flush;
1661
1662     TH1D *tmp = h->ProjectionY("tmp",i-rebinFactorX+1,i,"e");
1663
1664     // first fit
1665     fLandauGaus->SetParameters(0.2,0.,(Double_t)tmp->GetEntries(),1.);
1666     tmp->Fit("fLandauGaus","WWNQ");
1667
1668     // rebin histo
1669     Double_t fwhm = fLandauGaus->GetParameter(0);
1670     Double_t sigma = fLandauGaus->GetParameter(3);
1671     Double_t sigmaP = TMath::Sqrt(sigma*sigma + fwhm*fwhm/(8.*log(2.)));
1672     Int_t rebin = TMath::Max(static_cast<Int_t>(0.5*sigmaP/tmp->GetBinWidth(1)),1);
1673     while (tmp->GetNbinsX()%rebin!=0) rebin--;
1674     tmp->Rebin(rebin);
1675
1676     // second fit
1677     tmp->Fit("fLandauGaus","NQ");
1678
1679     // get fit results and fill histograms
1680     fwhm = fLandauGaus->GetParameter(0);
1681     sigma = fLandauGaus->GetParameter(3);
1682     sigmaP = TMath::Sqrt(sigma*sigma + fwhm*fwhm/(8.*log(2.)));
1683     Double_t fwhmErr = fLandauGaus->GetParError(0);
1684     Double_t sigmaErr = fLandauGaus->GetParError(3);
1685     Double_t sigmaPErr = TMath::Sqrt(sigma*sigma*sigmaErr*sigmaErr + fwhm*fwhm*fwhmErr*fwhmErr/(64.*log(2.)*log(2.))) / sigmaP;
1686     h->GetXaxis()->SetRange(i-rebinFactorX+1,i);
1687     Double_t p = (tmp->GetEntries() > 0) ? h->GetMean() : 0.5 * (h->GetBinLowEdge(i-rebinFactorX+1) + h->GetBinLowEdge(i+1));
1688     h->GetXaxis()->SetRange();
1689     Double_t pErr[2] = {p-h->GetBinLowEdge(i-rebinFactorX+1), h->GetBinLowEdge(i+1)-p};
1690     gMean->SetPoint(i/rebinFactorX-1, p, tmp->GetMean());
1691     gMean->SetPointError(i/rebinFactorX-1, pErr[0], pErr[1], tmp->GetMeanError(), tmp->GetMeanError());
1692     gMostProb->SetPoint(i/rebinFactorX-1, p, -fLandauGaus->GetParameter(1));
1693     gMostProb->SetPointError(i/rebinFactorX-1, pErr[0], pErr[1], fLandauGaus->GetParError(1), fLandauGaus->GetParError(1));
1694     gSigma->SetPoint(i/rebinFactorX-1, p, sigmaP);
1695     gSigma->SetPointError(i/rebinFactorX-1, pErr[0], pErr[1], sigmaPErr, sigmaPErr);
1696
1697     // clean memory
1698     delete tmp;
1699   }
1700
1701   cout<<Form("\rFitting %s... %d/%d",fitting,fNPBins,fNPBins)<<endl;
1702
1703 }
1704
1705 //________________________________________________________________________
1706 void AliAnalysisTaskMuonPerformance::FitGausResVsMom(TH2* h, const Double_t mean0,
1707                                                      const Double_t sigma0, const char* fitting,
1708                                                      TGraphAsymmErrors* gMean, TGraphAsymmErrors* gSigma)
1709 {
1710   /// generic function to fit residuals versus momentum with a gaussian
1711
1712   static TF1* fGaus = 0x0;
1713   if (!fGaus) fGaus = new TF1("fGaus","gaus");
1714
1715   Int_t rebinFactorX = TMath::Max(h->GetNbinsX()/fNPBins, 1);
1716   for (Int_t i = rebinFactorX; i <= h->GetNbinsX(); i+=rebinFactorX) {
1717
1718     cout<<Form("\rFitting %s... %d/%d",fitting,i/rebinFactorX,fNPBins)<<flush;
1719
1720     TH1D *tmp = h->ProjectionY("tmp",i-rebinFactorX+1,i,"e");
1721
1722     // first fit
1723     fGaus->SetParameters(tmp->GetEntries(), mean0, sigma0);
1724     tmp->Fit("fGaus","WWNQ");
1725
1726     // rebin histo
1727     Int_t rebin = TMath::Max(static_cast<Int_t>(0.5*fGaus->GetParameter(2)/tmp->GetBinWidth(1)),1);
1728     while (tmp->GetNbinsX()%rebin!=0) rebin--;
1729     tmp->Rebin(rebin);
1730
1731     // second fit
1732     tmp->Fit("fGaus","NQ");
1733
1734     // get fit results and fill histograms
1735     h->GetXaxis()->SetRange(i-rebinFactorX+1,i);
1736     Double_t p = (tmp->GetEntries() > 0) ? h->GetMean() : 0.5 * (h->GetBinLowEdge(i-rebinFactorX+1) + h->GetBinLowEdge(i+1));
1737     h->GetXaxis()->SetRange();
1738     Double_t pErr[2] = {p-h->GetBinLowEdge(i-rebinFactorX+1), h->GetBinLowEdge(i+1)-p};
1739     gMean->SetPoint(i/rebinFactorX-1, p, fGaus->GetParameter(1));
1740     gMean->SetPointError(i/rebinFactorX-1, pErr[0], pErr[1], fGaus->GetParError(1), fGaus->GetParError(1));
1741     gSigma->SetPoint(i/rebinFactorX-1, p, fGaus->GetParameter(2));
1742     gSigma->SetPointError(i/rebinFactorX-1, pErr[0], pErr[1], fGaus->GetParError(2), fGaus->GetParError(2));
1743
1744     // clean memory
1745     delete tmp;
1746   }
1747
1748   cout<<Form("\rFitting %s... %d/%d",fitting,fNPBins,fNPBins)<<endl;
1749
1750 }
1751
1752 //________________________________________________________________________
1753 void AliAnalysisTaskMuonPerformance::FitPDCAVsMom(TH2* h, const char* fitting,
1754                                                   TGraphAsymmErrors* gMean, TGraphAsymmErrors* gSigma)
1755 {
1756   /// generic function to fit p*DCA distributions
1757
1758   static TF1* fPGaus = 0x0;
1759   if (!fPGaus) fPGaus = new TF1("fPGaus","x*gaus");
1760
1761   Int_t rebinFactorX = TMath::Max(h->GetNbinsX()/fNPBins, 1);
1762   for (Int_t i = rebinFactorX; i <= h->GetNbinsX(); i+=rebinFactorX) {
1763
1764     cout<<Form("\rFitting %s... %d/%d",fitting,i/rebinFactorX,fNPBins)<<flush;
1765
1766     TH1D *tmp = h->ProjectionY("tmp",i-rebinFactorX+1,i,"e");
1767
1768     // rebin histo
1769     Int_t rebin = static_cast<Int_t>(25*(tmp->GetNbinsX()/(tmp->GetBinLowEdge(tmp->GetNbinsX()+1)-tmp->GetBinLowEdge(1))));
1770     while (tmp->GetNbinsX()%rebin!=0) rebin--;
1771     tmp->Rebin(rebin);
1772
1773     // fit
1774     fPGaus->SetParameters(1.,0.,80.);
1775     tmp->Fit("fPGaus","NQ");
1776
1777     // get fit results and fill histograms
1778     h->GetXaxis()->SetRange(i-rebinFactorX+1,i);
1779     Double_t p = (tmp->GetEntries() > 0) ? h->GetMean() : 0.5 * (h->GetBinLowEdge(i-rebinFactorX+1) + h->GetBinLowEdge(i+1));
1780     h->GetXaxis()->SetRange();
1781     Double_t pErr[2] = {p-h->GetBinLowEdge(i-rebinFactorX+1), h->GetBinLowEdge(i+1)-p};
1782     gMean->SetPoint(i/rebinFactorX-1, p, fPGaus->GetParameter(1));
1783     gMean->SetPointError(i/rebinFactorX-1, pErr[0], pErr[1], fPGaus->GetParError(1), fPGaus->GetParError(1));
1784     gSigma->SetPoint(i/rebinFactorX-1, p, fPGaus->GetParameter(2));
1785     gSigma->SetPointError(i/rebinFactorX-1, pErr[0], pErr[1], fPGaus->GetParError(2), fPGaus->GetParError(2));
1786
1787     // clean memory
1788     delete tmp;
1789   }
1790
1791   cout<<Form("\rFitting %s... %d/%d",fitting,fNPBins,fNPBins)<<endl;
1792
1793 }
1794
1795 //________________________________________________________________________
1796 void AliAnalysisTaskMuonPerformance::FitClusterResidual(TH1* h, Int_t i, Double_t& sigma,
1797                                                         TGraphErrors* gMean, TGraphErrors* gSigma)
1798 {
1799   /// fill graphs with residual mean and sigma
1800
1801   static TF1* fRGaus = 0x0;
1802   Double_t mean, meanErr, sigmaErr;
1803
1804   if (fFitResiduals) {
1805
1806     if (!fRGaus) fRGaus = new TF1("fRGaus","gaus");
1807
1808     // first fit
1809     fRGaus->SetParameters(h->GetEntries(), 0., 0.1);
1810     h->Fit("fRGaus", "WWNQ");
1811
1812     // rebin histo
1813     Int_t rebin = TMath::Max(static_cast<Int_t>(0.2*fRGaus->GetParameter(2)/h->GetBinWidth(1)),1);
1814     while (h->GetNbinsX()%rebin!=0) rebin--;
1815     h->Rebin(rebin);
1816
1817     // second fit
1818     h->Fit("fRGaus","NQ");
1819
1820     mean = fRGaus->GetParameter(1);
1821     meanErr = fRGaus->GetParError(1);
1822     sigma = fRGaus->GetParameter(2);
1823     sigmaErr = fRGaus->GetParError(2);
1824
1825   } else {
1826
1827     Zoom(h);
1828     mean = h->GetMean();
1829     meanErr = h->GetMeanError();
1830     sigma = h->GetRMS();
1831     sigmaErr = h->GetRMSError();
1832     h->GetXaxis()->SetRange(0,0);
1833
1834   }
1835
1836   gMean->SetPoint(i, i+1, mean);
1837   gMean->SetPointError(i, 0., meanErr);
1838
1839   if (fCorrectForSystematics) {
1840     Double_t s = TMath::Sqrt(mean*mean + sigma*sigma);
1841     sigmaErr = (s>0.) ? TMath::Sqrt(sigma*sigma*sigmaErr*sigmaErr + mean*mean*meanErr*meanErr) / s : 0.;
1842     sigma = s;
1843   }
1844
1845   gSigma->SetPoint(i, i+1, sigma);
1846   gSigma->SetPointError(i, 0., sigmaErr);
1847
1848 }
1849
1850 //________________________________________________________________________
1851 TCanvas* AliAnalysisTaskMuonPerformance::DrawVsAng(const char* name, const char* title, TH1* h1, TH2* h2)
1852 {
1853   /// generic function to draw histograms versus absorber angular region
1854   TCanvas* c = new TCanvas(name, title);
1855   c->cd();
1856   h1->DrawCopy();
1857   TH1D *proj1 = h2->ProjectionY(Form("%s_proj_0_2",h2->GetName()),1,2);
1858   proj1->SetLineColor(2);
1859   proj1->Draw("sames");
1860   TH1D *proj2 = h2->ProjectionY(Form("%s_proj_2_3",h2->GetName()),3,3);
1861   proj2->SetLineColor(4);
1862   proj2->Draw("sames");
1863   TH1D *proj3 = h2->ProjectionY(Form("%s__proj_3_10",h2->GetName()),4,10);
1864   proj3->SetLineColor(3);
1865   proj3->Draw("sames");
1866   return c;
1867 }
1868
1869 //________________________________________________________________________
1870 TCanvas* AliAnalysisTaskMuonPerformance::DrawVsPos(const char* name, const char* title, TH2* h1, TH2* h2, TH2* h3)
1871 {
1872   /// generic function to draw histograms versus position at absorber end
1873   TCanvas* c = new TCanvas(name, title);
1874   c->cd();
1875   h1->SetMarkerColor(2);
1876   h1->DrawCopy();
1877   h2->SetMarkerColor(4);
1878   h2->DrawCopy("sames");
1879   h3->SetMarkerColor(3);
1880   h3->DrawCopy("sames");
1881   return c;
1882 }
1883
1884 //________________________________________________________________________
1885 TCanvas* AliAnalysisTaskMuonPerformance::DrawFitLandauGausResPVsP(const char* name, const char* title,
1886                                                                   TH2* h, const Int_t nBins, const char* fitting)
1887 {
1888   /// generic function to draw and fit momentum residuals versus momentum
1889
1890   static TF1* fLandauGaus = 0x0;
1891   if (!fLandauGaus) fLandauGaus = new TF1("fLandauGaus",langaufun,h->GetYaxis()->GetBinLowEdge(1),h->GetYaxis()->GetBinLowEdge(h->GetNbinsY()+1),4);
1892
1893   TLegend* l = new TLegend(0.15,0.25,0.3,0.85);
1894   TCanvas* c = new TCanvas(name, title);
1895   c->cd();
1896
1897   h->Sumw2();
1898
1899   Int_t rebinFactorX = TMath::Max(h->GetNbinsX()/nBins, 1);
1900   for (Int_t i = rebinFactorX; i <= h->GetNbinsX(); i+=rebinFactorX) {
1901
1902     cout<<Form("\rFitting %s... %d/%d",fitting,i/rebinFactorX,nBins)<<flush;
1903
1904     // draw projection
1905     TH1D* proj = h->ProjectionY(Form("%s_%d",h->GetName(),i/rebinFactorX),i-rebinFactorX+1,i);
1906     if (proj->GetEntries() > 0) proj->Scale(1./proj->GetEntries());
1907     proj->Draw((i==rebinFactorX)?"hist":"histsames");
1908     proj->SetLineColor(i/rebinFactorX);
1909
1910     // first fit
1911     fLandauGaus->SetParameters(0.2,0.,1.,1.);
1912     fLandauGaus->SetLineColor(i/rebinFactorX);
1913     proj->Fit("fLandauGaus","WWNQ","sames");
1914
1915     // rebin histo
1916     Double_t fwhm = fLandauGaus->GetParameter(0);
1917     Double_t sigma = fLandauGaus->GetParameter(3);
1918     Double_t sigmaP = TMath::Sqrt(sigma*sigma + fwhm*fwhm/(8.*log(2.)));
1919     Int_t rebin = TMath::Max(static_cast<Int_t>(0.5*sigmaP/proj->GetBinWidth(1)),1);
1920     while (proj->GetNbinsX()%rebin!=0) rebin--;
1921     proj->Rebin(rebin);
1922     proj->Scale(1./rebin);
1923
1924     // second fit
1925     proj->Fit("fLandauGaus","Q","sames");
1926
1927     // set label
1928     Double_t p = 0.5 * (h->GetBinLowEdge(i-rebinFactorX+1) + h->GetBinLowEdge(i+1));
1929     l->AddEntry(proj,Form("%5.1f GeV",p));
1930
1931   }
1932
1933   cout<<Form("\rFitting %s... %d/%d",fitting,nBins,nBins)<<endl;
1934
1935   l->Draw("same");
1936
1937   return c;
1938 }
1939
1940 //________________________________________________________________________
1941 TCanvas* AliAnalysisTaskMuonPerformance::DrawResPVsP(const char* name, const char* title, TH2* h, const Int_t nBins)
1942 {
1943   /// generic function to draw momentum residuals versus momentum
1944
1945   TLegend* l = new TLegend(0.15,0.25,0.3,0.85);
1946   TCanvas* c = new TCanvas(name, title);
1947   c->cd();
1948
1949   h->Sumw2();
1950
1951   Int_t rebinFactorX = TMath::Max(h->GetNbinsX()/nBins, 1);
1952   for (Int_t i = rebinFactorX; i <= h->GetNbinsX(); i+=rebinFactorX) {
1953
1954     // draw projection
1955     TH1D* proj = h->ProjectionY(Form("%s_%d",h->GetName(),i/rebinFactorX),i-rebinFactorX+1,i);
1956     if (proj->GetEntries() > 0) proj->Scale(1./proj->GetEntries());
1957     proj->Draw((i==rebinFactorX)?"hist":"histsames");
1958     proj->SetLineColor(i/rebinFactorX);
1959     proj->SetLineWidth(2);
1960
1961     // set label
1962     Double_t p = 0.5 * (h->GetBinLowEdge(i-rebinFactorX+1) + h->GetBinLowEdge(i+1));
1963     l->AddEntry(proj,Form("%5.1f GeV",p));
1964
1965   }
1966
1967   l->Draw("same");
1968
1969   return c;
1970 }
1971
1972 //________________________________________________________________________
1973 void AliAnalysisTaskMuonPerformance::Zoom(TH1* h, Double_t fractionCut)
1974 {
1975   /// Reduce the range of the histogram by removing a given fration of the statistic at each edge
1976
1977   Double_t maxEventsCut = fractionCut * h->GetEntries();
1978   Int_t nBins = h->GetNbinsX();
1979
1980   // set low edge
1981   Int_t minBin;
1982   Double_t eventsCut = 0.;
1983   for (minBin = 1; minBin <= nBins; minBin++) {
1984     eventsCut += h->GetBinContent(minBin);
1985     if (eventsCut > maxEventsCut) break;
1986   }
1987
1988   // set high edge
1989   Int_t maxBin;
1990   eventsCut = 0.;
1991   for (maxBin = nBins; maxBin >= 1; maxBin--) {
1992     eventsCut += h->GetBinContent(maxBin);
1993     if (eventsCut > maxEventsCut) break;
1994   }
1995
1996   // set new axis range
1997   h->GetXaxis()->SetRange(--minBin, ++maxBin);
1998 }