TRD dEdx OADB framework: Nch and Ncls dependence added in the OADB object
[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 "AliVTrack.h"
40
41 #include "AliTRDPIDResponseObject.h"
42 #include "AliTRDPIDResponse.h"
43 //#include "AliTRDTKDInterpolator.h"
44 #include "AliTRDNDFast.h"
45 #include "AliTRDdEdxParams.h"
46
47 ClassImp(AliTRDPIDResponse)
48
49 //____________________________________________________________
50 AliTRDPIDResponse::AliTRDPIDResponse():
51   TObject()
52   ,fkPIDResponseObject(NULL)
53   ,fkTRDdEdxParams(NULL)
54   ,fGainNormalisationFactor(1.)
55 {
56   //
57   // Default constructor
58   //
59 }
60
61 //____________________________________________________________
62 AliTRDPIDResponse::AliTRDPIDResponse(const AliTRDPIDResponse &ref):
63   TObject(ref)
64   ,fkPIDResponseObject(NULL)
65   ,fkTRDdEdxParams(NULL)
66   ,fGainNormalisationFactor(ref.fGainNormalisationFactor)
67 {
68   //
69   // Copy constructor
70   //
71 }
72
73 //____________________________________________________________
74 AliTRDPIDResponse &AliTRDPIDResponse::operator=(const AliTRDPIDResponse &ref){
75   //
76   // Assignment operator
77   //
78   if(this == &ref) return *this;
79   
80   // Make copy
81   TObject::operator=(ref);
82   fGainNormalisationFactor = ref.fGainNormalisationFactor;
83   fkPIDResponseObject = ref.fkPIDResponseObject;
84   fkTRDdEdxParams = ref.fkTRDdEdxParams;
85
86   return *this;
87 }
88
89 //____________________________________________________________
90 AliTRDPIDResponse::~AliTRDPIDResponse(){
91   //
92   // Destructor
93   //
94   if(IsOwner()) {
95     delete fkPIDResponseObject;
96     delete fkTRDdEdxParams;
97   }
98 }
99
100 //____________________________________________________________
101 Bool_t AliTRDPIDResponse::Load(const Char_t * filename){
102   //
103   // Load References into the toolkit
104   //
105   AliDebug(1, "Loading reference histos from root file");
106   TDirectory *owd = gDirectory;// store old working directory
107   
108   if(!filename)
109     filename = Form("%s/STEER/LQ1dRef_v1.root",gSystem->ExpandPathName("$ALICE_ROOT"));
110   TFile *in = TFile::Open(filename);
111   if(!in){
112     AliError("Ref file not available.");
113     return kFALSE;
114   }
115   
116   gROOT->cd();
117   fkPIDResponseObject = dynamic_cast<const AliTRDPIDResponseObject *>(in->Get("TRDPIDResponse")->Clone());
118   in->Close(); delete in;
119   owd->cd();
120   SetBit(kIsOwner, kTRUE);
121   AliDebug(2, Form("Successfully loaded References for %d Momentum bins", fkPIDResponseObject->GetNumberOfMomentumBins()));
122   return kTRUE;
123 }
124
125
126
127 //____________________________________________________________
128 Double_t AliTRDPIDResponse::GetNumberOfSigmas(const AliVTrack *track, AliPID::EParticleType type) const
129 {
130   //
131   //calculate the TRD nSigma
132   //
133
134   const Double_t badval = -9999;
135   Double_t info[5]; for(int i=0; i<5; i++){info[i]=badval;}
136
137   const Double_t delta = GetSignalDelta(track, type, kFALSE, info);
138
139   const Double_t mean = info[0];
140   const Double_t res = info[1];
141   if(res<0){
142     return badval;
143   }
144
145   const Double_t sigma = mean*res;
146   const Double_t eps = 1e-12;
147   return delta/(sigma + eps);
148 }
149
150 //____________________________________________________________
151 Double_t AliTRDPIDResponse::GetSignalDelta( const AliVTrack* track, AliPID::EParticleType type, Bool_t ratio/*=kFALSE*/, Double_t *info/*=0x0*/) const
152 {
153   //
154   //calculate the TRD signal difference w.r.t. the expected
155   //output other information in info[]
156   //
157
158   const Double_t badval = -9999;
159   Double_t pTRD = -999;
160   for(Int_t ich=0; ich<6; ich++){
161     pTRD = track->GetTRDmomentum(ich);
162     if(pTRD>0)
163       break;
164   }
165   if(pTRD<0){
166     return badval;
167   }
168
169   if(!fkTRDdEdxParams){
170     AliError("fkTRDdEdxParams null");
171     return -99999;
172   }
173
174   const Double_t nch = track->GetTRDNchamberdEdx();
175   const Double_t ncls = track->GetTRDNclusterdEdx();
176
177   const TVectorF meanvec = fkTRDdEdxParams->GetMeanParameter(type, nch, ncls);
178   if(meanvec.GetNrows()==0){
179     return badval;
180   }
181
182   const TVectorF resvec  = fkTRDdEdxParams->GetSigmaParameter(type, nch, ncls);
183   if(resvec.GetNrows()==0){
184     return badval;
185   }
186
187   const Float_t *meanpar = meanvec.GetMatrixArray();
188   const Float_t *respar  = resvec.GetMatrixArray();
189
190   //============================================================================================<<<<<<<<<<<<<
191
192   const Double_t bg = pTRD/AliPID::ParticleMass(type);
193   const Double_t expsig = MeandEdxTR(&bg, meanpar);
194
195   if(info){
196     info[0]= expsig;
197     info[1]= ResolutiondEdxTR(&ncls, respar);
198   }
199
200   const Double_t eps = 1e-10;
201
202   if(ratio){
203     return track->GetTRDsignal()/(expsig +  eps);
204   }
205   else{
206     return track->GetTRDsignal() - expsig;
207   }
208 }
209
210
211 Double_t AliTRDPIDResponse::ResolutiondEdxTR(const Double_t * xx,  const Float_t * par)
212 {
213   //
214   //resolution 
215   //npar=3
216   //
217
218   const Double_t ncls = xx[0];
219
220   return par[0]+par[1]*TMath::Power(ncls, par[2]);
221 }
222
223 Double_t AliTRDPIDResponse::MeandEdxTR(const Double_t * xx,  const Float_t * pin)
224 {
225   //
226   //ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
227   //npar = 8 = 3+5
228   //
229
230   Float_t ptr[4]={0};
231   for(int ii=0; ii<3; ii++){
232     ptr[ii+1]=pin[ii];
233   }
234   return MeanTR(xx,ptr) + MeandEdx(xx,&(pin[3]));
235 }
236
237 Double_t AliTRDPIDResponse::MeanTR(const Double_t * xx,  const Float_t * par)
238 {
239   //
240   //LOGISTIC parametrization for TR, in unit of MIP
241   //npar = 4
242   //
243                                                                                                                                                                                                                                                         
244   const Double_t bg = xx[0];
245   const Double_t gamma = sqrt(1+bg*bg);
246
247   const Double_t p0 = TMath::Abs(par[1]);
248   const Double_t p1 = TMath::Abs(par[2]);
249   const Double_t p2 = TMath::Abs(par[3]);
250
251   const Double_t zz = TMath::Log(gamma);
252   const Double_t tryield = p0/( 1 + TMath::Exp(-p1*(zz-p2)) );
253
254   return par[0]+tryield;
255 }
256
257 Double_t AliTRDPIDResponse::MeandEdx(const Double_t * xx,  const Float_t * par)
258 {
259   //
260   //ALEPH parametrization for dEdx
261   //npar = 5
262   //
263
264   const Double_t bg = xx[0];
265   const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
266
267   const Double_t p0 = TMath::Abs(par[0]);
268   const Double_t p1 = TMath::Abs(par[1]);
269
270   const Double_t p2 = TMath::Abs(par[2]);
271
272   const Double_t p3 = TMath::Abs(par[3]);
273   const Double_t p4 = TMath::Abs(par[4]);
274
275   const Double_t aa = TMath::Power(beta, p3);
276
277   const Double_t bb = TMath::Log( p2 + TMath::Power(1./bg, p4) );
278
279   return (p1-aa-bb)*p0/aa;
280
281 }
282
283
284 //____________________________________________________________
285 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
286 {
287   //
288   // Calculate TRD likelihood values for the track based on dedx and 
289   // momentum values. The likelihoods are calculated by query the 
290   // reference data depending on the PID method selected
291   //
292   // Input parameter :
293   //   n - number of dedx slices/chamber
294   //   dedx - array of dedx slices organized layer wise
295   //   p - array of momentum measurements organized layer wise
296   // 
297   // Return parameters
298   //   prob - probabilities allocated by TRD for particle specis
299   //   kNorm - switch to normalize probabilities to 1. By default true. If false return not normalized prob.
300   // 
301   // Return value
302   //   number of tracklets used for PID, 0 if no PID
303     //
304     AliDebug(3,Form(" Response for PID method: %d",PIDmethod));
305
306
307     if(!fkPIDResponseObject){
308         AliDebug(3,"Missing reference data. PID calculation not possible.");
309         return 0;
310     }
311
312     for(Int_t is(AliPID::kSPECIES); is--;) prob[is]=.2;
313     Double_t prLayer[AliPID::kSPECIES];
314     Double_t dE[10], s(0.);
315     Int_t ntrackletsPID=0;
316     for(Int_t il(kNlayer); il--;){
317         memset(prLayer, 0, AliPID::kSPECIES*sizeof(Double_t));
318         if(!CookdEdx(n, &dedx[il*n], &dE[0],PIDmethod)) continue;
319         s=0.;
320         Bool_t filled=kTRUE;
321         for(Int_t is(AliPID::kSPECIES); is--;){
322             //if((PIDmethod==kLQ2D)&&(!(is==0||is==2)))continue;
323             if((dE[0] > 0.) && (p[il] > 0.)) prLayer[is] = GetProbabilitySingleLayer(is, p[il], &dE[0],PIDmethod);
324             AliDebug(3, Form("Probability for Species %d in Layer %d: %e", is, il, prLayer[is]));
325             if(prLayer[is]<1.e-30){
326                 AliDebug(2, Form("Null for species %d species prob layer[%d].",is,il));
327                 filled=kFALSE;
328                 break;
329             }
330             s+=prLayer[is];
331         }
332         if(!filled){
333             continue;
334         }
335         if(s<1.e-30){
336             AliDebug(2, Form("Null all species prob layer[%d].", il));
337             continue;
338         }
339         for(Int_t is(AliPID::kSPECIES); is--;){
340             if(kNorm) prLayer[is] /= s;
341             prob[is] *= prLayer[is];
342         }
343         ntrackletsPID++;
344     }
345     if(!kNorm) return ntrackletsPID;
346
347     s=0.;
348     for(Int_t is(AliPID::kSPECIES); is--;) s+=prob[is];
349     if(s<1.e-30){
350         AliDebug(2, "Null total prob.");
351         return 0;
352     }
353     for(Int_t is(AliPID::kSPECIES); is--;) prob[is]/=s;
354     return ntrackletsPID;
355 }
356
357 //____________________________________________________________
358 Double_t AliTRDPIDResponse::GetProbabilitySingleLayer(Int_t species, Double_t plocal, Double_t *dEdx,ETRDPIDMethod PIDmethod) const {
359   //
360   // Get the non-normalized probability for a certain particle species as coming
361   // from the reference histogram
362   // Interpolation between momentum bins
363   //
364   AliDebug(1, Form("Make Probability for Species %d with Momentum %f", species, plocal));
365
366   Double_t probLayer = 0.;
367
368   Float_t pLower, pUpper;
369         
370   AliTRDNDFast *refUpper = dynamic_cast<AliTRDNDFast *>(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,PIDmethod)),
371       *refLower = dynamic_cast<AliTRDNDFast *>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,PIDmethod));
372
373   // Do Interpolation exept for underflow and overflow
374   if(refLower && refUpper){
375       Double_t probLower = refLower->Eval(dEdx);
376       Double_t probUpper = refUpper->Eval(dEdx);
377
378       probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
379   } else if(refLower){
380      // underflow
381       probLayer = refLower->Eval(dEdx);
382   } else if(refUpper){
383       // overflow
384       probLayer = refUpper->Eval(dEdx);
385   } else {
386       AliError("No references available");
387   }
388   AliDebug(1, Form("Eval 1D dEdx %f Probability %e", dEdx[0],probLayer));
389
390   return probLayer;
391
392 /* old implementation
393
394 switch(PIDmethod){
395 case kNN: // NN
396       break;
397   case kLQ2D: // 2D LQ
398       {
399           if(species==0||species==2){ // references only for electrons and pions
400               Double_t error = 0.;
401               Double_t point[kNslicesLQ2D];
402               for(Int_t idim=0;idim<kNslicesLQ2D;idim++){point[idim]=dEdx[idim];}
403
404               AliTRDTKDInterpolator *refUpper = dynamic_cast<AliTRDTKDInterpolator *>(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,kLQ2D)),
405               *refLower = dynamic_cast<AliTRDTKDInterpolator *>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,kLQ2D));
406               // Do Interpolation exept for underflow and overflow
407               if(refLower && refUpper){
408                   Double_t probLower=0,probUpper=0;
409                   refLower->Eval(point,probLower,error);
410                   refUpper->Eval(point,probUpper,error);
411                   probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
412               } else if(refLower){
413                   // underflow
414                   refLower->Eval(point,probLayer,error);
415                   } else if(refUpper){
416                   // overflow
417                   refUpper->Eval(point,probLayer,error);
418               } else {
419                   AliError("No references available");
420               }
421               AliDebug(2,Form("Eval 2D Q0 %f Q1 %f P %e Err %e",point[0],point[1],probLayer,error));
422           }
423       }
424       break;
425   case kLQ1D: // 1D LQ
426       {
427           TH1 *refUpper = dynamic_cast<TH1 *>(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,kLQ1D)),
428               *refLower = dynamic_cast<TH1 *>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,kLQ1D));
429           // Do Interpolation exept for underflow and overflow
430           if(refLower && refUpper){
431               Double_t probLower = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx[0]));
432               Double_t probUpper = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx[0]));
433
434               probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
435           } else if(refLower){
436               // underflow
437               probLayer = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx[0]));
438           } else if(refUpper){
439               // overflow
440               probLayer = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx[0]));
441           } else {
442               AliError("No references available");
443           }
444           AliDebug(1, Form("Eval 1D dEdx %f Probability %e", dEdx[0],probLayer));
445       }
446       break;
447   default:
448       break;
449       }
450       return probLayer;
451       */
452
453 }
454
455 //____________________________________________________________
456 void AliTRDPIDResponse::SetOwner(){
457   //
458   // Make Deep Copy of the Reference Histograms
459   //
460     if(!fkPIDResponseObject || IsOwner()) return;
461     const AliTRDPIDResponseObject *tmp = fkPIDResponseObject;
462     fkPIDResponseObject = dynamic_cast<const AliTRDPIDResponseObject *>(tmp->Clone());
463     SetBit(kIsOwner, kTRUE);
464 }
465
466 //____________________________________________________________
467 Bool_t AliTRDPIDResponse::CookdEdx(Int_t nSlice, const Double_t * const in, Double_t *out,ETRDPIDMethod PIDmethod) const
468 {
469     //
470     // Recalculate dE/dx
471     //
472   switch(PIDmethod){
473   case kNN: // NN 
474       break;
475   case kLQ2D: // 2D LQ
476       out[0]=0;
477       out[1]=0;
478       for(Int_t islice = 0; islice < nSlice; islice++){
479           if(in[islice]<=0){out[0]=0;out[1]=0;return kFALSE;}  // Require that all slices are filled
480           if(islice<fkPIDResponseObject->GetNSlicesQ0())out[0]+= in[islice];
481           else out[1]+= in[islice];
482       }
483       // normalize signal to number of slices
484       out[0]*=1./Double_t(fkPIDResponseObject->GetNSlicesQ0());
485       out[1]*=1./Double_t(nSlice-fkPIDResponseObject->GetNSlicesQ0());
486       if(out[0] < 1e-6) return kFALSE;
487       AliDebug(3,Form("CookdEdx Q0 %f Q1 %f",out[0],out[1]));
488       break;
489   case kLQ1D: // 1D LQ
490       out[0]= 0.;
491       for(Int_t islice = 0; islice < nSlice; islice++)
492           if(in[islice] > 0) out[0] += in[islice] * fGainNormalisationFactor;   // Protect against negative values for slices having no dE/dx information
493       out[0]*=1./Double_t(nSlice);
494       if(out[0] < 1e-6) return kFALSE;
495       AliDebug(3,Form("CookdEdx dEdx %f",out[0]));
496       break;
497
498   default:
499       return kFALSE;
500   }
501   return kTRUE;
502 }
503
504 //____________________________________________________________
505 Bool_t AliTRDPIDResponse::IdentifiedAsElectron(Int_t nTracklets, const Double_t *like, Double_t p, Double_t level,Double_t centrality,ETRDPIDMethod PIDmethod) const {
506     //
507     // Check whether particle is identified as electron assuming a certain electron efficiency level
508     // Only electron and pion hypothesis is taken into account
509     //
510     // Inputs:
511     //         Number of tracklets
512     //         Likelihood values
513     //         Momentum
514     //         Electron efficiency level
515     //
516     // If the function fails when the params are not accessible, the function returns true
517     //
518   if(!fkPIDResponseObject){
519     AliDebug(3,"No PID Param object available");
520     return kTRUE;
521   } 
522   Double_t probEle = like[AliPID::kElectron]/(like[AliPID::kElectron] + like[AliPID::kPion]);
523   Double_t params[4];
524   if(!fkPIDResponseObject->GetThresholdParameters(nTracklets, level, params,centrality,PIDmethod)){
525     AliError("No Params found for the given configuration");
526     return kTRUE;
527   }
528   Double_t threshold = 1. - params[0] - params[1] * p - params[2] * TMath::Exp(-params[3] * p);
529   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
530   return kFALSE;
531 }
532
533 //____________________________________________________________
534 Bool_t AliTRDPIDResponse::SetPIDResponseObject(const AliTRDPIDResponseObject * obj){
535
536     fkPIDResponseObject = obj;
537     if((AliLog::GetDebugLevel("",IsA()->GetName()))>0)fkPIDResponseObject->Print("");
538     return kTRUE;
539 }
540
541
542 //____________________________________________________________    
543 Bool_t AliTRDPIDResponse::SetdEdxParams(const AliTRDdEdxParams * par)
544 {
545   fkTRDdEdxParams = par;
546   return kTRUE;
547 }