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