PWG2/RESONANCES/AliRsnPID.cxx

   1 //
   2 // Class AliRsnPID
   3 //
   4 // Processes the tracks to compute PID probabilities,
   5 // by combining the prior probabilities with the PID weights
   6 // stored in each AliRsnDaughter
   7 //
   8 // authors: A. Pulvirenti (alberto.pulvirenti@ct.infn.it)
   9 //          M. Vala (martin.vala@cern.ch)
  10 //
  11
  12
  13 #include <TMath.h>
  14 #include <TDirectory.h>
  15
  16 #include "AliLog.h"
  17 #include "AliRsnMCInfo.h"
  18 #include "AliRsnDaughter.h"
  19 #include "AliRsnEvent.h"
  20
  21 #include "AliRsnPID.h"
  22
  23 ClassImp(AliRsnPID)
  24
  25 const Double_t AliRsnPID::fgkParticleMass[AliRsnPID::kSpeciesAll + 1] =
  26 {
  27   0.00051099892,   // electron
  28   0.105658369,     // muon
  29   0.13957018,      // pion
  30   0.493677,        // kaon
  31   0.93827203,      // proton
  32   0.0,             // nothing
  33   1.019413,        // phi
  34   3.09693          // jpsi
  35 };
  36
  37 const char* AliRsnPID::fgkParticleNameLong[AliRsnPID::kSpeciesAll + 1] =
  38 {
  39   "electron",
  40   "muon",
  41   "pion",
  42   "kaon",
  43   "proton",
  44   "unknown",
  45   "phi",
  46   "jpsi"
  47 };
  48
  49 const char* AliRsnPID::fgkParticleNameShort[AliRsnPID::kSpeciesAll + 1] =
  50 {
  51   "e",
  52   "mu",
  53   "pi",
  54   "K",
  55   "p",
  56   "unk",
  57   "phi",
  58   "jpsi"
  59 };
  60
  61 const char* AliRsnPID::fgkParticleNameLatex[AliRsnPID::kSpeciesAll + 1] =
  62 {
  63   "e",
  64   "#mu",
  65   "#pi",
  66   "K",
  67   "p",
  68   "?",
  69   "#phi",
  70   "J/#psi"
  71 };
  72
  73 const Int_t AliRsnPID::fgkParticlePDG[AliRsnPID::kSpeciesAll + 1] =
  74 {
  75   11,
  76   13,
  77   211,
  78   321,
  79   2212,
  80   0,
  81   333,
  82   443
  83 };
  84
  85 //_____________________________________________________________________________
  86 AliRsnPID::AliRsnPID() :
  87   TNamed("RsnPID", ""),
  88   fMaxPt(100.0),
  89   fMinProb(0.0),
  90   fPIDArraysSize(1000)
  91 {
  92 //
  93 // Constructor.
  94 // Adds also the object in the default directory of the session.
  95 //
  96
  97   Int_t i;
  98   for (i = 0; i < kSpecies; i++) fPrior[i] = 0.2;
  99 }
 100
 101 //_____________________________________________________________________________
 102 AliRsnPID::EType AliRsnPID::InternalType(Int_t pdg)
 103 //
 104 // Return the internal enum value corresponding to the PDG
 105 // code passed as argument, if possible.
 106 // Otherwise, returns 'kUnknown' by default.
 107 //
 108 {
 109   EType value;
 110   Int_t absPDG = TMath::Abs(pdg);
 111
 112   switch (absPDG)
 113   {
 114     case 11:
 115       value = kElectron;
 116       break;
 117     case 13:
 118       value = kMuon;
 119       break;
 120     case 211:
 121       value = kPion;
 122       break;
 123     case 321:
 124       value = kKaon;
 125       break;
 126     case 2212:
 127       value = kProton;
 128       break;
 129     case 333:
 130       value = kPhi;
 131       break;
 132     case 443:
 133       value = kJPsi;
 134       break;
 135     default:
 136       value = kUnknown;
 137   }
 138   return value;
 139 }
 140
 141
 142 //_____________________________________________________________________________
 143 Int_t AliRsnPID::PDGCode(EType type)
 144 {
 145 //
 146 // Returns the PDG code of the particle type
 147 // specified as argument (w.r. to the internal enum)
 148 //
 149
 150   if (type >= kElectron && type <= kUnknown)
 151   {
 152     return fgkParticlePDG[type];
 153   }
 154   else
 155   {
 156     return 0;
 157   }
 158 }
 159
 160 //_____________________________________________________________________________
 161 const char* AliRsnPID::ParticleName(EType type, Bool_t shortName)
 162 {
 163 //
 164 // Returns the name of the particle type
 165 // specified as argument (w.r. to the internal enum)
 166 //
 167
 168   if (type >= 0 && type <= kSpeciesAll)
 169   {
 170     return shortName ? fgkParticleNameShort[type] : fgkParticleNameLong[type];
 171   }
 172   else
 173   {
 174     return shortName ? "unk" : "unknown";
 175   }
 176 }
 177
 178 //_____________________________________________________________________________
 179 const char* AliRsnPID::ParticleNameLatex(EType type)
 180 {
 181 //
 182 // Returns the name of the particle type
 183 // specified as argument (w.r. to the internal enum)
 184 //
 185
 186   if (type >= kElectron && type <= kSpeciesAll)
 187   {
 188     return fgkParticleNameLatex[type];
 189   }
 190   else
 191   {
 192     return "?";
 193   }
 194 }
 195
 196 //_____________________________________________________________________________
 197 Double_t AliRsnPID::ParticleMass(EType type)
 198 {
 199 //
 200 // Returns the mass corresponding to the particle type
 201 // specified as argument (w.r. to the internal enum)
 202 //
 203   /*
 204   TDatabasePDG *db = TDatabasePDG::Instance();
 205   Int_t pdg = PDGCode(type);
 206   return db->GetParticle(pdg)->Mass();
 207   */
 208   if (type >= kElectron && type < kSpeciesAll) return fgkParticleMass[type];
 209   return 0.0;
 210 }
 211
 212 //_____________________________________________________________________________
 213 Bool_t AliRsnPID::ComputeProbs(AliRsnDaughter *daughter)
 214 {
 215 //
 216 // Uses the Bayesian combination of prior probabilities
 217 // with the PID weights of the passed object to compute
 218 // the overall PID probabilities for each particle type.
 219 //
 220 // Once this computation is done, the argument is assigned
 221 // the PID corresponding to the largest probability,
 222 // and its data members are updated accordingly.
 223 // If the track Pt is larger than the cut defined (fMaxPt)
 224 // or the probability is smaller than the cut defined (fMinProb),
 225 // the track is considered unidentified.
 226 //
 227 // If the identification procedure encounters errors,
 228 // the return value will be "FALSE", otherwise it is "TRUE".
 229 //
 230
 231   // reset all PID probabilities to 0.0
 232   Int_t i;
 233   for (i = 0; i < kSpecies; i++) daughter->SetPIDProb(i, 1.0 / (Double_t) kSpecies);
 234
 235   // multiply weights and priors
 236   Double_t sum = 0.0, prob[kSpecies];
 237   for (i = 0; i < kSpecies; i++)
 238   {
 239     prob[i] = fPrior[i] * daughter->PID()[i];
 240     sum += prob[i];
 241   }
 242   if (sum <= (Double_t) 0.)
 243   {
 244     AliError(Form("Sum of weights = %f <= 0", sum));
 245     return kFALSE;
 246   }
 247
 248   // normalize
 249   for (i = 0; i < kSpecies; i++)
 250   {
 251     prob[i] /= sum;
 252     daughter->SetPIDProb(i, prob[i]);
 253   }
 254
 255   daughter->AssignRealisticPID();
 256   Double_t assprob;
 257   AliRsnDaughter::SetPIDMethod(AliRsnDaughter::kRealistic);
 258   AliRsnPID::EType type = daughter->PIDType(assprob);
 259   AliDebug(5, Form("Assigned PID: %s [%.2f %]", AliRsnPID::ParticleName(type), assprob*100.));
 260
 261   return kTRUE;
 262 }
 263
 264 //_____________________________________________________________________________
 265 Bool_t AliRsnPID::IdentifiedAs(AliRsnDaughter *d, EType type, Short_t charge)
 266 {
 267 //
 268 // Tells if a particle has can be identified to be of a given tipe and charge.
 269 // If the charge is zero, the check is done only on the PID type, otherwise
 270 // both charge and PID type are required to match.
 271 // If the track momentum is larger than the pt threshold passed to this object,
 272 // or the maximum probability is smaller than the prob thrashold, the return value
 273 // is kFALSE even when the type and charge are matched.
 274 //
 275
 276   EType dType = TrackType(d);
 277   if (dType != type) return kFALSE;
 278   if (charge == 0)
 279   {
 280     return kTRUE;
 281   }
 282   else if (charge > 0)
 283   {
 284     return (d->Charge() > 0);
 285   }
 286   else
 287   {
 288     return (d->Charge() < 0);
 289   }
 290 }
 291
 292 //_____________________________________________________________________________
 293 AliRsnPID::EType AliRsnPID::TrackType(AliRsnDaughter *d)
 294 {
 295 //
 296 // Returns the track type according to the object settings
 297 // and to the static settings in the AliRsnDaughter object.
 298 //
 299
 300   Double_t prob;
 301   EType type = d->PIDType(prob);
 302
 303   if (d->Pt() > fMaxPt) return kUnknown;
 304   if (prob < fMinProb) return kUnknown;
 305
 306   return type;
 307 }
 308
 309 //_____________________________________________________________________________
 310 Bool_t AliRsnPID::Process(AliRsnEvent *event)
 311 {
 312 //
 313 // Performs identification for all tracks in a given event.
 314 // Returns the logical AND of all PID operations.
 315 //
 316
 317   Bool_t check = kTRUE;
 318   if (!event) return check;
 319   if (!event->GetTracks()) return check;
 320   if (event->GetTracks()->IsEmpty()) return check;
 321
 322   AliRsnDaughter *daughter = 0;
 323   TObjArrayIter iter(event->GetTracks());
 324   while ((daughter = (AliRsnDaughter*) iter.Next()))
 325   {
 326     check = check && ComputeProbs(daughter);
 327   }
 328   event->FillPIDArrays(fPIDArraysSize);
 329
 330   return check;
 331 }
 332
 333
 334 //_____________________________________________________________________________
 335 void AliRsnPID::SetPriorProbability(EType type, Double_t p)
 336 {
 337 //
 338 // Sets the prior probability for Realistic PID, for a
 339 // given particle species.
 340 //
 341
 342   if (type >= kElectron && type < kSpecies)
 343   {
 344     fPrior[type] = p;
 345   }
 346 }
 347
 348 //_____________________________________________________________________________
 349 void AliRsnPID::DumpPriors()
 350 {
 351 //
 352 // Print all prior probabilities
 353 //
 354
 355   Int_t i;
 356   Char_t msg[200];
 357
 358   for (i = 0; i < kSpecies; i++)
 359   {
 360     sprintf(msg, "Prior probability for '%s' = %3.5f", fgkParticleNameLong[i], fPrior[i]);
 361     AliInfo(msg);
 362   }
 363 }