]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TPC/AliTPCClusterParam.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
bugfix
[u/mrichter/AliRoot.git] / TPC / AliTPCClusterParam.cxx
diff --git a/TPC/AliTPCClusterParam.cxx b/TPC/AliTPCClusterParam.cxx
index 9d023b70934801e7f662b5f5980970150e6c5b8b..1e652b98877b23ca7cfd4968da86f19166347c4b 100644 (file)
--- a/TPC/AliTPCClusterParam.cxx
+++ b/TPC/AliTPCClusterParam.cxx
@@ -97,6 +97,7 @@ AliTPCClusterParam::SetInstance(param);
 #include <TObjArray.h>
 #include "AliTPCcalibDB.h"
 #include "AliTPCParam.h"
+#include "THnBase.h"
 
 #include "AliMathBase.h"
 
@@ -153,7 +154,12 @@ AliTPCClusterParam::AliTPCClusterParam():
   fQNormCorr(0),
   fQNormHis(0),
   fQpadTnorm(0),           // q pad normalization - Total charge
-  fQpadMnorm(0)            // q pad normalization - Max charge
+  fQpadMnorm(0),           // q pad normalization - Max charge
+  fWaveCorrectionMap(0),
+  fWaveCorrectionMirroredPad( kFALSE ),
+  fWaveCorrectionMirroredZ( kFALSE ),
+  fWaveCorrectionMirroredAngle( kFALSE ),
+  fResolutionYMap(0)
   //
 {
   //
@@ -174,13 +180,16 @@ AliTPCClusterParam::AliTPCClusterParam(const AliTPCClusterParam& param):
   fQNormCorr(0),
   fQNormHis(0),
   fQpadTnorm(new TVectorD(*(param.fQpadTnorm))),           // q pad normalization - Total charge
-  fQpadMnorm(new TVectorD(*(param.fQpadMnorm)))            // q pad normalization - Max charge
-
+  fQpadMnorm(new TVectorD(*(param.fQpadMnorm))),           // q pad normalization - Max charge
+  fWaveCorrectionMap(0),
+  fWaveCorrectionMirroredPad( kFALSE ),
+  fWaveCorrectionMirroredZ( kFALSE ),
+  fWaveCorrectionMirroredAngle( kFALSE ),
+  fResolutionYMap(0)
 {
   //
   // copy constructor
   //
-  memcpy(this, &param,sizeof(AliTPCClusterParam));
   if (param.fQNorm) fQNorm = (TObjArray*) param.fQNorm->Clone();
   if (param.fQNormHis) fQNormHis = (TObjArray*) param.fQNormHis->Clone();
   //
@@ -202,7 +211,49 @@ AliTPCClusterParam::AliTPCClusterParam(const AliTPCClusterParam& param):
     fPosZcor[1] = new TVectorD(*(param.fPosZcor[1]));
     fPosZcor[2] = new TVectorD(*(param.fPosZcor[2]));
   }
-  
+
+  for (Int_t ii = 0; ii < 2; ++ii) {
+    for (Int_t jj = 0; jj < 3; ++jj) {
+      for (Int_t kk = 0; kk < 4; ++kk) {
+       fParamS0[ii][jj][kk] = param.fParamS0[ii][jj][kk];
+       fErrorS0[ii][jj][kk] = param.fErrorS0[ii][jj][kk];
+       fParamRMS0[ii][jj][kk] = param.fParamRMS0[ii][jj][kk];
+       fErrorRMS0[ii][jj][kk] = param.fErrorRMS0[ii][jj][kk];
+      }
+      for (Int_t kk = 0; kk < 7; ++kk) {
+       fParamS0Par[ii][jj][kk] = param.fParamS0Par[ii][jj][kk];
+       fErrorS0Par[ii][jj][kk] = param.fErrorS0Par[ii][jj][kk];
+      }
+      for (Int_t kk = 0; kk < 6; ++kk) {
+       fParamSQ[ii][jj][kk] = param.fParamSQ[ii][jj][kk];
+       fErrorSQ[ii][jj][kk] = param.fErrorSQ[ii][jj][kk];
+       fParamRMSQ[ii][jj][kk] = param.fParamRMSQ[ii][jj][kk];
+       fErrorRMSQ[ii][jj][kk] = param.fErrorRMSQ[ii][jj][kk];
+      }
+      for (Int_t kk = 0; kk < 9; ++kk) {
+       fParamSQPar[ii][jj][kk] = param.fParamSQPar[ii][jj][kk];
+       fErrorSQPar[ii][jj][kk] = param.fErrorSQPar[ii][jj][kk];
+      }
+      for (Int_t kk = 0; kk < 2; ++kk) {
+       fRMSSigmaFit[ii][jj][kk] = param.fRMSSigmaFit[ii][jj][kk];
+      }
+    }
+    for (Int_t jj = 0; jj < 4; ++jj) {
+      fParamS1[ii][jj] = param.fParamS1[ii][jj];
+      fErrorS1[ii][jj] = param.fErrorS1[ii][jj];
+    }
+    for (Int_t jj = 0; jj < 5; ++jj) {
+      fParamRMS1[ii][jj] = param.fParamRMS1[ii][jj];
+      fErrorRMS1[ii][jj] = param.fErrorRMS1[ii][jj];
+    }
+    fErrorRMSSys[ii] = param.fErrorRMSSys[ii];
+    for (Int_t jj = 0; jj < 2; ++jj){
+      fRMSSigmaRatio[ii][jj] = param.fRMSSigmaRatio[ii][jj];
+    }
+  }
+
+  SetWaveCorrectionMap( param.fWaveCorrectionMap );
+  SetResolutionYMap( param.fResolutionYMap );
 }
 
 
@@ -211,7 +262,6 @@ AliTPCClusterParam & AliTPCClusterParam::operator=(const AliTPCClusterParam& par
   // Assignment operator
   //
   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]){
@@ -232,6 +282,49 @@ AliTPCClusterParam & AliTPCClusterParam::operator=(const AliTPCClusterParam& par
       fPosZcor[1] = new TVectorD(*(param.fPosZcor[1]));
       fPosZcor[2] = new TVectorD(*(param.fPosZcor[2]));
     }
+    
+    for (Int_t ii = 0; ii < 2; ++ii) {
+      for (Int_t jj = 0; jj < 3; ++jj) {
+       for (Int_t kk = 0; kk < 4; ++kk) {
+         fParamS0[ii][jj][kk] = param.fParamS0[ii][jj][kk];
+         fErrorS0[ii][jj][kk] = param.fErrorS0[ii][jj][kk];
+         fParamRMS0[ii][jj][kk] = param.fParamRMS0[ii][jj][kk];
+         fErrorRMS0[ii][jj][kk] = param.fErrorRMS0[ii][jj][kk];
+       }
+       for (Int_t kk = 0; kk < 7; ++kk) {
+         fParamS0Par[ii][jj][kk] = param.fParamS0Par[ii][jj][kk];
+         fErrorS0Par[ii][jj][kk] = param.fErrorS0Par[ii][jj][kk];
+       }
+       for (Int_t kk = 0; kk < 6; ++kk) {
+         fParamSQ[ii][jj][kk] = param.fParamSQ[ii][jj][kk];
+         fErrorSQ[ii][jj][kk] = param.fErrorSQ[ii][jj][kk];
+         fParamRMSQ[ii][jj][kk] = param.fParamRMSQ[ii][jj][kk];
+         fErrorRMSQ[ii][jj][kk] = param.fErrorRMSQ[ii][jj][kk];
+       }
+       for (Int_t kk = 0; kk < 9; ++kk) {
+         fParamSQPar[ii][jj][kk] = param.fParamSQPar[ii][jj][kk];
+         fErrorSQPar[ii][jj][kk] = param.fErrorSQPar[ii][jj][kk];
+       }
+       for (Int_t kk = 0; kk < 2; ++kk) {
+         fRMSSigmaFit[ii][jj][kk] = param.fRMSSigmaFit[ii][jj][kk];
+       }
+      }
+      for (Int_t jj = 0; jj < 4; ++jj) {
+       fParamS1[ii][jj] = param.fParamS1[ii][jj];
+       fErrorS1[ii][jj] = param.fErrorS1[ii][jj];
+      }
+      for (Int_t jj = 0; jj < 5; ++jj) {
+       fParamRMS1[ii][jj] = param.fParamRMS1[ii][jj];
+       fErrorRMS1[ii][jj] = param.fErrorRMS1[ii][jj];
+      }
+      fErrorRMSSys[ii] = param.fErrorRMSSys[ii];
+      for (Int_t jj = 0; jj < 2; ++jj){
+       fRMSSigmaRatio[ii][jj] = param.fRMSSigmaRatio[ii][jj];
+      }
+    }
+    
+    SetWaveCorrectionMap( param.fWaveCorrectionMap );
+    SetResolutionYMap( param.fResolutionYMap );
   }
   return *this;
 }
@@ -264,6 +357,8 @@ AliTPCClusterParam::~AliTPCClusterParam(){
     delete fPosZcor[1];
     delete fPosZcor[2];
   }
+  delete fWaveCorrectionMap;
+  delete fResolutionYMap;
 }
 
 
@@ -1598,20 +1693,26 @@ Double_t  AliTPCClusterParam::GaussConvolution(Double_t x0, Double_t x1, Double_
   //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;
+  const Double_t kEpsilon = 0.0001;
+  const Double_t twoPi = TMath::TwoPi();
+  const Double_t hnorm = 0.5/TMath::Sqrt(twoPi);
+  const Double_t sqtwo = TMath::Sqrt(2.);
+
   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());
+    Double_t val = TMath::Gaus(x0,0,s0)*TMath::Gaus(x1,0,s1)/(s0*s1*twoPi);
     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 sigma = TMath::Sqrt(sigma2);
+  Double_t exp0 = TMath::Exp(-(k1*x0-k0*x1)*(k1*x0-k0*x1)/(2.*sigma2));
+  //
+  Double_t sigmaErf =  1./(2.*s0*s1*sqtwo*sigma);
+  Double_t k0s1s1 = 2.*k0*s1*s1;
+  Double_t k1s0s0 = 2.*k1*s0*s0;
+  Double_t erf0 = AliMathBase::ErfFast((sigma2-k0s1s1*x0-k1s0s0*x1)*sigmaErf);
+  Double_t erf1 = AliMathBase::ErfFast((sigma2+k0s1s1*x0+k1s0s0*x1)*sigmaErf);
+  Double_t norm = hnorm/sigma;
   Double_t val  = norm*exp0*(erf0+erf1);
   return val;
 }
@@ -1792,7 +1893,69 @@ Double_t  AliTPCClusterParam::QtotCorrection(Int_t sector, Int_t row, Float_t cp
 
 
 
+void AliTPCClusterParam::SetWaveCorrectionMap( THnBase *Map)
+{
+  //
+  // Set Correction Map for Y
+  //
+  delete fWaveCorrectionMap;
+  fWaveCorrectionMap = 0;
+  fWaveCorrectionMirroredPad = kFALSE;
+  fWaveCorrectionMirroredZ = kFALSE;
+  fWaveCorrectionMirroredAngle = kFALSE;
+  if( Map ){
+    fWaveCorrectionMap = dynamic_cast<THnBase*>( Map->Clone(Map->GetName()));
+    if( fWaveCorrectionMap ){
+      fWaveCorrectionMirroredPad   = ( fWaveCorrectionMap->GetAxis(3)->FindFixBin(0.5)<=1 );   // Pad axis is mirrored at 0.5
+      fWaveCorrectionMirroredZ     = ( fWaveCorrectionMap->GetAxis(1)->FindFixBin(0)<=1); // Z axis is mirrored at 0
+      fWaveCorrectionMirroredAngle = ( fWaveCorrectionMap->GetAxis(4)->FindFixBin(0.0)<=1 );   // Angle axis is mirrored at 0
+    }
+  }
+}
 
+void AliTPCClusterParam::SetResolutionYMap( THnBase *Map)
+{
+  //
+  // Set Resolution Map for Y
+  //
+  delete fResolutionYMap;
+  fResolutionYMap = 0;
+  if( Map ){
+    fResolutionYMap = dynamic_cast<THnBase*>( Map->Clone(Map->GetName()));
+  }
+}
 
+Float_t AliTPCClusterParam::GetWaveCorrection(Int_t Type, Float_t Z, Int_t QMax, Float_t Pad, Float_t angleY ) const
+{
+  //
+  // Correct Y cluster coordinate using a map
+  //
+
+  if( !fWaveCorrectionMap ) return 0;
+  Bool_t swapY = kFALSE;
+  Pad = Pad-(Int_t)Pad;
+
+  if( TMath::Abs(Pad-0.5)<1.e-8 ){// one pad clusters a stored in underflow bins
+    Pad = -1.; 
+  } else {
+    if( fWaveCorrectionMirroredPad && (Pad<0.5) ){ // cog axis is mirrored at 0.5
+       swapY = !swapY;
+       Pad = 1.0 - Pad;
+    }
+  }
 
+  if( fWaveCorrectionMirroredZ && (Z<0) ){ // Z axis is mirrored at 0
+    swapY = !swapY;
+    Z = -Z;
+  }
+  
+  if( fWaveCorrectionMirroredAngle && (angleY<0) ){ // Angle axis is mirrored at 0
+    angleY = -angleY;
+  }  
+  double var[5] = { Type, Z, QMax, Pad, angleY };
+  Long64_t bin = fWaveCorrectionMap->GetBin(var, kFALSE );
+  if( bin<0 ) return 0;
+  Double_t dY = fWaveCorrectionMap->GetBinContent(bin);
+  return (swapY ?-dY :dY);
+}