STEER/AliPIDCombined.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 //-----------------------------------------------------------------
  18 //        Base class for combining PID of different detectors    //
  19 //        (user selected) and compute Bayesian probabilities     //
  20 //                                                               //
  21 //                                                               //
  22 //   Origin: Pietro Antonioli, INFN-BO Pietro.Antonioli@cern.ch  //
  23 //                                                               //
  24 //-----------------------------------------------------------------
  25
  26 #include <TH1.h>
  27
  28 #include <AliVTrack.h>
  29 #include <AliLog.h>
  30 #include <AliPID.h>
  31 #include <AliPIDResponse.h>
  32
  33 #include "AliPIDCombined.h"
  34
  35 AliPIDCombined::AliPIDCombined():
  36         TNamed("CombinedPID","CombinedPID"),
  37         fDetectorMask(0),
  38         fRejectMismatchMask(0x7F),
  39         fEnablePriors(kTRUE),
  40         fSelectedSpecies(AliPID::kSPECIES)
  41 {
  42   //
  43   // default constructor
  44   //
  45   for (Int_t i=0;i<AliPID::kSPECIES+AliPID::kSPECIESLN;i++) fPriorsDistributions[i]=NULL;
  46   AliLog::SetClassDebugLevel("AliPIDCombined",10);
  47 }
  48
  49 //-------------------------------------------------------------------------------------------------
  50 AliPIDCombined::AliPIDCombined(const TString& name,const TString& title):
  51         TNamed(name,title),
  52         fDetectorMask(0),
  53         fRejectMismatchMask(0x7F),
  54         fEnablePriors(kTRUE),
  55         fSelectedSpecies(AliPID::kSPECIES)
  56 {
  57   //
  58   // default constructor with name and title
  59   //
  60   for (Int_t i=0;i<AliPID::kSPECIES+AliPID::kSPECIESLN;i++) fPriorsDistributions[i]=NULL;
  61   AliLog::SetClassDebugLevel("AliPIDCombined",10);
  62
  63 }
  64
  65 //-------------------------------------------------------------------------------------------------
  66 AliPIDCombined::~AliPIDCombined() {
  67
  68   for(Int_t i=0;i < AliPID::kSPECIES+AliPID::kSPECIESLN;i++){
  69     if(fPriorsDistributions[i])
  70       delete fPriorsDistributions[i];
  71   }
  72 }
  73
  74 //-------------------------------------------------------------------------------------------------
  75 void AliPIDCombined::SetPriorDistribution(AliPID::EParticleType type,TH1F *prior) {
  76   if ( (type < 0) || ( type >= (AliPID::kSPECIESN+AliPID::kSPECIESLN) ) ){
  77     AliError(Form("Invalid EParticleType setting prior (offending type: %d)",type));
  78     return;
  79   }
  80   if(prior) {
  81     Int_t i = (Int_t)type;
  82     if (i >= AliPID::kSPECIES ) {
  83       if (i < (Int_t)AliPID::kDeuteron) {
  84         AliError(Form("Invalid EParticleType setting prior. Type: %d (neutral) not supported",i));
  85         return;
  86       } else i -= (AliPID::kSPECIESN-AliPID::kSPECIES);
  87     }
  88     if (i>(AliPID::kSPECIES+AliPID::kSPECIESLN)) {             // coverity fix (to put it mildly....)
  89       AliError(Form("Unexpected EParticleType setting prior. Type: %d (neutral) not supported",i));
  90       return;
  91     }
  92     if (fPriorsDistributions[i]) {
  93       delete fPriorsDistributions[i];
  94     }
  95     fPriorsDistributions[i]=new TH1F(*prior);
  96   }
  97 }
  98
  99 //-------------------------------------------------------------------------------------------------
 100 UInt_t AliPIDCombined::ComputeProbabilities(const AliVTrack *track, const AliPIDResponse *response, Double_t* bayesProbabilities) const {
 101         Double_t pITS[fSelectedSpecies];
 102         Double_t pTPC[fSelectedSpecies];
 103         Double_t pTRD[fSelectedSpecies];
 104         Double_t pTOF[fSelectedSpecies];
 105         Double_t pHMPID[fSelectedSpecies];
 106         Double_t pEMCAL[fSelectedSpecies];
 107         Double_t pPHOS[fSelectedSpecies];
 108         for (Int_t i=0;i<fSelectedSpecies;i++) {
 109          pITS[i]=1.;
 110          pTPC[i]=1.;
 111          pTRD[i]=1.;
 112          pTOF[i]=1.;
 113          pHMPID[i]=1.;
 114          pEMCAL[i]=1.;
 115          pPHOS[i]=1.;
 116         }
 117         UInt_t usedMask=0;
 118         AliPIDResponse::EDetPidStatus status=AliPIDResponse::kDetNoSignal;
 119         Double_t p[fSelectedSpecies];
 120
 121         // getting probability distributions for selected detectors only
 122         if (fDetectorMask & AliPIDResponse::kDetITS) {
 123           status = response->ComputeITSProbability(track, fSelectedSpecies, pITS);
 124           SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetITS,pITS);
 125         }
 126
 127         if (fDetectorMask & AliPIDResponse::kDetTPC) {
 128           status = response->ComputeTPCProbability(track, fSelectedSpecies, pTPC);
 129           SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetTPC,pTPC);
 130         }
 131
 132
 133         if (fDetectorMask & AliPIDResponse::kDetTRD) {
 134           status = response->ComputeTRDProbability(track, fSelectedSpecies, pTRD);
 135           SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetTRD,pTRD);
 136         }
 137
 138         if (fDetectorMask &  AliPIDResponse::kDetTOF) {
 139           status = response->ComputeTOFProbability(track, fSelectedSpecies, pTOF);
 140           SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetTOF,pTOF);
 141         }
 142
 143         if (fDetectorMask & AliPIDResponse::kDetHMPID) {
 144           status = response->ComputeHMPIDProbability(track, fSelectedSpecies, pHMPID);
 145           SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetHMPID,pHMPID);
 146         }
 147
 148
 149         if (fDetectorMask & AliPIDResponse::kDetEMCAL) {
 150           status = response->ComputeEMCALProbability(track, fSelectedSpecies, pEMCAL);
 151           SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetEMCAL,pEMCAL);
 152         }
 153
 154
 155         if (fDetectorMask & AliPIDResponse::kDetPHOS) {
 156           status = response->ComputePHOSProbability(track, fSelectedSpecies, pPHOS);
 157           SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetPHOS,pPHOS);
 158         }
 159
 160
 161         for (Int_t i =0; i<fSelectedSpecies; i++){
 162           p[i] = pITS[i]*pTPC[i]*pTRD[i]*pTOF[i]*pHMPID[i]*pEMCAL[i]*pPHOS[i];
 163         }
 164         Double_t priors[fSelectedSpecies];
 165         if (fEnablePriors) GetPriors(track,priors);
 166         else { for (Int_t i=0;i<fSelectedSpecies;i++) priors[i]=1.;}
 167         ComputeBayesProbabilities(bayesProbabilities,p,priors);
 168         return usedMask;
 169 }
 170
 171
 172 //-------------------------------------------------------------------------------------------------
 173 void AliPIDCombined::GetPriors(const AliVTrack *track, Double_t* p) const {
 174
 175         //
 176         // get priors from distributions
 177         //
 178
 179         Double_t pt=TMath::Abs(track->Pt());
 180         Double_t sumPriors = 0;
 181         for (Int_t i=0;i<fSelectedSpecies;++i){
 182                 if (fPriorsDistributions[i]){
 183                         p[i]=fPriorsDistributions[i]->Interpolate(pt);
 184                 }
 185                 else {
 186                         p[i]=0.;
 187                 }
 188                 sumPriors+=p[i];
 189         }
 190
 191         // normalizing priors
 192         if (sumPriors == 0) return;
 193         for (Int_t i=0;i<fSelectedSpecies;++i){
 194            p[i]=p[i]/sumPriors;
 195         }
 196         return;
 197 }
 198
 199 //-------------------------------------------------------------------------------------------------
 200 void AliPIDCombined::ComputeBayesProbabilities(Double_t* probabilities, const Double_t* probDensity, const Double_t* prior) const {
 201
 202
 203   //
 204   // calculate Bayesian probabilities
 205   //
 206   Double_t sum = 0.;
 207   for (Int_t i = 0; i < fSelectedSpecies; i++) {
 208     sum += probDensity[i] * prior[i];
 209   }
 210   if (sum <= 0) {
 211     AliError("Invalid probability densities or priors");
 212     for (Int_t i = 0; i < fSelectedSpecies; i++) probabilities[i] = -1;
 213     return;
 214   }
 215   for (Int_t i = 0; i < fSelectedSpecies; i++) {
 216     probabilities[i] = probDensity[i] * prior[i] / sum;
 217   }
 218
 219
 220 }
 221
 222 //----------------------------------------------------------------------------------------
 223 void AliPIDCombined::SetCombinedStatus(AliPIDResponse::EDetPidStatus status, UInt_t *maskDetIn, AliPIDResponse::EDetCode bit, Double_t* p) const {
 224   switch (status) {
 225   case AliPIDResponse::kDetNoSignal:
 226     break;
 227   case AliPIDResponse::kDetPidOk:
 228     *maskDetIn = *maskDetIn | bit;
 229     break;
 230   case AliPIDResponse::kDetMismatch:
 231     if ( fRejectMismatchMask & bit) for (Int_t j=0;j<fSelectedSpecies;j++) p[j]=1./fSelectedSpecies;
 232     else *maskDetIn = *maskDetIn | bit;
 233     break;
 234   }
 235
 236 }
 237
 238 ClassImp(AliPIDCombined);
 239