PID updates for TRD
[u/mrichter/AliRoot.git] / STEER / STEERBase / AliTRDPIDResponse.cxx
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 // 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
20 //
21 // Authors:
22 //  Markus Fasel <M.Fasel@gsi.de>
23 //  Anton Andronic <A.Andronic@gsi.de>
24 //
25 #include <TAxis.h>
26 #include <TClass.h>
27 #include <TDirectory.h>
28 #include <TFile.h>
29 #include <TH1.h>
30 #include <TKey.h>
31 #include <TMath.h>
32 #include <TObjArray.h>
33 #include <TROOT.h> 
34 #include <TString.h>
35 #include <TSystem.h>
36 #include <TVectorT.h>
37
38 #include "AliLog.h"
39
40 #include "AliTRDPIDReference.h"
41 #include "AliTRDPIDResponse.h"
42
43 ClassImp(AliTRDPIDResponse)
44
45 //____________________________________________________________
46 AliTRDPIDResponse::AliTRDPIDResponse():
47   TObject()
48   ,fkPIDReference(NULL)
49   ,fkPIDParams(NULL)
50   ,fGainNormalisationFactor(1.)
51   ,fPIDmethod(kLQ1D)
52 {
53   //
54   // Default constructor
55   //
56 }
57
58 //____________________________________________________________
59 AliTRDPIDResponse::AliTRDPIDResponse(const AliTRDPIDResponse &ref):
60   TObject(ref)
61   ,fkPIDReference(ref.fkPIDReference)
62   ,fkPIDParams(ref.fkPIDParams)
63   ,fGainNormalisationFactor(ref.fGainNormalisationFactor)
64   ,fPIDmethod(ref.fPIDmethod)
65 {
66   //
67   // Copy constructor
68   //
69 }
70
71 //____________________________________________________________
72 AliTRDPIDResponse &AliTRDPIDResponse::operator=(const AliTRDPIDResponse &ref){
73   //
74   // Assignment operator
75   //
76   if(this == &ref) return *this;
77   
78   // Make copy
79   TObject::operator=(ref);
80   fGainNormalisationFactor = ref.fGainNormalisationFactor;
81   fkPIDReference = ref.fkPIDReference;
82   fkPIDParams = ref.fkPIDParams;
83   fPIDmethod = ref.fPIDmethod;
84   
85   return *this;
86 }
87
88 //____________________________________________________________
89 AliTRDPIDResponse::~AliTRDPIDResponse(){
90   //
91   // Destructor
92   //
93   if(IsOwner()) delete fkPIDReference;
94 }
95
96 //____________________________________________________________
97 Bool_t AliTRDPIDResponse::Load(const Char_t * filename, const Char_t *refName){
98   //
99   // Load References into the toolkit
100   //
101   AliDebug(1, "Loading reference histos from root file");
102   TDirectory *owd = gDirectory;// store old working directory
103   
104   if(!filename)
105     filename = Form("%s/STEER/LQ1dRef_v1.root",gSystem->ExpandPathName("$ALICE_ROOT"));
106   TFile *in = TFile::Open(filename);
107   if(!in){
108     AliError("Ref file not available.");
109     return kFALSE;
110   }
111   
112   gROOT->cd();
113   fkPIDReference = dynamic_cast<const AliTRDPIDReference *>(in->Get(refName)->Clone());
114   in->Close(); delete in;
115   owd->cd();
116   SetBit(kIsOwner, kTRUE);
117   AliDebug(2, Form("Successfully loaded References for %d Momentum bins", fkPIDReference->GetNumberOfMomentumBins()));
118   return kTRUE;
119 }
120
121 //____________________________________________________________
122 Bool_t AliTRDPIDResponse::GetResponse(Int_t n, const Double_t * const dedx, const Float_t * const p, Double_t prob[AliPID::kSPECIES], Bool_t kNorm) const
123 {
124   //
125   // Calculate TRD likelihood values for the track based on dedx and 
126   // momentum values. The likelihoods are calculated by query the 
127   // reference data depending on the PID method selected
128   //
129   // Input parameter :
130   //   n - number of dedx slices/chamber
131   //   dedx - array of dedx slices organized layer wise
132   //   p - array of momentum measurements organized layer wise
133   // 
134   // Return parameters
135   //   prob - probabilities allocated by TRD for particle specis
136   //   kNorm - switch to normalize probabilities to 1. By default true. If false return not normalized prob.
137   // 
138   // Return value
139   //   true if calculation success
140   // 
141
142   if(!fkPIDReference){
143     AliWarning("Missing reference data. PID calculation not possible.");
144     return kFALSE;
145   }
146
147   for(Int_t is(AliPID::kSPECIES); is--;) prob[is]=.2;
148   Double_t prLayer[AliPID::kSPECIES];
149   Double_t dE[10], s(0.);
150   for(Int_t il(kNlayer); il--;){
151     memset(prLayer, 0, AliPID::kSPECIES*sizeof(Double_t));
152     if(!CookdEdx(n, &dedx[il*n], &dE[0])) continue;
153
154     s=0.;
155     for(Int_t is(AliPID::kSPECIES); is--;){
156       if((dE[0] > 0.) && (p[il] > 0.)) prLayer[is] = GetProbabilitySingleLayer(is, p[il], dE[0]);
157       AliDebug(3, Form("Probability for Species %d in Layer %d: %f", is, il, prLayer[is]));
158       s+=prLayer[is];
159     }
160     if(s<1.e-30){
161       AliDebug(2, Form("Null all species prob layer[%d].", il));
162       continue;
163     }
164     for(Int_t is(AliPID::kSPECIES); is--;){
165       if(kNorm) prLayer[is] /= s;
166       prob[is] *= prLayer[is];
167     }
168   }
169   if(!kNorm) return kTRUE;
170
171   s=0.;
172   for(Int_t is(AliPID::kSPECIES); is--;) s+=prob[is];
173   if(s<1.e-30){
174     AliDebug(2, "Null total prob.");
175     return kFALSE;
176   }
177   for(Int_t is(AliPID::kSPECIES); is--;) prob[is]/=s;
178   return kTRUE;
179 }
180
181 //____________________________________________________________
182 Double_t AliTRDPIDResponse::GetProbabilitySingleLayer(Int_t species, Double_t plocal, Double_t dEdx) const {
183   //
184   // Get the non-normalized probability for a certain particle species as coming
185   // from the reference histogram
186   // Interpolation between momentum bins
187   //
188   AliDebug(1, Form("Make Probability for Species %d with Momentum %f", species, plocal));
189   Float_t pLower, pUpper;
190   Double_t probLayer = 0.;
191   TH1 *refUpper = dynamic_cast<TH1 *>(fkPIDReference->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper)),
192       *refLower = dynamic_cast<TH1 *>(fkPIDReference->GetLowerReference((AliPID::EParticleType)species, plocal, pLower));
193   // Do Interpolation exept for underflow and overflow
194   if(refLower && refUpper){
195     Double_t probLower = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx));
196     Double_t probUpper = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx));
197   
198     probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
199   } else if(refLower){
200     // underflow
201     probLayer = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx));
202   } else if(refUpper){
203     // overflow
204     probLayer = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx));
205   } else {
206     AliError("No references available");
207   }
208   AliDebug(1, Form("Probability %f", probLayer));
209   return probLayer;
210 }
211
212 //____________________________________________________________
213 void AliTRDPIDResponse::SetOwner(){
214   //
215   // Make Deep Copy of the Reference Histograms
216   //
217   if(!fkPIDReference || IsOwner()) return;
218   const AliTRDPIDReference *tmp = fkPIDReference;
219   fkPIDReference = dynamic_cast<const AliTRDPIDReference *>(tmp->Clone());
220   SetBit(kIsOwner, kTRUE);
221 }
222
223 //____________________________________________________________
224 Bool_t AliTRDPIDResponse::CookdEdx(Int_t nSlice, const Double_t * const in, Double_t *out) const
225 {
226         //
227         // Recalculate dE/dx
228         //
229   switch(fPIDmethod){
230   case kNN: // NN 
231     break;
232   case kLQ2D: // 2D LQ 
233     break;
234   case kLQ1D: // 1D LQ 
235     out[0]= 0.;
236     for(Int_t islice = 0; islice < nSlice; islice++) 
237       if(in[islice] > 0) out[0] += in[islice] * fGainNormalisationFactor;   // Protect against negative values for slices having no dE/dx information
238     if(out[0] < 1e-6) return kFALSE;
239     break;
240   default:
241     return kFALSE;
242   }
243   return kTRUE;
244 }
245
246 //____________________________________________________________
247 Bool_t AliTRDPIDResponse::IdentifiedAsElectron(Int_t nTracklets, const Double_t *like, Double_t p, Double_t level) const {
248   //
249   // Check whether particle is identified as electron assuming a certain electron efficiency level
250   // Only electron and pion hypothesis is taken into account
251   //
252   // Inputs:
253   //         Number of tracklets
254   //         Likelihood values
255   //         Momentum
256   //         Electron efficiency level
257   //
258   // If the function fails when the params are not accessible, the function returns true
259   //
260   if(!fkPIDParams){
261     AliError("No PID Param object available");
262     return kTRUE;
263   } 
264   Double_t probEle = like[AliPID::kElectron]/(like[AliPID::kElectron] + like[AliPID::kPion]);
265   const TVectorD *params = GetParams(nTracklets, level);
266   if(!params){
267     AliError("No Params found for the given configuration");
268     return kTRUE;
269   }
270   Double_t threshold = 1. - (*params)[0] - (*params)[1] * p - (*params)[2] * TMath::Exp(-(*params)[3] * p);
271   if(probEle > TMath::Max(TMath::Min(threshold, 0.99), 0.2)) return kTRUE; // truncate the threshold upperwards to 0.999 and lowerwards to 0.2 and exclude unphysical values
272   return kFALSE;
273 }
274
275 //____________________________________________________________
276 const TVectorD* AliTRDPIDResponse::GetParams(Int_t ntracklets, Double_t level) const {
277   //
278   // returns the threshold for a given number of tracklets and a given efficiency level
279   //tby definition the lower of step is given.
280   //
281   if(ntracklets > 6 || ntracklets <=0) return NULL;
282   TObjArray * entry = dynamic_cast<TObjArray *>(fkPIDParams->At(ntracklets - 1));
283   if(!entry) return NULL;
284   
285   TObjArray*cut = NULL;
286   TVectorF *effLevel = NULL; const TVectorD *parameters = NULL;
287   Float_t currentLower = 0.;
288   TIter cutIter(entry);
289   while((cut = dynamic_cast<TObjArray *>(cutIter()))){
290     effLevel = static_cast<TVectorF *>(cut->At(0));
291     if((*effLevel)[0] > currentLower && (*effLevel)[0] <= level){
292       // New Lower entry found
293       parameters = static_cast<const TVectorD *>(cut->At(1));
294       currentLower = (*effLevel)[0];
295     }
296   }  
297   AliDebug(2, Form("Taking params for %d tracklets and %f electron Efficiency\n", ntracklets, currentLower));
298
299   return parameters;
300 }