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