PWG2/RESONANCES/AliRsnPID.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //-------------------------------------------------------------------------
  17 //                      Class AliRsnPID
  18 //                     -------------------
  19 //           Simple collection of reconstructed tracks
  20 //           selected from an ESD event
  21 //           to be used for analysis.
  22 //           .........................................
  23 //
  24 // author: A. Pulvirenti             (email: alberto.pulvirenti@ct.infn.it)
  25 //-------------------------------------------------------------------------
  26
  27 #include <TMath.h>
  28 #include <TString.h>
  29 #include <TClonesArray.h>
  30 #include <TDatabasePDG.h>
  31
  32 #include "AliLog.h"
  33 #include "AliRsnMCInfo.h"
  34 #include "AliRsnDaughter.h"
  35 #include "AliRsnEvent.h"
  36
  37 #include "AliRsnPID.h"
  38
  39 ClassImp(AliRsnPID)
  40
  41 const char* AliRsnPID::fgkParticleNameLong[AliRsnPID::kSpecies + 1] =
  42 {
  43   "electron",
  44   "muon",
  45   "pion",
  46   "kaon",
  47   "proton",
  48   "unknown"
  49 };
  50
  51 const char* AliRsnPID::fgkParticleNameShort[AliRsnPID::kSpecies + 1] =
  52 {
  53   "e",
  54   "mu",
  55   "pi",
  56   "K",
  57   "p",
  58   "unk"
  59 };
  60
  61 const char* AliRsnPID::fgkParticleNameLatex[AliRsnPID::kSpecies + 1] =
  62 {
  63   "e",
  64   "#mu",
  65   "#pi",
  66   "K",
  67   "p",
  68   "?"
  69 };
  70
  71 const Int_t AliRsnPID::fgkParticlePDG[AliRsnPID::kSpecies + 1] =
  72 {
  73     11,
  74     13,
  75     211,
  76     321,
  77     2212,
  78     0
  79 };
  80
  81 //_____________________________________________________________________________
  82 AliRsnPID::AliRsnPID() :
  83   fMethod(kRealistic),
  84   fMaxPt(10.0),
  85   fMinProb(0.0)
  86 {
  87 //
  88 // Default constructor
  89 // Sets a default setup:
  90 //  - realistic PID
  91 //  - no lower limit for probability
  92 //  - upper limit of 10 GeV for Pt
  93 //  - suitable values for prior probabilities
  94 //
  95
  96     fPrior[kElectron] = 0.20;
  97         fPrior[kMuon    ] = 0.20;
  98         fPrior[kPion    ] = 0.83;
  99         fPrior[kKaon    ] = 0.07;
 100         fPrior[kProton  ] = 0.06;
 101 }
 102
 103 //_____________________________________________________________________________
 104 AliRsnPID::AliRsnPID(const AliRsnPID &event) :
 105   TObject((TObject)event),
 106   fMethod(event.fMethod),
 107   fMaxPt(event.fMaxPt),
 108   fMinProb(event.fMinProb)
 109 {
 110 //
 111 // Copy constructor.
 112 // Implemented to manage the array safely.
 113 //
 114
 115     Int_t i;
 116     for (i = 0; i < kSpecies; i++) fPrior[i] = event.fPrior[i];
 117 }
 118
 119 //_____________________________________________________________________________
 120 AliRsnPID& AliRsnPID::operator=(const AliRsnPID &event)
 121 {
 122 //
 123 // Assignment operator.
 124 // Implemented to manage the array safely.
 125 //
 126
 127     fMethod = event.fMethod;
 128     fMaxPt = event.fMaxPt;
 129     fMinProb = event.fMinProb;
 130
 131     Int_t i;
 132     for (i = 0; i < kSpecies; i++) fPrior[i] = event.fPrior[i];
 133
 134     // return the newly created object
 135     return (*this);
 136 }
 137
 138 //_____________________________________________________________________________
 139 AliRsnPID::EType AliRsnPID::InternalType(Int_t pdg)
 140 //
 141 // Return the internal enum value corresponding to the PDG
 142 // code passed as argument, if possible.
 143 // Otherwise, returns 'kUnknown' by default.
 144 //
 145 {
 146     EType value;
 147     Int_t absPDG = TMath::Abs(pdg);
 148
 149         switch (absPDG) {
 150         case 11:
 151             value = kElectron;
 152             break;
 153         case 13:
 154             value = kMuon;
 155             break;
 156         case 211:
 157             value = kPion;
 158             break;
 159         case 321:
 160             value = kKaon;
 161             break;
 162         case 2212:
 163             value = kProton;
 164             break;
 165         default:
 166             value = kUnknown;
 167     }
 168     return value;
 169 }
 170
 171
 172 //_____________________________________________________________________________
 173 Int_t AliRsnPID::PDGCode(EType type)
 174 {
 175 //
 176 // Returns the PDG code of the particle type
 177 // specified as argument (w.r. to the internal enum)
 178 //
 179
 180     if (type >= kElectron && type <= kUnknown) {
 181         return fgkParticlePDG[type];
 182     }
 183     else {
 184         return 0;
 185     }
 186 }
 187
 188 //_____________________________________________________________________________
 189 const char* AliRsnPID::ParticleName(EType type, Bool_t shortName)
 190 {
 191 //
 192 // Returns the name of the particle type
 193 // specified as argument (w.r. to the internal enum)
 194 //
 195
 196     if (type >= kElectron && type <= kUnknown) {
 197         return shortName ? fgkParticleNameShort[type] : fgkParticleNameLong[type];
 198     }
 199     else {
 200         return shortName ? "unk" : "unknown";
 201     }
 202 }
 203
 204 //_____________________________________________________________________________
 205 const char* AliRsnPID::ParticleNameLatex(EType type)
 206 {
 207 //
 208 // Returns the name of the particle type
 209 // specified as argument (w.r. to the internal enum)
 210 //
 211
 212     if (type >= kElectron && type <= kUnknown) {
 213         return fgkParticleNameLatex[type];
 214     }
 215     else {
 216         return "?";
 217     }
 218 }
 219
 220 //_____________________________________________________________________________
 221 Double_t AliRsnPID::ParticleMass(EType type)
 222 {
 223 //
 224 // Returns the mass corresponding to the particle type
 225 // specified as argument (w.r. to the internal enum)
 226 //
 227     TDatabasePDG *db = TDatabasePDG::Instance();
 228     Int_t pdg = PDGCode(type);
 229     return db->GetParticle(pdg)->Mass();
 230 }
 231
 232 //_____________________________________________________________________________
 233 Bool_t AliRsnPID::IdentifyRealistic(AliRsnDaughter *daughter)
 234 {
 235 //
 236 // Uses the Bayesian combination of prior probabilities
 237 // with the PID weights of the passed object to compute
 238 // the overall PID probabilities for each particle type.
 239 //
 240 // Once this computation is done, the argument is assigned
 241 // the PID corresponding to the largest probability,
 242 // and its data members are updated accordingly.
 243 // If the track Pt is larger than the cut defined (fMaxPt)
 244 // or the probability is smaller than the cut defined (fMinProb),
 245 // the track is considered unidentified.
 246 //
 247 // If the identification procedure encounters errors,
 248 // the return value will be "FALSE", otherwise it is "TRUE".
 249 //
 250
 251     // reset all PID probabilities to 0.0
 252     Unidentify(daughter, 0.0);
 253
 254     // check the transverse momentum:
 255     // if it is larger than the cut, the particle is unidentified
 256     // setting all its probabilities to 1.0
 257     if (daughter->Pt() > fMaxPt) {
 258         Unidentify(daughter, 1.0);
 259         return kTRUE;
 260     }
 261
 262     // 0 - retrieve PID weights from argument
 263     Int_t     i;
 264     Double_t  sum    = 0.0;
 265     Double_t *prob   = new Double_t[kSpecies];
 266     Double_t *weight = new Double_t[kSpecies];
 267     for (i = 0; i < kSpecies; i++) weight[i] = (daughter->PID())[i];
 268
 269     // step 1 - compute the normalization factor
 270     for (i = 0; i < kSpecies; i++) {
 271         prob[i] = fPrior[i] * weight[i];
 272         sum += prob[i];
 273     }
 274     if (sum <= (Double_t)0.) {
 275         AliError(Form("Sum of weights = %f < 0", sum));
 276         return kFALSE;
 277     }
 278
 279         // step 2 - normalize PID weights and update daughter data-member
 280     for (i = 0; i < kSpecies; i++) {
 281         prob[i] /= sum;
 282         daughter->SetPIDProb(i, prob[i]);
 283     }
 284
 285         // step 3 - find the maximum probability and update daughter data members
 286     Int_t    imax = 0;
 287     Double_t pmax = prob[0];
 288     for (i = 1; i < kSpecies; i++) {
 289         if (prob[i] > pmax) {
 290             imax = i;
 291             pmax = prob[i];
 292         }
 293     }
 294     EType type = (EType)imax;
 295     if (pmax >= fMinProb) {
 296         daughter->SetPIDType(type);
 297         daughter->SetM(ParticleMass(type));
 298     }
 299     else {
 300         daughter->SetPIDType(kUnknown);
 301     }
 302
 303     return kTRUE;
 304 }
 305
 306 //_____________________________________________________________________________
 307 Bool_t AliRsnPID::IdentifyPerfect(AliRsnDaughter *daughter)
 308 {
 309 //
 310 // Uses the true PDG code to make a perfect identification.
 311 // If the true PDG code does not correspond to any
 312 // of the expected PID types, gives a warning, sets the
 313 // PID to 'unknown' and returns kFALSE.
 314 // Otherwise, returns kTRUE.
 315 //
 316
 317     // reset all PID probabilities to 0.0
 318     Unidentify(daughter, 0.0);
 319
 320     // if the MCInfo is not present, the particle cannot be identified perfectly
 321     // in this case, to notice the error, the probs are maintained all to 0.0
 322     AliRsnMCInfo *particle = daughter->GetMCInfo();
 323     if (!particle) {
 324         AliWarning("Particle object not initialized: impossible to do perfect PID");
 325         return kFALSE;
 326     }
 327
 328     // convert the PDG into the internal enum
 329     Int_t pdgCode = particle->PDG();
 330     EType type = InternalType(pdgCode);
 331
 332     // if the type is one of the available ones in the PID enum
 333     // (e, mu, pi, K, p)
 334     // the corresponding probability is set to 1, and the other remain 0
 335     if (type >= 0 && type < kSpecies) {
 336         daughter->SetPIDType(type);
 337         daughter->SetPIDProb(type, 1.0);
 338         daughter->SetM(ParticleMass(type));
 339         return kTRUE;
 340     }
 341
 342     return kFALSE;
 343 }
 344
 345 //_____________________________________________________________________________
 346 Bool_t AliRsnPID::Unidentify(AliRsnDaughter *daughter, Double_t value)
 347 {
 348 //
 349 // Sets the PID to 'unknown' to every track.
 350 //
 351
 352     Int_t i;
 353     for (i = 0; i < kSpecies; i++) daughter->SetPIDProb(i, value);
 354     daughter->SetPIDType(kUnknown);
 355     return kTRUE;
 356 }
 357
 358 //_____________________________________________________________________________
 359 Bool_t AliRsnPID::Identify(AliRsnDaughter *daughter)
 360 {
 361 //
 362 // Recalls one of the above methods, according to the one
 363 // defined in the related data member.
 364 // If the method is not recognized, returns kFALSE and
 365 // gives an alert. Otherwise, returns kTRUE.
 366 //
 367
 368     switch (fMethod) {
 369         case kNone:
 370             Unidentify(daughter);
 371             return kTRUE;
 372         case kRealistic:
 373             IdentifyRealistic(daughter);
 374             return kTRUE;
 375         case kPerfect:
 376             IdentifyPerfect(daughter);
 377             return kTRUE;
 378         default:
 379             AliError(Form("PID method '%d' unrecognized. Nothing done.", fMethod));
 380             return kFALSE;
 381     }
 382 }
 383
 384 //_____________________________________________________________________________
 385 Bool_t AliRsnPID::IdentifiedAs(AliRsnDaughter *d, EType type, Short_t charge)
 386 {
 387 //
 388 // Tells if a particle has been identified to be of a given tipe and charge.
 389 // If the charge is zero, the check is done only on the PID type, otherwise
 390 // both charge and PID type are required to match
 391 //
 392     if (charge == 0) {
 393         return (d->PIDType() == type);
 394     }
 395     else if (charge > 0) {
 396         return (d->PIDType() == type && d->Charge() > 0);
 397     }
 398     else {
 399         return (d->PIDType() == type && d->Charge() < 0);
 400     }
 401 }
 402
 403 //_____________________________________________________________________________
 404 Bool_t AliRsnPID::Identify(AliRsnEvent *event)
 405 {
 406 //
 407 // Performs identification for all tracks in a given event.
 408 // Returns the logical AND of all PID operations.
 409 //
 410
 411     Bool_t check = kTRUE;
 412     if (!event) return check;
 413     if (!event->GetTracks()) return check;
 414     if (event->GetTracks()->IsEmpty()) return check;
 415
 416     AliRsnDaughter *daughter = 0;
 417     TObjArrayIter iter(event->GetTracks());
 418     while ( (daughter = (AliRsnDaughter*)iter.Next()) ) {
 419         check = check && Identify(daughter);
 420     }
 421     event->FillPIDArrays();
 422
 423     return check;
 424 }
 425
 426
 427 //_____________________________________________________________________________
 428 void AliRsnPID::SetPriorProbability(EType type, Double_t p)
 429 {
 430 //
 431 // Sets the prior probability for Realistic PID, for a
 432 // given particle species.
 433 //
 434
 435     if (type >= kElectron && type < kSpecies) {
 436         fPrior[type] = p;
 437     }
 438 }