[u/mrichter/AliRoot.git] / PWGCF / EBYE / PIDFluctuation / task / AliEbyEPidRatioEffContExtra.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: ALICE Offline.                                                 *
 * 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.                  *
 **************************************************************************/


//=========================================================================//
//             AliEbyE Analysis for Particle Ratio Fluctuation             //
//                   Deepika Rathee  | Satyajit Jena                       //
//                   drathee@cern.ch | sjena@cern.ch                       //
//                  Date: Wed Jul  9 18:38:30 CEST 2014                    //
//          New approch to find particle ratio to reduce memory            //
//                             (Test Only)                                 //
//        Copied from NetParticle Classes
//        Origin: Jochen Thaeder <jochen@thaeder.de>
//                Michael Weber <m.weber@cern.ch>
//=========================================================================//

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

#include "AliESDEvent.h"
#include "AliStack.h"
#include "AliMCEvent.h"
#include "AliESDtrackCuts.h"
#include "AliAODEvent.h"
#include "AliAODMCParticle.h"

#include "AliEbyEPidRatioEffContExtra.h"

using namespace std;

ClassImp(AliEbyEPidRatioEffContExtra)

//________________________________________________________________________
AliEbyEPidRatioEffContExtra::AliEbyEPidRatioEffContExtra() :
  AliEbyEPidRatioBase("EffCont", "EffCont"),
  fLabelsRec(NULL),
 
  fHnNchEMc(NULL),
  fHnNchERec(NULL),
  fHnNpiEMc(NULL),
  fHnNpiERec(NULL),
  fHnNkaEMc(NULL),
  fHnNkaERec(NULL),
  fHnNprEMc(NULL),
  fHnNprERec(NULL),

  fHnNchCMc(NULL),
  fHnNchCRec(NULL),
  fHnNpiCMc(NULL),
  fHnNpiCRec(NULL),
  fHnNkaCMc(NULL),
  fHnNkaCRec(NULL),
  fHnNprCMc(NULL),
  fHnNprCRec(NULL) {
  AliLog::SetClassDebugLevel("AliEbyEPidRatioEffContExtra",10);
}

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


  for (Int_t ii = 0; ii < 2; ++ii) {
    for (Int_t kk = 0; kk < 4; ++kk)
      if (fLabelsRec[ii][kk]) delete[] fLabelsRec[ii][kk];
    if (fLabelsRec[ii]) delete[] fLabelsRec[ii];
  }
  if (fLabelsRec) delete[] fLabelsRec;

}


//________________________________________________________________________
void AliEbyEPidRatioEffContExtra::Process() {
  for(Int_t i = 0; i < 4; i++)  {
    FillMCLabels(i);
    FillMCEffHist(i);
  }
  return;
}      

//________________________________________________________________________
void AliEbyEPidRatioEffContExtra::Init() {
  fLabelsRec = new Int_t**[2];
  for (Int_t ii = 0 ; ii < 2; ++ii) {
    fLabelsRec[ii] = new Int_t*[4];
    for (Int_t kk = 0 ; kk < 4; ++kk)
      fLabelsRec[ii][kk] = NULL;
  }
  Printf(" >>>> AliEbyEPidRatioEffContExtra - inside");
}

//________________________________________________________________________
void AliEbyEPidRatioEffContExtra::CreateHistograms() {
   Int_t    binHnEffPIDMC[8] = {AliEbyEPidRatioHelper::fgkfHistNBinsCent, 
			       AliEbyEPidRatioHelper::fgkfHistNBinsSign, 
			       2, 2, 2,
			       AliEbyEPidRatioHelper::fgkfHistNBinsRap,     
			       AliEbyEPidRatioHelper::fgkfHistNBinsPhi, 
			       35};
  
  Double_t minHnEffPIDMC[8] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[0], 
			       AliEbyEPidRatioHelper::fgkfHistRangeSign[0], 
			       -0.5,     -0.5,     -0.5,
			       AliEbyEPidRatioHelper::fgkfHistRangeRap[0],  
			       AliEbyEPidRatioHelper::fgkfHistRangePhi[0], 
			       0.2};  
  
  Double_t maxHnEffPIDMC[8] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[1], 
			       AliEbyEPidRatioHelper::fgkfHistRangeSign[1], 
			       1.5,      1.5,      1.5,
			       AliEbyEPidRatioHelper::fgkfHistRangeRap[1],  
			       AliEbyEPidRatioHelper::fgkfHistRangePhi[1], 
			       3.0};  
  
  
  TString titilemc       = "cent:signMC:findable:recStatus:pidStatus:yMC:phiMC:ptMC";
 
  TString tiltlelaxmc[8]  = {"Centrality", 
			     "sign", 
			     "findable",
			     "recStatus",
			     "recPid",
			     "#it{y}_{MC}", 
			     "#varphi_{MC} (rad)", 
			     "#it{p}_{T,MC} (GeV/#it{c})"
  };

 
  //eff
  fHnNpiEMc  = new THnSparseF("hmNpiEffMc",titilemc.Data(),8,binHnEffPIDMC,minHnEffPIDMC, maxHnEffPIDMC);
  fHnNkaEMc  = new THnSparseF("hmNkaEffMc",titilemc.Data(),8,binHnEffPIDMC,minHnEffPIDMC, maxHnEffPIDMC);
  fHnNprEMc  = new THnSparseF("hmNprEffMc",titilemc.Data(),8,binHnEffPIDMC,minHnEffPIDMC, maxHnEffPIDMC);
  fHnNchEMc  = new THnSparseF("hmNchEffMc",titilemc.Data(),8,binHnEffPIDMC,minHnEffPIDMC, maxHnEffPIDMC);

  for (Int_t i = 0; i < 8; i++) {  
    fHnNchEMc->GetAxis(i)->SetTitle(tiltlelaxmc[i].Data());
    fHnNpiEMc->GetAxis(i)->SetTitle(tiltlelaxmc[i].Data());
    fHnNkaEMc->GetAxis(i)->SetTitle(tiltlelaxmc[i].Data());
    fHnNprEMc->GetAxis(i)->SetTitle(tiltlelaxmc[i].Data());
  }
  
  
  Int_t    binHnEffPID[5] = {AliEbyEPidRatioHelper::fgkfHistNBinsCent, 
			     AliEbyEPidRatioHelper::fgkfHistNBinsSign, 
			     AliEbyEPidRatioHelper::fgkfHistNBinsRap, 
			     AliEbyEPidRatioHelper::fgkfHistNBinsPhi,     
			     35};

  Double_t minHnEffPID[5] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[0], 
			     AliEbyEPidRatioHelper::fgkfHistRangeSign[0],			  
			     AliEbyEPidRatioHelper::fgkfHistRangeRap[0], 
			     AliEbyEPidRatioHelper::fgkfHistRangePhi[0],  
			     0.2};

  
  Double_t maxHnEffPID[5] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[1], 
			     AliEbyEPidRatioHelper::fgkfHistRangeSign[1], 
			     AliEbyEPidRatioHelper::fgkfHistRangeRap[1], 
			     AliEbyEPidRatioHelper::fgkfHistRangePhi[1],  
			     3.0};


  TString titilerec      = "cent:signRec:yRec:phiRec:ptRec";
  TString tiltlelaxrec[5]  = {"Centrality", 
			      "sign", 
			      "#it{y}_{Rec}", 
			      "#varphi_{Rec} (rad)", 
			      "#it{p}_{T,Rec} (GeV/#it{c})"};

 
  fHnNpiERec = new THnSparseF("hmNpiEffRec",titilerec.Data(),5,binHnEffPID, minHnEffPID, maxHnEffPID);
  fHnNkaERec = new THnSparseF("hmNkaEffRec",titilerec.Data(),5,binHnEffPID, minHnEffPID, maxHnEffPID);
  fHnNprERec = new THnSparseF("hmNprEffRec",titilerec.Data(),5,binHnEffPID, minHnEffPID, maxHnEffPID);
  fHnNchERec = new THnSparseF("hmNchEffRec",titilerec.Data(),5,binHnEffPID,minHnEffPID, maxHnEffPID);

  
  for (Int_t i = 0; i < 5; i++) { 
    fHnNchERec->GetAxis(i)->SetTitle(tiltlelaxrec[i].Data());
    fHnNpiERec->GetAxis(i)->SetTitle(tiltlelaxrec[i].Data());
    fHnNkaERec->GetAxis(i)->SetTitle(tiltlelaxrec[i].Data());
    fHnNprERec->GetAxis(i)->SetTitle(tiltlelaxrec[i].Data());
  }  

  //----- cont

  //----
  Int_t    binHnCont[6] = {AliEbyEPidRatioHelper::fgkfHistNBinsCent, 
			   AliEbyEPidRatioHelper::fgkfHistNBinsSign, 
			   8,                                        
			   AliEbyEPidRatioHelper::fgkfHistNBinsRap,  
			   AliEbyEPidRatioHelper::fgkfHistNBinsPhi, 
			   35};  
  
  Double_t minHnCont[6] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[0], 
			   AliEbyEPidRatioHelper::fgkfHistRangeSign[0], 
			   0.5,                                            
			   AliEbyEPidRatioHelper::fgkfHistRangeRap[0],  
			   AliEbyEPidRatioHelper::fgkfHistRangePhi[0], 
			   0.2};
  
  
  Double_t maxHnCont[6] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[1], 
			   AliEbyEPidRatioHelper::fgkfHistRangeSign[1],
			   8.5,                        
			   AliEbyEPidRatioHelper::fgkfHistRangeRap[1],  
			   AliEbyEPidRatioHelper::fgkfHistRangePhi[1], 
			   3.0};   


  TString titilecont     = "cent:signMC:contStatus:yMC:phiMC:ptMC";
  TString tiltlelaxcont[6] 
    = {"Centrality","sign","contPart","#it{y}_{MC}","#varphi_{MC} (rad)","#it{p}_{T,MC} (GeV/#it{c})"};
  
  fHnNpiCMc  = new THnSparseF("hmNpiContMc",titilecont.Data(),6,binHnCont,minHnCont, maxHnCont);
  fHnNkaCMc  = new THnSparseF("hmNkaContMc",titilecont.Data(),6,binHnCont,minHnCont, maxHnCont);
  fHnNprCMc  = new THnSparseF("hmNprContMc",titilecont.Data(),6,binHnCont,minHnCont, maxHnCont);
  fHnNchCMc  = new THnSparseF("hmNchContMc",titilecont.Data(),6,binHnCont,minHnCont, maxHnCont);

 for (Int_t i = 0; i < 6; i++) {  
    fHnNchCMc->GetAxis(i)->SetTitle(tiltlelaxcont[i].Data());
    fHnNpiCMc->GetAxis(i)->SetTitle(tiltlelaxcont[i].Data());
    fHnNkaCMc->GetAxis(i)->SetTitle(tiltlelaxcont[i].Data());
    fHnNprCMc->GetAxis(i)->SetTitle(tiltlelaxcont[i].Data());
  }

 
 fHnNpiCRec = new THnSparseF("hmNpiContRec",titilerec.Data(),5,binHnEffPID,minHnEffPID, maxHnEffPID);
 fHnNkaCRec = new THnSparseF("hmNkaContRec",titilerec.Data(),5,binHnEffPID,minHnEffPID, maxHnEffPID);
 fHnNprCRec = new THnSparseF("hmNprContRec",titilerec.Data(),5,binHnEffPID,minHnEffPID, maxHnEffPID);
 fHnNchCRec = new THnSparseF("hmNchContRec",titilerec.Data(),5,binHnEffPID,minHnEffPID, maxHnEffPID);

 for (Int_t i = 0; i < 5; i++) {  
   fHnNchCRec->GetAxis(i)->SetTitle(tiltlelaxrec[i].Data());
   fHnNpiCRec->GetAxis(i)->SetTitle(tiltlelaxrec[i].Data());
   fHnNkaCRec->GetAxis(i)->SetTitle(tiltlelaxrec[i].Data());
   fHnNprCRec->GetAxis(i)->SetTitle(tiltlelaxrec[i].Data());
 }  

 
  fHelper->BinLogAxis(fHnNchERec,4);
  fHelper->BinLogAxis(fHnNpiERec,4);
  fHelper->BinLogAxis(fHnNkaERec,4);
  fHelper->BinLogAxis(fHnNprERec,4);
  
  fHelper->BinLogAxis(fHnNchEMc,7);
  fHelper->BinLogAxis(fHnNpiEMc,7);
  fHelper->BinLogAxis(fHnNkaEMc,7);
  fHelper->BinLogAxis(fHnNprEMc,7);

  fHelper->BinLogAxis(fHnNchCMc,5);
  fHelper->BinLogAxis(fHnNpiCMc,5);
  fHelper->BinLogAxis(fHnNkaCMc,5);
  fHelper->BinLogAxis(fHnNprCMc,5);

  fHelper->BinLogAxis(fHnNchCRec,4);
  fHelper->BinLogAxis(fHnNpiCRec,4);
  fHelper->BinLogAxis(fHnNkaCRec,4);
  fHelper->BinLogAxis(fHnNprCRec,4);


  return;
}

//________________________________________________________________________
Int_t AliEbyEPidRatioEffContExtra::Setup() {
  for(Int_t i = 0; i < 4; i++) {
    fLabelsRec[0][i] = new Int_t[fNTracks];
    if(!fLabelsRec[0][i]) {
      AliError("Cannot create fLabelsRec[0]");
      return -1;
    }
    
    fLabelsRec[1][i] = new Int_t[fNTracks];
    if(!fLabelsRec[1][i]) {
      AliError("Cannot create fLabelsRec[1] for PID");
      return -1;
    }
    
    for(Int_t ii = 0; ii < fNTracks; ++ii) {
      fLabelsRec[0][i][ii] = 0;
      fLabelsRec[1][i][ii] = 0;
    }
  }
 
  return 0;
}

//________________________________________________________________________
void AliEbyEPidRatioEffContExtra::Reset() {
  // -- Reset eventwise

  for(Int_t i = 0; i < 4; i++) {
    for (Int_t ii = 0; ii < 2 ; ++ii) {
      if (fLabelsRec[ii][i])
	delete[] fLabelsRec[ii][i];
      fLabelsRec[ii][i] = NULL;
    }
  }
}

//________________________________________________________________________
void AliEbyEPidRatioEffContExtra::FillMCLabels(Int_t ipid) {
  // 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

  // -- Get ranges for AOD particles
  Float_t etaRange[2];
  fESDTrackCuts->GetEtaRange(etaRange[0],etaRange[1]);

  Float_t ptRange[2];
  fESDTrackCuts->GetPtRange(ptRange[0],ptRange[1]);

  // -- Track Loop
  for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
    AliVTrack *track = (fESD) ? static_cast<AliVTrack*>(fESD->GetTrack(idxTrack)) : static_cast<AliVTrack*>(fAOD->GetTrack(idxTrack)); 

    // -- Check if track is accepted for basic parameters
    if (!fHelper->IsTrackAcceptedBasicCharged(track))
      continue;
    
    // -- Check if accepted - ESD
    if (fESD && !fESDTrackCuts->AcceptTrack(dynamic_cast<AliESDtrack*>(track)))
      continue;
    
    // -- Check if accepted - AOD
    if (fAOD){
      AliAODTrack * trackAOD = dynamic_cast<AliAODTrack*>(track);
      
      if (!trackAOD) {
	AliError("Pointer to dynamic_cast<AliAODTrack*>(track) = ZERO");
	continue;
      }
      if (!trackAOD->TestFilterBit(fAODtrackCutBit))
	continue;

      // -- Check if in pT and eta range (is done in ESDTrackCuts for ESDs)
      if(!(track->Pt() > ptRange[0] && track->Pt() <= ptRange[1] && TMath::Abs(track->Eta()) <= etaRange[1]))
	continue;
    }
   
    // -- Check if accepted in rapidity window -- for identified particles
    //    for charged particles -- eta check is done in the trackcuts
    Double_t yP;
    if (fHelper->GetUsePID(ipid) && !fHelper->IsTrackAcceptedRapidity(track, yP, ipid))
      continue;

 
    if (!fHelper->IsTrackAcceptedDCA(track))
      continue;

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

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

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

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

  return;
}

//________________________________________________________________________
void AliEbyEPidRatioEffContExtra::CheckContTrack(AliVTrack *track, Int_t ipid) {
  // Check if particle is contamination or correctly identified for ESDs and AODs
  // Check for missidentified primaries and secondaries
  // Fill contamination THnSparse

  Int_t label     = TMath::Abs(track->GetLabel()); 
  Float_t signRec = track->Charge();

  
  AliVParticle* particle = (fESD) ? fMCEvent->GetTrack(label) : static_cast<AliVParticle*>(fArrayMC->At(label));
  if (!particle)
    return;

  Bool_t isPhysicalPrimary = (fESD) ? fStack->IsPhysicalPrimary(label): (static_cast<AliAODMCParticle*>(particle))->IsPhysicalPrimary();

  // -- Check if correctly identified 
  //    > return if correctly identified -> all ok, no action neededin this method
  //    > if PID required check -> for the correct (signed pdgcode) particle
  //    > no PID just check for primary 
  if (fHelper->GetUsePID(ipid)) {
    if (particle->PdgCode() == (signRec*fHelper->GetPdg(ipid)))
      if (isPhysicalPrimary)
	return;
  }
  else {
    if (isPhysicalPrimary)
      return;
  }

  // -- Check if secondaries from material or weak decay
  Bool_t isSecondaryFromWeakDecay = (fESD) ? fStack->IsSecondaryFromWeakDecay(label) : (static_cast<AliAODMCParticle*>(particle))->IsSecondaryFromWeakDecay();
  Bool_t isSecondaryFromMaterial  = (fESD) ? fStack->IsSecondaryFromMaterial(label)  : (static_cast<AliAODMCParticle*>(particle))->IsSecondaryFromMaterial();

  // -- Get PDG Charge of contaminating particle
  Float_t signMC = 0.;
  if      (particle->Charge() == 0.) signMC =  0.;
  else if (particle->Charge() <  0.) signMC = -1.;	
  else if (particle->Charge() >  0.) signMC =  1.;	

  // -- Get contaminating particle
  Double_t contPart = 0;
  if        (isSecondaryFromWeakDecay)                contPart = 7; // probeParticle from WeakDecay
  else if   (isSecondaryFromMaterial)                 contPart = 8; // probeParticle from Material
  else {
    if      (TMath::Abs(particle->PdgCode()) ==  211) contPart = 1; // pion
    else if (TMath::Abs(particle->PdgCode()) ==  321) contPart = 2; // kaon
    else if (TMath::Abs(particle->PdgCode()) == 2212) contPart = 3; // proton
    else if (TMath::Abs(particle->PdgCode()) ==   11) contPart = 4; // electron
    else if (TMath::Abs(particle->PdgCode()) ==   13) contPart = 5; // muon
    else                                              contPart = 6; // other
  }
  
  // -- Get Reconstructed y
  //    yRec = y for identified particles | yRec = eta for charged particles
  Double_t yRec  = 0.;
  fHelper->IsTrackAcceptedRapidity(track, yRec, ipid); 

 
  Double_t yetapid = (ipid == 0 ) ? particle->Eta() : particle->Y();
  Double_t yeta    = (ipid == 0 ) ? track->Eta() : yRec;

  Double_t hnContMc[6]  = {fCentralityBin,signMC,contPart,yetapid,particle->Phi(),particle->Pt()};
  Double_t hnContRec[5] = {fCentralityBin,signRec, yeta,track->Phi(),track->Pt()};
  
  if (ipid == 0) { fHnNchCRec->Fill(hnContRec); fHnNchCMc->Fill(hnContMc); }
  else if (ipid == 1) { fHnNpiCRec->Fill(hnContRec); fHnNpiCMc->Fill(hnContMc); }
  else if (ipid == 2) { fHnNkaCRec->Fill(hnContRec); fHnNkaCMc->Fill(hnContMc); }
  else if (ipid == 3) { fHnNprCRec->Fill(hnContRec); fHnNprCMc->Fill(hnContMc); }
  
  
}

//________________________________________________________________________
void AliEbyEPidRatioEffContExtra::FillMCEffHist(Int_t ipid) {
  // Fill efficiency THnSparse for ESDs

  Float_t etaRange[2];
  fESDTrackCuts->GetEtaRange(etaRange[0],etaRange[1]);

  Int_t nPart  = (fESD) ? fStack->GetNprimary() : fArrayMC->GetEntriesFast();

  for (Int_t idxMC = 0; idxMC < nPart; ++idxMC) {
    AliVParticle* particle = (fESD) ? fMCEvent->GetTrack(idxMC) : static_cast<AliVParticle*>(fArrayMC->At(idxMC));

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

    // -- Check if accepted in rapidity window -- for identified particles
    Double_t yMC;
    if (fHelper->GetUsePID(ipid) && !fHelper->IsParticleAcceptedRapidity(particle, yMC, ipid))
      continue;

    // -- Check if accepted in eta window -- for charged particles
    if (!fHelper->GetUsePID(ipid) && TMath::Abs(particle->Eta()) > etaRange[1])
      continue;

    // -- Check if probeParticle / anti-probeParticle 
    //    > skip check if PID is not required
    if (fHelper->GetUsePID(ipid) && TMath::Abs(particle->PdgCode()) != fHelper->GetPdg(ipid))
      continue;
    
    // -- Get sign of particle
    Float_t signMC    = (particle->PdgCode() < 0) ? -1. : 1.;

    // -- Get if particle is findable --- not availible for AODs yet
    Float_t findable  = (fESD) ? Float_t(fHelper->IsParticleFindable(idxMC)) : 1.;

    // -- 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.;
    Double_t yRec   = 0.;
    Float_t signRec = 0.;

    // -- Loop over all labels
    for (Int_t idxRec=0; idxRec < fNTracks; ++idxRec) {
      if (idxMC == fLabelsRec[0][ipid][idxRec]) {
	recStatus = 1.;
	
	if (idxMC == fLabelsRec[1][ipid][idxRec]) recPid = 1.;
	
        AliVTrack *track = NULL;
        if(fESD)
          track = fESD->GetTrack(idxRec);
        else if(fAOD)
          track = fAOD->GetTrack(idxRec);
	
        if (track) {
	  etaRec  = track->Eta();
          phiRec  = track->Phi();         
          ptRec   = track->Pt();
	  signRec = track->Charge();
          fHelper->IsTrackAcceptedRapidity(track, yRec, ipid); 
	  Double_t yeta    = (ipid == 0 ) ? etaRec : yRec;
	  Double_t hneffRec[5] = {fCentralityBin,signRec, yeta,phiRec,ptRec};
	  if (ipid == 0) { fHnNchERec->Fill(hneffRec); }
	  else if (ipid == 1) { fHnNpiERec->Fill(hneffRec); }
	  else if (ipid == 2) { fHnNkaERec->Fill(hneffRec); }
	  else if (ipid == 3) { fHnNprERec->Fill(hneffRec); }

        }     
        break;
      }
    } // for (Int_t idxRec=0; idxRec < fNTracks; ++idxRec) {  
    /*
    Double_t deltaPhi = particle->Phi()-phiRec;
    if (TMath::Abs(deltaPhi) > TMath::TwoPi()) {
      if (deltaPhi < 0)
	deltaPhi += TMath::TwoPi();
      else
    	deltaPhi -= TMath::TwoPi();
    }
    */
  
    Double_t yetapid = (ipid == 0 ) ? particle->Eta() : particle->Y();
   

    // Printf("%2d  %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f", ipid, yetapid, yeta, particle->Eta(), particle->Y(),  etaRec, yRec);    

    Double_t hneffMc[8]  = {fCentralityBin,signMC,findable, recStatus, recPid,yetapid,particle->Phi(),particle->Pt()};
    

    if (ipid == 0) { fHnNchEMc->Fill(hneffMc); }
    else if (ipid == 1) { fHnNpiEMc->Fill(hneffMc); }
    else if (ipid == 2) { fHnNkaEMc->Fill(hneffMc); }
    else if (ipid == 3) { fHnNprEMc->Fill(hneffMc); }
   

  } // for (Int_t idxMC = 0; idxMC < nPart; ++idxMC) {
  
  return;
}
Commit	Line	Data
5f071dcf	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: ALICE Offline. *
	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	//=========================================================================//
	18	// AliEbyE Analysis for Particle Ratio Fluctuation //
	19	// Deepika Rathee \| Satyajit Jena //
	20	// drathee@cern.ch \| sjena@cern.ch //
	21	// Date: Wed Jul 9 18:38:30 CEST 2014 //
	22	// New approch to find particle ratio to reduce memory //
	23	// (Test Only) //
df4f740f	24	// Copied from NetParticle Classes
1063a5fa	25	// Origin: Jochen Thaeder <jochen@thaeder.de>
1063a5fa	26	// Michael Weber <m.weber@cern.ch>
5f071dcf	27	//=========================================================================//
	28
	29	#include "TMath.h"
	30	#include "TAxis.h"
	31
	32	#include "AliESDEvent.h"
	33	#include "AliStack.h"
	34	#include "AliMCEvent.h"
	35	#include "AliESDtrackCuts.h"
	36	#include "AliAODEvent.h"
	37	#include "AliAODMCParticle.h"
	38
	39	#include "AliEbyEPidRatioEffContExtra.h"
	40
	41	using namespace std;
	42
	43	ClassImp(AliEbyEPidRatioEffContExtra)
	44
	45	//________________________________________________________________________
	46	AliEbyEPidRatioEffContExtra::AliEbyEPidRatioEffContExtra() :
	47	AliEbyEPidRatioBase("EffCont", "EffCont"),
	48	fLabelsRec(NULL),
	49
	50	fHnNchEMc(NULL),
	51	fHnNchERec(NULL),
	52	fHnNpiEMc(NULL),
	53	fHnNpiERec(NULL),
	54	fHnNkaEMc(NULL),
	55	fHnNkaERec(NULL),
	56	fHnNprEMc(NULL),
	57	fHnNprERec(NULL),
	58
	59	fHnNchCMc(NULL),
	60	fHnNchCRec(NULL),
	61	fHnNpiCMc(NULL),
	62	fHnNpiCRec(NULL),
	63	fHnNkaCMc(NULL),
	64	fHnNkaCRec(NULL),
	65	fHnNprCMc(NULL),
	66	fHnNprCRec(NULL) {
	67	AliLog::SetClassDebugLevel("AliEbyEPidRatioEffContExtra",10);
	68	}
	69
	70	//________________________________________________________________________
	71	AliEbyEPidRatioEffContExtra::~AliEbyEPidRatioEffContExtra() {
	72	// Destructor
	73
	74
	75	for (Int_t ii = 0; ii < 2; ++ii) {
	76	for (Int_t kk = 0; kk < 4; ++kk)
	77	if (fLabelsRec[ii][kk]) delete[] fLabelsRec[ii][kk];
	78	if (fLabelsRec[ii]) delete[] fLabelsRec[ii];
	79	}
	80	if (fLabelsRec) delete[] fLabelsRec;
	81
5f071dcf	82	}
	83
	84
	85	//________________________________________________________________________
	86	void AliEbyEPidRatioEffContExtra::Process() {
5f071dcf	87	for(Int_t i = 0; i < 4; i++) {
5f071dcf	88	FillMCLabels(i);
5f071dcf	89	FillMCEffHist(i);
	90	}
	91	return;
	92	}
	93
	94	//________________________________________________________________________
	95	void AliEbyEPidRatioEffContExtra::Init() {
5f071dcf	96	fLabelsRec = new Int_t**[2];
	97	for (Int_t ii = 0 ; ii < 2; ++ii) {
	98	fLabelsRec[ii] = new Int_t*[4];
	99	for (Int_t kk = 0 ; kk < 4; ++kk)
	100	fLabelsRec[ii][kk] = NULL;
	101	}
5f071dcf	102	Printf(" >>>> AliEbyEPidRatioEffContExtra - inside");
	103	}
	104
	105	//________________________________________________________________________
	106	void AliEbyEPidRatioEffContExtra::CreateHistograms() {
d61393df	107	Int_t binHnEffPIDMC[8] = {AliEbyEPidRatioHelper::fgkfHistNBinsCent,
5f071dcf	108	AliEbyEPidRatioHelper::fgkfHistNBinsSign,
	109	2, 2, 2,
	110	AliEbyEPidRatioHelper::fgkfHistNBinsRap,
	111	AliEbyEPidRatioHelper::fgkfHistNBinsPhi,
eb231666	112	35};
5f071dcf	113
	114	Double_t minHnEffPIDMC[8] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[0],
	115	AliEbyEPidRatioHelper::fgkfHistRangeSign[0],
	116	-0.5, -0.5, -0.5,
	117	AliEbyEPidRatioHelper::fgkfHistRangeRap[0],
	118	AliEbyEPidRatioHelper::fgkfHistRangePhi[0],
eb231666	119	0.2};
5f071dcf	120
	121	Double_t maxHnEffPIDMC[8] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[1],
	122	AliEbyEPidRatioHelper::fgkfHistRangeSign[1],
	123	1.5, 1.5, 1.5,
	124	AliEbyEPidRatioHelper::fgkfHistRangeRap[1],
	125	AliEbyEPidRatioHelper::fgkfHistRangePhi[1],
eb231666	126	3.0};
5f071dcf	127
	128
	129
	130	TString titilemc = "cent:signMC:findable:recStatus:pidStatus:yMC:phiMC:ptMC";
	131
	132	TString tiltlelaxmc[8] = {"Centrality",
	133	"sign",
	134	"findable",
	135	"recStatus",
	136	"recPid",
	137	"#it{y}_{MC}",
	138	"#varphi_{MC} (rad)",
	139	"#it{p}_{T,MC} (GeV/#it{c})"
	140	};
	141
	142
	143	//eff
	144	fHnNpiEMc = new THnSparseF("hmNpiEffMc",titilemc.Data(),8,binHnEffPIDMC,minHnEffPIDMC, maxHnEffPIDMC);
	145	fHnNkaEMc = new THnSparseF("hmNkaEffMc",titilemc.Data(),8,binHnEffPIDMC,minHnEffPIDMC, maxHnEffPIDMC);
	146	fHnNprEMc = new THnSparseF("hmNprEffMc",titilemc.Data(),8,binHnEffPIDMC,minHnEffPIDMC, maxHnEffPIDMC);
	147	fHnNchEMc = new THnSparseF("hmNchEffMc",titilemc.Data(),8,binHnEffPIDMC,minHnEffPIDMC, maxHnEffPIDMC);
	148
	149	for (Int_t i = 0; i < 8; i++) {
	150	fHnNchEMc->GetAxis(i)->SetTitle(tiltlelaxmc[i].Data());
	151	fHnNpiEMc->GetAxis(i)->SetTitle(tiltlelaxmc[i].Data());
	152	fHnNkaEMc->GetAxis(i)->SetTitle(tiltlelaxmc[i].Data());
	153	fHnNprEMc->GetAxis(i)->SetTitle(tiltlelaxmc[i].Data());
	154	}
	155
	156
	157
	158
	159	Int_t binHnEffPID[5] = {AliEbyEPidRatioHelper::fgkfHistNBinsCent,
	160	AliEbyEPidRatioHelper::fgkfHistNBinsSign,
	161	AliEbyEPidRatioHelper::fgkfHistNBinsRap,
	162	AliEbyEPidRatioHelper::fgkfHistNBinsPhi,
ebe7e6bf	163	35};
5f071dcf	164
	165	Double_t minHnEffPID[5] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[0],
	166	AliEbyEPidRatioHelper::fgkfHistRangeSign[0],
	167	AliEbyEPidRatioHelper::fgkfHistRangeRap[0],
	168	AliEbyEPidRatioHelper::fgkfHistRangePhi[0],
ebe7e6bf	169	0.2};
5f071dcf	170
	171
	172	Double_t maxHnEffPID[5] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[1],
ebe7e6bf	173	AliEbyEPidRatioHelper::fgkfHistRangeSign[1],
	174	AliEbyEPidRatioHelper::fgkfHistRangeRap[1],
	175	AliEbyEPidRatioHelper::fgkfHistRangePhi[1],
	176	3.0};
5f071dcf	177
	178
	179	TString titilerec = "cent:signRec:yRec:phiRec:ptRec";
	180	TString tiltlelaxrec[5] = {"Centrality",
	181	"sign",
	182	"#it{y}_{Rec}",
	183	"#varphi_{Rec} (rad)",
	184	"#it{p}_{T,Rec} (GeV/#it{c})"};
	185
	186
	187	fHnNpiERec = new THnSparseF("hmNpiEffRec",titilerec.Data(),5,binHnEffPID, minHnEffPID, maxHnEffPID);
	188	fHnNkaERec = new THnSparseF("hmNkaEffRec",titilerec.Data(),5,binHnEffPID, minHnEffPID, maxHnEffPID);
	189	fHnNprERec = new THnSparseF("hmNprEffRec",titilerec.Data(),5,binHnEffPID, minHnEffPID, maxHnEffPID);
	190	fHnNchERec = new THnSparseF("hmNchEffRec",titilerec.Data(),5,binHnEffPID,minHnEffPID, maxHnEffPID);
	191
	192
	193
	194	for (Int_t i = 0; i < 5; i++) {
	195	fHnNchERec->GetAxis(i)->SetTitle(tiltlelaxrec[i].Data());
	196	fHnNpiERec->GetAxis(i)->SetTitle(tiltlelaxrec[i].Data());
	197	fHnNkaERec->GetAxis(i)->SetTitle(tiltlelaxrec[i].Data());
	198	fHnNprERec->GetAxis(i)->SetTitle(tiltlelaxrec[i].Data());
	199	}
	200
	201	//----- cont
	202
	203	//----
	204	Int_t binHnCont[6] = {AliEbyEPidRatioHelper::fgkfHistNBinsCent,
	205	AliEbyEPidRatioHelper::fgkfHistNBinsSign,
	206	8,
	207	AliEbyEPidRatioHelper::fgkfHistNBinsRap,
	208	AliEbyEPidRatioHelper::fgkfHistNBinsPhi,
eb231666	209	35};
5f071dcf	210
	211	Double_t minHnCont[6] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[0],
	212	AliEbyEPidRatioHelper::fgkfHistRangeSign[0],
	213	0.5,
	214	AliEbyEPidRatioHelper::fgkfHistRangeRap[0],
	215	AliEbyEPidRatioHelper::fgkfHistRangePhi[0],
eb231666	216	0.2};
5f071dcf	217
	218
	219
	220	Double_t maxHnCont[6] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[1],
	221	AliEbyEPidRatioHelper::fgkfHistRangeSign[1],
	222	8.5,
	223	AliEbyEPidRatioHelper::fgkfHistRangeRap[1],
	224	AliEbyEPidRatioHelper::fgkfHistRangePhi[1],
eb231666	225	3.0};
5f071dcf	226
	227
	228
	229
	230	TString titilecont = "cent:signMC:contStatus:yMC:phiMC:ptMC";
	231	TString tiltlelaxcont[6]
	232	= {"Centrality","sign","contPart","#it{y}_{MC}","#varphi_{MC} (rad)","#it{p}_{T,MC} (GeV/#it{c})"};
	233
61627432	234	fHnNpiCMc = new THnSparseF("hmNpiContMc",titilecont.Data(),6,binHnCont,minHnCont, maxHnCont);
	235	fHnNkaCMc = new THnSparseF("hmNkaContMc",titilecont.Data(),6,binHnCont,minHnCont, maxHnCont);
	236	fHnNprCMc = new THnSparseF("hmNprContMc",titilecont.Data(),6,binHnCont,minHnCont, maxHnCont);
	237	fHnNchCMc = new THnSparseF("hmNchContMc",titilecont.Data(),6,binHnCont,minHnCont, maxHnCont);
5f071dcf	238
	239	for (Int_t i = 0; i < 6; i++) {
	240	fHnNchCMc->GetAxis(i)->SetTitle(tiltlelaxcont[i].Data());
	241	fHnNpiCMc->GetAxis(i)->SetTitle(tiltlelaxcont[i].Data());
	242	fHnNkaCMc->GetAxis(i)->SetTitle(tiltlelaxcont[i].Data());
	243	fHnNprCMc->GetAxis(i)->SetTitle(tiltlelaxcont[i].Data());
	244	}
	245
	246
	247
	248	fHnNpiCRec = new THnSparseF("hmNpiContRec",titilerec.Data(),5,binHnEffPID,minHnEffPID, maxHnEffPID);
	249	fHnNkaCRec = new THnSparseF("hmNkaContRec",titilerec.Data(),5,binHnEffPID,minHnEffPID, maxHnEffPID);
	250	fHnNprCRec = new THnSparseF("hmNprContRec",titilerec.Data(),5,binHnEffPID,minHnEffPID, maxHnEffPID);
	251	fHnNchCRec = new THnSparseF("hmNchContRec",titilerec.Data(),5,binHnEffPID,minHnEffPID, maxHnEffPID);
	252
	253	for (Int_t i = 0; i < 5; i++) {
	254	fHnNchCRec->GetAxis(i)->SetTitle(tiltlelaxrec[i].Data());
	255	fHnNpiCRec->GetAxis(i)->SetTitle(tiltlelaxrec[i].Data());
	256	fHnNkaCRec->GetAxis(i)->SetTitle(tiltlelaxrec[i].Data());
	257	fHnNprCRec->GetAxis(i)->SetTitle(tiltlelaxrec[i].Data());
	258	}
	259
	260
	261	fHelper->BinLogAxis(fHnNchERec,4);
	262	fHelper->BinLogAxis(fHnNpiERec,4);
	263	fHelper->BinLogAxis(fHnNkaERec,4);
	264	fHelper->BinLogAxis(fHnNprERec,4);
	265
	266	fHelper->BinLogAxis(fHnNchEMc,7);
	267	fHelper->BinLogAxis(fHnNpiEMc,7);
	268	fHelper->BinLogAxis(fHnNkaEMc,7);
	269	fHelper->BinLogAxis(fHnNprEMc,7);
	270
	271	fHelper->BinLogAxis(fHnNchCMc,5);
	272	fHelper->BinLogAxis(fHnNpiCMc,5);
	273	fHelper->BinLogAxis(fHnNkaCMc,5);
	274	fHelper->BinLogAxis(fHnNprCMc,5);
	275
	276	fHelper->BinLogAxis(fHnNchCRec,4);
	277	fHelper->BinLogAxis(fHnNpiCRec,4);
	278	fHelper->BinLogAxis(fHnNkaCRec,4);
	279	fHelper->BinLogAxis(fHnNprCRec,4);
	280
	281
	282	return;
	283	}
	284
	285	//________________________________________________________________________
	286	Int_t AliEbyEPidRatioEffContExtra::Setup() {
5f071dcf	287	for(Int_t i = 0; i < 4; i++) {
	288	fLabelsRec[0][i] = new Int_t[fNTracks];
	289	if(!fLabelsRec[0][i]) {
	290	AliError("Cannot create fLabelsRec[0]");
	291	return -1;
	292	}
	293
	294	fLabelsRec[1][i] = new Int_t[fNTracks];
	295	if(!fLabelsRec[1][i]) {
	296	AliError("Cannot create fLabelsRec[1] for PID");
	297	return -1;
	298	}
	299
	300	for(Int_t ii = 0; ii < fNTracks; ++ii) {
	301	fLabelsRec[0][i][ii] = 0;
	302	fLabelsRec[1][i][ii] = 0;
	303	}
	304	}
d61393df	305
5f071dcf	306	return 0;
	307	}
	308
	309	//________________________________________________________________________
	310	void AliEbyEPidRatioEffContExtra::Reset() {
	311	// -- Reset eventwise
	312
	313	for(Int_t i = 0; i < 4; i++) {
	314	for (Int_t ii = 0; ii < 2 ; ++ii) {
	315	if (fLabelsRec[ii][i])
	316	delete[] fLabelsRec[ii][i];
	317	fLabelsRec[ii][i] = NULL;
	318	}
	319	}
5f071dcf	320	}
5f071dcf	321
5f071dcf	322	//________________________________________________________________________
	323	void AliEbyEPidRatioEffContExtra::FillMCLabels(Int_t ipid) {
	324	// Fill MC labels
	325	// Loop over ESD tracks and fill arrays with MC lables
	326	// fLabelsRec[0] : all Tracks
	327	// fLabelsRec[1] : all Tracks accepted by PID of TPC
	328	// Check every accepted track if correctly identified
	329	// otherwise check for contamination
	330
	331	// -- Get ranges for AOD particles
	332	Float_t etaRange[2];
	333	fESDTrackCuts->GetEtaRange(etaRange[0],etaRange[1]);
	334
	335	Float_t ptRange[2];
	336	fESDTrackCuts->GetPtRange(ptRange[0],ptRange[1]);
	337
	338	// -- Track Loop
	339	for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
	340	AliVTrack track = (fESD) ? static_cast<AliVTrack>(fESD->GetTrack(idxTrack)) : static_cast<AliVTrack*>(fAOD->GetTrack(idxTrack));
	341
	342	// -- Check if track is accepted for basic parameters
	343	if (!fHelper->IsTrackAcceptedBasicCharged(track))
	344	continue;
	345
	346	// -- Check if accepted - ESD
	347	if (fESD && !fESDTrackCuts->AcceptTrack(dynamic_cast<AliESDtrack*>(track)))
	348	continue;
	349
	350	// -- Check if accepted - AOD
	351	if (fAOD){
	352	AliAODTrack * trackAOD = dynamic_cast<AliAODTrack*>(track);
	353
	354	if (!trackAOD) {
	355	AliError("Pointer to dynamic_cast<AliAODTrack*>(track) = ZERO");
	356	continue;
	357	}
	358	if (!trackAOD->TestFilterBit(fAODtrackCutBit))
	359	continue;
	360
	361	// -- Check if in pT and eta range (is done in ESDTrackCuts for ESDs)
	362	if(!(track->Pt() > ptRange[0] && track->Pt() <= ptRange[1] && TMath::Abs(track->Eta()) <= etaRange[1]))
	363	continue;
	364	}
	365
	366	// -- Check if accepted in rapidity window -- for identified particles
	367	// for charged particles -- eta check is done in the trackcuts
	368	Double_t yP;
	369	if (fHelper->GetUsePID(ipid) && !fHelper->IsTrackAcceptedRapidity(track, yP, ipid))
	370	continue;
	371
	372
5f071dcf	373	if (!fHelper->IsTrackAcceptedDCA(track))
	374	continue;
	375
	376	Int_t label = TMath::Abs(track->GetLabel());
	377
	378	// continue;
	379	// -- Fill Label of all reconstructed
	380	fLabelsRec[0][ipid][idxTrack] = label;
	381
	382	// -- Check if accepted by PID from TPC or TPC+TOF
	383	Double_t pid[3];
	384	if (!fHelper->IsTrackAcceptedPID(track, pid, fHelper->GetParticleSpecies(ipid))) // check it
	385	continue;
	386
	387
	388
	389	// -- Fill Label of all reconstructed && recPid_TPC+TOF
	390	fLabelsRec[1][ipid][idxTrack] = label;
	391
	392	// -- Check for contamination and fill contamination THnSparse
	393	CheckContTrack(track, ipid);
	394
	395	} // for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
	396
	397	return;
	398	}
	399
	400	//________________________________________________________________________
	401	void AliEbyEPidRatioEffContExtra::CheckContTrack(AliVTrack *track, Int_t ipid) {
	402	// Check if particle is contamination or correctly identified for ESDs and AODs
	403	// Check for missidentified primaries and secondaries
	404	// Fill contamination THnSparse
	405
	406	Int_t label = TMath::Abs(track->GetLabel());
	407	Float_t signRec = track->Charge();
	408
	409
	410	AliVParticle* particle = (fESD) ? fMCEvent->GetTrack(label) : static_cast<AliVParticle*>(fArrayMC->At(label));
	411	if (!particle)
	412	return;
	413
	414	Bool_t isPhysicalPrimary = (fESD) ? fStack->IsPhysicalPrimary(label): (static_cast<AliAODMCParticle*>(particle))->IsPhysicalPrimary();
1063a5fa	415
5f071dcf	416	// -- Check if correctly identified
	417	// > return if correctly identified -> all ok, no action neededin this method
	418	// > if PID required check -> for the correct (signed pdgcode) particle
	419	// > no PID just check for primary
	420	if (fHelper->GetUsePID(ipid)) {
	421	if (particle->PdgCode() == (signRec*fHelper->GetPdg(ipid)))
	422	if (isPhysicalPrimary)
	423	return;
	424	}
	425	else {
	426	if (isPhysicalPrimary)
	427	return;
	428	}
	429
	430	// -- Check if secondaries from material or weak decay
	431	Bool_t isSecondaryFromWeakDecay = (fESD) ? fStack->IsSecondaryFromWeakDecay(label) : (static_cast<AliAODMCParticle*>(particle))->IsSecondaryFromWeakDecay();
	432	Bool_t isSecondaryFromMaterial = (fESD) ? fStack->IsSecondaryFromMaterial(label) : (static_cast<AliAODMCParticle*>(particle))->IsSecondaryFromMaterial();
	433
	434	// -- Get PDG Charge of contaminating particle
	435	Float_t signMC = 0.;
	436	if (particle->Charge() == 0.) signMC = 0.;
	437	else if (particle->Charge() < 0.) signMC = -1.;
	438	else if (particle->Charge() > 0.) signMC = 1.;
	439
	440	// -- Get contaminating particle
	441	Double_t contPart = 0;
	442	if (isSecondaryFromWeakDecay) contPart = 7; // probeParticle from WeakDecay
	443	else if (isSecondaryFromMaterial) contPart = 8; // probeParticle from Material
	444	else {
	445	if (TMath::Abs(particle->PdgCode()) == 211) contPart = 1; // pion
	446	else if (TMath::Abs(particle->PdgCode()) == 321) contPart = 2; // kaon
	447	else if (TMath::Abs(particle->PdgCode()) == 2212) contPart = 3; // proton
	448	else if (TMath::Abs(particle->PdgCode()) == 11) contPart = 4; // electron
	449	else if (TMath::Abs(particle->PdgCode()) == 13) contPart = 5; // muon
	450	else contPart = 6; // other
	451	}
	452
	453	// -- Get Reconstructed y
	454	// yRec = y for identified particles \| yRec = eta for charged particles
	455	Double_t yRec = 0.;
	456	fHelper->IsTrackAcceptedRapidity(track, yRec, ipid);
	457
5f071dcf	458
	459	Double_t yetapid = (ipid == 0 ) ? particle->Eta() : particle->Y();
	460	Double_t yeta = (ipid == 0 ) ? track->Eta() : yRec;
	461
	462	Double_t hnContMc[6] = {fCentralityBin,signMC,contPart,yetapid,particle->Phi(),particle->Pt()};
	463	Double_t hnContRec[5] = {fCentralityBin,signRec, yeta,track->Phi(),track->Pt()};
	464
	465	if (ipid == 0) { fHnNchCRec->Fill(hnContRec); fHnNchCMc->Fill(hnContMc); }
1063a5fa	466	else if (ipid == 1) { fHnNpiCRec->Fill(hnContRec); fHnNpiCMc->Fill(hnContMc); }
	467	else if (ipid == 2) { fHnNkaCRec->Fill(hnContRec); fHnNkaCMc->Fill(hnContMc); }
	468	else if (ipid == 3) { fHnNprCRec->Fill(hnContRec); fHnNprCMc->Fill(hnContMc); }
5f071dcf	469
	470
	471	}
	472
	473	//________________________________________________________________________
	474	void AliEbyEPidRatioEffContExtra::FillMCEffHist(Int_t ipid) {
	475	// Fill efficiency THnSparse for ESDs
	476
	477	Float_t etaRange[2];
	478	fESDTrackCuts->GetEtaRange(etaRange[0],etaRange[1]);
	479
	480	Int_t nPart = (fESD) ? fStack->GetNprimary() : fArrayMC->GetEntriesFast();
	481
	482	for (Int_t idxMC = 0; idxMC < nPart; ++idxMC) {
	483	AliVParticle* particle = (fESD) ? fMCEvent->GetTrack(idxMC) : static_cast<AliVParticle*>(fArrayMC->At(idxMC));
	484
	485	// -- Check basic MC properties -> charged physical primary
	486	if (!fHelper->IsParticleAcceptedBasicCharged(particle, idxMC))
	487	continue;
	488
	489	// -- Check if accepted in rapidity window -- for identified particles
	490	Double_t yMC;
	491	if (fHelper->GetUsePID(ipid) && !fHelper->IsParticleAcceptedRapidity(particle, yMC, ipid))
	492	continue;
	493
	494	// -- Check if accepted in eta window -- for charged particles
	495	if (!fHelper->GetUsePID(ipid) && TMath::Abs(particle->Eta()) > etaRange[1])
	496	continue;
	497
	498	// -- Check if probeParticle / anti-probeParticle
	499	// > skip check if PID is not required
	500	if (fHelper->GetUsePID(ipid) && TMath::Abs(particle->PdgCode()) != fHelper->GetPdg(ipid))
	501	continue;
	502
	503	// -- Get sign of particle
	504	Float_t signMC = (particle->PdgCode() < 0) ? -1. : 1.;
	505
	506	// -- Get if particle is findable --- not availible for AODs yet
	507	Float_t findable = (fESD) ? Float_t(fHelper->IsParticleFindable(idxMC)) : 1.;
	508
	509	// -- Get recStatus and pidStatus
	510	Float_t recStatus = 0.;
	511	Float_t recPid = 0.;
	512
	513	// -- Get Reconstructed values
	514	Float_t etaRec = 0.;
	515	Float_t phiRec = 0.;
	516	Float_t ptRec = 0.;
	517	Double_t yRec = 0.;
	518	Float_t signRec = 0.;
	519
	520	// -- Loop over all labels
	521	for (Int_t idxRec=0; idxRec < fNTracks; ++idxRec) {
d61393df	522	if (idxMC == fLabelsRec[0][ipid][idxRec]) {
5f071dcf	523	recStatus = 1.;
5f071dcf	524
1063a5fa	525	if (idxMC == fLabelsRec[1][ipid][idxRec]) recPid = 1.;
5f071dcf	526
	527	AliVTrack *track = NULL;
	528	if(fESD)
	529	track = fESD->GetTrack(idxRec);
	530	else if(fAOD)
	531	track = fAOD->GetTrack(idxRec);
	532
	533	if (track) {
1063a5fa	534	etaRec = track->Eta();
5f071dcf	535	phiRec = track->Phi();
	536	ptRec = track->Pt();
	537	signRec = track->Charge();
1063a5fa	538	fHelper->IsTrackAcceptedRapidity(track, yRec, ipid);
	539	Double_t yeta = (ipid == 0 ) ? etaRec : yRec;
	540	Double_t hneffRec[5] = {fCentralityBin,signRec, yeta,phiRec,ptRec};
	541	if (ipid == 0) { fHnNchERec->Fill(hneffRec); }
	542	else if (ipid == 1) { fHnNpiERec->Fill(hneffRec); }
	543	else if (ipid == 2) { fHnNkaERec->Fill(hneffRec); }
	544	else if (ipid == 3) { fHnNprERec->Fill(hneffRec); }
	545
d61393df	546	}
5f071dcf	547	break;
	548	}
	549	} // for (Int_t idxRec=0; idxRec < fNTracks; ++idxRec) {
	550	/*
	551	Double_t deltaPhi = particle->Phi()-phiRec;
	552	if (TMath::Abs(deltaPhi) > TMath::TwoPi()) {
	553	if (deltaPhi < 0)
	554	deltaPhi += TMath::TwoPi();
	555	else
	556	deltaPhi -= TMath::TwoPi();
	557	}
	558	*/
	559
	560	Double_t yetapid = (ipid == 0 ) ? particle->Eta() : particle->Y();
1063a5fa	561
5f071dcf	562
d61393df	563	// Printf("%2d %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f", ipid, yetapid, yeta, particle->Eta(), particle->Y(), etaRec, yRec);
5f071dcf	564
5f071dcf	565	Double_t hneffMc[8] = {fCentralityBin,signMC,findable, recStatus, recPid,yetapid,particle->Phi(),particle->Pt()};
1063a5fa	566
5f071dcf	567
1063a5fa	568	if (ipid == 0) { fHnNchEMc->Fill(hneffMc); }
	569	else if (ipid == 1) { fHnNpiEMc->Fill(hneffMc); }
	570	else if (ipid == 2) { fHnNkaEMc->Fill(hneffMc); }
	571	else if (ipid == 3) { fHnNprEMc->Fill(hneffMc); }
5f071dcf	572
	573
	574	} // for (Int_t idxMC = 0; idxMC < nPart; ++idxMC) {
	575
	576	return;
	577	}