PWG2/RESONANCES/AliRsnCutPID.cxx

   1 //
   2 // Class AliRsnCutPID
   3 //
   4 // General implementation of a single cut strategy, which can be:
   5 // - a value contained in a given interval  [--> IsBetween()   ]
   6 // - a value equal to a given reference     [--> MatchesValue()]
   7 //
   8 // In all cases, the reference value(s) is (are) given as data members
   9 // and each kind of cut requires a given value type (Int, UInt, Double),
  10 // but the cut check procedure is then automatized and chosen thanks to
  11 // an enumeration of the implemented cut types.
  12 // At the end, the user (or any other point which uses this object) has
  13 // to use the method IsSelected() to check if this cut has been passed.
  14 //
  15 // authors: Martin Vala (martin.vala@cern.ch)
  16 //          Alberto Pulvirenti (alberto.pulvirenti@ct.infn.it)
  17 //
  18 #include "TMath.h"
  19
  20 #include "AliExternalTrackParam.h"
  21
  22 #include "AliRsnDaughter.h"
  23 #include "AliRsnCutPID.h"
  24
  25 ClassImp(AliRsnCutPID)
  26
  27 //_________________________________________________________________________________________________
  28 AliRsnCutPID::AliRsnCutPID() :
  29   AliRsnCut(),
  30   fPerfect(kFALSE),
  31   fUseDefault(kTRUE)
  32 {
  33 //
  34 // Default constructor.
  35 // Sets the cut to realistic PID with default weights,
  36 // and defines the 'fMinI' value of the base class as the PID
  37 // to which we want to compare this object.
  38 //
  39
  40   Int_t i;
  41
  42   for (i = 0; i < kDetectors; i++)
  43   {
  44     fUseDetector[i] = kFALSE;
  45     fPtThreshold[i] = 0.0;
  46     fGoAboveThreshold[i] = kTRUE;
  47   }
  48
  49   for (i = 0; i < AliPID::kSPECIES; i++)
  50   {
  51     fWeight[i] = 0.0;
  52     fPrior[i] = 1.0;
  53   }
  54
  55   SetTargetType(AliRsnTarget::kDaughter);
  56 }
  57
  58 //_________________________________________________________________________________________________
  59 AliRsnCutPID::AliRsnCutPID(const char *name, AliPID::EParticleType pid, Double_t probMin, Bool_t perfectPID) :
  60   AliRsnCut(name, AliRsnCut::kDaughter, (Int_t)pid),
  61   fPerfect(perfectPID),
  62   fUseDefault(kTRUE)
  63 {
  64 //
  65 // Default constructor.
  66 // Sets the cut to realistic PID with default weights,
  67 // and defines the 'fMinI' value of the base class as the PID
  68 // to which we want to compare this object.
  69 //
  70
  71   Int_t i;
  72
  73   for (i = 0; i < kDetectors; i++)
  74   {
  75     fUseDetector[i] = kFALSE;
  76     fPtThreshold[i] = 0.0;
  77     fGoAboveThreshold[i] = kTRUE;
  78   }
  79
  80   for (i = 0; i < AliPID::kSPECIES; i++)
  81   {
  82     fWeight[i] = 0.0;
  83     fPrior[i] = 1.0;
  84   }
  85
  86   fMinD = probMin;
  87   fMaxD = 1.000001;
  88
  89   SetTargetType(AliRsnTarget::kDaughter);
  90 }
  91
  92 //_____________________________________________________________________________
  93 Bool_t AliRsnCutPID::CheckThreshold(EDetector det, Double_t value)
  94 {
  95 //
  96 // Checks if the passed value (track pT) stays in the
  97 // interval where the detector should be accepted
  98 //
  99
 100   if (!CheckBounds(det)) return kFALSE;
 101
 102   if (fGoAboveThreshold[det]) return (value >= fPtThreshold[det]);
 103   else return (value <= fPtThreshold[det]);
 104 }
 105
 106 //_________________________________________________________________________________________________
 107 Bool_t AliRsnCutPID::ComputeWeights(AliRsnDaughter *daughter)
 108 {
 109 //
 110 // Compute the PID weights using the class settings.
 111 // If the argument is an ESD track, customized weights can be computed
 112 // It the argument is a track (ESD or AOD), at least default weights
 113 // can be computed, otherwise, no weights can be defined.
 114 // The return value tells if the operation was successful.
 115 //
 116
 117   Int_t  i, j;
 118   Bool_t useDefault = fUseDefault;
 119   Bool_t perfectPID = fPerfect;
 120   if (perfectPID && !daughter->GetRefMC()) return kFALSE;
 121   if (!daughter->GetRefESDtrack()) useDefault = kTRUE;
 122   if (!daughter->GetRefESDtrack() && !daughter->GetRefAODtrack()) return kFALSE;
 123
 124   // if perfect PID ise required,
 125   // compare the PDG code of the type stored in 'fMinI' of the cut
 126   // and that of the particle which is checked, looking at its MC
 127   if (perfectPID)
 128   {
 129     i = TMath::Abs(AliPID::ParticleCode(fMinI));
 130     j = daughter->GetPDG();
 131     return (i == j);
 132   }
 133
 134   // if default weights are (or need to be) used,
 135   // they are taken directly and function exits
 136   if (useDefault)
 137   {
 138     if (daughter->GetRefESDtrack())
 139       daughter->GetRefESDtrack()->GetESDpid(fWeight);
 140     else
 141     {
 142       for (i = 0; i < AliPID::kSPECIES; i++)
 143         fWeight[i] = daughter->GetRefAODtrack()->PID()[i];
 144     }
 145     return kTRUE;
 146   }
 147
 148   // if we arrive here, this means that we have an ESD track
 149   // and we want to customize the PID
 150   AliESDtrack *track = daughter->GetRefESDtrack();
 151   Double_t     w[kDetectors][AliPID::kSPECIES];
 152   track->GetITSpid(w[kITS]);
 153   track->GetTPCpid(w[kTPC]);
 154   track->GetTRDpid(w[kTRD]);
 155   track->GetTOFpid(w[kTOF]);
 156   track->GetHMPIDpid(w[kHMPID]);
 157
 158   // if a detector is excluded or the track has not the right pT
 159   // all related weights are set to 1 in order not to contribute
 160   // to final multiplication
 161   for (i = 0; i < kDetectors; i++)
 162   {
 163     if (!fUseDetector[i] || !CheckThreshold((EDetector)i, track->Pt()))
 164     {
 165       for (j = 0; j < AliPID::kSPECIES; j++) {
 166         w[i][j] = 1.0;
 167       }
 168     }
 169   }
 170
 171   // multiplicate all weights to compute final one
 172   for (i = 0; i < AliPID::kSPECIES; i++)
 173   {
 174     fWeight[i] = w[kITS][i] * w[kTPC][i] * w[kTRD][i] * w[kTOF][i] * w[kHMPID][i];
 175   }
 176
 177   return kTRUE;
 178 }
 179
 180 //_________________________________________________________________________________________________
 181 Int_t AliRsnCutPID::RealisticPID(AliRsnDaughter * const daughter, Double_t &prob)
 182 {
 183 //
 184 // Combines the weights collected from chosen detectors with the priors
 185 // and gives the corresponding particle with the largest probability,
 186 // in terms of the AliPID particle type enumeration.
 187 // The argument, passed by reference, gives the corresponding probability,
 188 // since the cut could require that it is larger than a threshold.
 189 //
 190
 191   // try to compute the weights
 192   if (!ComputeWeights(daughter))
 193   {
 194     prob = -1.0;
 195     return AliPID::kUnknown;
 196   }
 197
 198   // combine with priors and normalize
 199   Int_t    i;
 200   Double_t sum = 0.0, w[AliPID::kSPECIES];
 201   for (i = 0; i < AliPID::kSPECIES; i++)
 202   {
 203     w[i] = fWeight[i] * fPrior[i];
 204     sum += w[i];
 205   }
 206   if (sum <= 0.0)
 207   {
 208     AliError("Sum = 0");
 209     prob = -1.0;
 210     return AliPID::kUnknown;
 211   }
 212   for (i = 0; i < AliPID::kSPECIES; i++) w[i] /= sum;
 213
 214   // find the largest probability and related PID
 215   Int_t ibest = 0;
 216   prob = w[0];
 217   for (i = 1; i < AliPID::kSPECIES; i++)
 218   {
 219     if (w[i] > prob)
 220     {
 221       prob = w[i];
 222       ibest = i;
 223     }
 224   }
 225
 226   // return the value, while the probability
 227   // will be consequentially set
 228   return ibest;
 229 }
 230
 231 //_________________________________________________________________________________________________
 232 Int_t AliRsnCutPID::PerfectPID(AliRsnDaughter * const daughter)
 233 {
 234 //
 235 // If MC is present, retrieve the particle corresponding to the used track
 236 // (using the fRefMC data member) and then return the true particle type
 237 // taken from the PDG code of the reference particle itself, converted
 238 // into the enumeration from AliPID object.
 239 //
 240
 241   // works only if the MC is present
 242   if (!daughter->GetRefMC()) return AliPID::kUnknown;
 243
 244   // get the PDG code of the particle
 245   Int_t pdg = daughter->GetPDG();
 246
 247   // loop over all species listed in AliPID to find the match
 248   Int_t i;
 249   for (i = 0; i < AliPID::kSPECIES; i++)
 250   {
 251     if (AliPID::ParticleCode(i) == pdg) return i;
 252   }
 253
 254   return AliPID::kUnknown;
 255 }
 256
 257 //_________________________________________________________________________________________________
 258 Bool_t AliRsnCutPID::IsSelected(TObject *object)
 259 {
 260 //
 261 // Cut checker.
 262 //
 263
 264   // convert the object into the unique correct type
 265
 266   if (!TargetOK(object))
 267   {
 268     AliError(Form("[%s]: this cut works only with AliRsnDaughter objects", GetName()));
 269     return kTRUE;
 270   }
 271
 272   // if target is OK, do a dynamic cast
 273   AliRsnDaughter *daughter = dynamic_cast<AliRsnDaughter*>(object);
 274
 275   // depending on the PID type, recalls the appropriate method:
 276   // in case of perfect PID, checks only if the PID type is
 277   // corresponding to the request in the cut (fMinI)
 278   // while in case of realistic PID checks also the probability
 279   // to be within the required interval
 280   if (fPerfect)
 281   {
 282     fCutValueI = PerfectPID(daughter);
 283     return OkValueI();
 284   }
 285   else
 286   {
 287     fCutValueI = RealisticPID(daughter, fCutValueD);
 288     return OkValueI() && OkRangeD();
 289   }
 290 }
 291
 292 //__________________________________________________________________________________________________
 293 void AliRsnCutPID::IncludeDetector(EDetector det, Double_t threshold, Bool_t goAbove)
 294 {
 295 //
 296 // Include a detector for a customized weight computing
 297 // and specify also its eventual threshold and if the detector
 298 // must be used above or below the threshold.
 299 // By default the threshold is zero and detector is always used.
 300 //
 301
 302   if (!CheckBounds(det)) return;
 303
 304   fUseDetector[det] = kTRUE;
 305   fPtThreshold[det] = threshold;
 306   fGoAboveThreshold[det] = goAbove;
 307 }