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       fUseDetector[i] = kFALSE;
  44       fPtThreshold[i] = 0.0;
  45       fGoAboveThreshold[i] = kTRUE;
  46    }
  47
  48    for (i = 0; i < AliPID::kSPECIES; i++) {
  49       fWeight[i] = 0.0;
  50       fPrior[i] = 1.0;
  51    }
  52
  53    SetTargetType(AliRsnTarget::kDaughter);
  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       fUseDetector[i] = kFALSE;
  73       fPtThreshold[i] = 0.0;
  74       fGoAboveThreshold[i] = kTRUE;
  75    }
  76
  77    for (i = 0; i < AliPID::kSPECIES; i++) {
  78       fWeight[i] = 0.0;
  79       fPrior[i] = 1.0;
  80    }
  81
  82    fMinD = probMin;
  83    fMaxD = 1.000001;
  84
  85    SetTargetType(AliRsnTarget::kDaughter);
  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       i = TMath::Abs(AliPID::ParticleCode(fMinI));
 125       j = daughter->GetPDG();
 126       return (i == j);
 127    }
 128
 129    // if default weights are (or need to be) used,
 130    // they are taken directly and function exits
 131    if (useDefault) {
 132       if (daughter->GetRefESDtrack())
 133          daughter->GetRefESDtrack()->GetESDpid(fWeight);
 134       else {
 135          for (i = 0; i < AliPID::kSPECIES; i++)
 136             fWeight[i] = daughter->GetRefAODtrack()->PID()[i];
 137       }
 138       return kTRUE;
 139    }
 140
 141    // if we arrive here, this means that we have an ESD track
 142    // and we want to customize the PID
 143    AliESDtrack *track = daughter->GetRefESDtrack();
 144    if (!track) return kFALSE;
 145    Double_t     w[kDetectors][AliPID::kSPECIES];
 146    track->GetITSpid(w[kITS]);
 147    track->GetTPCpid(w[kTPC]);
 148    track->GetTRDpid(w[kTRD]);
 149    track->GetTOFpid(w[kTOF]);
 150    track->GetHMPIDpid(w[kHMPID]);
 151
 152    // if a detector is excluded or the track has not the right pT
 153    // all related weights are set to 1 in order not to contribute
 154    // to final multiplication
 155    for (i = 0; i < kDetectors; i++) {
 156       if (!fUseDetector[i] || !CheckThreshold((EDetector)i, track->Pt())) {
 157          for (j = 0; j < AliPID::kSPECIES; j++) {
 158             w[i][j] = 1.0;
 159          }
 160       }
 161    }
 162
 163    // multiplicate all weights to compute final one
 164    for (i = 0; i < AliPID::kSPECIES; i++) {
 165       fWeight[i] = w[kITS][i] * w[kTPC][i] * w[kTRD][i] * w[kTOF][i] * w[kHMPID][i];
 166    }
 167
 168    return kTRUE;
 169 }
 170
 171 //_________________________________________________________________________________________________
 172 Int_t AliRsnCutPID::RealisticPID(AliRsnDaughter * const daughter, Double_t &prob)
 173 {
 174 //
 175 // Combines the weights collected from chosen detectors with the priors
 176 // and gives the corresponding particle with the largest probability,
 177 // in terms of the AliPID particle type enumeration.
 178 // The argument, passed by reference, gives the corresponding probability,
 179 // since the cut could require that it is larger than a threshold.
 180 //
 181
 182    // try to compute the weights
 183    if (!ComputeWeights(daughter)) {
 184       prob = -1.0;
 185       return AliPID::kUnknown;
 186    }
 187
 188    // combine with priors and normalize
 189    Int_t    i;
 190    Double_t sum = 0.0, w[AliPID::kSPECIES];
 191    for (i = 0; i < AliPID::kSPECIES; i++) {
 192       w[i] = fWeight[i] * fPrior[i];
 193       sum += w[i];
 194    }
 195    if (sum <= 0.0) {
 196       AliError("Sum = 0");
 197       prob = -1.0;
 198       return AliPID::kUnknown;
 199    }
 200    for (i = 0; i < AliPID::kSPECIES; i++) w[i] /= sum;
 201
 202    // find the largest probability and related PID
 203    Int_t ibest = 0;
 204    prob = w[0];
 205    for (i = 1; i < AliPID::kSPECIES; i++) {
 206       if (w[i] > prob) {
 207          prob = w[i];
 208          ibest = i;
 209       }
 210    }
 211
 212    // return the value, while the probability
 213    // will be consequentially set
 214    return ibest;
 215 }
 216
 217 //_________________________________________________________________________________________________
 218 Int_t AliRsnCutPID::PerfectPID(AliRsnDaughter * const daughter)
 219 {
 220 //
 221 // If MC is present, retrieve the particle corresponding to the used track
 222 // (using the fRefMC data member) and then return the true particle type
 223 // taken from the PDG code of the reference particle itself, converted
 224 // into the enumeration from AliPID object.
 225 //
 226
 227    // works only if the MC is present
 228    if (!daughter->GetRefMC()) return AliPID::kUnknown;
 229
 230    // get the PDG code of the particle
 231    Int_t pdg = daughter->GetPDG();
 232
 233    // loop over all species listed in AliPID to find the match
 234    Int_t i;
 235    for (i = 0; i < AliPID::kSPECIES; i++) {
 236       if (AliPID::ParticleCode(i) == pdg) return i;
 237    }
 238
 239    return AliPID::kUnknown;
 240 }
 241
 242 //_________________________________________________________________________________________________
 243 Bool_t AliRsnCutPID::IsSelected(TObject *object)
 244 {
 245 //
 246 // Cut checker.
 247 //
 248
 249    // convert the object into the unique correct type
 250
 251    if (!TargetOK(object)) {
 252       AliError(Form("[%s]: this cut works only with AliRsnDaughter objects", GetName()));
 253       return kTRUE;
 254    }
 255
 256    // if target is OK, do a dynamic cast
 257    AliRsnDaughter *daughter = fDaughter;
 258
 259    // depending on the PID type, recalls the appropriate method:
 260    // in case of perfect PID, checks only if the PID type is
 261    // corresponding to the request in the cut (fMinI)
 262    // while in case of realistic PID checks also the probability
 263    // to be within the required interval
 264    if (fPerfect && daughter) {
 265       fCutValueI = PerfectPID(daughter);
 266       return OkValueI();
 267    } else if (daughter) {
 268       fCutValueI = RealisticPID(daughter, fCutValueD);
 269       return OkValueI() && OkRangeD();
 270    } else
 271       return kFALSE;
 272 }
 273
 274 //__________________________________________________________________________________________________
 275 void AliRsnCutPID::IncludeDetector(EDetector det, Double_t threshold, Bool_t goAbove)
 276 {
 277 //
 278 // Include a detector for a customized weight computing
 279 // and specify also its eventual threshold and if the detector
 280 // must be used above or below the threshold.
 281 // By default the threshold is zero and detector is always used.
 282 //
 283
 284    if (!CheckBounds(det)) return;
 285
 286    fUseDetector[det] = kTRUE;
 287    fPtThreshold[det] = threshold;
 288    fGoAboveThreshold[det] = goAbove;
 289 }