o First Version of TRDnSigma implementation (Xianguo) o still requires some catching...
authorwiechula <Jens.Wiechula@cern.ch>
Thu, 3 Jul 2014 11:18:07 +0000 (13:18 +0200)
committermorsch <andreas.morsch@cern.ch>
Fri, 4 Jul 2014 20:44:52 +0000 (22:44 +0200)
STEER/STEERBase/AliPIDResponse.cxx
STEER/STEERBase/AliPIDResponse.h
STEER/STEERBase/AliTRDPIDResponse.cxx
STEER/STEERBase/AliTRDPIDResponse.h
STEER/STEERBase/AliVTrack.h

index 8f5c167..84f1cc5 100644 (file)
@@ -291,6 +291,15 @@ Float_t AliPIDResponse::NumberOfSigmasTPC( const AliVParticle *vtrack,
 }
 
 //______________________________________________________________________________
+Float_t AliPIDResponse::NumberOfSigmasTRD(const AliVParticle *vtrack, AliPID::EParticleType type) const
+{
+  //
+  // Calculate the number of sigmas in the TRD
+  //
+  return NumberOfSigmas(kTRD, vtrack, type);
+}
+
+//______________________________________________________________________________
 Float_t AliPIDResponse::NumberOfSigmasTOF(const AliVParticle *vtrack, AliPID::EParticleType type) const
 {
   //
@@ -399,6 +408,7 @@ AliPIDResponse::EDetPidStatus AliPIDResponse::GetSignalDelta(EDetector detector,
   switch (detector){
     case kITS:   return GetSignalDeltaITS(track,type,val,ratio); break;
     case kTPC:   return GetSignalDeltaTPC(track,type,val,ratio); break;
+    case kTRD:   return GetSignalDeltaTRD(track,type,val,ratio); break;
     case kTOF:   return GetSignalDeltaTOF(track,type,val,ratio); break;
     case kHMPID: return GetSignalDeltaHMPID(track,type,val,ratio); break;
     default: return kDetNoSignal;
@@ -2055,6 +2065,7 @@ Float_t AliPIDResponse::GetNumberOfSigmas(EDetector detector, const AliVParticle
   switch (detector){
     case kITS:   return GetNumberOfSigmasITS(track, type);   break;
     case kTPC:   return GetNumberOfSigmasTPC(track, type);   break;
+    case kTRD:   return GetNumberOfSigmasTRD(track, type);   break;
     case kTOF:   return GetNumberOfSigmasTOF(track, type);   break;
     case kHMPID: return GetNumberOfSigmasHMPID(track, type); break;
     case kEMCAL: return GetNumberOfSigmasEMCAL(track, type); break;
@@ -2101,6 +2112,21 @@ Float_t AliPIDResponse::GetNumberOfSigmasTPC(const AliVParticle *vtrack, AliPID:
 }
 
 //______________________________________________________________________________
+Float_t AliPIDResponse::GetNumberOfSigmasTRD(const AliVParticle *vtrack, AliPID::EParticleType type) const
+{
+  //
+  // Calculate the number of sigmas in the TRD
+  //
+  
+  AliVTrack *track=(AliVTrack*)vtrack;
+  
+  const EDetPidStatus pidStatus=GetTRDPIDStatus(track);
+  if (pidStatus!=kDetPidOk) return -999.;
+  
+  return fTRDResponse.GetNumberOfSigmas(track,type);
+}
+
+//______________________________________________________________________________
 Float_t AliPIDResponse::GetNumberOfSigmasTOF(const AliVParticle *vtrack, AliPID::EParticleType type) const
 {
   //
@@ -2185,6 +2211,18 @@ AliPIDResponse::EDetPidStatus AliPIDResponse::GetSignalDeltaTPC(const AliVPartic
 }
 
 //______________________________________________________________________________
+AliPIDResponse::EDetPidStatus AliPIDResponse::GetSignalDeltaTRD(const AliVParticle *vtrack, AliPID::EParticleType type, Double_t &val, Bool_t ratio/*=kFALSE*/) const
+{
+  //
+  // Signal minus expected Signal for TRD
+  //
+  AliVTrack *track=(AliVTrack*)vtrack;
+  val=fTRDResponse.GetSignalDelta(track,type,ratio);
+  
+  return GetTRDPIDStatus(track);
+}
+
+//______________________________________________________________________________
 AliPIDResponse::EDetPidStatus AliPIDResponse::GetSignalDeltaTOF(const AliVParticle *vtrack, AliPID::EParticleType type, Double_t &val, Bool_t ratio/*=kFALSE*/) const
 {
   //
index ac7ae84..280c68e 100644 (file)
@@ -95,6 +95,7 @@ public:
   virtual Float_t NumberOfSigmasITS  (const AliVParticle *track, AliPID::EParticleType type) const;
   virtual Float_t NumberOfSigmasTPC  (const AliVParticle *track, AliPID::EParticleType type) const;
   virtual Float_t NumberOfSigmasTPC  (const AliVParticle *track, AliPID::EParticleType type, AliTPCPIDResponse::ETPCdEdxSource dedxSource) const;
+  virtual Float_t NumberOfSigmasTRD  (const AliVParticle *track, AliPID::EParticleType type) const;
   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, Float_t /*timeZeroTOF*/) const { return NumberOfSigmasTOF(track,type); }
@@ -298,6 +299,7 @@ private:
   Float_t GetNumberOfSigmas(EDetector detCode, const AliVParticle *track, AliPID::EParticleType type) const;
   Float_t GetNumberOfSigmasITS  (const AliVParticle *track, AliPID::EParticleType type) const;
   Float_t GetNumberOfSigmasTPC  (const AliVParticle *track, AliPID::EParticleType type) const;
+  Float_t GetNumberOfSigmasTRD  (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;
@@ -308,6 +310,7 @@ private:
   // Signal deltas
   EDetPidStatus GetSignalDeltaITS(const AliVParticle *track, AliPID::EParticleType type, Double_t &val, Bool_t ratio=kFALSE) const;
   EDetPidStatus GetSignalDeltaTPC(const AliVParticle *track, AliPID::EParticleType type, Double_t &val, Bool_t ratio=kFALSE) const;
+  EDetPidStatus GetSignalDeltaTRD(const AliVParticle *track, AliPID::EParticleType type, Double_t &val, Bool_t ratio=kFALSE) const;
   EDetPidStatus GetSignalDeltaTOF(const AliVParticle *track, AliPID::EParticleType type, Double_t &val, Bool_t ratio=kFALSE) const;
   EDetPidStatus GetSignalDeltaHMPID(const AliVParticle *vtrack, AliPID::EParticleType type, Double_t &val, Bool_t ratio=kFALSE) const;
   
index 0db79d9..8b6cfbf 100644 (file)
@@ -35,6 +35,8 @@
 #include <TSystem.h>
 
 #include "AliLog.h"
+#include "AliVTrack.h"
 
 #include "AliTRDPIDResponseObject.h"
 #include "AliTRDPIDResponse.h"
@@ -113,6 +115,185 @@ Bool_t AliTRDPIDResponse::Load(const Char_t * filename){
   return kTRUE;
 }
 
+
+
+//____________________________________________________________
+Double_t AliTRDPIDResponse::GetNumberOfSigmas(const AliVTrack *track, AliPID::EParticleType type) const
+{
+  //
+  //calculate the TRD nSigma
+  //
+
+  const Double_t badval = -9999;
+  Double_t info[5]; for(int i=0; i<5; i++){info[i]=badval;}
+
+  const Double_t delta = GetSignalDelta(track, type, kFALSE, info);
+
+  const Double_t mean = info[0];
+  const Double_t res = info[1];
+  if(res<0){
+    return badval;
+  }
+
+  const Double_t sigma = mean*res;
+  return delta/sigma;
+}
+
+//____________________________________________________________
+Double_t AliTRDPIDResponse::GetSignalDelta( const AliVTrack* track, AliPID::EParticleType type, Bool_t ratio/*=kFALSE*/, Double_t *info/*=0x0*/) const
+{
+  //
+  //calculate the TRD signal difference w.r.t. the expected
+  //output other information in info[]
+  //
+
+  //============================================================================================>>>>>>>>>>>>
+  //Data production dependent, in the end should go to OADB
+  //============================================================================================
+  //~/.task/TRDdEdx/2013output/resolution$ r bgmean.root
+  // hfit;1  1.521e+00 , 7.294e-01 , 8.066e+00 , 2.146e-02 , 3.137e+01 , 7.160e-15 , 2.838e+00 , 1.100e+01 , (14290.45/67) conv_runlist_LHC13b.pass3.root trdncls/trdnch>=18 && trdnch>=6 && tpcncls>120 && pidV0>=1          
+  // hscn0_wid;1     -4.122e-01 , 1.019e+00 , -1.319e-01 , (15071.70/120) conv_runlist_LHC13b.pass3.root tpcncls>120 && pidV0>=1
+
+  //mean signal parametrization
+  Double_t meanpars[]={1.521e+00 , 7.294e-01 , 8.066e+00 , 2.146e-02 , 3.137e+01 , 7.160e-15 , 2.838e+00 , 1.100e+01};
+  const Int_t nmp = sizeof(meanpars)/sizeof(Double_t);
+
+  if(type==AliPID::kProton){
+    const Double_t tmpproton[]={2.026e+00 , -1.462e-04 , 1.202e+00 , 1.162e-01 , 2.092e+00 , -3.018e-02 , 3.665e+00 , 3.487e-07};
+    for(Int_t ii=0; ii<nmp; ii++){
+      meanpars[ii] = tmpproton[ii];
+    }
+  }
+  else if(type==AliPID::kPion || type==AliPID::kElectron){
+    const Double_t tmppe[]={1.508e+00 , 7.356e-01 , 8.002e+00 , 1.932e-02 , 3.058e+01 , 5.114e-18 , 3.561e+00 , 1.408e+01};
+    for(Int_t ii=0; ii<nmp; ii++){
+      meanpars[ii] = tmppe[ii];
+    }
+  }
+
+  //resolution parametrization
+  const Double_t respars[]={-4.122e-01 , 1.019e+00 , -1.319e-01};
+
+  //cut on number of chambers
+  const Double_t cutch = 6;
+
+  //cut on mean number of clusters per chamber
+  const Double_t cutclsperch = 18;
+  //============================================================================================<<<<<<<<<<<<<
+
+
+  const Double_t badval = -9999;
+
+  const Double_t nch = track->GetTRDNchamberdEdx();
+  if(nch < cutch){
+    return badval;
+  }
+
+  const Double_t ncls = track->GetTRDNclusterdEdx();
+  if( ncls/nch < cutclsperch){
+    return badval;
+  }
+
+  Double_t pTRD = -999;
+  for(Int_t ich=0; ich<6; ich++){
+    pTRD = track->GetTRDmomentum(ich);
+    if(pTRD>0)
+      break;
+  }
+  if(pTRD<0){
+    return badval;
+  }
+
+  const Double_t bg = pTRD/AliPID::ParticleMass(type);
+  const Double_t expsig = MeandEdxTR(&bg, meanpars);
+
+  if(info){
+    info[0]= expsig;
+    info[1]= ResolutiondEdxTR(&ncls, respars);
+  }
+
+  if(ratio){
+    return track->GetTRDsignal()/expsig;
+  }
+  else{
+    return track->GetTRDsignal() - expsig;
+  }
+}
+
+
+Double_t AliTRDPIDResponse::ResolutiondEdxTR(const Double_t * xx,  const Double_t * par)
+{
+  //
+  //resolution 
+  //npar=3
+  //
+
+  const Double_t ncls = xx[0];
+
+  return par[0]+par[1]*TMath::Power(ncls, par[2]);
+}
+
+Double_t AliTRDPIDResponse::MeandEdxTR(const Double_t * xx,  const Double_t * pin)
+{
+  //
+  //ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
+  //npar = 8 = 3+5
+  //
+
+  Double_t ptr[4]={0};
+  for(int ii=0; ii<3; ii++){
+    ptr[ii+1]=pin[ii];
+  }
+  return MeanTR(xx,ptr) + MeandEdx(xx,&(pin[3]));
+}
+
+Double_t AliTRDPIDResponse::MeanTR(const Double_t * xx,  const Double_t * par)
+{
+  //
+  //LOGISTIC parametrization for TR, in unit of MIP
+  //npar = 4
+  //
+                                                                                                                                                                                                                                                        
+  const Double_t bg = xx[0];
+  const Double_t gamma = sqrt(1+bg*bg);
+
+  const Double_t p0 = TMath::Abs(par[1]);
+  const Double_t p1 = TMath::Abs(par[2]);
+  const Double_t p2 = TMath::Abs(par[3]);
+
+  const Double_t zz = TMath::Log(gamma);
+  const Double_t tryield = p0/( 1 + TMath::Exp(-p1*(zz-p2)) );
+
+  return par[0]+tryield;
+}
+
+Double_t AliTRDPIDResponse::MeandEdx(const Double_t * xx,  const Double_t * par)
+{
+  //
+  //ALEPH parametrization for dEdx
+  //npar = 5
+  //
+
+  const Double_t bg = xx[0];
+  const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
+
+  const Double_t p0 = TMath::Abs(par[0]);
+  const Double_t p1 = TMath::Abs(par[1]);
+
+  const Double_t p2 = TMath::Abs(par[2]);
+
+  const Double_t p3 = TMath::Abs(par[3]);
+  const Double_t p4 = TMath::Abs(par[4]);
+
+  const Double_t aa = TMath::Power(beta, p3);
+
+  const Double_t bb = TMath::Log( p2 + TMath::Power(1./bg, p4) );
+
+  return (p1-aa-bb)*p0/aa;
+
+}
+
+
 //____________________________________________________________
 Int_t AliTRDPIDResponse::GetResponse(Int_t n, const Double_t * const dedx, const Float_t * const p, Double_t prob[AliPID::kSPECIES],ETRDPIDMethod PIDmethod,Bool_t kNorm) const
 {
index b9c1db6..d8c156f 100644 (file)
@@ -59,6 +59,13 @@ class AliTRDPIDResponse : public TObject {
     AliTRDPIDResponse& operator=(const AliTRDPIDResponse &ref);
     ~AliTRDPIDResponse();
     
+    Double_t GetNumberOfSigmas(const AliVTrack *track, AliPID::EParticleType type) const;
+    Double_t GetSignalDelta( const AliVTrack* track, AliPID::EParticleType type, Bool_t ratio=kFALSE, Double_t *info=0x0) const;
+    static Double_t MeandEdx(const Double_t * xx, const Double_t * par);
+    static Double_t MeanTR(const Double_t * xx, const Double_t * par);
+    static Double_t MeandEdxTR(const Double_t * xx, const Double_t * par);
+    static Double_t ResolutiondEdxTR(const Double_t * xx,  const Double_t * par);
+
     Int_t    GetResponse(Int_t n, const Double_t * const dedx, const Float_t * const p, Double_t prob[AliPID::kSPECIES],ETRDPIDMethod PIDmethod=kLQ1D, Bool_t kNorm=kTRUE) const;
     inline ETRDNslices  GetNumberOfSlices(ETRDPIDMethod PIDmethod=kLQ1D) const;
     
index b928d06..99f0df9 100644 (file)
@@ -121,6 +121,8 @@ public:
   virtual Double_t  GetTOFsignalTunedOnData() const {return 0.;}
   virtual Double_t  GetHMPIDsignal()     const {return 0.;}
   virtual Double_t  GetTRDsignal()       const {return 0.;}
+  virtual UChar_t GetTRDNchamberdEdx() const {return 0;}
+  virtual UChar_t GetTRDNclusterdEdx() const {return 0;}
 
   virtual Double_t  GetHMPIDoccupancy()  const {return 0.;}