//-----------------------------------------------------------------
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 27 Apr 2010 09:52:48 +0000 (09:52 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 27 Apr 2010 09:52:48 +0000 (09:52 +0000)
//           Implementation of the combined PID class
//           For the AOD Class
//           containing information on the particle identification
//      Origin: Rosa Romita, GSI, r.romita@gsi.de
//-----------------------------------------------------------------
Marco van Leeuwen

STEER/AODLinkDef.h
STEER/AliAODpidUtil.cxx [new file with mode: 0644]
STEER/AliAODpidUtil.h [new file with mode: 0644]
STEER/libAOD.pkg

index 3f74305..b37ae35 100644 (file)
@@ -44,6 +44,7 @@
 #pragma link C++ class AliAODPWG4Particle+;
 #pragma link C++ class AliAODPWG4ParticleCorrelation+;
 #pragma link C++ class AliAODDimuon+;
+#pragma link C++ class AliAODpidUtil+;
 
 
 #endif
diff --git a/STEER/AliAODpidUtil.cxx b/STEER/AliAODpidUtil.cxx
new file mode 100644 (file)
index 0000000..3d57da6
--- /dev/null
@@ -0,0 +1,203 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+/* $Id: AliAODpidUtil.cxx 38329 2010-01-17 19:17:24Z hristov $ */
+
+//-----------------------------------------------------------------
+//           Implementation of the combined PID class
+//           For the AOD Class
+//           containing information on the particle identification
+//      Origin: Rosa Romita, GSI, r.romita@gsi.de
+//-----------------------------------------------------------------
+
+#include "AliLog.h"
+#include "AliPID.h"
+#include "AliAODpidUtil.h"
+#include "AliAODEvent.h"
+#include "AliAODTrack.h"
+#include "AliAODPid.h"
+#include "AliTRDPIDResponse.h"
+
+ClassImp(AliAODpidUtil)
+
+  Int_t AliAODpidUtil::MakePID(AliAODTrack *track,Float_t TimeZeroTOF,Double_t *p) const {
+  //
+  //  Calculate probabilities for all detectors, except if TPConly==kTRUE
+  //  and combine PID
+  //  
+  //   Option TPConly==kTRUE is used during reconstruction, 
+  //  because ITS tracking uses TPC pid
+  //  HMPID and TRD pid are done in detector reconstructors
+  //
+
+  /*
+    Float_t TimeZeroTOF = 0;
+    if (subtractT0) 
+    TimeZeroTOF = event->GetT0();
+  */
+  Int_t ns=AliPID::kSPECIES;
+  Double_t tpcPid[AliPID::kSPECIES];
+  MakeTPCPID(track,tpcPid);
+  Double_t itsPid[AliPID::kSPECIES];
+  Double_t tofPid[AliPID::kSPECIES];
+  Double_t trdPid[AliPID::kSPECIES];
+  MakeITSPID(track,itsPid);
+  MakeTOFPID(track, TimeZeroTOF,tofPid);
+  //MakeHMPIDPID(track);
+  MakeTRDPID(track,trdPid);
+  for (Int_t j=0; j<ns; j++) {
+    p[j]=tpcPid[j]*itsPid[j]*tofPid[j]*trdPid[j];
+  }
+
+  return 0;
+}
+//_________________________________________________________________________
+void AliAODpidUtil::MakeTPCPID(AliAODTrack *track,Double_t *p) const
+{
+  //
+  //  TPC pid using bethe-bloch and gaussian response
+  //
+
+  Double_t mom = track->P();
+  AliAODPid *pidObj = track->GetDetPid();
+  if (pidObj) mom = pidObj->GetTPCmomentum();
+   
+  Double_t dedx=pidObj->GetTPCsignal(); 
+  Bool_t mismatch=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,50,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;
+    }
+
+  }
+
+  if (mismatch)
+    for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1/AliPID::kSPECIES;
+
+
+  return;
+}
+//_________________________________________________________________________
+void AliAODpidUtil::MakeITSPID(AliAODTrack *track,Double_t *p) const
+{
+  //
+  // ITS PID
+  //  1) Truncated mean method
+  //
+
+
+  Double_t mom=track->P();  
+  AliAODPid *pidObj = track->GetDetPid();
+
+  Double_t dedx=pidObj->GetITSsignal();
+  Bool_t mismatch=kTRUE;
+  for (Int_t j=0; j<AliPID::kSPECIES; j++) {
+    Double_t mass=AliPID::ParticleMass(j);//GeV/c^2
+    Double_t bethe=fITSResponse.Bethe(mom,mass);
+    Double_t sigma=fITSResponse.GetResolution(bethe);
+    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 (mismatch)
+    for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
+
+  return;
+
+}
+//_________________________________________________________________________
+void AliAODpidUtil::MakeTOFPID(AliAODTrack *track, Float_t TimeZeroTOF,Double_t *p) const
+{
+  //
+  //   TOF PID using gaussian response
+  //
+
+  Double_t time[AliPID::kSPECIESN];
+  Double_t sigma[AliPID::kSPECIESN];
+  AliAODPid *pidObj = track->GetDetPid();
+  pidObj->GetIntegratedTimes(time);
+
+  for (Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++) {
+    sigma[iPart] = fTOFResponse.GetExpectedSigma(track->P(),time[iPart],AliPID::ParticleMass(iPart));
+  }
+
+  AliDebugGeneral("AliESDpid::MakeTOFPID",2,
+                 Form("Expected TOF signals [ps]: %f %f %f %f %f",
+                      time[AliPID::kElectron],
+                      time[AliPID::kMuon],
+                      time[AliPID::kPion],
+                      time[AliPID::kKaon],
+                      time[AliPID::kProton]));
+
+  AliDebugGeneral("AliESDpid::MakeTOFPID",2,
+                 Form("Expected TOF std deviations [ps]: %f %f %f %f %f",
+                      sigma[AliPID::kElectron],
+                      sigma[AliPID::kMuon],
+                      sigma[AliPID::kPion],
+                      sigma[AliPID::kKaon],
+                      sigma[AliPID::kProton]
+                      ));
+
+  Double_t tof = pidObj->GetTOFsignal() - TimeZeroTOF;
+
+  Bool_t mismatch = kTRUE;
+  for (Int_t j=0; j<AliPID::kSPECIES; j++) {
+    Double_t sig = sigma[j];
+    if (TMath::Abs(tof-time[j]) > fRange*sig) {
+      p[j] = TMath::Exp(-0.5*fRange*fRange)/sig;
+    } else
+      p[j] = TMath::Exp(-0.5*(tof-time[j])*(tof-time[j])/(sig*sig))/sig;
+
+    // Check the mismatching
+    Double_t mass = AliPID::ParticleMass(j);
+    Double_t pm = fTOFResponse.GetMismatchProbability(track->P(),mass);
+    if (p[j]>pm) mismatch = kFALSE;
+
+    // Check for particles heavier than (AliPID::kSPECIES - 1)
+
+  }
+
+  if (mismatch)
+    for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1/AliPID::kSPECIES;
+
+  return;
+}
+//_________________________________________________________________________
+void AliAODpidUtil::MakeTRDPID(AliAODTrack *track,Double_t *p) const
+{
+  
+  // Method to recalculate the TRD PID probabilities
+  AliAODPid *pidObj = track->GetDetPid();
+  Float_t *mom=pidObj->GetTRDmomentum();
+  Double_t *dedx=pidObj->GetTRDsignal();
+  Bool_t norm=kTRUE;
+  AliTRDPIDResponse *trdResponse;
+  trdResponse->GetResponse(pidObj->GetTRDnSlices(),dedx,mom,p,norm);
+  return;
+}
diff --git a/STEER/AliAODpidUtil.h b/STEER/AliAODpidUtil.h
new file mode 100644 (file)
index 0000000..f8c4e97
--- /dev/null
@@ -0,0 +1,79 @@
+#ifndef ALIAODPIDUTIL_H
+#define ALIAODPIDUTIL_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id: AliAODpidUtil.h 38493 2010-01-26 16:33:03Z hristov $ */
+
+//-------------------------------------------------------
+//                    Combined PID class
+//                    for the AOD class
+//   Origin: Rosa Romita, GSI, r.romita@gsi.de 
+//-------------------------------------------------------
+#include <Rtypes.h>
+#include <TMatrixD.h>
+#include "AliAODTrack.h" // Needed for inline functions
+#include "AliAODPid.h" // Needed for inline functions
+#include "AliTPCPIDResponse.h"
+#include "AliITSPIDResponse.h"
+#include "AliTOFPIDResponse.h"
+#include "AliTRDPIDResponse.h"
+//#include "HMPID/AliHMPID.h"
+
+class AliAODEvent;
+
+class AliAODpidUtil {
+public:
+  AliAODpidUtil(): fRange(5.), fTPCResponse(), fITSResponse(), fTOFResponse(){;}
+  virtual ~AliAODpidUtil() {;}
+
+
+  Int_t MakePID(AliAODTrack *track,Float_t TimeZeroTOF,Double_t *p) const;
+  void MakeTPCPID(AliAODTrack *track,Double_t *p) const;
+  void MakeITSPID(AliAODTrack *track,Double_t *p) const;
+  void MakeTOFPID(AliAODTrack *track, Float_t TimeZeroTOF,Double_t *p) const;
+  //  void MakeHMPIDPID(AliESDtrack *track);
+  void MakeTRDPID(AliAODTrack *track,Double_t *p) const;
+
+  Float_t NumberOfSigmasTPC(const AliAODTrack *track, AliPID::EParticleType type) const;
+  Float_t NumberOfSigmasTOF(const AliAODTrack *track, AliPID::EParticleType type, const Float_t TimeZeroTOF) const;
+  Float_t NumberOfSigmasITS(const AliAODTrack *track, AliPID::EParticleType type) const;
+
+  AliITSPIDResponse &GetITSResponse() {return fITSResponse;}
+  AliTPCPIDResponse &GetTPCResponse() {return fTPCResponse;}
+  AliTOFPIDResponse &GetTOFResponse() {return fTOFResponse;}
+
+private:
+  Float_t           fRange;          // nSigma max in likelihood
+  AliTPCPIDResponse fTPCResponse;
+  AliITSPIDResponse fITSResponse;
+  AliTOFPIDResponse fTOFResponse;
+  // AliHMPIDPIDResponse fHMPIDResponse;
+ // AliTRDPIDResponse fTRDResponse;
+
+  ClassDef(AliAODpidUtil,1)  // PID calculation class
+};
+
+inline Float_t AliAODpidUtil::NumberOfSigmasTPC(const AliAODTrack *track, AliPID::EParticleType type) const {
+  
+  Double_t mom = track->P();
+  AliAODPid *pidObj = track->GetDetPid();
+  if (pidObj)
+    mom = pidObj->GetTPCmomentum();
+  return fTPCResponse.GetNumberOfSigmas(mom,pidObj->GetTPCsignal(),0,type); 
+}
+
+inline Float_t AliAODpidUtil::NumberOfSigmasTOF(const AliAODTrack *track, AliPID::EParticleType type, const Float_t TimeZeroTOF) const {
+  Double_t times[AliPID::kSPECIES];
+  AliAODPid *pidObj = track->GetDetPid();
+  pidObj->GetIntegratedTimes(times);
+  return (pidObj->GetTOFsignal() - TimeZeroTOF - times[type])/fTOFResponse.GetExpectedSigma(track->P(),times[type],AliPID::ParticleMass(type));
+}
+
+inline Float_t AliAODpidUtil::NumberOfSigmasITS(const AliAODTrack *track, AliPID::EParticleType type) const {
+  AliAODPid *pidObj = track->GetDetPid();
+  return fITSResponse.GetNumberOfSigmas(track->P(),pidObj->GetITSsignal(),type); 
+}
+#endif
+
+
index fb0f0dc..9ae69d0 100644 (file)
@@ -8,7 +8,7 @@ SRCS = AliAODEvent.cxx AliAODHeader.cxx \
        AliAODv0.cxx AliAODcascade.cxx AliAODCaloCells.cxx AliAODInputHandler.cxx \
        AliAODDiJet.cxx AliAODMCParticle.cxx AliAODMCHeader.cxx \
         AliAODPWG4Particle.cxx AliAODPWG4ParticleCorrelation.cxx \
-       AliAODDimuon.cxx
+       AliAODDimuon.cxx AliAODpidUtil.cxx
 
 HDRS:= $(SRCS:.cxx=.h)