[u/mrichter/AliRoot.git] / PWGHF / correlationHF / AliDxHFEParticleSelectionEl.cxx

// $Id$

//**************************************************************************
//* This file is property of and copyright by the ALICE Project            * 
//* ALICE Experiment at CERN, All rights reserved.                         *
//*                                                                        *
//* Primary Authors: Matthias Richter <Matthias.Richter@ift.uib.no>        *
//*                  Sedat Altinpinar <Sedat.Altinpinar@cern.ch>           *
//*                  Hege Erdal       <hege.erdal@gmail.com>               *
//*                                                                        *
//* 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.                  *
//**************************************************************************

/// @file   AliDxHFEParticleSelectionEl.cxx
/// @author Sedat Altinpinar, Hege Erdal, Matthias Richter
/// @date   2012-03-19
/// @brief  D0 selection for D0-HFE correlation
///
#include "AliDxHFEParticleSelectionEl.h"
#include "AliVParticle.h"
#include "AliVEvent.h"
#include "AliPID.h"
#include "AliPIDResponse.h"
#include "AliHFEcontainer.h"
#include "AliHFEpid.h"
#include "AliHFEpidBase.h"
#include "AliHFEtools.h"
#include "AliHFEcuts.h"
#include "AliAODTrack.h"
#include "AliAnalysisDataSlot.h"
#include "AliAnalysisDataContainer.h"
#include "AliAnalysisManager.h"
#include "AliCFManager.h"
#include "THnSparse.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TAxis.h"
#include "TObjArray.h"
#include <iostream>
#include <cerrno>
#include <memory>

using namespace std;

/// ROOT macro for the implementation of ROOT specific class methods
ClassImp(AliDxHFEParticleSelectionEl)

AliDxHFEParticleSelectionEl::AliDxHFEParticleSelectionEl(const char* opt)
  : AliDxHFEParticleSelection("electron", opt)
  , fPID(NULL)
  , fPIDTOF(NULL)
  , fElectronProperties(NULL)
  , fWhichCut(NULL)
  , fdEdx(NULL)
  , fdEdxCut(NULL)
  , fdEdxPid(NULL)
  , fdEdxPidTOF(NULL)
  , fnSigTPC(NULL)
  , fnSigTPCCut(NULL)
  , fnSigTPCPid(NULL)
  , fnSigTPCPidTOF(NULL)
  , fnSigTOF(NULL)
  , fnSigTOFCut(NULL)
  , fnSigTOFPid(NULL)
  , fnSigTOFPidTOF(NULL)
  , fPIDResponse(NULL)
  , fCuts(NULL)
  , fCFM(NULL)
{
  // constructor
  // 
  // 
  // 
  // 
}

AliDxHFEParticleSelectionEl::~AliDxHFEParticleSelectionEl()
{
  // destructor
  if (fElectronProperties) {
    delete fElectronProperties;
    fElectronProperties=NULL;
  }
  if(fCFM){
    delete fCFM;
    fCFM=NULL;
  }
  if(fWhichCut){
    delete fWhichCut;
    fWhichCut=NULL;
  }
  if(fdEdx){
    delete fdEdx;
    fdEdx=NULL;
  }
  if(fdEdxCut){
    delete fdEdxCut;
    fdEdxCut=NULL;
  }
  if(fdEdxPid){
    delete fdEdxPid;
    fdEdxPid=NULL;
  }
  if(fdEdxPidTOF){
    delete fdEdxPidTOF;
    fdEdxPidTOF=NULL;
  }
  if(fnSigTPC){
    delete fnSigTPC;
    fnSigTPC=NULL;
  }
  if(fnSigTPCCut){
    delete fnSigTPCCut;
    fnSigTPCCut=NULL;
  }
  if(fnSigTPCPid){
    delete fnSigTPCPid;
    fnSigTPCPid=NULL;
  }
  if(fnSigTPCPidTOF){
    delete fnSigTPCPidTOF;
    fnSigTPCPidTOF=NULL;
  }
  if(fnSigTOF){
    delete fnSigTOF;
    fnSigTOF=NULL;
  }
  if(fnSigTOFCut){
    delete fnSigTOFCut;
    fnSigTOFCut=NULL;
  }
  if(fnSigTOFPid){
    delete fnSigTOFPid;
    fnSigTOFPid=NULL;
  }

  if(fnSigTOFPidTOF){
    delete fnSigTOFPidTOF;
    fnSigTOFPidTOF=NULL;
  }

  if(fPIDResponse){
    delete fPIDResponse;
    fPIDResponse=NULL;
  }
  // NOTE: external objects fPID and fCuts are not deleted here
  fPID=NULL;
  fPIDTOF=NULL;
  fCuts=NULL;
}

const char* AliDxHFEParticleSelectionEl::fgkCutBinNames[]={
  "kRecKineITSTPC",
  "kRecPrim",
  "kHFEcuts",
  "kHFEcutsTOFPID",
  "kHFEcutsTPCPID",
  "kPIDTOF",
  "kPID",
  "Selected e"

};

int AliDxHFEParticleSelectionEl::Init()
{
  //
  // init function
  // 
  int iResult=0;

  // Implicit call to InitControlObjects() before setting up CFM and fCuts
  // (if not there)
  iResult=AliDxHFEParticleSelection::Init();
  if (iResult<0) return iResult;

  //--------Initialize correction Framework and Cuts-------------------------
  // Consider moving this, either to separate function or
  // add a set function for AliCFManager
  // Do we need this? Can we just call AliHFEcuts::CheckParticleCuts
  AliInfo("Setting up CFM");
  fCFM = new AliCFManager;
  // the setup of cut objects is done in AliHFEcuts::Initialize
  // the ids used by this class must be the same, the code might be broken if
  // the sequence in AliHFEcuts::Initialize is changed
  const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
  // reset pointers in the CF manager
  fCFM->SetNStepParticle(kNcutSteps);
  for(Int_t istep = 0; istep < kNcutSteps; istep++) {
    fCFM->SetParticleCutsList(istep, NULL);
  }
  if(!fCuts) {
    AliWarning("Cuts not available. Default cuts will be used");
    fCuts = new AliHFEcuts;
    fCuts->CreateStandardCuts();
  }
  // TODO: error handling?
  fCuts->Initialize(fCFM);

  return 0;

}

THnSparse* AliDxHFEParticleSelectionEl::DefineTHnSparse()
{
  //
  // Defines the THnSparse. For now, only calls CreatControlTHnSparse

  // here is the only place to change the dimension
  const int thnSize = 3;
  InitTHnSparseArray(thnSize);
  const double Pi=TMath::Pi();
  TString name;
  name.Form("%s info", GetName());

  //                                 0    1       2
  // 	 	                     Pt   Phi    Eta
  int         thnBins [thnSize] = { 1000,  200, 500};
  double      thnMin  [thnSize] = {    0,    0, -1.};
  double      thnMax  [thnSize] = {  100, 2*Pi,  1.};
  const char* thnNames[thnSize] = { "Pt","Phi","Eta"};

  return CreateControlTHnSparse(name,thnSize,thnBins,thnMin,thnMax,thnNames);
}

int AliDxHFEParticleSelectionEl::InitControlObjects()
{
  /// init control and monitoring objects
  AliInfo("Electron THnSparse");

  fElectronProperties=DefineTHnSparse();
  AddControlObject(fElectronProperties);

  //Look into VarManager to set up and fill histos
  fWhichCut= new TH1F("fWhichCut","effective cut for a rejected particle",kNCutLabels,-0.5,kNCutLabels-0.5);
  for (int iLabel=0; iLabel<kNCutLabels; iLabel++)
    fWhichCut->GetXaxis()->SetBinLabel(iLabel+1, fgkCutBinNames[iLabel]);
  AddControlObject(fWhichCut);
  //
  // dEdx plots, TPC signal vs momentum
  // 
  fdEdx = new TH2F("fdEdx", "dEdx before cuts", 1000,0.,10.,200,0.,200.);
  fdEdx->GetXaxis()->SetTitle("momentum (GeV/c)");
  fdEdx->GetYaxis()->SetTitle("dE/dx in TPC (a.u.)");
  AddControlObject(fdEdx);
  //
  fdEdxCut = new TH2F("fdEdxCut", "dEdx after cuts", 1000,0.,10.,200,0.,200.);
  fdEdxCut->GetXaxis()->SetTitle("momentum (GeV/c)");
  fdEdxCut->GetYaxis()->SetTitle("dE/dx in TPC (a.u.)");
  AddControlObject(fdEdxCut);
  //
  fdEdxPid = new TH2F("fdEdxPid", "dEdx after pid",  1000,0.,10.,200,0.,200.);
  fdEdxPid->GetXaxis()->SetTitle("momentum (GeV/c)");
  fdEdxPid->GetYaxis()->SetTitle("dE/dx in TPC (a.u.)");
  AddControlObject(fdEdxPid);
  //
  fdEdxPidTOF = new TH2F("fdEdxPidTOF", "dEdx after TOF pid",  1000,0.,10.,200,0.,200.);
  fdEdxPidTOF->GetXaxis()->SetTitle("momentum (GeV/c)");
  fdEdxPidTOF->GetYaxis()->SetTitle("dE/dx in TPC (a.u.)");
  AddControlObject(fdEdxPidTOF);

  //
  // nSigmaTPC vs momentum
  //

  fnSigTPC = new TH2F("fnSigTPC", "nSigTPC before cuts",  1000,0.,10.,200,-10.,10.);
  fnSigTPC->GetXaxis()->SetTitle("momentum (GeV/c)");
  fnSigTPC->GetYaxis()->SetTitle("nSigma in TPC (a.u.)");
  AddControlObject(fnSigTPC);


  fnSigTPCCut = new TH2F("fnSigTPCCut", "nSigmaTPC after cuts",  1000,0.,10.,200,-10.,10.);
  fnSigTPCCut->GetXaxis()->SetTitle("momentum (GeV/c)");
  fnSigTPCCut->GetYaxis()->SetTitle("nSigma in TPC (a.u.)");
  AddControlObject(fnSigTPCCut);


  fnSigTPCPid = new TH2F("fnSigTPCPid", "nSigmaTPC after PID",  1000,0.,10.,200,-10.,10.);
  fnSigTPCPid->GetXaxis()->SetTitle("momentum (GeV/c)");
  fnSigTPCPid->GetYaxis()->SetTitle("nSigma in TPC (a.u.)");
  AddControlObject(fnSigTPCPid);


  fnSigTPCPidTOF = new TH2F("fnSigTPCPidTOF", "nSigmaTPC after TOF PID",  1000,0.,10.,200,-10.,10.);
  fnSigTPCPidTOF->GetXaxis()->SetTitle("momentum (GeV/c)");
  fnSigTPCPidTOF->GetYaxis()->SetTitle("nSigma in TPC (a.u.)");
  AddControlObject(fnSigTPCPidTOF);

  //
  // nSigmaTOF vs momentum
  //

  fnSigTOF = new TH2F("fnSigTOF", "nSigmaTOF before cuts",  1000,0.,10.,200,-10.,10.);
  fnSigTOF->GetXaxis()->SetTitle("momentum (GeV/c)");
  fnSigTOF->GetYaxis()->SetTitle("nSigma in TOF (a.u.)");
  AddControlObject(fnSigTOF);


  fnSigTOFCut = new TH2F("fnSigTOFCut", "nSigmaTOF after cuts",  1000,0.,10.,200,-10.,10.);
  fnSigTOFCut->GetXaxis()->SetTitle("momentum (GeV/c)");
  fnSigTOFCut->GetYaxis()->SetTitle("nSigma in TOF (a.u.)");
  AddControlObject(fnSigTOFCut);


  fnSigTOFPid = new TH2F("fnSigTOFPid", "nSigmaTOF after PID",  1000,0.,10.,200,-10.,10.);
  fnSigTOFPid->GetXaxis()->SetTitle("momentum (GeV/c)");
  fnSigTOFPid->GetYaxis()->SetTitle("nSigma in TOF (a.u.)");
  AddControlObject(fnSigTOFPid);


  fnSigTOFPidTOF = new TH2F("fnSigTOFPidTOF", "nSigmaTOF after TOF PID",  1000,0.,10.,200,-10.,10.);
  fnSigTOFPidTOF->GetXaxis()->SetTitle("momentum (GeV/c)");
  fnSigTOFPidTOF->GetYaxis()->SetTitle("nSigma in TOF (a.u.)");
  AddControlObject(fnSigTOFPidTOF);


  return AliDxHFEParticleSelection::InitControlObjects();
}

int AliDxHFEParticleSelectionEl::HistogramParticleProperties(AliVParticle* p, int selectionCode)
{
  /// histogram particle properties
  if (!p) return -EINVAL;
  //if (!fControlObjects) return 0;
  if(selectionCode==0) return  0;
  
  if(fElectronProperties && ParticleProperties()) {
    FillParticleProperties(p, ParticleProperties(), GetDimTHnSparse());
    fElectronProperties->Fill(ParticleProperties());
  }

  return 0;
}

int AliDxHFEParticleSelectionEl::FillParticleProperties(AliVParticle* p, Double_t* data, int dimension) const
{
  // fill the data array from the particle data
  if (!data) return -EINVAL;
  AliAODTrack *track=(AliAODTrack*)p;
  if (!track) return -ENODATA;
  int i=0;
  if (dimension!=GetDimTHnSparse()) {
    // TODO: think about filling only the available data and throwing a warning
    return -ENOSPC;
  }
  data[i++]=track->Pt();
  data[i++]=track->Phi();
  data[i++]=track->Eta();
  return i;
}


int AliDxHFEParticleSelectionEl::IsSelected(AliVParticle* pEl, const AliVEvent* pEvent)
{
  /// select El candidates
  // TODO: How to handle MC? would be too much duplicate code if copy entire IsSelected. 

  AliAODTrack *track=(AliAODTrack*)pEl;
  fCFM->SetRecEventInfo(pEvent);
  
  fdEdx->Fill(track->GetTPCmomentum(), track->GetTPCsignal());
  fnSigTPC->Fill(track->GetTPCmomentum(), fPIDResponse->NumberOfSigmasTPC(track,AliPID::kElectron));
  fnSigTOF->Fill(track->P(), fPIDResponse->NumberOfSigmasTOF(track,AliPID::kElectron));
  
  //--------track cut selection-----------------------
  //Using AliHFECuts:
  // RecKine: ITSTPC cuts  
  if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)){
    fWhichCut->Fill(kRecKineITSTPC);
    return 0;
  }
  
  // RecPrim
  if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) {
    fWhichCut->Fill(kRecPrim);
    return 0;
  }
  
  // HFEcuts: ITS layers cuts
  if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) {
    fWhichCut->Fill(kHFEcutsITS);
    return 0;
  }
  
  // HFE cuts: TOF PID and mismatch flag
  if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTOF, track)) {
    fWhichCut->Fill(kHFEcutsTOF);
    return 0;
  }
  
  // HFE cuts: TPC PID cleanup
  if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)){
    fWhichCut->Fill(kHFEcutsTPC);
    return 0;
  } 
 
  // HFEcuts: Nb of tracklets TRD0
  //if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTRD, track)) continue;
  fdEdxCut->Fill(track->GetTPCmomentum(), track->GetTPCsignal());
  fnSigTPCCut->Fill(track->GetTPCmomentum(), fPIDResponse->NumberOfSigmasTPC(track,AliPID::kElectron));
  fnSigTOFCut->Fill(track->P(), fPIDResponse->NumberOfSigmasTOF(track,AliPID::kElectron));


  //--------PID selection-----------------------
  AliHFEpidObject hfetrack;
  hfetrack.SetAnalysisType(AliHFEpidObject::kAODanalysis);
  hfetrack.SetRecTrack(track);

  // TODO: configurable colliding system
  //if(IsPbPb()) hfetrack.SetPbPb();
  hfetrack.SetPP();

  if(fPIDTOF && fPIDTOF->IsSelected(&hfetrack)) {
    fdEdxPidTOF->Fill(track->GetTPCmomentum(), track->GetTPCsignal());
    fnSigTPCPidTOF->Fill(track->GetTPCmomentum(), fPIDResponse->NumberOfSigmasTPC(track,AliPID::kElectron));
    fnSigTOFPidTOF->Fill(track->P(), fPIDResponse->NumberOfSigmasTOF(track,AliPID::kElectron));
  }
  else{
    fWhichCut->Fill(kPIDTOF);
    return 0;
  }

  //Combined tof & tpc pid
  if(fPID && fPID->IsSelected(&hfetrack)) {
    AliDebug(3,"Inside FilldPhi, electron is selected");
    fWhichCut->Fill(kSelected);
    fdEdxPid->Fill(track->GetTPCmomentum(), track->GetTPCsignal());
    fnSigTPCPid->Fill(track->GetTPCmomentum(), fPIDResponse->NumberOfSigmasTPC(track,AliPID::kElectron));
    fnSigTOFPid->Fill(track->P(), fPIDResponse->NumberOfSigmasTOF(track,AliPID::kElectron));
    return 1;
  }
  else{
    fWhichCut->Fill(kPID);
    return 0;
  }


}

void AliDxHFEParticleSelectionEl::SetCuts(TObject* cuts, int level)
{
  /// set cut objects
  if (level==kCutHFE) {
    fCuts=dynamic_cast<AliHFEcuts*>(cuts);
    if (!fCuts && cuts) {
      AliError(Form("Cut object is not of required type AliHFEcuts but %s", cuts->ClassName()));
    }
    return;
  }

  if (level==kCutPID) {
    fPID=dynamic_cast<AliHFEpid*>(cuts);
    if (!fPID && cuts) {
      AliError(Form("cuts object is not of required type AliHFEpid but %s", cuts->ClassName()));
    }
    return;
  }

  if (level==kCutPIDTOF) {
    fPIDTOF=dynamic_cast<AliHFEpid*>(cuts);
    if (!fPIDTOF && cuts) {
      AliError(Form("cuts object is not of required type AliHFEpid but %s", cuts->ClassName()));
    }
    return;
  }
}

//________________________________________________________________________
Bool_t AliDxHFEParticleSelectionEl::ProcessCutStep(Int_t cutStep, AliVParticle *track)
{
  // Check single track cuts for a given cut step
  const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
  if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
  //if(!fCuts->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
  return kTRUE;
}
Commit	Line	Data
72c0a987	1	// $Id$
	2
	3	//**************************************************************************
	4	//* This file is property of and copyright by the ALICE Project *
	5	//* ALICE Experiment at CERN, All rights reserved. *
	6	//* *
	7	//* Primary Authors: Matthias Richter <Matthias.Richter@ift.uib.no> *
	8	//* Sedat Altinpinar <Sedat.Altinpinar@cern.ch> *
	9	//* Hege Erdal <hege.erdal@gmail.com> *
	10	//* *
	11	//* Permission to use, copy, modify and distribute this software and its *
	12	//* documentation strictly for non-commercial purposes is hereby granted *
	13	//* without fee, provided that the above copyright notice appears in all *
	14	//* copies and that both the copyright notice and this permission notice *
	15	//* appear in the supporting documentation. The authors make no claims *
	16	//* about the suitability of this software for any purpose. It is *
	17	//* provided "as is" without express or implied warranty. *
	18	//**************************************************************************
	19
	20	/// @file AliDxHFEParticleSelectionEl.cxx
	21	/// @author Sedat Altinpinar, Hege Erdal, Matthias Richter
	22	/// @date 2012-03-19
	23	/// @brief D0 selection for D0-HFE correlation
	24	///
72c0a987	25	#include "AliDxHFEParticleSelectionEl.h"
72c0a987	26	#include "AliVParticle.h"
9535cec9	27	#include "AliVEvent.h"
	28	#include "AliPID.h"
	29	#include "AliPIDResponse.h"
	30	#include "AliHFEcontainer.h"
	31	#include "AliHFEpid.h"
	32	#include "AliHFEpidBase.h"
	33	#include "AliHFEtools.h"
	34	#include "AliHFEcuts.h"
	35	#include "AliAODTrack.h"
	36	#include "AliAnalysisDataSlot.h"
	37	#include "AliAnalysisDataContainer.h"
	38	#include "AliAnalysisManager.h"
	39	#include "AliCFManager.h"
	40	#include "THnSparse.h"
	41	#include "TH1F.h"
b86451e1	42	#include "TH2F.h"
9535cec9	43	#include "TAxis.h"
72c0a987	44	#include "TObjArray.h"
9535cec9	45	#include <iostream>
	46	#include <cerrno>
	47	#include <memory>
	48
	49	using namespace std;
72c0a987	50
	51	/// ROOT macro for the implementation of ROOT specific class methods
	52	ClassImp(AliDxHFEParticleSelectionEl)
	53
	54	AliDxHFEParticleSelectionEl::AliDxHFEParticleSelectionEl(const char* opt)
9535cec9	55	: AliDxHFEParticleSelection("electron", opt)
9535cec9	56	, fPID(NULL)
b86451e1	57	, fPIDTOF(NULL)
9535cec9	58	, fElectronProperties(NULL)
b86451e1	59	, fWhichCut(NULL)
	60	, fdEdx(NULL)
	61	, fdEdxCut(NULL)
	62	, fdEdxPid(NULL)
	63	, fdEdxPidTOF(NULL)
	64	, fnSigTPC(NULL)
	65	, fnSigTPCCut(NULL)
	66	, fnSigTPCPid(NULL)
	67	, fnSigTPCPidTOF(NULL)
	68	, fnSigTOF(NULL)
	69	, fnSigTOFCut(NULL)
	70	, fnSigTOFPid(NULL)
	71	, fnSigTOFPidTOF(NULL)
	72	, fPIDResponse(NULL)
9535cec9	73	, fCuts(NULL)
9535cec9	74	, fCFM(NULL)
72c0a987	75	{
	76	// constructor
	77	//
	78	//
	79	//
	80	//
	81	}
	82
	83	AliDxHFEParticleSelectionEl::~AliDxHFEParticleSelectionEl()
	84	{
	85	// destructor
9535cec9	86	if (fElectronProperties) {
	87	delete fElectronProperties;
	88	fElectronProperties=NULL;
	89	}
	90	if(fCFM){
	91	delete fCFM;
	92	fCFM=NULL;
	93	}
d731501a	94	if(fWhichCut){
	95	delete fWhichCut;
	96	fWhichCut=NULL;
	97	}
b86451e1	98	if(fdEdx){
	99	delete fdEdx;
	100	fdEdx=NULL;
	101	}
	102	if(fdEdxCut){
	103	delete fdEdxCut;
	104	fdEdxCut=NULL;
	105	}
	106	if(fdEdxPid){
	107	delete fdEdxPid;
	108	fdEdxPid=NULL;
	109	}
	110	if(fdEdxPidTOF){
	111	delete fdEdxPidTOF;
	112	fdEdxPidTOF=NULL;
	113	}
	114	if(fnSigTPC){
	115	delete fnSigTPC;
	116	fnSigTPC=NULL;
	117	}
	118	if(fnSigTPCCut){
	119	delete fnSigTPCCut;
	120	fnSigTPCCut=NULL;
	121	}
	122	if(fnSigTPCPid){
	123	delete fnSigTPCPid;
	124	fnSigTPCPid=NULL;
	125	}
	126	if(fnSigTPCPidTOF){
	127	delete fnSigTPCPidTOF;
	128	fnSigTPCPidTOF=NULL;
	129	}
	130	if(fnSigTOF){
	131	delete fnSigTOF;
	132	fnSigTOF=NULL;
	133	}
	134	if(fnSigTOFCut){
	135	delete fnSigTOFCut;
	136	fnSigTOFCut=NULL;
	137	}
	138	if(fnSigTOFPid){
	139	delete fnSigTOFPid;
	140	fnSigTOFPid=NULL;
	141	}
9535cec9	142
b86451e1	143	if(fnSigTOFPidTOF){
	144	delete fnSigTOFPidTOF;
	145	fnSigTOFPidTOF=NULL;
	146	}
	147
	148	if(fPIDResponse){
	149	delete fPIDResponse;
	150	fPIDResponse=NULL;
	151	}
9535cec9	152	// NOTE: external objects fPID and fCuts are not deleted here
9535cec9	153	fPID=NULL;
b86451e1	154	fPIDTOF=NULL;
9535cec9	155	fCuts=NULL;
9535cec9	156	}
9535cec9	157
d731501a	158	const char* AliDxHFEParticleSelectionEl::fgkCutBinNames[]={
	159	"kRecKineITSTPC",
	160	"kRecPrim",
	161	"kHFEcuts",
	162	"kHFEcutsTOFPID",
	163	"kHFEcutsTPCPID",
b86451e1	164	"kPIDTOF",
d731501a	165	"kPID",
d731501a	166	"Selected e"
b86451e1	167
d731501a	168	};
9535cec9	169
d731501a	170	int AliDxHFEParticleSelectionEl::Init()
	171	{
	172	//
	173	// init function
	174	//
	175	int iResult=0;
9535cec9	176
d731501a	177	// Implicit call to InitControlObjects() before setting up CFM and fCuts
	178	// (if not there)
	179	iResult=AliDxHFEParticleSelection::Init();
	180	if (iResult<0) return iResult;
9535cec9	181
dfe96b90	182	//--------Initialize correction Framework and Cuts-------------------------
9535cec9	183	// Consider moving this, either to separate function or
	184	// add a set function for AliCFManager
	185	// Do we need this? Can we just call AliHFEcuts::CheckParticleCuts
	186	AliInfo("Setting up CFM");
	187	fCFM = new AliCFManager;
	188	// the setup of cut objects is done in AliHFEcuts::Initialize
	189	// the ids used by this class must be the same, the code might be broken if
	190	// the sequence in AliHFEcuts::Initialize is changed
	191	const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
	192	// reset pointers in the CF manager
	193	fCFM->SetNStepParticle(kNcutSteps);
	194	for(Int_t istep = 0; istep < kNcutSteps; istep++) {
	195	fCFM->SetParticleCutsList(istep, NULL);
	196	}
	197	if(!fCuts) {
	198	AliWarning("Cuts not available. Default cuts will be used");
	199	fCuts = new AliHFEcuts;
	200	fCuts->CreateStandardCuts();
	201	}
	202	// TODO: error handling?
	203	fCuts->Initialize(fCFM);
	204
	205	return 0;
d731501a	206
	207	}
	208
dcf83226	209	THnSparse* AliDxHFEParticleSelectionEl::DefineTHnSparse()
d731501a	210	{
	211	//
	212	// Defines the THnSparse. For now, only calls CreatControlTHnSparse
dcf83226	213
dcf83226	214	// here is the only place to change the dimension
d731501a	215	const int thnSize = 3;
dcf83226	216	InitTHnSparseArray(thnSize);
d731501a	217	const double Pi=TMath::Pi();
	218	TString name;
	219	name.Form("%s info", GetName());
	220
dcf83226	221	// 0 1 2
	222	// Pt Phi Eta
	223	int thnBins [thnSize] = { 1000, 200, 500};
	224	double thnMin [thnSize] = { 0, 0, -1.};
	225	double thnMax [thnSize] = { 100, 2*Pi, 1.};
	226	const char* thnNames[thnSize] = { "Pt","Phi","Eta"};
d731501a	227
dcf83226	228	return CreateControlTHnSparse(name,thnSize,thnBins,thnMin,thnMax,thnNames);
d731501a	229	}
	230
	231	int AliDxHFEParticleSelectionEl::InitControlObjects()
	232	{
	233	/// init control and monitoring objects
	234	AliInfo("Electron THnSparse");
	235
	236	fElectronProperties=DefineTHnSparse();
	237	AddControlObject(fElectronProperties);
	238
b86451e1	239	//Look into VarManager to set up and fill histos
d731501a	240	fWhichCut= new TH1F("fWhichCut","effective cut for a rejected particle",kNCutLabels,-0.5,kNCutLabels-0.5);
	241	for (int iLabel=0; iLabel<kNCutLabels; iLabel++)
	242	fWhichCut->GetXaxis()->SetBinLabel(iLabel+1, fgkCutBinNames[iLabel]);
	243	AddControlObject(fWhichCut);
b86451e1	244	//
	245	// dEdx plots, TPC signal vs momentum
	246	//
	247	fdEdx = new TH2F("fdEdx", "dEdx before cuts", 1000,0.,10.,200,0.,200.);
	248	fdEdx->GetXaxis()->SetTitle("momentum (GeV/c)");
	249	fdEdx->GetYaxis()->SetTitle("dE/dx in TPC (a.u.)");
	250	AddControlObject(fdEdx);
	251	//
	252	fdEdxCut = new TH2F("fdEdxCut", "dEdx after cuts", 1000,0.,10.,200,0.,200.);
	253	fdEdxCut->GetXaxis()->SetTitle("momentum (GeV/c)");
	254	fdEdxCut->GetYaxis()->SetTitle("dE/dx in TPC (a.u.)");
	255	AddControlObject(fdEdxCut);
	256	//
	257	fdEdxPid = new TH2F("fdEdxPid", "dEdx after pid", 1000,0.,10.,200,0.,200.);
	258	fdEdxPid->GetXaxis()->SetTitle("momentum (GeV/c)");
	259	fdEdxPid->GetYaxis()->SetTitle("dE/dx in TPC (a.u.)");
	260	AddControlObject(fdEdxPid);
	261	//
	262	fdEdxPidTOF = new TH2F("fdEdxPidTOF", "dEdx after TOF pid", 1000,0.,10.,200,0.,200.);
	263	fdEdxPidTOF->GetXaxis()->SetTitle("momentum (GeV/c)");
	264	fdEdxPidTOF->GetYaxis()->SetTitle("dE/dx in TPC (a.u.)");
	265	AddControlObject(fdEdxPidTOF);
	266
	267	//
	268	// nSigmaTPC vs momentum
	269	//
	270
	271	fnSigTPC = new TH2F("fnSigTPC", "nSigTPC before cuts", 1000,0.,10.,200,-10.,10.);
	272	fnSigTPC->GetXaxis()->SetTitle("momentum (GeV/c)");
	273	fnSigTPC->GetYaxis()->SetTitle("nSigma in TPC (a.u.)");
	274	AddControlObject(fnSigTPC);
	275
	276
	277	fnSigTPCCut = new TH2F("fnSigTPCCut", "nSigmaTPC after cuts", 1000,0.,10.,200,-10.,10.);
	278	fnSigTPCCut->GetXaxis()->SetTitle("momentum (GeV/c)");
	279	fnSigTPCCut->GetYaxis()->SetTitle("nSigma in TPC (a.u.)");
	280	AddControlObject(fnSigTPCCut);
	281
	282
	283	fnSigTPCPid = new TH2F("fnSigTPCPid", "nSigmaTPC after PID", 1000,0.,10.,200,-10.,10.);
	284	fnSigTPCPid->GetXaxis()->SetTitle("momentum (GeV/c)");
	285	fnSigTPCPid->GetYaxis()->SetTitle("nSigma in TPC (a.u.)");
	286	AddControlObject(fnSigTPCPid);
	287
	288
	289	fnSigTPCPidTOF = new TH2F("fnSigTPCPidTOF", "nSigmaTPC after TOF PID", 1000,0.,10.,200,-10.,10.);
	290	fnSigTPCPidTOF->GetXaxis()->SetTitle("momentum (GeV/c)");
	291	fnSigTPCPidTOF->GetYaxis()->SetTitle("nSigma in TPC (a.u.)");
	292	AddControlObject(fnSigTPCPidTOF);
	293
	294	//
	295	// nSigmaTOF vs momentum
	296	//
	297
	298	fnSigTOF = new TH2F("fnSigTOF", "nSigmaTOF before cuts", 1000,0.,10.,200,-10.,10.);
	299	fnSigTOF->GetXaxis()->SetTitle("momentum (GeV/c)");
	300	fnSigTOF->GetYaxis()->SetTitle("nSigma in TOF (a.u.)");
	301	AddControlObject(fnSigTOF);
	302
	303
	304	fnSigTOFCut = new TH2F("fnSigTOFCut", "nSigmaTOF after cuts", 1000,0.,10.,200,-10.,10.);
	305	fnSigTOFCut->GetXaxis()->SetTitle("momentum (GeV/c)");
	306	fnSigTOFCut->GetYaxis()->SetTitle("nSigma in TOF (a.u.)");
	307	AddControlObject(fnSigTOFCut);
308
309
310	fnSigTOFPid = new TH2F("fnSigTOFPid", "nSigmaTOF after PID", 1000,0.,10.,200,-10.,10.);
311	fnSigTOFPid->GetXaxis()->SetTitle("momentum (GeV/c)");
312	fnSigTOFPid->GetYaxis()->SetTitle("nSigma in TOF (a.u.)");
313	AddControlObject(fnSigTOFPid);
314
315
316	fnSigTOFPidTOF = new TH2F("fnSigTOFPidTOF", "nSigmaTOF after TOF PID", 1000,0.,10.,200,-10.,10.);
317	fnSigTOFPidTOF->GetXaxis()->SetTitle("momentum (GeV/c)");
318	fnSigTOFPidTOF->GetYaxis()->SetTitle("nSigma in TOF (a.u.)");
319	AddControlObject(fnSigTOFPidTOF);
320
d731501a	321
d731501a	322	return AliDxHFEParticleSelection::InitControlObjects();
9535cec9	323	}
	324
	325	int AliDxHFEParticleSelectionEl::HistogramParticleProperties(AliVParticle* p, int selectionCode)
	326	{
	327	/// histogram particle properties
	328	if (!p) return -EINVAL;
	329	//if (!fControlObjects) return 0;
	330	if(selectionCode==0) return 0;
dcf83226	331
	332	if(fElectronProperties && ParticleProperties()) {
	333	FillParticleProperties(p, ParticleProperties(), GetDimTHnSparse());
	334	fElectronProperties->Fill(ParticleProperties());
	335	}
	336
9535cec9	337	return 0;
	338	}
	339
dcf83226	340	int AliDxHFEParticleSelectionEl::FillParticleProperties(AliVParticle* p, Double_t* data, int dimension) const
d731501a	341	{
	342	// fill the data array from the particle data
	343	if (!data) return -EINVAL;
	344	AliAODTrack track=(AliAODTrack)p;
	345	if (!track) return -ENODATA;
	346	int i=0;
dcf83226	347	if (dimension!=GetDimTHnSparse()) {
d731501a	348	// TODO: think about filling only the available data and throwing a warning
	349	return -ENOSPC;
	350	}
	351	data[i++]=track->Pt();
	352	data[i++]=track->Phi();
	353	data[i++]=track->Eta();
	354	return i;
	355	}
	356
	357
b32d523b	358	int AliDxHFEParticleSelectionEl::IsSelected(AliVParticle* pEl, const AliVEvent* pEvent)
9535cec9	359	{
9535cec9	360	/// select El candidates
d731501a	361	// TODO: How to handle MC? would be too much duplicate code if copy entire IsSelected.
d731501a	362
9535cec9	363	AliAODTrack track=(AliAODTrack)pEl;
b32d523b	364	fCFM->SetRecEventInfo(pEvent);
b86451e1	365
	366	fdEdx->Fill(track->GetTPCmomentum(), track->GetTPCsignal());
	367	fnSigTPC->Fill(track->GetTPCmomentum(), fPIDResponse->NumberOfSigmasTPC(track,AliPID::kElectron));
	368	fnSigTOF->Fill(track->P(), fPIDResponse->NumberOfSigmasTOF(track,AliPID::kElectron));
	369
9535cec9	370	//--------track cut selection-----------------------
	371	//Using AliHFECuts:
	372	// RecKine: ITSTPC cuts
	373	if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)){
d731501a	374	fWhichCut->Fill(kRecKineITSTPC);
9535cec9	375	return 0;
	376	}
	377
	378	// RecPrim
	379	if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) {
d731501a	380	fWhichCut->Fill(kRecPrim);
9535cec9	381	return 0;
	382	}
	383
	384	// HFEcuts: ITS layers cuts
	385	if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) {
d731501a	386	fWhichCut->Fill(kHFEcutsITS);
9535cec9	387	return 0;
	388	}
	389
	390	// HFE cuts: TOF PID and mismatch flag
	391	if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTOF, track)) {
d731501a	392	fWhichCut->Fill(kHFEcutsTOF);
9535cec9	393	return 0;
	394	}
	395
	396	// HFE cuts: TPC PID cleanup
	397	if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)){
d731501a	398	fWhichCut->Fill(kHFEcutsTPC);
9535cec9	399	return 0;
	400	}
	401
	402	// HFEcuts: Nb of tracklets TRD0
	403	//if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTRD, track)) continue;
b86451e1	404	fdEdxCut->Fill(track->GetTPCmomentum(), track->GetTPCsignal());
	405	fnSigTPCCut->Fill(track->GetTPCmomentum(), fPIDResponse->NumberOfSigmasTPC(track,AliPID::kElectron));
	406	fnSigTOFCut->Fill(track->P(), fPIDResponse->NumberOfSigmasTOF(track,AliPID::kElectron));
9535cec9	407
	408
	409	//--------PID selection-----------------------
	410	AliHFEpidObject hfetrack;
	411	hfetrack.SetAnalysisType(AliHFEpidObject::kAODanalysis);
	412	hfetrack.SetRecTrack(track);
	413
	414	// TODO: configurable colliding system
	415	//if(IsPbPb()) hfetrack.SetPbPb();
	416	hfetrack.SetPP();
	417
b86451e1	418	if(fPIDTOF && fPIDTOF->IsSelected(&hfetrack)) {
	419	fdEdxPidTOF->Fill(track->GetTPCmomentum(), track->GetTPCsignal());
	420	fnSigTPCPidTOF->Fill(track->GetTPCmomentum(), fPIDResponse->NumberOfSigmasTPC(track,AliPID::kElectron));
	421	fnSigTOFPidTOF->Fill(track->P(), fPIDResponse->NumberOfSigmasTOF(track,AliPID::kElectron));
	422	}
	423	else{
	424	fWhichCut->Fill(kPIDTOF);
	425	return 0;
	426	}
	427
	428	//Combined tof & tpc pid
9535cec9	429	if(fPID && fPID->IsSelected(&hfetrack)) {
9535cec9	430	AliDebug(3,"Inside FilldPhi, electron is selected");
d731501a	431	fWhichCut->Fill(kSelected);
b86451e1	432	fdEdxPid->Fill(track->GetTPCmomentum(), track->GetTPCsignal());
	433	fnSigTPCPid->Fill(track->GetTPCmomentum(), fPIDResponse->NumberOfSigmasTPC(track,AliPID::kElectron));
	434	fnSigTOFPid->Fill(track->P(), fPIDResponse->NumberOfSigmasTOF(track,AliPID::kElectron));
9535cec9	435	return 1;
	436	}
	437	else{
d731501a	438	fWhichCut->Fill(kPID);
9535cec9	439	return 0;
9535cec9	440	}
b86451e1	441
b86451e1	442
9535cec9	443	}
	444
	445	void AliDxHFEParticleSelectionEl::SetCuts(TObject* cuts, int level)
	446	{
	447	/// set cut objects
	448	if (level==kCutHFE) {
	449	fCuts=dynamic_cast<AliHFEcuts*>(cuts);
	450	if (!fCuts && cuts) {
	451	AliError(Form("Cut object is not of required type AliHFEcuts but %s", cuts->ClassName()));
	452	}
	453	return;
	454	}
	455
	456	if (level==kCutPID) {
	457	fPID=dynamic_cast<AliHFEpid*>(cuts);
	458	if (!fPID && cuts) {
	459	AliError(Form("cuts object is not of required type AliHFEpid but %s", cuts->ClassName()));
	460	}
	461	return;
	462	}
b86451e1	463
	464	if (level==kCutPIDTOF) {
	465	fPIDTOF=dynamic_cast<AliHFEpid*>(cuts);
	466	if (!fPIDTOF && cuts) {
	467	AliError(Form("cuts object is not of required type AliHFEpid but %s", cuts->ClassName()));
	468	}
	469	return;
	470	}
72c0a987	471	}
72c0a987	472
9535cec9	473	//________________________________________________________________________
9535cec9	474	Bool_t AliDxHFEParticleSelectionEl::ProcessCutStep(Int_t cutStep, AliVParticle *track)
72c0a987	475	{
9535cec9	476	// Check single track cuts for a given cut step
	477	const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
	478	if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
	479	//if(!fCuts->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
	480	return kTRUE;
72c0a987	481	}