[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),
  fNEvent(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(),
  fNEvent(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=40;

  fNEvent = new TH1F("fNEvent","NEvent",1,-0.5,0.5);
  fHistList->Add(fNEvent);
  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.,20,-1.,1.);
  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.,40,-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.,20,-1.,1.);
  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.,40,-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"));  
  
  if (!fESD) {
    AliDebug(2,Form("ERROR: fESD not available"));
    return;
  }

 AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
  if (!eventHandler) {
    AliDebug(2,Form("ERROR: Could not retrieve MC event handler"));
    return;
  }

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

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

  if (!fESD) {
    AliDebug(2,Form("ERROR: fESD not available"));
    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)
    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
  fNEvent->Fill(0.);

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