- Code cleanup
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 18 Mar 2013 17:59:20 +0000 (17:59 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 18 Mar 2013 17:59:20 +0000 (17:59 +0000)
- Fix problem for light nuclei in 'TunedOnData' (Fancesco N.)
- Add HMPID PID response and NSigma calculation (Giacomo)

23 files changed:
STEER/AOD/AliAODEvent.cxx
STEER/AOD/AliAODEvent.h
STEER/AOD/AliAODHMPIDrings.cxx
STEER/AOD/AliAODHMPIDrings.h
STEER/AOD/AliAODTrack.cxx
STEER/AOD/AliAODTrack.h
STEER/AOD/AliAODpidUtil.cxx
STEER/AOD/AliAODpidUtil.h
STEER/CMakelibSTEERBase.pkg
STEER/ESD/AliESDpid.cxx
STEER/ESD/AliESDpid.h
STEER/STEERBase/AliHMPIDPIDParams.cxx [new file with mode: 0644]
STEER/STEERBase/AliHMPIDPIDParams.h [new file with mode: 0644]
STEER/STEERBase/AliHMPIDPIDResponse.cxx [new file with mode: 0644]
STEER/STEERBase/AliHMPIDPIDResponse.h [new file with mode: 0644]
STEER/STEERBase/AliITSPIDResponse.cxx
STEER/STEERBase/AliITSPIDResponse.h
STEER/STEERBase/AliPIDResponse.cxx
STEER/STEERBase/AliPIDResponse.h
STEER/STEERBase/AliTPCPIDResponse.cxx
STEER/STEERBase/AliTPCPIDResponse.h
STEER/STEERBase/AliVTrack.h
STEER/STEERBaseLinkDef.h

index df9ca0e..b1de8eb 100644 (file)
@@ -955,7 +955,7 @@ void  AliAODEvent::SetTOFHeader(const AliTOFHeader *header)
 
 }
 //------------------------------------------------------------
-AliAODHMPIDrings *AliAODEvent::GetHMPIDringForTrackID(Int_t trackID)
+AliAODHMPIDrings *AliAODEvent::GetHMPIDringForTrackID(Int_t trackID) const
 {
   //
   // Returns the HMPID object if any for a given track ID
@@ -970,7 +970,7 @@ AliAODHMPIDrings *AliAODEvent::GetHMPIDringForTrackID(Int_t trackID)
   return 0;
 }
 //------------------------------------------------------------
-Int_t AliAODEvent::GetNHMPIDrings()   
+Int_t AliAODEvent::GetNHMPIDrings() const   
 { 
   //
   // If there is a list of HMPID rings in the given AOD event, return their number
@@ -979,7 +979,7 @@ Int_t AliAODEvent::GetNHMPIDrings()
   else return -1;
 } 
 //------------------------------------------------------------
-AliAODHMPIDrings *AliAODEvent::GetHMPIDring(Int_t nRings)   
+AliAODHMPIDrings *AliAODEvent::GetHMPIDring(Int_t nRings) const
 { 
   //
   // If there is a list of HMPID rings in the given AOD event, return corresponding ring
index 0d5f2e9..3f5230c 100644 (file)
@@ -215,12 +215,12 @@ class AliAODEvent : public AliVEvent {
 
   // -- HMPID objects 
   TClonesArray *GetHMPIDrings()       const {return fHMPIDrings; } 
-  Int_t         GetNHMPIDrings();
-  AliAODHMPIDrings *GetHMPIDring(Int_t nRings);
+  Int_t         GetNHMPIDrings()      const;
+  AliAODHMPIDrings *GetHMPIDring(Int_t nRings) const;
   Int_t         AddHMPIDrings(const  AliAODHMPIDrings* ring) 
   {new((*fHMPIDrings)[fHMPIDrings->GetEntriesFast()]) AliAODHMPIDrings(*ring); return fHMPIDrings->GetEntriesFast()-1;}
   
-  AliAODHMPIDrings *GetHMPIDringForTrackID(Int_t trackID);
+  AliAODHMPIDrings *GetHMPIDringForTrackID(Int_t trackID) const;
   
   
   // -- Jet
index f80af83..78605e4 100644 (file)
@@ -138,13 +138,13 @@ AliAODHMPIDrings& AliAODHMPIDrings::operator=(const AliAODHMPIDrings& hmpidAOD)
                                   
 }
 //________________________________________________________________________________________________________________________________________________________
-void AliAODHMPIDrings::GetHmpPidProbs(Double32_t *pid)
+void AliAODHMPIDrings::GetHmpPidProbs(Double32_t *pid) const
 {
   // Gets probabilities of each particle type (in HMPID)
   for (Int_t i=0; i<AliPID::kSPECIES; i++) pid[i]=fHmpidAODpid[i];
 }
 //________________________________________________________________________________________________________________________________________________________
-void  AliAODHMPIDrings::GetHmpMom(Double32_t *mom)
+void  AliAODHMPIDrings::GetHmpMom(Double32_t *mom) const
 {
   for( Int_t ico = 0 ; ico < 3; ico++) mom[ico] = fHMPIDmom[ico];
 }
index 939295c..2b81253 100644 (file)
@@ -43,34 +43,36 @@ class AliAODHMPIDrings : public TObject {
   virtual ~AliAODHMPIDrings() {};
     
   //___ Getters
-  Double32_t GetHmpTrkID()                { return fHmpidAODtrkId; }
+  Double32_t GetHmpTrkID()                const { return fHmpidAODtrkId; }
   
-  Double32_t GetHmpMipCharge()            { return fHmpidAODqn%1000000; }
-  Double32_t GetHmpNumOfPhotonClusters()  { return fHmpidAODqn/1000000;}
+  Double32_t GetHmpMipCharge()            const { return fHmpidAODqn%1000000; }
+  Double32_t GetHmpNumOfPhotonClusters()  const { return fHmpidAODqn/1000000;}
  
-  Int_t      GetHmpChamber()              { return fHmpidAODcluIdx/1000000; }
+  Int_t      GetHmpChamber()              const { return fHmpidAODcluIdx/1000000; }
   
-  Double32_t GetHmpTrackTheta()           { return fHmpidAODtrkTheta;}
-  Double32_t GetHmpTrackPhi()             { return fHmpidAODtrkPhi;}
+  Int_t      GetHmpCluIdx()               const { return fHmpidAODcluIdx; }
   
-  Double32_t GetHmpSignal()               { return fHmpidAODsignal;}
-  Double32_t GetHmpOccupancy()            { return fHmpidAODocc;}
+  Double32_t GetHmpTrackTheta()           const { return fHmpidAODtrkTheta;}
+  Double32_t GetHmpTrackPhi()             const { return fHmpidAODtrkPhi;}
   
-  Double32_t GetHmpChi2()                 { return fHmpidAODchi2;}
+  Double32_t GetHmpSignal()               const { return fHmpidAODsignal;}
+  Double32_t GetHmpOccupancy()            const { return fHmpidAODocc;}
+  
+  Double32_t GetHmpChi2()                 const { return fHmpidAODchi2;}
 
-  Double32_t GetHmpTrackX()               { return fHmpidAODtrkX;}
-  Double32_t GetHmpTrackY()               { return fHmpidAODtrkY;}
+  Double32_t GetHmpTrackX()               const { return fHmpidAODtrkX;}
+  Double32_t GetHmpTrackY()               const { return fHmpidAODtrkY;}
 
-  Double32_t GetHmpMipX()                 { return fHmpidAODmipX;}
-  Double32_t GetHmpMipY()                 { return fHmpidAODmipY;}
+  Double32_t GetHmpMipX()                 const { return fHmpidAODmipX;}
+  Double32_t GetHmpMipY()                 const { return fHmpidAODmipY;}
 
-  Double32_t GetHmpDX()                   { return fHmpidAODmipX - fHmpidAODtrkX;}
-  Double32_t GetHmpDY()                   { return fHmpidAODmipY - fHmpidAODtrkY;}
-  Double32_t GetHmpDist()                 { return TMath::Sqrt((fHmpidAODmipX - fHmpidAODtrkX)*(fHmpidAODmipX - fHmpidAODtrkX) + (fHmpidAODmipY - fHmpidAODtrkY)*(fHmpidAODmipY - fHmpidAODtrkY));}
+  Double32_t GetHmpDX()                   const { return fHmpidAODmipX - fHmpidAODtrkX;}
+  Double32_t GetHmpDY()                   const { return fHmpidAODmipY - fHmpidAODtrkY;}
+  Double32_t GetHmpDist()                 const { return TMath::Sqrt((fHmpidAODmipX - fHmpidAODtrkX)*(fHmpidAODmipX - fHmpidAODtrkX) + (fHmpidAODmipY - fHmpidAODtrkY)*(fHmpidAODmipY - fHmpidAODtrkY));}
   
   
-  void GetHmpPidProbs(Double32_t *pid);   //defined in cxx
-  void GetHmpMom(Double32_t *mom);        //defined in cxx
+  void GetHmpPidProbs(Double32_t *pid) const;   //defined in cxx
+  void GetHmpMom(Double32_t *mom)      const;   //defined in cxx
   
   //___ Setters
   
index c044fc5..cddaff2 100644 (file)
 #include "AliLog.h"
 #include "AliExternalTrackParam.h"
 #include "AliVVertex.h"
-#include "AliAODTrack.h"
 #include "AliDetectorPID.h"
 #include "AliAODEvent.h"
+#include "AliAODHMPIDrings.h"
+
+#include "AliAODTrack.h"
 
 ClassImp(AliAODTrack)
 
@@ -795,6 +797,62 @@ void AliAODTrack::SetDetectorPID(const AliDetectorPID *pid)
 }
 
 //_____________________________________________________________________________
+Double_t AliAODTrack::GetHMPIDsignal() const
+{
+  if(fAODEvent->GetHMPIDringForTrackID(fID)) return fAODEvent->GetHMPIDringForTrackID(fID)->GetHmpSignal();
+  else return -999.;
+}
+
+//_____________________________________________________________________________
+Double_t AliAODTrack::GetHMPIDoccupancy() const
+{
+  if(fAODEvent->GetHMPIDringForTrackID(fID)) return fAODEvent->GetHMPIDringForTrackID(fID)->GetHmpOccupancy();
+  else return -999.;
+}
+
+//_____________________________________________________________________________
+Int_t AliAODTrack::GetHMPIDcluIdx() const
+{
+  if(fAODEvent->GetHMPIDringForTrackID(fID)) return fAODEvent->GetHMPIDringForTrackID(fID)->GetHmpCluIdx();
+  else return -999;
+}
+
+//_____________________________________________________________________________
+void AliAODTrack::GetHMPIDtrk(Float_t &x, Float_t &y, Float_t &th, Float_t &ph) const
+{
+  x = -999; y = -999.; th = -999.; ph = -999.;
+
+  const AliAODHMPIDrings *ring=fAODEvent->GetHMPIDringForTrackID(fID);
+  if(ring){
+    x  = ring->GetHmpTrackX();
+    y  = ring->GetHmpTrackY();
+    th = ring->GetHmpTrackTheta();
+    ph = ring->GetHmpTrackPhi();
+  }
+}
+
+//_____________________________________________________________________________
+void AliAODTrack::GetHMPIDmip(Float_t &x,Float_t &y,Int_t &q, Int_t &nph) const
+{
+  x = -999; y = -999.; q = -999; nph = -999;
+  
+  const AliAODHMPIDrings *ring=fAODEvent->GetHMPIDringForTrackID(fID);
+  if(ring){
+    x   = ring->GetHmpMipX();
+    y   = ring->GetHmpMipY();
+    q   = (Int_t)ring->GetHmpMipCharge();
+    nph = (Int_t)ring->GetHmpNumOfPhotonClusters();
+  }
+}
+
+//_____________________________________________________________________________
+Bool_t AliAODTrack::GetOuterHmpPxPyPz(Double_t *p) const 
+{ 
+ if(fAODEvent->GetHMPIDringForTrackID(fID)) {fAODEvent->GetHMPIDringForTrackID(fID)->GetHmpMom(p); return kTRUE;}
+ else return kFALSE;      
+}
+//_____________________________________________________________________________
 Bool_t AliAODTrack::GetXYZAt(Double_t x, Double_t b, Double_t *r) const
 {
   //---------------------------------------------------------------------
index c198cb3..d429f41 100644 (file)
@@ -284,7 +284,15 @@ class AliAODTrack : public AliVTrack {
   virtual AliTPCdEdxInfo* GetTPCdEdxInfo() const {return fDetPid?fDetPid->GetTPCdEdxInfo():0;}
   Double_t  GetTPCmomentum()     const { return fDetPid?fDetPid->GetTPCmomentum():0.;  }
   Double_t  GetTOFsignal()       const { return fDetPid?fDetPid->GetTOFsignal():0.;    }
-  Double_t  GetHMPIDsignal()     const { return 0.;  } // TODO: To be implemented properly with the new HMPID object
+  Double_t  GetHMPIDsignal()     const; 
+  Double_t  GetHMPIDoccupancy()  const;
+
+  Int_t     GetHMPIDcluIdx()     const;
+    
+  void GetHMPIDtrk(Float_t &x, Float_t &y, Float_t &th, Float_t &ph) const;  
+  void GetHMPIDmip(Float_t &x,Float_t &y,Int_t &q,Int_t &nph) const;
+  
+  Bool_t GetOuterHmpPxPyPz(Double_t *p) const;
   
   void      GetIntegratedTimes(Double_t *times) const {if (fDetPid) fDetPid->GetIntegratedTimes(times); }
   Double_t  GetTRDslice(Int_t plane, Int_t slice) const;
index 757973a..4cfc10d 100644 (file)
@@ -95,7 +95,7 @@ Float_t AliAODpidUtil::GetTPCsignalTunedOnData(const AliVTrack *t) const {
         Double_t sigma = fTPCResponse.GetExpectedSigma(track, type, AliTPCPIDResponse::kdEdxDefault, this->UseTPCEtaCorrection());
         dedx = gRandom->Gaus(bethe,sigma);
         
-           if(iS == AliPID::ParticleCode(AliPID::kHe3) || iS == AliPID::ParticleCode(AliPID::kAlpha)) dedx *= 5;
+//         if(iS == AliPID::ParticleCode(AliPID::kHe3) || iS == AliPID::ParticleCode(AliPID::kAlpha)) dedx *= 5;
        }
 
     }
@@ -103,6 +103,36 @@ Float_t AliAODpidUtil::GetTPCsignalTunedOnData(const AliVTrack *t) const {
     track->SetTPCsignalTunedOnData(dedx);
     return dedx;
 }
+
+//_________________________________________________________________________
+Float_t AliAODpidUtil::GetSignalDeltaTOFold(const AliVParticle *vtrack, AliPID::EParticleType type) const
+{
+  //
+  // Number of sigma implementation for the TOF
+  //
+  
+  AliAODTrack *track=(AliAODTrack*)vtrack;
+  
+  AliAODPid *pidObj = track->GetDetPid();
+  if (!pidObj) return -9999.;
+  Double_t tofTime=pidObj->GetTOFsignal();
+  const Double_t expTime=fTOFResponse.GetExpectedSignal((AliVTrack*)vtrack,type);
+  Double_t sigmaTOFPid[AliPID::kSPECIES];
+  pidObj->GetTOFpidResolution(sigmaTOFPid);
+  AliAODEvent *event=(AliAODEvent*)track->GetAODEvent();
+  if (event) {  // protection if the user didn't call GetTrack, which sets the internal pointer
+    AliTOFHeader* tofH=(AliTOFHeader*)event->GetTOFHeader();
+    if (tofH && (TMath::Abs(sigmaTOFPid[0]) <= 1.E-16) ) { // new AOD
+        tofTime -= fTOFResponse.GetStartTime(vtrack->P());
+    }
+  } else {
+    AliError("pointer to AliAODEvent not found, please call GetTrack to set it");
+    return -9999.;
+  }
+  
+  return tofTime - expTime;
+}
+
 //_________________________________________________________________________
 Float_t AliAODpidUtil::GetNumberOfSigmasTOFold(const AliVParticle *vtrack, AliPID::EParticleType type) const
 {
@@ -113,7 +143,7 @@ Float_t AliAODpidUtil::GetNumberOfSigmasTOFold(const AliVParticle *vtrack, AliPI
   AliAODTrack *track=(AliAODTrack*)vtrack;
 
   Bool_t oldAod=kTRUE;
-  Double_t sigTOF;
+  Double_t sigTOF=0.;
   AliAODPid *pidObj = track->GetDetPid();
   if (!pidObj) return -999.;
   Double_t tofTime=pidObj->GetTOFsignal();
index aadad04..bd13dd9 100644 (file)
@@ -36,6 +36,7 @@ public:
   Float_t GetTPCsignalTunedOnData(const AliVTrack *t) const;
 
 protected:
+  virtual Float_t GetSignalDeltaTOFold(const AliVParticle *track, AliPID::EParticleType type) const;
   virtual Float_t GetNumberOfSigmasTOFold(const AliVParticle *vtrack, AliPID::EParticleType type) const;
   
 private:
index e1d56ed..2366ac2 100644 (file)
@@ -71,6 +71,8 @@ set ( SRCS
     STEERBase/AliTOFPIDParams.cxx 
     STEERBase/AliTRDPIDResponse.cxx 
     STEERBase/AliEMCALPIDResponse.cxx 
+    STEERBase/AliHMPIDPIDResponse.cxx
+    STEERBase/AliHMPIDPIDParams.cxx
     STEERBase/AliPIDCombined.cxx
     STEERBase/AliPIDValues.cxx
     STEERBase/AliDetectorPID.cxx
index 1c3473f..8e3c077 100644 (file)
@@ -126,7 +126,7 @@ Float_t AliESDpid::GetTPCsignalTunedOnData(const AliVTrack *t) const {
         Double_t bethe = fTPCResponse.GetExpectedSignal(track, type, AliTPCPIDResponse::kdEdxDefault, this->UseTPCEtaCorrection());
         Double_t sigma = fTPCResponse.GetExpectedSigma(track, type, AliTPCPIDResponse::kdEdxDefault, this->UseTPCEtaCorrection());
                dedx = gRandom->Gaus(bethe,sigma);
-               if(iS == AliPID::ParticleCode(AliPID::kHe3) || iS == AliPID::ParticleCode(AliPID::kAlpha)) dedx *= 5;
+//             if(iS == AliPID::ParticleCode(AliPID::kHe3) || iS == AliPID::ParticleCode(AliPID::kAlpha)) dedx *= 5;
            }
        }
     }
@@ -410,6 +410,19 @@ Bool_t AliESDpid::CheckTOFMatching(AliESDtrack *track) const{
     return status;
 }
 
+//_________________________________________________________________________
+Float_t AliESDpid::GetSignalDeltaTOFold(const AliVParticle *track, AliPID::EParticleType type) const
+{
+  //
+  // TOF signal - expected
+  //
+  AliVTrack *vtrack=(AliVTrack*)track;
+  
+  Double_t expTime = fTOFResponse.GetExpectedSignal(vtrack,type);
+  return (vtrack->GetTOFsignal() - fTOFResponse.GetStartTime(vtrack->P()) - expTime);
+}
+
+//_________________________________________________________________________
 Float_t AliESDpid::GetNumberOfSigmasTOFold(const AliVParticle *track, AliPID::EParticleType type) const
 {
   //
index 3a50ebc..cd51b36 100644 (file)
@@ -50,6 +50,7 @@ AliESDpid(const AliESDpid&a): AliPIDResponse(a), fRangeTOFMismatch(a.fRangeTOFMi
 
   void SetEventHandler(AliVEventHandler *event){fEventHandler=event;};
 protected:
+  virtual Float_t GetSignalDeltaTOFold(const AliVParticle *track, AliPID::EParticleType type) const;
   virtual Float_t GetNumberOfSigmasTOFold(const AliVParticle *track, AliPID::EParticleType type) const;
 
 private:
diff --git a/STEER/STEERBase/AliHMPIDPIDParams.cxx b/STEER/STEERBase/AliHMPIDPIDParams.cxx
new file mode 100644 (file)
index 0000000..4ab5250
--- /dev/null
@@ -0,0 +1,70 @@
+/**************************************************************************
+ * Copyright(c) 1998-2010, 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.                  *
+ **************************************************************************/
+
+//***********************************************************
+// Class AliHMPIDPIDParams
+//
+// class to store PID parameters for HMPID in OADB
+//
+// Author: G. Volpe, giacomo.volpe@cern.ch
+//***********************************************************
+
+#include <TNamed.h>
+#include "AliHMPIDPIDParams.h"
+
+ClassImp(AliHMPIDPIDParams)
+
+//_____________________________________________________________________________
+AliHMPIDPIDParams::AliHMPIDPIDParams():
+ TNamed("default",""),
+ fHMPIDRefIndexArray(0x0)
+{
+  
+}
+//_____________________________________________________________________________
+AliHMPIDPIDParams::AliHMPIDPIDParams(Char_t *name):
+  TNamed(name,""),
+  fHMPIDRefIndexArray(0x0)
+{
+  
+}
+//___________________________________________________________________________
+AliHMPIDPIDParams& AliHMPIDPIDParams::operator=(const AliHMPIDPIDParams& c)
+{
+  //
+  // Assignment operator
+  //
+  if (this!=&c) {
+    fHMPIDRefIndexArray=c.fHMPIDRefIndexArray;
+  }
+  return *this;
+}
+
+//___________________________________________________________________________
+AliHMPIDPIDParams::AliHMPIDPIDParams(const AliHMPIDPIDParams& c) :
+ TNamed(c),
+ fHMPIDRefIndexArray(c.fHMPIDRefIndexArray)   
+{
+  //
+  // Copy Constructor
+  //
+  
+}
+//_____________________________________________________________________________
+AliHMPIDPIDParams::~AliHMPIDPIDParams(){
+}
+
+
diff --git a/STEER/STEERBase/AliHMPIDPIDParams.h b/STEER/STEERBase/AliHMPIDPIDParams.h
new file mode 100644 (file)
index 0000000..d385041
--- /dev/null
@@ -0,0 +1,36 @@
+#ifndef ALIHMPIDPIDPARAMS_H
+#define ALIHMPIDPIDPARAMS_H
+/* Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+//***********************************************************
+// Class AliHMPIDPIDparams
+// class to store PID parameters for HMPID in OADB
+// Author: G. Volpe, giacomo.volpe@cern.ch
+//***********************************************************
+
+#include <TNamed.h>
+
+class TObjArray;
+
+class AliHMPIDPIDParams : public TNamed {
+
+ public:
+  AliHMPIDPIDParams();
+  AliHMPIDPIDParams(Char_t * name);
+  AliHMPIDPIDParams& operator= (const AliHMPIDPIDParams& c);
+  AliHMPIDPIDParams(const AliHMPIDPIDParams& c);
+  virtual ~AliHMPIDPIDParams();
+
+  TObjArray*  GetHMPIDrefIndex() const     {return fHMPIDRefIndexArray;}  
+  void SetHMPIDrefIndex(TObjArray *array)  {fHMPIDRefIndexArray = array;}
+
+ private:
+  TObjArray   *fHMPIDRefIndexArray;                           // C6F14 refractive index
+
+  ClassDef(AliHMPIDPIDParams,1);
+
+};
+
+#endif
+
diff --git a/STEER/STEERBase/AliHMPIDPIDResponse.cxx b/STEER/STEERBase/AliHMPIDPIDResponse.cxx
new file mode 100644 (file)
index 0000000..45295f8
--- /dev/null
@@ -0,0 +1,586 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+//***********************************************************
+// Class AliHMPIDPIDResponse
+//
+// HMPID class to perfom particle identification
+// 
+// Author: G. Volpe, giacomo.volpe@cern.ch
+//***********************************************************
+
+
+#include "AliHMPIDPIDResponse.h"  //class header
+#include "AliPID.h"               //FindPid()
+#include "AliVTrack.h"            //FindPid()
+#include "AliLog.h"               //general
+#include <TRandom.h>              //Resolution()
+#include <TVector2.h>             //Resolution()
+#include <TRotation.h>
+#include <TF1.h>
+#include <TGeoManager.h>          //Instance()
+#include <TGeoMatrix.h>           //Instance()
+#include <TGeoPhysicalNode.h>     //ctor
+#include <TGeoBBox.h>
+#include <TObjArray.h>
+
+Float_t AliHMPIDPIDResponse::fgkMinPcX[]={0.,0.,0.,0.,0.,0.};
+Float_t AliHMPIDPIDResponse::fgkMaxPcX[]={0.,0.,0.,0.,0.,0.};
+Float_t AliHMPIDPIDResponse::fgkMinPcY[]={0.,0.,0.,0.,0.,0.};
+Float_t AliHMPIDPIDResponse::fgkMaxPcY[]={0.,0.,0.,0.,0.,0.};
+
+Float_t AliHMPIDPIDResponse::fgCellX=0.;
+Float_t AliHMPIDPIDResponse::fgCellY=0.;
+
+Float_t AliHMPIDPIDResponse::fgPcX=0;
+Float_t AliHMPIDPIDResponse::fgPcY=0;
+
+Float_t AliHMPIDPIDResponse::fgAllX=0;
+Float_t AliHMPIDPIDResponse::fgAllY=0;
+
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+AliHMPIDPIDResponse::AliHMPIDPIDResponse():
+  TNamed("HMPIDPIDResponseRec","HMPIDPIDResponsePid"),    
+  fRefIdx(1.28947),
+  fTrkDir(0,0,1),  
+  fTrkPos(30,40),  
+  fRefIndexArray(0x0)
+{
+  //
+  // ctor
+  //
+  
+  Float_t dead=2.6;// cm of the dead zones between PCs-> See 2CRC2099P1
+
+  fgCellX=0.8; fgCellY=0.84;
+  
+  fgPcX  = 80.*fgCellX;      fgPcY = 48.*fgCellY;
+  fgAllX = 2.*fgPcX+dead;
+  fgAllY = 3.*fgPcY+2.*dead;
+
+  fgkMinPcX[1]=fgPcX+dead; fgkMinPcX[3]=fgkMinPcX[1];  fgkMinPcX[5]=fgkMinPcX[3];
+  fgkMaxPcX[0]=fgPcX;      fgkMaxPcX[2]=fgkMaxPcX[0];  fgkMaxPcX[4]=fgkMaxPcX[2];
+  fgkMaxPcX[1]=fgAllX;     fgkMaxPcX[3]=fgkMaxPcX[1];  fgkMaxPcX[5]=fgkMaxPcX[3];
+
+  fgkMinPcY[2]=fgPcY+dead;       fgkMinPcY[3]=fgkMinPcY[2];  
+  fgkMinPcY[4]=2.*fgPcY+2.*dead; fgkMinPcY[5]=fgkMinPcY[4];
+  fgkMaxPcY[0]=fgPcY;            fgkMaxPcY[1]=fgkMaxPcY[0];  
+  fgkMaxPcY[2]=2.*fgPcY+dead;    fgkMaxPcY[3]=fgkMaxPcY[2]; 
+  fgkMaxPcY[4]=fgAllY;           fgkMaxPcY[5]=fgkMaxPcY[4];   
+    
+  for(Int_t i=kMinCh;i<=kMaxCh;i++)
+    if(gGeoManager && gGeoManager->IsClosed()) {
+      TGeoPNEntry* pne = gGeoManager->GetAlignableEntry(Form("/HMPID/Chamber%i",i));
+      if (!pne) {
+        AliErrorClass(Form("The symbolic volume %s does not correspond to any physical entry!",Form("HMPID_%i",i)));
+        fM[i]=new TGeoHMatrix;
+        IdealPosition(i,fM[i]);
+      } else {
+        TGeoPhysicalNode *pnode = pne->GetPhysicalNode();
+        if(pnode) fM[i]=new TGeoHMatrix(*(pnode->GetMatrix()));
+        else {
+          fM[i]=new TGeoHMatrix;
+          IdealPosition(i,fM[i]);
+        }
+      }
+    } else{
+      fM[i]=new TGeoHMatrix;
+      IdealPosition(i,fM[i]);
+    } 
+    
+}//ctor
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+AliHMPIDPIDResponse::AliHMPIDPIDResponse(const AliHMPIDPIDResponse& c):
+  TNamed(c), 
+  fRefIdx(c.fRefIdx),
+  fTrkDir(c.fTrkDir),  
+  fTrkPos(c.fTrkPos),  
+  fRefIndexArray(c.fRefIndexArray)
+ {
+   //
+   // copy ctor
+   //
+   
+   for(Int_t i=0; i<6; i++) {
+      
+      fgkMinPcX[i] = c.fgkMinPcX[i];
+      fgkMinPcY[i] = c.fgkMinPcY[i];
+      fgkMaxPcX[i] = c.fgkMaxPcX[i];
+      fgkMaxPcY[i] = c.fgkMaxPcY[i];
+     }
+   
+   for(Int_t i=0; i<7; i++) fM[i] = c.fM[i];
+ }
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+AliHMPIDPIDResponse& AliHMPIDPIDResponse::operator=(const AliHMPIDPIDResponse& c) {
+
+   //
+   // assignment operator
+   //       
+  if(this!=&c){
+     TNamed::operator=(c);
+     fgCellX = c.fgCellX;
+     fgCellY = c.fgCellY;
+     fgPcX   = c.fgPcX;
+     fgPcY   = c.fgPcY;
+     fgAllX  = c.fgAllX;
+     fgAllY  = c.fgAllY;
+     fRefIdx = c.fRefIdx;
+     fTrkDir = c.fTrkDir;  
+     fTrkPos = c.fTrkPos;  
+     fRefIndexArray = c.fRefIndexArray;
+     for(Int_t i=0; i<6; i++) {    
+      fgkMinPcX[i] = c.fgkMinPcX[i];
+      fgkMinPcY[i] = c.fgkMinPcY[i];
+      fgkMaxPcX[i] = c.fgkMaxPcX[i];
+      fgkMaxPcY[i] = c.fgkMaxPcY[i];
+     }   
+     for(Int_t i=0; i<7; i++) fM[i] = c.fM[i];                 
+    } 
+    
+  return *this; 
+}    
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+void AliHMPIDPIDResponse::IdealPosition(Int_t iCh, TGeoHMatrix *pMatrix) {
+
+// Construct ideal position matrix for a given chamber
+// Arguments: iCh- chamber ID; pMatrix- pointer to precreated unity matrix where to store the results
+// Returns: none
+
+  const Double_t kAngHor=19.5;        //  horizontal angle between chambers  19.5 grad
+  const Double_t kAngVer=20;          //  vertical angle between chambers    20   grad     
+  const Double_t kAngCom=30;          //  common HMPID rotation with respect to x axis  30   grad     
+  const Double_t kTrans[3]={490,0,0}; //  center of the chamber is on window-gap surface
+  pMatrix->RotateY(90);               //  rotate around y since initial position is in XY plane -> now in YZ plane
+  pMatrix->SetTranslation(kTrans);    //  now plane in YZ is shifted along x 
+  switch(iCh){
+    case 0:                pMatrix->RotateY(kAngHor);  pMatrix->RotateZ(-kAngVer);  break; //right and down 
+    case 1:                                            pMatrix->RotateZ(-kAngVer);  break; //down              
+    case 2:                pMatrix->RotateY(kAngHor);                               break; //right 
+    case 3:                                                                         break; //no rotation
+    case 4:                pMatrix->RotateY(-kAngHor);                              break; //left   
+    case 5:                                            pMatrix->RotateZ(kAngVer);   break; //up
+    case 6:                pMatrix->RotateY(-kAngHor); pMatrix->RotateZ(kAngVer);   break; //left and up 
+  }
+  pMatrix->RotateZ(kAngCom);     //apply common rotation  in XY plane    
+   
+}
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+Double_t AliHMPIDPIDResponse::GetExpectedSignal(const AliVTrack *vTrk, AliPID::EParticleType specie) const {
+  
+  // expected Cherenkov angle calculation
+  
+  const Double_t nmean = GetNMean(vTrk);
+  return ExpectedSignal(vTrk,nmean,specie);
+}
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+Double_t AliHMPIDPIDResponse::GetExpectedSigma(const AliVTrack *vTrk, AliPID::EParticleType specie) const {
+  
+  // expected resolution calculation
+  
+  const Double_t nmean = GetNMean(vTrk);
+  return ExpectedSigma(vTrk,nmean,specie);
+}
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+Double_t AliHMPIDPIDResponse::ExpectedSignal(const AliVTrack *vTrk, Double_t nmean, AliPID::EParticleType specie) const {
+  
+  // expected Cherenkov angle calculation
+  
+  Double_t thetaTheor = -999.;
+  
+  Double_t p[3] = {0}, mom = 0;
+  if(vTrk->GetOuterHmpPxPyPz(p))  mom = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);  // Momentum of the charged particle    
+  else return thetaTheor;
+    
+  const Double_t mass = AliPID::ParticleMass(specie); 
+  const Double_t cosTheta = TMath::Sqrt(mass*mass+mom*mom)/(nmean*mom);
+   
+  if(cosTheta>1) return thetaTheor;
+                  
+  else thetaTheor = TMath::ACos(cosTheta);
+  
+  return thetaTheor;                                                                                          //  evaluate the theor. Theta Cherenkov
+}
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+Double_t AliHMPIDPIDResponse::ExpectedSigma(const AliVTrack *vTrk, Double_t nmean, AliPID::EParticleType specie) const {
+  
+  // expected resolution calculation
+  
+  Float_t x=0., y=0.;
+  Int_t q=0, nph=0;
+  Float_t xPc=0.,yPc=0.,thRa=0.,phRa=0.;
+  
+  vTrk->GetHMPIDmip(x,y,q,nph);  
+  vTrk->GetHMPIDtrk(xPc,yPc,thRa,phRa);
+      
+  const Double_t xRa = xPc - (RadThick()+WinThick()+GapThick())*TMath::Cos(phRa)*TMath::Tan(thRa);  //just linear extrapolation back to RAD
+  const Double_t yRa = yPc - (RadThick()+WinThick()+GapThick())*TMath::Sin(phRa)*TMath::Tan(thRa);  //just linear extrapolation back to RAD
+  
+  const Double_t thetaCerTh = ExpectedSignal(vTrk,nmean,specie);
+  const Double_t occupancy  = vTrk->GetHMPIDoccupancy();
+  const Double_t thetaMax   = TMath::ACos(1./nmean);
+  const Int_t    nPhotsTh   = (Int_t)(12.*TMath::Sin(thetaCerTh)*TMath::Sin(thetaCerTh)/(TMath::Sin(thetaMax)*TMath::Sin(thetaMax))+0.01);
+
+  Double_t sigmatot = 0;
+  Int_t nTrks = 20;
+  for(Int_t iTrk=0;iTrk<nTrks;iTrk++) {
+    Double_t invSigma = 0;
+    Int_t nPhotsAcc = 0;
+    
+    Int_t nPhots = 0; 
+    if(nph<nPhotsTh+TMath::Sqrt(nPhotsTh) && nph>nPhotsTh-TMath::Sqrt(nPhotsTh)) nPhots = nph;
+    else nPhots = gRandom->Poisson(nPhotsTh);
+    
+    for(Int_t j=0;j<nPhots;j++){
+      Double_t phi = gRandom->Rndm()*TMath::TwoPi();
+      TVector2 pos; pos = TracePhot(xRa,yRa,thRa,phRa,thetaCerTh,phi);
+      if(!IsInside(pos.X(),pos.Y())) continue;
+      if(IsInDead(pos.X(),pos.Y()))  continue;
+      Double_t sigma2 = Sigma2(thRa,thRa,thetaCerTh,phi); //photon candidate sigma^2
+      
+      if(sigma2!=0) {
+        invSigma += 1./sigma2;
+        nPhotsAcc++;
+      }
+    }      
+    if(invSigma!=0) sigmatot += 1./TMath::Sqrt(invSigma);  
+  }
+    
+  return (sigmatot/nTrks)*SigmaCorrFact(specie,occupancy);
+}
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+Double_t AliHMPIDPIDResponse::GetNumberOfSigmas(const AliVTrack *vTrk, AliPID::EParticleType specie) const {
+
+  // Number of sigmas calculation
+    
+  Double_t nSigmas = -999.;
+
+  if(vTrk->GetHMPIDsignal()<0.) return nSigmas;
+    
+  const Double_t nmean = GetNMean(vTrk);
+  
+  const Double_t expSigma = ExpectedSigma(vTrk, nmean, specie);
+  
+  if(expSigma > 0.) nSigmas = (vTrk->GetHMPIDsignal() - ExpectedSignal(vTrk,nmean,specie))/expSigma;
+                         
+  return nSigmas;
+    
+}
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+void AliHMPIDPIDResponse::GetProbability(const AliVTrack *vTrk,Int_t nSpecies,Double_t *prob) const {
+
+// Calculates probability to be a electron-muon-pion-kaon-proton with the "amplitude" method
+// from the given Cerenkov angle and momentum assuming no initial particle composition
+  
+  const Double_t thetaCerExp = vTrk->GetHMPIDsignal();                                                                           
+
+  const Double_t nmean = GetNMean(vTrk);
+    
+  if(thetaCerExp<=0){                                                                                     // HMPID does not find anything reasonable for this track, assign 0.2 for all species
+    for(Int_t iPart=0;iPart<nSpecies;iPart++) prob[iPart]=1.0/(Float_t)nSpecies;
+    return;
+  } 
+  
+  Double_t p[3] = {0,0,0};
+  
+  if(!(vTrk->GetOuterHmpPxPyPz(p))) for(Int_t iPart=0;iPart<nSpecies;iPart++) prob[iPart]=1.0/(Float_t)nSpecies;
+  
+  Double_t hTot=0;                                                                                        // Initialize the total height of the amplitude method
+  Double_t *h = new Double_t [nSpecies];                                                                  // number of charged particles to be considered
+
+  Bool_t desert = kTRUE;                                                                                  // Flag to evaluate if ThetaC is far ("desert") from the given Gaussians
+  
+  for(Int_t iPart=0;iPart<nSpecies;iPart++){                                                              // for each particle
+
+        
+    h[iPart] = 0;                                                                                         // reset the height
+    Double_t thetaCerTh = ExpectedSignal(vTrk,nmean,(AliPID::EParticleType)iPart);                        // theoretical Theta Cherenkov
+    if(thetaCerTh>900.) continue;                                                                         // no light emitted, zero height
+    Double_t sigmaRing = ExpectedSigma(vTrk,nmean,(AliPID::EParticleType)iPart);
+    
+    if(sigmaRing==0) continue;
+      
+    if(TMath::Abs(thetaCerExp-thetaCerTh)<4*sigmaRing) desert = kFALSE;                                                               
+    h[iPart] =TMath::Gaus(thetaCerTh,thetaCerExp,sigmaRing,kTRUE);
+    hTot    +=h[iPart];                                                                                   // total height of all theoretical heights for normalization
+    
+  }//species loop
+
+  for(Int_t iPart=0;iPart<nSpecies;iPart++) {                                                             // species loop to assign probabilities
+     
+    if(!desert) prob[iPart]=h[iPart]/hTot;
+    else        prob[iPart]=1.0/(Float_t)nSpecies;                                                        // all theoretical values are far away from experemental one
+    
+  }
+  
+  delete [] h;
+}
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+Double_t AliHMPIDPIDResponse::GetSignalDelta(const AliVTrack *vTrk, AliPID::EParticleType specie) const {
+  
+  //
+  // calculation of Experimental Cherenkov angle - Theoretical Cherenkov angle  
+  //
+  const Double_t delta = vTrk->GetHMPIDsignal() - GetExpectedSignal(vTrk,specie);
+  
+  return delta;
+  
+}
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+TVector2 AliHMPIDPIDResponse::TracePhot(Double_t xRa, Double_t yRa, Double_t thRa, Double_t phRa, Double_t ckovThe,Double_t ckovPhi) const {
+
+// Trace a single Ckov photon from emission point somewhere in radiator up to photocathode taking into account ref indexes of materials it travereses
+// Returns: distance between photon point on PC and track projection  
+  
+  Double_t theta=0.,phi=0.;
+  TVector3  dirTRS,dirLORS;
+  dirTRS.SetMagThetaPhi(1,ckovThe,ckovPhi);                     //photon in TRS
+  Trs2Lors(thRa,phRa,dirTRS,theta,phi);
+  dirLORS.SetMagThetaPhi(1,theta,phi);                          //photon in LORS
+  return TraceForward(xRa,yRa,dirLORS);                                 //now foward tracing
+}//TracePhot()
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+TVector2 AliHMPIDPIDResponse::TraceForward(Double_t xRa, Double_t yRa, TVector3 dirCkov) const {
+
+// Trace forward a photon from (x,y) up to PC
+// Returns: pos of traced photon at PC
+  
+  TVector2 pos(-999,-999);
+  Double_t thetaCer = dirCkov.Theta();
+  if(thetaCer > TMath::ASin(1./GetRefIdx())) return pos;          //total refraction on WIN-GAP boundary
+  Double_t zRad= -0.5*RadThick()-0.5*WinThick();          //z position of middle of RAD
+  TVector3  posCkov(xRa,yRa,zRad);                        //RAD: photon position is track position @ middle of RAD 
+  Propagate(dirCkov,posCkov,           -0.5*WinThick());          //go to RAD-WIN boundary  
+  Refract  (dirCkov,         GetRefIdx(),WinIdx());       //RAD-WIN refraction
+  Propagate(dirCkov,posCkov,            0.5*WinThick());          //go to WIN-GAP boundary
+  Refract  (dirCkov,         WinIdx(),GapIdx());          //WIN-GAP refraction
+  Propagate(dirCkov,posCkov,0.5*WinThick()+GapThick());   //go to PC
+  pos.Set(posCkov.X(),posCkov.Y());
+  return pos;
+}//TraceForward()
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+void AliHMPIDPIDResponse::Propagate(const TVector3 dir,TVector3 &pos,Double_t z) const {
+  
+// Finds an intersection point between a line and XY plane shifted along Z.
+// Arguments:  dir,pos   - vector along the line and any point of the line
+//             z         - z coordinate of plain 
+// Returns:  none
+// On exit:  pos is the position if this intesection if any
+
+  static TVector3 nrm(0,0,1); 
+         TVector3 pnt(0,0,z);
+  
+  TVector3 diff=pnt-pos;
+  Double_t sint=(nrm*diff)/(nrm*dir);
+  pos+=sint*dir;
+}//Propagate()
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+void AliHMPIDPIDResponse::Refract(TVector3 &dir,Double_t n1,Double_t n2) const {
+
+// Refract direction vector according to Snell law
+// Arguments: 
+//            n1 - ref idx of first substance
+//            n2 - ref idx of second substance
+//   Returns: none
+//   On exit: dir is new direction
+  
+  Double_t sinref=(n1/n2)*TMath::Sin(dir.Theta());
+  if(TMath::Abs(sinref)>1.) dir.SetXYZ(-999,-999,-999);
+  else             dir.SetTheta(TMath::ASin(sinref));
+}//Refract()
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+void AliHMPIDPIDResponse::Trs2Lors(Double_t thRa, Double_t phRa, TVector3 dirCkov,Double_t &thetaCer,Double_t &phiCer) const {
+
+  // Theta Cerenkov reconstruction 
+  // Returns: thetaCer of photon in LORS
+  //          phiCer of photon in LORS
+  
+  TRotation mtheta;   mtheta.RotateY(thRa);
+  TRotation mphi;       mphi.RotateZ(phRa);
+  TRotation mrot=mphi*mtheta;
+  TVector3 dirCkovLORS;
+  dirCkovLORS=mrot*dirCkov;
+  phiCer  = dirCkovLORS.Phi();                                          //actual value of the phi of the photon
+  thetaCer= dirCkovLORS.Theta();                                        //actual value of thetaCerenkov of the photon
+}
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+Bool_t AliHMPIDPIDResponse::IsInDead(Float_t x,Float_t y)  {
+  
+// Check is the current point is outside of sensitive area or in dead zones
+// Arguments: x,y -position
+// Returns: 1 if not in sensitive zone
+             
+  for(Int_t iPc=0;iPc<6;iPc++)
+    if(x>=fgkMinPcX[iPc] && x<=fgkMaxPcX[iPc] && y>=fgkMinPcY[iPc] && y<=fgkMaxPcY [iPc]) return kFALSE; //in current pc
+  
+  return kTRUE;
+}
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+Double_t AliHMPIDPIDResponse::Sigma2(Double_t trkTheta,Double_t trkPhi,Double_t ckovTh, Double_t ckovPh) const {
+  
+// Analithical calculation of total error (as a sum of localization, geometrical and chromatic errors) on Cerenkov angle for a given Cerenkov photon 
+// created by a given MIP. Formules according to CERN-EP-2000-058 
+// Arguments: Cerenkov and azimuthal angles for Cerenkov photon, [radians]
+//            dip and azimuthal angles for MIP taken at the entrance to radiator, [radians]        
+//            MIP beta
+// Returns: absolute error on Cerenkov angle, [radians]    
+  
+  TVector3 v(-999,-999,-999);
+  Double_t trkBeta = 1./(TMath::Cos(ckovTh)*GetRefIdx());
+  
+  if(trkBeta > 1) trkBeta = 1;                 //protection against bad measured thetaCer  
+  if(trkBeta < 0) trkBeta = 0.0001;            //
+
+  v.SetX(SigLoc (trkTheta,trkPhi,ckovTh,ckovPh,trkBeta));
+  v.SetY(SigGeom(trkTheta,trkPhi,ckovTh,ckovPh,trkBeta));
+  v.SetZ(SigCrom(trkTheta,ckovTh,ckovPh,trkBeta));
+
+  return v.Mag2();
+}
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+Double_t AliHMPIDPIDResponse::SigLoc(Double_t trkTheta,Double_t trkPhi,Double_t thetaC, Double_t phiC,Double_t betaM) const {
+  
+// Analitical calculation of localization error (due to finite segmentation of PC) on Cerenkov angle for a given Cerenkov photon 
+// created by a given MIP. Fromulae according to CERN-EP-2000-058 
+// Arguments: Cerenkov and azimuthal angles for Cerenkov photon, [radians]
+//            dip and azimuthal angles for MIP taken at the entrance to radiator, [radians]        
+//            MIP beta
+// Returns: absolute error on Cerenkov angle, [radians]    
+  
+  Double_t phiDelta = phiC;
+
+  Double_t sint     = TMath::Sin(trkTheta);
+  Double_t cost     = TMath::Cos(trkTheta);
+  Double_t sinf     = TMath::Sin(trkPhi);
+  Double_t cosf     = TMath::Cos(trkPhi);
+  Double_t sinfd    = TMath::Sin(phiDelta);
+  Double_t cosfd    = TMath::Cos(phiDelta);
+  Double_t tantheta = TMath::Tan(thetaC);
+  
+  Double_t alpha =cost-tantheta*cosfd*sint;                                                    // formula (11)
+  Double_t k = 1.-GetRefIdx()*GetRefIdx()+alpha*alpha/(betaM*betaM);                           // formula (after 8 in the text)
+  if (k<0) return 1e10;
+  Double_t mu =sint*sinf+tantheta*(cost*cosfd*sinf+sinfd*cosf);                                // formula (10)
+  Double_t e  =sint*cosf+tantheta*(cost*cosfd*cosf-sinfd*sinf);                                // formula (9)
+
+  Double_t kk = betaM*TMath::Sqrt(k)/(GapThick()*alpha);                                       // formula (6) and (7)
+  Double_t dtdxc = kk*(k*(cosfd*cosf-cost*sinfd*sinf)-(alpha*mu/(betaM*betaM))*sint*sinfd);    // formula (6)           
+  Double_t dtdyc = kk*(k*(cosfd*sinf+cost*sinfd*cosf)+(alpha* e/(betaM*betaM))*sint*sinfd);    // formula (7)            pag.4
+  
+  Double_t errX = 0.2,errY=0.25;                                                                //end of page 7
+  return  TMath::Sqrt(errX*errX*dtdxc*dtdxc + errY*errY*dtdyc*dtdyc); 
+}
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+Double_t AliHMPIDPIDResponse::SigCrom(Double_t trkTheta,Double_t thetaC, Double_t phiC,Double_t betaM) const {
+
+// Analitical calculation of chromatic error (due to lack of knowledge of Cerenkov photon energy) on Cerenkov angle for a given Cerenkov photon 
+// created by a given MIP. Fromulae according to CERN-EP-2000-058 
+// Arguments: Cerenkov and azimuthal angles for Cerenkov photon, [radians]
+//            dip and azimuthal angles for MIP taken at the entrance to radiator, [radians]        
+//            MIP beta
+//   Returns: absolute error on Cerenkov angle, [radians]    
+  
+  Double_t phiDelta = phiC;
+
+  Double_t sint     = TMath::Sin(trkTheta);
+  Double_t cost     = TMath::Cos(trkTheta);
+  Double_t cosfd    = TMath::Cos(phiDelta);
+  Double_t tantheta = TMath::Tan(thetaC);
+  
+  Double_t alpha =cost-tantheta*cosfd*sint;                                         // formula (11)
+  Double_t dtdn = cost*GetRefIdx()*betaM*betaM/(alpha*tantheta);                    // formula (12)
+            
+  Double_t f = 0.0172*(7.75-5.635)/TMath::Sqrt(24.);
+
+  return f*dtdn;
+}//SigCrom()
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+Double_t AliHMPIDPIDResponse::SigGeom(Double_t trkTheta,Double_t trkPhi,Double_t thetaC, Double_t phiC,Double_t betaM) const {
+  
+// Analitical calculation of geometric error (due to lack of knowledge of creation point in radiator) on Cerenkov angle for a given Cerenkov photon 
+// created by a given MIP. Formulae according to CERN-EP-2000-058 
+// Arguments: Cerenkov and azimuthal angles for Cerenkov photon, [radians]
+//            dip and azimuthal angles for MIP taken at the entrance to radiator, [radians]        
+//            MIP beta
+//   Returns: absolute error on Cerenkov angle, [radians]    
+
+  Double_t phiDelta = phiC;
+
+  Double_t sint     = TMath::Sin(trkTheta);
+  Double_t cost     = TMath::Cos(trkTheta);
+  Double_t sinf     = TMath::Sin(trkPhi);
+  Double_t cosfd    = TMath::Cos(phiDelta);
+  Double_t costheta = TMath::Cos(thetaC);
+  Double_t tantheta = TMath::Tan(thetaC);
+  
+  Double_t alpha =cost-tantheta*cosfd*sint;                                                // formula (11)
+  
+  Double_t k = 1.-GetRefIdx()*GetRefIdx()+alpha*alpha/(betaM*betaM);                       // formula (after 8 in the text)
+  if (k<0) return 1e10;
+
+  Double_t eTr = 0.5*RadThick()*betaM*TMath::Sqrt(k)/(GapThick()*alpha);                    // formula (14)
+  Double_t lambda = (1.-sint*sinf)*(1.+sint*sinf);                                          // formula (15)
+
+  Double_t c1 = 1./(1.+ eTr*k/(alpha*alpha*costheta*costheta));                             // formula (13.a)
+  Double_t c2 = betaM*TMath::Power(k,1.5)*tantheta*lambda/(GapThick()*alpha*alpha);         // formula (13.b)
+  Double_t c3 = (1.+eTr*k*betaM*betaM)/((1+eTr)*alpha*alpha);                               // formula (13.c)
+  Double_t c4 = TMath::Sqrt(k)*tantheta*(1-lambda)/(GapThick()*betaM);                      // formula (13.d)
+  Double_t dtdT = c1 * (c2+c3*c4);
+  Double_t trErr = RadThick()/(TMath::Sqrt(12.)*cost);
+
+  return trErr*dtdT;
+}//SigGeom()
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+Double_t AliHMPIDPIDResponse::GetNMean(const AliVTrack *vTrk) const {
+
+  // 
+  // mean refractive index calculation
+  //
+  Double_t nmean = -999.; 
+  
+  Float_t xPc=0.,yPc=0.,thRa=0.,phRa=0.;
+  vTrk->GetHMPIDtrk(xPc,yPc,thRa,phRa);
+  
+  const Int_t ch = vTrk->GetHMPIDcluIdx()/1000000;
+  
+  const Double_t yRa = yPc - (RadThick()+WinThick()+GapThick())*TMath::Sin(phRa)*TMath::Tan(thRa);  //just linear extrapolation back to RAD
+      
+  TF1 *RefIndex=0x0;
+  
+  if(GetRefIndexArray()) RefIndex = (TF1*)(GetRefIndexArray()->At(ch));
+  else return nmean;
+  
+  if(RefIndex) nmean = RefIndex->Eval(yRa);
+  else return nmean;
+  
+  return nmean;   
+}  
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+Double_t AliHMPIDPIDResponse::SigmaCorrFact  (Int_t iPart, Double_t occupancy) {
+
+// calculation of sigma correction factor
+  
+  Double_t corr = 1.0;
+       
+  switch(iPart) {
+    case 0: corr = 0.115*occupancy + 1.166; break; 
+    case 1: corr = 0.115*occupancy + 1.166; break;
+    case 2: corr = 0.115*occupancy + 1.166; break;
+    case 3: corr = 0.065*occupancy + 1.137; break;
+    case 4: corr = 0.048*occupancy + 1.202; break;
+  }                                                                                                                           
+ return corr; 
+}
+
diff --git a/STEER/STEERBase/AliHMPIDPIDResponse.h b/STEER/STEERBase/AliHMPIDPIDResponse.h
new file mode 100644 (file)
index 0000000..566118e
--- /dev/null
@@ -0,0 +1,96 @@
+#ifndef ALIHMPIDPIDRESPONSE_H
+#define ALIHMPIDPIDRESPONSE_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+//***********************************************************
+// Class AliHMPIDPIDResponse
+//
+// HMPID class to perfom particle identification
+// 
+// Author: G. Volpe, giacomo.volpe@cern.ch
+//***********************************************************
+
+
+#include <TNamed.h>          //base class
+#include <TVector3.h>
+#include <TVector2.h>
+
+#include "AliPID.h"
+        
+class AliVTrack;
+class TObjArray;
+class TGeoHMatrix;
+
+class AliHMPIDPIDResponse : public TNamed 
+{
+public : 
+             AliHMPIDPIDResponse();    //ctor
+             AliHMPIDPIDResponse(const AliHMPIDPIDResponse& c);                //copy constructor
+             AliHMPIDPIDResponse &operator=(const AliHMPIDPIDResponse& c);     //dummy assignment operator
+    virtual ~AliHMPIDPIDResponse() {;} //dtor
+    
+    enum EChamberData{kMinCh=0,kMaxCh=6,kMinPc=0,kMaxPc=5};      //Segmenation
+    enum EPadxData{kPadPcX=80,kMinPx=0,kMaxPx=79,kMaxPcx=159};   //Segmentation structure along x
+    enum EPadyData{kPadPcY=48,kMinPy=0,kMaxPy=47,kMaxPcy=143};   //Segmentation structure along y 
+    
+    Double_t GetExpectedSignal  (const AliVTrack *vTrk, AliPID::EParticleType specie                 ) const;
+    Double_t GetExpectedSigma   (const AliVTrack *vTrk, AliPID::EParticleType specie                 ) const;                                                                     //Find the sigma for a given ThetaCerTh
+    Double_t GetNumberOfSigmas  (const AliVTrack *vTrk, AliPID::EParticleType specie                 ) const;                                                                     //Find the expected Cherenkov angle for a given track
+    void     GetProbability     (const AliVTrack *vTrk, Int_t nSpecies,Double_t *prob                ) const;                                                                     //Find the PID probability array
+    Double_t GetSignalDelta     (const AliVTrack *vTrk, AliPID::EParticleType specie                 ) const;    
+    void     Propagate          (const TVector3  dir,   TVector3 &pos,  Double_t z                   ) const;                                                                     //propagate photon alogn the line  
+    void     Refract            (TVector3 &dir,         Double_t n1,    Double_t n2                  ) const;                                                                     //refract photon on the boundary
+    TVector2 TracePhot          (Double_t xRa, Double_t yRa,  Double_t thRa, Double_t phRa, Double_t ckovThe,Double_t ckovPhi) const;                                             //trace photon created by track to PC 
+    void     Trs2Lors           (Double_t thRa, Double_t phRa, TVector3 dirCkov,      Double_t &thetaCer,Double_t &phiCer) const;                                                 //TRS to LORS
+    TVector2 TraceForward       (Double_t xRa, Double_t yRa, TVector3 dirCkov                        ) const;                                                                     //tracing forward a photon from (x,y) to PC
+    void     SetTrack           (Double_t xRad,         Double_t yRad,  Double_t theta,Double_t phi  )       {fTrkDir.SetMagThetaPhi(1,theta,phi);  fTrkPos.Set(xRad,yRad);}      //set track parameter at RAD
+    Double_t RadThick           (                                                                    ) const {return 1.5;}                                                        //Radiator thickness
+    Double_t WinThick           (                                                                    ) const {return 0.5;}                                                        //Window thickness
+    Double_t GapThick           (                                                                    ) const {return 8.0;}                                                        //Proximity gap thicknes
+    Double_t GetRefIdx          (                                                                    ) const {return fRefIdx;}                                                    //running refractive index
+    Double_t WinIdx             (                                                                    ) const {return 1.5787;}                                                     //Mean refractive index of WIN material (SiO2) 
+    Double_t GapIdx             (                                                                    ) const {return 1.0005;}                                                     //Mean refractive index of GAP material (CH4)
+    static Bool_t  IsInside     (Float_t x,Float_t y,Float_t d=0                                     )       {return  x>-d&&y>-d&&x<fgkMaxPcX[kMaxPc]+d&&y<fgkMaxPcY[kMaxPc]+d; } //is point inside chamber boundaries?
+    static Bool_t  IsInDead     (Float_t x,Float_t y                                                 );                                                                           //is the point in a dead area?
+    static Float_t SizeAllX     (                                                                    )       {return fgAllX;}                                                     //all PCs size x, [cm]        
+    static Float_t SizeAllY     (                                                                    )       {return fgAllY;}                                                     //all PCs size y, [cm]    
+    static void    IdealPosition(Int_t iCh,TGeoHMatrix *m                                            );                                                                           //ideal position of given chamber 
+        
+    Double_t SigLoc             (Double_t trkTheta,Double_t trkPhi,Double_t ckovTh,Double_t ckovPh,Double_t beta) const;                                                          //error due to cathode segmetation
+    Double_t SigGeom            (Double_t trkTheta,Double_t trkPhi,Double_t ckovTh,Double_t ckovPh,Double_t beta) const;                                                          //error due to unknown photon origin
+    Double_t SigCrom            (Double_t trkTheta,Double_t ckovTh,Double_t ckovPh,Double_t beta                ) const;                                                          //error due to unknonw photon energy
+    Double_t Sigma2             (Double_t trkTheta,Double_t trkPhi,Double_t ckovTh,Double_t ckovPh              ) const;                                                          //photon candidate sigma^2  
+    Double_t GetNMean           (const AliVTrack *vTrk                                                          ) const;
+    static Double_t SigmaCorrFact(Int_t iPart, Double_t occupancy                                               )      ;                                                          //correction factor for theoretical resolution
+    
+    void SetRefIndexArray       (TObjArray *array                                                     )       {fRefIndexArray = array;}
+    TObjArray* GetRefIndexArray (                                                                     ) const {return fRefIndexArray;}
+
+//
+private:
+        
+    Double_t ExpectedSignal     (const AliVTrack *vTrk, Double_t nmean, AliPID::EParticleType specie ) const;
+    Double_t ExpectedSigma      (const AliVTrack *vTrk, Double_t nmean, AliPID::EParticleType specie ) const;                                                                     //Find the sigma for a given ThetaCerTh
+    
+protected:
+
+  static /*const*/ Float_t fgkMinPcX[6];                               //limits PC
+  static /*const*/ Float_t fgkMinPcY[6];                               //limits PC
+  static /*const*/ Float_t fgkMaxPcX[6];                               //limits PC
+  static /*const*/ Float_t fgkMaxPcY[6]; 
+          
+  static Float_t fgCellX, fgCellY, fgPcX, fgPcY, fgAllX, fgAllY;       //definition of HMPID geometric parameters 
+        
+  TGeoHMatrix *fM[7];                                                  //pointers to matrices defining HMPID chambers rotations-translations
+  
+  Double_t  fRefIdx;                                                   //running refractive index of C6F14
+  TVector3  fTrkDir;                                                   //track direction in LORS at RAD
+  TVector2  fTrkPos;                                                   //track positon in LORS at RAD
+  TObjArray *fRefIndexArray;                                           //array of refracive index funxtion;
+  
+  ClassDef(AliHMPIDPIDResponse,1)
+};
+#endif // #ifdef AliHMPIDPIDResponse_cxx
+
index e3828fb..52430c6 100644 (file)
@@ -22,6 +22,7 @@
 //      Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
 //-----------------------------------------------------------------
 #include "TMath.h"
+#include "AliVTrack.h"
 #include "AliITSPIDResponse.h"
 #include "AliITSPidParams.h"
 #include "AliExternalTrackParam.h"
@@ -273,6 +274,52 @@ void AliITSPIDResponse::GetITSProbabilities(Float_t mom, Double_t qclu[4], Doubl
 }
 
 //_________________________________________________________________________
+Double_t AliITSPIDResponse::GetNumberOfSigmas( const AliVTrack* track, AliPID::EParticleType type) const
+{
+  //
+  // number of sigmas
+  //
+  UChar_t clumap=track->GetITSClusterMap();
+  Int_t nPointsForPid=0;
+  for(Int_t i=2; i<6; i++){
+    if(clumap&(1<<i)) ++nPointsForPid;
+  }
+  Float_t mom=track->P();
+  
+  //check for ITS standalone tracks
+  Bool_t isSA=kTRUE;
+  if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE;
+  
+  const Float_t dEdx=track->GetITSsignal();
+  
+  //TODO: in case of the electron, use the SA parametrisation,
+  //      this needs to be changed if ITS provides a parametrisation
+  //      for electrons also for ITS+TPC tracks
+  return GetNumberOfSigmas(mom,dEdx,type,nPointsForPid,isSA || (type==AliPID::kElectron));
+}
+
+//_________________________________________________________________________
+Double_t AliITSPIDResponse::GetSignalDelta( const AliVTrack* track, AliPID::EParticleType type) const
+{
+  //
+  // Signal - expected
+  //
+  const Float_t mom=track->P();
+  const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(type),2.);
+  Bool_t isSA=kTRUE;
+  if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE;
+  
+  const Float_t dEdx=track->GetITSsignal();
+  
+  //TODO: in case of the electron, use the SA parametrisation,
+  //      this needs to be changed if ITS provides a parametrisation
+  //      for electrons also for ITS+TPC tracks
+  
+  Float_t bethe = Bethe(mom,AliPID::ParticleMassZ(type), isSA || (type==AliPID::kElectron))*chargeFactor;
+  return (dEdx - bethe);
+}
+
+//_________________________________________________________________________
 Int_t AliITSPIDResponse::GetParticleIdFromdEdxVsP(Float_t mom, Float_t signal, Bool_t isSA) const{
   // method to get particle identity with simple cuts on dE/dx vs. momentum
 
index 325ee66..6dd6543 100644 (file)
@@ -14,6 +14,8 @@
 #include <TObject.h>
 #include "AliPID.h"
 
+class AliVTrack;
+
 class AliITSPIDResponse : public TObject {
 
 public:
@@ -40,6 +42,11 @@ public:
  Double_t BetheITSsaHybrid(Double_t p, Double_t mass) const;
  Double_t GetResolution(Double_t bethe, Int_t nPtsForPid=4, Bool_t isSA=kFALSE) const;
  void GetITSProbabilities(Float_t mom, Double_t qclu[4], Double_t condprobfun[AliPID::kSPECIES],Bool_t isMC=kFALSE) const;
+
+ Double_t GetNumberOfSigmas( const AliVTrack* track, AliPID::EParticleType species) const;
+
+ Double_t GetSignalDelta( const AliVTrack* track, AliPID::EParticleType species) const;
  Float_t GetNumberOfSigmas(Float_t mom, Float_t signal, AliPID::EParticleType type, Int_t nPtsForPid=4, Bool_t isSA=kFALSE) const {
    const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(type),2.);
    Float_t bethe = Bethe(mom,AliPID::ParticleMassZ(type),isSA)*chargeFactor;
index 5da7852..6f71e1c 100644 (file)
@@ -41,6 +41,7 @@
 #include <AliOADBContainer.h>
 #include <AliTRDPIDResponseObject.h>
 #include <AliTOFPIDParams.h>
+#include <AliHMPIDPIDParams.h>
 
 #include "AliPIDResponse.h"
 #include "AliDetectorPID.h"
@@ -55,6 +56,7 @@ fITSResponse(isMC),
 fTPCResponse(),
 fTRDResponse(),
 fTOFResponse(),
+fHMPIDResponse(),
 fEMCALResponse(),
 fRange(5.),
 fITSPIDmethod(kITSTruncMean),
@@ -81,6 +83,7 @@ fUseTPCEtaCorrection(kFALSE),//TODO: In future, default kTRUE
 fTRDPIDResponseObject(NULL),
 fTOFtail(1.1),
 fTOFPIDParams(NULL),
+fHMPIDPIDParams(NULL),
 fEMCALPIDParams(NULL),
 fCurrentEvent(NULL),
 fCurrCentrality(0.0),
@@ -113,6 +116,7 @@ fITSResponse(other.fITSResponse),
 fTPCResponse(other.fTPCResponse),
 fTRDResponse(other.fTRDResponse),
 fTOFResponse(other.fTOFResponse),
+fHMPIDResponse(other.fHMPIDResponse),
 fEMCALResponse(other.fEMCALResponse),
 fRange(other.fRange),
 fITSPIDmethod(other.fITSPIDmethod),
@@ -139,6 +143,7 @@ fUseTPCEtaCorrection(other.fUseTPCEtaCorrection),
 fTRDPIDResponseObject(NULL),
 fTOFtail(1.1),
 fTOFPIDParams(NULL),
+fHMPIDPIDParams(NULL),
 fEMCALPIDParams(NULL),
 fCurrentEvent(NULL),
 fCurrCentrality(0.0),
@@ -162,6 +167,7 @@ AliPIDResponse& AliPIDResponse::operator=(const AliPIDResponse &other)
     fTPCResponse=other.fTPCResponse;
     fTRDResponse=other.fTRDResponse;
     fTOFResponse=other.fTOFResponse;
+    fHMPIDResponse=other.fHMPIDResponse;
     fEMCALResponse=other.fEMCALResponse;
     fRange=other.fRange;
     fITSPIDmethod=other.fITSPIDmethod;
@@ -189,6 +195,7 @@ AliPIDResponse& AliPIDResponse::operator=(const AliPIDResponse &other)
     fEMCALPIDParams=NULL;
     fTOFtail=1.1;
     fTOFPIDParams=NULL;
+    fHMPIDPIDParams=NULL;
     fCurrentEvent=other.fCurrentEvent;
 
   }
@@ -278,6 +285,16 @@ Float_t AliPIDResponse::NumberOfSigmasTOF(const AliVParticle *vtrack, AliPID::EP
 }
 
 //______________________________________________________________________________
+Float_t AliPIDResponse::NumberOfSigmasHMPID(const AliVParticle *vtrack, AliPID::EParticleType type) const
+{
+  //
+  // Calculate the number of sigmas in the EMCAL
+  //
+  
+  return NumberOfSigmas(kHMPID, vtrack, type);
+}
+
+//______________________________________________________________________________
 Float_t AliPIDResponse::NumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type) const
 {
   //
@@ -357,6 +374,35 @@ Float_t  AliPIDResponse::NumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID:
 }
 
 //______________________________________________________________________________
+AliPIDResponse::EDetPidStatus AliPIDResponse::GetSignalDelta(EDetector detector, const AliVParticle *track, AliPID::EParticleType type, Double_t &val) const
+{
+  //
+  //
+  //
+  val=-9999.;
+  switch (detector){
+    case kITS:   return GetSignalDeltaITS(track,type,val); break;
+    case kTPC:   return GetSignalDeltaTPC(track,type,val); break;
+    case kTOF:   return GetSignalDeltaTOF(track,type,val); break;
+    case kHMPID: return GetSignalDeltaHMPID(track,type,val); break;
+    default: return kDetNoSignal;
+  }
+  return kDetNoSignal;
+}
+
+//______________________________________________________________________________
+Double_t AliPIDResponse::GetSignalDelta(EDetector detCode, const AliVParticle *track, AliPID::EParticleType type) const
+{
+  //
+  //
+  //
+  Double_t val=-9999.;
+  EDetPidStatus stat=GetSignalDelta(detCode, track, type, val);
+  if ( stat==kDetNoSignal ) val=-9999.;
+  return val;
+}
+
+//______________________________________________________________________________
 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability  (EDetCode  detCode, const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
 {
   // Compute PID response of 'detCode'
@@ -532,6 +578,9 @@ void AliPIDResponse::ExecNewRun()
   SetTOFPidResponseMaster();
   InitializeTOFResponse();
 
+  SetHMPIDPidResponseMaster();
+  InitializeHMPIDResponse();
+
   if (fCurrentEvent) fTPCResponse.SetMagField(fCurrentEvent->GetMagneticField());
 }
 
@@ -1322,6 +1371,37 @@ void AliPIDResponse::InitializeTOFResponse(){
 
 }
 
+//______________________________________________________________________________
+void AliPIDResponse::SetHMPIDPidResponseMaster()
+{
+  //
+  // Load the HMPID pid params from the OADB
+  //
+
+  if (fHMPIDPIDParams) delete fHMPIDPIDParams;
+  fHMPIDPIDParams=NULL;
+
+  TFile *oadbf = new TFile(Form("%s/COMMON/PID/data/HMPIDPIDParams.root",fOADBPath.Data()));
+  if (oadbf && oadbf->IsOpen()) {
+    AliInfo(Form("Loading HMPID Params from %s/COMMON/PID/data/HMPIDPIDParams.root", fOADBPath.Data()));
+    AliOADBContainer *oadbc = (AliOADBContainer *)oadbf->Get("HMPoadb");
+    if (oadbc) fHMPIDPIDParams = dynamic_cast<AliHMPIDPIDParams *>(oadbc->GetObject(fRun,"HMPparams"));
+    oadbf->Close();
+    delete oadbc;
+  }
+  delete oadbf;
+
+  if (!fHMPIDPIDParams) AliFatal("HMPIDPIDParams could not be retrieved");
+}
+
+//______________________________________________________________________________
+void AliPIDResponse::InitializeHMPIDResponse(){
+  //
+  // Set PID Params to the HMPID PID response
+  //
+
+  fHMPIDResponse.SetRefIndexArray(fHMPIDPIDParams->GetHMPIDrefIndex());
+}
 
 //______________________________________________________________________________
 Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack, Double_t efficiencyLevel,Double_t centrality,AliTRDPIDResponse::ETRDPIDMethod PIDmethod) const {
@@ -1681,9 +1761,10 @@ Float_t AliPIDResponse::GetNumberOfSigmas(EDetector detector, const AliVParticle
   const AliVTrack *track=static_cast<const AliVTrack*>(vtrack);
   
   switch (detector){
-    case kITS:   return GetNumberOfSigmasITS(track, type); break;
-    case kTPC:   return GetNumberOfSigmasTPC(track, type); break;
-    case kTOF:   return GetNumberOfSigmasTOF(track, type); break;
+    case kITS:   return GetNumberOfSigmasITS(track, type);   break;
+    case kTPC:   return GetNumberOfSigmasTPC(track, type);   break;
+    case kTOF:   return GetNumberOfSigmasTOF(track, type);   break;
+    case kHMPID: return GetNumberOfSigmasHMPID(track, type); break;
     case kEMCAL: return GetNumberOfSigmasEMCAL(track, type); break;
     default: return -999.;
   }
@@ -1702,24 +1783,8 @@ Float_t AliPIDResponse::GetNumberOfSigmasITS(const AliVParticle *vtrack, AliPID:
 
   const EDetPidStatus pidStatus=GetITSPIDStatus(track);
   if (pidStatus!=kDetPidOk) return -999.;
-    
-  UChar_t clumap=track->GetITSClusterMap();
-  Int_t nPointsForPid=0;
-  for(Int_t i=2; i<6; i++){
-    if(clumap&(1<<i)) ++nPointsForPid;
-  }
-  Float_t mom=track->P();
-  
-  //check for ITS standalone tracks
-  Bool_t isSA=kTRUE;
-  if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE;
-  
-  const Float_t dEdx=track->GetITSsignal();
 
-  //TODO: in case of the electron, use the SA parametrisation,
-  //      this needs to be changed if ITS provides a parametrisation
-  //      for electrons also for ITS+TPC tracks
-  return fITSResponse.GetNumberOfSigmas(mom,dEdx,type,nPointsForPid,isSA || (type==AliPID::kElectron));
+  return fITSResponse.GetNumberOfSigmas(track,type);
 }
 
 //______________________________________________________________________________
@@ -1734,14 +1799,7 @@ Float_t AliPIDResponse::GetNumberOfSigmasTPC(const AliVParticle *vtrack, AliPID:
   const EDetPidStatus pidStatus=GetTPCPIDStatus(track);
   if (pidStatus!=kDetPidOk) return -999.;
   
-  Double_t nSigma = -999.;
-  
-  if (fTuneMConData)
-    this->GetTPCsignalTunedOnData(track);
-  
-  nSigma = fTPCResponse.GetNumberOfSigmas(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
-  
-  return nSigma;
+  return fTPCResponse.GetNumberOfSigmas(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
 }
 
 //______________________________________________________________________________
@@ -1755,10 +1813,23 @@ Float_t AliPIDResponse::GetNumberOfSigmasTOF(const AliVParticle *vtrack, AliPID:
 
   const EDetPidStatus pidStatus=GetTOFPIDStatus(track);
   if (pidStatus!=kDetPidOk) return -999.;
-
   
   return GetNumberOfSigmasTOFold(vtrack, type);
 }
+//______________________________________________________________________________
+
+Float_t AliPIDResponse::GetNumberOfSigmasHMPID(const AliVParticle *vtrack, AliPID::EParticleType type) const
+{
+  //
+  // Calculate the number of sigmas in the HMPID
+  //  
+  AliVTrack *track=(AliVTrack*)vtrack;
+    
+  const EDetPidStatus pidStatus=GetHMPIDPIDStatus(track);
+  if (pidStatus!=kDetPidOk) return -999.; 
+  
+  return fHMPIDResponse.GetNumberOfSigmas(track, type);
+}
 
 //______________________________________________________________________________
 Float_t AliPIDResponse::GetNumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type) const
@@ -1784,6 +1855,53 @@ Float_t AliPIDResponse::GetNumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPI
   return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge);
 }
 
+//______________________________________________________________________________
+AliPIDResponse::EDetPidStatus AliPIDResponse::GetSignalDeltaITS(const AliVParticle *vtrack, AliPID::EParticleType type, Double_t &val) const
+{
+  //
+  // Signal minus expected Signal for ITS
+  //
+  AliVTrack *track=(AliVTrack*)vtrack;
+  val=fITSResponse.GetSignalDelta(track,type);
+  
+  return GetITSPIDStatus(track);
+}
+
+//______________________________________________________________________________
+AliPIDResponse::EDetPidStatus AliPIDResponse::GetSignalDeltaTPC(const AliVParticle *vtrack, AliPID::EParticleType type, Double_t &val) const
+{
+  //
+  // Signal minus expected Signal for TPC
+  //
+  AliVTrack *track=(AliVTrack*)vtrack;
+  val=fTPCResponse.GetSignalDelta(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
+  
+  return GetTPCPIDStatus(track);
+}
+
+//______________________________________________________________________________
+AliPIDResponse::EDetPidStatus AliPIDResponse::GetSignalDeltaTOF(const AliVParticle *vtrack, AliPID::EParticleType type, Double_t &val) const
+{
+  //
+  // Signal minus expected Signal for TOF
+  //
+  AliVTrack *track=(AliVTrack*)vtrack;
+  val=GetSignalDeltaTOFold(track, type);
+  
+  return GetTOFPIDStatus(track);
+}
+
+//______________________________________________________________________________
+AliPIDResponse::EDetPidStatus AliPIDResponse::GetSignalDeltaHMPID(const AliVParticle *vtrack, AliPID::EParticleType type, Double_t &val) const
+{
+  //
+  // Signal minus expected Signal for HMPID
+  //
+  AliVTrack *track=(AliVTrack*)vtrack;
+  val=fHMPIDResponse.GetSignalDelta(track, type);
+  
+  return GetHMPIDPIDStatus(track);
+}
 
 //______________________________________________________________________________
 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputePIDProbability  (EDetector detCode,  const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
@@ -2016,8 +2134,8 @@ AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeHMPIDProbability(const A
   const EDetPidStatus pidStatus=GetHMPIDPIDStatus(track);
   if (pidStatus!=kDetPidOk) return pidStatus;
   
-  track->GetHMPIDpid(p);
-  
+  fHMPIDResponse.GetProbability(track,nSpecies,p);
+    
   return kDetPidOk;
 }
 
@@ -2118,13 +2236,16 @@ Float_t AliPIDResponse::GetTOFMismatchProbability(const AliVTrack *track) const
   return 0.;
 }
 
-
-
 //______________________________________________________________________________
 AliPIDResponse::EDetPidStatus AliPIDResponse:: GetHMPIDPIDStatus(const AliVTrack *track) const
 {
   // compute HMPID pid status
-  if((track->GetStatus()&AliVTrack::kHMPIDpid)==0) return kDetNoSignal;
+  
+  Int_t ch = track->GetHMPIDcluIdx()/1000000;
+  Double_t HMPIDsignal = track->GetHMPIDsignal(); 
+  
+  if((track->GetStatus()&AliVTrack::kHMPIDpid)==0 || ch<0 || ch>6 || HMPIDsignal<0) return kDetNoSignal;
+  
   return kDetPidOk;
 }
 
index 1ecf0aa..85897d7 100644 (file)
@@ -15,6 +15,7 @@
 #include "AliTPCPIDResponse.h"
 #include "AliTRDPIDResponse.h"
 #include "AliTOFPIDResponse.h"
+#include "AliHMPIDPIDResponse.h"
 #include "AliEMCALPIDResponse.h"
 
 
@@ -30,6 +31,7 @@ class TLinearFitter;
 class AliVEvent;
 class AliTRDPIDResponseObject;
 class AliTOFPIDParams;
+class AliHMPIDPIDParams;
 class AliOADBContainer;
 
 class AliPIDResponse : public TNamed {
@@ -90,10 +92,15 @@ public:
   virtual Float_t NumberOfSigmasEMCAL(const AliVParticle *track, AliPID::EParticleType type, Double_t &eop, Double_t showershape[4]) const;
   virtual Float_t NumberOfSigmasTOF  (const AliVParticle *track, AliPID::EParticleType type) const;
   virtual Float_t NumberOfSigmasTOF  (const AliVParticle *track, AliPID::EParticleType type, const Float_t /*timeZeroTOF*/) const { return NumberOfSigmasTOF(track,type); }
+  virtual Float_t NumberOfSigmasHMPID(const AliVParticle *track, AliPID::EParticleType type) const;
   virtual Float_t NumberOfSigmasEMCAL(const AliVParticle *track, AliPID::EParticleType type) const;
 
   Bool_t IdentifiedAsElectronTRD(const AliVTrack *track, Double_t efficiencyLevel,Double_t centrality=-1,AliTRDPIDResponse::ETRDPIDMethod PIDmethod=AliTRDPIDResponse::kLQ1D) const;
 
+  // Signal delta
+  EDetPidStatus GetSignalDelta(EDetector detCode, const AliVParticle *track, AliPID::EParticleType type, Double_t &val) const;
+  Double_t GetSignalDelta(EDetector detCode, const AliVParticle *track, AliPID::EParticleType type) const;
+  
   // Probabilities
   EDetPidStatus ComputePIDProbability  (EDetCode  detCode, const AliVTrack *track, Int_t nSpecies, Double_t p[]) const;
   EDetPidStatus ComputePIDProbability  (EDetector detCode, const AliVTrack *track, Int_t nSpecies, Double_t p[]) const;
@@ -156,10 +163,11 @@ public:
 
   
 protected:
-  AliITSPIDResponse fITSResponse;    //PID response function of the ITS
-  AliTPCPIDResponse fTPCResponse;    //PID response function of the TPC
-  AliTRDPIDResponse fTRDResponse;    //PID response function of the TRD
-  AliTOFPIDResponse fTOFResponse;    //PID response function of the TOF
+  AliITSPIDResponse   fITSResponse;    //PID response function of the ITS
+  AliTPCPIDResponse   fTPCResponse;    //PID response function of the TPC
+  AliTRDPIDResponse   fTRDResponse;    //PID response function of the TRD
+  AliTOFPIDResponse   fTOFResponse;    //PID response function of the TOF
+  AliHMPIDPIDResponse fHMPIDResponse;  //PID response function of the HMPID
   AliEMCALPIDResponse fEMCALResponse;  //PID response function of the EMCAL
 
   Float_t           fRange;          // nSigma max in likelihood
@@ -167,6 +175,8 @@ protected:
 
   //unbuffered PID calculation
   virtual Float_t GetNumberOfSigmasTOFold  (const AliVParticle */*track*/, AliPID::EParticleType /*type*/) const {return 0;}
+  virtual Float_t GetSignalDeltaTOFold(const AliVParticle */*track*/, AliPID::EParticleType /*type*/) const {return -9999.;}
+  
   EDetPidStatus GetComputeTRDProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[],AliTRDPIDResponse::ETRDPIDMethod PIDmethod=AliTRDPIDResponse::kLQ1D) const;
   EDetPidStatus GetTOFPIDStatus(const AliVTrack *track) const;
 
@@ -199,6 +209,8 @@ private:
 
   Float_t fTOFtail;                    //! TOF tail effect used in TOF probability
   AliTOFPIDParams *fTOFPIDParams;      //! TOF PID Params - period depending (OADB loaded)
+  
+  AliHMPIDPIDParams *fHMPIDPIDParams;  //! HMPID PID Params (OADB loaded)
 
   TObjArray *fEMCALPIDParams;             //! EMCAL PID Params
 
@@ -237,6 +249,10 @@ private:
   void SetTOFPidResponseMaster();
   void InitializeTOFResponse();
 
+  //HMPID
+  void SetHMPIDPidResponseMaster();
+  void InitializeHMPIDResponse();
+
   //EMCAL
   void SetEMCALPidResponseMaster();
   void InitializeEMCALResponse();
@@ -253,10 +269,17 @@ private:
   Float_t GetNumberOfSigmasITS  (const AliVParticle *track, AliPID::EParticleType type) const;
   Float_t GetNumberOfSigmasTPC  (const AliVParticle *track, AliPID::EParticleType type) const;
   Float_t GetNumberOfSigmasTOF  (const AliVParticle *track, AliPID::EParticleType type) const;
+  Float_t GetNumberOfSigmasHMPID(const AliVParticle *track, AliPID::EParticleType type) const;
   Float_t GetNumberOfSigmasEMCAL(const AliVParticle *track, AliPID::EParticleType type, Double_t &eop, Double_t showershape[4]) const;
   Float_t GetNumberOfSigmasEMCAL(const AliVParticle *track, AliPID::EParticleType type) const;
 
   Float_t GetBufferedNumberOfSigmas(EDetector detCode, const AliVParticle *track, AliPID::EParticleType type) const;
+
+  // Signal deltas
+  EDetPidStatus GetSignalDeltaITS(const AliVParticle *track, AliPID::EParticleType type, Double_t &val) const;
+  EDetPidStatus GetSignalDeltaTPC(const AliVParticle *track, AliPID::EParticleType type, Double_t &val) const;
+  EDetPidStatus GetSignalDeltaTOF(const AliVParticle *track, AliPID::EParticleType type, Double_t &val) const;
+  EDetPidStatus GetSignalDeltaHMPID(const AliVParticle *vtrack, AliPID::EParticleType type, Double_t &val) const;
   
   // Probabilities
   EDetPidStatus GetComputePIDProbability  (EDetector detCode,  const AliVTrack *track, Int_t nSpecies, Double_t p[]) const;
index ec59c6e..2de82b0 100644 (file)
@@ -514,7 +514,7 @@ Float_t AliTPCPIDResponse::GetNumberOfSigmas(const AliVTrack* track,
   //Calculates the number of sigmas of the PID signal from the expected value
   //for a given particle species in the presence of multiple gain scenarios
   //inside the TPC
-                             
+  
   Double_t dEdx = -1;
   Int_t nPoints = -1;
   ETPCgainScenario gainScenario = kGainScenarioInvalid;
@@ -522,7 +522,7 @@ Float_t AliTPCPIDResponse::GetNumberOfSigmas(const AliVTrack* track,
   
   if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
     return -999; //TODO: Better handling!
-                             
+    
   Double_t bethe = GetExpectedSignal(track, species, dEdx, responseFunction, correctEta);
   Double_t sigma = GetExpectedSigma(track, species, gainScenario, dEdx, nPoints, responseFunction, correctEta);
   // 999 will be returned by GetExpectedSigma e.g. in case of 0 dEdx clusters
@@ -533,6 +533,29 @@ Float_t AliTPCPIDResponse::GetNumberOfSigmas(const AliVTrack* track,
 }
 
 //_________________________________________________________________________
+Float_t AliTPCPIDResponse::GetSignalDelta(const AliVTrack* track,
+                                          AliPID::EParticleType species,
+                                          ETPCdEdxSource dedxSource,
+                                          Bool_t correctEta) const
+{
+  //Calculates the number of sigmas of the PID signal from the expected value
+  //for a given particle species in the presence of multiple gain scenarios
+  //inside the TPC
+
+  Double_t dEdx = -1;
+  Int_t nPoints = -1;
+  ETPCgainScenario gainScenario = kGainScenarioInvalid;
+  TSpline3* responseFunction = 0x0;
+
+  if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
+    return -9999.; //TODO: Better handling!
+
+  Double_t bethe = GetExpectedSignal(track, species, dEdx, responseFunction, correctEta);
+  // 999 will be returned by GetExpectedSigma e.g. in case of 0 dEdx clusters
+  return dEdx-bethe;
+}
+
+//_________________________________________________________________________
 Bool_t AliTPCPIDResponse::ResponseFunctiondEdxN( const AliVTrack* track, 
                                                  AliPID::EParticleType species,
                                                  ETPCdEdxSource dedxSource,
index 1816551..b78c751 100644 (file)
@@ -113,7 +113,12 @@ public:
                              AliPID::EParticleType species,
                              ETPCdEdxSource dedxSource = kdEdxDefault,
                              Bool_t correctEta = kFALSE) const;//TODO: In future, default kTRUE
-
+  
+  Float_t GetSignalDelta( const AliVTrack* track,
+                          AliPID::EParticleType species,
+                          ETPCdEdxSource dedxSource = kdEdxDefault,
+                          Bool_t correctEta = kFALSE) const;
+  
   void SetResponseFunction(TObject* o,
                            AliPID::EParticleType type,
                            ETPCgainScenario gainScenario);
index 19ec538..c954912 100644 (file)
@@ -94,6 +94,15 @@ public:
   virtual Double_t  GetHMPIDsignal()     const {return 0.;}
   virtual Double_t  GetTRDsignal()       const {return 0.;}
 
+  virtual Double_t  GetHMPIDoccupancy()  const {return 0.;}
+  
+  virtual Int_t     GetHMPIDcluIdx()     const {return 0;}
+  
+  virtual void GetHMPIDtrk(Float_t &/*&x*/, Float_t &/*y*/, Float_t &/*th*/, Float_t &/*ph*/) const {;}  
+  virtual void GetHMPIDmip(Float_t &/*x*/, Float_t &/*y*/, Int_t &/*q*/,Int_t &/*nph*/) const {;}
+  
+  virtual Bool_t GetOuterHmpPxPyPz(Double_t */*p*/) const {return kFALSE;}
+  
   virtual void      GetIntegratedTimes(Double_t */*times*/) const { return; }
   virtual Double_t  GetTRDmomentum(Int_t /*plane*/, Double_t */*sp*/=0x0) const {return 0.;}
   virtual void      GetHMPIDpid(Double_t */*p*/) const {;}
index 8346351..2865bbd 100644 (file)
@@ -99,6 +99,8 @@
 #pragma link C++ class AliTOFPIDResponse+;
 #pragma link C++ class AliTRDPIDResponse+;
 #pragma link C++ class AliEMCALPIDResponse+;
+#pragma link C++ class AliHMPIDPIDResponse+;
+#pragma link C++ class AliHMPIDPIDParams+;
 #pragma link C++ class AliPIDCombined+;
 #pragma link C++ class AliPIDValues+;
 #pragma link C++ class AliDetectorPID+;