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