cd2b72b2d1818823b113a32b7a93c275272d9a39
[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
160   if(!track){
161     return badval;
162   }
163
164   Double_t pTRD = -999;
165   for(Int_t ich=0; ich<6; ich++){
166     pTRD = track->GetTRDmomentum(ich);
167     if(pTRD>0)
168       break;
169   }
170   if(pTRD<0){
171     return badval;
172   }
173
174   if(!fkTRDdEdxParams){
175     AliError("fkTRDdEdxParams null");
176     return -99999;
177   }
178
179   const Double_t nch = track->GetTRDNchamberdEdx();
180   const Double_t ncls = track->GetTRDNclusterdEdx();
181
182   const TVectorF meanvec = fkTRDdEdxParams->GetMeanParameter(type, nch, ncls);
183   if(meanvec.GetNrows()==0){
184     return badval;
185   }
186
187   const TVectorF resvec  = fkTRDdEdxParams->GetSigmaParameter(type, nch, ncls);
188   if(resvec.GetNrows()==0){
189     return badval;
190   }
191
192   const Float_t *meanpar = meanvec.GetMatrixArray();
193   const Float_t *respar  = resvec.GetMatrixArray();
194
195   //============================================================================================<<<<<<<<<<<<<
196
197   const Double_t bg = pTRD/AliPID::ParticleMass(type);
198   const Double_t expsig = MeandEdxTR(&bg, meanpar);
199
200   if(info){
201     info[0]= expsig;
202     info[1]= ResolutiondEdxTR(&ncls, respar);
203   }
204
205   const Double_t eps = 1e-10;
206
207   if(ratio){
208     return track->GetTRDsignal()/(expsig +  eps);
209   }
210   else{
211     return track->GetTRDsignal() - expsig;
212   }
213 }
214
215
216 Double_t AliTRDPIDResponse::ResolutiondEdxTR(const Double_t * xx,  const Float_t * par)
217 {
218   //
219   //resolution 
220   //npar=3
221   //
222
223   const Double_t ncls = xx[0];
224
225   return par[0]+par[1]*TMath::Power(ncls, par[2]);
226 }
227
228 Double_t AliTRDPIDResponse::MeandEdxTR(const Double_t * xx,  const Float_t * pin)
229 {
230   //
231   //ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
232   //npar = 8 = 3+5
233   //
234
235   Float_t ptr[4]={0};
236   for(int ii=0; ii<3; ii++){
237     ptr[ii+1]=pin[ii];
238   }
239   return MeanTR(xx,ptr) + MeandEdx(xx,&(pin[3]));
240 }
241
242 Double_t AliTRDPIDResponse::MeanTR(const Double_t * xx,  const Float_t * par)
243 {
244   //
245   //LOGISTIC parametrization for TR, in unit of MIP
246   //npar = 4
247   //
248                                                                                                                                                                                                                                                         
249   const Double_t bg = xx[0];
250   const Double_t gamma = sqrt(1+bg*bg);
251
252   const Double_t p0 = TMath::Abs(par[1]);
253   const Double_t p1 = TMath::Abs(par[2]);
254   const Double_t p2 = TMath::Abs(par[3]);
255
256   const Double_t zz = TMath::Log(gamma);
257   const Double_t tryield = p0/( 1 + TMath::Exp(-p1*(zz-p2)) );
258
259   return par[0]+tryield;
260 }
261
262 Double_t AliTRDPIDResponse::MeandEdx(const Double_t * xx,  const Float_t * par)
263 {
264   //
265   //ALEPH parametrization for dEdx
266   //npar = 5
267   //
268
269   const Double_t bg = xx[0];
270   const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
271
272   const Double_t p0 = TMath::Abs(par[0]);
273   const Double_t p1 = TMath::Abs(par[1]);
274
275   const Double_t p2 = TMath::Abs(par[2]);
276
277   const Double_t p3 = TMath::Abs(par[3]);
278   const Double_t p4 = TMath::Abs(par[4]);
279
280   const Double_t aa = TMath::Power(beta, p3);
281
282   const Double_t bb = TMath::Log( p2 + TMath::Power(1./bg, p4) );
283
284   return (p1-aa-bb)*p0/aa;
285
286 }
287
288
289 //____________________________________________________________
290 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
291 {
292   //
293   // Calculate TRD likelihood values for the track based on dedx and 
294   // momentum values. The likelihoods are calculated by query the 
295   // reference data depending on the PID method selected
296   //
297   // Input parameter :
298   //   n - number of dedx slices/chamber
299   //   dedx - array of dedx slices organized layer wise
300   //   p - array of momentum measurements organized layer wise
301   // 
302   // Return parameters
303   //   prob - probabilities allocated by TRD for particle specis
304   //   kNorm - switch to normalize probabilities to 1. By default true. If false return not normalized prob.
305   // 
306   // Return value
307   //   number of tracklets used for PID, 0 if no PID
308     //
309     AliDebug(3,Form(" Response for PID method: %d",PIDmethod));
310
311
312     if(!fkPIDResponseObject){
313         AliDebug(3,"Missing reference data. PID calculation not possible.");
314         return 0;
315     }
316
317     for(Int_t is(AliPID::kSPECIES); is--;) prob[is]=.2;
318     Double_t prLayer[AliPID::kSPECIES];
319     Double_t dE[10], s(0.);
320     Int_t ntrackletsPID=0;
321     for(Int_t il(kNlayer); il--;){
322         memset(prLayer, 0, AliPID::kSPECIES*sizeof(Double_t));
323         if(!CookdEdx(n, &dedx[il*n], &dE[0],PIDmethod)) continue;
324         s=0.;
325         Bool_t filled=kTRUE;
326         for(Int_t is(AliPID::kSPECIES); is--;){
327             //if((PIDmethod==kLQ2D)&&(!(is==0||is==2)))continue;
328             if((dE[0] > 0.) && (p[il] > 0.)) prLayer[is] = GetProbabilitySingleLayer(is, p[il], &dE[0],PIDmethod);
329             AliDebug(3, Form("Probability for Species %d in Layer %d: %e", is, il, prLayer[is]));
330             if(prLayer[is]<1.e-30){
331                 AliDebug(2, Form("Null for species %d species prob layer[%d].",is,il));
332                 filled=kFALSE;
333                 break;
334             }
335             s+=prLayer[is];
336         }
337         if(!filled){
338             continue;
339         }
340         if(s<1.e-30){
341             AliDebug(2, Form("Null all species prob layer[%d].", il));
342             continue;
343         }
344         for(Int_t is(AliPID::kSPECIES); is--;){
345             if(kNorm) prLayer[is] /= s;
346             prob[is] *= prLayer[is];
347         }
348         ntrackletsPID++;
349     }
350     if(!kNorm) return ntrackletsPID;
351
352     s=0.;
353     for(Int_t is(AliPID::kSPECIES); is--;) s+=prob[is];
354     if(s<1.e-30){
355         AliDebug(2, "Null total prob.");
356         return 0;
357     }
358     for(Int_t is(AliPID::kSPECIES); is--;) prob[is]/=s;
359     return ntrackletsPID;
360 }
361
362 //____________________________________________________________
363 Double_t AliTRDPIDResponse::GetProbabilitySingleLayer(Int_t species, Double_t plocal, Double_t *dEdx,ETRDPIDMethod PIDmethod) const {
364   //
365   // Get the non-normalized probability for a certain particle species as coming
366   // from the reference histogram
367   // Interpolation between momentum bins
368   //
369   AliDebug(1, Form("Make Probability for Species %d with Momentum %f", species, plocal));
370
371   Double_t probLayer = 0.;
372
373   Float_t pLower, pUpper;
374         
375   AliTRDNDFast *refUpper = dynamic_cast<AliTRDNDFast *>(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,PIDmethod)),
376       *refLower = dynamic_cast<AliTRDNDFast *>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,PIDmethod));
377
378   // Do Interpolation exept for underflow and overflow
379   if(refLower && refUpper){
380       Double_t probLower = refLower->Eval(dEdx);
381       Double_t probUpper = refUpper->Eval(dEdx);
382
383       probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
384   } else if(refLower){
385      // underflow
386       probLayer = refLower->Eval(dEdx);
387   } else if(refUpper){
388       // overflow
389       probLayer = refUpper->Eval(dEdx);
390   } else {
391       AliError("No references available");
392   }
393   AliDebug(1, Form("Eval 1D dEdx %f Probability %e", dEdx[0],probLayer));
394
395   return probLayer;
396
397 /* old implementation
398
399 switch(PIDmethod){
400 case kNN: // NN
401       break;
402   case kLQ2D: // 2D LQ
403       {
404           if(species==0||species==2){ // references only for electrons and pions
405               Double_t error = 0.;
406               Double_t point[kNslicesLQ2D];
407               for(Int_t idim=0;idim<kNslicesLQ2D;idim++){point[idim]=dEdx[idim];}
408
409               AliTRDTKDInterpolator *refUpper = dynamic_cast<AliTRDTKDInterpolator *>(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,kLQ2D)),
410               *refLower = dynamic_cast<AliTRDTKDInterpolator *>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,kLQ2D));
411               // Do Interpolation exept for underflow and overflow
412               if(refLower && refUpper){
413                   Double_t probLower=0,probUpper=0;
414                   refLower->Eval(point,probLower,error);
415                   refUpper->Eval(point,probUpper,error);
416                   probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
417               } else if(refLower){
418                   // underflow
419                   refLower->Eval(point,probLayer,error);
420                   } else if(refUpper){
421                   // overflow
422                   refUpper->Eval(point,probLayer,error);
423               } else {
424                   AliError("No references available");
425               }
426               AliDebug(2,Form("Eval 2D Q0 %f Q1 %f P %e Err %e",point[0],point[1],probLayer,error));
427           }
428       }
429       break;
430   case kLQ1D: // 1D LQ
431       {
432           TH1 *refUpper = dynamic_cast<TH1 *>(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,kLQ1D)),
433               *refLower = dynamic_cast<TH1 *>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,kLQ1D));
434           // Do Interpolation exept for underflow and overflow
435           if(refLower && refUpper){
436               Double_t probLower = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx[0]));
437               Double_t probUpper = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx[0]));
438
439               probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
440           } else if(refLower){
441               // underflow
442               probLayer = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx[0]));
443           } else if(refUpper){
444               // overflow
445               probLayer = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx[0]));
446           } else {
447               AliError("No references available");
448           }
449           AliDebug(1, Form("Eval 1D dEdx %f Probability %e", dEdx[0],probLayer));
450       }
451       break;
452   default:
453       break;
454       }
455       return probLayer;
456       */
457
458 }
459
460 //____________________________________________________________
461 void AliTRDPIDResponse::SetOwner(){
462   //
463   // Make Deep Copy of the Reference Histograms
464   //
465     if(!fkPIDResponseObject || IsOwner()) return;
466     const AliTRDPIDResponseObject *tmp = fkPIDResponseObject;
467     fkPIDResponseObject = dynamic_cast<const AliTRDPIDResponseObject *>(tmp->Clone());
468     SetBit(kIsOwner, kTRUE);
469 }
470
471 //____________________________________________________________
472 Bool_t AliTRDPIDResponse::CookdEdx(Int_t nSlice, const Double_t * const in, Double_t *out,ETRDPIDMethod PIDmethod) const
473 {
474     //
475     // Recalculate dE/dx
476     //
477   switch(PIDmethod){
478   case kNN: // NN 
479       break;
480   case kLQ2D: // 2D LQ
481       out[0]=0;
482       out[1]=0;
483       for(Int_t islice = 0; islice < nSlice; islice++){
484           if(in[islice]<=0){out[0]=0;out[1]=0;return kFALSE;}  // Require that all slices are filled
485           if(islice<fkPIDResponseObject->GetNSlicesQ0())out[0]+= in[islice];
486           else out[1]+= in[islice];
487       }
488       // normalize signal to number of slices
489       out[0]*=1./Double_t(fkPIDResponseObject->GetNSlicesQ0());
490       out[1]*=1./Double_t(nSlice-fkPIDResponseObject->GetNSlicesQ0());
491       if(out[0] < 1e-6) return kFALSE;
492       AliDebug(3,Form("CookdEdx Q0 %f Q1 %f",out[0],out[1]));
493       break;
494   case kLQ1D: // 1D LQ
495       out[0]= 0.;
496       for(Int_t islice = 0; islice < nSlice; islice++)
497           if(in[islice] > 0) out[0] += in[islice] * fGainNormalisationFactor;   // Protect against negative values for slices having no dE/dx information
498       out[0]*=1./Double_t(nSlice);
499       if(out[0] < 1e-6) return kFALSE;
500       AliDebug(3,Form("CookdEdx dEdx %f",out[0]));
501       break;
502
503   default:
504       return kFALSE;
505   }
506   return kTRUE;
507 }
508
509 //____________________________________________________________
510 Bool_t AliTRDPIDResponse::IdentifiedAsElectron(Int_t nTracklets, const Double_t *like, Double_t p, Double_t level,Double_t centrality,ETRDPIDMethod PIDmethod) const {
511     //
512     // Check whether particle is identified as electron assuming a certain electron efficiency level
513     // Only electron and pion hypothesis is taken into account
514     //
515     // Inputs:
516     //         Number of tracklets
517     //         Likelihood values
518     //         Momentum
519     //         Electron efficiency level
520     //
521     // If the function fails when the params are not accessible, the function returns true
522     //
523   if(!fkPIDResponseObject){
524     AliDebug(3,"No PID Param object available");
525     return kTRUE;
526   } 
527   Double_t probEle = like[AliPID::kElectron]/(like[AliPID::kElectron] + like[AliPID::kPion]);
528   Double_t params[4];
529   if(!fkPIDResponseObject->GetThresholdParameters(nTracklets, level, params,centrality,PIDmethod)){
530     AliError("No Params found for the given configuration");
531     return kTRUE;
532   }
533   Double_t threshold = 1. - params[0] - params[1] * p - params[2] * TMath::Exp(-params[3] * p);
534   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
535   return kFALSE;
536 }
537
538 //____________________________________________________________
539 Bool_t AliTRDPIDResponse::SetPIDResponseObject(const AliTRDPIDResponseObject * obj){
540
541     fkPIDResponseObject = obj;
542     if((AliLog::GetDebugLevel("",IsA()->GetName()))>0)fkPIDResponseObject->Print("");
543     return kTRUE;
544 }
545
546
547 //____________________________________________________________    
548 Bool_t AliTRDPIDResponse::SetdEdxParams(const AliTRDdEdxParams * par)
549 {
550   fkTRDdEdxParams = par;
551   return kTRUE;
552 }