[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliPWG4HighPtQAMC.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

//-----------------------------------------------------------------------
// This class compares the global reconstruction with the MC information
// Momentum resolution is stored as function of track cuts and pt.
// Output: Histograms for different set of cuts
//-----------------------------------------------------------------------
// Author : Marta Verweij - UU
//-----------------------------------------------------------------------

#ifndef ALIPWG4HighPtQAMC_CXX
#define ALIPWG4HighPtQAMC_CXX

#include "AliPWG4HighPtQAMC.h"

#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TList.h"
#include "TChain.h"
#include "TH3F.h"
#include "AliAnalysisManager.h"
#include "AliESDInputHandler.h"
#include "AliMCEvent.h"
#include "AliMCEventHandler.h"
#include "AliStack.h"
#include "AliESDtrack.h"
#include "AliESDtrackCuts.h"
#include "AliExternalTrackParam.h"
#include "AliLog.h"

using namespace std; //required for resolving the 'cout' symbol

ClassImp(AliPWG4HighPtQAMC)

AliPWG4HighPtQAMC::AliPWG4HighPtQAMC(): AliAnalysisTask("AliPWG4HighPtQAMC", ""), 
  fESD(0), 
  fTrackCuts(0), 
  fTrackCutsITS(0),
  fNEventAll(0),
  fNEventSel(0),
  fPtAll(0),  
  fPtSel(0),  
  fPtAllminPtMCvsPtAll(0),
  fPtAllminPtMCvsPtAllNPointTPC(0),
  fPtAllminPtMCvsPtAllDCAR(0),
  fPtAllminPtMCvsPtAllDCAZ(0),
  fPtAllminPtMCvsPtAllPhi(0),
  fPtAllminPtMCvsPtAllNPointITS(0),
  fPtAllminPtMCvsPtAllNSigmaToVertex(0),
  fPtAllminPtMCvsPtAllChi2C(0),
  fPtAllminPtMCvsPtAllRel1PtUncertainty(0),
  fPtAllMC(0),
  fPtSelMC(0),
  fPtSelMCITS(0),
  fHistList(0),
  fPtSelITS(0),
  fPtITSminPtMCvsPtITS(0),
  fPtITSminPtMCvsPtITSNPointTPC(0),
  fPtITSminPtMCvsPtITSDCAR(0),
  fPtITSminPtMCvsPtITSDCAZ(0),
  fPtITSminPtMCvsPtITSPhi(0),
  fPtITSminPtMCvsPtITSNPointITS(0),
  fPtITSminPtMCvsPtITSNSigmaToVertex(0),
  fPtITSminPtMCvsPtITSChi2C(0),
  fPtITSminPtMCvsPtITSRel1PtUncertainty(0),
  fHistListITS(0)
{

}
//________________________________________________________________________
AliPWG4HighPtQAMC::AliPWG4HighPtQAMC(const char *name): 
  AliAnalysisTask(name, ""), 
  fESD(0),
  fTrackCuts(),
  fTrackCutsITS(),
  fNEventAll(0),
  fNEventSel(0),
  fPtAll(0),
  fPtSel(0),
  fPtAllminPtMCvsPtAll(0),
  fPtAllminPtMCvsPtAllNPointTPC(0),
  fPtAllminPtMCvsPtAllDCAR(0),
  fPtAllminPtMCvsPtAllDCAZ(0),
  fPtAllminPtMCvsPtAllPhi(0),
  fPtAllminPtMCvsPtAllNPointITS(0),
  fPtAllminPtMCvsPtAllNSigmaToVertex(0),
  fPtAllminPtMCvsPtAllChi2C(0),
  fPtAllminPtMCvsPtAllRel1PtUncertainty(0),
  fPtAllMC(0),
  fPtSelMC(0),
  fPtSelMCITS(0),
  fHistList(0),
  fPtSelITS(0),
  fPtITSminPtMCvsPtITS(0),
  fPtITSminPtMCvsPtITSNPointTPC(0),
  fPtITSminPtMCvsPtITSDCAR(0),
  fPtITSminPtMCvsPtITSDCAZ(0),
  fPtITSminPtMCvsPtITSPhi(0),
  fPtITSminPtMCvsPtITSNPointITS(0),
  fPtITSminPtMCvsPtITSNSigmaToVertex(0),
  fPtITSminPtMCvsPtITSChi2C(0),
  fPtITSminPtMCvsPtITSRel1PtUncertainty(0),
  fHistListITS(0)
{
  //
  // Constructor. Initialization of Inputs and Outputs
  //
  AliDebug(2,Form("AliPWG4HighPtQAMC","Calling Constructor"));
  // Input slot #0 works with a TChain ESD
  DefineInput(0, TChain::Class());
  // Output slot #0 writes into a TList
  DefineOutput(0, TList::Class());
  // Output slot #1 writes into a TList
  DefineOutput(1, TList::Class());
  // Output slot #2 writes into a TList
  DefineOutput(2, TList::Class());
}

//________________________________________________________________________
void AliPWG4HighPtQAMC::ConnectInputData(Option_t *) 
{
  // Connect ESD here
  // Called once
  AliDebug(2,Form(">> AliPWG4HighPtSpectra::ConnectInputData \n"));
  TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
  if (!tree) {
    AliDebug(2,Form("ERROR: Could not read chain from input slot 0"));
  } else {
    
    AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
    
    if (!esdH) {
      AliDebug(2,Form("ERROR: Could not get ESDInputHandler"));
    } else
      fESD = esdH->GetEvent();
  }
}

//________________________________________________________________________
void AliPWG4HighPtQAMC::CreateOutputObjects() {
  //Create output objects
  AliDebug(2,Form(">> AliPWG4HighPtQATPConly::CreateOutputObjects \n"));

  Bool_t oldStatus = TH1::AddDirectoryStatus();
  TH1::AddDirectory(kFALSE); 

  OpenFile(0);
  fHistList = new TList();
  OpenFile(1);
  fHistListITS = new TList();

  Int_t fgkNPhiBins=18;
  Float_t kMinPhi = 0.;
  Float_t kMaxPhi = 2.*TMath::Pi();
  
  Int_t fgkNPtBins=98;
  Float_t fgkPtMin=2.;
  Float_t fgkPtMax=100.;
  Int_t fgkResPtBins=80;

  fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
  fHistList->Add(fNEventAll);
  fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
  fHistList->Add(fNEventSel);
  fPtAll = new TH1F("fPtAll","PtAll",fgkNPtBins, fgkPtMin, fgkPtMax);
  fHistList->Add(fPtAll);
  fPtSel = new TH1F("fPtSel","PtSel",fgkNPtBins, fgkPtMin, fgkPtMax);
  fHistList->Add(fPtSel);
  
  fPtAllminPtMCvsPtAll = new TH2F("fPtAllminPtMCvsPtAll","PtAllminPtMCvsPtAll",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.);
  fPtAllminPtMCvsPtAll->SetXTitle("p_{t}^{MC}");
  fPtAllminPtMCvsPtAll->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
  fHistList->Add(fPtAllminPtMCvsPtAll);
  
  fPtAllminPtMCvsPtAllNPointTPC = new TH3F("fPtAllminPtMCvsPtAllNPointTPC","PtAllminPtMCvsPtAllNPointTPC",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,160,0.5,160.5);
  fPtAllminPtMCvsPtAllNPointTPC->SetXTitle("p_{t}^{MC}");
  fPtAllminPtMCvsPtAllNPointTPC->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
  fPtAllminPtMCvsPtAllNPointTPC->SetZTitle("N_{point,TPC}");
  fHistList->Add(fPtAllminPtMCvsPtAllNPointTPC);

  fPtAllminPtMCvsPtAllDCAR = new TH3F("fPtAllminPtMCvsPtAllDCAR","PtAllminPtMCvsPtAllDCAR",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,80,-0.2,0.2);
  fPtAllminPtMCvsPtAllDCAR->SetXTitle("p_{t}^{MC}");
  fPtAllminPtMCvsPtAllDCAR->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
  fPtAllminPtMCvsPtAllDCAR->SetZTitle("DCA_{R}");
  fHistList->Add(fPtAllminPtMCvsPtAllDCAR);

  fPtAllminPtMCvsPtAllDCAZ = new TH3F("fPtAllminPtMCvsPtAllDCAZ","PtAllminPtMCvsPtAllDCAZ",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,80,-2.,2.);
  fPtAllminPtMCvsPtAllDCAZ->SetXTitle("p_{t}^{MC}");
  fPtAllminPtMCvsPtAllDCAZ->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
  fPtAllminPtMCvsPtAllDCAZ->SetZTitle("DCA_{Z}");
  fHistList->Add(fPtAllminPtMCvsPtAllDCAZ);

  fPtAllminPtMCvsPtAllPhi = new TH3F("fPtAllminPtMCvsPtAllPhi","PtAllminPtMCvsPtAllPhi",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,fgkNPhiBins,kMinPhi,kMaxPhi);
  fPtAllminPtMCvsPtAllPhi->SetXTitle("p_{t}^{MC}");
  fPtAllminPtMCvsPtAllPhi->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
  fPtAllminPtMCvsPtAllPhi->SetZTitle("#phi");
  fHistList->Add(fPtAllminPtMCvsPtAllPhi);

  fPtAllminPtMCvsPtAllNPointITS = new TH3F("fPtAllminPtMCvsPtAllNPointITS","PtAllminPtMCvsPtAllNPointITS",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,9,-0.5,8.5);
  fPtAllminPtMCvsPtAllNPointITS->SetXTitle("p_{t}^{MC}");
  fPtAllminPtMCvsPtAllNPointITS->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
  fPtAllminPtMCvsPtAllNPointITS->SetZTitle("N_{point,ITS}}");
  fHistList->Add(fPtAllminPtMCvsPtAllNPointITS);
  
  fPtAllminPtMCvsPtAllNSigmaToVertex = new TH3F("fPtAllminPtMCvsPtAllNSigmaToVertex","PtAllminPtMCvsPtAllNSigmaToVertex",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,40,0.,8.);
  fPtAllminPtMCvsPtAllNSigmaToVertex->SetXTitle("p_{t}^{MC}");
  fPtAllminPtMCvsPtAllNSigmaToVertex->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
  fPtAllminPtMCvsPtAllNSigmaToVertex->SetZTitle("N#sigma to vertex");
  fHistList->Add(fPtAllminPtMCvsPtAllNSigmaToVertex);

  fPtAllminPtMCvsPtAllChi2C = new TH3F("fPtAllminPtMCvsPtAllChi2C","PtAllminPtMCvsPtAllChi2C",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,20,0.,10.);
  fPtAllminPtMCvsPtAllChi2C->SetXTitle("p_{t}^{MC}");
  fPtAllminPtMCvsPtAllChi2C->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
  fPtAllminPtMCvsPtAllChi2C->SetZTitle("Constrained #chi^{2}");
  fHistList->Add(fPtAllminPtMCvsPtAllChi2C);

  fPtAllminPtMCvsPtAllRel1PtUncertainty = new TH3F("fPtAllminPtMCvsPtAllRel1PtUncertainty","PtAllminPtMCvsPtAllRel1PtUncertainty",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,30,0.,0.3);
  fPtAllminPtMCvsPtAllRel1PtUncertainty->SetXTitle("p_{t}^{MC}");
  fPtAllminPtMCvsPtAllRel1PtUncertainty->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
  fPtAllminPtMCvsPtAllRel1PtUncertainty->SetZTitle("Rel1PtUncertainty");
  fHistList->Add(fPtAllminPtMCvsPtAllRel1PtUncertainty);

  //ITSrefit
  fPtSelITS = new TH1F("fPtSelITSrefit","PtSel",fgkNPtBins, fgkPtMin, fgkPtMax);
  fHistListITS->Add(fPtSelITS);
  
  fPtITSminPtMCvsPtITS = new TH2F("fPtITSminPtMCvsPtITS","PtITSminPtMCvsPtITS",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.);
  fPtITSminPtMCvsPtITS->SetXTitle("p_{t}^{MC}");
  fPtITSminPtMCvsPtITS->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
  fHistListITS->Add(fPtITSminPtMCvsPtITS);
  
  fPtITSminPtMCvsPtITSNPointTPC = new TH3F("fPtITSminPtMCvsPtITSNPointTPC","PtITSminPtMCvsPtITSNPointTPC",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,160,0.5,160.5);
  fPtITSminPtMCvsPtITSNPointTPC->SetXTitle("p_{t}^{MC}");
  fPtITSminPtMCvsPtITSNPointTPC->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
  fPtITSminPtMCvsPtITSNPointTPC->SetZTitle("N_{point,TPC}");
  fHistListITS->Add(fPtITSminPtMCvsPtITSNPointTPC);
    
  fPtITSminPtMCvsPtITSDCAR = new TH3F("fPtITSminPtMCvsPtITSDCAR","PtITSminPtMCvsPtITSDCAR",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,80,-0.2,0.2);
  fPtITSminPtMCvsPtITSDCAR->SetXTitle("p_{t}^{MC}");
  fPtITSminPtMCvsPtITSDCAR->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
  fPtITSminPtMCvsPtITSDCAR->SetZTitle("DCA_{R}");
  fHistListITS->Add(fPtITSminPtMCvsPtITSDCAR);
  
  fPtITSminPtMCvsPtITSDCAZ = new TH3F("fPtITSminPtMCvsPtITSDCAZ","PtITSminPtMCvsPtITSDCAZ",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,80,-2.,2.);
  fPtITSminPtMCvsPtITSDCAZ->SetXTitle("p_{t}^{MC}");
  fPtITSminPtMCvsPtITSDCAZ->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
  fPtITSminPtMCvsPtITSDCAZ->SetZTitle("DCA_{Z}");
  fHistListITS->Add(fPtITSminPtMCvsPtITSDCAZ);
  
  fPtITSminPtMCvsPtITSPhi = new TH3F("fPtITSminPtMCvsPtITSPhi","PtITSminPtMCvsPtITSPhi",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,fgkNPhiBins,kMinPhi,kMaxPhi);
  fPtITSminPtMCvsPtITSPhi->SetXTitle("p_{t}^{MC}");
  fPtITSminPtMCvsPtITSPhi->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
  fPtITSminPtMCvsPtITSPhi->SetZTitle("#phi");
  fHistListITS->Add(fPtITSminPtMCvsPtITSPhi);
  
  fPtITSminPtMCvsPtITSNPointITS = new TH3F("fPtITSminPtMCvsPtITSNPointITS","PtITSminPtMCvsPtITSNPointITS",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,9,-0.5,8.5);
  fPtITSminPtMCvsPtITSNPointITS->SetXTitle("p_{t}^{MC}");
  fPtITSminPtMCvsPtITSNPointITS->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})}");
  fPtITSminPtMCvsPtITSNPointITS->SetZTitle("N_{point,ITS}}");
  fHistListITS->Add(fPtITSminPtMCvsPtITSNPointITS); 
  
  fPtITSminPtMCvsPtITSNSigmaToVertex = new TH3F("fPtITSminPtMCvsPtITSNSigmaToVertex","PtITSminPtMCvsPtITSNSigmaToVertex",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,40,0.,8.);
  fPtITSminPtMCvsPtITSNSigmaToVertex->SetXTitle("p_{t}^{MC}");
  fPtITSminPtMCvsPtITSNSigmaToVertex->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
  fPtITSminPtMCvsPtITSNSigmaToVertex->SetZTitle("N#sigma to vertex");
  fHistListITS->Add(fPtITSminPtMCvsPtITSNSigmaToVertex);

  fPtITSminPtMCvsPtITSChi2C = new TH3F("fPtITSminPtMCvsPtITSChi2C","PtITSminPtMCvsPtITSChi2C",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,20,0.,10.);
  fPtITSminPtMCvsPtITSChi2C->SetXTitle("p_{t}^{MC}");
  fPtITSminPtMCvsPtITSChi2C->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
  fPtITSminPtMCvsPtITSChi2C->SetZTitle("Constrained #chi^{2}");
  fHistListITS->Add(fPtITSminPtMCvsPtITSChi2C);

  fPtITSminPtMCvsPtITSRel1PtUncertainty = new TH3F("fPtITSminPtMCvsPtITSRel1PtUncertainty","PtITSminPtMCvsPtITSRel1PtUncertainty",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,30,0.,0.3);
  fPtITSminPtMCvsPtITSRel1PtUncertainty->SetXTitle("p_{t}^{MC}");
  fPtITSminPtMCvsPtITSRel1PtUncertainty->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
  fPtITSminPtMCvsPtITSRel1PtUncertainty->SetZTitle("Rel1PtUncertainty");
  fHistListITS->Add(fPtITSminPtMCvsPtITSRel1PtUncertainty);

  fPtAllMC = new TH1F("fPtAllMC","PtAll",fgkNPtBins, fgkPtMin, fgkPtMax);
  fHistList->Add(fPtAllMC);
  fPtSelMC = new TH1F("fPtSelMC","PtSel",fgkNPtBins, fgkPtMin, fgkPtMax);
  fHistList->Add(fPtSelMC);
  fPtSelMCITS = new TH1F("fPtSelMCITS","PtSel",fgkNPtBins, fgkPtMin, fgkPtMax);
  fHistList->Add(fPtSelMCITS);
  
  TH1::AddDirectory(oldStatus); 

}
//________________________________________________________________________
void AliPWG4HighPtQAMC::Exec(Option_t *) {  
  // Main loop
  // Called for each event
  AliDebug(2,Form(">> AliPWG4HighPtQATPConly::Exec \n"));  
  
  // All events without selection
  fNEventAll->Fill(0.);

  if (!fESD) {
    AliDebug(2,Form("ERROR: fESD not available"));
    PostData(0, fHistList);
    PostData(1, fHistListITS);
    return;
  }

  Bool_t isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
  if(!isSelected) { //Select collison candidates
    AliDebug(2,Form(" Trigger Selection: event REJECTED ... "));
    // Post output data
    PostData(0, fHistList);
    PostData(1, fHistListITS);
    return;
  }
  
 AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
  if (!eventHandler) {
    AliDebug(2,Form("ERROR: Could not retrieve MC event handler"));
    PostData(0, fHistList);
    PostData(1, fHistListITS);
    return;
  }

  AliMCEvent* mcEvent = eventHandler->MCEvent();
  if (!mcEvent) {
    AliDebug(2,Form("ERROR: Could not retrieve MC event"));
    PostData(0, fHistList);
    PostData(1, fHistListITS);
    return;
  }

  AliDebug(2,Form("MC particles: %d", mcEvent->GetNumberOfTracks()));

  if (!fESD) {
    AliDebug(2,Form("ERROR: fESD not available"));
    PostData(0, fHistList);
    PostData(1, fHistListITS);
    return;
  }

  AliStack* stack = mcEvent->Stack();                //Particles Stack

  AliDebug(2,Form("MC particles stack: %d", stack->GetNtrack()));

  const AliESDVertex *vtx = fESD->GetPrimaryVertex();
  // Need vertex cut
  if (vtx->GetNContributors() < 2) {
    PostData(0, fHistList);
    PostData(1, fHistListITS);
    return;
  }
 double primVtx[3];
  vtx->GetXYZ(primVtx);
  //  printf("primVtx: %g  %g  %g \n",primVtx[0],primVtx[1],primVtx[2]);
  if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
    // Post output data
    PostData(0, fHistList);
    PostData(1, fHistListITS);
    return;
  }
  
  AliDebug(2,Form("Vertex title %s, status %d, nCont %d\n",vtx->GetTitle(), vtx->GetStatus(), vtx->GetNContributors()));

  // Need to keep track of evts without vertex
  fNEventSel->Fill(0.);

  if(!fESD->GetNumberOfTracks() || fESD->GetNumberOfTracks()<2)  {
    PostData(0, fHistList);
    PostData(1, fHistListITS);
    return;
  }

  Int_t nTracks = fESD->GetNumberOfTracks();
  AliDebug(2,Form("nTracks %d", nTracks));

  for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
    
    AliESDtrack *track = fESD->GetTrack(iTrack);
    if(!track) continue;
    Int_t label = TMath::Abs(track->GetLabel());
    TParticle *particle = stack->Particle(label) ;
    if(!particle) continue;

    Float_t pt = track->Pt();
    Float_t ptMC = particle->Pt();
    Float_t phi = track->Phi();
    Float_t dca2D, dcaZ;
    track->GetImpactParameters(dca2D,dcaZ);
    UChar_t itsMap = track->GetITSClusterMap();
    Int_t nPointITS = 0;
    for (Int_t i=0; i < 6; i++) {
      if (itsMap & (1 << i))
	nPointITS ++;
    }
    Float_t nSigmaToVertex = fTrackCuts->GetSigmaToVertex(track);// Calculates the number of sigma to the vertex for a track.
    Float_t chi2C = track->GetConstrainedChi2();
    Float_t relUncertainty1Pt = TMath::Sqrt(track->GetSigma1Pt2())*pt;

    fPtAll->Fill(pt);
    fPtAllMC->Fill(ptMC);
    
    if (fTrackCuts->AcceptTrack(track)) {

      fPtSel->Fill(pt);
      
      fPtSelMC->Fill(ptMC);
      fPtAllminPtMCvsPtAll->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC) );
      fPtAllminPtMCvsPtAllNPointTPC->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),track->GetTPCNcls());
      fPtAllminPtMCvsPtAllDCAR->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dca2D);
      fPtAllminPtMCvsPtAllDCAZ->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dcaZ);
      fPtAllminPtMCvsPtAllPhi->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),phi);
      fPtAllminPtMCvsPtAllNPointITS->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),nPointITS);
      fPtAllminPtMCvsPtAllNSigmaToVertex->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),nSigmaToVertex);
      fPtAllminPtMCvsPtAllChi2C->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),chi2C);
      fPtAllminPtMCvsPtAllRel1PtUncertainty->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),relUncertainty1Pt);
    }//fTrackCuts selection
    
    
    //ITSrefit selection
    if (fTrackCutsITS->AcceptTrack(track)) {
      
      fPtSelITS->Fill(pt);
      fPtSelMCITS->Fill(ptMC);
      fPtITSminPtMCvsPtITS->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC) );
      fPtITSminPtMCvsPtITSNPointTPC->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),track->GetTPCNcls());
      fPtITSminPtMCvsPtITSDCAR->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dca2D);
      fPtITSminPtMCvsPtITSDCAZ->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dcaZ);
      fPtITSminPtMCvsPtITSPhi->Fill(ptMC,(pt-ptMC)/(pt),phi);
      fPtITSminPtMCvsPtITSNPointITS->Fill(ptMC,(pt-ptMC)/(pt),nPointITS);
      fPtITSminPtMCvsPtITSNSigmaToVertex->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),nSigmaToVertex);
      fPtITSminPtMCvsPtITSChi2C->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),chi2C);
      fPtITSminPtMCvsPtITSRel1PtUncertainty->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),relUncertainty1Pt);
    }//fTrackCutsITS loop
      
  }//ESD track loop
   
  // Post output data
  PostData(0, fHistList);
  PostData(1, fHistListITS);

}
//________________________________________________________________________
void AliPWG4HighPtQAMC::Terminate(Option_t *)
{
  // The Terminate() function is the last function to be called during
  // a query. It always runs on the client, it can be used to present
  // the results graphically or save the results to file.

}

#endif
Commit	Line	Data
fdceab34	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
df943115	24	#ifndef ALIPWG4HighPtQAMC_CXX
	25	#define ALIPWG4HighPtQAMC_CXX
	26
fdceab34	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"
df943115	43	#include "AliLog.h"
fdceab34	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),
8f0faa80	53	fNEventAll(0),
8f0faa80	54	fNEventSel(0),
df943115	55	fPtAll(0),
df943115	56	fPtSel(0),
fdceab34	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	{
df943115	82
fdceab34	83	}
	84	//________________________________________________________________________
	85	AliPWG4HighPtQAMC::AliPWG4HighPtQAMC(const char *name):
	86	AliAnalysisTask(name, ""),
	87	fESD(0),
	88	fTrackCuts(),
	89	fTrackCutsITS(),
8f0faa80	90	fNEventAll(0),
8f0faa80	91	fNEventSel(0),
df943115	92	fPtAll(0),
df943115	93	fPtSel(0),
fdceab34	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	//
df943115	122	AliDebug(2,Form("AliPWG4HighPtQAMC","Calling Constructor"));
fdceab34	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());
df943115	131	}
fdceab34	132
	133	//________________________________________________________________________
	134	void AliPWG4HighPtQAMC::ConnectInputData(Option_t *)
	135	{
	136	// Connect ESD here
	137	// Called once
df943115	138	AliDebug(2,Form(">> AliPWG4HighPtSpectra::ConnectInputData \n"));
fdceab34	139	TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
fdceab34	140	if (!tree) {
df943115	141	AliDebug(2,Form("ERROR: Could not read chain from input slot 0"));
fdceab34	142	} else {
	143
	144	AliESDInputHandler esdH = dynamic_cast<AliESDInputHandler> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
	145
	146	if (!esdH) {
df943115	147	AliDebug(2,Form("ERROR: Could not get ESDInputHandler"));
fdceab34	148	} else
	149	fESD = esdH->GetEvent();
	150	}
	151	}
	152
fdceab34	153	//________________________________________________________________________
	154	void AliPWG4HighPtQAMC::CreateOutputObjects() {
	155	//Create output objects
df943115	156	AliDebug(2,Form(">> AliPWG4HighPtQATPConly::CreateOutputObjects \n"));
	157
	158	Bool_t oldStatus = TH1::AddDirectoryStatus();
	159	TH1::AddDirectory(kFALSE);
	160
fdceab34	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.;
8f0faa80	173	Int_t fgkResPtBins=80;
fdceab34	174
8f0faa80	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);
fdceab34	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
8f0faa80	195	fPtAllminPtMCvsPtAllDCAR = new TH3F("fPtAllminPtMCvsPtAllDCAR","PtAllminPtMCvsPtAllDCAR",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,80,-0.2,0.2);
fdceab34	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
8f0faa80	201	fPtAllminPtMCvsPtAllDCAZ = new TH3F("fPtAllminPtMCvsPtAllDCAZ","PtAllminPtMCvsPtAllDCAZ",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,80,-2.,2.);
fdceab34	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
8f0faa80	252	fPtITSminPtMCvsPtITSDCAR = new TH3F("fPtITSminPtMCvsPtITSDCAR","PtITSminPtMCvsPtITSDCAR",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,80,-0.2,0.2);
fdceab34	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
8f0faa80	258	fPtITSminPtMCvsPtITSDCAZ = new TH3F("fPtITSminPtMCvsPtITSDCAZ","PtITSminPtMCvsPtITSDCAZ",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,80,-2.,2.);
fdceab34	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
df943115	301	TH1::AddDirectory(oldStatus);
df943115	302
fdceab34	303	}
	304	//________________________________________________________________________
	305	void AliPWG4HighPtQAMC::Exec(Option_t *) {
	306	// Main loop
	307	// Called for each event
df943115	308	AliDebug(2,Form(">> AliPWG4HighPtQATPConly::Exec \n"));
fdceab34	309
8f0faa80	310	// All events without selection
	311	fNEventAll->Fill(0.);
	312
fdceab34	313	if (!fESD) {
df943115	314	AliDebug(2,Form("ERROR: fESD not available"));
8f0faa80	315	PostData(0, fHistList);
8f0faa80	316	PostData(1, fHistListITS);
fdceab34	317	return;
	318	}
	319
8f0faa80	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
fdceab34	329	AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
fdceab34	330	if (!eventHandler) {
df943115	331	AliDebug(2,Form("ERROR: Could not retrieve MC event handler"));
8f0faa80	332	PostData(0, fHistList);
8f0faa80	333	PostData(1, fHistListITS);
fdceab34	334	return;
	335	}
	336
	337	AliMCEvent* mcEvent = eventHandler->MCEvent();
	338	if (!mcEvent) {
df943115	339	AliDebug(2,Form("ERROR: Could not retrieve MC event"));
8f0faa80	340	PostData(0, fHistList);
8f0faa80	341	PostData(1, fHistListITS);
fdceab34	342	return;
	343	}
	344
df943115	345	AliDebug(2,Form("MC particles: %d", mcEvent->GetNumberOfTracks()));
fdceab34	346
fdceab34	347	if (!fESD) {
df943115	348	AliDebug(2,Form("ERROR: fESD not available"));
8f0faa80	349	PostData(0, fHistList);
8f0faa80	350	PostData(1, fHistListITS);
fdceab34	351	return;
	352	}
	353
	354	AliStack* stack = mcEvent->Stack(); //Particles Stack
	355
df943115	356	AliDebug(2,Form("MC particles stack: %d", stack->GetNtrack()));
fdceab34	357
fdceab34	358	const AliESDVertex *vtx = fESD->GetPrimaryVertex();
fdceab34	359	// Need vertex cut
8f0faa80	360	if (vtx->GetNContributors() < 2) {
	361	PostData(0, fHistList);
	362	PostData(1, fHistListITS);
fdceab34	363	return;
8f0faa80	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	}
df943115	374
df943115	375	AliDebug(2,Form("Vertex title %s, status %d, nCont %d\n",vtx->GetTitle(), vtx->GetStatus(), vtx->GetNContributors()));
fdceab34	376
fdceab34	377	// Need to keep track of evts without vertex
8f0faa80	378	fNEventSel->Fill(0.);
	379
	380	if(!fESD->GetNumberOfTracks() \|\| fESD->GetNumberOfTracks()<2) {
	381	PostData(0, fHistList);
	382	PostData(1, fHistListITS);
	383	return;
	384	}
fdceab34	385
fdceab34	386	Int_t nTracks = fESD->GetNumberOfTracks();
df943115	387	AliDebug(2,Form("nTracks %d", nTracks));
df943115	388
fdceab34	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	{
df943115	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.
fdceab34	461
fdceab34	462	}
df943115	463
df943115	464	#endif