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 "TProfile.h"
  33 #include "TList.h"
  34 #include "TFile.h"
  35 #include "TChain.h"
  36 #include "TH3F.h"
  37 #include "TKey.h"
  38 #include "TSystem.h"
  39
  40 #include "AliAnalysisTask.h"
  41 #include "AliAnalysisManager.h"
  42 #include "AliESDInputHandler.h"
  43 #include "AliMCEvent.h"
  44 #include "AliMCEventHandler.h"
  45 #include "AliStack.h"
  46 #include "AliESDtrack.h"
  47 #include "AliESDtrackCuts.h"
  48 #include "AliExternalTrackParam.h"
  49 #include "AliLog.h"
  50 #include "AliGenPythiaEventHeader.h"
  51 #include "AliGenCocktailEventHeader.h"
  52 //#include "AliAnalysisHelperJetTasks.h"
  53
  54 using namespace std; //required for resolving the 'cout' symbol
  55
  56 ClassImp(AliPWG4HighPtQAMC)
  57
  58 AliPWG4HighPtQAMC::AliPWG4HighPtQAMC()
  59 : AliAnalysisTask("AliPWG4HighPtQAMC", ""),
  60   fESD(0),
  61   fMC(0),
  62   fTrackCuts(0),
  63   fTrackCutsITS(0),
  64   fTrackType(0),
  65   fPtMax(100.),
  66   fAvgTrials(1),
  67   fNEventAll(0),
  68   fNEventSel(0),
  69   fh1Xsec(0),
  70   fh1Trials(0),
  71   fh1PtHard(0),
  72   fh1PtHardTrials(0),
  73   fPtAll(0),
  74   fPtSel(0),
  75   fPtAllminPtMCvsPtAll(0),
  76   fPtAllminPtMCvsPtAllNPointTPC(0),
  77   fPtAllminPtMCvsPtAllDCAR(0),
  78   fPtAllminPtMCvsPtAllDCAZ(0),
  79   fPtAllminPtMCvsPtAllPhi(0),
  80   fPtAllminPtMCvsPtAllNPointITS(0),
  81   fPtAllminPtMCvsPtAllNSigmaToVertex(0),
  82   fPtAllminPtMCvsPtAllChi2C(0),
  83   fPtAllminPtMCvsPtAllRel1PtUncertainty(0),
  84   fPtAllMC(0),
  85   fPtSelMC(0),
  86   fPtSelMCITS(0),
  87   fHistList(0),
  88   fPtSelITS(0),
  89   fPtITSminPtMCvsPtITS(0),
  90   fPtITSminPtMCvsPtITSNPointTPC(0),
  91   fPtITSminPtMCvsPtITSDCAR(0),
  92   fPtITSminPtMCvsPtITSDCAZ(0),
  93   fPtITSminPtMCvsPtITSPhi(0),
  94   fPtITSminPtMCvsPtITSNPointITS(0),
  95   fPtITSminPtMCvsPtITSNSigmaToVertex(0),
  96   fPtITSminPtMCvsPtITSChi2C(0),
  97   fPtITSminPtMCvsPtITSRel1PtUncertainty(0),
  98   fHistListITS(0)
  99 {
 100
 101 }
 102 //________________________________________________________________________
 103 AliPWG4HighPtQAMC::AliPWG4HighPtQAMC(const char *name):
 104   AliAnalysisTask(name,""),
 105   fESD(0),
 106   fMC(0),
 107   fTrackCuts(),
 108   fTrackCutsITS(),
 109   fTrackType(0),
 110   fPtMax(100.),
 111   fAvgTrials(1),
 112   fNEventAll(0),
 113   fNEventSel(0),
 114   fh1Xsec(0),
 115   fh1Trials(0),
 116   fh1PtHard(0),
 117   fh1PtHardTrials(0),
 118   fPtAll(0),
 119   fPtSel(0),
 120   fPtAllminPtMCvsPtAll(0),
 121   fPtAllminPtMCvsPtAllNPointTPC(0),
 122   fPtAllminPtMCvsPtAllDCAR(0),
 123   fPtAllminPtMCvsPtAllDCAZ(0),
 124   fPtAllminPtMCvsPtAllPhi(0),
 125   fPtAllminPtMCvsPtAllNPointITS(0),
 126   fPtAllminPtMCvsPtAllNSigmaToVertex(0),
 127   fPtAllminPtMCvsPtAllChi2C(0),
 128   fPtAllminPtMCvsPtAllRel1PtUncertainty(0),
 129   fPtAllMC(0),
 130   fPtSelMC(0),
 131   fPtSelMCITS(0),
 132   fHistList(0),
 133   fPtSelITS(0),
 134   fPtITSminPtMCvsPtITS(0),
 135   fPtITSminPtMCvsPtITSNPointTPC(0),
 136   fPtITSminPtMCvsPtITSDCAR(0),
 137   fPtITSminPtMCvsPtITSDCAZ(0),
 138   fPtITSminPtMCvsPtITSPhi(0),
 139   fPtITSminPtMCvsPtITSNPointITS(0),
 140   fPtITSminPtMCvsPtITSNSigmaToVertex(0),
 141   fPtITSminPtMCvsPtITSChi2C(0),
 142   fPtITSminPtMCvsPtITSRel1PtUncertainty(0),
 143   fHistListITS(0)
 144 {
 145   //
 146   // Constructor. Initialization of Inputs and Outputs
 147   //
 148   AliDebug(2,Form("AliPWG4HighPtQAMC Calling Constructor"));
 149   // Input slot #0 works with a TChain ESD
 150   DefineInput(0, TChain::Class());
 151   // Output slot #0, #1 write into a TList
 152   DefineOutput(0, TList::Class());
 153   DefineOutput(1, TList::Class());
 154 }
 155
 156 //________________________________________________________________________
 157 void AliPWG4HighPtQAMC::ConnectInputData(Option_t *)
 158 {
 159   // Connect ESD and MC event handler here
 160   // Called once
 161   AliDebug(2,Form(">> AliPWG4HighPtSpectra::ConnectInputData \n"));
 162   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
 163   if (!tree) {
 164     AliDebug(2,Form( "ERROR: Could not read chain from input slot 0 \n"));
 165     return;
 166   }
 167
 168   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 169
 170   if (!esdH) {
 171     AliDebug(2,Form("ERROR: Could not get ESDInputHandler"));
 172     return;
 173   } else
 174     fESD = esdH->GetEvent();
 175
 176  AliMCEventHandler *eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 177  if (!eventHandler) {
 178    AliDebug(2,Form( "ERROR: Could not retrieve MC event handler \n"));
 179  }
 180   else
 181     fMC = eventHandler->MCEvent();
 182
 183 }
 184
 185
 186 //________________________________________________________________________
 187 void AliPWG4HighPtQAMC::CreateOutputObjects() {
 188   //Create output objects
 189   AliDebug(2,Form(">> AliPWG4HighPtQATPConly::CreateOutputObjects \n"));
 190
 191   Bool_t oldStatus = TH1::AddDirectoryStatus();
 192   TH1::AddDirectory(kFALSE);
 193
 194   OpenFile(0);
 195   fHistList = new TList();
 196   fHistList->SetOwner(kTRUE);
 197   OpenFile(1);
 198   fHistListITS = new TList();
 199   fHistListITS->SetOwner(kTRUE);
 200
 201   Int_t fgkNPhiBins=18;
 202   Float_t kMinPhi = 0.;
 203   Float_t kMaxPhi = 2.*TMath::Pi();
 204
 205   Int_t fgkNPtBins=98;
 206   Float_t fgkPtMin=2.;
 207   Float_t fgkPtMax=fPtMax;
 208   Int_t fgkResPtBins=80;
 209
 210   fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
 211   fHistList->Add(fNEventAll);
 212   fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
 213   fHistList->Add(fNEventSel);
 214
 215   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 216   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 217   fHistList->Add(fh1Xsec);
 218
 219   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 220   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 221   fHistList->Add(fh1Trials);
 222
 223   fh1PtHard       = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
 224   fHistList->Add(fh1PtHard);
 225   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
 226   fHistList->Add(fh1PtHardTrials);
 227
 228   fPtAll = new TH1F("fPtAll","PtAll",fgkNPtBins, fgkPtMin, fgkPtMax);
 229   fHistList->Add(fPtAll);
 230   fPtSel = new TH1F("fPtSel","PtSel",fgkNPtBins, fgkPtMin, fgkPtMax);
 231   fHistList->Add(fPtSel);
 232
 233   fPtAllminPtMCvsPtAll = new TH2F("fPtAllminPtMCvsPtAll","PtAllminPtMCvsPtAll",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.);
 234   fPtAllminPtMCvsPtAll->SetXTitle("p_{t}^{MC}");
 235   fPtAllminPtMCvsPtAll->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 236   fHistList->Add(fPtAllminPtMCvsPtAll);
 237
 238   fPtAllminPtMCvsPtAllNPointTPC = new TH3F("fPtAllminPtMCvsPtAllNPointTPC","PtAllminPtMCvsPtAllNPointTPC",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,160,0.5,160.5);
 239   fPtAllminPtMCvsPtAllNPointTPC->SetXTitle("p_{t}^{MC}");
 240   fPtAllminPtMCvsPtAllNPointTPC->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 241   fPtAllminPtMCvsPtAllNPointTPC->SetZTitle("N_{point,TPC}");
 242   fHistList->Add(fPtAllminPtMCvsPtAllNPointTPC);
 243
 244   fPtAllminPtMCvsPtAllDCAR = new TH3F("fPtAllminPtMCvsPtAllDCAR","PtAllminPtMCvsPtAllDCAR",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,80,-0.2,0.2);
 245   fPtAllminPtMCvsPtAllDCAR->SetXTitle("p_{t}^{MC}");
 246   fPtAllminPtMCvsPtAllDCAR->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 247   fPtAllminPtMCvsPtAllDCAR->SetZTitle("DCA_{R}");
 248   fHistList->Add(fPtAllminPtMCvsPtAllDCAR);
 249
 250   fPtAllminPtMCvsPtAllDCAZ = new TH3F("fPtAllminPtMCvsPtAllDCAZ","PtAllminPtMCvsPtAllDCAZ",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,80,-2.,2.);
 251   fPtAllminPtMCvsPtAllDCAZ->SetXTitle("p_{t}^{MC}");
 252   fPtAllminPtMCvsPtAllDCAZ->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 253   fPtAllminPtMCvsPtAllDCAZ->SetZTitle("DCA_{Z}");
 254   fHistList->Add(fPtAllminPtMCvsPtAllDCAZ);
 255
 256   fPtAllminPtMCvsPtAllPhi = new TH3F("fPtAllminPtMCvsPtAllPhi","PtAllminPtMCvsPtAllPhi",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,fgkNPhiBins,kMinPhi,kMaxPhi);
 257   fPtAllminPtMCvsPtAllPhi->SetXTitle("p_{t}^{MC}");
 258   fPtAllminPtMCvsPtAllPhi->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 259   fPtAllminPtMCvsPtAllPhi->SetZTitle("#phi");
 260   fHistList->Add(fPtAllminPtMCvsPtAllPhi);
 261
 262   fPtAllminPtMCvsPtAllNPointITS = new TH3F("fPtAllminPtMCvsPtAllNPointITS","PtAllminPtMCvsPtAllNPointITS",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,9,-0.5,8.5);
 263   fPtAllminPtMCvsPtAllNPointITS->SetXTitle("p_{t}^{MC}");
 264   fPtAllminPtMCvsPtAllNPointITS->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 265   fPtAllminPtMCvsPtAllNPointITS->SetZTitle("N_{point,ITS}}");
 266   fHistList->Add(fPtAllminPtMCvsPtAllNPointITS);
 267
 268   fPtAllminPtMCvsPtAllNSigmaToVertex = new TH3F("fPtAllminPtMCvsPtAllNSigmaToVertex","PtAllminPtMCvsPtAllNSigmaToVertex",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,40,0.,8.);
 269   fPtAllminPtMCvsPtAllNSigmaToVertex->SetXTitle("p_{t}^{MC}");
 270   fPtAllminPtMCvsPtAllNSigmaToVertex->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 271   fPtAllminPtMCvsPtAllNSigmaToVertex->SetZTitle("N#sigma to vertex");
 272   fHistList->Add(fPtAllminPtMCvsPtAllNSigmaToVertex);
 273
 274   fPtAllminPtMCvsPtAllChi2C = new TH3F("fPtAllminPtMCvsPtAllChi2C","PtAllminPtMCvsPtAllChi2C",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,20,0.,10.);
 275   fPtAllminPtMCvsPtAllChi2C->SetXTitle("p_{t}^{MC}");
 276   fPtAllminPtMCvsPtAllChi2C->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 277   fPtAllminPtMCvsPtAllChi2C->SetZTitle("Constrained #chi^{2}");
 278   fHistList->Add(fPtAllminPtMCvsPtAllChi2C);
 279
 280   fPtAllminPtMCvsPtAllRel1PtUncertainty = new TH3F("fPtAllminPtMCvsPtAllRel1PtUncertainty","PtAllminPtMCvsPtAllRel1PtUncertainty",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,30,0.,0.3);
 281   fPtAllminPtMCvsPtAllRel1PtUncertainty->SetXTitle("p_{t}^{MC}");
 282   fPtAllminPtMCvsPtAllRel1PtUncertainty->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 283   fPtAllminPtMCvsPtAllRel1PtUncertainty->SetZTitle("Rel1PtUncertainty");
 284   fHistList->Add(fPtAllminPtMCvsPtAllRel1PtUncertainty);
 285
 286   //ITSrefit
 287   fPtSelITS = new TH1F("fPtSelITSrefit","PtSel",fgkNPtBins, fgkPtMin, fgkPtMax);
 288   fHistListITS->Add(fPtSelITS);
 289
 290   fPtITSminPtMCvsPtITS = new TH2F("fPtITSminPtMCvsPtITS","PtITSminPtMCvsPtITS",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.);
 291   fPtITSminPtMCvsPtITS->SetXTitle("p_{t}^{MC}");
 292   fPtITSminPtMCvsPtITS->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 293   fHistListITS->Add(fPtITSminPtMCvsPtITS);
 294
 295   fPtITSminPtMCvsPtITSNPointTPC = new TH3F("fPtITSminPtMCvsPtITSNPointTPC","PtITSminPtMCvsPtITSNPointTPC",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,160,0.5,160.5);
 296   fPtITSminPtMCvsPtITSNPointTPC->SetXTitle("p_{t}^{MC}");
 297   fPtITSminPtMCvsPtITSNPointTPC->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 298   fPtITSminPtMCvsPtITSNPointTPC->SetZTitle("N_{point,TPC}");
 299   fHistListITS->Add(fPtITSminPtMCvsPtITSNPointTPC);
 300
 301   fPtITSminPtMCvsPtITSDCAR = new TH3F("fPtITSminPtMCvsPtITSDCAR","PtITSminPtMCvsPtITSDCAR",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,80,-0.2,0.2);
 302   fPtITSminPtMCvsPtITSDCAR->SetXTitle("p_{t}^{MC}");
 303   fPtITSminPtMCvsPtITSDCAR->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 304   fPtITSminPtMCvsPtITSDCAR->SetZTitle("DCA_{R}");
 305   fHistListITS->Add(fPtITSminPtMCvsPtITSDCAR);
 306
 307   fPtITSminPtMCvsPtITSDCAZ = new TH3F("fPtITSminPtMCvsPtITSDCAZ","PtITSminPtMCvsPtITSDCAZ",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,80,-2.,2.);
 308   fPtITSminPtMCvsPtITSDCAZ->SetXTitle("p_{t}^{MC}");
 309   fPtITSminPtMCvsPtITSDCAZ->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 310   fPtITSminPtMCvsPtITSDCAZ->SetZTitle("DCA_{Z}");
 311   fHistListITS->Add(fPtITSminPtMCvsPtITSDCAZ);
 312
 313   fPtITSminPtMCvsPtITSPhi = new TH3F("fPtITSminPtMCvsPtITSPhi","PtITSminPtMCvsPtITSPhi",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,fgkNPhiBins,kMinPhi,kMaxPhi);
 314   fPtITSminPtMCvsPtITSPhi->SetXTitle("p_{t}^{MC}");
 315   fPtITSminPtMCvsPtITSPhi->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 316   fPtITSminPtMCvsPtITSPhi->SetZTitle("#phi");
 317   fHistListITS->Add(fPtITSminPtMCvsPtITSPhi);
 318
 319   fPtITSminPtMCvsPtITSNPointITS = new TH3F("fPtITSminPtMCvsPtITSNPointITS","PtITSminPtMCvsPtITSNPointITS",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,9,-0.5,8.5);
 320   fPtITSminPtMCvsPtITSNPointITS->SetXTitle("p_{t}^{MC}");
 321   fPtITSminPtMCvsPtITSNPointITS->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})}");
 322   fPtITSminPtMCvsPtITSNPointITS->SetZTitle("N_{point,ITS}}");
 323   fHistListITS->Add(fPtITSminPtMCvsPtITSNPointITS);
 324
 325   fPtITSminPtMCvsPtITSNSigmaToVertex = new TH3F("fPtITSminPtMCvsPtITSNSigmaToVertex","PtITSminPtMCvsPtITSNSigmaToVertex",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,40,0.,8.);
 326   fPtITSminPtMCvsPtITSNSigmaToVertex->SetXTitle("p_{t}^{MC}");
 327   fPtITSminPtMCvsPtITSNSigmaToVertex->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 328   fPtITSminPtMCvsPtITSNSigmaToVertex->SetZTitle("N#sigma to vertex");
 329   fHistListITS->Add(fPtITSminPtMCvsPtITSNSigmaToVertex);
 330
 331   fPtITSminPtMCvsPtITSChi2C = new TH3F("fPtITSminPtMCvsPtITSChi2C","PtITSminPtMCvsPtITSChi2C",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,20,0.,10.);
 332   fPtITSminPtMCvsPtITSChi2C->SetXTitle("p_{t}^{MC}");
 333   fPtITSminPtMCvsPtITSChi2C->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 334   fPtITSminPtMCvsPtITSChi2C->SetZTitle("Constrained #chi^{2}");
 335   fHistListITS->Add(fPtITSminPtMCvsPtITSChi2C);
 336
 337   fPtITSminPtMCvsPtITSRel1PtUncertainty = new TH3F("fPtITSminPtMCvsPtITSRel1PtUncertainty","PtITSminPtMCvsPtITSRel1PtUncertainty",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,30,0.,0.3);
 338   fPtITSminPtMCvsPtITSRel1PtUncertainty->SetXTitle("p_{t}^{MC}");
 339   fPtITSminPtMCvsPtITSRel1PtUncertainty->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 340   fPtITSminPtMCvsPtITSRel1PtUncertainty->SetZTitle("Rel1PtUncertainty");
 341   fHistListITS->Add(fPtITSminPtMCvsPtITSRel1PtUncertainty);
 342
 343   fPtAllMC = new TH1F("fPtAllMC","PtAll",fgkNPtBins, fgkPtMin, fgkPtMax);
 344   fHistList->Add(fPtAllMC);
 345   fPtSelMC = new TH1F("fPtSelMC","PtSel",fgkNPtBins, fgkPtMin, fgkPtMax);
 346   fHistList->Add(fPtSelMC);
 347   fPtSelMCITS = new TH1F("fPtSelMCITS","PtSel",fgkNPtBins, fgkPtMin, fgkPtMax);
 348   fHistList->Add(fPtSelMCITS);
 349
 350   TH1::AddDirectory(oldStatus);
 351
 352 }
 353 //________________________________________________________________________
 354 void AliPWG4HighPtQAMC::Exec(Option_t *) {
 355   // Main loop
 356   // Called for each event
 357   AliDebug(2,Form(">> AliPWG4HighPtQATPConly::Exec \n"));
 358   // All events without selection
 359   fNEventAll->Fill(0.);
 360
 361   if (!fESD) {
 362     AliDebug(2,Form("ERROR: fInputEvent not available\n"));
 363     PostData(0, fHistList);
 364     PostData(1, fHistListITS);
 365     return;
 366   }
 367
 368   UInt_t isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 369   if(!(isSelected&AliVEvent::kMB)) { //Select collison candidates
 370     AliDebug(2,Form(" Trigger Selection: event REJECTED ... "));
 371     // Post output data
 372     PostData(0, fHistList);
 373     PostData(1, fHistListITS);
 374     return;
 375   }
 376
 377   AliStack* stack = 0x0;
 378
 379   if(fMC) {
 380     AliDebug(2,Form("MC particles: %d", fMC->GetNumberOfTracks()));
 381     stack = fMC->Stack();                //Particles Stack
 382     AliDebug(2,Form("MC particles stack: %d", stack->GetNtrack()));
 383   } else {
 384     AliDebug(2,Form("ERROR: Could not retrieve MC eventHandler"));
 385     PostData(0, fHistList);
 386     PostData(1, fHistListITS);
 387     return;
 388   }
 389
 390   //___ get MC information __________________________________________________________________
 391
 392   Double_t ptHard = 0.;
 393   Double_t nTrials = 1; // trials for MC trigger weight for real data
 394
 395   if(fMC){
 396     AliGenPythiaEventHeader*  pythiaGenHeader = GetPythiaEventHeader(fMC);
 397      if(!pythiaGenHeader){
 398        AliDebug(2,Form("ERROR: Could not retrieve AliGenPythiaEventHeader"));
 399        PostData(0, fHistList);
 400        PostData(1, fHistListITS);
 401        return;
 402      } else {
 403         nTrials = pythiaGenHeader->Trials();
 404         ptHard  = pythiaGenHeader->GetPtHard();
 405
 406         fh1PtHard->Fill(ptHard);
 407         fh1PtHardTrials->Fill(ptHard,nTrials);
 408
 409         fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
 410      }
 411    }
 412
 413   const AliESDVertex *vtx = fESD->GetPrimaryVertex();
 414   // Need vertex cut
 415   TString vtxName(vtx->GetName());
 416   if(vtx->GetNContributors() < 2 || (vtxName.Contains("TPCVertex")) ) {
 417     // SPD vertex
 418     vtx = fESD->GetPrimaryVertexSPD();
 419     if(vtx->GetNContributors()<2) {
 420       vtx = 0x0;
 421       // Post output data
 422       PostData(0, fHistList);
 423       PostData(1, fHistListITS);
 424       return;
 425     }
 426   }
 427
 428   double primVtx[3];
 429   vtx->GetXYZ(primVtx);
 430   if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
 431     // Post output data
 432     PostData(0, fHistList);
 433     PostData(1, fHistListITS);
 434     return;
 435   }
 436
 437   AliDebug(2,Form("Vertex title %s, status %d, nCont %d\n",vtx->GetTitle(), vtx->GetStatus(), vtx->GetNContributors()));
 438
 439   if(!fESD->GetNumberOfTracks() || fESD->GetNumberOfTracks()<2)  {
 440     PostData(0, fHistList);
 441     PostData(1, fHistListITS);
 442     return;
 443   }
 444
 445   //Need to keep track of selected events
 446   fNEventSel->Fill(0.);
 447
 448   Int_t nTracks = fESD->GetNumberOfTracks();
 449   AliDebug(2,Form("nTracks ESD%d", nTracks));
 450
 451   int nMCtracks = stack->GetNtrack();
 452
 453   Float_t pt      = 0.;
 454   Float_t ptMC    = 0.;
 455   Float_t phi     = 0.;
 456   Float_t dca2D   = 0.;
 457   Float_t dcaZ    = 0.;
 458   Int_t nPointITS = 0;
 459   Float_t chi2C   = 0.;
 460   Float_t nSigmaToVertex    = 0.;
 461   Float_t relUncertainty1Pt = 0.;
 462
 463   for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
 464
 465     AliESDtrack *track;
 466     AliESDtrack *esdtrack = fESD->GetTrack(iTrack);
 467     if(!esdtrack) continue;
 468     AliExternalTrackParam *trackTPC = (AliExternalTrackParam *)esdtrack->GetTPCInnerParam();
 469
 470     if(fTrackType==1)
 471       track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
 472     else
 473       track = esdtrack;
 474
 475
 476     if(!track) continue;
 477
 478     Int_t label = TMath::Abs(track->GetLabel());
 479     if(label>=nMCtracks)continue;
 480     TParticle *particle = stack->Particle(label) ;
 481     if(!particle) continue;
 482
 483     ptMC = particle->Pt();
 484
 485     if(fTrackType==0) {       //Global
 486       pt  = track->Pt();
 487       phi = track->Phi();
 488       track->GetImpactParameters(dca2D,dcaZ);
 489     }
 490     else if(fTrackType==1) {  //TPConly
 491       pt  = trackTPC->Pt();
 492       phi = trackTPC->Phi();
 493       track->GetImpactParametersTPC(dca2D,dcaZ);
 494     }
 495     else {continue;}
 496
 497
 498     UChar_t itsMap = track->GetITSClusterMap();
 499     for (Int_t i=0; i < 6; i++) {
 500       if (itsMap & (1 << i))
 501         nPointITS ++;
 502     }
 503     nSigmaToVertex = fTrackCuts->GetSigmaToVertex(track);// Calculates the number of sigma to the vertex for a track.
 504     chi2C = track->GetConstrainedChi2();
 505     relUncertainty1Pt = TMath::Sqrt(track->GetSigma1Pt2())*pt;
 506
 507     fPtAll->Fill(pt);
 508     fPtAllMC->Fill(ptMC);
 509
 510     if (fTrackCuts->AcceptTrack(track)) {
 511
 512       fPtSel->Fill(pt);
 513
 514       fPtSelMC->Fill(ptMC);
 515       fPtAllminPtMCvsPtAll->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC) );
 516       fPtAllminPtMCvsPtAllNPointTPC->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),track->GetTPCNcls());
 517       fPtAllminPtMCvsPtAllDCAR->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dca2D);
 518       fPtAllminPtMCvsPtAllDCAZ->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dcaZ);
 519       fPtAllminPtMCvsPtAllPhi->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),phi);
 520       fPtAllminPtMCvsPtAllNPointITS->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),nPointITS);
 521       fPtAllminPtMCvsPtAllNSigmaToVertex->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),nSigmaToVertex);
 522       fPtAllminPtMCvsPtAllChi2C->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),chi2C);
 523       fPtAllminPtMCvsPtAllRel1PtUncertainty->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),relUncertainty1Pt);
 524     }//fTrackCuts selection
 525
 526
 527
 528     //ITSrefit selection
 529     if (fTrackCutsITS->AcceptTrack(track) && fTrackType==0) {
 530
 531       fPtSelITS->Fill(pt);
 532       fPtSelMCITS->Fill(ptMC);
 533       fPtITSminPtMCvsPtITS->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC) );
 534       fPtITSminPtMCvsPtITSNPointTPC->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),track->GetTPCNcls());
 535       fPtITSminPtMCvsPtITSDCAR->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dca2D);
 536       fPtITSminPtMCvsPtITSDCAZ->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dcaZ);
 537       fPtITSminPtMCvsPtITSPhi->Fill(ptMC,(pt-ptMC)/(pt),phi);
 538       fPtITSminPtMCvsPtITSNPointITS->Fill(ptMC,(pt-ptMC)/(pt),nPointITS);
 539       fPtITSminPtMCvsPtITSNSigmaToVertex->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),nSigmaToVertex);
 540       fPtITSminPtMCvsPtITSChi2C->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),chi2C);
 541       fPtITSminPtMCvsPtITSRel1PtUncertainty->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),relUncertainty1Pt);
 542     }//fTrackCutsITS loop
 543
 544   }//ESD track loop
 545
 546   // Post output data
 547   PostData(0, fHistList);
 548   PostData(1, fHistListITS);
 549
 550 }
 551 //________________________________________________________________________
 552 Bool_t AliPWG4HighPtQAMC::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
 553   //
 554   // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
 555   // This is to called in Notify and should provide the path to the AOD/ESD file
 556   // Copied from AliAnalysisTaskJetSpectrum2
 557   //
 558
 559   TString file(currFile);
 560   fXsec = 0;
 561   fTrials = 1;
 562
 563   if(file.Contains("root_archive.zip#")){
 564     Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
 565     Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
 566     file.Replace(pos+1,20,"");
 567   }
 568   else {
 569     // not an archive take the basename....
 570     file.ReplaceAll(gSystem->BaseName(file.Data()),"");
 571   }
 572   Printf("%s",file.Data());
 573
 574
 575   TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")); // problem that we cannot really test the existance of a file in a archive so we have to lvie with open error message from root
 576   if(!fxsec){
 577     // next trial fetch the histgram file
 578     fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
 579     if(!fxsec){
 580         // not a severe condition but inciate that we have no information
 581       return kFALSE;
 582     }
 583     else{
 584       // find the tlist we want to be independtent of the name so use the Tkey
 585       TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
 586       if(!key){
 587         fxsec->Close();
 588         return kFALSE;
 589       }
 590       TList *list = dynamic_cast<TList*>(key->ReadObj());
 591       if(!list){
 592         fxsec->Close();
 593         return kFALSE;
 594       }
 595       fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
 596       fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
 597       fxsec->Close();
 598     }
 599   } // no tree pyxsec.root
 600   else {
 601     TTree *xtree = (TTree*)fxsec->Get("Xsection");
 602     if(!xtree){
 603       fxsec->Close();
 604       return kFALSE;
 605     }
 606     UInt_t   ntrials  = 0;
 607     Double_t  xsection  = 0;
 608     xtree->SetBranchAddress("xsection",&xsection);
 609     xtree->SetBranchAddress("ntrials",&ntrials);
 610     xtree->GetEntry(0);
 611     fTrials = ntrials;
 612     fXsec = xsection;
 613     fxsec->Close();
 614   }
 615   return kTRUE;
 616 }
 617 //________________________________________________________________________
 618 Bool_t AliPWG4HighPtQAMC::Notify()
 619 {
 620   //
 621   // Implemented Notify() to read the cross sections
 622   // and number of trials from pyxsec.root
 623   // Copied from AliAnalysisTaskJetSpectrum2
 624   //
 625
 626   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 627   Float_t xsection = 0;
 628   Float_t ftrials  = 1;
 629
 630   fAvgTrials = 1;
 631   if(tree){
 632     TFile *curfile = tree->GetCurrentFile();
 633     if (!curfile) {
 634       Error("Notify","No current file");
 635       return kFALSE;
 636     }
 637     if(!fh1Xsec||!fh1Trials){
 638       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
 639       return kFALSE;
 640     }
 641     PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
 642     fh1Xsec->Fill("<#sigma>",xsection);
 643     // construct a poor man average trials
 644     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
 645     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
 646   }
 647   return kTRUE;
 648 }
 649
 650 //________________________________________________________________________
 651 AliGenPythiaEventHeader*  AliPWG4HighPtQAMC::GetPythiaEventHeader(AliMCEvent *mcEvent){
 652
 653   if(!mcEvent)return 0;
 654   AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
 655   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
 656   if(!pythiaGenHeader){
 657     // cocktail ??
 658     AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
 659
 660     if (!genCocktailHeader) {
 661       AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
 662       //      AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
 663       return 0;
 664     }
 665     TList* headerList = genCocktailHeader->GetHeaders();
 666     for (Int_t i=0; i<headerList->GetEntries(); i++) {
 667       pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
 668       if (pythiaGenHeader)
 669         break;
 670     }
 671     if(!pythiaGenHeader){
 672       AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
 673       return 0;
 674     }
 675   }
 676   return pythiaGenHeader;
 677
 678 }
 679
 680 //________________________________________________________________________
 681 void AliPWG4HighPtQAMC::Terminate(Option_t *)
 682 {
 683   // The Terminate() function is the last function to be called during
 684   // a query. It always runs on the client, it can be used to present
 685   // the results graphically or save the results to file.
 686
 687 }
 688
 689 #endif