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 Bool_t AliRsnCutPID::IsSelected(TObject *obj1, TObject* /*obj2*/)
 178 {
 179 //
 180 // Cut checker.
 181 //
 182
 183   // coherence check
 184   if (!AliRsnCut::TargetOK(obj1))
 185   {
 186     AliError(Form("Wrong target. Skipping cut", GetName()));
 187     return kTRUE;
 188   }
 189
 190   // convert the object into the unique correct type
 191   AliRsnDaughter *daughter = dynamic_cast<AliRsnDaughter*>(obj1);
 192
 193   // try to compute the weights
 194   if (!ComputeWeights(daughter)) return kFALSE;
 195
 196   // combine with priors and get the majority
 197   Int_t    i;
 198   Double_t sum = 0.0, w[AliPID::kSPECIES];
 199   for (i = 0; i < AliPID::kSPECIES; i++)
 200   {
 201     w[i] = fWeight[i] * fPrior[i];
 202     sum += w[i];
 203   }
 204   if (sum <= 0.0)
 205   {
 206     AliError("Sum = 0");
 207     return kFALSE;
 208   }
 209   for (i = 0; i < AliPID::kSPECIES; i++) w[i] /= sum;
 210
 211   // find the largest probability and related PID
 212   // and assign them to the mother class members which
 213   // are checked for the cut
 214   fCutValueI = 0;
 215   fCutValueD = w[0];
 216   for (i = 1; i < AliPID::kSPECIES; i++)
 217   {
 218     if (w[i] > fCutValueD)
 219     {
 220       fCutValueD = w[i];
 221       fCutValueI = i;
 222     }
 223   }
 224
 225   // if the best probability is too small, the cut is failed anyway
 226   if (!OkRangeD()) return kFALSE;
 227
 228   // if the best probability is OK, the cut is passed
 229   // if it correspond to the right particle
 230   return OkValue();
 231 }
 232
 233 void AliRsnCutPID::IncludeDetector(EDetector det, Double_t threshold, Bool_t goAbove)
 234 {
 235 //
 236 // Include a detector for a customized weight computing
 237 // and specify also its eventual threshold and if the detector
 238 // must be used above or below the threshold.
 239 // By default the threshold is zero and detector is always used.
 240 //
 241
 242   if (!CheckBounds(det)) return;
 243
 244   fUseDetector[det] = kTRUE;
 245   fPtThreshold[det] = threshold;
 246   fGoAboveThreshold[det] = goAbove;
 247 }