o First Version of TRDnSigma implementation (Xianguo) o still requires some catching...

diff --git a/STEER/STEERBase/AliPIDResponse.cxx b/STEER/STEERBase/AliPIDResponse.cxx
index 8f5c167347d9ba471fccc0b4b60cd975b7dc83b4..84f1cc57e23282e4d38d661050ee8107f0db6a0c 100644 (file)
--- a/STEER/STEERBase/AliPIDResponse.cxx
+++ b/STEER/STEERBase/AliPIDResponse.cxx
@@ -290,6 +290,15 @@ Float_t AliPIDResponse::NumberOfSigmasTPC( const AliVParticle *vtrack,
   return nSigma;
 }
 
+//______________________________________________________________________________
+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;
@@ -2100,6 +2111,21 @@ Float_t AliPIDResponse::GetNumberOfSigmasTPC(const AliVParticle *vtrack, AliPID:
   return fTPCResponse.GetNumberOfSigmas(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection, fUseTPCMultiplicityCorrection);
 }
 
+//______________________________________________________________________________
+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
 {
@@ -2184,6 +2210,18 @@ AliPIDResponse::EDetPidStatus AliPIDResponse::GetSignalDeltaTPC(const AliVPartic
   return GetTPCPIDStatus(track);
 }
 
+//______________________________________________________________________________
+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
 {

diff --git a/STEER/STEERBase/AliPIDResponse.h b/STEER/STEERBase/AliPIDResponse.h
index ac7ae843710dd9c4d950bbe9342da5ffc0be108a..280c68e72f437c72da080d02313b53a70bf08756 100644 (file)
--- a/STEER/STEERBase/AliPIDResponse.h
+++ b/STEER/STEERBase/AliPIDResponse.h
@@ -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;
   

diff --git a/STEER/STEERBase/AliTRDPIDResponse.cxx b/STEER/STEERBase/AliTRDPIDResponse.cxx
index 0db79d9ea8786ef18d529e3583a6298ac9a9bed4..8b6cfbf7e44924861d504a4be39b31c42b5c9af2 100644 (file)
--- a/STEER/STEERBase/AliTRDPIDResponse.cxx
+++ b/STEER/STEERBase/AliTRDPIDResponse.cxx
@@ -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
 {

diff --git a/STEER/STEERBase/AliTRDPIDResponse.h b/STEER/STEERBase/AliTRDPIDResponse.h
index b9c1db6e6d43e5d71f72a0bab75634b5bcf3ca58..d8c156fe7cdbad393ed99aaee7c28d3ab6eed8e2 100644 (file)
--- a/STEER/STEERBase/AliTRDPIDResponse.h
+++ b/STEER/STEERBase/AliTRDPIDResponse.h
@@ -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;
     

diff --git a/STEER/STEERBase/AliVTrack.h b/STEER/STEERBase/AliVTrack.h
index b928d06855855d33d2ee77feab3ce326820adb02..99f0df918d1238bc1f90bd11a67fd70d257f356c 100644 (file)
--- a/STEER/STEERBase/AliVTrack.h
+++ b/STEER/STEERBase/AliVTrack.h
@@ -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.;}