X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TPC%2FAliTPCClusterParam.cxx;h=703ffb0d05b7f5da8fa8f167ae259a78baefcf81;hb=b9908d0b8e9e96b9f6cf84fd98a9907fd2a4e5c6;hp=9ff841554a49381fa28acd75d7d7cf0d44fd527a;hpb=db2fdcfb8dd86cb3f07f3cbd56451d08216f76bd;p=u%2Fmrichter%2FAliRoot.git

diff --git a/TPC/AliTPCClusterParam.cxx b/TPC/AliTPCClusterParam.cxx
index 9ff841554a4..703ffb0d05b 100644
--- a/TPC/AliTPCClusterParam.cxx
+++ b/TPC/AliTPCClusterParam.cxx
@@ -61,7 +61,7 @@
 //
 //  Example how to retrieve the paramterization:
 /*    
-      AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT");
+      AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
       AliCDBManager::Instance()->SetRun(0) 
       AliTPCClusterParam * param = AliTPCcalibDB::Instance()->GetClusterParam();
 
@@ -91,10 +91,14 @@ AliTPCClusterParam::SetInstance(param);
 #include <TVectorF.h>
 #include <TLinearFitter.h>
 #include <TH1F.h>
+#include <TH3F.h>
 #include <TProfile2D.h>
 #include <TVectorD.h>
 #include <TObjArray.h>
 #include "AliTPCcalibDB.h"
+#include "AliTPCParam.h"
+
+#include "AliMathBase.h"
 
 ClassImp(AliTPCClusterParam)
 
@@ -146,6 +150,8 @@ AliTPCClusterParam::AliTPCClusterParam():
   TObject(),
   fRatio(0),
   fQNorm(0),
+  fQNormCorr(0),
+  fQNormHis(0),
   fQpadTnorm(0),           // q pad normalization - Total charge
   fQpadMnorm(0)            // q pad normalization - Max charge
   //
@@ -164,6 +170,8 @@ AliTPCClusterParam::AliTPCClusterParam(const AliTPCClusterParam& param):
   TObject(param),
   fRatio(0),
   fQNorm(0),
+  fQNormCorr(0),
+  fQNormHis(0),
   fQpadTnorm(new TVectorD(*(param.fQpadTnorm))),           // q pad normalization - Total charge
   fQpadMnorm(new TVectorD(*(param.fQpadMnorm)))            // q pad normalization - Max charge
 
@@ -173,6 +181,7 @@ AliTPCClusterParam::AliTPCClusterParam(const AliTPCClusterParam& param):
   //
   memcpy(this, &param,sizeof(AliTPCClusterParam));
   if (param.fQNorm) fQNorm = (TObjArray*) param.fQNorm->Clone();
+  if (param.fQNormHis) fQNormHis = (TObjArray*) param.fQNormHis->Clone();
   //
   if (param.fPosQTnorm[0]){
     fPosQTnorm[0] = new TVectorD(*(param.fPosQTnorm[0]));
@@ -203,6 +212,7 @@ AliTPCClusterParam & AliTPCClusterParam::operator=(const AliTPCClusterParam& par
   if (this != &param) {
     memcpy(this, &param,sizeof(AliTPCClusterParam));
     if (param.fQNorm) fQNorm = (TObjArray*) param.fQNorm->Clone();
+    if (param.fQNormHis) fQNormHis = (TObjArray*) param.fQNormHis->Clone();
     if (param.fPosQTnorm[0]){
       fPosQTnorm[0] = new TVectorD(*(param.fPosQTnorm[0]));
       fPosQTnorm[1] = new TVectorD(*(param.fPosQTnorm[1]));
@@ -231,7 +241,10 @@ AliTPCClusterParam::~AliTPCClusterParam(){
   // destructor
   //
   if (fQNorm) fQNorm->Delete();
+  if (fQNormCorr) delete fQNormCorr;
+  if (fQNormHis) fQNormHis->Delete();
   delete fQNorm;
+  delete fQNormHis;
   if (fPosQTnorm[0]){
     delete fPosQTnorm[0];
     delete fPosQTnorm[1];
@@ -779,7 +792,7 @@ void AliTPCClusterParam::FitRMSSigma(TTree * tree, Int_t dim, Int_t type, Float_
 
 
 
-Float_t  AliTPCClusterParam::GetError0(Int_t dim, Int_t type, Float_t z, Float_t angle){
+Float_t  AliTPCClusterParam::GetError0(Int_t dim, Int_t type, Float_t z, Float_t angle) const {
   //
   //
   //
@@ -792,7 +805,7 @@ Float_t  AliTPCClusterParam::GetError0(Int_t dim, Int_t type, Float_t z, Float_t
 }
 
 
-Float_t  AliTPCClusterParam::GetError0Par(Int_t dim, Int_t type, Float_t z, Float_t angle){
+Float_t  AliTPCClusterParam::GetError0Par(Int_t dim, Int_t type, Float_t z, Float_t angle) const {
   //
   //
   //
@@ -809,7 +822,7 @@ Float_t  AliTPCClusterParam::GetError0Par(Int_t dim, Int_t type, Float_t z, Floa
 
 
 
-Float_t  AliTPCClusterParam::GetError1(Int_t dim, Int_t type, Float_t z, Float_t angle){
+Float_t  AliTPCClusterParam::GetError1(Int_t dim, Int_t type, Float_t z, Float_t angle) const {
   //
   //
   //
@@ -824,7 +837,7 @@ Float_t  AliTPCClusterParam::GetError1(Int_t dim, Int_t type, Float_t z, Float_t
   return value;
 }
 
-Float_t  AliTPCClusterParam::GetErrorQ(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){
+Float_t  AliTPCClusterParam::GetErrorQ(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean) const {
   //
   //
   //
@@ -840,7 +853,7 @@ Float_t  AliTPCClusterParam::GetErrorQ(Int_t dim, Int_t type, Float_t z, Float_t
 
 }
 
-Float_t  AliTPCClusterParam::GetErrorQPar(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){
+Float_t  AliTPCClusterParam::GetErrorQPar(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean) const {
   //
   //
   //
@@ -859,7 +872,7 @@ Float_t  AliTPCClusterParam::GetErrorQPar(Int_t dim, Int_t type, Float_t z, Floa
 
 }
 
-Float_t  AliTPCClusterParam::GetErrorQParScaled(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){
+Float_t  AliTPCClusterParam::GetErrorQParScaled(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean) const {
   //
   //
   //
@@ -880,7 +893,7 @@ Float_t  AliTPCClusterParam::GetErrorQParScaled(Int_t dim, Int_t type, Float_t z
 
 }
 
-Float_t  AliTPCClusterParam::GetRMS0(Int_t dim, Int_t type, Float_t z, Float_t angle){
+Float_t  AliTPCClusterParam::GetRMS0(Int_t dim, Int_t type, Float_t z, Float_t angle) const {
   //
   // calculate mean RMS of cluster - z,angle - parameters for each pad and dimension separatelly
   //
@@ -892,7 +905,7 @@ Float_t  AliTPCClusterParam::GetRMS0(Int_t dim, Int_t type, Float_t z, Float_t a
   return value;
 }
 
-Float_t  AliTPCClusterParam::GetRMS1(Int_t dim, Int_t type, Float_t z, Float_t angle){
+Float_t  AliTPCClusterParam::GetRMS1(Int_t dim, Int_t type, Float_t z, Float_t angle) const {
   //
   // calculate mean RMS of cluster - z,angle - pad length scalling
   //
@@ -911,7 +924,7 @@ Float_t  AliTPCClusterParam::GetRMS1(Int_t dim, Int_t type, Float_t z, Float_t a
   return value;
 }
 
-Float_t  AliTPCClusterParam::GetRMSQ(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){
+Float_t  AliTPCClusterParam::GetRMSQ(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean) const {
   //
   // calculate mean RMS of cluster - z,angle, Q dependence
   //
@@ -925,7 +938,7 @@ Float_t  AliTPCClusterParam::GetRMSQ(Int_t dim, Int_t type, Float_t z, Float_t a
   return value;
 }
 
-Float_t  AliTPCClusterParam::GetRMSSigma(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){
+Float_t  AliTPCClusterParam::GetRMSSigma(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean) const {
   //
   // calculates RMS of signal shape fluctuation
   //
@@ -935,7 +948,7 @@ Float_t  AliTPCClusterParam::GetRMSSigma(Int_t dim, Int_t type, Float_t z, Float
   return value;
 }
 
-Float_t  AliTPCClusterParam::GetShapeFactor(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean, Float_t rmsL, Float_t rmsM){
+Float_t  AliTPCClusterParam::GetShapeFactor(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean, Float_t rmsL, Float_t rmsM) const {
   //
   // calculates vallue - sigma distortion contribution
   //
@@ -1317,7 +1330,7 @@ Float_t AliTPCClusterParam::Qnorm(Int_t ipad, Int_t itype, Float_t dr, Float_t t
   // type - 0 Qtot 1 Qmax
   // ipad - 0 (0.75 cm) ,1 (1 cm), 2 (1.5 cm)
   //
-  //expession formula - TString *strq0 = toolkit.FitPlane(chain,"dedxQ.fElements[2]","dr++ty++tz++dr*ty++dr*tz++ty*tz++ty^2++tz^2","IPad==0",chi2,npoints,param,covar,0,100000);
+  //expession formula - TString *strq0 = toolkit.FitPlane(chain,"dedxQ.fElements[2]","dr++ty++tz++dr*ty++dr*tz++++dr*dr++ty*tz++ty^2++tz^2","IPad==0",chi2,npoints,param,covar,0,100000);
 
   if (fQNorm==0) return 0;
   TVectorD * norm = (TVectorD*)fQNorm->At(3*itype+ipad);
@@ -1340,7 +1353,34 @@ Float_t AliTPCClusterParam::Qnorm(Int_t ipad, Int_t itype, Float_t dr, Float_t t
 
 
 
-void AliTPCClusterParam::SetQnorm(Int_t ipad, Int_t itype, TVectorD * norm){
+Float_t AliTPCClusterParam::QnormHis(Int_t ipad, Int_t itype, Float_t dr, Float_t p2, Float_t p3){
+  // get Q normalization
+  // type - 0 Qtot 1 Qmax
+  // ipad - 0 (0.75 cm) ,1 (1 cm), 2 (1.5 cm)
+  //
+
+  if (fQNormHis==0) return 0;
+  TH3F * norm = (TH3F*)fQNormHis->At(4*itype+ipad);
+  if (!norm) return 1;
+  p2=TMath::Abs(p2);
+  dr=TMath::Min(dr,Float_t(norm->GetXaxis()->GetXmax()-norm->GetXaxis()->GetBinWidth(0)));
+  dr=TMath::Max(dr,Float_t(norm->GetXaxis()->GetXmin()+norm->GetXaxis()->GetBinWidth(0)));
+  //
+  p2=TMath::Min(p2,Float_t(norm->GetYaxis()->GetXmax()-norm->GetYaxis()->GetBinWidth(0)));
+  p2=TMath::Max(p2,Float_t(norm->GetYaxis()->GetXmin()+norm->GetYaxis()->GetBinWidth(0)));
+  //
+  p3=TMath::Min(p3,Float_t(norm->GetZaxis()->GetXmax()-norm->GetZaxis()->GetBinWidth(0)));
+  p3=TMath::Max(p3,Float_t(norm->GetZaxis()->GetXmin()+norm->GetZaxis()->GetBinWidth(0)));
+  //
+  Double_t res = norm->GetBinContent(norm->FindBin(dr,p2,p3));
+  if (res==0) res = norm->GetBinContent(norm->FindBin(0.5,0.5,0.5));  // This is just hack - to be fixed entries without 
+
+  return res;
+}
+
+
+
+void AliTPCClusterParam::SetQnorm(Int_t ipad, Int_t itype, const TVectorD *const norm){
   //
   // set normalization
   //
@@ -1352,6 +1392,50 @@ void AliTPCClusterParam::SetQnorm(Int_t ipad, Int_t itype, TVectorD * norm){
   fQNorm->AddAt(new TVectorD(*norm), itype*3+ipad);
 }
 
+void AliTPCClusterParam::ResetQnormCorr(){
+  //
+  //
+  //
+  if (!fQNormCorr) fQNormCorr= new TMatrixD(12,6);
+  for (Int_t irow=0;irow<12; irow++)
+    for (Int_t icol=0;icol<6; icol++){
+      (*fQNormCorr)(irow,icol)=1.;             // default - no correction
+      if (irow>5) (*fQNormCorr)(irow,icol)=0.; // default - no correction
+    } 
+}
+
+void AliTPCClusterParam::SetQnormCorr(Int_t ipad, Int_t itype, Int_t corrType, Float_t val){
+  //
+  // ipad        - pad type
+  // itype       - 0- qtot 1-qmax
+  // corrType    - 0 - s0y corr     - eff. PRF corr
+  //             - 1 - s0z corr     - eff. TRF corr
+  //             - 2 - d0y          - eff. diffusion correction y
+  //             - 3 - d0z          - eff. diffusion correction
+  //             - 4 - eff length   - eff.length - wire pitch + x diffsion
+  //             - 5 - pad type normalization
+  if (!fQNormCorr) {
+    ResetQnormCorr();
+  }
+  //
+  // eff shap parameterization matrix
+  //
+  // rows
+  // itype*3+ipad  - itype=0 qtot itype=1 qmax ipad=0
+  // 
+  if (itype<2) (*fQNormCorr)(itype*3+ipad, corrType) *= val;  // multiplicative correction
+  if (itype>=2) (*fQNormCorr)(itype*3+ipad, corrType)+= val;  // additive       correction  
+}
+
+Double_t  AliTPCClusterParam::GetQnormCorr(Int_t ipad, Int_t itype, Int_t corrType) const{
+  //
+  // see AliTPCClusterParam::SetQnormCorr
+  //
+  if (!fQNormCorr) return 0;
+  return  (*fQNormCorr)(itype*3+ipad, corrType);
+}
+
+
 Float_t AliTPCClusterParam::QnormPos(Int_t ipad,Bool_t isMax, Float_t pad, Float_t time, Float_t z, Float_t sy2, Float_t sz2, Float_t qm, Float_t qt){
   //
   // Make Q normalization as function of following parameters
@@ -1504,3 +1588,210 @@ Float_t AliTPCClusterParam::PosCorrection(Int_t type, Int_t ipad,  Float_t pad,
 
 
 
+Double_t  AliTPCClusterParam::GaussConvolution(Double_t x0, Double_t x1, Double_t k0, Double_t k1, Double_t s0, Double_t s1){
+  //
+  // 2 D gaus convoluted with angular effect
+  // See in mathematica: 
+  //Simplify[Integrate[Exp[-(x0-k0*xd)*(x0-k0*xd)/(2*s0*s0)-(x1-k1*xd)*(x1-k1*xd)/(2*s1*s1)]/(s0*s1),{xd,-1/2,1/2}]]
+  // 
+  //TF1 f1("f1","AliTPCClusterParam::GaussConvolution(x,0,1,0,0.1,0.1)",-2,2)
+  //TF2 f2("f2","AliTPCClusterParam::GaussConvolution(x,y,1,1,0.1,0.1)",-2,2,-2,2)
+  //
+  const Float_t kEpsilon = 0.0001;
+  if ((TMath::Abs(k0)+TMath::Abs(k1))<kEpsilon*(s0+s1)){
+    // small angular effect
+    Double_t val = (TMath::Gaus(x0,0,s0)*TMath::Gaus(x1,0,s1))/(s0*s1*2.*TMath::Pi());
+    return val;
+  }
+  Double_t sigma2 = k1*k1*s0*s0+k0*k0*s1*s1;
+  Double_t exp0 = TMath::Exp(-(k1*x0-k0*x1)*(k1*x0-k0*x1)/(2*sigma2));
+  //
+  Double_t sigmaErf =  2*s0*s1*TMath::Sqrt(2*sigma2);			     
+  Double_t erf0 = AliMathBase::ErfFast( (k0*s1*s1*(k0-2*x0)+k1*s0*s0*(k1-2*x1))/sigmaErf);
+  Double_t erf1 = AliMathBase::ErfFast( (k0*s1*s1*(k0+2*x0)+k1*s0*s0*(k1+2*x1))/sigmaErf);
+  Double_t norm = 1./TMath::Sqrt(sigma2);
+  norm/=2.*TMath::Sqrt(2.*TMath::Pi());
+  Double_t val  = norm*exp0*(erf0+erf1);
+  return val;
+}
+
+
+Double_t  AliTPCClusterParam::GaussConvolutionTail(Double_t x0, Double_t x1, Double_t k0, Double_t k1, Double_t s0, Double_t s1, Double_t tau){
+  //
+  // 2 D gaus convoluted with angular effect and exponential tail in z-direction
+  // tail integrated numerically 
+  // Integral normalized to one
+  // Mean at 0
+  // 
+  // TF1 f1t("f1t","AliTPCClusterParam::GaussConvolutionTail(0,x,0,0,0.5,0.5,0.9)",-5,5)
+  Double_t sum =1, mean=0;
+  // the COG of exponent
+  for (Float_t iexp=0;iexp<5;iexp+=0.2){
+    mean+=iexp*TMath::Exp(-iexp/tau);
+    sum +=TMath::Exp(-iexp/tau);
+  }
+  mean/=sum;
+  //
+  sum = 1;
+  Double_t val = GaussConvolution(x0,x1+mean, k0, k1 , s0,s1);
+  for (Float_t iexp=0;iexp<5;iexp+=0.2){
+    val+=GaussConvolution(x0,x1+mean-iexp, k0, k1 , s0,s1)*TMath::Exp(-iexp/tau);
+    sum+=TMath::Exp(-iexp/tau);
+  }
+  return val/sum;
+}
+
+Double_t  AliTPCClusterParam::GaussConvolutionGamma4(Double_t x0, Double_t x1, Double_t k0, Double_t k1, Double_t s0, Double_t s1, Double_t tau){
+  //
+  // 2 D gaus convoluted with angular effect and exponential tail in z-direction
+  // tail integrated numerically 
+  // Integral normalized to one
+  // Mean at 0
+  // 
+  // TF1 f1g4("f1g4","AliTPCClusterParam::GaussConvolutionGamma4(0,x,0,0,0.5,0.2,1.6)",-5,5)
+  // TF2 f2g4("f2g4","AliTPCClusterParam::GaussConvolutionGamma4(y,x,0,0,0.5,0.2,1.6)",-5,5,-5,5)
+  Double_t sum =0, mean=0;
+  // the COG of G4
+  for (Float_t iexp=0;iexp<5;iexp+=0.2){
+    Double_t g4 = TMath::Exp(-4.*iexp/tau)*TMath::Power(iexp/tau,4.);
+    mean+=iexp*g4;
+    sum +=g4;
+  }
+  mean/=sum;
+  //
+  sum = 0;
+  Double_t val = 0;
+  for (Float_t iexp=0;iexp<5;iexp+=0.2){ 
+    Double_t g4 = TMath::Exp(-4.*iexp/tau)*TMath::Power(iexp/tau,4.);
+    val+=GaussConvolution(x0,x1+mean-iexp, k0, k1 , s0,s1)*g4;
+    sum+=g4;
+  }
+  return val/sum;
+}
+
+Double_t  AliTPCClusterParam::QmaxCorrection(Int_t sector, Int_t row, Float_t cpad, Float_t ctime, Float_t ky, Float_t kz, Float_t rmsy0, Float_t rmsz0, Float_t effPad, Float_t effDiff){
+  //
+  //
+  // cpad      - pad (y) coordinate
+  // ctime     - time(z) coordinate
+  // ky        - dy/dx
+  // kz        - dz/dx
+  // rmsy0     - RF width in pad units
+  // rmsz0     - RF width in time bin  units
+  // effLength - contibution of PRF and diffusion
+  // effDiff   - overwrite diffusion
+
+  // Response function aproximated by convolution of gaussian with angular effect (integral=1)
+  //  
+  // Gaus width sy and sz is determined by RF width and diffusion 
+  // Integral of Q is equal 1
+  // Q max is calculated at position cpad, ctime
+  // Example function:         
+  //  TF1 f1("f1", "AliTPCClusterParam::QmaxCorrection(0,0.5,x,0,0,0.5,0.6)",0,1000) 
+  //
+  AliTPCParam * param   = AliTPCcalibDB::Instance()->GetParameters(); 
+  Double_t padLength= param->GetPadPitchLength(sector,row);
+  Double_t padWidth = param->GetPadPitchWidth(sector);
+  Double_t zwidth   = param->GetZWidth();
+  Double_t effLength= padLength+(param->GetWWPitch(0)+TMath::Sqrt(ctime*zwidth)*param->GetDiffT())*effPad;
+
+  // diffusion in pad, time bin  units
+  Double_t diffT=TMath::Sqrt(ctime*zwidth)*param->GetDiffT()/padWidth;
+  Double_t diffL=TMath::Sqrt(ctime*zwidth)*param->GetDiffL()/zwidth;
+  diffT*=effDiff;  //
+  diffL*=effDiff;  //
+  //
+  // transform angular effect to pad units
+  //
+  Double_t pky   = ky*effLength/padWidth;
+  Double_t pkz   = kz*effLength/zwidth;
+  // position in pad unit
+  Double_t py = (cpad+0.5)-TMath::Nint(cpad+0.5);
+  Double_t pz = (ctime+0.5)-TMath::Nint(ctime+0.5);
+  //
+  //
+  Double_t sy = TMath::Sqrt(rmsy0*rmsy0+diffT*diffT);
+  Double_t sz = TMath::Sqrt(rmsz0*rmsz0+diffL*diffL); 
+  //return GaussConvolutionGamma4(py,pz, pky,pkz,sy,sz,tau);
+  Double_t length = padLength*TMath::Sqrt(1+ky*ky+kz*kz);
+  return GaussConvolution(py,pz, pky,pkz,sy,sz)*length;
+}
+
+Double_t  AliTPCClusterParam::QtotCorrection(Int_t sector, Int_t row, Float_t cpad, Float_t ctime, Float_t ky, Float_t kz, Float_t rmsy0, Float_t rmsz0,  Float_t qtot, Float_t thr, Float_t effPad, Float_t effDiff){
+  //
+  //
+  // cpad      - pad (y) coordinate
+  // ctime     - time(z) coordinate
+  // ky        - dy/dx
+  // kz        - dz/dx
+  // rmsy0     - RF width in pad units
+  // rmsz0     - RF width in time bin  units
+  // qtot      - the sum of signal in cluster - without thr correction
+  // thr       - threshold
+  // effLength - contibution of PRF and diffusion
+  // effDiff   - overwrite diffusion
+
+  // Response function aproximated by convolution of gaussian with angular effect (integral=1)
+  //  
+  // Gaus width sy and sz is determined by RF width and diffusion 
+  // Integral of Q is equal 1
+  // Q max is calculated at position cpad, ctime
+  //          
+  //  
+  //
+  AliTPCParam * param   = AliTPCcalibDB::Instance()->GetParameters(); 
+  Double_t padLength= param->GetPadPitchLength(sector,row);
+  Double_t padWidth = param->GetPadPitchWidth(sector);
+  Double_t zwidth   = param->GetZWidth();
+  Double_t effLength= padLength+(param->GetWWPitch(0)+TMath::Sqrt(ctime*zwidth)*param->GetDiffT())*effPad;
+  //
+  // diffusion in pad units
+  Double_t diffT=TMath::Sqrt(ctime*zwidth)*param->GetDiffT()/padWidth;
+  Double_t diffL=TMath::Sqrt(ctime*zwidth)*param->GetDiffL()/zwidth;
+  diffT*=effDiff;  //
+  diffL*=effDiff;  //
+  //
+  // transform angular effect to pad units 
+  Double_t pky   = ky*effLength/padWidth;
+  Double_t pkz   = kz*effLength/zwidth;
+  // position in pad unit
+  //  
+  Double_t py = (cpad+0.5)-TMath::Nint(cpad+0.5);
+  Double_t pz = (ctime+0.5)-TMath::Nint(ctime+0.5);
+  //
+  Double_t sy = TMath::Sqrt(rmsy0*rmsy0+diffT*diffT);
+  Double_t sz = TMath::Sqrt(rmsz0*rmsz0+diffL*diffL); 
+  //
+  //
+  //
+  Double_t sumAll=0,sumThr=0;
+  //
+  Double_t corr =1;
+  Double_t qnorm=qtot;
+  for (Float_t iy=-3;iy<=3;iy+=1.)
+    for (Float_t iz=-4;iz<=4;iz+=1.){
+      //      Double_t val = GaussConvolutionGamma4(py-iy,pz-iz, pky,pkz, sy,sz,tau);      
+      Double_t val = GaussConvolution(py-iy,pz-iz, pky,pkz, sy,sz);      
+      Double_t qlocal =qnorm*val;
+      if (TMath::Abs(iy)<1.5&&TMath::Abs(iz)<1.5){
+      	sumThr+=qlocal;   // Virtual charge used in cluster finder
+      }
+      else{
+	if (qlocal>thr && TMath::Abs(iz)<2.5&&TMath::Abs(iy)<2.5) sumThr+=qlocal;
+      }
+      sumAll+=qlocal;
+    }
+  if (sumAll>0&&sumThr>0) {
+    corr=(sumThr)/sumAll;
+  }
+  //
+  Double_t length = padLength*TMath::Sqrt(1+ky*ky+kz*kz);
+  return corr*length;
+}
+
+
+
+
+
+
+