As part of the efforts of the combined PID subgroup we are continuing to develop...
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 17 Jun 2011 18:52:18 +0000 (18:52 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 17 Jun 2011 18:52:18 +0000 (18:52 +0000)
it is possible ESD/AOD "neutral" classes using AliVTrack and pushing detectors to provide
code for the AliPIDResponse class.

This update introduces also a new class AliPIDCombined which uses AliPIDResponse interfaces to obtain
PID probabilities and then computes Bayesian probabilities within the framework discussed in the sub-group.
Note the updated .pkg and STEERBaseLinkDef.h

Pietro Antonioli <Pietro.Antonioli@bo.infn.it>

STEER/AliAODTrack.cxx
STEER/AliAODTrack.h
STEER/AliESDpid.h
STEER/AliESDtrack.h
STEER/AliPIDCombined.cxx [new file with mode: 0644]
STEER/AliPIDCombined.h [new file with mode: 0644]
STEER/AliPIDResponse.cxx
STEER/AliPIDResponse.h
STEER/AliVTrack.h
STEER/CMakelibSTEERBase.pkg

index f2a95e7..47a23c1 100644 (file)
@@ -620,3 +620,47 @@ Float_t AliAODTrack::GetTPCClusterInfo(Int_t nNeighbours/*=3*/, Int_t type/*=0*/
   }
   return 0;  // undefined type - default value
 }
+
+//______________________________________________________________________________
+Double_t  AliAODTrack::GetTRDslice(Int_t plane, Int_t slice) const {
+  //
+  // return TRD Pid information
+  //
+  if (!fDetPid) return -1;
+  Double32_t *trdSlices=fDetPid->GetTRDsignal();
+  if (!trdSlices) return -1;
+  if ((plane<0) || (plane>=kTRDnPlanes)) {
+    return -1.;
+  }
+
+  Int_t ns=fDetPid->GetTRDnSlices();
+  if ((slice<-1) || (slice>=ns)) {
+    return -1.;
+  }
+
+  if(slice>=0) return trdSlices[plane*ns + slice];
+
+  // return average of the dEdx measurements
+  Double_t q=0.; Double32_t *s = &trdSlices[plane*ns];
+  for (Int_t i=0; i<ns; i++, s++) if((*s)>0.) q+=(*s);
+  return q/ns;
+}
+
+//______________________________________________________________________________
+Double_t AliAODTrack::GetTRDmomentum(Int_t plane, Double_t */*sp*/) const
+{
+  //Returns momentum estimation
+  // in TRD layer "plane".
+
+  if (!fDetPid) return -1;
+  Float_t *trdMomentum=fDetPid->GetTRDmomentum();
+
+  if (!trdMomentum) {
+    return -1.;
+  }
+  if ((plane<0) || (plane>=kTRDnPlanes)) {
+    return -1.;
+  }
+
+  return trdMomentum[plane];
+}
index 915af2c..52a3a8e 100644 (file)
@@ -225,7 +225,10 @@ class AliAODTrack : public AliVTrack {
   UShort_t  GetTPCsignalN()      const { return fDetPid?fDetPid->GetTPCsignalN():0;    }
   Double_t  GetTPCmomentum()     const { return fDetPid?fDetPid->GetTPCmomentum():0.;  }
   Double_t  GetTOFsignal()       const { return fDetPid?fDetPid->GetTOFsignal():0.; }
-  
+  void GetIntegratedTimes(Double_t *times) const {if (fDetPid) fDetPid->GetIntegratedTimes(times); }
+  Double_t  GetTRDslice(Int_t plane, Int_t slice) const;
+  Double_t GetTRDmomentum(Int_t plane, Double_t */*sp*/=0x0) const;
+
   
   AliAODPid    *GetDetPid() const { return fDetPid; }
   AliAODVertex *GetProdVertex() const { return (AliAODVertex*)fProdVertex.GetObject(); }
index e1f15d6..5096afc 100644 (file)
@@ -25,7 +25,7 @@ class AliVParticle;
 
 class AliESDpid : public AliPIDResponse  {
 public:
-  AliESDpid(Bool_t forMC=kFALSE): AliPIDResponse(forMC), fRange(5.), fRangeTOFMismatch(5.), fITSPIDmethod(kITSTruncMean) {;}
+  AliESDpid(Bool_t forMC=kFALSE): AliPIDResponse(forMC), fRangeTOFMismatch(5.) {;}
   virtual ~AliESDpid() {}
   
   Int_t MakePID(AliESDEvent *event, Bool_t TPCOnly = kFALSE, Float_t timeZeroTOF=9999) const;
@@ -36,9 +36,6 @@ public:
   //  void MakeHMPIDPID(AliESDtrack *track);
   void MakeTRDPID(AliESDtrack *track) const;
   void CombinePID(AliESDtrack *track) const;
-
-  enum ITSPIDmethod { kITSTruncMean, kITSLikelihood };
-  void SetITSPIDmethod(ITSPIDmethod pmeth) { fITSPIDmethod = pmeth; }
   
   virtual Float_t NumberOfSigmasTOF(const AliVParticle *vtrack, AliPID::EParticleType type, const Float_t timeZeroTOF) const;
   virtual Float_t NumberOfSigmasTOF(const AliVParticle *vtrack, AliPID::EParticleType type) const {return NumberOfSigmasTOF(vtrack,type,0); }
@@ -49,9 +46,7 @@ public:
   Float_t GetNMaxSigmaTOFTPCMismatch() const {return fRangeTOFMismatch;}
 
 private:
-  Float_t           fRange;          // nSigma max in likelihood
   Float_t           fRangeTOFMismatch; // nSigma max for TOF matching with TPC
-  ITSPIDmethod      fITSPIDmethod;   // 0 = trunc mean; 1 = likelihood 
 
   ClassDef(AliESDpid,6)  // PID calculation class
 };
index 0a94c47..835b39f 100644 (file)
@@ -42,29 +42,6 @@ class TPolyMarker3D;
 
 class AliESDtrack : public AliExternalTrackParam {
 public:
-  enum {
-    kITSin=0x0001,kITSout=0x0002,kITSrefit=0x0004,kITSpid=0x0008,
-    kTPCin=0x0010,kTPCout=0x0020,kTPCrefit=0x0040,kTPCpid=0x0080,
-    kTRDin=0x0100,kTRDout=0x0200,kTRDrefit=0x0400,kTRDpid=0x0800,
-    kTOFin=0x1000,kTOFout=0x2000,kTOFrefit=0x4000,kTOFpid=0x8000,
-    kTOFmismatch=0x100000,
-    kHMPIDout=0x10000,kHMPIDpid=0x20000,
-    kEMCALmatch=0x40000,
-    kPHOSmatch=0x200000,
-    kTRDbackup =0x80000,
-    kTRDStop=0x20000000,
-    kESDpid=0x40000000,
-    kTIME=0x80000000,
-    kGlobalMerge=0x08000000,
-    kITSpureSA=0x10000000,
-    kMultInV0 =0x2000000,    //BIT(25): assumed to be belong to V0 in multiplicity estimates
-    kMultSec  =0x4000000,     //BIT(26): assumed to be secondary (due to the DCA) in multiplicity estimates
-    kEmbedded =0x8000000     // BIT(27), 1<<27: Is a track that has been embedded into the event      
-  }; 
-  enum {
-    kTRDnPlanes = 6,
-    kEMCALNoMatch = -4096
-  };
   AliESDtrack();
   AliESDtrack(const AliESDtrack& track);
   AliESDtrack(const AliVTrack* track);
diff --git a/STEER/AliPIDCombined.cxx b/STEER/AliPIDCombined.cxx
new file mode 100644 (file)
index 0000000..1db2c2d
--- /dev/null
@@ -0,0 +1,228 @@
+/**************************************************************************\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
diff --git a/STEER/AliPIDCombined.h b/STEER/AliPIDCombined.h
new file mode 100644 (file)
index 0000000..446449f
--- /dev/null
@@ -0,0 +1,63 @@
+#ifndef ALIPIDCOMBINED_H\r
+#define ALIPIDCOMBINED_H\r
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
+ * See cxx source for full Copyright notice                               */\r
+\r
+//---------------------------------------------------------------//\r
+//        Base class for combining PID response of               //\r
+//        of different detectors                                 //\r
+//        and compute Bayesian probabilities                     //\r
+//                                                               //\r
+//   Origin: Pietro Antonioli, INFN-BO, pietro.antonioli@cern.ch //\r
+//                                                               //\r
+//---------------------------------------------------------------//\r
+\r
+\r
+\r
+#include <TNamed.h>\r
+#include <AliPID.h>\r
+#include <AliPIDResponse.h>\r
+#include <TH1F.h>\r
+\r
+//class TH1;\r
+class AliPIDResponse;\r
+\r
+class AliPIDCombined : public TNamed {\r
+public:\r
+  AliPIDCombined();\r
+  AliPIDCombined(const TString& name, const TString& title);\r
+  virtual ~AliPIDCombined();\r
+\r
+  void SetDetectorMask(Int_t mask) {fDetectorMask=mask;}\r
+  Int_t GetDetectorMask() const {return fDetectorMask;}\r
+  void SetRejectMismatchMask(Int_t mask) {fRejectMismatchMask=mask;}\r
+  Int_t GetRejectMismatchMask() const {return fRejectMismatchMask;}\r
+  void SetEnablePriors(Bool_t flag) {fEnablePriors=flag;}\r
+  Bool_t GetEnablePriors() const {return fEnablePriors;}\r
+  void SetPriorDistribution(AliPID::EParticleType type,TH1F *prior);\r
+  //  const TH1* GetPriorDistribution(AliPID::EParticleType type) const {return (TH1*)fPriorsDistributions[type];}\r
+  TH1* GetPriorDistribution(AliPID::EParticleType type)  const {return (TH1*)fPriorsDistributions[type];}\r
+       \r
+  UInt_t ComputeProbabilities(const AliVTrack *track, const AliPIDResponse *response, Double_t* bayesProbabilities) const;\r
+  void SetSelectedSpecies(Int_t selectedSpecies) {fSelectedSpecies = selectedSpecies;}\r
+  Int_t GetSelectedSpecies() const {return fSelectedSpecies;}\r
+\r
+protected:\r
+  void GetPriors(const AliVTrack *track,Double_t* priors) const;\r
+  void ComputeBayesProbabilities(Double_t* bayesProbabilities,const Double_t* probDensity, const Double_t* priors) const;\r
+  void SetCombinedStatus(const AliPIDResponse::EDetPidStatus status,UInt_t *mask, const AliPIDResponse::EDetCode bit, Double_t* p) const;\r
+\r
+private:\r
+  AliPIDCombined(const AliPIDCombined&);\r
+  AliPIDCombined &operator=(const AliPIDCombined&);\r
+\r
+  Int_t fDetectorMask;       // Detectors included in combined pid\r
+  Int_t fRejectMismatchMask; // Detectors set return flat prob. if mismatch detected \r
+  Bool_t fEnablePriors;      // Enable bayesian PID (if kFALSE priors set flat)\r
+  Int_t fSelectedSpecies;    // Number of selected species to study\r
+  TH1F *fPriorsDistributions[AliPID::kSPECIES+AliPID::kSPECIESLN]; // priors\r
+\r
+  ClassDef(AliPIDCombined,1);\r
+};\r
+\r
+#endif\r
index 89a8862..3e068cd 100644 (file)
@@ -31,6 +31,7 @@
 #include <TFile.h>
 
 #include <AliVEvent.h>
+#include <AliVTrack.h>
 #include <AliLog.h>
 #include <AliPID.h>
 
@@ -44,12 +45,16 @@ fITSResponse(isMC),
 fTPCResponse(),
 fTRDResponse(),
 fTOFResponse(),
+fRange(5.),
+fITSPIDmethod(kITSTruncMean),
 fIsMC(isMC),
 fOADBPath(),
 fBeamType("PP"),
 fLHCperiod(),
 fMCperiodTPC(),
+fMCperiodUser(),
 fRecoPass(0),
+fRecoPassUser(-1),
 fRun(0),
 fOldRun(0),
 fArrPidResponseMaster(0x0),
@@ -81,12 +86,16 @@ fITSResponse(other.fITSResponse),
 fTPCResponse(other.fTPCResponse),
 fTRDResponse(other.fTRDResponse),
 fTOFResponse(other.fTOFResponse),
+fRange(other.fRange),
+fITSPIDmethod(other.fITSPIDmethod),
 fIsMC(other.fIsMC),
 fOADBPath(other.fOADBPath),
 fBeamType("PP"),
 fLHCperiod(),
 fMCperiodTPC(),
+fMCperiodUser(other.fMCperiodUser),
 fRecoPass(0),
+fRecoPassUser(other.fRecoPassUser),
 fRun(0),
 fOldRun(0),
 fArrPidResponseMaster(0x0),
@@ -112,12 +121,16 @@ AliPIDResponse& AliPIDResponse::operator=(const AliPIDResponse &other)
     fTPCResponse=other.fTPCResponse;
     fTRDResponse=other.fTRDResponse;
     fTOFResponse=other.fTOFResponse;
+    fRange=other.fRange;
+    fITSPIDmethod=other.fITSPIDmethod;
     fOADBPath=other.fOADBPath;
     fIsMC=other.fIsMC;
     fBeamType="PP";
     fLHCperiod="";
     fMCperiodTPC="";
+    fMCperiodUser=other.fMCperiodUser;
     fRecoPass=0;
+    fRecoPassUser=other.fRecoPassUser;
     fRun=0;
     fOldRun=0;
     fArrPidResponseMaster=0x0;
@@ -129,6 +142,234 @@ AliPIDResponse& AliPIDResponse::operator=(const AliPIDResponse &other)
 }
 
 //______________________________________________________________________________
+Float_t AliPIDResponse::NumberOfSigmas(EDetCode detCode, const AliVParticle *track, AliPID::EParticleType type) const
+{
+  //
+  // NumberOfSigmas for 'detCode'
+  //
+
+  switch (detCode){
+    case kDetITS: return NumberOfSigmasITS(track, type); break;
+    case kDetTPC: return NumberOfSigmasTPC(track, type); break;
+    case kDetTOF: return NumberOfSigmasTOF(track, type); break;
+//     case kDetTRD: return ComputeTRDProbability(track, type); break;
+//     case kDetPHOS: return ComputePHOSProbability(track, type); break;
+//     case kDetEMCAL: return ComputeEMCALProbability(track, type); break;
+//     case kDetHMPID: return ComputeHMPIDProbability(track, type); break;
+    default: return -999.;
+  }
+
+}
+
+//______________________________________________________________________________
+AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability  (EDetCode detCode,  const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
+{
+  //
+  // Compute PID response of 'detCode'
+  //
+
+  switch (detCode){
+    case kDetITS: return ComputeITSProbability(track, nSpecies, p); break;
+    case kDetTPC: return ComputeTPCProbability(track, nSpecies, p); break;
+    case kDetTOF: return ComputeTOFProbability(track, nSpecies, p); break;
+    case kDetTRD: return ComputeTRDProbability(track, nSpecies, p); break;
+    case kDetPHOS: return ComputePHOSProbability(track, nSpecies, p); break;
+    case kDetEMCAL: return ComputeEMCALProbability(track, nSpecies, p); break;
+    case kDetHMPID: return ComputeHMPIDProbability(track, nSpecies, p); break;
+    default: return kDetNoSignal;
+  }
+}
+
+//______________________________________________________________________________
+AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeITSProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
+{
+  //
+  // Compute PID response for the ITS
+  //
+
+  // set flat distribution (no decision)
+  for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
+
+  if ((track->GetStatus()&AliVTrack::kITSin)==0 &&
+    (track->GetStatus()&AliVTrack::kITSout)==0) return kDetNoSignal;
+  
+  Double_t mom=track->P();
+  Double_t dedx=track->GetITSsignal();
+  Bool_t isSA=kTRUE;
+  Double_t momITS=mom;
+  ULong_t trStatus=track->GetStatus();
+  if(trStatus&AliVTrack::kTPCin) isSA=kFALSE;
+  UChar_t clumap=track->GetITSClusterMap();
+  Int_t nPointsForPid=0;
+  for(Int_t i=2; i<6; i++){
+    if(clumap&(1<<i)) ++nPointsForPid;
+  }
+  
+  if(nPointsForPid<3) { // track not to be used for combined PID purposes
+    //       track->ResetStatus(AliVTrack::kITSpid);
+    return kDetNoSignal;
+  }
+
+  Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
+  for (Int_t j=0; j<AliPID::kSPECIES; j++) {
+    Double_t mass=AliPID::ParticleMass(j);//GeV/c^2
+    Double_t bethe=fITSResponse.Bethe(momITS,mass);
+    Double_t sigma=fITSResponse.GetResolution(bethe,nPointsForPid,isSA);
+    if (TMath::Abs(dedx-bethe) > fRange*sigma) {
+      p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
+    } else {
+      p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
+      mismatch=kFALSE;
+    }
+
+    // Check for particles heavier than (AliPID::kSPECIES - 1)
+    //       if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
+
+  }
+
+  if (mismatch){
+    for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
+    return kDetNoSignal;
+  }
+
+    
+  return kDetPidOk;
+}
+//______________________________________________________________________________
+AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTPCProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
+{
+  //
+  // Compute PID response for the TPC
+  //
+
+  // set flat distribution (no decision)
+  for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
+
+  // check quality of the track
+  if ( (track->GetStatus()&AliVTrack::kTPCin )==0 && (track->GetStatus()&AliVTrack::kTPCout)==0 ) return kDetNoSignal;
+
+  Double_t mom = track->GetTPCmomentum();
+
+  Double_t dedx=track->GetTPCsignal();
+  Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
+
+  for (Int_t j=0; j<AliPID::kSPECIES; j++) {
+    AliPID::EParticleType type=AliPID::EParticleType(j);
+    Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type);
+    Double_t sigma=fTPCResponse.GetExpectedSigma(mom,track->GetTPCsignalN(),type);
+    if (TMath::Abs(dedx-bethe) > fRange*sigma) {
+      p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
+    } else {
+      p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
+      mismatch=kFALSE;
+    }
+
+    // TODO: Light nuclei, also in TPC pid response
+    
+    // Check for particles heavier than (AliPID::kSPECIES - 1)
+//     if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
+
+  }
+
+  if (mismatch){
+    for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
+    return kDetNoSignal;
+  }
+
+  return kDetPidOk;
+}
+//______________________________________________________________________________
+AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTOFProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
+{
+  //
+  // Compute PID response for the
+  //
+
+  // set flat distribution (no decision)
+  for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
+  
+  if ((track->GetStatus()&AliVTrack::kTOFout)==0) return kDetNoSignal;
+  if ((track->GetStatus()&AliVTrack::kTIME)==0) return kDetNoSignal;
+  
+  Double_t time[AliPID::kSPECIESN];
+  track->GetIntegratedTimes(time);
+  
+  Double_t sigma[AliPID::kSPECIES];
+  for (Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++) {
+    sigma[iPart] = fTOFResponse.GetExpectedSigma(track->P(),time[iPart],AliPID::ParticleMass(iPart));
+  }
+  
+  Bool_t mismatch = kTRUE/*, heavy = kTRUE*/;
+  for (Int_t j=0; j<AliPID::kSPECIES; j++) {
+    AliPID::EParticleType type=AliPID::EParticleType(j);
+    Double_t nsigmas=NumberOfSigmasTOF(track,type);
+    
+    Double_t sig = sigma[j];
+    if (TMath::Abs(nsigmas) > (fRange+2)) {
+      p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
+    } else
+      p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
+
+    if (TMath::Abs(nsigmas)<5.){
+      Double_t nsigmasTPC=NumberOfSigmasTPC(track,type);
+      if (TMath::Abs(nsigmasTPC)<5.) mismatch=kFALSE;
+    }
+  }
+
+  if (mismatch){
+    return kDetMismatch;    
+  }
+
+    // TODO: Light nuclei
+    
+  return kDetPidOk;
+}
+//______________________________________________________________________________
+AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTRDProbability  (const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const
+{
+  //
+  // Compute PID response for the
+  //
+
+  // set flat distribution (no decision)
+  for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
+  return kDetNoSignal;
+}
+//______________________________________________________________________________
+AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeEMCALProbability(const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const
+{
+  //
+  // Compute PID response for the EMCAL
+  //
+
+  // set flat distribution (no decision)
+  for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
+  return kDetNoSignal;
+}
+//______________________________________________________________________________
+AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePHOSProbability (const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const
+{
+  //
+  // Compute PID response for the PHOS
+  //
+
+  // set flat distribution (no decision)
+  for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
+  return kDetNoSignal;
+}
+//______________________________________________________________________________
+AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeHMPIDProbability(const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const
+{
+  //
+  // Compute PID response for the HMPID
+  //
+
+  // set flat distribution (no decision)
+  for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
+  return kDetNoSignal;
+}
+
+//______________________________________________________________________________
 void AliPIDResponse::InitialiseEvent(AliVEvent *event, Int_t pass)
 {
   //
index 7d89dee..c769da0 100644 (file)
@@ -29,22 +29,58 @@ public:
   AliPIDResponse(Bool_t isMC=kFALSE);
   virtual ~AliPIDResponse();
 
+  enum EDetCode {
+    kDetITS = 0x1,
+    kDetTPC = 0x2,
+    kDetTRD = 0x4,
+    kDetTOF = 0x8,
+    kDetHMPID = 0x10,
+    kDetEMCAL = 0x20,
+    kDetPHOS = 0x40
+  };
+
   enum EStartTimeType_t {kFILL_T0,kTOF_T0, kT0_T0, kBest_T0};
+
+  enum ITSPIDmethod { kITSTruncMean, kITSLikelihood };
+
+  enum EDetPidStatus {
+    kDetNoSignal=0,
+    kDetPidOk=1,
+    kDetMismatch=2
+  };
   
   AliITSPIDResponse &GetITSResponse() {return fITSResponse;}
   AliTPCPIDResponse &GetTPCResponse() {return fTPCResponse;}
   AliTOFPIDResponse &GetTOFResponse() {return fTOFResponse;}
   AliTRDPIDResponse &GetTRDResponse() {return fTRDResponse;}
+
+  Float_t NumberOfSigmas(EDetCode detCode, const AliVParticle *track, AliPID::EParticleType type) const;
   
   virtual Float_t NumberOfSigmasITS(const AliVParticle *track, AliPID::EParticleType type) const;
   virtual Float_t NumberOfSigmasTPC(const AliVParticle *track, AliPID::EParticleType type) const;
   virtual Float_t NumberOfSigmasTOF(const AliVParticle *track, AliPID::EParticleType type) const = 0;
 
+  EDetPidStatus ComputePIDProbability  (EDetCode detCode, const AliVTrack *track, Int_t nSpecies, Double_t p[]) const;
+  
+  EDetPidStatus ComputeITSProbability  (const AliVTrack *track, Int_t nSpecies, Double_t  p[]) const;
+  EDetPidStatus ComputeTPCProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const;
+  EDetPidStatus ComputeTOFProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const;
+  EDetPidStatus ComputeTRDProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const;
+  EDetPidStatus ComputeEMCALProbability(const AliVTrack *track, Int_t nSpecies, Double_t p[]) const;
+  EDetPidStatus ComputePHOSProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const;
+  EDetPidStatus ComputeHMPIDProbability(const AliVTrack *track, Int_t nSpecies, Double_t p[]) const;
+
+
+  void SetITSPIDmethod(ITSPIDmethod pmeth) { fITSPIDmethod = pmeth; }
   virtual void SetTOFResponse(AliVEvent */*event*/,EStartTimeType_t /*option*/) {;}
 
   void SetOADBPath(const char* path) {fOADBPath=path;}
   void InitialiseEvent(AliVEvent *event, Int_t pass);
 
+  // User settings for the MC period and reco pass
+  void SetMCperiod(const char *mcPeriod) {fMCperiodUser=mcPeriod;}
+  void SetRecoPass(Int_t recoPass)       {fRecoPassUser=recoPass;}
+  
   AliPIDResponse(const AliPIDResponse &other);
   AliPIDResponse& operator=(const AliPIDResponse &other);
   
@@ -54,6 +90,9 @@ protected:
   AliTRDPIDResponse fTRDResponse;    //PID response function of the TRD
   AliTOFPIDResponse fTOFResponse;    //PID response function of the TOF
 
+  Float_t           fRange;          // nSigma max in likelihood
+  ITSPIDmethod      fITSPIDmethod;   // 0 = trunc mean; 1 = likelihood
+
 private:
   Bool_t fIsMC;                        //  If we run on MC data
 
@@ -62,10 +101,12 @@ private:
   TString fBeamType;                   //! beam type (PP) or (PBPB)
   TString fLHCperiod;                  //! LHC period
   TString fMCperiodTPC;                //! corresponding MC period to use for the TPC splines
+  TString fMCperiodUser;               //  MC prodution requested by the user
   Int_t   fRecoPass;                   //! reconstruction pass
+  Int_t   fRecoPassUser;               //  reconstruction pass explicitly set by the user
   Int_t   fRun;                        //! current run number
   Int_t   fOldRun;                     //! current run number
-
+  
   TObjArray *fArrPidResponseMaster;    //!  TPC pid splines
   TF1       *fResolutionCorrection;    //! TPC resolution correction
 
index faba22c..80636f6 100644 (file)
@@ -19,6 +19,30 @@ class AliExternalTrackParam;
 class AliVTrack: public AliVParticle {
 
 public:
+  enum {
+    kITSin=0x0001,kITSout=0x0002,kITSrefit=0x0004,kITSpid=0x0008,
+    kTPCin=0x0010,kTPCout=0x0020,kTPCrefit=0x0040,kTPCpid=0x0080,
+    kTRDin=0x0100,kTRDout=0x0200,kTRDrefit=0x0400,kTRDpid=0x0800,
+    kTOFin=0x1000,kTOFout=0x2000,kTOFrefit=0x4000,kTOFpid=0x8000,
+    kTOFmismatch=0x100000,
+    kHMPIDout=0x10000,kHMPIDpid=0x20000,
+    kEMCALmatch=0x40000,
+    kPHOSmatch=0x200000,
+    kTRDbackup =0x80000,
+    kTRDStop=0x20000000,
+    kESDpid=0x40000000,
+    kTIME=0x80000000,
+    kGlobalMerge=0x08000000,
+    kITSpureSA=0x10000000,
+    kMultInV0 =0x2000000,    //BIT(25): assumed to be belong to V0 in multiplicity estimates
+    kMultSec  =0x4000000,     //BIT(26): assumed to be secondary (due to the DCA) in multiplicity estimates
+    kEmbedded =0x8000000     // BIT(27), 1<<27: Is a track that has been embedded into the event
+  };
+  enum {
+    kTRDnPlanes = 6,
+    kEMCALNoMatch = -4096
+  };
+
   AliVTrack() { }
   virtual ~AliVTrack() { }
   AliVTrack(const AliVTrack& vTrack); 
@@ -29,6 +53,7 @@ public:
   virtual Float_t  GetTPCClusterInfo(Int_t /*nNeighbours*/, Int_t /*type*/, Int_t /*row0*/=0, Int_t /*row1*/=159) const {return 0.;}
   virtual UShort_t GetTPCNcls() const { return 0;}
   virtual UShort_t GetTPCNclsF() const { return 0;}
+  virtual Double_t GetTRDslice(Int_t /*plane*/, Int_t /*slice*/) const { return -1.; }
   
   //pid info
   virtual Double_t  GetITSsignal()       const {return 0.;}
@@ -36,6 +61,8 @@ public:
   virtual UShort_t  GetTPCsignalN()      const {return 0 ;}
   virtual Double_t  GetTPCmomentum()     const {return 0.;}
   virtual Double_t  GetTOFsignal()       const {return 0.;}
+  virtual void GetIntegratedTimes(Double_t */*times*/) const { return; }
+  virtual Double_t GetTRDmomentum(Int_t /*plane*/, Double_t */*sp*/=0x0) const {return 0.;}
   
   virtual ULong_t  GetStatus() const = 0;
   virtual Bool_t   GetXYZ(Double_t *p) const = 0;
index 79d9a8e..ac04108 100644 (file)
@@ -65,6 +65,7 @@ set ( SRCS
     AliTPCPIDResponse.cxx 
     AliTOFPIDResponse.cxx 
     AliTRDPIDResponse.cxx 
+    AliPIDCombined.cxx
     AliDAQ.cxx 
     AliRefArray.cxx
     )