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   fMC(0),
  52   fTrackCuts(0),
  53   fTrackCutsITS(0),
  54   fTrackType(0),
  55   fPtMax(100.),
  56   fNEventAll(0),
  57   fNEventSel(0),
  58   fPtAll(0),
  59   fPtSel(0),
  60   fPtAllminPtMCvsPtAll(0),
  61   fPtAllminPtMCvsPtAllNPointTPC(0),
  62   fPtAllminPtMCvsPtAllDCAR(0),
  63   fPtAllminPtMCvsPtAllDCAZ(0),
  64   fPtAllminPtMCvsPtAllPhi(0),
  65   fPtAllminPtMCvsPtAllNPointITS(0),
  66   fPtAllminPtMCvsPtAllNSigmaToVertex(0),
  67   fPtAllminPtMCvsPtAllChi2C(0),
  68   fPtAllminPtMCvsPtAllRel1PtUncertainty(0),
  69   fPtAllMC(0),
  70   fPtSelMC(0),
  71   fPtSelMCITS(0),
  72   fHistList(0),
  73   fPtSelITS(0),
  74   fPtITSminPtMCvsPtITS(0),
  75   fPtITSminPtMCvsPtITSNPointTPC(0),
  76   fPtITSminPtMCvsPtITSDCAR(0),
  77   fPtITSminPtMCvsPtITSDCAZ(0),
  78   fPtITSminPtMCvsPtITSPhi(0),
  79   fPtITSminPtMCvsPtITSNPointITS(0),
  80   fPtITSminPtMCvsPtITSNSigmaToVertex(0),
  81   fPtITSminPtMCvsPtITSChi2C(0),
  82   fPtITSminPtMCvsPtITSRel1PtUncertainty(0),
  83   fHistListITS(0)
  84 {
  85
  86 }
  87 //________________________________________________________________________
  88 AliPWG4HighPtQAMC::AliPWG4HighPtQAMC(const char *name):
  89   AliAnalysisTask(name, ""),
  90   fESD(0),
  91   fMC(0),
  92   fTrackCuts(),
  93   fTrackCutsITS(),
  94   fTrackType(0),
  95   fPtMax(100.),
  96   fNEventAll(0),
  97   fNEventSel(0),
  98   fPtAll(0),
  99   fPtSel(0),
 100   fPtAllminPtMCvsPtAll(0),
 101   fPtAllminPtMCvsPtAllNPointTPC(0),
 102   fPtAllminPtMCvsPtAllDCAR(0),
 103   fPtAllminPtMCvsPtAllDCAZ(0),
 104   fPtAllminPtMCvsPtAllPhi(0),
 105   fPtAllminPtMCvsPtAllNPointITS(0),
 106   fPtAllminPtMCvsPtAllNSigmaToVertex(0),
 107   fPtAllminPtMCvsPtAllChi2C(0),
 108   fPtAllminPtMCvsPtAllRel1PtUncertainty(0),
 109   fPtAllMC(0),
 110   fPtSelMC(0),
 111   fPtSelMCITS(0),
 112   fHistList(0),
 113   fPtSelITS(0),
 114   fPtITSminPtMCvsPtITS(0),
 115   fPtITSminPtMCvsPtITSNPointTPC(0),
 116   fPtITSminPtMCvsPtITSDCAR(0),
 117   fPtITSminPtMCvsPtITSDCAZ(0),
 118   fPtITSminPtMCvsPtITSPhi(0),
 119   fPtITSminPtMCvsPtITSNPointITS(0),
 120   fPtITSminPtMCvsPtITSNSigmaToVertex(0),
 121   fPtITSminPtMCvsPtITSChi2C(0),
 122   fPtITSminPtMCvsPtITSRel1PtUncertainty(0),
 123   fHistListITS(0)
 124 {
 125   //
 126   // Constructor. Initialization of Inputs and Outputs
 127   //
 128   AliDebug(2,Form("AliPWG4HighPtQAMC Calling Constructor"));
 129   // Input slot #0 works with a TChain ESD
 130   DefineInput(0, TChain::Class());
 131   // Output slot #0 writes into a TList
 132   DefineOutput(0, TList::Class());
 133   // Output slot #1 writes into a TList
 134   DefineOutput(1, TList::Class());
 135   // Output slot #2 writes into a TList
 136   DefineOutput(2, TList::Class());
 137 }
 138
 139 //________________________________________________________________________
 140 void AliPWG4HighPtQAMC::ConnectInputData(Option_t *)
 141 {
 142   // Connect ESD and MC event handler here
 143   // Called once
 144   AliDebug(2,Form(">> AliPWG4HighPtSpectra::ConnectInputData \n"));
 145   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
 146   if (!tree) {
 147     AliDebug(2,Form( "ERROR: Could not read chain from input slot 0 \n"));
 148     return;
 149   }
 150
 151   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 152
 153   if (!esdH) {
 154     AliDebug(2,Form("ERROR: Could not get ESDInputHandler"));
 155     return;
 156   } else
 157     fESD = esdH->GetEvent();
 158
 159  AliMCEventHandler *eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 160  if (!eventHandler) {
 161    AliDebug(2,Form( "ERROR: Could not retrieve MC event handler \n"));
 162  }
 163   else
 164     fMC = eventHandler->MCEvent();
 165
 166 }
 167
 168 //________________________________________________________________________
 169 void AliPWG4HighPtQAMC::CreateOutputObjects() {
 170   //Create output objects
 171   AliDebug(2,Form(">> AliPWG4HighPtQATPConly::CreateOutputObjects \n"));
 172
 173   Bool_t oldStatus = TH1::AddDirectoryStatus();
 174   TH1::AddDirectory(kFALSE);
 175
 176   OpenFile(0);
 177   fHistList = new TList();
 178   OpenFile(1);
 179   fHistListITS = new TList();
 180
 181   Int_t fgkNPhiBins=18;
 182   Float_t kMinPhi = 0.;
 183   Float_t kMaxPhi = 2.*TMath::Pi();
 184
 185   Int_t fgkNPtBins=98;
 186   Float_t fgkPtMin=2.;
 187   Float_t fgkPtMax=fPtMax;
 188   Int_t fgkResPtBins=80;
 189
 190   fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
 191   fHistList->Add(fNEventAll);
 192   fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
 193   fHistList->Add(fNEventSel);
 194   fPtAll = new TH1F("fPtAll","PtAll",fgkNPtBins, fgkPtMin, fgkPtMax);
 195   fHistList->Add(fPtAll);
 196   fPtSel = new TH1F("fPtSel","PtSel",fgkNPtBins, fgkPtMin, fgkPtMax);
 197   fHistList->Add(fPtSel);
 198
 199   fPtAllminPtMCvsPtAll = new TH2F("fPtAllminPtMCvsPtAll","PtAllminPtMCvsPtAll",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.);
 200   fPtAllminPtMCvsPtAll->SetXTitle("p_{t}^{MC}");
 201   fPtAllminPtMCvsPtAll->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 202   fHistList->Add(fPtAllminPtMCvsPtAll);
 203
 204   fPtAllminPtMCvsPtAllNPointTPC = new TH3F("fPtAllminPtMCvsPtAllNPointTPC","PtAllminPtMCvsPtAllNPointTPC",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,160,0.5,160.5);
 205   fPtAllminPtMCvsPtAllNPointTPC->SetXTitle("p_{t}^{MC}");
 206   fPtAllminPtMCvsPtAllNPointTPC->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 207   fPtAllminPtMCvsPtAllNPointTPC->SetZTitle("N_{point,TPC}");
 208   fHistList->Add(fPtAllminPtMCvsPtAllNPointTPC);
 209
 210   fPtAllminPtMCvsPtAllDCAR = new TH3F("fPtAllminPtMCvsPtAllDCAR","PtAllminPtMCvsPtAllDCAR",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,80,-0.2,0.2);
 211   fPtAllminPtMCvsPtAllDCAR->SetXTitle("p_{t}^{MC}");
 212   fPtAllminPtMCvsPtAllDCAR->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 213   fPtAllminPtMCvsPtAllDCAR->SetZTitle("DCA_{R}");
 214   fHistList->Add(fPtAllminPtMCvsPtAllDCAR);
 215
 216   fPtAllminPtMCvsPtAllDCAZ = new TH3F("fPtAllminPtMCvsPtAllDCAZ","PtAllminPtMCvsPtAllDCAZ",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,80,-2.,2.);
 217   fPtAllminPtMCvsPtAllDCAZ->SetXTitle("p_{t}^{MC}");
 218   fPtAllminPtMCvsPtAllDCAZ->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 219   fPtAllminPtMCvsPtAllDCAZ->SetZTitle("DCA_{Z}");
 220   fHistList->Add(fPtAllminPtMCvsPtAllDCAZ);
 221
 222   fPtAllminPtMCvsPtAllPhi = new TH3F("fPtAllminPtMCvsPtAllPhi","PtAllminPtMCvsPtAllPhi",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,fgkNPhiBins,kMinPhi,kMaxPhi);
 223   fPtAllminPtMCvsPtAllPhi->SetXTitle("p_{t}^{MC}");
 224   fPtAllminPtMCvsPtAllPhi->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 225   fPtAllminPtMCvsPtAllPhi->SetZTitle("#phi");
 226   fHistList->Add(fPtAllminPtMCvsPtAllPhi);
 227
 228   fPtAllminPtMCvsPtAllNPointITS = new TH3F("fPtAllminPtMCvsPtAllNPointITS","PtAllminPtMCvsPtAllNPointITS",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,9,-0.5,8.5);
 229   fPtAllminPtMCvsPtAllNPointITS->SetXTitle("p_{t}^{MC}");
 230   fPtAllminPtMCvsPtAllNPointITS->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 231   fPtAllminPtMCvsPtAllNPointITS->SetZTitle("N_{point,ITS}}");
 232   fHistList->Add(fPtAllminPtMCvsPtAllNPointITS);
 233
 234   fPtAllminPtMCvsPtAllNSigmaToVertex = new TH3F("fPtAllminPtMCvsPtAllNSigmaToVertex","PtAllminPtMCvsPtAllNSigmaToVertex",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,40,0.,8.);
 235   fPtAllminPtMCvsPtAllNSigmaToVertex->SetXTitle("p_{t}^{MC}");
 236   fPtAllminPtMCvsPtAllNSigmaToVertex->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 237   fPtAllminPtMCvsPtAllNSigmaToVertex->SetZTitle("N#sigma to vertex");
 238   fHistList->Add(fPtAllminPtMCvsPtAllNSigmaToVertex);
 239
 240   fPtAllminPtMCvsPtAllChi2C = new TH3F("fPtAllminPtMCvsPtAllChi2C","PtAllminPtMCvsPtAllChi2C",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,20,0.,10.);
 241   fPtAllminPtMCvsPtAllChi2C->SetXTitle("p_{t}^{MC}");
 242   fPtAllminPtMCvsPtAllChi2C->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 243   fPtAllminPtMCvsPtAllChi2C->SetZTitle("Constrained #chi^{2}");
 244   fHistList->Add(fPtAllminPtMCvsPtAllChi2C);
 245
 246   fPtAllminPtMCvsPtAllRel1PtUncertainty = new TH3F("fPtAllminPtMCvsPtAllRel1PtUncertainty","PtAllminPtMCvsPtAllRel1PtUncertainty",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,30,0.,0.3);
 247   fPtAllminPtMCvsPtAllRel1PtUncertainty->SetXTitle("p_{t}^{MC}");
 248   fPtAllminPtMCvsPtAllRel1PtUncertainty->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 249   fPtAllminPtMCvsPtAllRel1PtUncertainty->SetZTitle("Rel1PtUncertainty");
 250   fHistList->Add(fPtAllminPtMCvsPtAllRel1PtUncertainty);
 251
 252   //ITSrefit
 253   fPtSelITS = new TH1F("fPtSelITSrefit","PtSel",fgkNPtBins, fgkPtMin, fgkPtMax);
 254   fHistListITS->Add(fPtSelITS);
 255
 256   fPtITSminPtMCvsPtITS = new TH2F("fPtITSminPtMCvsPtITS","PtITSminPtMCvsPtITS",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.);
 257   fPtITSminPtMCvsPtITS->SetXTitle("p_{t}^{MC}");
 258   fPtITSminPtMCvsPtITS->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 259   fHistListITS->Add(fPtITSminPtMCvsPtITS);
 260
 261   fPtITSminPtMCvsPtITSNPointTPC = new TH3F("fPtITSminPtMCvsPtITSNPointTPC","PtITSminPtMCvsPtITSNPointTPC",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,160,0.5,160.5);
 262   fPtITSminPtMCvsPtITSNPointTPC->SetXTitle("p_{t}^{MC}");
 263   fPtITSminPtMCvsPtITSNPointTPC->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 264   fPtITSminPtMCvsPtITSNPointTPC->SetZTitle("N_{point,TPC}");
 265   fHistListITS->Add(fPtITSminPtMCvsPtITSNPointTPC);
 266
 267   fPtITSminPtMCvsPtITSDCAR = new TH3F("fPtITSminPtMCvsPtITSDCAR","PtITSminPtMCvsPtITSDCAR",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,80,-0.2,0.2);
 268   fPtITSminPtMCvsPtITSDCAR->SetXTitle("p_{t}^{MC}");
 269   fPtITSminPtMCvsPtITSDCAR->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 270   fPtITSminPtMCvsPtITSDCAR->SetZTitle("DCA_{R}");
 271   fHistListITS->Add(fPtITSminPtMCvsPtITSDCAR);
 272
 273   fPtITSminPtMCvsPtITSDCAZ = new TH3F("fPtITSminPtMCvsPtITSDCAZ","PtITSminPtMCvsPtITSDCAZ",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,80,-2.,2.);
 274   fPtITSminPtMCvsPtITSDCAZ->SetXTitle("p_{t}^{MC}");
 275   fPtITSminPtMCvsPtITSDCAZ->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 276   fPtITSminPtMCvsPtITSDCAZ->SetZTitle("DCA_{Z}");
 277   fHistListITS->Add(fPtITSminPtMCvsPtITSDCAZ);
 278
 279   fPtITSminPtMCvsPtITSPhi = new TH3F("fPtITSminPtMCvsPtITSPhi","PtITSminPtMCvsPtITSPhi",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,fgkNPhiBins,kMinPhi,kMaxPhi);
 280   fPtITSminPtMCvsPtITSPhi->SetXTitle("p_{t}^{MC}");
 281   fPtITSminPtMCvsPtITSPhi->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 282   fPtITSminPtMCvsPtITSPhi->SetZTitle("#phi");
 283   fHistListITS->Add(fPtITSminPtMCvsPtITSPhi);
 284
 285   fPtITSminPtMCvsPtITSNPointITS = new TH3F("fPtITSminPtMCvsPtITSNPointITS","PtITSminPtMCvsPtITSNPointITS",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,9,-0.5,8.5);
 286   fPtITSminPtMCvsPtITSNPointITS->SetXTitle("p_{t}^{MC}");
 287   fPtITSminPtMCvsPtITSNPointITS->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})}");
 288   fPtITSminPtMCvsPtITSNPointITS->SetZTitle("N_{point,ITS}}");
 289   fHistListITS->Add(fPtITSminPtMCvsPtITSNPointITS);
 290
 291   fPtITSminPtMCvsPtITSNSigmaToVertex = new TH3F("fPtITSminPtMCvsPtITSNSigmaToVertex","PtITSminPtMCvsPtITSNSigmaToVertex",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,40,0.,8.);
 292   fPtITSminPtMCvsPtITSNSigmaToVertex->SetXTitle("p_{t}^{MC}");
 293   fPtITSminPtMCvsPtITSNSigmaToVertex->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 294   fPtITSminPtMCvsPtITSNSigmaToVertex->SetZTitle("N#sigma to vertex");
 295   fHistListITS->Add(fPtITSminPtMCvsPtITSNSigmaToVertex);
 296
 297   fPtITSminPtMCvsPtITSChi2C = new TH3F("fPtITSminPtMCvsPtITSChi2C","PtITSminPtMCvsPtITSChi2C",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,20,0.,10.);
 298   fPtITSminPtMCvsPtITSChi2C->SetXTitle("p_{t}^{MC}");
 299   fPtITSminPtMCvsPtITSChi2C->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 300   fPtITSminPtMCvsPtITSChi2C->SetZTitle("Constrained #chi^{2}");
 301   fHistListITS->Add(fPtITSminPtMCvsPtITSChi2C);
 302
 303   fPtITSminPtMCvsPtITSRel1PtUncertainty = new TH3F("fPtITSminPtMCvsPtITSRel1PtUncertainty","PtITSminPtMCvsPtITSRel1PtUncertainty",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,30,0.,0.3);
 304   fPtITSminPtMCvsPtITSRel1PtUncertainty->SetXTitle("p_{t}^{MC}");
 305   fPtITSminPtMCvsPtITSRel1PtUncertainty->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 306   fPtITSminPtMCvsPtITSRel1PtUncertainty->SetZTitle("Rel1PtUncertainty");
 307   fHistListITS->Add(fPtITSminPtMCvsPtITSRel1PtUncertainty);
 308
 309   fPtAllMC = new TH1F("fPtAllMC","PtAll",fgkNPtBins, fgkPtMin, fgkPtMax);
 310   fHistList->Add(fPtAllMC);
 311   fPtSelMC = new TH1F("fPtSelMC","PtSel",fgkNPtBins, fgkPtMin, fgkPtMax);
 312   fHistList->Add(fPtSelMC);
 313   fPtSelMCITS = new TH1F("fPtSelMCITS","PtSel",fgkNPtBins, fgkPtMin, fgkPtMax);
 314   fHistList->Add(fPtSelMCITS);
 315
 316   TH1::AddDirectory(oldStatus);
 317
 318 }
 319 //________________________________________________________________________
 320 void AliPWG4HighPtQAMC::Exec(Option_t *) {
 321   // Main loop
 322   // Called for each event
 323   AliDebug(2,Form(">> AliPWG4HighPtQATPConly::Exec \n"));
 324   // All events without selection
 325   fNEventAll->Fill(0.);
 326
 327   if (!fESD) {
 328     AliDebug(2,Form("ERROR: fESD not available\n"));
 329     PostData(0, fHistList);
 330     PostData(1, fHistListITS);
 331     return;
 332   }
 333
 334   UInt_t isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 335   if(!(isSelected&AliVEvent::kMB)) { //Select collison candidates
 336     AliDebug(2,Form(" Trigger Selection: event REJECTED ... "));
 337     // Post output data
 338     PostData(0, fHistList);
 339     PostData(1, fHistListITS);
 340     return;
 341   }
 342
 343   AliStack* stack = 0x0;
 344
 345   if(fMC) {
 346     AliDebug(2,Form("MC particles: %d", fMC->GetNumberOfTracks()));
 347     stack = fMC->Stack();                //Particles Stack
 348     AliDebug(2,Form("MC particles stack: %d", stack->GetNtrack()));
 349   } else {
 350     AliDebug(2,Form("ERROR: Could not retrieve MC eventHandler"));
 351     PostData(0, fHistList);
 352     PostData(1, fHistListITS);
 353     return;
 354   }
 355
 356   const AliESDVertex *vtx = fESD->GetPrimaryVertex();
 357   // Need vertex cut
 358   TString vtxName(vtx->GetName());
 359   if(vtx->GetNContributors() < 2 || (vtxName.Contains("TPCVertex")) ) {
 360     // SPD vertex
 361     vtx = fESD->GetPrimaryVertexSPD();
 362     if(vtx->GetNContributors()<2) {
 363       vtx = 0x0;
 364       // Post output data
 365       PostData(0, fHistList);
 366       PostData(1, fHistListITS);
 367       return;
 368     }
 369   }
 370
 371   double primVtx[3];
 372   vtx->GetXYZ(primVtx);
 373   if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
 374     // Post output data
 375     PostData(0, fHistList);
 376     PostData(1, fHistListITS);
 377     return;
 378   }
 379
 380   AliDebug(2,Form("Vertex title %s, status %d, nCont %d\n",vtx->GetTitle(), vtx->GetStatus(), vtx->GetNContributors()));
 381
 382   // Need to keep track of evts without vertex
 383   fNEventSel->Fill(0.);
 384
 385   if(!fESD->GetNumberOfTracks() || fESD->GetNumberOfTracks()<2)  {
 386     PostData(0, fHistList);
 387     PostData(1, fHistListITS);
 388     return;
 389   }
 390
 391   Int_t nTracks = fESD->GetNumberOfTracks();
 392   AliDebug(2,Form("nTracks ESD%d", nTracks));
 393
 394   int nMCtracks = stack->GetNtrack();
 395
 396   Float_t pt      = 0.;
 397   Float_t ptMC    = 0.;
 398   Float_t phi     = 0.;
 399   Float_t dca2D   = 0.;
 400   Float_t dcaZ    = 0.;
 401   Int_t nPointITS = 0;
 402   Float_t chi2C   = 0.;
 403   Float_t nSigmaToVertex    = 0.;
 404   Float_t relUncertainty1Pt = 0.;
 405
 406   for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
 407
 408     AliESDtrack *track;
 409     AliESDtrack *esdtrack = fESD->GetTrack(iTrack);
 410     if(!esdtrack) continue;
 411     AliExternalTrackParam *trackTPC = (AliExternalTrackParam *)esdtrack->GetTPCInnerParam();
 412
 413     if(fTrackType==1)
 414       track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
 415     else
 416       track = esdtrack;
 417
 418
 419     if(!track) continue;
 420
 421     Int_t label = TMath::Abs(track->GetLabel());
 422     if(label>=nMCtracks)continue;
 423     TParticle *particle = stack->Particle(label) ;
 424     if(!particle) continue;
 425
 426     ptMC = particle->Pt();
 427
 428     if(fTrackType==0) {       //Global
 429       pt  = track->Pt();
 430       phi = track->Phi();
 431       track->GetImpactParameters(dca2D,dcaZ);
 432     }
 433     else if(fTrackType==1) {  //TPConly
 434       pt  = trackTPC->Pt();
 435       phi = trackTPC->Phi();
 436       track->GetImpactParametersTPC(dca2D,dcaZ);
 437     }
 438     else {continue;}
 439
 440
 441     UChar_t itsMap = track->GetITSClusterMap();
 442     for (Int_t i=0; i < 6; i++) {
 443       if (itsMap & (1 << i))
 444         nPointITS ++;
 445     }
 446     nSigmaToVertex = fTrackCuts->GetSigmaToVertex(track);// Calculates the number of sigma to the vertex for a track.
 447     chi2C = track->GetConstrainedChi2();
 448     relUncertainty1Pt = TMath::Sqrt(track->GetSigma1Pt2())*pt;
 449
 450     fPtAll->Fill(pt);
 451     fPtAllMC->Fill(ptMC);
 452
 453     if (fTrackCuts->AcceptTrack(track)) {
 454
 455       fPtSel->Fill(pt);
 456
 457       fPtSelMC->Fill(ptMC);
 458       fPtAllminPtMCvsPtAll->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC) );
 459       fPtAllminPtMCvsPtAllNPointTPC->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),track->GetTPCNcls());
 460       fPtAllminPtMCvsPtAllDCAR->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dca2D);
 461       fPtAllminPtMCvsPtAllDCAZ->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dcaZ);
 462       fPtAllminPtMCvsPtAllPhi->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),phi);
 463       fPtAllminPtMCvsPtAllNPointITS->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),nPointITS);
 464       fPtAllminPtMCvsPtAllNSigmaToVertex->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),nSigmaToVertex);
 465       fPtAllminPtMCvsPtAllChi2C->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),chi2C);
 466       fPtAllminPtMCvsPtAllRel1PtUncertainty->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),relUncertainty1Pt);
 467     }//fTrackCuts selection
 468
 469
 470
 471     //ITSrefit selection
 472     if (fTrackCutsITS->AcceptTrack(track) && fTrackType==0) {
 473
 474       fPtSelITS->Fill(pt);
 475       fPtSelMCITS->Fill(ptMC);
 476       fPtITSminPtMCvsPtITS->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC) );
 477       fPtITSminPtMCvsPtITSNPointTPC->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),track->GetTPCNcls());
 478       fPtITSminPtMCvsPtITSDCAR->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dca2D);
 479       fPtITSminPtMCvsPtITSDCAZ->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dcaZ);
 480       fPtITSminPtMCvsPtITSPhi->Fill(ptMC,(pt-ptMC)/(pt),phi);
 481       fPtITSminPtMCvsPtITSNPointITS->Fill(ptMC,(pt-ptMC)/(pt),nPointITS);
 482       fPtITSminPtMCvsPtITSNSigmaToVertex->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),nSigmaToVertex);
 483       fPtITSminPtMCvsPtITSChi2C->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),chi2C);
 484       fPtITSminPtMCvsPtITSRel1PtUncertainty->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),relUncertainty1Pt);
 485     }//fTrackCutsITS loop
 486
 487   }//ESD track loop
 488
 489   // Post output data
 490   PostData(0, fHistList);
 491   PostData(1, fHistListITS);
 492
 493 }
 494 //________________________________________________________________________
 495 void AliPWG4HighPtQAMC::Terminate(Option_t *)
 496 {
 497   // The Terminate() function is the last function to be called during
 498   // a query. It always runs on the client, it can be used to present
 499   // the results graphically or save the results to file.
 500
 501 }
 502
 503 #endif