Updates in order to enable the '2D' PID for the TRD developed by Daniel Lohner.
[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
37 #include "AliLog.h"
38
39 #include "AliTRDPIDResponseObject.h"
40 //#include "AliTRDPIDParams.h"
41 //#include "AliTRDPIDReference.h"
42 #include "AliTRDPIDResponse.h"
43 #include "AliTRDTKDInterpolator.h"
44
45 ClassImp(AliTRDPIDResponse)
46
47 //____________________________________________________________
48 AliTRDPIDResponse::AliTRDPIDResponse():
49   TObject()
50   ,fkPIDResponseObject(NULL)
51   ,fGainNormalisationFactor(1.)
52   ,fPIDmethod(kLQ1D)
53 {
54   //
55   // Default constructor
56   //
57 }
58
59 //____________________________________________________________
60 AliTRDPIDResponse::AliTRDPIDResponse(const AliTRDPIDResponse &ref):
61   TObject(ref)
62   ,fkPIDResponseObject(NULL)
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   fkPIDResponseObject = ref.fkPIDResponseObject;
82   fPIDmethod = ref.fPIDmethod;
83   
84   return *this;
85 }
86
87 //____________________________________________________________
88 AliTRDPIDResponse::~AliTRDPIDResponse(){
89   //
90   // Destructor
91   //
92     if(IsOwner()) delete fkPIDResponseObject;
93 }
94
95 //____________________________________________________________
96 Bool_t AliTRDPIDResponse::Load(const Char_t * filename){
97   //
98   // Load References into the toolkit
99   //
100   AliDebug(1, "Loading reference histos from root file");
101   TDirectory *owd = gDirectory;// store old working directory
102   
103   if(!filename)
104     filename = Form("%s/STEER/LQ1dRef_v1.root",gSystem->ExpandPathName("$ALICE_ROOT"));
105   TFile *in = TFile::Open(filename);
106   if(!in){
107     AliError("Ref file not available.");
108     return kFALSE;
109   }
110   
111   gROOT->cd();
112   fkPIDResponseObject = dynamic_cast<const AliTRDPIDResponseObject *>(in->Get("TRDPIDResponse")->Clone());
113   in->Close(); delete in;
114   owd->cd();
115   SetBit(kIsOwner, kTRUE);
116   AliDebug(2, Form("Successfully loaded References for %d Momentum bins", fkPIDResponseObject->GetNumberOfMomentumBins()));
117   return kTRUE;
118 }
119
120 //____________________________________________________________
121 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
122 {
123   //
124   // Calculate TRD likelihood values for the track based on dedx and 
125   // momentum values. The likelihoods are calculated by query the 
126   // reference data depending on the PID method selected
127   //
128   // Input parameter :
129   //   n - number of dedx slices/chamber
130   //   dedx - array of dedx slices organized layer wise
131   //   p - array of momentum measurements organized layer wise
132   // 
133   // Return parameters
134   //   prob - probabilities allocated by TRD for particle specis
135   //   kNorm - switch to normalize probabilities to 1. By default true. If false return not normalized prob.
136   // 
137   // Return value
138   //   true if calculation success
139   // 
140
141     if(!fkPIDResponseObject){
142         AliWarning("Missing reference data. PID calculation not possible.");
143         return kFALSE;
144     }
145
146     for(Int_t is(AliPID::kSPECIES); is--;) prob[is]=.2;
147     Double_t prLayer[AliPID::kSPECIES];
148     Double_t dE[10], s(0.);
149     for(Int_t il(kNlayer); il--;){
150         memset(prLayer, 0, AliPID::kSPECIES*sizeof(Double_t));
151         if(!CookdEdx(n, &dedx[il*n], &dE[0])) continue;
152
153         s=0.;
154         Bool_t filled=kTRUE;
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             if(prLayer[is]<1.e-30){
159                 AliDebug(2, Form("Null for species %d species prob layer[%d].",is,il));
160                 filled=kFALSE;
161                 break;
162             }
163             s+=prLayer[is];
164         }
165         if(!filled){
166             continue;
167         }
168         if(s<1.e-30){
169             AliDebug(2, Form("Null all species prob layer[%d].", il));
170             continue;
171         }
172         for(Int_t is(AliPID::kSPECIES); is--;){
173             if(kNorm) prLayer[is] /= s;
174             prob[is] *= prLayer[is];
175         }
176     }
177     if(!kNorm) return kTRUE;
178
179     s=0.;
180     for(Int_t is(AliPID::kSPECIES); is--;) s+=prob[is];
181     if(s<1.e-30){
182         AliDebug(2, "Null total prob.");
183         return kFALSE;
184     }
185     for(Int_t is(AliPID::kSPECIES); is--;) prob[is]/=s;
186     return kTRUE;
187 }
188
189 //____________________________________________________________
190 Double_t AliTRDPIDResponse::GetProbabilitySingleLayer(Int_t species, Double_t plocal, Double_t *dEdx) const {
191   //
192   // Get the non-normalized probability for a certain particle species as coming
193   // from the reference histogram
194   // Interpolation between momentum bins
195   //
196   AliDebug(1, Form("Make Probability for Species %d with Momentum %f", species, plocal));
197
198   Double_t probLayer = 0.;
199   Float_t pLower, pUpper;
200         
201   switch(fPIDmethod){
202   case kNN: // NN
203       break;
204   case kLQ2D: // 2D LQ
205       {
206           Double_t error;
207           Double_t point[kNslicesLQ2D];
208           for(Int_t idim=0;idim<kNslicesLQ2D;idim++){point[idim]=dEdx[idim];}
209
210           AliTRDTKDInterpolator *refLower = dynamic_cast<AliTRDTKDInterpolator*>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,kLQ2D));
211
212           if(refLower){
213               refLower->Eval(point,probLayer,error);
214           }
215           else {
216               AliError("No references available");
217           }
218       }
219       break;
220   case kLQ1D: // 1D LQ
221       {
222           TH1 *refUpper = dynamic_cast<TH1 *>(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,kLQ1D)),
223               *refLower = dynamic_cast<TH1 *>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,kLQ1D));
224           // Do Interpolation exept for underflow and overflow
225           if(refLower && refUpper){
226               Double_t probLower = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx[0]));
227               Double_t probUpper = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx[0]));
228
229               probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
230           } else if(refLower){
231               // underflow
232               probLayer = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx[0]));
233           } else if(refUpper){
234               // overflow
235               probLayer = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx[0]));
236           } else {
237               AliError("No references available");
238           }
239           AliDebug(1, Form("Probability %f", probLayer));
240       }
241       break;
242   default:
243       break;
244   }
245   return probLayer;
246 }
247
248 //____________________________________________________________
249 void AliTRDPIDResponse::SetOwner(){
250   //
251   // Make Deep Copy of the Reference Histograms
252   //
253     if(!fkPIDResponseObject || IsOwner()) return;
254     const AliTRDPIDResponseObject *tmp = fkPIDResponseObject;
255     fkPIDResponseObject = dynamic_cast<const AliTRDPIDResponseObject *>(tmp->Clone());
256     SetBit(kIsOwner, kTRUE);
257 }
258
259 //____________________________________________________________
260 Bool_t AliTRDPIDResponse::CookdEdx(Int_t nSlice, const Double_t * const in, Double_t *out) const
261 {
262         //
263         // Recalculate dE/dx
264         //
265   switch(fPIDmethod){
266   case kNN: // NN 
267       break;
268   case kLQ2D: // 2D LQ
269       out[0]=0;
270       out[1]=0;
271       for(Int_t islice = 0; islice < nSlice; islice++){
272           if(islice<fkPIDResponseObject->GetNSlicesQ0())out[0]+= in[islice];
273           else out[1]+= in[islice];
274       }
275       break;
276   case kLQ1D: // 1D LQ
277       out[0]= 0.;
278       for(Int_t islice = 0; islice < nSlice; islice++)
279           if(in[islice] > 0) out[0] += in[islice] * fGainNormalisationFactor;   // Protect against negative values for slices having no dE/dx information
280       if(out[0] < 1e-6) return kFALSE;
281       break;
282
283   default:
284       return kFALSE;
285   }
286   return kTRUE;
287 }
288
289 //____________________________________________________________
290 Bool_t AliTRDPIDResponse::IdentifiedAsElectron(Int_t nTracklets, const Double_t *like, Double_t p, Double_t level) const {
291     //
292     // Check whether particle is identified as electron assuming a certain electron efficiency level
293     // Only electron and pion hypothesis is taken into account
294     //
295     // Inputs:
296     //         Number of tracklets
297     //         Likelihood values
298     //         Momentum
299     //         Electron efficiency level
300     //
301     // If the function fails when the params are not accessible, the function returns true
302     //
303   if(!fkPIDResponseObject){
304     AliError("No PID Param object available");
305     return kTRUE;
306   } 
307   Double_t probEle = like[AliPID::kElectron]/(like[AliPID::kElectron] + like[AliPID::kPion]);
308   Double_t params[4];
309   if(!fkPIDResponseObject->GetThresholdParameters(nTracklets, level, params,fPIDmethod)){
310     AliError("No Params found for the given configuration");
311     return kTRUE;
312   }
313   Double_t threshold = 1. - params[0] - params[1] * p - params[2] * TMath::Exp(-params[3] * p);
314   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
315   return kFALSE;
316 }
317