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