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 }