PWG4/JetTasks/AliPWG4HighPtQAMC.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //-----------------------------------------------------------------------
  17 // This class compares the global reconstruction with the MC information
  18 // Momentum resolution is stored as function of track cuts and pt.
  19 // Output: Histograms for different set of cuts
  20 //-----------------------------------------------------------------------
  21 // Author : Marta Verweij - UU
  22 //-----------------------------------------------------------------------
  23
  24 #ifndef ALIPWG4HighPtQAMC_CXX
  25 #define ALIPWG4HighPtQAMC_CXX
  26
  27 #include "AliPWG4HighPtQAMC.h"
  28
  29 #include "TH1.h"
  30 #include "TH2.h"
  31 #include "TH3.h"
  32 #include "TList.h"
  33 #include "TChain.h"
  34 #include "TH3F.h"
  35 #include "AliAnalysisManager.h"
  36 #include "AliESDInputHandler.h"
  37 #include "AliMCEvent.h"
  38 #include "AliMCEventHandler.h"
  39 #include "AliStack.h"
  40 #include "AliESDtrack.h"
  41 #include "AliESDtrackCuts.h"
  42 #include "AliExternalTrackParam.h"
  43 #include "AliLog.h"
  44
  45 using namespace std; //required for resolving the 'cout' symbol
  46
  47 ClassImp(AliPWG4HighPtQAMC)
  48
  49 AliPWG4HighPtQAMC::AliPWG4HighPtQAMC(): AliAnalysisTask("AliPWG4HighPtQAMC", ""),
  50   fESD(0),
  51   fTrackCuts(0),
  52   fTrackCutsITS(0),
  53   fNEventAll(0),
  54   fNEventSel(0),
  55   fPtAll(0),
  56   fPtSel(0),
  57   fPtAllminPtMCvsPtAll(0),
  58   fPtAllminPtMCvsPtAllNPointTPC(0),
  59   fPtAllminPtMCvsPtAllDCAR(0),
  60   fPtAllminPtMCvsPtAllDCAZ(0),
  61   fPtAllminPtMCvsPtAllPhi(0),
  62   fPtAllminPtMCvsPtAllNPointITS(0),
  63   fPtAllminPtMCvsPtAllNSigmaToVertex(0),
  64   fPtAllminPtMCvsPtAllChi2C(0),
  65   fPtAllminPtMCvsPtAllRel1PtUncertainty(0),
  66   fPtAllMC(0),
  67   fPtSelMC(0),
  68   fPtSelMCITS(0),
  69   fHistList(0),
  70   fPtSelITS(0),
  71   fPtITSminPtMCvsPtITS(0),
  72   fPtITSminPtMCvsPtITSNPointTPC(0),
  73   fPtITSminPtMCvsPtITSDCAR(0),
  74   fPtITSminPtMCvsPtITSDCAZ(0),
  75   fPtITSminPtMCvsPtITSPhi(0),
  76   fPtITSminPtMCvsPtITSNPointITS(0),
  77   fPtITSminPtMCvsPtITSNSigmaToVertex(0),
  78   fPtITSminPtMCvsPtITSChi2C(0),
  79   fPtITSminPtMCvsPtITSRel1PtUncertainty(0),
  80   fHistListITS(0)
  81 {
  82
  83 }
  84 //________________________________________________________________________
  85 AliPWG4HighPtQAMC::AliPWG4HighPtQAMC(const char *name):
  86   AliAnalysisTask(name, ""),
  87   fESD(0),
  88   fTrackCuts(),
  89   fTrackCutsITS(),
  90   fNEventAll(0),
  91   fNEventSel(0),
  92   fPtAll(0),
  93   fPtSel(0),
  94   fPtAllminPtMCvsPtAll(0),
  95   fPtAllminPtMCvsPtAllNPointTPC(0),
  96   fPtAllminPtMCvsPtAllDCAR(0),
  97   fPtAllminPtMCvsPtAllDCAZ(0),
  98   fPtAllminPtMCvsPtAllPhi(0),
  99   fPtAllminPtMCvsPtAllNPointITS(0),
 100   fPtAllminPtMCvsPtAllNSigmaToVertex(0),
 101   fPtAllminPtMCvsPtAllChi2C(0),
 102   fPtAllminPtMCvsPtAllRel1PtUncertainty(0),
 103   fPtAllMC(0),
 104   fPtSelMC(0),
 105   fPtSelMCITS(0),
 106   fHistList(0),
 107   fPtSelITS(0),
 108   fPtITSminPtMCvsPtITS(0),
 109   fPtITSminPtMCvsPtITSNPointTPC(0),
 110   fPtITSminPtMCvsPtITSDCAR(0),
 111   fPtITSminPtMCvsPtITSDCAZ(0),
 112   fPtITSminPtMCvsPtITSPhi(0),
 113   fPtITSminPtMCvsPtITSNPointITS(0),
 114   fPtITSminPtMCvsPtITSNSigmaToVertex(0),
 115   fPtITSminPtMCvsPtITSChi2C(0),
 116   fPtITSminPtMCvsPtITSRel1PtUncertainty(0),
 117   fHistListITS(0)
 118 {
 119   //
 120   // Constructor. Initialization of Inputs and Outputs
 121   //
 122   AliDebug(2,Form("AliPWG4HighPtQAMC","Calling Constructor"));
 123   // Input slot #0 works with a TChain ESD
 124   DefineInput(0, TChain::Class());
 125   // Output slot #0 writes into a TList
 126   DefineOutput(0, TList::Class());
 127   // Output slot #1 writes into a TList
 128   DefineOutput(1, TList::Class());
 129   // Output slot #2 writes into a TList
 130   DefineOutput(2, TList::Class());
 131 }
 132
 133 //________________________________________________________________________
 134 void AliPWG4HighPtQAMC::ConnectInputData(Option_t *)
 135 {
 136   // Connect ESD here
 137   // Called once
 138   AliDebug(2,Form(">> AliPWG4HighPtSpectra::ConnectInputData \n"));
 139   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
 140   if (!tree) {
 141     AliDebug(2,Form("ERROR: Could not read chain from input slot 0"));
 142   } else {
 143
 144     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 145
 146     if (!esdH) {
 147       AliDebug(2,Form("ERROR: Could not get ESDInputHandler"));
 148     } else
 149       fESD = esdH->GetEvent();
 150   }
 151 }
 152
 153 //________________________________________________________________________
 154 void AliPWG4HighPtQAMC::CreateOutputObjects() {
 155   //Create output objects
 156   AliDebug(2,Form(">> AliPWG4HighPtQATPConly::CreateOutputObjects \n"));
 157
 158   Bool_t oldStatus = TH1::AddDirectoryStatus();
 159   TH1::AddDirectory(kFALSE);
 160
 161   OpenFile(0);
 162   fHistList = new TList();
 163   OpenFile(1);
 164   fHistListITS = new TList();
 165
 166   Int_t fgkNPhiBins=18;
 167   Float_t kMinPhi = 0.;
 168   Float_t kMaxPhi = 2.*TMath::Pi();
 169
 170   Int_t fgkNPtBins=98;
 171   Float_t fgkPtMin=2.;
 172   Float_t fgkPtMax=100.;
 173   Int_t fgkResPtBins=80;
 174
 175   fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
 176   fHistList->Add(fNEventAll);
 177   fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
 178   fHistList->Add(fNEventSel);
 179   fPtAll = new TH1F("fPtAll","PtAll",fgkNPtBins, fgkPtMin, fgkPtMax);
 180   fHistList->Add(fPtAll);
 181   fPtSel = new TH1F("fPtSel","PtSel",fgkNPtBins, fgkPtMin, fgkPtMax);
 182   fHistList->Add(fPtSel);
 183
 184   fPtAllminPtMCvsPtAll = new TH2F("fPtAllminPtMCvsPtAll","PtAllminPtMCvsPtAll",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.);
 185   fPtAllminPtMCvsPtAll->SetXTitle("p_{t}^{MC}");
 186   fPtAllminPtMCvsPtAll->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 187   fHistList->Add(fPtAllminPtMCvsPtAll);
 188
 189   fPtAllminPtMCvsPtAllNPointTPC = new TH3F("fPtAllminPtMCvsPtAllNPointTPC","PtAllminPtMCvsPtAllNPointTPC",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,160,0.5,160.5);
 190   fPtAllminPtMCvsPtAllNPointTPC->SetXTitle("p_{t}^{MC}");
 191   fPtAllminPtMCvsPtAllNPointTPC->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 192   fPtAllminPtMCvsPtAllNPointTPC->SetZTitle("N_{point,TPC}");
 193   fHistList->Add(fPtAllminPtMCvsPtAllNPointTPC);
 194
 195   fPtAllminPtMCvsPtAllDCAR = new TH3F("fPtAllminPtMCvsPtAllDCAR","PtAllminPtMCvsPtAllDCAR",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,80,-0.2,0.2);
 196   fPtAllminPtMCvsPtAllDCAR->SetXTitle("p_{t}^{MC}");
 197   fPtAllminPtMCvsPtAllDCAR->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 198   fPtAllminPtMCvsPtAllDCAR->SetZTitle("DCA_{R}");
 199   fHistList->Add(fPtAllminPtMCvsPtAllDCAR);
 200
 201   fPtAllminPtMCvsPtAllDCAZ = new TH3F("fPtAllminPtMCvsPtAllDCAZ","PtAllminPtMCvsPtAllDCAZ",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,80,-2.,2.);
 202   fPtAllminPtMCvsPtAllDCAZ->SetXTitle("p_{t}^{MC}");
 203   fPtAllminPtMCvsPtAllDCAZ->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 204   fPtAllminPtMCvsPtAllDCAZ->SetZTitle("DCA_{Z}");
 205   fHistList->Add(fPtAllminPtMCvsPtAllDCAZ);
 206
 207   fPtAllminPtMCvsPtAllPhi = new TH3F("fPtAllminPtMCvsPtAllPhi","PtAllminPtMCvsPtAllPhi",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,fgkNPhiBins,kMinPhi,kMaxPhi);
 208   fPtAllminPtMCvsPtAllPhi->SetXTitle("p_{t}^{MC}");
 209   fPtAllminPtMCvsPtAllPhi->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 210   fPtAllminPtMCvsPtAllPhi->SetZTitle("#phi");
 211   fHistList->Add(fPtAllminPtMCvsPtAllPhi);
 212
 213   fPtAllminPtMCvsPtAllNPointITS = new TH3F("fPtAllminPtMCvsPtAllNPointITS","PtAllminPtMCvsPtAllNPointITS",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,9,-0.5,8.5);
 214   fPtAllminPtMCvsPtAllNPointITS->SetXTitle("p_{t}^{MC}");
 215   fPtAllminPtMCvsPtAllNPointITS->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 216   fPtAllminPtMCvsPtAllNPointITS->SetZTitle("N_{point,ITS}}");
 217   fHistList->Add(fPtAllminPtMCvsPtAllNPointITS);
 218
 219   fPtAllminPtMCvsPtAllNSigmaToVertex = new TH3F("fPtAllminPtMCvsPtAllNSigmaToVertex","PtAllminPtMCvsPtAllNSigmaToVertex",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,40,0.,8.);
 220   fPtAllminPtMCvsPtAllNSigmaToVertex->SetXTitle("p_{t}^{MC}");
 221   fPtAllminPtMCvsPtAllNSigmaToVertex->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 222   fPtAllminPtMCvsPtAllNSigmaToVertex->SetZTitle("N#sigma to vertex");
 223   fHistList->Add(fPtAllminPtMCvsPtAllNSigmaToVertex);
 224
 225   fPtAllminPtMCvsPtAllChi2C = new TH3F("fPtAllminPtMCvsPtAllChi2C","PtAllminPtMCvsPtAllChi2C",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,20,0.,10.);
 226   fPtAllminPtMCvsPtAllChi2C->SetXTitle("p_{t}^{MC}");
 227   fPtAllminPtMCvsPtAllChi2C->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 228   fPtAllminPtMCvsPtAllChi2C->SetZTitle("Constrained #chi^{2}");
 229   fHistList->Add(fPtAllminPtMCvsPtAllChi2C);
 230
 231   fPtAllminPtMCvsPtAllRel1PtUncertainty = new TH3F("fPtAllminPtMCvsPtAllRel1PtUncertainty","PtAllminPtMCvsPtAllRel1PtUncertainty",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,30,0.,0.3);
 232   fPtAllminPtMCvsPtAllRel1PtUncertainty->SetXTitle("p_{t}^{MC}");
 233   fPtAllminPtMCvsPtAllRel1PtUncertainty->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 234   fPtAllminPtMCvsPtAllRel1PtUncertainty->SetZTitle("Rel1PtUncertainty");
 235   fHistList->Add(fPtAllminPtMCvsPtAllRel1PtUncertainty);
 236
 237   //ITSrefit
 238   fPtSelITS = new TH1F("fPtSelITSrefit","PtSel",fgkNPtBins, fgkPtMin, fgkPtMax);
 239   fHistListITS->Add(fPtSelITS);
 240
 241   fPtITSminPtMCvsPtITS = new TH2F("fPtITSminPtMCvsPtITS","PtITSminPtMCvsPtITS",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.);
 242   fPtITSminPtMCvsPtITS->SetXTitle("p_{t}^{MC}");
 243   fPtITSminPtMCvsPtITS->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 244   fHistListITS->Add(fPtITSminPtMCvsPtITS);
 245
 246   fPtITSminPtMCvsPtITSNPointTPC = new TH3F("fPtITSminPtMCvsPtITSNPointTPC","PtITSminPtMCvsPtITSNPointTPC",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,160,0.5,160.5);
 247   fPtITSminPtMCvsPtITSNPointTPC->SetXTitle("p_{t}^{MC}");
 248   fPtITSminPtMCvsPtITSNPointTPC->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 249   fPtITSminPtMCvsPtITSNPointTPC->SetZTitle("N_{point,TPC}");
 250   fHistListITS->Add(fPtITSminPtMCvsPtITSNPointTPC);
 251
 252   fPtITSminPtMCvsPtITSDCAR = new TH3F("fPtITSminPtMCvsPtITSDCAR","PtITSminPtMCvsPtITSDCAR",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,80,-0.2,0.2);
 253   fPtITSminPtMCvsPtITSDCAR->SetXTitle("p_{t}^{MC}");
 254   fPtITSminPtMCvsPtITSDCAR->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 255   fPtITSminPtMCvsPtITSDCAR->SetZTitle("DCA_{R}");
 256   fHistListITS->Add(fPtITSminPtMCvsPtITSDCAR);
 257
 258   fPtITSminPtMCvsPtITSDCAZ = new TH3F("fPtITSminPtMCvsPtITSDCAZ","PtITSminPtMCvsPtITSDCAZ",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,80,-2.,2.);
 259   fPtITSminPtMCvsPtITSDCAZ->SetXTitle("p_{t}^{MC}");
 260   fPtITSminPtMCvsPtITSDCAZ->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 261   fPtITSminPtMCvsPtITSDCAZ->SetZTitle("DCA_{Z}");
 262   fHistListITS->Add(fPtITSminPtMCvsPtITSDCAZ);
 263
 264   fPtITSminPtMCvsPtITSPhi = new TH3F("fPtITSminPtMCvsPtITSPhi","PtITSminPtMCvsPtITSPhi",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,fgkNPhiBins,kMinPhi,kMaxPhi);
 265   fPtITSminPtMCvsPtITSPhi->SetXTitle("p_{t}^{MC}");
 266   fPtITSminPtMCvsPtITSPhi->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 267   fPtITSminPtMCvsPtITSPhi->SetZTitle("#phi");
 268   fHistListITS->Add(fPtITSminPtMCvsPtITSPhi);
 269
 270   fPtITSminPtMCvsPtITSNPointITS = new TH3F("fPtITSminPtMCvsPtITSNPointITS","PtITSminPtMCvsPtITSNPointITS",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,9,-0.5,8.5);
 271   fPtITSminPtMCvsPtITSNPointITS->SetXTitle("p_{t}^{MC}");
 272   fPtITSminPtMCvsPtITSNPointITS->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})}");
 273   fPtITSminPtMCvsPtITSNPointITS->SetZTitle("N_{point,ITS}}");
 274   fHistListITS->Add(fPtITSminPtMCvsPtITSNPointITS);
 275
 276   fPtITSminPtMCvsPtITSNSigmaToVertex = new TH3F("fPtITSminPtMCvsPtITSNSigmaToVertex","PtITSminPtMCvsPtITSNSigmaToVertex",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,40,0.,8.);
 277   fPtITSminPtMCvsPtITSNSigmaToVertex->SetXTitle("p_{t}^{MC}");
 278   fPtITSminPtMCvsPtITSNSigmaToVertex->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 279   fPtITSminPtMCvsPtITSNSigmaToVertex->SetZTitle("N#sigma to vertex");
 280   fHistListITS->Add(fPtITSminPtMCvsPtITSNSigmaToVertex);
 281
 282   fPtITSminPtMCvsPtITSChi2C = new TH3F("fPtITSminPtMCvsPtITSChi2C","PtITSminPtMCvsPtITSChi2C",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,20,0.,10.);
 283   fPtITSminPtMCvsPtITSChi2C->SetXTitle("p_{t}^{MC}");
 284   fPtITSminPtMCvsPtITSChi2C->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 285   fPtITSminPtMCvsPtITSChi2C->SetZTitle("Constrained #chi^{2}");
 286   fHistListITS->Add(fPtITSminPtMCvsPtITSChi2C);
 287
 288   fPtITSminPtMCvsPtITSRel1PtUncertainty = new TH3F("fPtITSminPtMCvsPtITSRel1PtUncertainty","PtITSminPtMCvsPtITSRel1PtUncertainty",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,30,0.,0.3);
 289   fPtITSminPtMCvsPtITSRel1PtUncertainty->SetXTitle("p_{t}^{MC}");
 290   fPtITSminPtMCvsPtITSRel1PtUncertainty->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 291   fPtITSminPtMCvsPtITSRel1PtUncertainty->SetZTitle("Rel1PtUncertainty");
 292   fHistListITS->Add(fPtITSminPtMCvsPtITSRel1PtUncertainty);
 293
 294   fPtAllMC = new TH1F("fPtAllMC","PtAll",fgkNPtBins, fgkPtMin, fgkPtMax);
 295   fHistList->Add(fPtAllMC);
 296   fPtSelMC = new TH1F("fPtSelMC","PtSel",fgkNPtBins, fgkPtMin, fgkPtMax);
 297   fHistList->Add(fPtSelMC);
 298   fPtSelMCITS = new TH1F("fPtSelMCITS","PtSel",fgkNPtBins, fgkPtMin, fgkPtMax);
 299   fHistList->Add(fPtSelMCITS);
 300
 301   TH1::AddDirectory(oldStatus);
 302
 303 }
 304 //________________________________________________________________________
 305 void AliPWG4HighPtQAMC::Exec(Option_t *) {
 306   // Main loop
 307   // Called for each event
 308   AliDebug(2,Form(">> AliPWG4HighPtQATPConly::Exec \n"));
 309
 310   // All events without selection
 311   fNEventAll->Fill(0.);
 312
 313   if (!fESD) {
 314     AliDebug(2,Form("ERROR: fESD not available"));
 315     PostData(0, fHistList);
 316     PostData(1, fHistListITS);
 317     return;
 318   }
 319
 320   Bool_t isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 321   if(!isSelected) { //Select collison candidates
 322     AliDebug(2,Form(" Trigger Selection: event REJECTED ... "));
 323     // Post output data
 324     PostData(0, fHistList);
 325     PostData(1, fHistListITS);
 326     return;
 327   }
 328
 329  AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 330   if (!eventHandler) {
 331     AliDebug(2,Form("ERROR: Could not retrieve MC event handler"));
 332     PostData(0, fHistList);
 333     PostData(1, fHistListITS);
 334     return;
 335   }
 336
 337   AliMCEvent* mcEvent = eventHandler->MCEvent();
 338   if (!mcEvent) {
 339     AliDebug(2,Form("ERROR: Could not retrieve MC event"));
 340     PostData(0, fHistList);
 341     PostData(1, fHistListITS);
 342     return;
 343   }
 344
 345   AliDebug(2,Form("MC particles: %d", mcEvent->GetNumberOfTracks()));
 346
 347   if (!fESD) {
 348     AliDebug(2,Form("ERROR: fESD not available"));
 349     PostData(0, fHistList);
 350     PostData(1, fHistListITS);
 351     return;
 352   }
 353
 354   AliStack* stack = mcEvent->Stack();                //Particles Stack
 355
 356   AliDebug(2,Form("MC particles stack: %d", stack->GetNtrack()));
 357
 358   const AliESDVertex *vtx = fESD->GetPrimaryVertex();
 359   // Need vertex cut
 360   if (vtx->GetNContributors() < 2) {
 361     PostData(0, fHistList);
 362     PostData(1, fHistListITS);
 363     return;
 364   }
 365  double primVtx[3];
 366   vtx->GetXYZ(primVtx);
 367   //  printf("primVtx: %g  %g  %g \n",primVtx[0],primVtx[1],primVtx[2]);
 368   if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
 369     // Post output data
 370     PostData(0, fHistList);
 371     PostData(1, fHistListITS);
 372     return;
 373   }
 374
 375   AliDebug(2,Form("Vertex title %s, status %d, nCont %d\n",vtx->GetTitle(), vtx->GetStatus(), vtx->GetNContributors()));
 376
 377   // Need to keep track of evts without vertex
 378   fNEventSel->Fill(0.);
 379
 380   if(!fESD->GetNumberOfTracks() || fESD->GetNumberOfTracks()<2)  {
 381     PostData(0, fHistList);
 382     PostData(1, fHistListITS);
 383     return;
 384   }
 385
 386   Int_t nTracks = fESD->GetNumberOfTracks();
 387   AliDebug(2,Form("nTracks %d", nTracks));
 388
 389   for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
 390
 391     AliESDtrack *track = fESD->GetTrack(iTrack);
 392     if(!track) continue;
 393     Int_t label = TMath::Abs(track->GetLabel());
 394     TParticle *particle = stack->Particle(label) ;
 395     if(!particle) continue;
 396
 397     Float_t pt = track->Pt();
 398     Float_t ptMC = particle->Pt();
 399     Float_t phi = track->Phi();
 400     Float_t dca2D, dcaZ;
 401     track->GetImpactParameters(dca2D,dcaZ);
 402     UChar_t itsMap = track->GetITSClusterMap();
 403     Int_t nPointITS = 0;
 404     for (Int_t i=0; i < 6; i++) {
 405       if (itsMap & (1 << i))
 406         nPointITS ++;
 407     }
 408     Float_t nSigmaToVertex = fTrackCuts->GetSigmaToVertex(track);// Calculates the number of sigma to the vertex for a track.
 409     Float_t chi2C = track->GetConstrainedChi2();
 410     Float_t relUncertainty1Pt = TMath::Sqrt(track->GetSigma1Pt2())*pt;
 411
 412     fPtAll->Fill(pt);
 413     fPtAllMC->Fill(ptMC);
 414
 415     if (fTrackCuts->AcceptTrack(track)) {
 416
 417       fPtSel->Fill(pt);
 418
 419       fPtSelMC->Fill(ptMC);
 420       fPtAllminPtMCvsPtAll->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC) );
 421       fPtAllminPtMCvsPtAllNPointTPC->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),track->GetTPCNcls());
 422       fPtAllminPtMCvsPtAllDCAR->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dca2D);
 423       fPtAllminPtMCvsPtAllDCAZ->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dcaZ);
 424       fPtAllminPtMCvsPtAllPhi->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),phi);
 425       fPtAllminPtMCvsPtAllNPointITS->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),nPointITS);
 426       fPtAllminPtMCvsPtAllNSigmaToVertex->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),nSigmaToVertex);
 427       fPtAllminPtMCvsPtAllChi2C->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),chi2C);
 428       fPtAllminPtMCvsPtAllRel1PtUncertainty->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),relUncertainty1Pt);
 429     }//fTrackCuts selection
 430
 431
 432     //ITSrefit selection
 433     if (fTrackCutsITS->AcceptTrack(track)) {
 434
 435       fPtSelITS->Fill(pt);
 436       fPtSelMCITS->Fill(ptMC);
 437       fPtITSminPtMCvsPtITS->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC) );
 438       fPtITSminPtMCvsPtITSNPointTPC->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),track->GetTPCNcls());
 439       fPtITSminPtMCvsPtITSDCAR->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dca2D);
 440       fPtITSminPtMCvsPtITSDCAZ->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dcaZ);
 441       fPtITSminPtMCvsPtITSPhi->Fill(ptMC,(pt-ptMC)/(pt),phi);
 442       fPtITSminPtMCvsPtITSNPointITS->Fill(ptMC,(pt-ptMC)/(pt),nPointITS);
 443       fPtITSminPtMCvsPtITSNSigmaToVertex->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),nSigmaToVertex);
 444       fPtITSminPtMCvsPtITSChi2C->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),chi2C);
 445       fPtITSminPtMCvsPtITSRel1PtUncertainty->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),relUncertainty1Pt);
 446     }//fTrackCutsITS loop
 447
 448   }//ESD track loop
 449
 450   // Post output data
 451   PostData(0, fHistList);
 452   PostData(1, fHistListITS);
 453
 454 }
 455 //________________________________________________________________________
 456 void AliPWG4HighPtQAMC::Terminate(Option_t *)
 457 {
 458   // The Terminate() function is the last function to be called during
 459   // a query. It always runs on the client, it can be used to present
 460   // the results graphically or save the results to file.
 461
 462 }
 463
 464 #endif