[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 "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)
  , fElectronProperties(NULL)
  , fWhichCut(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;
  }

  // NOTE: external objects fPID and fCuts are not deleted here
  fPID=NULL;
  fCuts=NULL;
}

const char* AliDxHFEParticleSelectionEl::fgkCutBinNames[]={
  "kRecKineITSTPC",
  "kRecPrim",
  "kHFEcuts",
  "kHFEcutsTOFPID",
  "kHFEcutsTPCPID",
  "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);

  //
  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);

  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);

  //--------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;


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

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

  if(fPID && fPID->IsSelected(&hfetrack)) {
    AliDebug(3,"Inside FilldPhi, electron is selected");
    fWhichCut->Fill(kSelected);
    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;
  }
}

//________________________________________________________________________
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"
	42	#include "TAxis.h"
72c0a987	43	#include "TObjArray.h"
9535cec9	44	#include <iostream>
	45	#include <cerrno>
	46	#include <memory>
	47
	48	using namespace std;
72c0a987	49
	50	/// ROOT macro for the implementation of ROOT specific class methods
	51	ClassImp(AliDxHFEParticleSelectionEl)
	52
	53	AliDxHFEParticleSelectionEl::AliDxHFEParticleSelectionEl(const char* opt)
9535cec9	54	: AliDxHFEParticleSelection("electron", opt)
	55	, fPID(NULL)
	56	, fElectronProperties(NULL)
	57	, fWhichCut(NULL)
	58	, fCuts(NULL)
	59	, fCFM(NULL)
72c0a987	60	{
	61	// constructor
	62	//
	63	//
	64	//
	65	//
	66	}
	67
	68	AliDxHFEParticleSelectionEl::~AliDxHFEParticleSelectionEl()
	69	{
	70	// destructor
9535cec9	71	if (fElectronProperties) {
	72	delete fElectronProperties;
	73	fElectronProperties=NULL;
	74	}
	75	if(fCFM){
	76	delete fCFM;
	77	fCFM=NULL;
	78	}
d731501a	79	if(fWhichCut){
	80	delete fWhichCut;
	81	fWhichCut=NULL;
	82	}
9535cec9	83
	84	// NOTE: external objects fPID and fCuts are not deleted here
	85	fPID=NULL;
	86	fCuts=NULL;
	87	}
9535cec9	88
d731501a	89	const char* AliDxHFEParticleSelectionEl::fgkCutBinNames[]={
	90	"kRecKineITSTPC",
	91	"kRecPrim",
	92	"kHFEcuts",
	93	"kHFEcutsTOFPID",
	94	"kHFEcutsTPCPID",
	95	"kPID",
	96	"Selected e"
	97	};
9535cec9	98
d731501a	99	int AliDxHFEParticleSelectionEl::Init()
	100	{
	101	//
	102	// init function
	103	//
	104	int iResult=0;
9535cec9	105
d731501a	106	// Implicit call to InitControlObjects() before setting up CFM and fCuts
	107	// (if not there)
	108	iResult=AliDxHFEParticleSelection::Init();
	109	if (iResult<0) return iResult;
9535cec9	110
dfe96b90	111	//--------Initialize correction Framework and Cuts-------------------------
9535cec9	112	// Consider moving this, either to separate function or
	113	// add a set function for AliCFManager
	114	// Do we need this? Can we just call AliHFEcuts::CheckParticleCuts
	115	AliInfo("Setting up CFM");
	116	fCFM = new AliCFManager;
	117	// the setup of cut objects is done in AliHFEcuts::Initialize
	118	// the ids used by this class must be the same, the code might be broken if
	119	// the sequence in AliHFEcuts::Initialize is changed
	120	const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
	121	// reset pointers in the CF manager
	122	fCFM->SetNStepParticle(kNcutSteps);
	123	for(Int_t istep = 0; istep < kNcutSteps; istep++) {
	124	fCFM->SetParticleCutsList(istep, NULL);
	125	}
	126	if(!fCuts) {
	127	AliWarning("Cuts not available. Default cuts will be used");
	128	fCuts = new AliHFEcuts;
	129	fCuts->CreateStandardCuts();
	130	}
	131	// TODO: error handling?
	132	fCuts->Initialize(fCFM);
	133
	134	return 0;
d731501a	135
	136	}
	137
dcf83226	138	THnSparse* AliDxHFEParticleSelectionEl::DefineTHnSparse()
d731501a	139	{
	140	//
	141	// Defines the THnSparse. For now, only calls CreatControlTHnSparse
dcf83226	142
dcf83226	143	// here is the only place to change the dimension
d731501a	144	const int thnSize = 3;
dcf83226	145	InitTHnSparseArray(thnSize);
d731501a	146	const double Pi=TMath::Pi();
	147	TString name;
	148	name.Form("%s info", GetName());
	149
dcf83226	150	// 0 1 2
	151	// Pt Phi Eta
	152	int thnBins [thnSize] = { 1000, 200, 500};
	153	double thnMin [thnSize] = { 0, 0, -1.};
	154	double thnMax [thnSize] = { 100, 2*Pi, 1.};
	155	const char* thnNames[thnSize] = { "Pt","Phi","Eta"};
d731501a	156
dcf83226	157	return CreateControlTHnSparse(name,thnSize,thnBins,thnMin,thnMax,thnNames);
d731501a	158	}
	159
	160	int AliDxHFEParticleSelectionEl::InitControlObjects()
	161	{
	162	/// init control and monitoring objects
	163	AliInfo("Electron THnSparse");
	164
	165	fElectronProperties=DefineTHnSparse();
	166	AddControlObject(fElectronProperties);
	167
	168	//
	169	fWhichCut= new TH1F("fWhichCut","effective cut for a rejected particle",kNCutLabels,-0.5,kNCutLabels-0.5);
	170	for (int iLabel=0; iLabel<kNCutLabels; iLabel++)
	171	fWhichCut->GetXaxis()->SetBinLabel(iLabel+1, fgkCutBinNames[iLabel]);
	172	AddControlObject(fWhichCut);
	173
	174	return AliDxHFEParticleSelection::InitControlObjects();
9535cec9	175	}
	176
	177	int AliDxHFEParticleSelectionEl::HistogramParticleProperties(AliVParticle* p, int selectionCode)
	178	{
	179	/// histogram particle properties
	180	if (!p) return -EINVAL;
	181	//if (!fControlObjects) return 0;
	182	if(selectionCode==0) return 0;
dcf83226	183
	184	if(fElectronProperties && ParticleProperties()) {
	185	FillParticleProperties(p, ParticleProperties(), GetDimTHnSparse());
	186	fElectronProperties->Fill(ParticleProperties());
	187	}
	188
9535cec9	189	return 0;
	190	}
	191
dcf83226	192	int AliDxHFEParticleSelectionEl::FillParticleProperties(AliVParticle* p, Double_t* data, int dimension) const
d731501a	193	{
	194	// fill the data array from the particle data
	195	if (!data) return -EINVAL;
	196	AliAODTrack track=(AliAODTrack)p;
	197	if (!track) return -ENODATA;
	198	int i=0;
dcf83226	199	if (dimension!=GetDimTHnSparse()) {
d731501a	200	// TODO: think about filling only the available data and throwing a warning
	201	return -ENOSPC;
	202	}
	203	data[i++]=track->Pt();
	204	data[i++]=track->Phi();
	205	data[i++]=track->Eta();
	206	return i;
	207	}
	208
	209
b32d523b	210	int AliDxHFEParticleSelectionEl::IsSelected(AliVParticle* pEl, const AliVEvent* pEvent)
9535cec9	211	{
9535cec9	212	/// select El candidates
d731501a	213	// TODO: How to handle MC? would be too much duplicate code if copy entire IsSelected.
d731501a	214
9535cec9	215	AliAODTrack track=(AliAODTrack)pEl;
b32d523b	216	fCFM->SetRecEventInfo(pEvent);
9535cec9	217
	218	//--------track cut selection-----------------------
	219	//Using AliHFECuts:
	220	// RecKine: ITSTPC cuts
	221	if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)){
d731501a	222	fWhichCut->Fill(kRecKineITSTPC);
9535cec9	223	return 0;
	224	}
	225
	226	// RecPrim
	227	if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) {
d731501a	228	fWhichCut->Fill(kRecPrim);
9535cec9	229	return 0;
	230	}
	231
	232	// HFEcuts: ITS layers cuts
	233	if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) {
d731501a	234	fWhichCut->Fill(kHFEcutsITS);
9535cec9	235	return 0;
	236	}
	237
	238	// HFE cuts: TOF PID and mismatch flag
	239	if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTOF, track)) {
d731501a	240	fWhichCut->Fill(kHFEcutsTOF);
9535cec9	241	return 0;
	242	}
	243
	244	// HFE cuts: TPC PID cleanup
	245	if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)){
d731501a	246	fWhichCut->Fill(kHFEcutsTPC);
9535cec9	247	return 0;
	248	}
	249
	250	// HFEcuts: Nb of tracklets TRD0
	251	//if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTRD, track)) continue;
	252
	253
	254	//--------PID selection-----------------------
	255	AliHFEpidObject hfetrack;
	256	hfetrack.SetAnalysisType(AliHFEpidObject::kAODanalysis);
	257	hfetrack.SetRecTrack(track);
	258
	259	// TODO: configurable colliding system
	260	//if(IsPbPb()) hfetrack.SetPbPb();
	261	hfetrack.SetPP();
	262
	263	if(fPID && fPID->IsSelected(&hfetrack)) {
	264	AliDebug(3,"Inside FilldPhi, electron is selected");
d731501a	265	fWhichCut->Fill(kSelected);
9535cec9	266	return 1;
	267	}
	268	else{
d731501a	269	fWhichCut->Fill(kPID);
9535cec9	270	return 0;
	271	}
	272	}
	273
	274	void AliDxHFEParticleSelectionEl::SetCuts(TObject* cuts, int level)
	275	{
	276	/// set cut objects
	277	if (level==kCutHFE) {
	278	fCuts=dynamic_cast<AliHFEcuts*>(cuts);
	279	if (!fCuts && cuts) {
	280	AliError(Form("Cut object is not of required type AliHFEcuts but %s", cuts->ClassName()));
	281	}
	282	return;
	283	}
	284
	285	if (level==kCutPID) {
	286	fPID=dynamic_cast<AliHFEpid*>(cuts);
	287	if (!fPID && cuts) {
	288	AliError(Form("cuts object is not of required type AliHFEpid but %s", cuts->ClassName()));
	289	}
	290	return;
	291	}
72c0a987	292	}
72c0a987	293
9535cec9	294	//________________________________________________________________________
9535cec9	295	Bool_t AliDxHFEParticleSelectionEl::ProcessCutStep(Int_t cutStep, AliVParticle *track)
72c0a987	296	{
9535cec9	297	// Check single track cuts for a given cut step
	298	const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
	299	if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
	300	//if(!fCuts->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
	301	return kTRUE;
72c0a987	302	}