1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 // PID Response class for the TRD detector
17 // Based on 1D Likelihood approach
18 // Calculation of probabilities using Bayesian approach
19 // Attention: This method is only used to separate electrons from pions
22 // Markus Fasel <M.Fasel@gsi.de>
23 // Anton Andronic <A.Andronic@gsi.de>
29 #include <TObjArray.h>
35 #include "AliTRDPIDResponse.h"
37 ClassImp(AliTRDPIDResponse)
39 const Double_t AliTRDPIDResponse::fgkPBins[kNPBins] = {1, 2, 3, 4, 5, 6};
41 //____________________________________________________________
42 AliTRDPIDResponse::AliTRDPIDResponse(const Char_t * filename):
48 // Default constructor
50 memset(fMapRefHists, 0, sizeof(Int_t) * 10);
54 //____________________________________________________________
55 AliTRDPIDResponse::AliTRDPIDResponse(const AliTRDPIDResponse &ref):
57 ,fReferences(ref.fReferences)
58 ,fPIDmethod(ref.fPIDmethod)
62 // Flat copy of the reference histos
63 // For creating a deep copy call SetOwner
65 memcpy(fMapRefHists, ref.fMapRefHists, sizeof(Double_t) * 10);
68 //____________________________________________________________
69 AliTRDPIDResponse &AliTRDPIDResponse::operator=(const AliTRDPIDResponse &ref){
71 // Assignment operator
72 // Performs a flat copy of the reference histos
74 if(this == &ref) return *this;
77 TObject::operator=(ref);
80 fReferences->Delete();
83 SetBit(kIsOwner, kFALSE);
85 fReferences = ref.fReferences;
86 memcpy(fMapRefHists, ref.fMapRefHists, sizeof(Double_t) * 10);
87 fPIDmethod = ref.fPIDmethod;
92 //____________________________________________________________
93 AliTRDPIDResponse::~AliTRDPIDResponse(){
97 if(fReferences && IsOwner()){
98 fReferences->Delete();
103 //____________________________________________________________
104 Bool_t AliTRDPIDResponse::Load(const Char_t * filename){
106 // Load References into the toolkit
109 filename = Form("%s/STEER/LQ1dRef_v1.root",gSystem->ExpandPathName("$ALICE_ROOT"));
110 TFile *in = TFile::Open(filename);
112 AliError("Ref file not available.");
116 fReferences = new TObjArray();
117 fReferences->SetOwner();
118 TIter iter(in->GetListOfKeys());
120 Int_t arrayPos = 0, pbin = 0;
121 AliPID::EParticleType species;
123 while((tmp = iter())){
124 histname = tmp->GetName();
125 if(histname.BeginsWith("fHQel")){ // Electron histogram
126 histname.ReplaceAll("fHQel_p","");
127 species = AliPID::kElectron;
128 } else { // Pion histogram
129 histname.ReplaceAll("fHQpi_p","");
130 species = AliPID::kPion;
132 pbin = histname.Atoi();
133 fMapRefHists[species][pbin] = arrayPos;
134 fReferences->AddAt(tmp->Clone(), arrayPos);
139 AliDebug(2, Form("Successfully loaded %d Reference Histograms", arrayPos));
144 //____________________________________________________________
145 Bool_t AliTRDPIDResponse::GetResponse(Int_t n, Double_t *dedx, Float_t *p, Double_t prob[AliPID::kSPECIES], Bool_t kNorm)
148 // Calculate TRD likelihood values for the track based on dedx and
149 // momentum values. The likelihoods are calculated by query the
150 // reference data depending on the PID method selected
153 // n - number of dedx slices/chamber
154 // dedx - array of dedx slices organized layer wise
155 // p - array of momentum measurements organized layer wise
158 // prob - probabilities allocated by TRD for particle specis
159 // kNorm - switch to normalize probabilities to 1. By default true. If false return not normalized prob.
162 // true if calculation success
166 AliWarning("Missing reference data. PID calculation not possible.");
170 for(Int_t is(AliPID::kSPECIES); is--;) prob[is]=.2;
171 Double_t prLayer[AliPID::kSPECIES];
172 memset(prLayer, 0, AliPID::kSPECIES*sizeof(Double_t));
173 Double_t DE[10], s(0.);
174 for(Int_t il(kNlayer); il--;){
175 if(!CookdEdx(&dedx[il*n], &DE[0])) continue;
178 for(Int_t is(AliPID::kSPECIES); is--;){
179 prLayer[is] = GetProbabilitySingleLayer(is, p[il], DE[0]);
183 AliDebug(2, Form("Null all species prob layer[%d].", il));
186 for(Int_t is(AliPID::kSPECIES); is--;){
187 if(kNorm) prLayer[is] /= s;
188 prob[is] *= prLayer[is];
191 if(!kNorm) return kTRUE;
194 for(Int_t is(AliPID::kSPECIES); is--;) s+=prob[is];
196 AliDebug(2, "Null total prob.");
199 for(Int_t is(AliPID::kSPECIES); is--;) prob[is]/=s;
204 //____________________________________________________________
205 Double_t AliTRDPIDResponse::GetProbabilitySingleLayer(Int_t species, Double_t plocal, Double_t dEdx){
207 // Get the non-normalized probability for a certain particle species as coming
208 // from the reference histogram
209 // Interpolation between momentum bins
211 Int_t pbin = GetLowerMomentumBin(plocal);
212 Double_t pLayer = 0.;
213 // Do Interpolation exept for underflow and overflow
214 if(pbin >= 0 && pbin < kNPBins-1){
215 TH1 *refHistUpper = NULL, *refHistLower = NULL;
216 refHistLower = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[species][pbin]));
217 refHistUpper = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[species][pbin+1]));
219 Double_t pLower = refHistLower->GetBinContent(refHistLower->GetXaxis()->FindBin(dEdx));
220 Double_t pUpper = refHistUpper->GetBinContent(refHistUpper->GetXaxis()->FindBin(dEdx));
222 pLayer = pLower + (pUpper - pLower)/(fgkPBins[pbin+1]-fgkPBins[pbin]) * plocal;
228 refHist = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[0][0]));
231 refHist = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[species][kNPBins]));
233 pLayer = refHist->GetBinContent(refHist->GetXaxis()->FindBin(dEdx));
238 //____________________________________________________________
239 Int_t AliTRDPIDResponse::GetLowerMomentumBin(Double_t p){
241 // Get the momentum bin for a given momentum value
244 if(p > fgkPBins[kNPBins-1]) return kNPBins-1;
245 for(Int_t ibin = 0; ibin < kNPBins - 1; ibin++){
246 if(p >= fgkPBins[ibin] && p < fgkPBins[ibin+1]){
254 //____________________________________________________________
255 void AliTRDPIDResponse::SetOwner(){
257 // Make Deep Copy of the Reference Histograms
259 if(!fReferences || IsOwner()) return;
260 TObjArray *ctmp = new TObjArray();
261 for(Int_t ien = 0; ien < fReferences->GetEntriesFast(); ien++)
262 ctmp->AddAt(fReferences->UncheckedAt(ien)->Clone(), ien);
264 SetBit(kIsOwner, kTRUE);
267 //____________________________________________________________
268 Bool_t AliTRDPIDResponse::CookdEdx(Double_t *in, Double_t *out)
277 for(Int_t islice = 0; islice < 8; islice++) out[0] += in[islice];
278 if(out[0] < 1e-6) return kFALSE;