]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
Bug fixed (out-of-bounds for light nuclei) in SetPriorDistribution
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 20 Jun 2011 14:00:38 +0000 (14:00 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 20 Jun 2011 14:00:38 +0000 (14:00 +0000)
P. Antonioli <Pietro.Antonioli@bo.infn.it>

STEER/AliPIDCombined.cxx

index 1db2c2dea80885ef9e36a0857bfcad430c92457f..fc1218645d4c558824dbbb48e4fca4bde32bb35e 100644 (file)
-/**************************************************************************\r
- * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
- *                                                                        *\r
- * Author: The ALICE Off-line Project.                                    *\r
- * Contributors are mentioned in the code where appropriate.              *\r
- *                                                                        *\r
- * Permission to use, copy, modify and distribute this software and its   *\r
- * documentation strictly for non-commercial purposes is hereby granted   *\r
- * without fee, provided that the above copyright notice appears in all   *\r
- * copies and that both the copyright notice and this permission notice   *\r
- * appear in the supporting documentation. The authors make no claims     *\r
- * about the suitability of this software for any purpose. It is          *\r
- * provided "as is" without express or implied warranty.                  *\r
- **************************************************************************/\r
-\r
-\r
-//-----------------------------------------------------------------\r
-//        Base class for combining PID of different detectors    //\r
-//        (user selected) and compute Bayesian probabilities     //\r
-//                                                               //\r
-//                                                               //\r
-//   Origin: Pietro Antonioli, INFN-BO Pietro.Antonioli@cern.ch  //\r
-//                                                               //\r
-//-----------------------------------------------------------------\r
-\r
-#include <TH1.h>\r
-\r
-#include <AliVTrack.h>\r
-#include <AliLog.h>\r
-#include <AliPID.h>\r
-#include <AliPIDResponse.h>\r
-\r
-#include "AliPIDCombined.h"\r
-\r
-AliPIDCombined::AliPIDCombined():\r
-       TNamed("CombinedPID","CombinedPID"),\r
-       fDetectorMask(0),\r
-       fRejectMismatchMask(0x7F),\r
-       fEnablePriors(kTRUE),\r
-       fSelectedSpecies(AliPID::kSPECIES)\r
-{      \r
-  //\r
-  // default constructor\r
-  //   \r
-  for (Int_t i=0;i<AliPID::kSPECIES+AliPID::kSPECIESLN;i++) fPriorsDistributions[i]=NULL;\r
-  AliLog::SetClassDebugLevel("AliPIDCombined",10);\r
-}\r
-\r
-//-------------------------------------------------------------------------------------------------    \r
-AliPIDCombined::AliPIDCombined(const TString& name,const TString& title):\r
-       TNamed(name,title),\r
-       fDetectorMask(0),\r
-       fRejectMismatchMask(0x7F),\r
-       fEnablePriors(kTRUE),\r
-       fSelectedSpecies(AliPID::kSPECIES)\r
-{\r
-  //\r
-  // default constructor with name and title\r
-  //\r
-  for (Int_t i=0;i<AliPID::kSPECIES+AliPID::kSPECIESLN;i++) fPriorsDistributions[i]=NULL;\r
-  AliLog::SetClassDebugLevel("AliPIDCombined",10);\r
-\r
-}\r
-\r
-//-------------------------------------------------------------------------------------------------    \r
-AliPIDCombined::~AliPIDCombined() {\r
-\r
-  for(Int_t i=0;i < AliPID::kSPECIES+AliPID::kSPECIESLN;i++){\r
-    if(fPriorsDistributions[i])\r
-      delete fPriorsDistributions[i];\r
-  }\r
-}\r
-\r
-//-------------------------------------------------------------------------------------------------    \r
-void AliPIDCombined::SetPriorDistribution(AliPID::EParticleType type,TH1F *prior) {\r
-  if ( (type < 0) || (type>AliPID::kSPECIES+AliPID::kSPECIESLN) ) {\r
-    AliError(Form("Invalid EParticleType setting prior (type: %d)",type));\r
-    return;\r
-  }\r
-  if(prior) {\r
-    if (fPriorsDistributions[type]) {\r
-      delete fPriorsDistributions[type]; \r
-    }\r
-    fPriorsDistributions[type]=new TH1F(*prior);\r
-  }\r
-}\r
-\r
-//-------------------------------------------------------------------------------------------------    \r
-UInt_t AliPIDCombined::ComputeProbabilities(const AliVTrack *track, const AliPIDResponse *response, Double_t* bayesProbabilities) const {\r
-       Double_t pITS[fSelectedSpecies];\r
-       Double_t pTPC[fSelectedSpecies];\r
-       Double_t pTRD[fSelectedSpecies];\r
-       Double_t pTOF[fSelectedSpecies];\r
-       Double_t pHMPID[fSelectedSpecies];\r
-       Double_t pEMCAL[fSelectedSpecies];\r
-       Double_t pPHOS[fSelectedSpecies];\r
-       for (Int_t i=0;i<fSelectedSpecies;i++) {\r
-        pITS[i]=1.;\r
-        pTPC[i]=1.;\r
-        pTRD[i]=1.;\r
-        pTOF[i]=1.;\r
-        pHMPID[i]=1.;\r
-        pEMCAL[i]=1.;\r
-        pPHOS[i]=1.;\r
-       }\r
-       UInt_t usedMask=0;\r
-       AliPIDResponse::EDetPidStatus status=AliPIDResponse::kDetNoSignal;\r
-       Double_t p[fSelectedSpecies];\r
-\r
-       // getting probability distributions for selected detectors only\r
-       if (fDetectorMask & AliPIDResponse::kDetITS) {\r
-         status = response->ComputeITSProbability(track, fSelectedSpecies, pITS);\r
-         SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetITS,pITS);\r
-       }\r
-\r
-       if (fDetectorMask & AliPIDResponse::kDetTPC) { \r
-         status = response->ComputeTPCProbability(track, fSelectedSpecies, pTPC);\r
-         SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetTPC,pTPC);\r
-       }\r
-\r
-\r
-       if (fDetectorMask & AliPIDResponse::kDetTRD) { \r
-         status = response->ComputeTRDProbability(track, fSelectedSpecies, pTRD);\r
-         SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetTRD,pTRD);\r
-       }\r
-\r
-       if (fDetectorMask &  AliPIDResponse::kDetTOF) { \r
-         status = response->ComputeTOFProbability(track, fSelectedSpecies, pTOF);\r
-         SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetTOF,pTOF);\r
-       }\r
-\r
-       if (fDetectorMask & AliPIDResponse::kDetHMPID) { \r
-         status = response->ComputeHMPIDProbability(track, fSelectedSpecies, pHMPID);\r
-         SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetHMPID,pHMPID);\r
-       }\r
-\r
-\r
-       if (fDetectorMask & AliPIDResponse::kDetEMCAL) { \r
-         status = response->ComputeEMCALProbability(track, fSelectedSpecies, pEMCAL);\r
-         SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetEMCAL,pEMCAL);\r
-       }\r
-\r
-\r
-       if (fDetectorMask & AliPIDResponse::kDetPHOS) { \r
-         status = response->ComputePHOSProbability(track, fSelectedSpecies, pPHOS);\r
-         SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetPHOS,pPHOS);\r
-       }\r
-\r
-\r
-       for (Int_t i =0; i<fSelectedSpecies; i++){\r
-         p[i] = pITS[i]*pTPC[i]*pTRD[i]*pTOF[i]*pHMPID[i]*pEMCAL[i]*pPHOS[i];\r
-       }\r
-       Double_t priors[fSelectedSpecies];\r
-       if (fEnablePriors) GetPriors(track,priors);\r
-       else { for (Int_t i=0;i<fSelectedSpecies;i++) priors[i]=1.;}\r
-       ComputeBayesProbabilities(bayesProbabilities,p,priors);\r
-       return usedMask;\r
-}\r
-\r
-\r
-//-------------------------------------------------------------------------------------------------\r
-void AliPIDCombined::GetPriors(const AliVTrack *track, Double_t* p) const {\r
-       \r
-       //\r
-       // get priors from distributions\r
-       //\r
-       \r
-       Double_t pt=TMath::Abs(track->Pt());\r
-       Double_t sumPriors = 0;\r
-       for (Int_t i=0;i<fSelectedSpecies;++i){\r
-               if (fPriorsDistributions[i]){\r
-                       p[i]=fPriorsDistributions[i]->Interpolate(pt);\r
-               }\r
-               else {\r
-                       p[i]=0.;\r
-               }\r
-               sumPriors+=p[i];                \r
-       }\r
-\r
-       // normalizing priors\r
-       if (sumPriors == 0) return;\r
-       for (Int_t i=0;i<AliPID::kSPECIES;++i){\r
-          p[i]=p[i]/sumPriors;\r
-       }\r
-       return;\r
-}\r
-\r
-//-------------------------------------------------------------------------------------------------    \r
-void AliPIDCombined::ComputeBayesProbabilities(Double_t* probabilities, const Double_t* probDensity, const Double_t* prior) const {\r
-\r
-\r
-  //\r
-  // calculate Bayesian probabilities\r
-  //\r
-  Double_t sum = 0.;\r
-  for (Int_t i = 0; i < fSelectedSpecies; i++) {\r
-    sum += probDensity[i] * prior[i];\r
-  }\r
-  if (sum <= 0) {\r
-    AliError("Invalid probability densities or priors");\r
-    for (Int_t i = 0; i < fSelectedSpecies; i++) probabilities[i] = -1;\r
-    return;\r
-  }\r
-  for (Int_t i = 0; i < fSelectedSpecies; i++) {\r
-    probabilities[i] = probDensity[i] * prior[i] / sum;\r
-  }\r
-\r
-\r
-}\r
-\r
-//----------------------------------------------------------------------------------------\r
-void AliPIDCombined::SetCombinedStatus(AliPIDResponse::EDetPidStatus status, UInt_t *maskDetIn, AliPIDResponse::EDetCode bit, Double_t* p) const {\r
-  switch (status) {\r
-  case AliPIDResponse::kDetNoSignal:\r
-    break;\r
-  case AliPIDResponse::kDetPidOk:\r
-    *maskDetIn = *maskDetIn | bit;\r
-    break;\r
-  case AliPIDResponse::kDetMismatch:\r
-    if ( fRejectMismatchMask & bit) for (Int_t j=0;j<fSelectedSpecies;j++) p[j]=1./fSelectedSpecies;\r
-    else *maskDetIn = *maskDetIn | bit;\r
-    break;\r
-  }\r
-\r
-}\r
-\r
-ClassImp(AliPIDCombined);\r
-\r
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+
+//-----------------------------------------------------------------
+//        Base class for combining PID of different detectors    //
+//        (user selected) and compute Bayesian probabilities     //
+//                                                               //
+//                                                               //
+//   Origin: Pietro Antonioli, INFN-BO Pietro.Antonioli@cern.ch  //
+//                                                               //
+//-----------------------------------------------------------------
+
+#include <TH1.h>
+
+#include <AliVTrack.h>
+#include <AliLog.h>
+#include <AliPID.h>
+#include <AliPIDResponse.h>
+
+#include "AliPIDCombined.h"
+
+AliPIDCombined::AliPIDCombined():
+       TNamed("CombinedPID","CombinedPID"),
+       fDetectorMask(0),
+       fRejectMismatchMask(0x7F),
+       fEnablePriors(kTRUE),
+       fSelectedSpecies(AliPID::kSPECIES)
+{      
+  //
+  // default constructor
+  //   
+  for (Int_t i=0;i<AliPID::kSPECIES+AliPID::kSPECIESLN;i++) fPriorsDistributions[i]=NULL;
+  AliLog::SetClassDebugLevel("AliPIDCombined",10);
+}
+
+//-------------------------------------------------------------------------------------------------    
+AliPIDCombined::AliPIDCombined(const TString& name,const TString& title):
+       TNamed(name,title),
+       fDetectorMask(0),
+       fRejectMismatchMask(0x7F),
+       fEnablePriors(kTRUE),
+       fSelectedSpecies(AliPID::kSPECIES)
+{
+  //
+  // default constructor with name and title
+  //
+  for (Int_t i=0;i<AliPID::kSPECIES+AliPID::kSPECIESLN;i++) fPriorsDistributions[i]=NULL;
+  AliLog::SetClassDebugLevel("AliPIDCombined",10);
+
+}
+
+//-------------------------------------------------------------------------------------------------    
+AliPIDCombined::~AliPIDCombined() {
+
+  for(Int_t i=0;i < AliPID::kSPECIES+AliPID::kSPECIESLN;i++){
+    if(fPriorsDistributions[i])
+      delete fPriorsDistributions[i];
+  }
+}
+
+//-------------------------------------------------------------------------------------------------    
+void AliPIDCombined::SetPriorDistribution(AliPID::EParticleType type,TH1F *prior) {
+  if ( (type < 0) || ( type >= (AliPID::kSPECIESN+AliPID::kSPECIESLN) ) ){
+    AliError(Form("Invalid EParticleType setting prior (offending type: %d)",type));
+    return;
+  }
+  if(prior) {
+    Int_t i = type;
+    if (type >= (AliPID::EParticleType)AliPID::kSPECIES ) {
+      if (type < AliPID::kDeuteron) {
+       AliError(Form("Invalid EParticleType setting prior. Type: %d (neutral) not supported",type));
+       return;
+      } else i = (Int_t)type - (AliPID::kSPECIESN-AliPID::kSPECIES);
+    }
+    if (fPriorsDistributions[i]) {
+      delete fPriorsDistributions[i]; 
+    }
+    fPriorsDistributions[i]=new TH1F(*prior);
+  }
+}
+
+//-------------------------------------------------------------------------------------------------    
+UInt_t AliPIDCombined::ComputeProbabilities(const AliVTrack *track, const AliPIDResponse *response, Double_t* bayesProbabilities) const {
+       Double_t pITS[fSelectedSpecies];
+       Double_t pTPC[fSelectedSpecies];
+       Double_t pTRD[fSelectedSpecies];
+       Double_t pTOF[fSelectedSpecies];
+       Double_t pHMPID[fSelectedSpecies];
+       Double_t pEMCAL[fSelectedSpecies];
+       Double_t pPHOS[fSelectedSpecies];
+       for (Int_t i=0;i<fSelectedSpecies;i++) {
+        pITS[i]=1.;
+        pTPC[i]=1.;
+        pTRD[i]=1.;
+        pTOF[i]=1.;
+        pHMPID[i]=1.;
+        pEMCAL[i]=1.;
+        pPHOS[i]=1.;
+       }
+       UInt_t usedMask=0;
+       AliPIDResponse::EDetPidStatus status=AliPIDResponse::kDetNoSignal;
+       Double_t p[fSelectedSpecies];
+
+       // getting probability distributions for selected detectors only
+       if (fDetectorMask & AliPIDResponse::kDetITS) {
+         status = response->ComputeITSProbability(track, fSelectedSpecies, pITS);
+         SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetITS,pITS);
+       }
+
+       if (fDetectorMask & AliPIDResponse::kDetTPC) { 
+         status = response->ComputeTPCProbability(track, fSelectedSpecies, pTPC);
+         SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetTPC,pTPC);
+       }
+
+
+       if (fDetectorMask & AliPIDResponse::kDetTRD) { 
+         status = response->ComputeTRDProbability(track, fSelectedSpecies, pTRD);
+         SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetTRD,pTRD);
+       }
+
+       if (fDetectorMask &  AliPIDResponse::kDetTOF) { 
+         status = response->ComputeTOFProbability(track, fSelectedSpecies, pTOF);
+         SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetTOF,pTOF);
+       }
+
+       if (fDetectorMask & AliPIDResponse::kDetHMPID) { 
+         status = response->ComputeHMPIDProbability(track, fSelectedSpecies, pHMPID);
+         SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetHMPID,pHMPID);
+       }
+
+
+       if (fDetectorMask & AliPIDResponse::kDetEMCAL) { 
+         status = response->ComputeEMCALProbability(track, fSelectedSpecies, pEMCAL);
+         SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetEMCAL,pEMCAL);
+       }
+
+
+       if (fDetectorMask & AliPIDResponse::kDetPHOS) { 
+         status = response->ComputePHOSProbability(track, fSelectedSpecies, pPHOS);
+         SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetPHOS,pPHOS);
+       }
+
+
+       for (Int_t i =0; i<fSelectedSpecies; i++){
+         p[i] = pITS[i]*pTPC[i]*pTRD[i]*pTOF[i]*pHMPID[i]*pEMCAL[i]*pPHOS[i];
+       }
+       Double_t priors[fSelectedSpecies];
+       if (fEnablePriors) GetPriors(track,priors);
+       else { for (Int_t i=0;i<fSelectedSpecies;i++) priors[i]=1.;}
+       ComputeBayesProbabilities(bayesProbabilities,p,priors);
+       return usedMask;
+}
+
+
+//-------------------------------------------------------------------------------------------------
+void AliPIDCombined::GetPriors(const AliVTrack *track, Double_t* p) const {
+       
+       //
+       // get priors from distributions
+       //
+       
+       Double_t pt=TMath::Abs(track->Pt());
+       Double_t sumPriors = 0;
+       for (Int_t i=0;i<fSelectedSpecies;++i){
+               if (fPriorsDistributions[i]){
+                       p[i]=fPriorsDistributions[i]->Interpolate(pt);
+               }
+               else {
+                       p[i]=0.;
+               }
+               sumPriors+=p[i];                
+       }
+
+       // normalizing priors
+       if (sumPriors == 0) return;
+       for (Int_t i=0;i<fSelectedSpecies;++i){
+          p[i]=p[i]/sumPriors;
+       }
+       return;
+}
+
+//-------------------------------------------------------------------------------------------------    
+void AliPIDCombined::ComputeBayesProbabilities(Double_t* probabilities, const Double_t* probDensity, const Double_t* prior) const {
+
+
+  //
+  // calculate Bayesian probabilities
+  //
+  Double_t sum = 0.;
+  for (Int_t i = 0; i < fSelectedSpecies; i++) {
+    sum += probDensity[i] * prior[i];
+  }
+  if (sum <= 0) {
+    AliError("Invalid probability densities or priors");
+    for (Int_t i = 0; i < fSelectedSpecies; i++) probabilities[i] = -1;
+    return;
+  }
+  for (Int_t i = 0; i < fSelectedSpecies; i++) {
+    probabilities[i] = probDensity[i] * prior[i] / sum;
+  }
+
+
+}
+
+//----------------------------------------------------------------------------------------
+void AliPIDCombined::SetCombinedStatus(AliPIDResponse::EDetPidStatus status, UInt_t *maskDetIn, AliPIDResponse::EDetCode bit, Double_t* p) const {
+  switch (status) {
+  case AliPIDResponse::kDetNoSignal:
+    break;
+  case AliPIDResponse::kDetPidOk:
+    *maskDetIn = *maskDetIn | bit;
+    break;
+  case AliPIDResponse::kDetMismatch:
+    if ( fRejectMismatchMask & bit) for (Int_t j=0;j<fSelectedSpecies;j++) p[j]=1./fSelectedSpecies;
+    else *maskDetIn = *maskDetIn | bit;
+    break;
+  }
+
+}
+
+ClassImp(AliPIDCombined);
+