- fast TRD PID implementation with muon constraints (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 "AliTRDPIDResponse.h"
41 //#include "AliTRDTKDInterpolator.h"
42 #include "AliTRDNDFast.h"
43
44 ClassImp(AliTRDPIDResponse)
45
46 //____________________________________________________________
47 AliTRDPIDResponse::AliTRDPIDResponse():
48   TObject()
49   ,fkPIDResponseObject(NULL)
50   ,fGainNormalisationFactor(1.)
51 {
52   //
53   // Default constructor
54   //
55 }
56
57 //____________________________________________________________
58 AliTRDPIDResponse::AliTRDPIDResponse(const AliTRDPIDResponse &ref):
59   TObject(ref)
60   ,fkPIDResponseObject(NULL)
61   ,fGainNormalisationFactor(ref.fGainNormalisationFactor)
62 {
63   //
64   // Copy constructor
65   //
66 }
67
68 //____________________________________________________________
69 AliTRDPIDResponse &AliTRDPIDResponse::operator=(const AliTRDPIDResponse &ref){
70   //
71   // Assignment operator
72   //
73   if(this == &ref) return *this;
74   
75   // Make copy
76   TObject::operator=(ref);
77   fGainNormalisationFactor = ref.fGainNormalisationFactor;
78   fkPIDResponseObject = ref.fkPIDResponseObject;
79   
80   return *this;
81 }
82
83 //____________________________________________________________
84 AliTRDPIDResponse::~AliTRDPIDResponse(){
85   //
86   // Destructor
87   //
88     if(IsOwner()) delete fkPIDResponseObject;
89 }
90
91 //____________________________________________________________
92 Bool_t AliTRDPIDResponse::Load(const Char_t * filename){
93   //
94   // Load References into the toolkit
95   //
96   AliDebug(1, "Loading reference histos from root file");
97   TDirectory *owd = gDirectory;// store old working directory
98   
99   if(!filename)
100     filename = Form("%s/STEER/LQ1dRef_v1.root",gSystem->ExpandPathName("$ALICE_ROOT"));
101   TFile *in = TFile::Open(filename);
102   if(!in){
103     AliError("Ref file not available.");
104     return kFALSE;
105   }
106   
107   gROOT->cd();
108   fkPIDResponseObject = dynamic_cast<const AliTRDPIDResponseObject *>(in->Get("TRDPIDResponse")->Clone());
109   in->Close(); delete in;
110   owd->cd();
111   SetBit(kIsOwner, kTRUE);
112   AliDebug(2, Form("Successfully loaded References for %d Momentum bins", fkPIDResponseObject->GetNumberOfMomentumBins()));
113   return kTRUE;
114 }
115
116 //____________________________________________________________
117 Int_t AliTRDPIDResponse::GetResponse(Int_t n, const Double_t * const dedx, const Float_t * const p, Double_t prob[AliPID::kSPECIES],ETRDPIDMethod PIDmethod,Bool_t kNorm) const
118 {
119   //
120   // Calculate TRD likelihood values for the track based on dedx and 
121   // momentum values. The likelihoods are calculated by query the 
122   // reference data depending on the PID method selected
123   //
124   // Input parameter :
125   //   n - number of dedx slices/chamber
126   //   dedx - array of dedx slices organized layer wise
127   //   p - array of momentum measurements organized layer wise
128   // 
129   // Return parameters
130   //   prob - probabilities allocated by TRD for particle specis
131   //   kNorm - switch to normalize probabilities to 1. By default true. If false return not normalized prob.
132   // 
133   // Return value
134   //   number of tracklets used for PID, 0 if no PID
135     //
136     AliDebug(3,Form(" Response for PID method: %d",PIDmethod));
137
138
139     if(!fkPIDResponseObject){
140         AliDebug(3,"Missing reference data. PID calculation not possible.");
141         return 0;
142     }
143
144     for(Int_t is(AliPID::kSPECIES); is--;) prob[is]=.2;
145     Double_t prLayer[AliPID::kSPECIES];
146     Double_t dE[10], s(0.);
147     Int_t ntrackletsPID=0;
148     for(Int_t il(kNlayer); il--;){
149         memset(prLayer, 0, AliPID::kSPECIES*sizeof(Double_t));
150         if(!CookdEdx(n, &dedx[il*n], &dE[0],PIDmethod)) continue;
151         s=0.;
152         Bool_t filled=kTRUE;
153         for(Int_t is(AliPID::kSPECIES); is--;){
154             //if((PIDmethod==kLQ2D)&&(!(is==0||is==2)))continue;
155             if((dE[0] > 0.) && (p[il] > 0.)) prLayer[is] = GetProbabilitySingleLayer(is, p[il], &dE[0],PIDmethod);
156             AliDebug(3, Form("Probability for Species %d in Layer %d: %e", is, il, prLayer[is]));
157             if(prLayer[is]<1.e-30){
158                 AliDebug(2, Form("Null for species %d species prob layer[%d].",is,il));
159                 filled=kFALSE;
160                 break;
161             }
162             s+=prLayer[is];
163         }
164         if(!filled){
165             continue;
166         }
167         if(s<1.e-30){
168             AliDebug(2, Form("Null all species prob layer[%d].", il));
169             continue;
170         }
171         for(Int_t is(AliPID::kSPECIES); is--;){
172             if(kNorm) prLayer[is] /= s;
173             prob[is] *= prLayer[is];
174         }
175         ntrackletsPID++;
176     }
177     if(!kNorm) return ntrackletsPID;
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 0;
184     }
185     for(Int_t is(AliPID::kSPECIES); is--;) prob[is]/=s;
186     return ntrackletsPID;
187 }
188
189 //____________________________________________________________
190 Double_t AliTRDPIDResponse::GetProbabilitySingleLayer(Int_t species, Double_t plocal, Double_t *dEdx,ETRDPIDMethod PIDmethod) 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
200   Float_t pLower, pUpper;
201         
202   AliTRDNDFast *refUpper = dynamic_cast<AliTRDNDFast *>(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,PIDmethod)),
203       *refLower = dynamic_cast<AliTRDNDFast *>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,PIDmethod));
204
205   // Do Interpolation exept for underflow and overflow
206   if(refLower && refUpper){
207       Double_t probLower = refLower->Eval(dEdx);
208       Double_t probUpper = refUpper->Eval(dEdx);
209
210       probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
211   } else if(refLower){
212      // underflow
213       probLayer = refLower->Eval(dEdx);
214   } else if(refUpper){
215       // overflow
216       probLayer = refUpper->Eval(dEdx);
217   } else {
218       AliError("No references available");
219   }
220   AliDebug(1, Form("Eval 1D dEdx %f Probability %e", dEdx[0],probLayer));
221
222   return probLayer;
223
224 /* old implementation
225
226 switch(PIDmethod){
227 case kNN: // NN
228       break;
229   case kLQ2D: // 2D LQ
230       {
231           if(species==0||species==2){ // references only for electrons and pions
232               Double_t error = 0.;
233               Double_t point[kNslicesLQ2D];
234               for(Int_t idim=0;idim<kNslicesLQ2D;idim++){point[idim]=dEdx[idim];}
235
236               AliTRDTKDInterpolator *refUpper = dynamic_cast<AliTRDTKDInterpolator *>(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,kLQ2D)),
237               *refLower = dynamic_cast<AliTRDTKDInterpolator *>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,kLQ2D));
238               // Do Interpolation exept for underflow and overflow
239               if(refLower && refUpper){
240                   Double_t probLower=0,probUpper=0;
241                   refLower->Eval(point,probLower,error);
242                   refUpper->Eval(point,probUpper,error);
243                   probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
244               } else if(refLower){
245                   // underflow
246                   refLower->Eval(point,probLayer,error);
247                   } else if(refUpper){
248                   // overflow
249                   refUpper->Eval(point,probLayer,error);
250               } else {
251                   AliError("No references available");
252               }
253               AliDebug(2,Form("Eval 2D Q0 %f Q1 %f P %e Err %e",point[0],point[1],probLayer,error));
254           }
255       }
256       break;
257   case kLQ1D: // 1D LQ
258       {
259           TH1 *refUpper = dynamic_cast<TH1 *>(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,kLQ1D)),
260               *refLower = dynamic_cast<TH1 *>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,kLQ1D));
261           // Do Interpolation exept for underflow and overflow
262           if(refLower && refUpper){
263               Double_t probLower = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx[0]));
264               Double_t probUpper = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx[0]));
265
266               probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
267           } else if(refLower){
268               // underflow
269               probLayer = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx[0]));
270           } else if(refUpper){
271               // overflow
272               probLayer = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx[0]));
273           } else {
274               AliError("No references available");
275           }
276           AliDebug(1, Form("Eval 1D dEdx %f Probability %e", dEdx[0],probLayer));
277       }
278       break;
279   default:
280       break;
281       }
282       return probLayer;
283       */
284
285 }
286
287 //____________________________________________________________
288 void AliTRDPIDResponse::SetOwner(){
289   //
290   // Make Deep Copy of the Reference Histograms
291   //
292     if(!fkPIDResponseObject || IsOwner()) return;
293     const AliTRDPIDResponseObject *tmp = fkPIDResponseObject;
294     fkPIDResponseObject = dynamic_cast<const AliTRDPIDResponseObject *>(tmp->Clone());
295     SetBit(kIsOwner, kTRUE);
296 }
297
298 //____________________________________________________________
299 Bool_t AliTRDPIDResponse::CookdEdx(Int_t nSlice, const Double_t * const in, Double_t *out,ETRDPIDMethod PIDmethod) const
300 {
301     //
302     // Recalculate dE/dx
303     //
304   switch(PIDmethod){
305   case kNN: // NN 
306       break;
307   case kLQ2D: // 2D LQ
308       out[0]=0;
309       out[1]=0;
310       for(Int_t islice = 0; islice < nSlice; islice++){
311           if(in[islice]<=0){out[0]=0;out[1]=0;return kFALSE;}  // Require that all slices are filled
312           if(islice<fkPIDResponseObject->GetNSlicesQ0())out[0]+= in[islice];
313           else out[1]+= in[islice];
314       }
315       // normalize signal to number of slices
316       out[0]*=1./Double_t(fkPIDResponseObject->GetNSlicesQ0());
317       out[1]*=1./Double_t(nSlice-fkPIDResponseObject->GetNSlicesQ0());
318       if(out[0] < 1e-6) return kFALSE;
319       AliDebug(3,Form("CookdEdx Q0 %f Q1 %f",out[0],out[1]));
320       break;
321   case kLQ1D: // 1D LQ
322       out[0]= 0.;
323       for(Int_t islice = 0; islice < nSlice; islice++)
324           if(in[islice] > 0) out[0] += in[islice] * fGainNormalisationFactor;   // Protect against negative values for slices having no dE/dx information
325       out[0]*=1./Double_t(nSlice);
326       if(out[0] < 1e-6) return kFALSE;
327       AliDebug(3,Form("CookdEdx dEdx %f",out[0]));
328       break;
329
330   default:
331       return kFALSE;
332   }
333   return kTRUE;
334 }
335
336 //____________________________________________________________
337 Bool_t AliTRDPIDResponse::IdentifiedAsElectron(Int_t nTracklets, const Double_t *like, Double_t p, Double_t level,Double_t centrality,ETRDPIDMethod PIDmethod) const {
338     //
339     // Check whether particle is identified as electron assuming a certain electron efficiency level
340     // Only electron and pion hypothesis is taken into account
341     //
342     // Inputs:
343     //         Number of tracklets
344     //         Likelihood values
345     //         Momentum
346     //         Electron efficiency level
347     //
348     // If the function fails when the params are not accessible, the function returns true
349     //
350   if(!fkPIDResponseObject){
351     AliDebug(3,"No PID Param object available");
352     return kTRUE;
353   } 
354   Double_t probEle = like[AliPID::kElectron]/(like[AliPID::kElectron] + like[AliPID::kPion]);
355   Double_t params[4];
356   if(!fkPIDResponseObject->GetThresholdParameters(nTracklets, level, params,centrality,PIDmethod)){
357     AliError("No Params found for the given configuration");
358     return kTRUE;
359   }
360   Double_t threshold = 1. - params[0] - params[1] * p - params[2] * TMath::Exp(-params[3] * p);
361   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
362   return kFALSE;
363 }
364
365 //____________________________________________________________
366 Bool_t AliTRDPIDResponse::SetPIDResponseObject(const AliTRDPIDResponseObject * obj){
367
368     fkPIDResponseObject = obj;
369     if((AliLog::GetDebugLevel("",IsA()->GetName()))>0)fkPIDResponseObject->Print("");
370     return kTRUE;
371 }