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