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   fStack(0),
  63   fVtx(0x0),
  64   fTrackCuts(0),
  65   fTrackCutsITS(0),
  66   fTrackType(0),
  67   fPtMax(100.),
  68   fAvgTrials(1),
  69   fNEventAll(0),
  70   fNEventSel(0),
  71   fNEventReject(0),
  72   fh1Xsec(0),
  73   fh1Trials(0),
  74   fh1PtHard(0),
  75   fh1PtHardTrials(0),
  76   fPtAll(0),
  77   fPtSel(0),
  78   fPtSelFakes(0),
  79   fNPointTPCFakes(0),
  80   fPtSelLargeLabel(0),
  81   fMultRec(0),
  82   fNPointTPCMultRec(0),
  83   fDeltaPtMultRec(0),
  84   fPtAllminPtMCvsPtAll(0),
  85   fPtAllminPtMCvsPtAllNPointTPC(0),
  86   fPtAllminPtMCvsPtAllDCAR(0),
  87   fPtAllminPtMCvsPtAllDCAZ(0),
  88   fPtAllminPtMCvsPtAllPhi(0),
  89   fPtAllminPtMCvsPtAllNPointITS(0),
  90   fPtAllminPtMCvsPtAllNSigmaToVertex(0),
  91   fPtAllminPtMCvsPtAllChi2C(0),
  92   fPtAllminPtMCvsPtAllRel1PtUncertainty(0),
  93   fPtAllMC(0),
  94   fPtSelMC(0),
  95   fPtSelMCITS(0),
  96   fHistList(0),
  97   fPtSelITS(0),
  98   fPtITSminPtMCvsPtITS(0),
  99   fPtITSminPtMCvsPtITSNPointTPC(0),
 100   fPtITSminPtMCvsPtITSDCAR(0),
 101   fPtITSminPtMCvsPtITSDCAZ(0),
 102   fPtITSminPtMCvsPtITSPhi(0),
 103   fPtITSminPtMCvsPtITSNPointITS(0),
 104   fPtITSminPtMCvsPtITSNSigmaToVertex(0),
 105   fPtITSminPtMCvsPtITSChi2C(0),
 106   fPtITSminPtMCvsPtITSRel1PtUncertainty(0),
 107   fHistListITS(0)
 108 {
 109
 110 }
 111 //________________________________________________________________________
 112 AliPWG4HighPtQAMC::AliPWG4HighPtQAMC(const char *name):
 113   AliAnalysisTask(name,""),
 114   fESD(0),
 115   fMC(0),
 116   fStack(0),
 117   fVtx(0x0),
 118   fTrackCuts(),
 119   fTrackCutsITS(),
 120   fTrackType(0),
 121   fPtMax(100.),
 122   fAvgTrials(1),
 123   fNEventAll(0),
 124   fNEventSel(0),
 125   fNEventReject(0),
 126   fh1Xsec(0),
 127   fh1Trials(0),
 128   fh1PtHard(0),
 129   fh1PtHardTrials(0),
 130   fPtAll(0),
 131   fPtSel(0),
 132   fPtSelFakes(0),
 133   fNPointTPCFakes(0),
 134   fPtSelLargeLabel(0),
 135   fMultRec(0),
 136   fNPointTPCMultRec(0),
 137   fDeltaPtMultRec(0),
 138   fPtAllminPtMCvsPtAll(0),
 139   fPtAllminPtMCvsPtAllNPointTPC(0),
 140   fPtAllminPtMCvsPtAllDCAR(0),
 141   fPtAllminPtMCvsPtAllDCAZ(0),
 142   fPtAllminPtMCvsPtAllPhi(0),
 143   fPtAllminPtMCvsPtAllNPointITS(0),
 144   fPtAllminPtMCvsPtAllNSigmaToVertex(0),
 145   fPtAllminPtMCvsPtAllChi2C(0),
 146   fPtAllminPtMCvsPtAllRel1PtUncertainty(0),
 147   fPtAllMC(0),
 148   fPtSelMC(0),
 149   fPtSelMCITS(0),
 150   fHistList(0),
 151   fPtSelITS(0),
 152   fPtITSminPtMCvsPtITS(0),
 153   fPtITSminPtMCvsPtITSNPointTPC(0),
 154   fPtITSminPtMCvsPtITSDCAR(0),
 155   fPtITSminPtMCvsPtITSDCAZ(0),
 156   fPtITSminPtMCvsPtITSPhi(0),
 157   fPtITSminPtMCvsPtITSNPointITS(0),
 158   fPtITSminPtMCvsPtITSNSigmaToVertex(0),
 159   fPtITSminPtMCvsPtITSChi2C(0),
 160   fPtITSminPtMCvsPtITSRel1PtUncertainty(0),
 161   fHistListITS(0)
 162 {
 163   //
 164   // Constructor. Initialization of Inputs and Outputs
 165   //
 166   AliDebug(2,Form("AliPWG4HighPtQAMC Calling Constructor"));
 167   // Input slot #0 works with a TChain ESD
 168   DefineInput(0, TChain::Class());
 169   // Output slot #0, #1 write into a TList
 170   DefineOutput(0, TList::Class());
 171   DefineOutput(1, TList::Class());
 172 }
 173
 174 //________________________________________________________________________
 175 void AliPWG4HighPtQAMC::ConnectInputData(Option_t *)
 176 {
 177   // Connect ESD and MC event handler here
 178   // Called once
 179   AliDebug(2,Form(">> AliPWG4HighPtSpectra::ConnectInputData \n"));
 180   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
 181   if (!tree) {
 182     AliDebug(2,Form( "ERROR: Could not read chain from input slot 0 \n"));
 183     return;
 184   }
 185
 186   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 187
 188   if (!esdH) {
 189     AliDebug(2,Form("ERROR: Could not get ESDInputHandler"));
 190     return;
 191   } else
 192     fESD = esdH->GetEvent();
 193
 194   AliMCEventHandler *eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 195   if (!eventHandler) {
 196     AliDebug(2,Form( "ERROR: Could not retrieve MC event handler \n"));
 197   }
 198   else
 199     fMC = eventHandler->MCEvent();
 200
 201 }
 202
 203
 204 //________________________________________________________________________
 205 void AliPWG4HighPtQAMC::CreateOutputObjects() {
 206   //Create output objects
 207   AliDebug(2,Form(">> AliPWG4HighPtQATPConly::CreateOutputObjects \n"));
 208
 209   Bool_t oldStatus = TH1::AddDirectoryStatus();
 210   TH1::AddDirectory(kFALSE);
 211
 212   OpenFile(0);
 213   fHistList = new TList();
 214   fHistList->SetOwner(kTRUE);
 215   OpenFile(1);
 216   fHistListITS = new TList();
 217   fHistListITS->SetOwner(kTRUE);
 218
 219   Int_t fgkNPhiBins=18;
 220   Float_t kMinPhi = 0.;
 221   Float_t kMaxPhi = 2.*TMath::Pi();
 222
 223   Int_t fgkNPtBins=98;
 224   Float_t fgkPtMin=2.;
 225   Float_t fgkPtMax=fPtMax;
 226   Int_t fgkResPtBins=80;
 227
 228   fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
 229   fHistList->Add(fNEventAll);
 230   fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
 231   fHistList->Add(fNEventSel);
 232   fNEventReject = new TH1F("fNEventReject","Reason events are rejectected for analysis",20,0,20);
 233   //Set labels
 234   fNEventReject->Fill("noESD",0);
 235   fNEventReject->Fill("Trigger",0);
 236   fNEventReject->Fill("noMCEvent",0);
 237   fNEventReject->Fill("NTracks<2",0);
 238   fNEventReject->Fill("noVTX",0);
 239   fNEventReject->Fill("VtxStatus",0);
 240   fNEventReject->Fill("NCont<2",0);
 241   fNEventReject->Fill("ZVTX>10",0);
 242   fHistList->Add(fNEventReject);
 243
 244   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 245   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 246   fHistList->Add(fh1Xsec);
 247
 248   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 249   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 250   fHistList->Add(fh1Trials);
 251
 252   fh1PtHard       = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
 253   fHistList->Add(fh1PtHard);
 254   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
 255   fHistList->Add(fh1PtHardTrials);
 256
 257   fPtAll = new TH1F("fPtAll","PtAll",fgkNPtBins, fgkPtMin, fgkPtMax);
 258   fHistList->Add(fPtAll);
 259   fPtSel = new TH1F("fPtSel","PtSel",fgkNPtBins, fgkPtMin, fgkPtMax);
 260   fHistList->Add(fPtSel);
 261
 262   fPtSelFakes = new TH1F("fPtSelFakes","PtSelFakes",fgkNPtBins, fgkPtMin, fgkPtMax);
 263   fHistList->Add(fPtSelFakes);
 264
 265   fNPointTPCFakes = new TH1F("fNPointTPCFakes","fNPointTPCFakes",160,0.5,160.5);
 266   fHistList->Add(fNPointTPCFakes);
 267
 268   fPtSelLargeLabel = new TH1F("fPtSelLargeLabel","PtSelLargeLabel",fgkNPtBins, fgkPtMin, fgkPtMax);
 269   fHistList->Add(fPtSelLargeLabel);
 270
 271   fMultRec = new TH1F("fMultRec","Multiple reconstruction of tracks",20,0,20);
 272   fHistList->Add(fMultRec);
 273
 274   fNPointTPCMultRec = new TH1F("fNPointTPCMultRec","Multiple reconstruction of tracks",160,0.5,160.5);
 275   fHistList->Add(fNPointTPCMultRec);
 276
 277   fDeltaPtMultRec = new TH2F("fDeltaPtMultRec","Delta pT vs pT for multiple reconstructed tracks",100,0.,50.,100,-20.,20.);
 278   fHistList->Add(fDeltaPtMultRec);
 279
 280   fPtAllminPtMCvsPtAll = new TH2F("fPtAllminPtMCvsPtAll","PtAllminPtMCvsPtAll",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.);
 281   fPtAllminPtMCvsPtAll->SetXTitle("p_{t}^{MC}");
 282   fPtAllminPtMCvsPtAll->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 283   fHistList->Add(fPtAllminPtMCvsPtAll);
 284
 285   fPtAllminPtMCvsPtAllNPointTPC = new TH3F("fPtAllminPtMCvsPtAllNPointTPC","PtAllminPtMCvsPtAllNPointTPC",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,160,0.5,160.5);
 286   fPtAllminPtMCvsPtAllNPointTPC->SetXTitle("p_{t}^{MC}");
 287   fPtAllminPtMCvsPtAllNPointTPC->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 288   fPtAllminPtMCvsPtAllNPointTPC->SetZTitle("N_{point,TPC}");
 289   fHistList->Add(fPtAllminPtMCvsPtAllNPointTPC);
 290
 291   fPtAllminPtMCvsPtAllDCAR = new TH3F("fPtAllminPtMCvsPtAllDCAR","PtAllminPtMCvsPtAllDCAR",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,80,-0.2,0.2);
 292   fPtAllminPtMCvsPtAllDCAR->SetXTitle("p_{t}^{MC}");
 293   fPtAllminPtMCvsPtAllDCAR->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 294   fPtAllminPtMCvsPtAllDCAR->SetZTitle("DCA_{R}");
 295   fHistList->Add(fPtAllminPtMCvsPtAllDCAR);
 296
 297   fPtAllminPtMCvsPtAllDCAZ = new TH3F("fPtAllminPtMCvsPtAllDCAZ","PtAllminPtMCvsPtAllDCAZ",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,80,-2.,2.);
 298   fPtAllminPtMCvsPtAllDCAZ->SetXTitle("p_{t}^{MC}");
 299   fPtAllminPtMCvsPtAllDCAZ->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 300   fPtAllminPtMCvsPtAllDCAZ->SetZTitle("DCA_{Z}");
 301   fHistList->Add(fPtAllminPtMCvsPtAllDCAZ);
 302
 303   fPtAllminPtMCvsPtAllPhi = new TH3F("fPtAllminPtMCvsPtAllPhi","PtAllminPtMCvsPtAllPhi",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,fgkNPhiBins,kMinPhi,kMaxPhi);
 304   fPtAllminPtMCvsPtAllPhi->SetXTitle("p_{t}^{MC}");
 305   fPtAllminPtMCvsPtAllPhi->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 306   fPtAllminPtMCvsPtAllPhi->SetZTitle("#phi");
 307   fHistList->Add(fPtAllminPtMCvsPtAllPhi);
 308
 309   fPtAllminPtMCvsPtAllNPointITS = new TH3F("fPtAllminPtMCvsPtAllNPointITS","PtAllminPtMCvsPtAllNPointITS",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,9,-0.5,8.5);
 310   fPtAllminPtMCvsPtAllNPointITS->SetXTitle("p_{t}^{MC}");
 311   fPtAllminPtMCvsPtAllNPointITS->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 312   fPtAllminPtMCvsPtAllNPointITS->SetZTitle("N_{point,ITS}}");
 313   fHistList->Add(fPtAllminPtMCvsPtAllNPointITS);
 314
 315   fPtAllminPtMCvsPtAllNSigmaToVertex = new TH3F("fPtAllminPtMCvsPtAllNSigmaToVertex","PtAllminPtMCvsPtAllNSigmaToVertex",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,40,0.,8.);
 316   fPtAllminPtMCvsPtAllNSigmaToVertex->SetXTitle("p_{t}^{MC}");
 317   fPtAllminPtMCvsPtAllNSigmaToVertex->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 318   fPtAllminPtMCvsPtAllNSigmaToVertex->SetZTitle("N#sigma to vertex");
 319   fHistList->Add(fPtAllminPtMCvsPtAllNSigmaToVertex);
 320
 321   fPtAllminPtMCvsPtAllChi2C = new TH3F("fPtAllminPtMCvsPtAllChi2C","PtAllminPtMCvsPtAllChi2C",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,20,0.,10.);
 322   fPtAllminPtMCvsPtAllChi2C->SetXTitle("p_{t}^{MC}");
 323   fPtAllminPtMCvsPtAllChi2C->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 324   fPtAllminPtMCvsPtAllChi2C->SetZTitle("Constrained #chi^{2}");
 325   fHistList->Add(fPtAllminPtMCvsPtAllChi2C);
 326
 327   fPtAllminPtMCvsPtAllRel1PtUncertainty = new TH3F("fPtAllminPtMCvsPtAllRel1PtUncertainty","PtAllminPtMCvsPtAllRel1PtUncertainty",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,30,0.,0.3);
 328   fPtAllminPtMCvsPtAllRel1PtUncertainty->SetXTitle("p_{t}^{MC}");
 329   fPtAllminPtMCvsPtAllRel1PtUncertainty->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 330   fPtAllminPtMCvsPtAllRel1PtUncertainty->SetZTitle("Rel1PtUncertainty");
 331   fHistList->Add(fPtAllminPtMCvsPtAllRel1PtUncertainty);
 332
 333   //ITSrefit
 334   fPtSelITS = new TH1F("fPtSelITSrefit","PtSel",fgkNPtBins, fgkPtMin, fgkPtMax);
 335   fHistListITS->Add(fPtSelITS);
 336
 337   fPtITSminPtMCvsPtITS = new TH2F("fPtITSminPtMCvsPtITS","PtITSminPtMCvsPtITS",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.);
 338   fPtITSminPtMCvsPtITS->SetXTitle("p_{t}^{MC}");
 339   fPtITSminPtMCvsPtITS->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 340   fHistListITS->Add(fPtITSminPtMCvsPtITS);
 341
 342   fPtITSminPtMCvsPtITSNPointTPC = new TH3F("fPtITSminPtMCvsPtITSNPointTPC","PtITSminPtMCvsPtITSNPointTPC",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,160,0.5,160.5);
 343   fPtITSminPtMCvsPtITSNPointTPC->SetXTitle("p_{t}^{MC}");
 344   fPtITSminPtMCvsPtITSNPointTPC->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 345   fPtITSminPtMCvsPtITSNPointTPC->SetZTitle("N_{point,TPC}");
 346   fHistListITS->Add(fPtITSminPtMCvsPtITSNPointTPC);
 347
 348   fPtITSminPtMCvsPtITSDCAR = new TH3F("fPtITSminPtMCvsPtITSDCAR","PtITSminPtMCvsPtITSDCAR",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,80,-0.2,0.2);
 349   fPtITSminPtMCvsPtITSDCAR->SetXTitle("p_{t}^{MC}");
 350   fPtITSminPtMCvsPtITSDCAR->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 351   fPtITSminPtMCvsPtITSDCAR->SetZTitle("DCA_{R}");
 352   fHistListITS->Add(fPtITSminPtMCvsPtITSDCAR);
 353
 354   fPtITSminPtMCvsPtITSDCAZ = new TH3F("fPtITSminPtMCvsPtITSDCAZ","PtITSminPtMCvsPtITSDCAZ",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,80,-2.,2.);
 355   fPtITSminPtMCvsPtITSDCAZ->SetXTitle("p_{t}^{MC}");
 356   fPtITSminPtMCvsPtITSDCAZ->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 357   fPtITSminPtMCvsPtITSDCAZ->SetZTitle("DCA_{Z}");
 358   fHistListITS->Add(fPtITSminPtMCvsPtITSDCAZ);
 359
 360   fPtITSminPtMCvsPtITSPhi = new TH3F("fPtITSminPtMCvsPtITSPhi","PtITSminPtMCvsPtITSPhi",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,fgkNPhiBins,kMinPhi,kMaxPhi);
 361   fPtITSminPtMCvsPtITSPhi->SetXTitle("p_{t}^{MC}");
 362   fPtITSminPtMCvsPtITSPhi->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 363   fPtITSminPtMCvsPtITSPhi->SetZTitle("#phi");
 364   fHistListITS->Add(fPtITSminPtMCvsPtITSPhi);
 365
 366   fPtITSminPtMCvsPtITSNPointITS = new TH3F("fPtITSminPtMCvsPtITSNPointITS","PtITSminPtMCvsPtITSNPointITS",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,9,-0.5,8.5);
 367   fPtITSminPtMCvsPtITSNPointITS->SetXTitle("p_{t}^{MC}");
 368   fPtITSminPtMCvsPtITSNPointITS->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})}");
 369   fPtITSminPtMCvsPtITSNPointITS->SetZTitle("N_{point,ITS}}");
 370   fHistListITS->Add(fPtITSminPtMCvsPtITSNPointITS);
 371
 372   fPtITSminPtMCvsPtITSNSigmaToVertex = new TH3F("fPtITSminPtMCvsPtITSNSigmaToVertex","PtITSminPtMCvsPtITSNSigmaToVertex",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,40,0.,8.);
 373   fPtITSminPtMCvsPtITSNSigmaToVertex->SetXTitle("p_{t}^{MC}");
 374   fPtITSminPtMCvsPtITSNSigmaToVertex->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 375   fPtITSminPtMCvsPtITSNSigmaToVertex->SetZTitle("N#sigma to vertex");
 376   fHistListITS->Add(fPtITSminPtMCvsPtITSNSigmaToVertex);
 377
 378   fPtITSminPtMCvsPtITSChi2C = new TH3F("fPtITSminPtMCvsPtITSChi2C","PtITSminPtMCvsPtITSChi2C",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,20,0.,10.);
 379   fPtITSminPtMCvsPtITSChi2C->SetXTitle("p_{t}^{MC}");
 380   fPtITSminPtMCvsPtITSChi2C->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 381   fPtITSminPtMCvsPtITSChi2C->SetZTitle("Constrained #chi^{2}");
 382   fHistListITS->Add(fPtITSminPtMCvsPtITSChi2C);
 383
 384   fPtITSminPtMCvsPtITSRel1PtUncertainty = new TH3F("fPtITSminPtMCvsPtITSRel1PtUncertainty","PtITSminPtMCvsPtITSRel1PtUncertainty",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,30,0.,0.3);
 385   fPtITSminPtMCvsPtITSRel1PtUncertainty->SetXTitle("p_{t}^{MC}");
 386   fPtITSminPtMCvsPtITSRel1PtUncertainty->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 387   fPtITSminPtMCvsPtITSRel1PtUncertainty->SetZTitle("Rel1PtUncertainty");
 388   fHistListITS->Add(fPtITSminPtMCvsPtITSRel1PtUncertainty);
 389
 390   fPtAllMC = new TH1F("fPtAllMC","PtAll",fgkNPtBins, fgkPtMin, fgkPtMax);
 391   fHistList->Add(fPtAllMC);
 392   fPtSelMC = new TH1F("fPtSelMC","PtSel",fgkNPtBins, fgkPtMin, fgkPtMax);
 393   fHistList->Add(fPtSelMC);
 394   fPtSelMCITS = new TH1F("fPtSelMCITS","PtSel",fgkNPtBins, fgkPtMin, fgkPtMax);
 395   fHistList->Add(fPtSelMCITS);
 396
 397   TH1::AddDirectory(oldStatus);
 398
 399 }
 400
 401 //________________________________________________________________________
 402 Bool_t AliPWG4HighPtQAMC::SelectEvent() {
 403   //
 404   // Decide if event should be selected for analysis
 405   //
 406
 407   // Checks following requirements:
 408   // - fESD available
 409   // - trigger info from AliPhysicsSelection
 410   // - MCevent available
 411   // - number of reconstructed tracks > 1
 412   // - primary vertex reconstructed
 413   // - z-vertex < 10 cm
 414
 415   Bool_t selectEvent = kTRUE;
 416
 417   //fESD object available?
 418   if (!fESD) {
 419     AliDebug(2,Form("ERROR: fInputEvent not available\n"));
 420     fNEventReject->Fill("noESD",1);
 421     selectEvent = kFALSE;
 422     return selectEvent;
 423   }
 424
 425   //Trigger
 426   UInt_t isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 427   if(!(isSelected&AliVEvent::kMB)) { //Select collison candidates
 428     AliDebug(2,Form(" Trigger Selection: event REJECTED ... "));
 429     fNEventReject->Fill("Trigger",1);
 430     selectEvent = kFALSE;
 431     return selectEvent;
 432   }
 433
 434   //MCEvent available?
 435   //if yes: get stack
 436   if(fMC) {
 437     AliDebug(2,Form("MC particles: %d", fMC->GetNumberOfTracks()));
 438     fStack = fMC->Stack();                //Particles Stack
 439     AliDebug(2,Form("MC particles stack: %d", fStack->GetNtrack()));
 440   } else {
 441     AliDebug(2,Form("ERROR: Could not retrieve MC eventHandler"));
 442     fNEventReject->Fill("noMCEvent",1);
 443     selectEvent = kFALSE;
 444     return selectEvent;
 445   }
 446
 447   //Check if number of reconstructed tracks is larger than 1
 448   if(!fESD->GetNumberOfTracks() || fESD->GetNumberOfTracks()<2)  {
 449     fNEventReject->Fill("NTracks<2",1);
 450     selectEvent = kFALSE;
 451     return selectEvent;
 452   }
 453
 454   //Check if vertex is reconstructed
 455   if(fTrackType==1) fVtx = fESD->GetPrimaryVertexTPC();
 456   else              fVtx = fESD->GetPrimaryVertexSPD();
 457
 458   if(!fVtx) {
 459     fNEventReject->Fill("noVTX",1);
 460     selectEvent = kFALSE;
 461     return selectEvent;
 462   }
 463
 464   if(!fVtx->GetStatus()) {
 465     fNEventReject->Fill("VtxStatus",1);
 466     selectEvent = kFALSE;
 467     return selectEvent;
 468   }
 469
 470   // Need vertex cut
 471   //  TString vtxName(fVtx->GetName());
 472   if(fVtx->GetNContributors()<2) {
 473     fNEventReject->Fill("NCont<2",1);
 474     selectEvent = kFALSE;
 475     return selectEvent;
 476   }
 477
 478   //Check if z-vertex < 10 cm
 479   double primVtx[3];
 480   fVtx->GetXYZ(primVtx);
 481   if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
 482     fNEventReject->Fill("ZVTX>10",1);
 483     selectEvent = kFALSE;
 484     return selectEvent;
 485   }
 486
 487   AliDebug(2,Form("Vertex title %s, status %d, nCont %d\n",fVtx->GetTitle(), fVtx->GetStatus(), fVtx->GetNContributors()));
 488
 489   return selectEvent;
 490
 491 }
 492
 493 //________________________________________________________________________
 494 void AliPWG4HighPtQAMC::Exec(Option_t *) {
 495   // Main loop
 496   // Called for each event
 497   AliDebug(2,Form(">> AliPWG4HighPtQATPConly::Exec \n"));
 498   // All events without selection
 499   fNEventAll->Fill(0.);
 500
 501   if(!SelectEvent()) {
 502     // Post output data
 503     PostData(0, fHistList);
 504     PostData(1, fHistListITS);
 505     return;
 506   }
 507
 508   // ---- Get MC Header information (for MC productions in pThard bins) ----
 509   Double_t ptHard = 0.;
 510   Double_t nTrials = 1; // trials for MC trigger weight for real data
 511
 512   if(fMC){
 513     AliGenPythiaEventHeader*  pythiaGenHeader = GetPythiaEventHeader(fMC);
 514      if(pythiaGenHeader){
 515        nTrials = pythiaGenHeader->Trials();
 516        ptHard  = pythiaGenHeader->GetPtHard();
 517
 518        fh1PtHard->Fill(ptHard);
 519        fh1PtHardTrials->Fill(ptHard,nTrials);
 520
 521        fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
 522      }
 523   }
 524
 525   //Need to keep track of selected events
 526   fNEventSel->Fill(0.);
 527
 528   Int_t nTracks = fESD->GetNumberOfTracks();
 529   AliDebug(2,Form("nTracks ESD%d", nTracks));
 530
 531   int nMCtracks = fStack->GetNtrack();
 532
 533   Float_t pt      = 0.;
 534   Float_t ptMC    = 0.;
 535   Float_t phi     = 0.;
 536   Float_t dca2D   = 0.;
 537   Float_t dcaZ    = 0.;
 538   Int_t nPointITS = 0;
 539   Float_t chi2C   = 0.;
 540   Float_t nSigmaToVertex    = 0.;
 541   Float_t relUncertainty1Pt = 0.;
 542
 543   for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
 544
 545     AliESDtrack *track;
 546     AliESDtrack *esdtrack = fESD->GetTrack(iTrack);
 547     if(!esdtrack) continue;
 548
 549     if(fTrackType==1)
 550       track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
 551     else if(fTrackType==2) {
 552       track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
 553       if(!track) {
 554         delete track;
 555         continue;
 556       }
 557       AliExternalTrackParam exParam;
 558       Bool_t relate = track->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
 559       if( !relate ) {
 560         delete track;
 561         continue;
 562       }
 563       track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
 564     }
 565     else
 566       track = esdtrack;
 567
 568     if(!track) {
 569       if(fTrackType==1 || fTrackType==2) delete track;
 570       continue;
 571     }
 572
 573     Int_t label = TMath::Abs(track->GetLabel());
 574     if(label>=nMCtracks) {
 575       if (fTrackCuts->AcceptTrack(track)) { fPtSelLargeLabel->Fill(pt); }
 576       if(fTrackType==1 || fTrackType==2) delete track;
 577       continue;
 578     }
 579     TParticle *particle = fStack->Particle(label) ;
 580     if(!particle) {
 581       if(fTrackType==1 || fTrackType==2) delete track;
 582       continue;
 583     }
 584
 585     ptMC = particle->Pt();
 586
 587     if(fTrackType==0) {       //Global
 588       pt  = track->Pt();
 589       phi = track->Phi();
 590       track->GetImpactParameters(dca2D,dcaZ);
 591     }
 592     else if(fTrackType==1 || fTrackType==2) {  //TPConly
 593       pt  = track->Pt();
 594       phi = track->Phi();
 595       track->GetImpactParametersTPC(dca2D,dcaZ);
 596     }
 597     else {continue;}
 598
 599
 600     UChar_t itsMap = track->GetITSClusterMap();
 601     for (Int_t i=0; i < 6; i++) {
 602       if (itsMap & (1 << i))
 603         nPointITS ++;
 604     }
 605     nSigmaToVertex = fTrackCuts->GetSigmaToVertex(track);// Calculates the number of sigma to the vertex for a track.
 606     chi2C = track->GetConstrainedChi2();
 607     relUncertainty1Pt = TMath::Sqrt(track->GetSigma1Pt2())*pt;
 608
 609     fPtAll->Fill(pt);
 610     fPtAllMC->Fill(ptMC);
 611
 612     if (fTrackCuts->AcceptTrack(track)) {
 613
 614       fPtSel->Fill(pt);
 615       if(track->GetLabel()<0) {
 616         fPtSelFakes->Fill(pt);
 617         fNPointTPCFakes->Fill(track->GetTPCNcls());
 618       }
 619       fPtSelMC->Fill(ptMC);
 620       fPtAllminPtMCvsPtAll->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC) );
 621       fPtAllminPtMCvsPtAllNPointTPC->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),track->GetTPCNcls());
 622       fPtAllminPtMCvsPtAllDCAR->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dca2D);
 623       fPtAllminPtMCvsPtAllDCAZ->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dcaZ);
 624       fPtAllminPtMCvsPtAllPhi->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),phi);
 625       fPtAllminPtMCvsPtAllNPointITS->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),nPointITS);
 626       fPtAllminPtMCvsPtAllNSigmaToVertex->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),nSigmaToVertex);
 627       fPtAllminPtMCvsPtAllChi2C->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),chi2C);
 628       fPtAllminPtMCvsPtAllRel1PtUncertainty->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),relUncertainty1Pt);
 629
 630
 631       //Check if track is reconstructed multiple times
 632       int multCounter = 1;
 633       for (Int_t iTrack2 = iTrack+1; iTrack2 < nTracks; iTrack2++) {
 634         //   AliESDtrack *track2 = GetTrackForAnalysis(iTrack2);
 635         AliESDtrack *track2;
 636         AliESDtrack *esdtrack2 = fESD->GetTrack(iTrack2);
 637         if(!esdtrack2) continue;
 638
 639         if(fTrackType==1)
 640           track2 = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack2->GetID());
 641         else if(fTrackType==2) {
 642           track2 = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack2->GetID());
 643           if(!track2) {
 644             delete track2;
 645             continue;
 646           }
 647           AliExternalTrackParam exParam2;
 648           Bool_t relate = track2->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam2);
 649           if( !relate ) {
 650             delete track2;
 651             continue;
 652           }
 653           track2->Set(exParam2.GetX(),exParam2.GetAlpha(),exParam2.GetParameter(),exParam2.GetCovariance());
 654         }
 655         else
 656           track2 = esdtrack2;
 657         if(!track2) {
 658           if(fTrackType==1 || fTrackType==2) delete track2;
 659           continue;
 660         }
 661
 662         if (fTrackCuts->AcceptTrack(track2)) {
 663           Int_t label2 = TMath::Abs(track2->GetLabel());
 664           if(label==label2) {
 665             fNPointTPCMultRec->Fill(track->GetTPCNcls());
 666             fNPointTPCMultRec->Fill(track2->GetTPCNcls());
 667             fDeltaPtMultRec->Fill(track->Pt(),track->Pt()-track2->Pt());
 668             multCounter++;
 669           }
 670         }
 671         if(fTrackType==1 || fTrackType==2) delete track2;
 672       }//track2 loop
 673       if(multCounter>1) fMultRec->Fill(multCounter);
 674
 675     }//fTrackCuts selection
 676
 677     //ITSrefit selection
 678     if (fTrackCutsITS->AcceptTrack(track) && fTrackType==0) {
 679
 680       fPtSelITS->Fill(pt);
 681       fPtSelMCITS->Fill(ptMC);
 682       fPtITSminPtMCvsPtITS->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC) );
 683       fPtITSminPtMCvsPtITSNPointTPC->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),track->GetTPCNcls());
 684       fPtITSminPtMCvsPtITSDCAR->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dca2D);
 685       fPtITSminPtMCvsPtITSDCAZ->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dcaZ);
 686       fPtITSminPtMCvsPtITSPhi->Fill(ptMC,(pt-ptMC)/(pt),phi);
 687       fPtITSminPtMCvsPtITSNPointITS->Fill(ptMC,(pt-ptMC)/(pt),nPointITS);
 688       fPtITSminPtMCvsPtITSNSigmaToVertex->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),nSigmaToVertex);
 689       fPtITSminPtMCvsPtITSChi2C->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),chi2C);
 690       fPtITSminPtMCvsPtITSRel1PtUncertainty->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),relUncertainty1Pt);
 691     }//fTrackCutsITS loop
 692
 693     if(fTrackType==1 || fTrackType==2) delete track;
 694
 695   }//ESD track loop
 696
 697   // Post output data
 698   PostData(0, fHistList);
 699   PostData(1, fHistListITS);
 700
 701 }
 702 //________________________________________________________________________
 703 Bool_t AliPWG4HighPtQAMC::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
 704   //
 705   // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
 706   // This is to called in Notify and should provide the path to the AOD/ESD file
 707   // Copied from AliAnalysisTaskJetSpectrum2
 708   //
 709
 710   TString file(currFile);
 711   fXsec = 0;
 712   fTrials = 1;
 713
 714   if(file.Contains("root_archive.zip#")){
 715     Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
 716     Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
 717     file.Replace(pos+1,20,"");
 718   }
 719   else {
 720     // not an archive take the basename....
 721     file.ReplaceAll(gSystem->BaseName(file.Data()),"");
 722   }
 723   //  Printf("%s",file.Data());
 724
 725
 726   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
 727   if(!fxsec){
 728     // next trial fetch the histgram file
 729     fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
 730     if(!fxsec){
 731         // not a severe condition but inciate that we have no information
 732       return kFALSE;
 733     }
 734     else{
 735       // find the tlist we want to be independtent of the name so use the Tkey
 736       TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
 737       if(!key){
 738         fxsec->Close();
 739         return kFALSE;
 740       }
 741       TList *list = dynamic_cast<TList*>(key->ReadObj());
 742       if(!list){
 743         fxsec->Close();
 744         return kFALSE;
 745       }
 746       fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
 747       fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
 748       fxsec->Close();
 749     }
 750   } // no tree pyxsec.root
 751   else {
 752     TTree *xtree = (TTree*)fxsec->Get("Xsection");
 753     if(!xtree){
 754       fxsec->Close();
 755       return kFALSE;
 756     }
 757     UInt_t   ntrials  = 0;
 758     Double_t  xsection  = 0;
 759     xtree->SetBranchAddress("xsection",&xsection);
 760     xtree->SetBranchAddress("ntrials",&ntrials);
 761     xtree->GetEntry(0);
 762     fTrials = ntrials;
 763     fXsec = xsection;
 764     fxsec->Close();
 765   }
 766   return kTRUE;
 767 }
 768 //________________________________________________________________________
 769 Bool_t AliPWG4HighPtQAMC::Notify()
 770 {
 771   //
 772   // Implemented Notify() to read the cross sections
 773   // and number of trials from pyxsec.root
 774   // Copied from AliAnalysisTaskJetSpectrum2
 775   //
 776
 777   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 778   Float_t xsection = 0;
 779   Float_t ftrials  = 1;
 780
 781   fAvgTrials = 1;
 782   if(tree){
 783     TFile *curfile = tree->GetCurrentFile();
 784     if (!curfile) {
 785       Error("Notify","No current file");
 786       return kFALSE;
 787     }
 788     if(!fh1Xsec||!fh1Trials){
 789       //      Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
 790       return kFALSE;
 791     }
 792     PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
 793     fh1Xsec->Fill("<#sigma>",xsection);
 794     // construct a poor man average trials
 795     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
 796     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
 797   }
 798   return kTRUE;
 799 }
 800
 801 //________________________________________________________________________
 802 AliGenPythiaEventHeader*  AliPWG4HighPtQAMC::GetPythiaEventHeader(AliMCEvent *mcEvent){
 803
 804   if(!mcEvent)return 0;
 805   AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
 806   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
 807   if(!pythiaGenHeader){
 808     // cocktail ??
 809     AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
 810
 811     if (!genCocktailHeader) {
 812       //      AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
 813       //      AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
 814       return 0;
 815     }
 816     TList* headerList = genCocktailHeader->GetHeaders();
 817     for (Int_t i=0; i<headerList->GetEntries(); i++) {
 818       pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
 819       if (pythiaGenHeader)
 820         break;
 821     }
 822     if(!pythiaGenHeader){
 823       AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
 824       return 0;
 825     }
 826   }
 827   return pythiaGenHeader;
 828
 829 }
 830
 831 //________________________________________________________________________
 832 void AliPWG4HighPtQAMC::Terminate(Option_t *)
 833 {
 834   // The Terminate() function is the last function to be called during
 835   // a query. It always runs on the client, it can be used to present
 836   // the results graphically or save the results to file.
 837
 838 }
 839
 840 #endif