PWG2/RESONANCES/AliRsnCutBetheBloch.cxx

   1 //
   2 // Class AliRsnCutBetheBloch
   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 "AliRsnDaughter.h"
  21 #include "AliRsnCutBetheBloch.h"
  22
  23 ClassImp(AliRsnCutBetheBloch)
  24
  25 //_________________________________________________________________________________________________
  26 AliRsnCutBetheBloch::AliRsnCutBetheBloch() :
  27     AliRsnCut(),
  28     fCorrect(kTRUE),
  29     fMIP(50.0),
  30     fType(AliPID::kUnknown)
  31 {
  32 //
  33 // Default constructor.
  34 //
  35
  36   fConst[0] = fConst[1] = fConst[2] = fConst[3] = fConst[4] = 0.0;
  37 }
  38
  39 //_________________________________________________________________________________________________
  40 AliRsnCutBetheBloch::AliRsnCutBetheBloch
  41 (const char *name, Double_t fractionRange, AliPID::EParticleType type, Double_t mip, Bool_t correct) :
  42     AliRsnCut(name, 0.0, fractionRange),
  43     fCorrect(correct),
  44     fMIP(mip),
  45     fType(type)
  46 {
  47 //
  48 // Main constructor.
  49 // the cut range is the relative fraction of the value:
  50 // BB*(1-fraction) < TPC < BB*(1+fraction)
  51 // which means:
  52 // -fraction < (TPC - BB)/BB < fraction
  53 //
  54
  55   fConst[0] = fConst[1] = fConst[2] = fConst[3] = fConst[4] = 0.0;
  56 }
  57
  58 //_____________________________________________________________________________
  59 Double_t AliRsnCutBetheBloch::BetheBloch(AliRsnDaughter * const track)
  60 {
  61 //
  62 // Computes the theoretical dE/dx according to
  63 // a given mass hypothesis, from which betaGamma is computed
  64 //
  65 // This is the empirical ALEPH parameterization of the Bethe-Bloch formula.
  66 // It is normalized to 1 at the minimum.
  67 //
  68 // The default values for the kp* parameters are for ALICE TPC.
  69 // The value is computed in MIP units, multiplied by 50 to have it in energy.
  70 //
  71
  72   AliPID pid;
  73   Double_t mass = pid.ParticleMass(fType);
  74
  75   Double_t betaGamma = track->P() / mass;
  76   Double_t beta = betaGamma / TMath::Sqrt(1.0 + betaGamma * betaGamma);
  77   Double_t aa = TMath::Power(beta, fConst[3]);
  78   Double_t bb = TMath::Power(1.0 / betaGamma, fConst[4]);
  79
  80   bb = TMath::Log(fConst[2] + bb);
  81
  82   Double_t out = (fConst[1] - aa - bb) * fConst[0] / aa;
  83
  84   if (fCorrect) {
  85     Double_t kMeanCorr = 0.1;
  86     Double_t meanCorr = (1 + (out - 1) * kMeanCorr);
  87     out *= meanCorr;
  88   }
  89
  90   return out * fMIP;
  91 }
  92
  93 //_____________________________________________________________________________
  94 Double_t AliRsnCutBetheBloch::RelDiff(AliRsnDaughter *track)
  95 {
  96 //
  97 // Relative difference between BB value and TPC signal
  98 //
  99
 100   if (!track->GetRefESD()) return -99999.9;
 101
 102   // compute Bethe-Bloch with the given mass hypothesis
 103   Double_t bb = BetheBloch(track);
 104   return TMath::Abs((track->GetRefESD()->GetTPCsignal() - bb) / bb);
 105 }
 106
 107 //_________________________________________________________________________________________________
 108 Bool_t AliRsnCutBetheBloch::IsSelected(ETarget tgt, AliRsnDaughter *track)
 109 {
 110 //
 111 // Cut checker.
 112 //
 113
 114   // coherence check
 115   if (tgt != AliRsnCut::kParticle) {
 116     AliError(Form("Wrong target. Skipping cut", GetName()));
 117     return kTRUE;
 118   }
 119
 120   // if the required PID of the track is not the same as the
 121   // reference of the cut, the cut is automatically skipped
 122   if (track->RequiredPID() != fType) return kTRUE;
 123
 124   // retrieve the TPC signal
 125   AliESDtrack *esd = track->GetRefESD();
 126   if (!esd) {
 127     AliError("ESD information unavailable");
 128     return kTRUE;
 129   }
 130
 131   // the cut range is the relative fraction of the value:
 132   // BB*(1-fraction) < TPC < BB*(1+fraction)
 133   // which means:
 134   // -fraction < (TPC - BB)/BB < fraction
 135   // so we must compute the cut value accordingly
 136   fCutValueD = RelDiff(track);
 137
 138   // then, this cut is checked inside the range
 139   return OkRange();
 140 }
 141
 142 //_________________________________________________________________________________________________
 143 Bool_t AliRsnCutBetheBloch::IsSelected(ETarget /*tgt*/, AliRsnPairParticle* /*pair*/)
 144 {
 145 //
 146 // Cut checker
 147 //
 148
 149   AliWarning("Cannot apply this cut to pairs");
 150   return kTRUE;
 151 }
 152
 153 //_________________________________________________________________________________________________
 154 Bool_t AliRsnCutBetheBloch::IsSelected(ETarget /*tgt*/, AliRsnEvent* /*event*/)
 155 {
 156 //
 157 // Cut checker
 158 //
 159
 160   AliWarning("Cannot apply this cut to events");
 161   return kTRUE;
 162 }
 163
 164 //_________________________________________________________________________________________________
 165 Bool_t AliRsnCutBetheBloch::IsSelected(ETarget /*tgt*/, AliRsnEvent* /*ev1*/, AliRsnEvent* /*ev2*/)
 166 {
 167 //
 168 // Cut checker
 169 //
 170
 171   AliWarning("Cannot apply this cut to event mixing");
 172   return kTRUE;
 173 }