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