[u/mrichter/AliRoot.git] / PWG3 / hfe / AliHFEpidTRD.cxx

/**************************************************************************
* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
*                                                                        *
* Author: The ALICE Off-line Project.                                    *
* Contributors are mentioned in the code where appropriate.              *
*                                                                        *
* Permission to use, copy, modify and distribute this software and its   *
* documentation strictly for non-commercial purposes is hereby granted   *
* without fee, provided that the above copyright notice appears in all   *
* copies and that both the copyright notice and this permission notice   *
* appear in the supporting documentation. The authors make no claims     *
* about the suitability of this software for any purpose. It is          *
* provided "as is" without express or implied warranty.                  *
**************************************************************************/
/************************************************************************
 *                                                                      *
 * Class for TRD PID                                                    *
 * Implements the abstract base class AliHFEpidbase        *
 * Make PID does the PID decision                                       *
 * Class further contains TRD specific cuts and QA histograms           *
 *                                                                      *
 * Authors:                                                             *
 *   Markus Fasel <M.Fasel@gsi.de>                                      *
 *                                                                      *
 ************************************************************************/
#include <TAxis.h>
#include <TFile.h>
#include <TH1F.h>
#include <TIterator.h>
#include <TKey.h>
#include <TMap.h>
#include <TObjArray.h>
#include <TObjString.h>
#include <TString.h>
#include <TROOT.h>

#include "AliESDtrack.h"
#include "AliPID.h"

#include "AliHFEpidTRD.h"

ClassImp(AliHFEpidTRD)

//___________________________________________________________________
AliHFEpidTRD::AliHFEpidTRD(const char* name) :
    AliHFEpidBase(name)
  , fThresholdFile("TRD.PIDthresholds.root")
  , fPIDMethod(kNN)
  , fTRDthresholds(0x0)
  , fTRDelectronEfficiencies(0x0)
{
  //
  // default  constructor
  // 
  fTRDthresholds = new TMap();
}

//___________________________________________________________________
AliHFEpidTRD::AliHFEpidTRD(const AliHFEpidTRD &ref):
    AliHFEpidBase("")
  , fThresholdFile("")
  , fPIDMethod(kLQ)
  , fTRDthresholds(0x0)
  , fTRDelectronEfficiencies(0x0)
{
  //
  // Copy constructor
  //
  ref.Copy(*this);
}

//___________________________________________________________________
AliHFEpidTRD &AliHFEpidTRD::operator=(const AliHFEpidTRD &ref){
  //
  // Assignment operator
  //
  if(this != &ref){
    ref.Copy(*this);
  }
  return *this;
}

//___________________________________________________________________
void AliHFEpidTRD::Copy(TObject &ref) const {
  //
  // Performs the copying of the object
  //
  AliHFEpidTRD &target = dynamic_cast<AliHFEpidTRD &>(ref);

  target.fThresholdFile = fThresholdFile;
  target.fPIDMethod = fPIDMethod;
  if(fTRDthresholds){
    target.fTRDthresholds = dynamic_cast<TMap *>(fTRDthresholds->Clone());
  }
  if(fTRDelectronEfficiencies){
    target.fTRDelectronEfficiencies = new TAxis(*fTRDelectronEfficiencies);
  }
  AliHFEpidBase::Copy(ref);
}

//___________________________________________________________________
AliHFEpidTRD::~AliHFEpidTRD(){
  //
  // Destructor
  //
  if(fTRDthresholds){
    fTRDthresholds->Delete();
    delete fTRDthresholds;
  }
  if(fTRDelectronEfficiencies) delete fTRDelectronEfficiencies;
}

//______________________________________________________
Bool_t AliHFEpidTRD::InitializePID(){
  //
  // InitializePID: Load TRD thresholds and create the electron efficiency axis
  // to navigate 
  //
  LoadTRDthresholds();
  // Fill the electron efficiencies axis for the TRD alone PID
  Int_t nEffs = fTRDthresholds->GetEntries() + 1;
  TIterator* thres =fTRDthresholds->MakeIterator();
  TObjString *key = 0x0;
  Double_t *tmp = new Double_t[nEffs], *titer = tmp;
  while((key = dynamic_cast<TObjString *>((*thres)()))){
    (*titer++) = static_cast<Double_t>(key->String().Atoi())/100.;
  }
  delete thres;
  *titer = 1.;
  // Sort the electron efficiencies and put them into the TAxis for later navigation
  Int_t *ind = new Int_t[nEffs], *iiter = ind;
  TMath::Sort(nEffs, tmp, ind, kFALSE);
  Double_t *eleffs = new Double_t[nEffs], *eiter = eleffs;
  while(eiter < &eleffs[nEffs]) *(eiter++) = tmp[*(iiter++)];
  // print the content
  Int_t cnt = 0;
  if(GetDebugLevel() > 1){
    printf("Printing electron efficiency bins:\n");
    eiter = eleffs;
    thres = fTRDthresholds->MakeIterator();
    TObject *object = 0x0;
    while(eiter < &eleffs[nEffs - 1]){
      printf("eleffs[%d] = %f", cnt++, *(eiter++));
      key = dynamic_cast<TObjString *>(thres->Next());
      object = fTRDthresholds->GetValue(key->String().Data());
      printf(", Content: %p\n", (void *)object);
    }
  }
  delete[] tmp; delete[] ind;
  fTRDelectronEfficiencies = new TAxis(nEffs - 1, eleffs);
  if(GetDebugLevel() > 1){
    printf("Printing axis content:\n");
    for(Int_t ibin = fTRDelectronEfficiencies->GetFirst(); ibin <= fTRDelectronEfficiencies->GetLast(); ibin++)
      printf("%d.) minimum: %f, maximum? %f\n", ibin, fTRDelectronEfficiencies->GetBinLowEdge(ibin), fTRDelectronEfficiencies->GetBinUpEdge(ibin));
  }
  delete[] eleffs;
  return kTRUE;
}

//______________________________________________________
Int_t AliHFEpidTRD::IsSelected(AliVParticle *track){
  //
  // Does PID for TRD alone:
  // PID thresholds based on 90% Electron Efficiency level approximated by a linear 
  // step function
  //
  AliESDtrack *esdTrack = 0x0;
  if(!(esdTrack = dynamic_cast<AliESDtrack *>(track))) return kFALSE;
  Double_t p = esdTrack->GetOuterParam() ? esdTrack->GetOuterParam()->P() : esdTrack->P();
  if(p < 0.6) return 0;

  // Get the Histograms
  TH1 *threshist = GetTRDthresholds(0.91);
  Int_t bin = 0;
  if(p > threshist->GetXaxis()->GetXmax()) 
    bin = threshist->GetXaxis()->GetLast();
  else if(p < threshist->GetXaxis()->GetXmin()) 
    bin = threshist->GetXaxis()->GetFirst();
  else
    bin = threshist->GetXaxis()->FindBin(p);

  Double_t pidProbs[AliPID::kSPECIES];
  esdTrack->GetTRDpid(pidProbs);
  if(pidProbs[AliPID::kElectron] > threshist->GetBinContent(bin)) return 11;
  return 0;
}

//___________________________________________________________________
void AliHFEpidTRD::LoadTRDthresholds(){
  //
  // Load TRD threshold histograms from File
  //
  TFile *mythresholds = TFile::Open(fThresholdFile);
  TKey *object = 0x0;
  TString electron_eff;
  TObjArray *histos = 0x0;
  Float_t eff;
  TH1F *refhist = 0x0;
  TIterator *keyIterator = mythresholds->GetListOfKeys()->MakeIterator();
  TString histnames[2] = {"fHistThreshLQ", "fHistThreshNN"};
  gROOT->cd();
  while((object = dynamic_cast<TKey *>((*keyIterator)()))){
    // Get the electron efficiency bin this histogram was taken with
    electron_eff = object->GetName();
    electron_eff = electron_eff.Remove(0,3);
    eff = static_cast<Float_t>(electron_eff.Atoi())/100.;

    // Get the threshold according to the selected 
    histos = dynamic_cast<TObjArray *>(object->ReadObj());
    refhist = dynamic_cast<TH1F *>(histos->FindObject(histnames[fPIDMethod].Data()));
    SetTRDthresholds(refhist, eff);
    histos->Delete();
    delete histos;
  }
  delete keyIterator;
  mythresholds->Close();
  delete mythresholds;
}

//___________________________________________________________________
void AliHFEpidTRD::SetTRDthresholds(TH1F *thresholds, Float_t electronEff){
  //
  // Set the threshold histogram for the TRD pid
  //
  fTRDthresholds->Add(new TObjString(Form("%d", TMath::Nint(electronEff * 100.))), new TH1F(*thresholds));
}

//___________________________________________________________________
TH1F *AliHFEpidTRD::GetTRDthresholds(Float_t electronEff){ 
  Int_t bin = 0;
  if(electronEff < fTRDelectronEfficiencies->GetXmin()) bin = fTRDelectronEfficiencies->GetFirst();
  else if(electronEff > fTRDelectronEfficiencies->GetXmax()) bin = fTRDelectronEfficiencies->GetLast();
  else bin = fTRDelectronEfficiencies->FindBin(electronEff);
 TObjString keyname = Form("%d", TMath::Nint(fTRDelectronEfficiencies->GetBinLowEdge(bin)* 100.));
/*  printf("Key: %s\n", keyname.String().Data());*/
  TH1F *thresholds = dynamic_cast<TH1F *>((dynamic_cast<TPair *>(fTRDthresholds->FindObject(&keyname)))->Value());
/*  printf("thresholds: %p\n", thresholds);*/
  return thresholds;
}
Commit	Line	Data
809a4336	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15	/************************************************************************
	16	* *
	17	* Class for TRD PID *
	18	* Implements the abstract base class AliHFEpidbase *
	19	* Make PID does the PID decision *
	20	* Class further contains TRD specific cuts and QA histograms *
	21	* *
	22	* Authors: *
	23	* Markus Fasel <M.Fasel@gsi.de> *
	24	* *
	25	************************************************************************/
	26	#include <TAxis.h>
	27	#include <TFile.h>
	28	#include <TH1F.h>
	29	#include <TIterator.h>
	30	#include <TKey.h>
	31	#include <TMap.h>
	32	#include <TObjArray.h>
	33	#include <TObjString.h>
	34	#include <TString.h>
	35	#include <TROOT.h>
	36
	37	#include "AliESDtrack.h"
	38	#include "AliPID.h"
	39
	40	#include "AliHFEpidTRD.h"
	41
	42	ClassImp(AliHFEpidTRD)
	43
	44	//___________________________________________________________________
	45	AliHFEpidTRD::AliHFEpidTRD(const char* name) :
	46	AliHFEpidBase(name)
	47	, fThresholdFile("TRD.PIDthresholds.root")
	48	, fPIDMethod(kNN)
	49	, fTRDthresholds(0x0)
	50	, fTRDelectronEfficiencies(0x0)
	51	{
	52	//
	53	// default constructor
	54	//
	55	fTRDthresholds = new TMap();
	56	}
	57
	58	//___________________________________________________________________
	59	AliHFEpidTRD::AliHFEpidTRD(const AliHFEpidTRD &ref):
	60	AliHFEpidBase("")
	61	, fThresholdFile("")
	62	, fPIDMethod(kLQ)
	63	, fTRDthresholds(0x0)
	64	, fTRDelectronEfficiencies(0x0)
65	{
66	//
67	// Copy constructor
68	//
69	ref.Copy(*this);
70	}
71
72	//___________________________________________________________________
73	AliHFEpidTRD &AliHFEpidTRD::operator=(const AliHFEpidTRD &ref){
74	//
75	// Assignment operator
76	//
77	if(this != &ref){
78	ref.Copy(*this);
79	}
80	return *this;
81	}
82
83	//___________________________________________________________________
84	void AliHFEpidTRD::Copy(TObject &ref) const {
85	//
86	// Performs the copying of the object
87	//
88	AliHFEpidTRD &target = dynamic_cast<AliHFEpidTRD &>(ref);
89
90	target.fThresholdFile = fThresholdFile;
91	target.fPIDMethod = fPIDMethod;
92	if(fTRDthresholds){
93	target.fTRDthresholds = dynamic_cast<TMap *>(fTRDthresholds->Clone());
94	}
95	if(fTRDelectronEfficiencies){
96	target.fTRDelectronEfficiencies = new TAxis(*fTRDelectronEfficiencies);
97	}
98	AliHFEpidBase::Copy(ref);
99	}
100
101	//___________________________________________________________________
102	AliHFEpidTRD::~AliHFEpidTRD(){
103	//
104	// Destructor
105	//
106	if(fTRDthresholds){
107	fTRDthresholds->Delete();
108	delete fTRDthresholds;
109	}
110	if(fTRDelectronEfficiencies) delete fTRDelectronEfficiencies;
111	}
112
113	//______________________________________________________
114	Bool_t AliHFEpidTRD::InitializePID(){
115	//
116	// InitializePID: Load TRD thresholds and create the electron efficiency axis
117	// to navigate
118	//
119	LoadTRDthresholds();
120	// Fill the electron efficiencies axis for the TRD alone PID
121	Int_t nEffs = fTRDthresholds->GetEntries() + 1;
122	TIterator* thres =fTRDthresholds->MakeIterator();
123	TObjString *key = 0x0;
124	Double_t tmp = new Double_t[nEffs], titer = tmp;
125	while((key = dynamic_cast<TObjString >((thres)()))){
126	(*titer++) = static_cast<Double_t>(key->String().Atoi())/100.;
127	}
128	delete thres;
129	*titer = 1.;
130	// Sort the electron efficiencies and put them into the TAxis for later navigation
131	Int_t ind = new Int_t[nEffs], iiter = ind;
132	TMath::Sort(nEffs, tmp, ind, kFALSE);
133	Double_t eleffs = new Double_t[nEffs], eiter = eleffs;
134	while(eiter < &eleffs[nEffs]) (eiter++) = tmp[(iiter++)];
135	// print the content
136	Int_t cnt = 0;
137	if(GetDebugLevel() > 1){
138	printf("Printing electron efficiency bins:\n");
139	eiter = eleffs;
140	thres = fTRDthresholds->MakeIterator();
141	TObject *object = 0x0;
142	while(eiter < &eleffs[nEffs - 1]){
143	printf("eleffs[%d] = %f", cnt++, *(eiter++));
144	key = dynamic_cast<TObjString *>(thres->Next());
145	object = fTRDthresholds->GetValue(key->String().Data());
146	printf(", Content: %p\n", (void *)object);
147	}
148	}
149	delete[] tmp; delete[] ind;
150	fTRDelectronEfficiencies = new TAxis(nEffs - 1, eleffs);
151	if(GetDebugLevel() > 1){
152	printf("Printing axis content:\n");
153	for(Int_t ibin = fTRDelectronEfficiencies->GetFirst(); ibin <= fTRDelectronEfficiencies->GetLast(); ibin++)
154	printf("%d.) minimum: %f, maximum? %f\n", ibin, fTRDelectronEfficiencies->GetBinLowEdge(ibin), fTRDelectronEfficiencies->GetBinUpEdge(ibin));
155	}
156	delete[] eleffs;
157	return kTRUE;
158	}
159
160	//______________________________________________________
161	Int_t AliHFEpidTRD::IsSelected(AliVParticle *track){
162	//
163	// Does PID for TRD alone:
164	// PID thresholds based on 90% Electron Efficiency level approximated by a linear
165	// step function
166	//
167	AliESDtrack *esdTrack = 0x0;
168	if(!(esdTrack = dynamic_cast<AliESDtrack *>(track))) return kFALSE;
169	Double_t p = esdTrack->GetOuterParam() ? esdTrack->GetOuterParam()->P() : esdTrack->P();
170	if(p < 0.6) return 0;
171
172	// Get the Histograms
173	TH1 *threshist = GetTRDthresholds(0.91);
174	Int_t bin = 0;
175	if(p > threshist->GetXaxis()->GetXmax())
176	bin = threshist->GetXaxis()->GetLast();
177	else if(p < threshist->GetXaxis()->GetXmin())
178	bin = threshist->GetXaxis()->GetFirst();
179	else
180	bin = threshist->GetXaxis()->FindBin(p);
181
182	Double_t pidProbs[AliPID::kSPECIES];
183	esdTrack->GetTRDpid(pidProbs);
184	if(pidProbs[AliPID::kElectron] > threshist->GetBinContent(bin)) return 11;
185	return 0;
186	}
187
188	//___________________________________________________________________
189	void AliHFEpidTRD::LoadTRDthresholds(){
190	//
191	// Load TRD threshold histograms from File
192	//
193	TFile *mythresholds = TFile::Open(fThresholdFile);
194	TKey *object = 0x0;
195	TString electron_eff;
196	TObjArray *histos = 0x0;
197	Float_t eff;
198	TH1F *refhist = 0x0;
199	TIterator *keyIterator = mythresholds->GetListOfKeys()->MakeIterator();
200	TString histnames[2] = {"fHistThreshLQ", "fHistThreshNN"};
201	gROOT->cd();
202	while((object = dynamic_cast<TKey >((keyIterator)()))){
203	// Get the electron efficiency bin this histogram was taken with
204	electron_eff = object->GetName();
205	electron_eff = electron_eff.Remove(0,3);
206	eff = static_cast<Float_t>(electron_eff.Atoi())/100.;
207
208	// Get the threshold according to the selected
209	histos = dynamic_cast<TObjArray *>(object->ReadObj());
210	refhist = dynamic_cast<TH1F *>(histos->FindObject(histnames[fPIDMethod].Data()));
211	SetTRDthresholds(refhist, eff);
212	histos->Delete();
213	delete histos;
214	}
215	delete keyIterator;
216	mythresholds->Close();
217	delete mythresholds;
218	}
219
220	//___________________________________________________________________
221	void AliHFEpidTRD::SetTRDthresholds(TH1F *thresholds, Float_t electronEff){
222	//
223	// Set the threshold histogram for the TRD pid
224	//
225	fTRDthresholds->Add(new TObjString(Form("%d", TMath::Nint(electronEff * 100.))), new TH1F(*thresholds));
226	}
227
228	//___________________________________________________________________
229	TH1F *AliHFEpidTRD::GetTRDthresholds(Float_t electronEff){
230	Int_t bin = 0;
231	if(electronEff < fTRDelectronEfficiencies->GetXmin()) bin = fTRDelectronEfficiencies->GetFirst();
232	else if(electronEff > fTRDelectronEfficiencies->GetXmax()) bin = fTRDelectronEfficiencies->GetLast();
233	else bin = fTRDelectronEfficiencies->FindBin(electronEff);
234	TObjString keyname = Form("%d", TMath::Nint(fTRDelectronEfficiencies->GetBinLowEdge(bin)* 100.));
235	/* printf("Key: %s\n", keyname.String().Data());*/
236	TH1F thresholds = dynamic_cast<TH1F >((dynamic_cast<TPair *>(fTRDthresholds->FindObject(&keyname)))->Value());
237	/* printf("thresholds: %p\n", thresholds);*/
238	return thresholds;
239	}