Line data Source code
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 :
44 : #include "AliTRDPIDResponseObject.h"
45 : #include "AliTRDPIDResponse.h"
46 : //#include "AliTRDTKDInterpolator.h"
47 : #include "AliTRDNDFast.h"
48 : #include "AliTRDdEdxParams.h"
49 : #include "AliExternalTrackParam.h"
50 :
51 :
52 :
53 176 : ClassImp(AliTRDPIDResponse)
54 :
55 : //____________________________________________________________
56 : AliTRDPIDResponse::AliTRDPIDResponse():
57 12 : TObject()
58 12 : ,fkPIDResponseObject(NULL)
59 12 : ,fkTRDdEdxParams(NULL)
60 12 : ,fGainNormalisationFactor(1.)
61 12 : ,fCorrectEta(kFALSE)
62 12 : ,fMagField(0.)
63 60 : {
64 : //
65 : // Default constructor
66 : //
67 24 : }
68 :
69 : //____________________________________________________________
70 : AliTRDPIDResponse::AliTRDPIDResponse(const AliTRDPIDResponse &ref):
71 0 : TObject(ref)
72 0 : ,fkPIDResponseObject(NULL)
73 0 : ,fkTRDdEdxParams(NULL)
74 0 : ,fGainNormalisationFactor(ref.fGainNormalisationFactor)
75 0 : ,fCorrectEta(kFALSE)
76 0 : ,fMagField(0.)
77 0 : {
78 : //
79 : // Copy constructor
80 : //
81 0 : }
82 :
83 : //____________________________________________________________
84 : AliTRDPIDResponse &AliTRDPIDResponse::operator=(const AliTRDPIDResponse &ref){
85 : //
86 : // Assignment operator
87 : //
88 0 : if(this == &ref) return *this;
89 :
90 : // Make copy
91 0 : TObject::operator=(ref);
92 0 : fGainNormalisationFactor = ref.fGainNormalisationFactor;
93 0 : fkPIDResponseObject = ref.fkPIDResponseObject;
94 0 : fkTRDdEdxParams = ref.fkTRDdEdxParams;
95 0 : fCorrectEta = ref.fCorrectEta;
96 0 : fMagField = ref.fMagField;
97 :
98 0 : return *this;
99 0 : }
100 :
101 : //____________________________________________________________
102 40 : AliTRDPIDResponse::~AliTRDPIDResponse(){
103 : //
104 : // Destructor
105 : //
106 20 : if(IsOwner()) {
107 0 : delete fkPIDResponseObject;
108 0 : delete fkTRDdEdxParams;
109 0 : delete fhEtaCorr[0];
110 : }
111 20 : }
112 :
113 : //____________________________________________________________
114 : Bool_t AliTRDPIDResponse::Load(const Char_t * filename){
115 : //
116 : // Load References into the toolkit
117 : //
118 0 : AliDebug(1, "Loading reference histos from root file");
119 0 : TDirectory *owd = gDirectory;// store old working directory
120 :
121 0 : if(!filename)
122 0 : filename = Form("%s/STEER/LQ1dRef_v1.root",gSystem->ExpandPathName("$ALICE_ROOT"));
123 0 : TFile *in = TFile::Open(filename);
124 0 : if(!in){
125 0 : AliError("Ref file not available.");
126 0 : return kFALSE;
127 : }
128 :
129 0 : gROOT->cd();
130 0 : fkPIDResponseObject = dynamic_cast<const AliTRDPIDResponseObject *>(in->Get("TRDPIDResponse")->Clone());
131 0 : in->Close(); delete in;
132 0 : owd->cd();
133 0 : SetBit(kIsOwner, kTRUE);
134 0 : AliDebug(2, Form("Successfully loaded References for %d Momentum bins", fkPIDResponseObject->GetNumberOfMomentumBins()));
135 0 : return kTRUE;
136 0 : }
137 :
138 : //_________________________________________________________________________
139 : Bool_t AliTRDPIDResponse::SetEtaCorrMap(Int_t i,TH2D* hMap)
140 : {
141 : //
142 : // Load map for TRD eta correction (a copy is stored and will be deleted automatically).
143 : // If hMap is 0x0,the eta correction will be disabled and kFALSE is returned.
144 : // If the map can be set, kTRUE is returned.
145 : //
146 :
147 : // delete fhEtaCorr[0];
148 :
149 0 : if (!hMap) {
150 0 : fhEtaCorr[0] = 0x0;
151 :
152 0 : return kFALSE;
153 : }
154 :
155 0 : fhEtaCorr[0] = (TH2D*)(hMap->Clone());
156 :
157 0 : return kTRUE;
158 0 : }
159 :
160 : //____________________________________________________________
161 : Double_t AliTRDPIDResponse::GetEtaCorrection(const AliVTrack *track, Double_t bg) const
162 : {
163 : //
164 : // eta correction
165 : //
166 :
167 :
168 0 : if (!fhEtaCorr[0]) {
169 0 : AliError(Form("Eta correction requested, but map not initialised for iterator:%i (usually via AliPIDResponse). Returning eta correction factor 1!",1));
170 0 : return 1.;
171 :
172 : }
173 :
174 : Double_t fEtaCorFactor=1;
175 :
176 :
177 0 : Int_t nch = track->GetTRDNchamberdEdx();
178 :
179 : const Int_t iter = 0;
180 :
181 : if (iter < 0) {
182 : AliError(Form("Eta correction requested for track with = %i, no map available. Returning default eta correction factor = 1!", nch));
183 : return 1.;
184 : }
185 :
186 :
187 : // Get eta from propagation of outer params to TRD
188 : const AliExternalTrackParam *tempparam = NULL;
189 : Double_t etalayer=-99;
190 0 : if(track->GetOuterParam()) tempparam = track->GetOuterParam();
191 0 : else if(track->GetInnerParam()) tempparam = track->GetInnerParam();
192 :
193 0 : if(tempparam) {
194 0 : AliExternalTrackParam param(*tempparam);
195 0 : param.PropagateTo(288.43,fMagField); // hardwired number where TRD begins
196 0 : etalayer= param.Eta();
197 0 : } else return 1.;
198 :
199 :
200 :
201 0 : if((TMath::Abs(etalayer)>0.9)||(bg<0.3)||(bg>1e4)){
202 0 : return 1;
203 : } else {
204 0 : if ((fhEtaCorr[iter]->GetBinContent(fhEtaCorr[iter]->FindBin(etalayer,bg)) != 0)) {
205 0 : fEtaCorFactor= fhEtaCorr[iter]->GetBinContent(fhEtaCorr[iter]->FindBin(etalayer,bg));
206 0 : return fEtaCorFactor;
207 : } else
208 : {
209 0 : return 1;
210 : }
211 : }
212 0 : }
213 :
214 :
215 :
216 : //____________________________________________________________
217 : Double_t AliTRDPIDResponse::GetNumberOfSigmas(const AliVTrack *track, AliPID::EParticleType type, Bool_t fCorrectEta) const
218 : {
219 : //
220 : //calculate the TRD nSigma
221 : //
222 :
223 : const Double_t badval = -9999;
224 0 : Double_t info[5]; for(int i=0; i<5; i++){info[i]=badval;}
225 :
226 0 : const Double_t delta = GetSignalDelta(track, type, kFALSE, fCorrectEta, info);
227 :
228 0 : const Double_t mean = info[0];
229 0 : const Double_t res = info[1];
230 0 : if(res<0){
231 0 : return badval;
232 : }
233 :
234 0 : const Double_t sigma = mean*res;
235 : const Double_t eps = 1e-12;
236 0 : return delta/(sigma + eps);
237 0 : }
238 :
239 : //____________________________________________________________
240 : Double_t AliTRDPIDResponse::GetSignalDelta( const AliVTrack* track, AliPID::EParticleType type, Bool_t ratio/*=kFALSE*/, Bool_t fCorrectEta, Double_t *info/*=0x0*/) const
241 : {
242 : //
243 : //calculate the TRD signal difference w.r.t. the expected
244 : //output other information in info[]
245 : //
246 :
247 : const Double_t badval = -9999;
248 :
249 0 : if(!track){
250 0 : return badval;
251 : }
252 :
253 : Double_t pTRD = 0;
254 : Int_t pTRDNorm =0 ;
255 0 : for(Int_t ich=0; ich<6; ich++){
256 0 : if(track->GetTRDmomentum(ich)>0)
257 : {
258 0 : pTRD += track->GetTRDmomentum(ich);
259 0 : pTRDNorm++;
260 0 : }
261 : }
262 :
263 0 : if(pTRDNorm>0)
264 : {
265 0 : pTRD/=pTRDNorm;
266 : }
267 0 : else return badval;
268 :
269 0 : if(!fkTRDdEdxParams){
270 0 : AliDebug(3,"fkTRDdEdxParams null");
271 0 : return -99999;
272 : }
273 :
274 0 : const Double_t nch = track->GetTRDNchamberdEdx();
275 0 : const Double_t ncls = track->GetTRDNclusterdEdx();
276 :
277 : // fkTRDdEdxParams->Print();
278 :
279 0 : const TVectorF meanvec = fkTRDdEdxParams->GetMeanParameter(type, nch, ncls,fCorrectEta);
280 0 : if(meanvec.GetNrows()==0){
281 0 : return badval;
282 : }
283 :
284 0 : const TVectorF resvec = fkTRDdEdxParams->GetSigmaParameter(type, nch, ncls,fCorrectEta);
285 0 : if(resvec.GetNrows()==0){
286 0 : return badval;
287 : }
288 :
289 0 : const Float_t *meanpar = meanvec.GetMatrixArray();
290 0 : const Float_t *respar = resvec.GetMatrixArray();
291 :
292 : //============================================================================================<<<<<<<<<<<<<
293 :
294 :
295 :
296 0 : const Double_t bg = pTRD/AliPID::ParticleMass(type);
297 0 : const Double_t expsig = MeandEdxTR(&bg, meanpar);
298 :
299 0 : if(info){
300 0 : info[0]= expsig;
301 0 : info[1]= ResolutiondEdxTR(&ncls, respar);
302 0 : }
303 :
304 :
305 :
306 : const Double_t eps = 1e-10;
307 :
308 : // eta asymmetry correction
309 : Double_t corrFactorEta = 1.0;
310 :
311 0 : if (fCorrectEta) {
312 0 : corrFactorEta = GetEtaCorrection(track,bg);
313 0 : }
314 :
315 0 : AliDebug(3,Form("TRD trunc PID expected signal %f exp. resolution %f bg %f nch %f ncls %f etcoron/off %i nsigma %f ratio %f \n",expsig,ResolutiondEdxTR(&ncls, respar),bg,nch,ncls,fCorrectEta,(corrFactorEta*track->GetTRDsignal())/(expsig + eps),(corrFactorEta*track->GetTRDsignal()) - expsig));
316 :
317 :
318 0 : if(ratio){
319 0 : return (corrFactorEta*track->GetTRDsignal())/(expsig + eps);
320 : }
321 : else{
322 0 : return (corrFactorEta*track->GetTRDsignal()) - expsig;
323 : }
324 :
325 :
326 :
327 0 : }
328 :
329 :
330 : Double_t AliTRDPIDResponse::ResolutiondEdxTR(const Double_t * xx, const Float_t * par)
331 : {
332 : //
333 : //resolution
334 : //npar=3
335 : //
336 :
337 0 : const Double_t ncls = xx[0];
338 : // return par[0]+par[1]*TMath::Power(ncls, par[2]);
339 0 : return TMath::Sqrt(par[0]*par[0]+par[1]*par[1]/ncls);
340 : }
341 :
342 : Double_t AliTRDPIDResponse::MeandEdxTR(const Double_t * xx, const Float_t * pin)
343 : {
344 : //
345 : //ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
346 : //npar = 8 = 3+5
347 : //
348 :
349 0 : Float_t ptr[4]={0};
350 0 : for(int ii=0; ii<3; ii++){
351 0 : ptr[ii+1]=pin[ii];
352 : }
353 0 : return MeanTR(xx,ptr) + MeandEdx(xx,&(pin[3]));
354 0 : }
355 :
356 : Double_t AliTRDPIDResponse::MeanTR(const Double_t * xx, const Float_t * par)
357 : {
358 : //
359 : //LOGISTIC parametrization for TR, in unit of MIP
360 : //npar = 4
361 : //
362 :
363 0 : const Double_t bg = xx[0];
364 0 : const Double_t gamma = sqrt(1+bg*bg);
365 :
366 0 : const Double_t p0 = TMath::Abs(par[1]);
367 0 : const Double_t p1 = TMath::Abs(par[2]);
368 0 : const Double_t p2 = TMath::Abs(par[3]);
369 :
370 0 : const Double_t zz = TMath::Log(gamma);
371 0 : const Double_t tryield = p0/( 1 + TMath::Exp(-p1*(zz-p2)) );
372 :
373 0 : return par[0]+tryield;
374 : }
375 :
376 : Double_t AliTRDPIDResponse::MeandEdx(const Double_t * xx, const Float_t * par)
377 : {
378 : //
379 : //ALEPH parametrization for dEdx
380 : //npar = 5
381 : //
382 :
383 0 : const Double_t bg = xx[0];
384 0 : const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
385 :
386 0 : const Double_t p0 = TMath::Abs(par[0]);
387 0 : const Double_t p1 = TMath::Abs(par[1]);
388 :
389 0 : const Double_t p2 = TMath::Abs(par[2]);
390 :
391 0 : const Double_t p3 = TMath::Abs(par[3]);
392 0 : const Double_t p4 = TMath::Abs(par[4]);
393 :
394 0 : const Double_t aa = TMath::Power(beta, p3);
395 :
396 0 : const Double_t bb = TMath::Log( p2 + TMath::Power(1./bg, p4) );
397 :
398 0 : return (p1-aa-bb)*p0/aa;
399 :
400 : }
401 :
402 :
403 : //____________________________________________________________
404 : 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
405 : {
406 : //
407 : // Calculate TRD likelihood values for the track based on dedx and
408 : // momentum values. The likelihoods are calculated by query the
409 : // reference data depending on the PID method selected
410 : //
411 : // Input parameter :
412 : // n - number of dedx slices/chamber
413 : // dedx - array of dedx slices organized layer wise
414 : // p - array of momentum measurements organized layer wise
415 : //
416 : // Return parameters
417 : // prob - probabilities allocated by TRD for particle specis
418 : // kNorm - switch to normalize probabilities to 1. By default true. If false return not normalized prob.
419 : //
420 : // Return value
421 : // number of tracklets used for PID, 0 if no PID
422 : //
423 312 : AliDebug(3,Form(" Response for PID method: %d",PIDmethod));
424 :
425 :
426 104 : if(!fkPIDResponseObject){
427 0 : AliDebug(3,"Missing reference data. PID calculation not possible.");
428 0 : return 0;
429 : }
430 :
431 1248 : for(Int_t is(AliPID::kSPECIES); is--;) prob[is]=.2;
432 104 : Double_t prLayer[AliPID::kSPECIES];
433 104 : Double_t dE[10], s(0.);
434 : Int_t ntrackletsPID=0;
435 832 : for(Int_t il(kNlayer); il--;){
436 624 : memset(prLayer, 0, AliPID::kSPECIES*sizeof(Double_t));
437 624 : if(!CookdEdx(n, &dedx[il*n], &dE[0],PIDmethod)) continue;
438 : s=0.;
439 : Bool_t filled=kTRUE;
440 412 : for(Int_t is(AliPID::kSPECIES); is--;){
441 : //if((PIDmethod==kLQ2D)&&(!(is==0||is==2)))continue;
442 618 : if((dE[0] > 0.) && (p[il] > 0.)) prLayer[is] = GetProbabilitySingleLayer(is, p[il], &dE[0],PIDmethod);
443 618 : AliDebug(3, Form("Probability for Species %d in Layer %d: %e", is, il, prLayer[is]));
444 206 : if(prLayer[is]<1.e-30){
445 618 : AliDebug(2, Form("Null for species %d species prob layer[%d].",is,il));
446 : filled=kFALSE;
447 206 : break;
448 : }
449 0 : s+=prLayer[is];
450 : }
451 206 : if(!filled){
452 206 : continue;
453 : }
454 0 : if(s<1.e-30){
455 0 : AliDebug(2, Form("Null all species prob layer[%d].", il));
456 0 : continue;
457 : }
458 0 : for(Int_t is(AliPID::kSPECIES); is--;){
459 0 : if(kNorm) prLayer[is] /= s; // probability in each layer for each particle species normalized to the sum of probabilities for given layer
460 0 : prob[is] *= prLayer[is]; // multiply single layer probabilities to get probability for complete track
461 : }
462 0 : ntrackletsPID++;
463 0 : }
464 104 : if(!kNorm) return ntrackletsPID;
465 :
466 : s=0.;
467 : // sum probabilities for all considered particle species
468 1248 : for(Int_t is(AliPID::kSPECIES); is--;) { s+=prob[is];
469 : }
470 104 : if(s<1.e-30){
471 0 : AliDebug(2, "Null total prob.");
472 0 : return 0;
473 : }
474 : // norm to the summed probability (default values: s=1 prob[is]=0.2)
475 1248 : for(Int_t is(AliPID::kSPECIES); is--;){ prob[is]/=s; }
476 104 : return ntrackletsPID;
477 208 : }
478 :
479 : //____________________________________________________________
480 : Double_t AliTRDPIDResponse::GetProbabilitySingleLayer(Int_t species, Double_t plocal, Double_t *dEdx,ETRDPIDMethod PIDmethod) const {
481 : //
482 : // Get the non-normalized probability for a certain particle species as coming
483 : // from the reference histogram
484 : // Interpolation between momentum bins
485 : //
486 824 : AliDebug(1, Form("Make Probability for Species %d with Momentum %f", species, plocal));
487 :
488 : Double_t probLayer = 0.;
489 :
490 206 : Float_t pLower, pUpper;
491 :
492 618 : AliTRDNDFast *refUpper = dynamic_cast<AliTRDNDFast *>(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,PIDmethod)),
493 488 : *refLower = dynamic_cast<AliTRDNDFast *>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,PIDmethod));
494 :
495 : // Do Interpolation exept for underflow and overflow
496 206 : if(refLower && refUpper){
497 0 : Double_t probLower = refLower->Eval(dEdx);
498 0 : Double_t probUpper = refUpper->Eval(dEdx);
499 :
500 0 : probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
501 206 : } else if(refLower){
502 : // underflow
503 0 : probLayer = refLower->Eval(dEdx);
504 206 : } else if(refUpper){
505 : // overflow
506 0 : probLayer = refUpper->Eval(dEdx);
507 0 : } else {
508 618 : AliDebug(3,"No references available");
509 : }
510 :
511 412 : switch(PIDmethod){
512 : case kLQ2D: // 2D LQ
513 : {
514 0 : AliDebug(1,Form("Eval 2D Q0 %f Q1 %f P %e ",dEdx[0],dEdx[1],probLayer));
515 : }
516 : break;
517 : case kLQ1D: // 1D LQ
518 : {
519 618 : AliDebug(1, Form("Eval 1D dEdx %f Probability %e", dEdx[0],probLayer));
520 : }
521 : break;
522 : case kLQ3D: // 3D LQ
523 : {
524 0 : AliDebug(1, Form("Eval 1D dEdx %f %f %f Probability %e", dEdx[0],dEdx[1],dEdx[2],probLayer));
525 : }
526 : break;
527 : case kLQ7D: // 7D LQ
528 : {
529 0 : 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));
530 : }
531 : break;
532 : default:
533 : break;
534 : }
535 :
536 206 : return probLayer;
537 :
538 : /* old implementation
539 :
540 : switch(PIDmethod){
541 : case kNN: // NN
542 : break;
543 : case kLQ2D: // 2D LQ
544 : {
545 : if(species==0||species==2){ // references only for electrons and pions
546 : Double_t error = 0.;
547 : Double_t point[kNslicesLQ2D];
548 : for(Int_t idim=0;idim<kNslicesLQ2D;idim++){point[idim]=dEdx[idim];}
549 :
550 : AliTRDTKDInterpolator *refUpper = dynamic_cast<AliTRDTKDInterpolator *>(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,kLQ2D)),
551 : *refLower = dynamic_cast<AliTRDTKDInterpolator *>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,kLQ2D));
552 : // Do Interpolation exept for underflow and overflow
553 : if(refLower && refUpper){
554 : Double_t probLower=0,probUpper=0;
555 : refLower->Eval(point,probLower,error);
556 : refUpper->Eval(point,probUpper,error);
557 : probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
558 : } else if(refLower){
559 : // underflow
560 : refLower->Eval(point,probLayer,error);
561 : } else if(refUpper){
562 : // overflow
563 : refUpper->Eval(point,probLayer,error);
564 : } else {
565 : AliError("No references available");
566 : }
567 : AliDebug(2,Form("Eval 2D Q0 %f Q1 %f P %e Err %e",point[0],point[1],probLayer,error));
568 : }
569 : }
570 : break;
571 : case kLQ1D: // 1D LQ
572 : {
573 : TH1 *refUpper = dynamic_cast<TH1 *>(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,kLQ1D)),
574 : *refLower = dynamic_cast<TH1 *>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,kLQ1D));
575 : // Do Interpolation exept for underflow and overflow
576 : if(refLower && refUpper){
577 : Double_t probLower = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx[0]));
578 : Double_t probUpper = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx[0]));
579 :
580 : probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
581 : } else if(refLower){
582 : // underflow
583 : probLayer = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx[0]));
584 : } else if(refUpper){
585 : // overflow
586 : probLayer = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx[0]));
587 : } else {
588 : AliError("No references available");
589 : }
590 : AliDebug(1, Form("Eval 1D dEdx %f Probability %e", dEdx[0],probLayer));
591 : }
592 : break;
593 : default:
594 : break;
595 : }
596 : return probLayer;
597 : */
598 :
599 206 : }
600 :
601 : //____________________________________________________________
602 : void AliTRDPIDResponse::SetOwner(){
603 : //
604 : // Make Deep Copy of the Reference Histograms
605 : //
606 0 : if(!fkPIDResponseObject || IsOwner()) return;
607 0 : const AliTRDPIDResponseObject *tmp = fkPIDResponseObject;
608 0 : fkPIDResponseObject = dynamic_cast<const AliTRDPIDResponseObject *>(tmp->Clone());
609 0 : SetBit(kIsOwner, kTRUE);
610 0 : }
611 :
612 : //____________________________________________________________
613 : Bool_t AliTRDPIDResponse::CookdEdx(Int_t nSlice, const Double_t * const in, Double_t *out,ETRDPIDMethod PIDmethod) const
614 : {
615 : //
616 : // Recalculate dE/dx
617 : // 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
618 : //
619 :
620 1454 : switch(PIDmethod){
621 : case kNN: // NN
622 : break;
623 : case kLQ2D: // 2D LQ
624 0 : out[0]=0;
625 0 : out[1]=0;
626 0 : for(Int_t islice = 0; islice < nSlice; islice++){
627 : // if(in[islice]<=0){out[0]=0;out[1]=0;return kFALSE;} // Require that all slices are filled
628 :
629 0 : if(islice<kNsliceQ0LQ2D)out[0]+= in[islice];
630 0 : else out[1]+= in[islice];
631 : }
632 : // normalize signal to number of slices
633 0 : out[0]*=1./Double_t(kNsliceQ0LQ2D);
634 0 : out[1]*=1./Double_t(nSlice-kNsliceQ0LQ2D);
635 0 : if(out[0] <= 0) return kFALSE;
636 0 : if(out[1] <= 0) return kFALSE;
637 0 : AliDebug(3,Form("CookdEdx Q0 %f Q1 %f",out[0],out[1]));
638 : break;
639 : case kLQ1D: // 1D LQ
640 624 : out[0]= 0.;
641 2496 : for(Int_t islice = 0; islice < nSlice; islice++) {
642 : // if(in[islice] > 0) out[0] += in[islice] * fGainNormalisationFactor; // Protect against negative values for slices having no dE/dx information
643 830 : if(in[islice] > 0) out[0] += in[islice]; // no neg dE/dx values
644 : }
645 624 : out[0]*=1./Double_t(kNsliceQ0LQ1D);
646 1042 : if(out[0] <= 0) return kFALSE;
647 618 : AliDebug(3,Form("CookdEdx dEdx %f",out[0]));
648 : break;
649 : case kLQ3D: // 3D LQ
650 0 : out[0]=0;
651 0 : out[1]=0;
652 0 : out[2]=0;
653 0 : for(Int_t islice = 0; islice < nSlice; islice++){
654 : // if(in[islice]<=0){out[0]=0;out[1]=0;out[2]=0;return kFALSE;} // Require that all slices are filled
655 0 : if(islice<kNsliceQ0LQ3D)out[0]+= in[islice];
656 0 : out[1]=(in[3]+in[4]);
657 0 : out[2]=(in[5]+in[6]);
658 : }
659 : // normalize signal to number of slices
660 0 : out[0]*=1./Double_t(kNsliceQ0LQ3D);
661 0 : out[1]*=1./2.;
662 0 : out[2]*=1./2.;
663 0 : if(out[0] <= 0) return kFALSE;
664 0 : if(out[1] <= 0) return kFALSE;
665 0 : if(out[2] <= 0) return kFALSE;
666 0 : AliDebug(3,Form("CookdEdx Q0 %f Q1 %f Q2 %f",out[0],out[1],out[2]));
667 : break;
668 : case kLQ7D: // 7D LQ
669 0 : for(Int_t i=0;i<nSlice;i++) {out[i]=0;}
670 0 : for(Int_t islice = 0; islice < nSlice; islice++){
671 0 : if(in[islice]<=0){
672 0 : for(Int_t i=0;i<8;i++){
673 0 : out[i]=0;
674 : }
675 0 : return kFALSE;} // Require that all slices are filled
676 0 : out[islice]=in[islice];
677 : }
678 0 : for(Int_t i=0;i<nSlice;i++) {if(out[i]<=0) return kFALSE; }
679 0 : 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]));
680 : break;
681 :
682 : default:
683 0 : return kFALSE;
684 : }
685 206 : return kTRUE;
686 624 : }
687 :
688 : //____________________________________________________________
689 : Bool_t AliTRDPIDResponse::IdentifiedAsElectron(Int_t nTracklets, const Double_t *like, Double_t p, Double_t level,Double_t centrality,ETRDPIDMethod PIDmethod) const {
690 : //
691 : // Check whether particle is identified as electron assuming a certain electron efficiency level
692 : // Only electron and pion hypothesis is taken into account
693 : //
694 : // Inputs:
695 : // Number of tracklets
696 : // Likelihood values
697 : // Momentum
698 : // Electron efficiency level
699 : //
700 : // If the function fails when the params are not accessible, the function returns true
701 : //
702 0 : if(!fkPIDResponseObject){
703 0 : AliDebug(3,"No PID Param object available");
704 0 : return kTRUE;
705 : }
706 0 : Double_t probEle = like[AliPID::kElectron]/(like[AliPID::kElectron] + like[AliPID::kPion]);
707 0 : AliDebug(3,Form("probabilities like %f %f %f \n",probEle,like[AliPID::kElectron],like[AliPID::kPion]));
708 0 : Double_t params[4];
709 0 : if(!fkPIDResponseObject->GetThresholdParameters(nTracklets, level, params,centrality,PIDmethod)){
710 0 : AliError("No Params found for the given configuration");
711 0 : return kTRUE;
712 : }
713 :
714 :
715 :
716 0 : Double_t threshold = 1. - params[0] - params[1] * p - params[2] * TMath::Exp(-params[3] * p);
717 0 : 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)));
718 0 : 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
719 0 : return kFALSE;
720 0 : }
721 :
722 : //____________________________________________________________
723 : Bool_t AliTRDPIDResponse::SetPIDResponseObject(const AliTRDPIDResponseObject * obj){
724 :
725 4 : fkPIDResponseObject = obj;
726 2 : if((AliLog::GetDebugLevel("",IsA()->GetName()))>0 && obj) fkPIDResponseObject->Print("");
727 2 : return kTRUE;
728 : }
729 :
730 :
731 : //____________________________________________________________
732 : Bool_t AliTRDPIDResponse::SetdEdxParams(const AliTRDdEdxParams * par)
733 : {
734 0 : fkTRDdEdxParams = par;
735 0 : return kTRUE;
736 : }
|
