PWGDQ/dielectron/AliDielectronTRDpidCut.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 // Class for TRD PID
  17 // Class further contains TRD specific cuts and QA histograms
  18 //
  19 // Authors:
  20 //   Jens Wiechula <Jens.Wiechula@cern.ch>
  21 //   Markus Fasel <M.Fasel@gsi.de>
  22 //
  23 #include <TList.h>
  24 #include <TString.h>
  25
  26 #include <AliVParticle.h>
  27 #include <AliESDtrack.h>
  28 #include <AliLog.h>
  29 #include <AliPID.h>
  30
  31 #include "AliDielectronTRDpidCut.h"
  32
  33 ClassImp(AliDielectronTRDpidCut)
  34
  35 const Double_t AliDielectronTRDpidCut::fgkVerySmall = 1e-12;
  36
  37 //___________________________________________________________________
  38 AliDielectronTRDpidCut::AliDielectronTRDpidCut() :
  39     AliAnalysisCuts()
  40   , fMinP(1.)
  41   , fElectronEfficiency(0.91)
  42   , fPIDMethod(kNN)
  43   , fRequirePIDbit(0)
  44 {
  45   //
  46   // default  constructor
  47   //
  48   memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
  49   SetUseDefaultParameters();
  50 }
  51
  52 //___________________________________________________________________
  53 AliDielectronTRDpidCut::AliDielectronTRDpidCut(const char* name) :
  54     AliAnalysisCuts(name,name)
  55   , fMinP(1.)
  56   , fElectronEfficiency(0.91)
  57   , fPIDMethod(kNN)
  58   , fRequirePIDbit(0)
  59 {
  60   //
  61   // default  constructor
  62   //
  63   memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
  64   SetUseDefaultParameters();
  65 }
  66
  67 //___________________________________________________________________
  68 AliDielectronTRDpidCut::AliDielectronTRDpidCut(const AliDielectronTRDpidCut &ref):
  69     AliAnalysisCuts(ref.GetName(),ref.GetTitle())
  70   , fMinP(ref.fMinP)
  71   , fElectronEfficiency(ref.fElectronEfficiency)
  72   , fPIDMethod(ref.fPIDMethod)
  73   , fRequirePIDbit(ref.fRequirePIDbit)
  74 {
  75   //
  76   // Copy constructor
  77   //
  78   memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
  79   ref.Copy(*this);
  80 }
  81
  82 //___________________________________________________________________
  83 AliDielectronTRDpidCut &AliDielectronTRDpidCut::operator=(const AliDielectronTRDpidCut &ref){
  84   //
  85   // Assignment operator
  86   //
  87   if(this != &ref){
  88     ref.Copy(*this);
  89   }
  90   return *this;
  91 }
  92
  93 //___________________________________________________________________
  94 void AliDielectronTRDpidCut::Copy(TObject &ref) const {
  95   //
  96   // Performs the copying of the object
  97   //
  98   AliDielectronTRDpidCut &target = dynamic_cast<AliDielectronTRDpidCut &>(ref);
  99
 100   Bool_t defaultParameters = UseDefaultParameters();
 101   target.SetUseDefaultParameters(defaultParameters);
 102   target.fMinP = fMinP;
 103   target.fPIDMethod = fPIDMethod;
 104   target.fRequirePIDbit = fRequirePIDbit;
 105   target.fElectronEfficiency = fElectronEfficiency;
 106   memcpy(target.fThreshParams, fThreshParams, sizeof(Double_t) * kThreshParams);
 107   AliAnalysisCuts::Copy(ref);
 108 }
 109
 110 //___________________________________________________________________
 111 AliDielectronTRDpidCut::~AliDielectronTRDpidCut(){
 112   //
 113   // Destructor
 114   //
 115 }
 116
 117 //______________________________________________________
 118 Bool_t AliDielectronTRDpidCut::InitializePID(){
 119   //
 120   // InitializePID: Load TRD thresholds and create the electron efficiency axis
 121   // to navigate
 122   //
 123   if(UseDefaultParameters()){
 124     if(fPIDMethod == kLQ)
 125       InitParameters1DLQ();
 126     else
 127       InitParameters();
 128   }
 129   return kTRUE;
 130 }
 131
 132 //______________________________________________________
 133 Bool_t AliDielectronTRDpidCut::IsSelected(TObject *track) {
 134   //
 135   // Does PID for TRD alone:
 136   // PID thresholds based on 90% Electron Efficiency level approximated by a linear
 137   // step function
 138   //
 139   const AliESDtrack *part = dynamic_cast<const AliESDtrack *>(track);
 140   if (!part) return kFALSE;
 141
 142   if (fRequirePIDbit==AliDielectronTRDpidCut::kRequire&&!(part->GetStatus()&AliESDtrack::kTRDpid)) return kFALSE;
 143   if (fRequirePIDbit==AliDielectronTRDpidCut::kIfAvailable&&!(part->GetStatus()&AliESDtrack::kTRDpid)) return kTRUE;
 144
 145   Double_t p = GetP((AliVParticle*)track);
 146   if(p < fMinP){
 147     AliDebug(1, Form("Track momentum below %f", fMinP));
 148     if(fRequirePIDbit==AliDielectronTRDpidCut::kRequire) return kFALSE;
 149         if(fRequirePIDbit==AliDielectronTRDpidCut::kIfAvailable) return kTRUE;
 150   }
 151
 152   Double_t electronLike = GetElectronLikelihood((AliVParticle*)track);
 153   Double_t threshold = GetTRDthresholds(fElectronEfficiency, p);
 154   AliDebug(1, Form("Threshold: %f\n", threshold));
 155   if(electronLike > threshold){
 156     return kTRUE;
 157   }
 158   return kFALSE;
 159
 160 }
 161
 162 //___________________________________________________________________
 163 Double_t AliDielectronTRDpidCut::GetTRDthresholds(Double_t electronEff, Double_t p) const {
 164   //
 165   // Return momentum dependent and electron efficiency dependent TRD thresholds
 166   //
 167   Double_t params[4];
 168   GetParameters(electronEff, params);
 169   Double_t threshold = 1. - params[0] - params[1] * p - params[2] * TMath::Exp(-params[3] * p);
 170   return TMath::Max(TMath::Min(threshold, 0.99), 0.2); // truncate the threshold upperwards to 0.999 and lowerwards to 0.2 and exclude unphysical values
 171 }
 172
 173 //___________________________________________________________________
 174 void AliDielectronTRDpidCut::SetThresholdParameters(Double_t electronEff, Double_t *params){
 175   //
 176   // Set threshold parameters for the given bin
 177   //
 178   if(electronEff >= 1. || electronEff < 0.7) return;
 179   Int_t effbin = static_cast<Int_t>((electronEff - 0.7)/0.05);
 180   memcpy(&fThreshParams[effbin * 4], params, sizeof(Double_t) * 4);
 181   SetUseDefaultParameters(kFALSE);
 182 }
 183
 184 //___________________________________________________________________
 185 void AliDielectronTRDpidCut::InitParameters(){
 186   //
 187   // Fill the Parameters into an array
 188   //
 189
 190   AliDebug(2, "Loading threshold Parameter");
 191   // Parameters for 6 Layers
 192   fThreshParams[0] = -0.001839; // 0.7 electron eff
 193   fThreshParams[1] = 0.000276;
 194   fThreshParams[2] = 0.044902;
 195   fThreshParams[3] = 1.726751;
 196   fThreshParams[4] = -0.002405; // 0.75 electron eff
 197   fThreshParams[5] = 0.000372;
 198   fThreshParams[6] = 0.061775;
 199   fThreshParams[7] = 1.739371;
 200   fThreshParams[8] = -0.003178; // 0.8 electron eff
 201   fThreshParams[9] = 0.000521;
 202   fThreshParams[10] = 0.087585;
 203   fThreshParams[11] = 1.749154;
 204   fThreshParams[12] = -0.004058; // 0.85 electron eff
 205   fThreshParams[13] = 0.000748;
 206   fThreshParams[14] = 0.129583;
 207   fThreshParams[15] = 1.782323;
 208   fThreshParams[16] = -0.004967; // 0.9 electron eff
 209   fThreshParams[17] = 0.001216;
 210   fThreshParams[18] = 0.210128;
 211   fThreshParams[19] = 1.807665;
 212   fThreshParams[20] = -0.000996; // 0.95 electron eff
 213   fThreshParams[21] = 0.002627;
 214   fThreshParams[22] = 0.409099;
 215   fThreshParams[23] = 1.787076;
 216 }
 217
 218 //___________________________________________________________________
 219 void AliDielectronTRDpidCut::InitParameters1DLQ(){
 220   //
 221   // Init Parameters for 1DLQ PID (M. Fasel, Sept. 6th, 2010)
 222   //
 223
 224   // Parameters for 6 Layers
 225   AliDebug(2, Form("Loading threshold parameter for Method 1DLQ"));
 226   fThreshParams[0] = -0.02241; // 0.7 electron eff
 227   fThreshParams[1] = 0.05043;
 228   fThreshParams[2] = 0.7925;
 229   fThreshParams[3] = 2.625;
 230   fThreshParams[4] = 0.07438; // 0.75 electron eff
 231   fThreshParams[5] = 0.05158;
 232   fThreshParams[6] = 2.864;
 233   fThreshParams[7] = 4.356;
 234   fThreshParams[8] = 0.1977; // 0.8 electron eff
 235   fThreshParams[9] = 0.05956;
 236   fThreshParams[10] = 2.853;
 237   fThreshParams[11] = 3.713;
 238   fThreshParams[12] = 0.5206; // 0.85 electron eff
 239   fThreshParams[13] = 0.03077;
 240   fThreshParams[14] = 2.966;
 241   fThreshParams[15] = 4.07;
 242   fThreshParams[16] = 0.8808; // 0.9 electron eff
 243   fThreshParams[17] = 0.002092;
 244   fThreshParams[18] = 1.17;
 245   fThreshParams[19] = 4.506;
 246   fThreshParams[20] = 1.; // 0.95 electron eff
 247   fThreshParams[21] = 0.;
 248   fThreshParams[22] = 0.;
 249   fThreshParams[23] = 0.;
 250
 251 }
 252
 253 //___________________________________________________________________
 254 void AliDielectronTRDpidCut::RenormalizeElPi(const Double_t *likein, Double_t *likeout) const {
 255   //
 256   // Renormalize likelihoods for electrons and pions neglecting the
 257   // likelihoods for protons, kaons and muons
 258   //
 259   memset(likeout, 0, sizeof(Double_t) * AliPID::kSPECIES);
 260   Double_t norm = likein[AliPID::kElectron] + likein[AliPID::kPion];
 261   if(norm == 0.) norm = 1.;   // Safety
 262   likeout[AliPID::kElectron] = likein[AliPID::kElectron] / norm;
 263   likeout[AliPID::kPion] = likein[AliPID::kPion] / norm;
 264 }
 265
 266 //___________________________________________________________________
 267 void AliDielectronTRDpidCut::GetParameters(Double_t electronEff, Double_t *parameters) const {
 268   //
 269   // return parameter set for the given efficiency bin
 270   //
 271   Int_t effbin = static_cast<Int_t>((electronEff - 0.7)/0.05);
 272   memcpy(parameters, fThreshParams + effbin * 4, sizeof(Double_t) * 4);
 273 }
 274
 275 //___________________________________________________________________
 276 Double_t AliDielectronTRDpidCut::GetElectronLikelihood(const AliVParticle *track) const {
 277   //
 278   // Get TRD likelihoods for ESD respectively AOD tracks
 279   //
 280   Double_t pidProbs[AliPID::kSPECIES]; memset(pidProbs, 0, sizeof(Double_t) * AliPID::kSPECIES);
 281   const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
 282   if(esdtrack) esdtrack->GetTRDpid(pidProbs);
 283   if(!IsRenormalizeElPi()) return pidProbs[AliPID::kElectron];
 284   Double_t probsNew[AliPID::kSPECIES];
 285   RenormalizeElPi(pidProbs, probsNew);
 286   return probsNew[AliPID::kElectron];
 287 }
 288
 289 //___________________________________________________________________
 290 Double_t AliDielectronTRDpidCut::GetP(const AliVParticle *track) const {
 291   //
 292   // Get the Momentum in the TRD
 293   //
 294   Double_t p = 0.;
 295   const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
 296   if(esdtrack) p = esdtrack->GetOuterParam() ? esdtrack->GetOuterParam()->P() : esdtrack->P();
 297   return p;
 298 }
 299
 300 //___________________________________________________________________
 301 Double_t AliDielectronTRDpidCut::GetChargeLayer(const AliVParticle *track, UInt_t layer) const {
 302   //
 303   // Get the Charge in a single TRD layer
 304   //
 305   if(layer >= 6) return 0.;
 306   Double_t charge = 0.;
 307   const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
 308   if(esdtrack)
 309     for(Int_t islice = 0; islice < esdtrack->GetNumberOfTRDslices(); islice++)
 310       charge += esdtrack->GetTRDslice(static_cast<UInt_t>(layer), islice);
 311   return charge;
 312 }
 313
 314 //___________________________________________________________________
 315 void AliDielectronTRDpidCut::GetTRDmomenta(const AliVParticle *track, Double_t *mom) const {
 316   //
 317   // Fill Array with momentum information at the TRD tracklet
 318   //
 319   const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
 320   if(esdtrack)
 321     for(Int_t itl = 0; itl < 6; itl++)
 322       mom[itl] = esdtrack->GetTRDmomentum(itl);
 323 }
 324
 325 //___________________________________________________________________
 326 Double_t AliDielectronTRDpidCut::GetTRDSignalV1(const AliESDtrack *track, Float_t truncation) const {
 327   //
 328   // Calculation of the TRD Signal via truncated mean
 329   // Method 1: Take all Slices available
 330   // cut out 0s
 331   // Order them in increasing order
 332   // Cut out the upper third
 333   // Calculate mean over the last 2/3 slices
 334   //
 335   const Int_t kNSlices = 48;
 336   const Int_t kSlicePerLayer = 7;
 337   // Weight the slice to equalize the MPV of the dQ/dl-distribution per slice to the one in the first slice
 338   // Pions are used as reference for the equalization
 339   const Double_t kWeightSlice[8] = {1., 2.122, 1.8, 1.635, 1.595, 1.614, 1.16, 7.0};
 340   AliDebug(3, Form("Number of Tracklets: %d\n", track->GetTRDntrackletsPID()));
 341   Double_t trdSlices[kNSlices], tmp[kNSlices];
 342   Int_t indices[48];
 343   Int_t icnt = 0;
 344   for(Int_t idet = 0; idet < 6; idet++)
 345     for(Int_t islice = 0; islice < kSlicePerLayer; islice++){
 346       AliDebug(2, Form("Chamber[%d], Slice[%d]: TRDSlice = %f", idet, islice, track->GetTRDslice(idet, islice)));
 347       if(TMath::Abs(track->GetTRDslice(idet, islice)) < fgkVerySmall) continue;;
 348       trdSlices[icnt++] = track->GetTRDslice(idet, islice) * kWeightSlice[islice];
 349     }
 350   AliDebug(1, Form("Number of Slices: %d\n", icnt));
 351   if(icnt < 6) return 0.;   // We need at least 6 Slices for the truncated mean
 352   TMath::Sort(icnt, trdSlices, indices, kFALSE);
 353   memcpy(tmp, trdSlices, sizeof(Double_t) * icnt);
 354   for(Int_t ien = 0; ien < icnt; ien++)
 355     trdSlices[ien] = tmp[indices[ien]];
 356   Double_t trdSignal = TMath::Mean(static_cast<Int_t>(static_cast<Float_t>(icnt) * truncation), trdSlices);
 357   Double_t mom = track->GetOuterParam() ? track->GetOuterParam()->P() : -1;
 358   AliDebug(3, Form("PID Meth. 1: p[%f], TRDSignal[%f]", mom, trdSignal));
 359   return trdSignal;
 360 }
 361
 362 //___________________________________________________________________
 363 Double_t AliDielectronTRDpidCut::GetTRDSignalV2(const AliESDtrack *track, Float_t truncation) const {
 364   //
 365   // Calculation of the TRD Signal via truncated mean
 366   // Method 2: Take only first 5 slices per chamber
 367   // Order them in increasing order
 368   // Cut out upper half
 369   // Now do mean with the reamining 3 slices per chamber
 370   //
 371   const Double_t kWeightSlice[8] = {1., 2.122, 1.8, 1.635, 1.595, 1.614, 1.16, 7.0};
 372   const Int_t kLayers = 6;
 373   const Int_t kSlicesLow = 6;
 374   const Int_t kSlicesHigh = 1;
 375   Double_t trdSlicesLowTime[kLayers*kSlicesLow], trdSlicesRemaining[kLayers*(kSlicesHigh + kSlicesLow)];
 376   Int_t indices[kLayers*kSlicesLow];
 377   Int_t cntLowTime=0, cntRemaining = 0;
 378   for(Int_t idet = 0; idet < 6; idet++)
 379     for(Int_t islice = 0; islice < kSlicesLow+kSlicesHigh; islice++){
 380       if(TMath::Abs(track->GetTRDslice(idet, islice)) < fgkVerySmall) continue;;
 381       if(islice < kSlicesLow){
 382         AliDebug(3, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice)));
 383         trdSlicesLowTime[cntLowTime++] = track->GetTRDslice(idet, islice) * kWeightSlice[islice];
 384       } else{
 385         AliDebug(3, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice)));
 386         trdSlicesRemaining[cntRemaining++] = track->GetTRDslice(idet, islice) * kWeightSlice[islice];
 387       }
 388     }
 389   if(cntLowTime < 4 || cntRemaining < 2) return 0.; // Min. Number of Slices at high time is 2 (matches with 1 layer), for the truncated mean we need at least 4 Slices
 390   TMath::Sort(cntLowTime, trdSlicesLowTime, indices, kFALSE);
 391   // Fill the second array with the lower half of the first time bins
 392   for(Int_t ien = 0; ien < static_cast<Int_t>(static_cast<Float_t>(cntLowTime) * truncation); ien++)
 393     trdSlicesRemaining[cntRemaining++] = trdSlicesLowTime[indices[ien]];
 394   Double_t trdSignal = TMath::Mean(cntRemaining, trdSlicesRemaining);
 395   Double_t mom = track->GetOuterParam() ? track->GetOuterParam()->P() : -1;
 396   AliDebug(3, Form("PID Meth. 2: p[%f], TRDSignal[%f]", mom, trdSignal));
 397   return trdSignal;
 398 }