[u/mrichter/AliRoot.git] / PWGCF / EBYE / NetParticle / AliAnalysisNetParticleEffCont.cxx

//-*- Mode: C++ -*-

#include "TMath.h"
#include "TAxis.h"

#include "AliESDEvent.h"
#include "AliESDInputHandler.h"
#include "AliStack.h"
#include "AliMCEvent.h"

#include "AliESDtrackCuts.h"

#include "AliAnalysisNetParticleEffCont.h"

using namespace std;

// Task for NetParticle checks
// Author: Jochen Thaeder <jochen@thaeder.de>

ClassImp(AliAnalysisNetParticleEffCont)

/*
 * ---------------------------------------------------------------------------------
 *                            Constructor / Destructor
 * ---------------------------------------------------------------------------------
 */

//________________________________________________________________________
AliAnalysisNetParticleEffCont::AliAnalysisNetParticleEffCont() :
  fHelper(NULL),
  fPdgCode(2212),

  fESD(NULL), 
  fESDTrackCuts(NULL),

  fCentralityBin(-1.),
  fNTracks(0),

  fStack(NULL),
  fMCEvent(NULL),

  fLabelsRec(NULL),

  fHnEff(NULL),
  fHnCont(NULL) {
  // Constructor   

  AliLog::SetClassDebugLevel("AliAnalysisNetParticleEffCont",10);
}

//________________________________________________________________________
AliAnalysisNetParticleEffCont::~AliAnalysisNetParticleEffCont() {
  // Destructor

  if (fLabelsRec[0])
    delete[] (fLabelsRec[0]);
  if (fLabelsRec[1])
    delete[] (fLabelsRec[1]);
  if (fLabelsRec)
    delete[] fLabelsRec;
}

/*
 * ---------------------------------------------------------------------------------
 *                                 Public Methods
 * ---------------------------------------------------------------------------------
 */

//________________________________________________________________________
void AliAnalysisNetParticleEffCont::Initialize(AliESDtrackCuts *cuts , AliAnalysisNetParticleHelper* helper) {

  // -- ESD track cuts
  // -------------------
  fESDTrackCuts     = cuts;

  // -- Get Helper
  // ---------------
  fHelper           = helper;

  // -- Get particle species / pdgCode
  // -------------------------
  fPdgCode          = AliPID::ParticleCode(fHelper->GetParticleSpecies());

  // -- MC Labels for efficiency
  // -----------------------------
  fLabelsRec        = new Int_t*[2];
  for (Int_t ii = 0; ii < 2 ; ++ii)
    fLabelsRec[ii] = NULL;

  // -- Create THnSparse Histograms
  // --------------------------------
  CreateHistograms();
  
  return;
}

/*
 * ---------------------------------------------------------------------------------
 *                            Setup/Reset Methods - private
 * ---------------------------------------------------------------------------------
 */

//________________________________________________________________________
Int_t AliAnalysisNetParticleEffCont::SetupEvent(AliESDInputHandler *esdHandler, AliMCEvent *mcEvent) {
  // -- Setup event

  // -- Reset of event
  // -------------------
  ResetEvent();

  // -- Setup of event
  // -------------------
  fESD           = esdHandler->GetEvent();
  fNTracks       = fESD->GetNumberOfTracks();

  fMCEvent       = mcEvent;
  fStack         = fMCEvent->Stack();

  fCentralityBin = fHelper->GetCentralityBin();

  // -- Create label arrays
  // ------------------------
  fLabelsRec[0] = new Int_t[fNTracks];
  if(!fLabelsRec[0]) {
    AliError("Cannot create fLabelsRec[0]");
    return -1;
  }

  fLabelsRec[1] = new Int_t[fNTracks];
  if(!fLabelsRec[1]) {
    AliError("Cannot create fLabelsRec[1] for TPC");
    return -1;
  }

  for(Int_t ii = 0; ii < fNTracks; ++ii) {
    fLabelsRec[0][ii] = 0;
    fLabelsRec[1][ii] = 0;
  }

  return 0;
}

//________________________________________________________________________
void AliAnalysisNetParticleEffCont::ResetEvent() {
  // -- Reset event
  
  for (Int_t ii = 0; ii < 2 ; ++ii) {
    if (fLabelsRec[ii])
      delete[] fLabelsRec[ii];
    fLabelsRec[ii] = NULL;
  }

  return;
}

//________________________________________________________________________
void AliAnalysisNetParticleEffCont::Process() {
  // -- Process event

  // -- Setup (clean, create and fill) MC labels
  // ---------------------------------------------
  FillMCLabels();

  // -- Fill  MC histograms for efficiency studies
  // -----------------------------------------------
  FillMCEffHist();

  return;
}      

/*
 * ---------------------------------------------------------------------------------
 *                                 Private Methods
 * ---------------------------------------------------------------------------------
 */

//________________________________________________________________________
void AliAnalysisNetParticleEffCont::CreateHistograms() {
  // -- Create histograms

  Double_t dCent[2] = {-0.5, 8.5};
  Int_t iCent       = 9;
  
  Double_t dEta[2]  = {-0.9, 0.9}; 
  Int_t iEta        = Int_t((dEta[1]-dEta[0]) / 0.075) +1 ; 

  Double_t dRap[2]  = {-0.5, 0.5}; 
  Int_t iRap        = Int_t((dRap[1]-dRap[0]) / 0.075) +1 ; 

  Double_t dPhi[2]  = {0.0, TMath::TwoPi()};
  Int_t iPhi        = 42;

  Double_t dPt[2]   = {0.1, 3.0};
  Int_t iPt         = Int_t((dPt[1]-dPt[0]) / 0.075);

  Double_t dSign[2] = {-1.5, 1.5};
  Int_t iSign       = 3;

  // ------------------------------------------------------------------
  // -- Create THnSparseF - Eff
  // ------------------------------------------------------------------

  //                           cent:   etaMC:     yMC:  phiMC:   ptMC:     sign:findable:recStatus:pidStatus:   etaRec:  phiRec:  ptRec:deltaEta: deltaPhi: deltaPt
  Int_t    binHnEff[15] = {   iCent,    iEta,    iRap,    iPhi,    iPt,    iSign,       2,      2  ,      2  ,    iEta,    iPhi,    iPt,    iEta, 2*iPhi+1, 2*iPt+1};
  Double_t minHnEff[15] = {dCent[0], dEta[0], dRap[0], dPhi[0], dPt[0], dSign[0],    -0.5,     -0.5,     -0.5, dEta[0], dPhi[0], dPt[0], dEta[0], -dPhi[1], -dPt[1]};
  Double_t maxHnEff[15] = {dCent[1], dEta[1], dRap[0], dPhi[1], dPt[1], dSign[1],     1.5,      1.5,      1.5, dEta[1], dPhi[1], dPt[1], dEta[1],  dPhi[1],  dPt[1]};

  fHnEff = new THnSparseF("fHnEff", "cent:etaMC:yMC:phiMC:ptMC:sign:findable:recStatus:pidStatus:etaRec:phiRec:ptRec:deltaEta:deltaPhi:deltaPt", 
			  15, binHnEff, minHnEff, maxHnEff);

  fHnEff->Sumw2();    

  fHnEff->GetAxis(0)->SetTitle("centrality");                   //  0-5|5-10|10-20|20-30|30-40|40-50|50-60|60-70|70-80|80-90 --> 10 bins
  fHnEff->GetAxis(1)->SetTitle("#eta_{MC}");                    //  eta  [-0.9, 0.9]
  fHnEff->GetAxis(2)->SetTitle("#it{y}_{MC}");                  //  rapidity  [-0.5, 0.5]
  fHnEff->GetAxis(3)->SetTitle("#varphi_{MC} (rad)");           //  phi  [ 0. , 2Pi]
  fHnEff->GetAxis(4)->SetTitle("#it{p}_{T,MC} (GeV/#it{c})");   //  pt   [ 0.1, 3.0]
  
  fHnEff->GetAxis(5)->SetTitle("sign");                         //  -1 | 0 | +1 
  fHnEff->GetAxis(6)->SetTitle("findable");                     //  0 not findable      |  1 findable
  fHnEff->GetAxis(7)->SetTitle("recStatus");                    //  0 not reconstructed |  1 reconstructed
  fHnEff->GetAxis(8)->SetTitle("recPid");                       //  0 not accepted      |  1 accepted

  fHnEff->GetAxis(9)->SetTitle("#eta_{Rec}");                   //  eta  [-0.9, 0.9]
  fHnEff->GetAxis(10)->SetTitle("#varphi_{Rec} (rad)");         //  phi  [ 0. , 2Pi]
  fHnEff->GetAxis(11)->SetTitle("#it{p}_{T,Rec} (GeV/#it{c})"); //  pt   [ 0.1, 3.0]

  fHnEff->GetAxis(12)->SetTitle("#eta_{MC}-#eta_{Rec}");                      //  eta  [-0.9,0.9]
  fHnEff->GetAxis(13)->SetTitle("#varphi_{MC}-#varphi_{Rec} (rad)");          //  phi  [ 0. ,2Pi]
  fHnEff->GetAxis(14)->SetTitle("#it{p}_{T,MC}-#it{p}_{T,Rec} (GeV/#it{c})"); //  pt   [ 3.0,3.0]

  fHelper->BinLogAxis(fHnEff,  4);
  fHelper->BinLogAxis(fHnEff, 11);

  /* 
     >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
     
     creation -> findable -> reconstructed -> pid_TPC+TOF
        (1)         (2)          (3)            (4)      
                                                                      ||   findable | recStatus | recPid 
     1) all primary probeParticles_MC                                 ||      -           -         -
     2) all findable primary probeParticles_MC                        ||      x           -         -
     3) all reconstructed primary probeParticles_MC                   ||      x           x         -
     4) all reconstructed primary probeParticles_MC & recPid_TPC+TOF  ||      x           x         x
     
     <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
  */ 

  // ------------------------------------------------------------------
  // -- Create THnSparseF - Cont
  // ------------------------------------------------------------------

  //                            cent:   etaMC:    yMC:   phiMC:   ptMC:     sign:contPart: contSign:  etaRec:  phiRec:  ptRec:deltaEta: deltaPhi: deltaPt
  Int_t    binHnCont[14] = {   iCent,    iEta,    iRap,    iPhi,    iPt,    iSign,       8,    iSign,    iEta,    iPhi,    iPt,    iEta, 2*iPhi+1, 2*iPt+1};
  Double_t minHnCont[14] = {dCent[0], dEta[0], dRap[0], dPhi[0], dPt[0], dSign[0],     0.5, dSign[0], dEta[0], dPhi[0], dPt[0], dEta[0], -dPhi[1], -dPt[1]};
  Double_t maxHnCont[14] = {dCent[1], dEta[1], dRap[1], dPhi[1], dPt[1], dSign[1],     8.5, dSign[1], dEta[1], dPhi[1], dPt[1], dEta[1],  dPhi[1],  dPt[1]};
  fHnCont = new THnSparseF("fHnCont", "cent:etaMC:yMC:phiMC:ptMC:sign:contPart:contSign:etaRec:phiRec:ptRec:deltaEta:deltaPhi:deltaPt",
			   14, binHnCont, minHnCont, maxHnCont);

  fHnCont->Sumw2();    
  
  fHnCont->GetAxis(0)->SetTitle("centrality");                   //  0-5|5-10|10-20|20-30|30-40|40-50|50-60|60-70|70-80|80-90 --> 10 bins
  fHnCont->GetAxis(1)->SetTitle("#eta_{MC}");                    //  eta  [-0.9,0.9]
  fHnCont->GetAxis(2)->SetTitle("#it{y}_{MC}");                  //  rapidity  [-0.5, 0.5]
  fHnCont->GetAxis(3)->SetTitle("#varphi_{MC} (rad)");           //  phi  [ 0. ,2Pi]
  fHnCont->GetAxis(4)->SetTitle("#it{p}_{T,MC} (GeV/#it{c})");   //  pT   [ 0.1,1.3]
  fHnCont->GetAxis(5)->SetTitle("sign");                         //  -1 | 0 | +1 
  
  fHnCont->GetAxis(6)->SetTitle("contPart");                     //  1 pi | 2 K | 3 p | 4 e | 5 mu | 6 other | 7 p from WeakDecay | 8 p from Material
  fHnCont->GetAxis(7)->SetTitle("contSign");                     //  -1 | 0 | +1 
   
  fHnCont->GetAxis(8)->SetTitle("#eta_{Rec}");                   //  eta  [-0.9, 0.9]
  fHnCont->GetAxis(9)->SetTitle("#varphi_{Rec} (rad)");          //  phi  [ 0. , 2Pi]
  fHnCont->GetAxis(10)->SetTitle("#it{p}_{T,Rec} (GeV/#it{c})"); //  pt   [ 0.1, 3.0]

  fHnCont->GetAxis(11)->SetTitle("#eta_{MC}-#eta_{Rec}");                      //  eta  [-0.9,0.9]
  fHnCont->GetAxis(12)->SetTitle("#varphi_{MC}-#varphi_{Rec} (rad)");          //  phi  [ 0. ,2Pi]
  fHnCont->GetAxis(13)->SetTitle("#it{p}_{T,MC}-#it{p}_{T,Rec} (GeV/#it{c})"); //  pt   [ 3.0,3.0]

  fHelper->BinLogAxis(fHnCont,  4);
  fHelper->BinLogAxis(fHnCont, 10);

  // ------------------------------------------------------------------
  
  return;
}

//________________________________________________________________________
void AliAnalysisNetParticleEffCont::FillMCLabels() {
  // Fill MC labels
  // Loop over ESD tracks and fill arrays with MC lables
  //  fLabelsRec[0] : all Tracks
  //  fLabelsRec[1] : all Tracks accepted by PID of TPC
  // Check every accepted track if correctly identified
  //  otherwise check for contamination

  for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
    AliESDtrack *track = fESD->GetTrack(idxTrack); 
    
    // -- Check if track is accepted for basic parameters
    if (!fHelper->IsTrackAcceptedBasicCharged(track))
      continue;
    
    // -- Check if accepted
    if (!fESDTrackCuts->AcceptTrack(track)) 
      continue;

    // -- Check if accepted in rapidity window
    Double_t yP;
    if (!fHelper->IsTrackAcceptedRapidity(track, yP))
      continue;

    // -- Check if accepted with thighter DCA cuts
    if (!fHelper->IsTrackAcceptedDCA(track))
      continue;

    Int_t label  = TMath::Abs(track->GetLabel()); 
    
    // -- Fill Label of all reconstructed
    fLabelsRec[0][idxTrack] = label;

    // -- Check if accepted by PID from TPC or TPC+TOF
    Double_t pid[2];
    if (!fHelper->IsTrackAcceptedPID(track, pid))
      continue;

    // -- Fill Label of all reconstructed && recPid_TPC+TOF    
    fLabelsRec[1][idxTrack] = label;    
    
    // -- Check for contamination and fill contamination THnSparse
    CheckContTrack(label, track->GetSign(), idxTrack);

  } // for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {

  return;
}

//________________________________________________________________________
void AliAnalysisNetParticleEffCont::CheckContTrack(Int_t label, Float_t sign, Int_t idxTrack) {
  // Check if particle is contamination or correctly identified
  // Fill contamination THnSparse

  TParticle* particle = fStack->Particle(label);
  if (!particle)
    return;

  Bool_t isSecondaryFromWeakDecay = kFALSE;
  Bool_t isSecondaryFromMaterial  = kFALSE;
  
  // -- Check if correctly identified 
  if (particle->GetPdgCode() == (sign*fPdgCode)) {

    // Check if is physical primary -> all ok 
    if (fStack->IsPhysicalPrimary(label))
      return;
    
    // -- Check if secondaries from material or weak decay
    isSecondaryFromWeakDecay = fStack->IsSecondaryFromWeakDecay(label);
    isSecondaryFromMaterial  = fStack->IsSecondaryFromMaterial(label);
  } 
  
  // -- Get MC pdg
  Float_t contSign = 0.;
  if      (particle->GetPDG()->Charge() == 0.) contSign =  0.;
  else if (particle->GetPDG()->Charge() <  0.) contSign = -1.;	
  else if (particle->GetPDG()->Charge() >  0.) contSign =  1.;	

  // -- Get contaminating particle
  Float_t contPart = 0;
  if        (isSecondaryFromWeakDecay)                   contPart = 7; // probeParticle from WeakDecay
  else if   (isSecondaryFromMaterial)                    contPart = 8; // probeParticle from Material
  else {
    if      (TMath::Abs(particle->GetPdgCode()) ==  211) contPart = 1; // pion
    else if (TMath::Abs(particle->GetPdgCode()) ==  321) contPart = 2; // kaon
    else if (TMath::Abs(particle->GetPdgCode()) == 2212) contPart = 3; // proton
    else if (TMath::Abs(particle->GetPdgCode()) ==   11) contPart = 4; // electron
    else if (TMath::Abs(particle->GetPdgCode()) ==   13) contPart = 5; // muon
    else                                                 contPart = 6; // other
  }
  
  // -- Get Reconstructed values 
  Float_t etaRec = 0.;
  Float_t phiRec = 0.;
  Float_t ptRec  = 0.;

  AliESDtrack *track = fESD->GetTrack(idxTrack);
  if (track) {
    // if no track present (which should not happen)
    // -> pt = 0. , which is not in the looked at range
    
    // -- Get Reconstructed eta/phi/pt
    etaRec = track->Eta();
    phiRec = track->Phi();	  
    ptRec  = track->Pt();
  } 	

  // -- Fill THnSparse
  Double_t hnCont[14] = {fCentralityBin, particle->Eta(), particle->Y(), particle->Phi(), particle->Pt(), sign, 
			 contPart, contSign, etaRec, phiRec, ptRec, 
			 particle->Eta()-etaRec, particle->Phi()-phiRec, particle->Pt()-ptRec};
  fHnCont->Fill(hnCont);
}

//________________________________________________________________________
void AliAnalysisNetParticleEffCont::FillMCEffHist() {
  // Fill efficiency THnSparse
  
  Int_t nPart  = fStack->GetNprimary();

  for (Int_t idxMC = 0; idxMC < nPart; ++idxMC) {
    TParticle* particle = fStack->Particle(idxMC);

    // -- Check basic MC properties -> charged physical primary
    if (!fHelper->IsParticleAcceptedBasicCharged(particle, idxMC))
      continue;

    // -- Check rapidity window
    Double_t yMC;
    if (!fHelper->IsParticleAcceptedRapidity(particle, yMC))
      continue;

    // -- Check if probeParticle / anti-probeParticle
    if (TMath::Abs(particle->GetPdgCode()) != fPdgCode)
      continue;
    
    // -- Get sign of particle
    Float_t sign      = (particle->GetPdgCode() < 0) ? -1. : 1.;

    // -- Get if particle is findable 
    Float_t findable  = Float_t(fHelper->IsParticleFindable(idxMC));

    // -- Get recStatus and pidStatus
    Float_t recStatus = 0.;
    Float_t recPid    = 0.;

    // -- Get Reconstructed values 
    Float_t etaRec = 0.;
    Float_t phiRec = 0.;
    Float_t ptRec  = 0.;

    // -- Loop over all labels
    for (Int_t idxRec=0; idxRec < fNTracks; ++idxRec) {
      if (idxMC == fLabelsRec[0][idxRec]) {
	recStatus = 1.;
	
	if (idxMC == fLabelsRec[1][idxRec])
	  recPid = 1.;

	AliESDtrack *track = fESD->GetTrack(idxRec);
	if (track) {
	  // if no track present (which should not happen)
	  // -> pt = 0. , which is not in the looked at range

	  // -- Get Reconstructed eta/phi/pt
	  etaRec = track->Eta();
	  phiRec = track->Phi();	  
	  ptRec  = track->Pt();
	} 	
	break;
      }
    } // for (Int_t idxRec=0; idxRec < fNTracks; ++idxRec) {  
    
    // -- Fill THnSparse
    Double_t hnEff[15] = {fCentralityBin, particle->Eta(), particle->Y(), particle->Phi(), particle->Pt(), sign, 
			  findable, recStatus, recPid, etaRec, phiRec, ptRec, 
			  particle->Eta()-etaRec, particle->Phi()-phiRec, particle->Pt()-ptRec};
    fHnEff->Fill(hnEff);

  } // for (Int_t idxMC = 0; idxMC < nPart; ++idxMC) {
  
  return;
}
Commit	Line	Data
cb68eb1d	1	//-- Mode: C++ --
	2
	3	#include "TMath.h"
	4	#include "TAxis.h"
	5
	6	#include "AliESDEvent.h"
	7	#include "AliESDInputHandler.h"
	8	#include "AliStack.h"
	9	#include "AliMCEvent.h"
	10
	11	#include "AliESDtrackCuts.h"
	12
	13	#include "AliAnalysisNetParticleEffCont.h"
	14
	15	using namespace std;
	16
	17	// Task for NetParticle checks
	18	// Author: Jochen Thaeder <jochen@thaeder.de>
	19
	20	ClassImp(AliAnalysisNetParticleEffCont)
	21
	22	/*
	23	* ---------------------------------------------------------------------------------
	24	* Constructor / Destructor
	25	* ---------------------------------------------------------------------------------
	26	*/
	27
	28	//________________________________________________________________________
	29	AliAnalysisNetParticleEffCont::AliAnalysisNetParticleEffCont() :
	30	fHelper(NULL),
	31	fPdgCode(2212),
	32
	33	fESD(NULL),
	34	fESDTrackCuts(NULL),
	35
	36	fCentralityBin(-1.),
	37	fNTracks(0),
	38
	39	fStack(NULL),
	40	fMCEvent(NULL),
	41
	42	fLabelsRec(NULL),
	43
	44	fHnEff(NULL),
	45	fHnCont(NULL) {
	46	// Constructor
	47
	48	AliLog::SetClassDebugLevel("AliAnalysisNetParticleEffCont",10);
	49	}
	50
	51	//________________________________________________________________________
	52	AliAnalysisNetParticleEffCont::~AliAnalysisNetParticleEffCont() {
	53	// Destructor
	54
	55	if (fLabelsRec[0])
	56	delete[] (fLabelsRec[0]);
	57	if (fLabelsRec[1])
	58	delete[] (fLabelsRec[1]);
	59	if (fLabelsRec)
	60	delete[] fLabelsRec;
	61	}
	62
	63	/*
	64	* ---------------------------------------------------------------------------------
65	* Public Methods
66	* ---------------------------------------------------------------------------------
67	*/
68
69	//________________________________________________________________________
70	void AliAnalysisNetParticleEffCont::Initialize(AliESDtrackCuts cuts , AliAnalysisNetParticleHelper helper) {
71
72	// -- ESD track cuts
73	// -------------------
74	fESDTrackCuts = cuts;
75
76	// -- Get Helper
77	// ---------------
78	fHelper = helper;
79
80	// -- Get particle species / pdgCode
81	// -------------------------
82	fPdgCode = AliPID::ParticleCode(fHelper->GetParticleSpecies());
83
84	// -- MC Labels for efficiency
85	// -----------------------------
86	fLabelsRec = new Int_t*[2];
87	for (Int_t ii = 0; ii < 2 ; ++ii)
88	fLabelsRec[ii] = NULL;
89
90	// -- Create THnSparse Histograms
91	// --------------------------------
92	CreateHistograms();
93
94	return;
95	}
96
97	/*
98	* ---------------------------------------------------------------------------------
99	* Setup/Reset Methods - private
100	* ---------------------------------------------------------------------------------
101	*/
102
103	//________________________________________________________________________
104	Int_t AliAnalysisNetParticleEffCont::SetupEvent(AliESDInputHandler esdHandler, AliMCEvent mcEvent) {
105	// -- Setup event
106
107	// -- Reset of event
108	// -------------------
109	ResetEvent();
110
111	// -- Setup of event
112	// -------------------
113	fESD = esdHandler->GetEvent();
114	fNTracks = fESD->GetNumberOfTracks();
115
116	fMCEvent = mcEvent;
117	fStack = fMCEvent->Stack();
118
119	fCentralityBin = fHelper->GetCentralityBin();
120
121	// -- Create label arrays
122	// ------------------------
123	fLabelsRec[0] = new Int_t[fNTracks];
124	if(!fLabelsRec[0]) {
125	AliError("Cannot create fLabelsRec[0]");
126	return -1;
127	}
128
129	fLabelsRec[1] = new Int_t[fNTracks];
130	if(!fLabelsRec[1]) {
131	AliError("Cannot create fLabelsRec[1] for TPC");
132	return -1;
133	}
134
135	for(Int_t ii = 0; ii < fNTracks; ++ii) {
136	fLabelsRec[0][ii] = 0;
137	fLabelsRec[1][ii] = 0;
138	}
139
140	return 0;
141	}
142
143	//________________________________________________________________________
144	void AliAnalysisNetParticleEffCont::ResetEvent() {
145	// -- Reset event
146
147	for (Int_t ii = 0; ii < 2 ; ++ii) {
148	if (fLabelsRec[ii])
149	delete[] fLabelsRec[ii];
150	fLabelsRec[ii] = NULL;
151	}
152
153	return;
154	}
155
156	//________________________________________________________________________
157	void AliAnalysisNetParticleEffCont::Process() {
158	// -- Process event
159
160	// -- Setup (clean, create and fill) MC labels
161	// ---------------------------------------------
162	FillMCLabels();
163
164	// -- Fill MC histograms for efficiency studies
165	// -----------------------------------------------
166	FillMCEffHist();
167
168	return;
169	}
170
171	/*
172	* ---------------------------------------------------------------------------------
173	* Private Methods
174	* ---------------------------------------------------------------------------------
175	*/
176
177	//________________________________________________________________________
178	void AliAnalysisNetParticleEffCont::CreateHistograms() {
179	// -- Create histograms
180
181	Double_t dCent[2] = {-0.5, 8.5};
182	Int_t iCent = 9;
183
184	Double_t dEta[2] = {-0.9, 0.9};
185	Int_t iEta = Int_t((dEta[1]-dEta[0]) / 0.075) +1 ;
186
187	Double_t dRap[2] = {-0.5, 0.5};
188	Int_t iRap = Int_t((dRap[1]-dRap[0]) / 0.075) +1 ;
189
190	Double_t dPhi[2] = {0.0, TMath::TwoPi()};
191	Int_t iPhi = 42;
192
193	Double_t dPt[2] = {0.1, 3.0};
194	Int_t iPt = Int_t((dPt[1]-dPt[0]) / 0.075);
195
196	Double_t dSign[2] = {-1.5, 1.5};
197	Int_t iSign = 3;
198
199	// ------------------------------------------------------------------
200	// -- Create THnSparseF - Eff
201	// ------------------------------------------------------------------
202
203	// cent: etaMC: yMC: phiMC: ptMC: sign:findable:recStatus:pidStatus: etaRec: phiRec: ptRec:deltaEta: deltaPhi: deltaPt
204	Int_t binHnEff[15] = { iCent, iEta, iRap, iPhi, iPt, iSign, 2, 2 , 2 , iEta, iPhi, iPt, iEta, 2iPhi+1, 2iPt+1};
205	Double_t minHnEff[15] = {dCent[0], dEta[0], dRap[0], dPhi[0], dPt[0], dSign[0], -0.5, -0.5, -0.5, dEta[0], dPhi[0], dPt[0], dEta[0], -dPhi[1], -dPt[1]};
206	Double_t maxHnEff[15] = {dCent[1], dEta[1], dRap[0], dPhi[1], dPt[1], dSign[1], 1.5, 1.5, 1.5, dEta[1], dPhi[1], dPt[1], dEta[1], dPhi[1], dPt[1]};
207
208	fHnEff = new THnSparseF("fHnEff", "cent:etaMC:yMC:phiMC:ptMC:sign:findable:recStatus:pidStatus:etaRec:phiRec:ptRec:deltaEta:deltaPhi:deltaPt",
209	15, binHnEff, minHnEff, maxHnEff);
210
211	fHnEff->Sumw2();
212
213	fHnEff->GetAxis(0)->SetTitle("centrality"); // 0-5\|5-10\|10-20\|20-30\|30-40\|40-50\|50-60\|60-70\|70-80\|80-90 --> 10 bins
214	fHnEff->GetAxis(1)->SetTitle("#eta_{MC}"); // eta [-0.9, 0.9]
215	fHnEff->GetAxis(2)->SetTitle("#it{y}_{MC}"); // rapidity [-0.5, 0.5]
216	fHnEff->GetAxis(3)->SetTitle("#varphi_{MC} (rad)"); // phi [ 0. , 2Pi]
217	fHnEff->GetAxis(4)->SetTitle("#it{p}_{T,MC} (GeV/#it{c})"); // pt [ 0.1, 3.0]
218
219	fHnEff->GetAxis(5)->SetTitle("sign"); // -1 \| 0 \| +1
220	fHnEff->GetAxis(6)->SetTitle("findable"); // 0 not findable \| 1 findable
221	fHnEff->GetAxis(7)->SetTitle("recStatus"); // 0 not reconstructed \| 1 reconstructed
222	fHnEff->GetAxis(8)->SetTitle("recPid"); // 0 not accepted \| 1 accepted
223
224	fHnEff->GetAxis(9)->SetTitle("#eta_{Rec}"); // eta [-0.9, 0.9]
225	fHnEff->GetAxis(10)->SetTitle("#varphi_{Rec} (rad)"); // phi [ 0. , 2Pi]
226	fHnEff->GetAxis(11)->SetTitle("#it{p}_{T,Rec} (GeV/#it{c})"); // pt [ 0.1, 3.0]
227
228	fHnEff->GetAxis(12)->SetTitle("#eta_{MC}-#eta_{Rec}"); // eta [-0.9,0.9]
229	fHnEff->GetAxis(13)->SetTitle("#varphi_{MC}-#varphi_{Rec} (rad)"); // phi [ 0. ,2Pi]
230	fHnEff->GetAxis(14)->SetTitle("#it{p}_{T,MC}-#it{p}_{T,Rec} (GeV/#it{c})"); // pt [ 3.0,3.0]
231
232	fHelper->BinLogAxis(fHnEff, 4);
233	fHelper->BinLogAxis(fHnEff, 11);
234
235	/*
236	>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
237
238	creation -> findable -> reconstructed -> pid_TPC+TOF
239	(1) (2) (3) (4)
240	\|\| findable \| recStatus \| recPid
241	1) all primary probeParticles_MC \|\| - - -
242	2) all findable primary probeParticles_MC \|\| x - -
243	3) all reconstructed primary probeParticles_MC \|\| x x -
244	4) all reconstructed primary probeParticles_MC & recPid_TPC+TOF \|\| x x x
245
246	<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
247	*/
248
249	// ------------------------------------------------------------------
250	// -- Create THnSparseF - Cont
251	// ------------------------------------------------------------------
252
253	// cent: etaMC: yMC: phiMC: ptMC: sign:contPart: contSign: etaRec: phiRec: ptRec:deltaEta: deltaPhi: deltaPt
254	Int_t binHnCont[14] = { iCent, iEta, iRap, iPhi, iPt, iSign, 8, iSign, iEta, iPhi, iPt, iEta, 2iPhi+1, 2iPt+1};
255	Double_t minHnCont[14] = {dCent[0], dEta[0], dRap[0], dPhi[0], dPt[0], dSign[0], 0.5, dSign[0], dEta[0], dPhi[0], dPt[0], dEta[0], -dPhi[1], -dPt[1]};
256	Double_t maxHnCont[14] = {dCent[1], dEta[1], dRap[1], dPhi[1], dPt[1], dSign[1], 8.5, dSign[1], dEta[1], dPhi[1], dPt[1], dEta[1], dPhi[1], dPt[1]};
257	fHnCont = new THnSparseF("fHnCont", "cent:etaMC:yMC:phiMC:ptMC:sign:contPart:contSign:etaRec:phiRec:ptRec:deltaEta:deltaPhi:deltaPt",
258	14, binHnCont, minHnCont, maxHnCont);
259
260	fHnCont->Sumw2();
261
262	fHnCont->GetAxis(0)->SetTitle("centrality"); // 0-5\|5-10\|10-20\|20-30\|30-40\|40-50\|50-60\|60-70\|70-80\|80-90 --> 10 bins
263	fHnCont->GetAxis(1)->SetTitle("#eta_{MC}"); // eta [-0.9,0.9]
264	fHnCont->GetAxis(2)->SetTitle("#it{y}_{MC}"); // rapidity [-0.5, 0.5]
265	fHnCont->GetAxis(3)->SetTitle("#varphi_{MC} (rad)"); // phi [ 0. ,2Pi]
266	fHnCont->GetAxis(4)->SetTitle("#it{p}_{T,MC} (GeV/#it{c})"); // pT [ 0.1,1.3]
267	fHnCont->GetAxis(5)->SetTitle("sign"); // -1 \| 0 \| +1
268
269	fHnCont->GetAxis(6)->SetTitle("contPart"); // 1 pi \| 2 K \| 3 p \| 4 e \| 5 mu \| 6 other \| 7 p from WeakDecay \| 8 p from Material
270	fHnCont->GetAxis(7)->SetTitle("contSign"); // -1 \| 0 \| +1
271
272	fHnCont->GetAxis(8)->SetTitle("#eta_{Rec}"); // eta [-0.9, 0.9]
273	fHnCont->GetAxis(9)->SetTitle("#varphi_{Rec} (rad)"); // phi [ 0. , 2Pi]
274	fHnCont->GetAxis(10)->SetTitle("#it{p}_{T,Rec} (GeV/#it{c})"); // pt [ 0.1, 3.0]
275
276	fHnCont->GetAxis(11)->SetTitle("#eta_{MC}-#eta_{Rec}"); // eta [-0.9,0.9]
277	fHnCont->GetAxis(12)->SetTitle("#varphi_{MC}-#varphi_{Rec} (rad)"); // phi [ 0. ,2Pi]
278	fHnCont->GetAxis(13)->SetTitle("#it{p}_{T,MC}-#it{p}_{T,Rec} (GeV/#it{c})"); // pt [ 3.0,3.0]
279
280	fHelper->BinLogAxis(fHnCont, 4);
281	fHelper->BinLogAxis(fHnCont, 10);
282
283	// ------------------------------------------------------------------
284
285	return;
286	}
287
288	//________________________________________________________________________
289	void AliAnalysisNetParticleEffCont::FillMCLabels() {
290	// Fill MC labels
291	// Loop over ESD tracks and fill arrays with MC lables
292	// fLabelsRec[0] : all Tracks
293	// fLabelsRec[1] : all Tracks accepted by PID of TPC
294	// Check every accepted track if correctly identified
295	// otherwise check for contamination
296
297	for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
298	AliESDtrack *track = fESD->GetTrack(idxTrack);
299
300	// -- Check if track is accepted for basic parameters
301	if (!fHelper->IsTrackAcceptedBasicCharged(track))
302	continue;
303
304	// -- Check if accepted
305	if (!fESDTrackCuts->AcceptTrack(track))
306	continue;
307
308	// -- Check if accepted in rapidity window
309	Double_t yP;
310	if (!fHelper->IsTrackAcceptedRapidity(track, yP))
311	continue;
312
313	// -- Check if accepted with thighter DCA cuts
314	if (!fHelper->IsTrackAcceptedDCA(track))
315	continue;
316
317	Int_t label = TMath::Abs(track->GetLabel());
318
319	// -- Fill Label of all reconstructed
320	fLabelsRec[0][idxTrack] = label;
321
322	// -- Check if accepted by PID from TPC or TPC+TOF
323	Double_t pid[2];
324	if (!fHelper->IsTrackAcceptedPID(track, pid))
325	continue;
326
327	// -- Fill Label of all reconstructed && recPid_TPC+TOF
328	fLabelsRec[1][idxTrack] = label;
329
330	// -- Check for contamination and fill contamination THnSparse
331	CheckContTrack(label, track->GetSign(), idxTrack);
332
333	} // for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
334
335	return;
336	}
337
338	//________________________________________________________________________
339	void AliAnalysisNetParticleEffCont::CheckContTrack(Int_t label, Float_t sign, Int_t idxTrack) {
340	// Check if particle is contamination or correctly identified
341	// Fill contamination THnSparse
342
343	TParticle* particle = fStack->Particle(label);
344	if (!particle)
345	return;
346
347	Bool_t isSecondaryFromWeakDecay = kFALSE;
348	Bool_t isSecondaryFromMaterial = kFALSE;
349
350	// -- Check if correctly identified
351	if (particle->GetPdgCode() == (sign*fPdgCode)) {
352
353	// Check if is physical primary -> all ok
354	if (fStack->IsPhysicalPrimary(label))
355	return;
356
357	// -- Check if secondaries from material or weak decay
358	isSecondaryFromWeakDecay = fStack->IsSecondaryFromWeakDecay(label);
359	isSecondaryFromMaterial = fStack->IsSecondaryFromMaterial(label);
360	}
361
362	// -- Get MC pdg
363	Float_t contSign = 0.;
364	if (particle->GetPDG()->Charge() == 0.) contSign = 0.;
365	else if (particle->GetPDG()->Charge() < 0.) contSign = -1.;
366	else if (particle->GetPDG()->Charge() > 0.) contSign = 1.;
367
368	// -- Get contaminating particle
369	Float_t contPart = 0;
370	if (isSecondaryFromWeakDecay) contPart = 7; // probeParticle from WeakDecay
371	else if (isSecondaryFromMaterial) contPart = 8; // probeParticle from Material
372	else {
373	if (TMath::Abs(particle->GetPdgCode()) == 211) contPart = 1; // pion
374	else if (TMath::Abs(particle->GetPdgCode()) == 321) contPart = 2; // kaon
375	else if (TMath::Abs(particle->GetPdgCode()) == 2212) contPart = 3; // proton
376	else if (TMath::Abs(particle->GetPdgCode()) == 11) contPart = 4; // electron
377	else if (TMath::Abs(particle->GetPdgCode()) == 13) contPart = 5; // muon
378	else contPart = 6; // other
379	}
380
381	// -- Get Reconstructed values
382	Float_t etaRec = 0.;
383	Float_t phiRec = 0.;
384	Float_t ptRec = 0.;
385
386	AliESDtrack *track = fESD->GetTrack(idxTrack);
387	if (track) {
388	// if no track present (which should not happen)
389	// -> pt = 0. , which is not in the looked at range
390
391	// -- Get Reconstructed eta/phi/pt
392	etaRec = track->Eta();
393	phiRec = track->Phi();
394	ptRec = track->Pt();
395	}
396
397	// -- Fill THnSparse
398	Double_t hnCont[14] = {fCentralityBin, particle->Eta(), particle->Y(), particle->Phi(), particle->Pt(), sign,
399	contPart, contSign, etaRec, phiRec, ptRec,
400	particle->Eta()-etaRec, particle->Phi()-phiRec, particle->Pt()-ptRec};
401	fHnCont->Fill(hnCont);
402	}
403
404	//________________________________________________________________________
405	void AliAnalysisNetParticleEffCont::FillMCEffHist() {
406	// Fill efficiency THnSparse
407
408	Int_t nPart = fStack->GetNprimary();
409
410	for (Int_t idxMC = 0; idxMC < nPart; ++idxMC) {
411	TParticle* particle = fStack->Particle(idxMC);
412
413	// -- Check basic MC properties -> charged physical primary
414	if (!fHelper->IsParticleAcceptedBasicCharged(particle, idxMC))
415	continue;
416
417	// -- Check rapidity window
418	Double_t yMC;
419	if (!fHelper->IsParticleAcceptedRapidity(particle, yMC))
420	continue;
421
422	// -- Check if probeParticle / anti-probeParticle
423	if (TMath::Abs(particle->GetPdgCode()) != fPdgCode)
424	continue;
425
426	// -- Get sign of particle
427	Float_t sign = (particle->GetPdgCode() < 0) ? -1. : 1.;
428
429	// -- Get if particle is findable
430	Float_t findable = Float_t(fHelper->IsParticleFindable(idxMC));
431
432	// -- Get recStatus and pidStatus
433	Float_t recStatus = 0.;
434	Float_t recPid = 0.;
435
436	// -- Get Reconstructed values
437	Float_t etaRec = 0.;
438	Float_t phiRec = 0.;
439	Float_t ptRec = 0.;
440
441	// -- Loop over all labels
442	for (Int_t idxRec=0; idxRec < fNTracks; ++idxRec) {
443	if (idxMC == fLabelsRec[0][idxRec]) {
444	recStatus = 1.;
445
446	if (idxMC == fLabelsRec[1][idxRec])
447	recPid = 1.;
448
449	AliESDtrack *track = fESD->GetTrack(idxRec);
450	if (track) {
451	// if no track present (which should not happen)
452	// -> pt = 0. , which is not in the looked at range
453
454	// -- Get Reconstructed eta/phi/pt
455	etaRec = track->Eta();
456	phiRec = track->Phi();
457	ptRec = track->Pt();
458	}
459	break;
460	}
461	} // for (Int_t idxRec=0; idxRec < fNTracks; ++idxRec) {
462
463	// -- Fill THnSparse
464	Double_t hnEff[15] = {fCentralityBin, particle->Eta(), particle->Y(), particle->Phi(), particle->Pt(), sign,
465	findable, recStatus, recPid, etaRec, phiRec, ptRec,
466	particle->Eta()-etaRec, particle->Phi()-phiRec, particle->Pt()-ptRec};
467	fHnEff->Fill(hnEff);
468
469	} // for (Int_t idxMC = 0; idxMC < nPart; ++idxMC) {
470
471	return;
472	}