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>
30 #include <TObjArray.h>
34 #include <TDirectory.h>
38 #include "AliTRDPIDResponse.h"
40 ClassImp(AliTRDPIDResponse)
42 const Double_t AliTRDPIDResponse::fgkPBins[kNPBins] = {1, 2, 3, 4, 5, 6};
44 //____________________________________________________________
45 AliTRDPIDResponse::AliTRDPIDResponse(const Char_t * filename):
48 ,fGainNormalisationFactor(1.)
52 // Default constructor
54 for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++)
55 for(Int_t ipbin = 0; ipbin < kNPBins; ipbin++)
56 fMapRefHists[ispec][ipbin] = -1;
60 //____________________________________________________________
61 AliTRDPIDResponse::AliTRDPIDResponse(const AliTRDPIDResponse &ref):
63 ,fReferences(ref.fReferences)
64 ,fGainNormalisationFactor(ref.fGainNormalisationFactor)
65 ,fPIDmethod(ref.fPIDmethod)
69 // Flat copy of the reference histos
70 // For creating a deep copy call SetOwner
72 Int_t size = (AliPID::kSPECIES)*(kNPBins);
73 memcpy(fMapRefHists, ref.fMapRefHists, sizeof(Double_t) * size);
76 //____________________________________________________________
77 AliTRDPIDResponse &AliTRDPIDResponse::operator=(const AliTRDPIDResponse &ref){
79 // Assignment operator
80 // Performs a flat copy of the reference histos
82 if(this == &ref) return *this;
85 TObject::operator=(ref);
88 fReferences->Delete();
91 SetBit(kIsOwner, kFALSE);
92 } else if(fReferences) delete fReferences;
93 fReferences = ref.fReferences;
94 Int_t size = (AliPID::kSPECIES)*(kNPBins);
95 memcpy(fMapRefHists, ref.fMapRefHists, sizeof(Double_t) * size);
96 fPIDmethod = ref.fPIDmethod;
101 //____________________________________________________________
102 AliTRDPIDResponse::~AliTRDPIDResponse(){
108 if(fReferences) fReferences->Delete();
110 if(fReferences) delete fReferences;
113 //____________________________________________________________
114 Bool_t AliTRDPIDResponse::Load(const Char_t * filename){
116 // Load References into the toolkit
119 TDirectory *owd = gDirectory;// store old working directory
122 filename = Form("%s/STEER/LQ1dRef_v1.root",gSystem->ExpandPathName("$ALICE_ROOT"));
123 TFile *in = TFile::Open(filename);
125 AliError("Ref file not available.");
130 fReferences = new TObjArray();
131 fReferences->SetOwner();
132 TIter iter(in->GetListOfKeys());
133 TKey *histkey = NULL;
135 Int_t arrayPos = 0, pbin = 0;
136 AliPID::EParticleType species;
138 while((histkey = dynamic_cast<TKey *>(iter()))){
139 tmp = histkey->ReadObj();
140 histname = tmp->GetName();
141 if(histname.BeginsWith("fHQel")){ // Electron histogram
142 histname.ReplaceAll("fHQel_p","");
143 species = AliPID::kElectron;
144 } else { // Pion histogram
145 histname.ReplaceAll("fHQpi_p","");
146 species = AliPID::kPion;
148 pbin = histname.Atoi() - 1;
149 AliDebug(1, Form("Species %d, Histname %s, Pbin %d, Position in container %d", species, histname.Data(), pbin, arrayPos));
150 fMapRefHists[species][pbin] = arrayPos;
151 fReferences->AddAt(new TH1F(*dynamic_cast<TH1F *>(tmp)), arrayPos);
156 AliDebug(2, Form("Successfully loaded %d Reference Histograms", arrayPos));
157 SetBit(kIsOwner, kTRUE);
162 //____________________________________________________________
163 Bool_t AliTRDPIDResponse::GetResponse(Int_t n, Double_t *dedx, Float_t *p, Double_t prob[AliPID::kSPECIES], Bool_t kNorm) const
166 // Calculate TRD likelihood values for the track based on dedx and
167 // momentum values. The likelihoods are calculated by query the
168 // reference data depending on the PID method selected
171 // n - number of dedx slices/chamber
172 // dedx - array of dedx slices organized layer wise
173 // p - array of momentum measurements organized layer wise
176 // prob - probabilities allocated by TRD for particle specis
177 // kNorm - switch to normalize probabilities to 1. By default true. If false return not normalized prob.
180 // true if calculation success
184 AliWarning("Missing reference data. PID calculation not possible.");
188 for(Int_t is(AliPID::kSPECIES); is--;) prob[is]=.2;
189 Double_t prLayer[AliPID::kSPECIES];
190 Double_t DE[10], s(0.);
191 for(Int_t il(kNlayer); il--;){
192 memset(prLayer, 0, AliPID::kSPECIES*sizeof(Double_t));
193 if(!CookdEdx(&dedx[il*n], &DE[0])) continue;
196 for(Int_t is(AliPID::kSPECIES); is--;){
197 if((DE[0] > 0.) && (p[il] > 0.)) prLayer[is] = GetProbabilitySingleLayer(is, p[il], DE[0]);
198 AliDebug(3, Form("Probability for Species %d in Layer %d: %f", is, il, prLayer[is]));
202 AliDebug(2, Form("Null all species prob layer[%d].", il));
205 for(Int_t is(AliPID::kSPECIES); is--;){
206 if(kNorm) prLayer[is] /= s;
207 prob[is] *= prLayer[is];
210 if(!kNorm) return kTRUE;
213 for(Int_t is(AliPID::kSPECIES); is--;) s+=prob[is];
215 AliDebug(2, "Null total prob.");
218 for(Int_t is(AliPID::kSPECIES); is--;) prob[is]/=s;
223 //____________________________________________________________
224 Double_t AliTRDPIDResponse::GetProbabilitySingleLayer(Int_t species, Double_t plocal, Double_t dEdx) const {
226 // Get the non-normalized probability for a certain particle species as coming
227 // from the reference histogram
228 // Interpolation between momentum bins
230 AliDebug(1, Form("Make Probability for Species %d with Momentum %f", species, plocal));
231 Int_t pbin = GetLowerMomentumBin(plocal);
232 AliDebug(1, Form("Bin %d", pbin));
233 Double_t pLayer = 0.;
234 // Do Interpolation exept for underflow and overflow
235 if(pbin >= 0 && pbin < kNPBins-1){
236 TH1 *refHistUpper = NULL, *refHistLower = NULL;
237 if(fMapRefHists[species][pbin] >= 0)
238 refHistLower = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[species][pbin]));
239 if(fMapRefHists[species][pbin+1] >= 0)
240 refHistUpper = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[species][pbin+1]));
241 AliDebug(1, Form("Reference Histos (Upper/Lower): [%p|%p]", refHistUpper, refHistLower));
243 if (refHistLower && refHistUpper ) {
244 Double_t pLower = refHistLower->GetBinContent(refHistLower->GetXaxis()->FindBin(dEdx));
245 Double_t pUpper = refHistUpper->GetBinContent(refHistUpper->GetXaxis()->FindBin(dEdx));
247 pLayer = pLower + (pUpper - pLower)/(fgkPBins[pbin+1]-fgkPBins[pbin]) * (plocal - fgkPBins[pbin]);
254 if(fMapRefHists[species][0] >= 0){
255 refHist = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[species][0]));
259 if(fMapRefHists[species][kNPBins-1] > 0)
260 refHist = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[species][kNPBins-1]));
263 pLayer = refHist->GetBinContent(refHist->GetXaxis()->FindBin(dEdx));
265 AliDebug(1, Form("Probability %f", pLayer));
269 //____________________________________________________________
270 Int_t AliTRDPIDResponse::GetLowerMomentumBin(Double_t p) const {
272 // Get the momentum bin for a given momentum value
275 if(p > fgkPBins[kNPBins-1]) return kNPBins-1;
276 for(Int_t ibin = 0; ibin < kNPBins - 1; ibin++){
277 if(p >= fgkPBins[ibin] && p < fgkPBins[ibin+1]){
285 //____________________________________________________________
286 void AliTRDPIDResponse::SetOwner(){
288 // Make Deep Copy of the Reference Histograms
290 if(!fReferences || IsOwner()) return;
291 TObjArray *ctmp = new TObjArray();
292 for(Int_t ien = 0; ien < fReferences->GetEntriesFast(); ien++)
293 ctmp->AddAt(fReferences->UncheckedAt(ien)->Clone(), ien);
296 SetBit(kIsOwner, kTRUE);
299 //____________________________________________________________
300 Bool_t AliTRDPIDResponse::CookdEdx(Double_t *in, Double_t *out) const
309 for(Int_t islice = 0; islice < 8; islice++) out[0] += in[islice] * fGainNormalisationFactor;
310 if(out[0] < 1e-6) return kFALSE;