AliRoot Core  ee782a0 (ee782a0)
AliTRDPIDResponse.cxx
Go to the documentation of this file.
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 // modifications 29.10. Yvonne Pachmayer <pachmay@physi.uni-heidelberg.de>
26 //
27 #include <TAxis.h>
28 #include <TClass.h>
29 #include <TDirectory.h>
30 #include <TFile.h>
31 #include <TH1.h>
32 #include <TH2.h>
33 #include <TKey.h>
34 #include <TMath.h>
35 #include <TObjArray.h>
36 #include <TROOT.h>
37 #include <TString.h>
38 #include <TSystem.h>
39 
40 #include "AliLog.h"
41 
42 #include "AliVTrack.h"
43 
45 #include "AliTRDPIDResponse.h"
46 //#include "AliTRDTKDInterpolator.h"
47 #include "AliTRDNDFast.h"
48 #include "AliTRDdEdxParams.h"
49 #include "AliExternalTrackParam.h"
50 
51 
52 
53 ClassImp(AliTRDPIDResponse)
54 
55 //____________________________________________________________
57  TObject()
58  ,fkPIDResponseObject(NULL)
59  ,fkTRDdEdxParams(NULL)
60  ,fGainNormalisationFactor(1.)
61  ,fCorrectEta(kFALSE)
62  ,fCorrectCluster(kFALSE)
63  ,fCorrectCentrality(kFALSE)
64  ,fCurrCentrality(-1.)
65  ,fMagField(0.)
66 {
67  //
68  // Default constructor
69  //
70 }
71 
72 //____________________________________________________________
74  TObject(ref)
75  ,fkPIDResponseObject(NULL)
76  ,fkTRDdEdxParams(NULL)
77  ,fGainNormalisationFactor(ref.fGainNormalisationFactor)
78  ,fCorrectEta(kFALSE)
79  ,fCorrectCluster(kFALSE)
80  ,fCorrectCentrality(kFALSE)
81  ,fCurrCentrality(-1.)
82  ,fMagField(0.)
83 {
84  //
85  // Copy constructor
86  //
87 }
88 
89 //____________________________________________________________
91  //
92  // Assignment operator
93  //
94  if(this == &ref) return *this;
95 
96  // Make copy
97  TObject::operator=(ref);
101  fCorrectEta = ref.fCorrectEta;
105  fMagField = ref.fMagField;
106 
107  return *this;
108 }
109 
110 //____________________________________________________________
112  //
113  // Destructor
114  //
115  if(IsOwner()) {
116  delete fkPIDResponseObject;
117  delete fkTRDdEdxParams;
118  delete fhEtaCorr[0];
119  for (Int_t i=0;i<3;i++) delete fhClusterCorr[i];
120  delete fhCentralityCorr[0];
121  }
122 }
123 
124 //____________________________________________________________
125 Bool_t AliTRDPIDResponse::Load(const Char_t * filename){
126  //
127  // Load References into the toolkit
128  //
129  AliDebug(1, "Loading reference histos from root file");
130  TDirectory *owd = gDirectory;// store old working directory
131 
132  if(!filename)
133  filename = Form("%s/STEER/LQ1dRef_v1.root",gSystem->ExpandPathName("$ALICE_ROOT"));
134  TFile *in = TFile::Open(filename);
135  if(!in){
136  AliError("Ref file not available.");
137  return kFALSE;
138  }
139 
140  gROOT->cd();
141  fkPIDResponseObject = dynamic_cast<const AliTRDPIDResponseObject *>(in->Get("TRDPIDResponse")->Clone());
142  in->Close(); delete in;
143  owd->cd();
144  SetBit(kIsOwner, kTRUE);
145  AliDebug(2, Form("Successfully loaded References for %d Momentum bins", fkPIDResponseObject->GetNumberOfMomentumBins()));
146  return kTRUE;
147 }
148 
149 //_________________________________________________________________________
150 Bool_t AliTRDPIDResponse::SetEtaCorrMap(Int_t i,TH2D* hMap)
151 {
152  //
153  // Load map for TRD eta correction (a copy is stored and will be deleted automatically).
154  // If hMap is 0x0,the eta correction will be disabled and kFALSE is returned.
155  // If the map can be set, kTRUE is returned.
156  //
157 
158  // delete fhEtaCorr[0];
159 
160  if (!hMap) {
161  fhEtaCorr[0] = 0x0;
162 
163  return kFALSE;
164  }
165 
166  fhEtaCorr[0] = (TH2D*)(hMap->Clone());
167 
168  return kTRUE;
169 }
170 
171 //_________________________________________________________________________
172 Bool_t AliTRDPIDResponse::SetClusterCorrMap(Int_t i,TH2D* hMap)
173 {
174  //
175  // Load map for TRD cluster correction (a copy is stored and will be deleted automatically).
176  // If hMap is 0x0,the cluster correction will be disabled and kFALSE is returned.
177  // If the map can be set, kTRUE is returned.
178  //
179 
180 
181  if (!hMap) {
182  fhClusterCorr[i] = 0x0;
183 
184  return kFALSE;
185  }
186 
187  fhClusterCorr[i] = (TH2D*)(hMap->Clone());
188 
189  return kTRUE;
190 }
191 
192 //_________________________________________________________________________
193 Bool_t AliTRDPIDResponse::SetCentralityCorrMap(Int_t i,TH2D* hMap)
194 {
195  //
196  // Load map for TRD centrality correction (a copy is stored and will be deleted automatically).
197  // If hMap is 0x0,the centrality correction will be disabled and kFALSE is returned.
198  // If the map can be set, kTRUE is returned.
199  //
200 
201  if (!hMap) {
202  fhCentralityCorr[0] = 0x0;
203 
204  return kFALSE;
205  }
206 
207  fhCentralityCorr[0] = (TH2D*)(hMap->Clone());
208 
209  return kTRUE;
210 }
211 
212 //____________________________________________________________
213 Double_t AliTRDPIDResponse::GetEtaCorrection(const AliVTrack *track, Double_t bg) const
214 {
215  //
216  // eta correction
217  //
218 
219 
220  if (!fhEtaCorr[0]) {
221  AliError(Form("Eta correction requested, but map not initialised for iterator:%i (usually via AliPIDResponse). Returning eta correction factor 1!",1));
222  return 1.;
223 
224  }
225 
226  Double_t fEtaCorFactor=1;
227 
228  Int_t nch = track->GetTRDNchamberdEdx();
229  Int_t iter=0;
230 
231  if (nch < 4) {
232  AliError(Form("Eta correction requested for track with = %i, no map available. Returning default eta correction factor = 1!", nch));
233  return 1.;
234  }
235 
236  Double_t tpctgl= 1.1;
237  tpctgl=track->GetTPCTgl();
238 
239 
240 
241  if ((fhEtaCorr[iter]->GetBinContent(fhEtaCorr[iter]->FindBin(tpctgl,bg)) != 0)) {
242  fEtaCorFactor= fhEtaCorr[iter]->GetBinContent(fhEtaCorr[iter]->FindBin(tpctgl,bg));
243  return fEtaCorFactor;
244  } else
245  {
246  return 1;
247  }
248 
249 
250 }
251 
252 
253 //____________________________________________________________
255 {
256  //
257  // eta correction
258  //
259 
260  Int_t nch = track->GetTRDNchamberdEdx();
261 
262  Double_t fClusterCorFactor=1;
263 
264  if (nch < 4) {
265  AliError(Form("Cluster correction requested for track with = %i, no map available. Returning default cluster correction factor = 1!", nch));
266  return 1.;
267  }
268 
269  Int_t offset =4;
270  const Int_t iter = nch-offset;
271  const Int_t ncls = track->GetTRDNclusterdEdx();
272 
273  if (!fhClusterCorr[iter]) {
274  AliError(Form("Cluster correction requested, but map not initialised for iterator:%i (usually via AliPIDResponse). Returning cluster correction factor 1!",1));
275  return 1.;
276  }
277 
278  if ((fhClusterCorr[iter]->GetBinContent(fhClusterCorr[iter]->FindBin(ncls,bg)) != 0)) {
279  fClusterCorFactor= fhClusterCorr[iter]->GetBinContent(fhClusterCorr[iter]->FindBin(ncls,bg));
280  return fClusterCorFactor;
281  } else
282  {
283  return 1;
284  }
285 
286 }
287 
288 //____________________________________________________________
290 {
291  //
292  // centrality correction
293  //
294 
295 
296  if (!fhCentralityCorr[0]) {
297  // AliInfo(Form("centrality correction requested, but map not initialised for iterator:%i (usually via AliPIDResponse). Returning centrality correction factor 1!",1));
298  return 1.;
299 
300  }
301 
302  Double_t fCentralityCorFactor=1;
303 
304  Int_t nch = track->GetTRDNchamberdEdx();
305  Int_t iter=0;
306 
307  if (nch < 4) {
308  // Ali(Form("Centrality correction requested for track with = %i, no map available. Returning default centrality correction factor = 1!", nch));
309  return 1.;
310  }
311 
312 
313  if(fCurrCentrality<0) return 1.;
314 
315 
316  if ((fhCentralityCorr[iter]->GetBinContent(fhCentralityCorr[iter]->FindBin(fCurrCentrality,bg)) != 0)) {
317  fCentralityCorFactor= fhCentralityCorr[iter]->GetBinContent(fhCentralityCorr[iter]->FindBin(fCurrCentrality,bg));
318  return fCentralityCorFactor;
319  } else
320  {
321  return 1;
322  }
323 
324 
325 }
326 
327 
328 //____________________________________________________________
330 {
331  //
332  //calculate the TRD nSigma
333  //
334 
335  const Double_t badval = -9999;
336  Double_t info[5]; for(Int_t i=0; i<5; i++){info[i]=badval;}
337  const Double_t delta = GetSignalDelta(track, type, kFALSE, fCorrectEta, fCorrectCluster, fCorrectCentrality, info);
338 
339  const Double_t mean = info[0];
340  const Double_t res = info[1];
341  if(res<0){
342  return badval;
343  }
344 
345  const Double_t sigma = mean*res;
346  const Double_t eps = 1e-12;
347  return delta/(sigma + eps);
348 }
349 
350 //____________________________________________________________
351 Double_t AliTRDPIDResponse::GetSignalDelta( const AliVTrack* track, AliPID::EParticleType type, Bool_t ratio/*=kFALSE*/, Bool_t fCorrectEta, Bool_t fCorrectCluster, Bool_t fCorrectCentrality, Double_t *info/*=0x0*/) const
352 {
353  //
354  //calculate the TRD signal difference w.r.t. the expected
355  //output other information in info[]
356  //
357 
358  const Double_t badval = -9999;
359 
360  if(!track){
361  return badval;
362  }
363 
364  Double_t pTRD = 0;
365  Int_t pTRDNorm =0 ;
366  for(Int_t ich=0; ich<6; ich++){
367  if(track->GetTRDmomentum(ich)>0)
368  {
369  pTRD += track->GetTRDmomentum(ich);
370  pTRDNorm++;
371  }
372  }
373 
374  if(pTRDNorm>0)
375  {
376  pTRD/=pTRDNorm;
377  }
378  else return badval;
379 
380  if(!fkTRDdEdxParams){
381  AliDebug(3,"fkTRDdEdxParams null");
382  return -99999;
383  }
384 
385  const Double_t nch = track->GetTRDNchamberdEdx();
386  const Double_t ncls = track->GetTRDNclusterdEdx();
387 
388  // fkTRDdEdxParams->Print();
389 
390  const TVectorF meanvec = fkTRDdEdxParams->GetMeanParameter(type, nch, ncls,fCorrectEta);
391  if(meanvec.GetNrows()==0){
392  return badval;
393  }
394 
395  const TVectorF resvec = fkTRDdEdxParams->GetSigmaParameter(type, nch, ncls,fCorrectEta);
396  if(resvec.GetNrows()==0){
397  return badval;
398  }
399 
400  const Float_t *meanpar = meanvec.GetMatrixArray();
401  const Float_t *respar = resvec.GetMatrixArray();
402 
403  //============================================================================================<<<<<<<<<<<<<
404 
405 
406 
407  const Double_t bg = pTRD/AliPID::ParticleMass(type);
408  const Double_t expsig = MeandEdxTR(&bg, meanpar);
409 
410  if(info){
411  info[0]= expsig;
412  info[1]= ResolutiondEdxTR(&ncls, respar);
413  }
414 
415 
416 
417  const Double_t eps = 1e-10;
418 
419  // eta asymmetry correction
420  Double_t corrFactorEta = 1.0;
421 
422  if (fCorrectEta) {
423  corrFactorEta = GetEtaCorrection(track,bg);
424  }
425 
426  // cluster correction
427  Double_t corrFactorCluster = 1.0;
428  if (fCorrectCluster) {
429  corrFactorCluster = GetClusterCorrection(track,bg);
430  }
431 
432 
433  // centrality correction
434  Double_t corrFactorCentrality = 1.0;
435  if (fCorrectCentrality) {
436  corrFactorCentrality = GetCentralityCorrection(track,bg);
437  }
438 
439 
440  AliDebug(3,Form("TRD trunc PID expected signal %f exp. resolution %f bg %f nch %f ncls %f etcoron/off %i clustercoron/off %i centralitycoron/off %i nsigma %f ratio %f \n",
441  expsig,ResolutiondEdxTR(&ncls, respar),bg,nch,ncls,fCorrectEta,fCorrectCluster,fCorrectCentrality,(corrFactorEta*corrFactorCluster*corrFactorCentrality*track->GetTRDsignal())/(expsig + eps),
442  (corrFactorEta*corrFactorCluster*corrFactorCentrality*track->GetTRDsignal()) - expsig));
443 
444 
445  if(ratio){
446  return (corrFactorEta*corrFactorCluster*corrFactorCentrality*track->GetTRDsignal())/(expsig + eps);
447  }
448  else{
449  return (corrFactorEta*corrFactorCluster*corrFactorCentrality*track->GetTRDsignal()) - expsig;
450  }
451 
452 
453 
454 }
455 
456 
457 Double_t AliTRDPIDResponse::ResolutiondEdxTR(const Double_t * xx, const Float_t * par)
458 {
459  //
460  //resolution
461  //npar=3
462  //
463 
464  const Double_t ncls = xx[0];
465  // return par[0]+par[1]*TMath::Power(ncls, par[2]);
466  return TMath::Sqrt(par[0]*par[0]+par[1]*par[1]/ncls);
467 }
468 
469 Double_t AliTRDPIDResponse::MeandEdxTR(const Double_t * xx, const Float_t * pin)
470 {
471  //
472  //ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
473  //npar = 8 = 3+5
474  //
475 
476  Float_t ptr[4]={0};
477  for(Int_t ii=0; ii<3; ii++){
478  ptr[ii+1]=pin[ii];
479  }
480  return MeanTR(xx,ptr) + MeandEdx(xx,&(pin[3]));
481 }
482 
483 Double_t AliTRDPIDResponse::MeanTR(const Double_t * xx, const Float_t * par)
484 {
485  //
486  //LOGISTIC parametrization for TR, in unit of MIP
487  //npar = 4
488  //
489 
490  const Double_t bg = xx[0];
491  const Double_t gamma = sqrt(1+bg*bg);
492 
493  const Double_t p0 = TMath::Abs(par[1]);
494  const Double_t p1 = TMath::Abs(par[2]);
495  const Double_t p2 = TMath::Abs(par[3]);
496 
497  const Double_t zz = TMath::Log(gamma);
498  const Double_t tryield = p0/( 1 + TMath::Exp(-p1*(zz-p2)) );
499 
500  return par[0]+tryield;
501 }
502 
503 Double_t AliTRDPIDResponse::MeandEdx(const Double_t * xx, const Float_t * par)
504 {
505  //
506  //ALEPH parametrization for dEdx
507  //npar = 5
508  //
509 
510  const Double_t bg = xx[0];
511  const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
512 
513  const Double_t p0 = TMath::Abs(par[0]);
514  const Double_t p1 = TMath::Abs(par[1]);
515 
516  const Double_t p2 = TMath::Abs(par[2]);
517 
518  const Double_t p3 = TMath::Abs(par[3]);
519  const Double_t p4 = TMath::Abs(par[4]);
520 
521  const Double_t aa = TMath::Power(beta, p3);
522 
523  const Double_t bb = TMath::Log( p2 + TMath::Power(1./bg, p4) );
524 
525  return (p1-aa-bb)*p0/aa;
526 
527 }
528 
529 
530 //____________________________________________________________
531 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 AliVTrack *track) const
532 {
533  //
534  // Calculate TRD likelihood values for the track based on dedx and
535  // momentum values. The likelihoods are calculated by query the
536  // reference data depending on the PID method selected
537  //
538  // Input parameter :
539  // n - number of dedx slices/chamber
540  // dedx - array of dedx slices organized layer wise
541  // p - array of momentum measurements organized layer wise
542  //
543  // Return parameters
544  // prob - probabilities allocated by TRD for particle specis
545  // kNorm - switch to normalize probabilities to 1. By default true. If false return not normalized prob.
546  //
547  // Return value
548  // number of tracklets used for PID, 0 if no PID
549  //
550  AliDebug(3,Form(" Response for PID method: %d",PIDmethod));
551  Int_t iCharge=AliPID::kNoCharge;
552  if (track!=NULL){
553  if (track->Charge()>0){
554  iCharge=AliPID::kPosCharge;
555  }
556  else {
557  iCharge=AliPID::kNegCharge;
558  }
559  }
560  if(!fkPIDResponseObject){
561  AliDebug(3,"Missing reference data. PID calculation not possible.");
562  return 0;
563  }
564 
565  for(Int_t is(AliPID::kSPECIES); is--;) prob[is]=.2;
566  Double_t prLayer[AliPID::kSPECIES];
567  Double_t dE[10], s(0.);
568  Int_t ntrackletsPID=0;
569  for(Int_t il(kNlayer); il--;){
570  memset(prLayer, 0, AliPID::kSPECIES*sizeof(Double_t));
571  if(!CookdEdx(n, &dedx[il*n], &dE[0],PIDmethod)) continue;
572  s=0.;
573  Bool_t filled=kTRUE;
574  for(Int_t is(AliPID::kSPECIES); is--;){
575  //if((PIDmethod==kLQ2D)&&(!(is==0||is==2)))continue;
576  if((dE[0] > 0.) && (p[il] > 0.)) prLayer[is] = GetProbabilitySingleLayer(is, p[il], &dE[0],PIDmethod,iCharge);
577  AliDebug(3, Form("Probability for Species %d in Layer %d: %e", is, il, prLayer[is]));
578  if(prLayer[is]<1.e-30){
579  AliDebug(2, Form("Null for species %d species prob layer[%d].",is,il));
580  filled=kFALSE;
581  break;
582  }
583  s+=prLayer[is];
584  }
585  if(!filled){
586  continue;
587  }
588  if(s<1.e-30){
589  AliDebug(2, Form("Null all species prob layer[%d].", il));
590  continue;
591  }
592  for(Int_t is(AliPID::kSPECIES); is--;){
593  if(kNorm) prLayer[is] /= s; // probability in each layer for each particle species normalized to the sum of probabilities for given layer
594  prob[is] *= prLayer[is]; // multiply single layer probabilities to get probability for complete track
595  }
596  ntrackletsPID++;
597  }
598  if(!kNorm) return ntrackletsPID;
599 
600  s=0.;
601  // sum probabilities for all considered particle species
602  for(Int_t is(AliPID::kSPECIES); is--;) { s+=prob[is];
603  }
604  if(s<1.e-30){
605  AliDebug(2, "Null total prob.");
606  return 0;
607  }
608  // norm to the summed probability (default values: s=1 prob[is]=0.2)
609  for(Int_t is(AliPID::kSPECIES); is--;){ prob[is]/=s; }
610  return ntrackletsPID;
611 }
612 
613 //____________________________________________________________
614 Double_t AliTRDPIDResponse::GetProbabilitySingleLayer(Int_t species, Double_t plocal, Double_t *dEdx,ETRDPIDMethod PIDmethod, Int_t Charge) const {
615  //
616  // Get the non-normalized probability for a certain particle species as coming
617  // from the reference histogram
618  // Interpolation between momentum bins
619  //
620  AliDebug(1, Form("Make Probability for Species %d with Momentum %f", species, plocal));
621 
622  Double_t probLayer = 0.;
623 
624  Float_t pLower, pUpper;
625 
626  AliTRDNDFast *refUpper = dynamic_cast<AliTRDNDFast *>(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,PIDmethod,Charge)),
627  *refLower = dynamic_cast<AliTRDNDFast *>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,PIDmethod, Charge));
628  if ((!refLower)&&(!refUpper)&&(Charge!=AliPID::kNoCharge)){
629  AliDebug(3,Form("No references available for Charge; References for both Charges used"));
630  refUpper = dynamic_cast<AliTRDNDFast *>(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,PIDmethod,AliPID::kNoCharge));
631  refLower = dynamic_cast<AliTRDNDFast *>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,PIDmethod,AliPID::kNoCharge));
632  }
633  // Do Interpolation exept for underflow and overflow
634  if(refLower && refUpper){
635  Double_t probLower = refLower->Eval(dEdx);
636  Double_t probUpper = refUpper->Eval(dEdx);
637 
638  probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
639  } else if(refLower){
640  // underflow
641  probLayer = refLower->Eval(dEdx);
642  } else if(refUpper){
643  // overflow
644  probLayer = refUpper->Eval(dEdx);
645  } else {
646  AliDebug(3,"No references available");
647  }
648 
649  switch(PIDmethod){
650  case kLQ2D: // 2D LQ
651  {
652  AliDebug(1,Form("Eval 2D Q0 %f Q1 %f P %e ",dEdx[0],dEdx[1],probLayer));
653  }
654  break;
655  case kLQ1D: // 1D LQ
656  {
657  AliDebug(1, Form("Eval 1D dEdx %f Probability %e", dEdx[0],probLayer));
658  }
659  break;
660  case kLQ3D: // 3D LQ
661  {
662  AliDebug(1, Form("Eval 1D dEdx %f %f %f Probability %e", dEdx[0],dEdx[1],dEdx[2],probLayer));
663  }
664  break;
665  case kLQ7D: // 7D LQ
666  {
667  AliDebug(1, Form("Eval 1D dEdx %f %f %f %f %f %f %f Probability %e", dEdx[0],dEdx[1],dEdx[2],dEdx[3],dEdx[4],dEdx[5],dEdx[6],probLayer));
668  }
669  break;
670  default:
671  break;
672  }
673 
674  return probLayer;
675 
676  /* old implementation
677 
678 switch(PIDmethod){
679 case kNN: // NN
680  break;
681  case kLQ2D: // 2D LQ
682  {
683  if(species==0||species==2){ // references only for electrons and pions
684  Double_t error = 0.;
685  Double_t point[kNslicesLQ2D];
686  for(Int_t idim=0;idim<kNslicesLQ2D;idim++){point[idim]=dEdx[idim];}
687 
688  AliTRDTKDInterpolator *refUpper = dynamic_cast<AliTRDTKDInterpolator *>(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,kLQ2D)),
689  *refLower = dynamic_cast<AliTRDTKDInterpolator *>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,kLQ2D));
690  // Do Interpolation exept for underflow and overflow
691  if(refLower && refUpper){
692  Double_t probLower=0,probUpper=0;
693  refLower->Eval(point,probLower,error);
694  refUpper->Eval(point,probUpper,error);
695  probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
696  } else if(refLower){
697  // underflow
698  refLower->Eval(point,probLayer,error);
699  } else if(refUpper){
700  // overflow
701  refUpper->Eval(point,probLayer,error);
702  } else {
703  AliError("No references available");
704  }
705  AliDebug(2,Form("Eval 2D Q0 %f Q1 %f P %e Err %e",point[0],point[1],probLayer,error));
706  }
707  }
708  break;
709  case kLQ1D: // 1D LQ
710  {
711  TH1 *refUpper = dynamic_cast<TH1 *>(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,kLQ1D)),
712  *refLower = dynamic_cast<TH1 *>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,kLQ1D));
713  // Do Interpolation exept for underflow and overflow
714  if(refLower && refUpper){
715  Double_t probLower = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx[0]));
716  Double_t probUpper = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx[0]));
717 
718  probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
719  } else if(refLower){
720  // underflow
721  probLayer = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx[0]));
722  } else if(refUpper){
723  // overflow
724  probLayer = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx[0]));
725  } else {
726  AliError("No references available");
727  }
728  AliDebug(1, Form("Eval 1D dEdx %f Probability %e", dEdx[0],probLayer));
729  }
730  break;
731  default:
732  break;
733  }
734  return probLayer;
735  */
736 
737 }
738 
739 //____________________________________________________________
741  //
742  // Make Deep Copy of the Reference Histograms
743  //
744  if(!fkPIDResponseObject || IsOwner()) return;
746  fkPIDResponseObject = dynamic_cast<const AliTRDPIDResponseObject *>(tmp->Clone());
747  SetBit(kIsOwner, kTRUE);
748 }
749 
750 //____________________________________________________________
751 Bool_t AliTRDPIDResponse::CookdEdx(Int_t nSlice, const Double_t * const in, Double_t *out,ETRDPIDMethod PIDmethod) const
752 {
753  //
754  // Recalculate dE/dx
755  // removed missing slices cut, detailed checks see presentation in TRD meeting: https://indico.cern.ch/event/506345/contribution/3/attachments/1239069/1821088/TRDPID_missingslices_charge.pdf
756  //
757 
758  switch(PIDmethod){
759  case kNN: // NN
760  break;
761  case kLQ2D: // 2D LQ
762  out[0]=0;
763  out[1]=0;
764  for(Int_t islice = 0; islice < nSlice; islice++){
765  // if(in[islice]<=0){out[0]=0;out[1]=0;return kFALSE;} // Require that all slices are filled
766 
767  if(islice<kNsliceQ0LQ2D)out[0]+= in[islice];
768  else out[1]+= in[islice];
769  }
770  // normalize signal to number of slices
771  out[0]*=1./Double_t(kNsliceQ0LQ2D);
772  out[1]*=1./Double_t(nSlice-kNsliceQ0LQ2D);
773  if(out[0] <= 0) return kFALSE;
774  if(out[1] <= 0) return kFALSE;
775  AliDebug(3,Form("CookdEdx Q0 %f Q1 %f",out[0],out[1]));
776  break;
777  case kLQ1D: // 1D LQ
778  out[0]= 0.;
779  for(Int_t islice = 0; islice < nSlice; islice++) {
780  // if(in[islice] > 0) out[0] += in[islice] * fGainNormalisationFactor; // Protect against negative values for slices having no dE/dx information
781  if(in[islice] > 0) out[0] += in[islice]; // no neg dE/dx values
782  }
783  out[0]*=1./Double_t(kNsliceQ0LQ1D);
784  if(out[0] <= 0) return kFALSE;
785  AliDebug(3,Form("CookdEdx dEdx %f",out[0]));
786  break;
787  case kLQ3D: // 3D LQ
788  out[0]=0;
789  out[1]=0;
790  out[2]=0;
791  for(Int_t islice = 0; islice < nSlice; islice++){
792  // if(in[islice]<=0){out[0]=0;out[1]=0;out[2]=0;return kFALSE;} // Require that all slices are filled
793  if(islice<kNsliceQ0LQ3D)out[0]+= in[islice];
794  out[1]=(in[3]+in[4]);
795  out[2]=(in[5]+in[6]);
796  }
797  // normalize signal to number of slices
798  out[0]*=1./Double_t(kNsliceQ0LQ3D);
799  out[1]*=1./2.;
800  out[2]*=1./2.;
801  if(out[0] <= 0) return kFALSE;
802  if(out[1] <= 0) return kFALSE;
803  if(out[2] <= 0) return kFALSE;
804  AliDebug(3,Form("CookdEdx Q0 %f Q1 %f Q2 %f",out[0],out[1],out[2]));
805  break;
806  case kLQ7D: // 7D LQ
807  for(Int_t i=0;i<nSlice;i++) {out[i]=0;}
808  for(Int_t islice = 0; islice < nSlice; islice++){
809  if(in[islice]<=0){
810  for(Int_t i=0;i<8;i++){
811  out[i]=0;
812  }
813  return kFALSE;} // Require that all slices are filled
814  out[islice]=in[islice];
815  }
816  for(Int_t i=0;i<nSlice;i++) {if(out[i]<=0) return kFALSE; }
817  AliDebug(3,Form("CookdEdx Q0 %f Q1 %f Q2 %f Q3 %f Q4 %f Q5 %f Q6 %f Q7 %f",out[0],out[1],out[2],out[3],out[4],out[5],out[6],out[7]));
818  break;
819 
820  default:
821  return kFALSE;
822  }
823  return kTRUE;
824 }
825 
826 //____________________________________________________________
827 Bool_t AliTRDPIDResponse::IdentifiedAsElectron(Int_t nTracklets, const Double_t *like, Double_t p, Double_t level,Double_t centrality,ETRDPIDMethod PIDmethod, const AliVTrack *vtrack) const {
828  //
829  // Check whether particle is identified as electron assuming a certain electron efficiency level
830  // Only electron and pion hypothesis is taken into account
831  //
832  // Inputs:
833  // Number of tracklets
834  // Likelihood values
835  // Momentum
836  // Electron efficiency level
837  //
838  // If the function fails when the params are not accessible, the function returns true
839  //
840  Int_t iCharge=AliPID::kNoCharge;
841  if (vtrack!=NULL){
842  Int_t vTrCharge=vtrack->Charge();
843  if (vTrCharge>0){
844  iCharge=AliPID::kPosCharge;
845  }
846  else{
847  iCharge=AliPID::kNegCharge;
848  }
849  }
850  if(!fkPIDResponseObject){
851  AliDebug(3,"No PID Param object available");
852  return kTRUE;
853  }
854  Double_t probEle = like[AliPID::kElectron]/(like[AliPID::kElectron] + like[AliPID::kPion]);
855  AliDebug(3,Form("probabilities like %f %f %f \n",probEle,like[AliPID::kElectron],like[AliPID::kPion]));
856  Double_t params[4];
857  if(!fkPIDResponseObject->GetThresholdParameters(nTracklets, level, params,centrality,PIDmethod,iCharge)){
858  AliError("No Params found for the given configuration with chosen Charge");
859  AliError("Using Parameters for both charges");
860  if((iCharge!=AliPID::kNoCharge)&&(!fkPIDResponseObject->GetThresholdParameters(nTracklets, level, params,centrality,PIDmethod,AliPID::kNoCharge))){
861  AliError("No Params found for the given configuration with charge 0");
862  return kTRUE;
863  }
864  if(iCharge==AliPID::kNoCharge){
865  return kTRUE;
866  }
867  }
868 
869  Double_t threshold = 1. - params[0] - params[1] * p - params[2] * TMath::Exp(-params[3] * p);
870  AliDebug(3,Form("is ident details %i %f %f %i %f %f %f %f \n",nTracklets, level, centrality,PIDmethod,probEle, threshold,TMath::Min(threshold, 0.99),TMath::Max(TMath::Min(threshold, 0.99), 0.2)));
871  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
872  return kFALSE;
873 }
874 
875 //____________________________________________________________
877 
878  fkPIDResponseObject = obj;
879  if((AliLog::GetDebugLevel("",IsA()->GetName()))>0 && obj) fkPIDResponseObject->Print("");
880  return kTRUE;
881 }
882 
883 
884 //____________________________________________________________
886 {
887  fkTRDdEdxParams = par;
888  return kTRUE;
889 }
TFile * Open(const char *filename, Long64_t &nevents)
virtual UChar_t GetTRDNchamberdEdx() const
Definition: AliVTrack.h:178
virtual Double_t GetTPCTgl() const
Definition: AliVTrack.h:170
virtual Short_t Charge() const =0
Bool_t SetPIDResponseObject(const AliTRDPIDResponseObject *obj)
Bool_t IdentifiedAsElectron(Int_t nTracklets, const Double_t *like, Double_t p, Double_t level, Double_t centrality=-1, ETRDPIDMethod PIDmethod=kLQ1D, const AliVTrack *vtrack=NULL) const
TH2D * fhClusterCorr[3]
switch for cluster correction
static Float_t ParticleMass(Int_t iType)
Definition: AliPID.h:56
TObject * GetUpperReference(AliPID::EParticleType spec, Float_t p, Float_t &pUpper, AliTRDPIDResponse::ETRDPIDMethod method=AliTRDPIDResponse::kLQ1D, Int_t Charge=0) const
const TVectorF & GetSigmaParameter(const Int_t itype, const Int_t nch, const Int_t ncls, const Bool_t etaCorrection) const
EParticleType
Definition: AliPID.h:27
TROOT * gROOT
virtual void Print(Option_t *opt) const
Float_t p[]
Definition: kNNTest.C:133
Bool_t fCorrectCentrality
Map for TRD eta correction.
Double_t GetClusterCorrection(const AliVTrack *track, Double_t bg) const
Double_t Eval(Double_t *point) const
AliTPCfastTrack * track
Bool_t fCorrectEta
Gain normalisation factor.
Bool_t SetClusterCorrMap(Int_t n, TH2D *hMapn)
Double_t GetEtaCorrection(const AliVTrack *track, Double_t bg) const
AliExternalInfo info
Bool_t bg
Definition: RunAnaESD.C:6
static Double_t MeandEdx(const Double_t *xx, const Float_t *par)
TObject * GetLowerReference(AliPID::EParticleType spec, Float_t p, Float_t &pLower, AliTRDPIDResponse::ETRDPIDMethod method=AliTRDPIDResponse::kLQ1D, Int_t Charge=0) const
AliTRDPIDResponse & operator=(const AliTRDPIDResponse &ref)
Bool_t GetThresholdParameters(Int_t ntracklets, Double_t efficiency, Double_t *params, Double_t centrality=-1, AliTRDPIDResponse::ETRDPIDMethod method=AliTRDPIDResponse::kLQ1D, Int_t charge=0) const
Bool_t Load(const Char_t *filename=NULL)
Bool_t SetCentralityCorrMap(Int_t n, TH2D *hMapn)
const AliTRDPIDResponseObject * fkPIDResponseObject
Double_t GetSignalDelta(const AliVTrack *track, AliPID::EParticleType type, Bool_t ratio=kFALSE, Bool_t fCorrectEta=kFALSE, Bool_t fCorrectCluster=kFALSE, Bool_t fCorrectCentrality=kFALSE, Double_t *info=0x0) const
TH2D * fhEtaCorr[1]
switch for eta correction
Double_t GetCentralityCorrection(const AliVTrack *track, Double_t bg) const
Double_t fGainNormalisationFactor
parametrisation for truncated mean
Double_t GetProbabilitySingleLayer(Int_t species, Double_t plocal, Double_t *dEdx, ETRDPIDMethod PIDmethod=kLQ1D, Int_t Charge=0) const
Int_t GetResponse(Int_t n, const Double_t *const dedx, const Float_t *const p, Double_t prob[AliPID::kSPECIES], ETRDPIDMethod PIDmethod=kLQ1D, Bool_t kNorm=kTRUE, const AliVTrack *track=NULL) const
static Double_t MeanTR(const Double_t *xx, const Float_t *par)
Bool_t SetdEdxParams(const AliTRDdEdxParams *par)
Double_t GetNumberOfSigmas(const AliVTrack *track, AliPID::EParticleType type, Bool_t fCorrectEta, Bool_t fCorrectCluster, Bool_t fCorrectCentrality) const
static Double_t ResolutiondEdxTR(const Double_t *xx, const Float_t *par)
#define AliDebug(logLevel, message)
Definition: AliLog.h:300
virtual UChar_t GetTRDNclusterdEdx() const
Definition: AliVTrack.h:179
Bool_t SetEtaCorrMap(Int_t n, TH2D *hMapn)
const TVectorF & GetMeanParameter(const Int_t itype, const Int_t nch, const Int_t ncls, const Bool_t etaCorrection) const
const AliTRDdEdxParams * fkTRDdEdxParams
PID References and Params.
Bool_t CookdEdx(Int_t nSlice, const Double_t *const in, Double_t *out, ETRDPIDMethod PIDmethod=kLQ1D) const
virtual Double_t GetTRDsignal() const
Definition: AliVTrack.h:177
static Double_t MeandEdxTR(const Double_t *xx, const Float_t *par)
Double_t fCurrCentrality
Map for TRD centrality correction.
static Int_t GetDebugLevel(const char *module, const char *className)
Definition: AliLog.cxx:843
#define AliError(message)
Definition: AliLog.h:591
Int_t GetNumberOfMomentumBins(AliTRDPIDResponse::ETRDPIDMethod method=AliTRDPIDResponse::kLQ1D) const
Bool_t fCorrectCluster
Map for TRD eta correction.
void res(Char_t i)
Definition: Resolution.C:2
Bool_t IsOwner() const
TH2D * fhCentralityCorr[1]
switch for centrality correction
virtual Double_t GetTRDmomentum(Int_t, Double_t *=0x0) const
Definition: AliVTrack.h:196