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()) perfectPID = kFALSE;
 117   if (!daughter->GetRefESDtrack()) useDefault = kTRUE;
 118   if (!daughter->GetRefESDtrack() && !daughter->GetRefAODtrack()) return kFALSE;
 119
 120   // if perfect PID ise required, this overcomes all
 121   // in this case the weight of the correct species is set to 1
 122   // and the others to 0
 123   // of course this happens only if there is a reference MC
 124   if (perfectPID)
 125   {
 126     j = TMath::Abs(daughter->GetRefMC()->Particle()->GetPdgCode());
 127     for (i = 0; i < AliPID::kSPECIES; i++)
 128     {
 129       if (AliPID::ParticleCode((AliPID::EParticleType)i) == j) fWeight[i] = 1.0;
 130       else fWeight[i] = 0.0;
 131     }
 132     return kTRUE;
 133   }
 134
 135   // if default weights are (or need to be) used,
 136   // they are taken directly and function exits
 137   if (useDefault)
 138   {
 139     if (daughter->GetRefESDtrack())
 140       daughter->GetRefESDtrack()->GetESDpid(fWeight);
 141     else
 142     {
 143       for (i = 0; i < AliPID::kSPECIES; i++)
 144         fWeight[i] = daughter->GetRefAODtrack()->PID()[i];
 145     }
 146     return kTRUE;
 147   }
 148
 149   // if we arrive here, this means that we have an ESD track
 150   // and we want to customize the PID
 151   AliESDtrack *track = daughter->GetRefESDtrack();
 152   Double_t     w[kDetectors][AliPID::kSPECIES];
 153   track->GetITSpid(w[kITS]);
 154   track->GetTPCpid(w[kTPC]);
 155   track->GetTRDpid(w[kTRD]);
 156   track->GetTOFpid(w[kTOF]);
 157   track->GetHMPIDpid(w[kHMPID]);
 158
 159   // if a detector is excluded or the track has not the right pT
 160   // all related weights are set to 1 in order not to contribute
 161   // to final multiplication
 162   for (i = 0; i < kDetectors; i++)
 163   {
 164     if (!fUseDetector[i] || !CheckThreshold((EDetector)i, track->Pt()))
 165     {
 166       for (j = 0; j < AliPID::kSPECIES; j++) {
 167         w[i][j] = 1.0;
 168       }
 169     }
 170   }
 171
 172   // multiplicate all weights to compute final one
 173   for (i = 0; i < AliPID::kSPECIES; i++)
 174   {
 175     fWeight[i] = w[kITS][i] * w[kTPC][i] * w[kTRD][i] * w[kTOF][i] * w[kHMPID][i];
 176   }
 177
 178   return kTRUE;
 179 }
 180
 181 //_________________________________________________________________________________________________
 182 Bool_t AliRsnCutPID::IsSelected(TObject *obj1, TObject* /*obj2*/)
 183 {
 184 //
 185 // Cut checker.
 186 //
 187
 188   // coherence check
 189   if (!AliRsnCut::TargetOK(obj1))
 190   {
 191     AliError(Form("Wrong target. Skipping cut", GetName()));
 192     return kTRUE;
 193   }
 194
 195   // convert the object into the unique correct type
 196   AliRsnDaughter *daughter = dynamic_cast<AliRsnDaughter*>(obj1);
 197
 198   // try to compute the weights
 199   if (!ComputeWeights(daughter)) return kFALSE;
 200
 201   // combine with priors and get the majority
 202   Int_t    i;
 203   Double_t sum = 0.0, w[AliPID::kSPECIES];
 204   for (i = 0; i < AliPID::kSPECIES; i++)
 205   {
 206     w[i] = fWeight[i] * fPrior[i];
 207     sum += w[i];
 208   }
 209   if (sum <= 0.0)
 210   {
 211     AliError("Sum = 0");
 212     return kFALSE;
 213   }
 214   for (i = 0; i < AliPID::kSPECIES; i++) w[i] /= sum;
 215
 216   // find the largest probability and related PID
 217   // and assign them to the mother class members which
 218   // are checked for the cut
 219   fCutValueI = 0;
 220   fCutValueD = w[0];
 221   for (i = 1; i < AliPID::kSPECIES; i++)
 222   {
 223     if (w[i] > fCutValueD)
 224     {
 225       fCutValueD = w[i];
 226       fCutValueI = i;
 227     }
 228   }
 229
 230   // if the best probability is too small, the cut is failed anyway
 231   if (!OkRangeD()) return kFALSE;
 232
 233   // if the best probability is OK, the cut is passed
 234   // if it correspond to the right particle
 235   return OkValue();
 236 }
 237
 238 void AliRsnCutPID::IncludeDetector(EDetector det, Double_t threshold, Bool_t goAbove)
 239 {
 240 //
 241 // Include a detector for a customized weight computing
 242 // and specify also its eventual threshold and if the detector
 243 // must be used above or below the threshold.
 244 // By default the threshold is zero and detector is always used.
 245 //
 246
 247   if (!CheckBounds(det)) return;
 248
 249   fUseDetector[det] = kTRUE;
 250   fPtThreshold[det] = threshold;
 251   fGoAboveThreshold[det] = goAbove;
 252 }