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