]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TPC/AliTPCPointCorrection.cxx
Avoiding crash in crash in GetELossRandomKFastR
[u/mrichter/AliRoot.git] / TPC / AliTPCPointCorrection.cxx
index 1b0fb3ea17a4d228a0ce6fda96df413b148be25f..28afd3c17fd0302e38231e31769a753e4bc25737 100644 (file)
@@ -28,7 +28,7 @@
     
 */
 
-
+#include "AliTPCcalibDB.h"
 #include "TLinearFitter.h"
 #include "Riostream.h"
 #include "TList.h"
@@ -39,6 +39,7 @@
 #include "TVectorD.h"
 #include "AliLog.h"
 #include "AliTPCROC.h"
+#include "AliTPCClusterParam.h"
 #include "AliTPCPointCorrection.h"
 
 AliTPCPointCorrection* AliTPCPointCorrection::fgInstance = 0;
@@ -52,11 +53,27 @@ AliTPCPointCorrection::AliTPCPointCorrection():
   fParamOutRVersion(0),
   fErrorsOutR(38),
   fErrorsOutZ(38),
-  fParamOutZVersion(0)
+  fParamOutZVersion(0),
+  //
+  fXIO(0),
+  fXmiddle(0),
+  fXquadrant(0),
+  fArraySectorIntParam(36), // array of sector alignment parameters
+  fArraySectorIntCovar(36), // array of sector alignment covariances 
+  //
+  // Kalman filter for global alignment
+  //
+  fSectorParam(0),     // Kalman parameter   
+  fSectorCovar(0)     // Kalman covariance  
+
 {
   //
   // Default constructor
   //
+  AliTPCROC * roc = AliTPCROC::Instance();
+  fXquadrant = roc->GetPadRowRadii(36,53);
+  fXmiddle   = ( roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
+  fXIO       = ( roc->GetPadRowRadii(0,roc->GetNRows(0)-1)+roc->GetPadRowRadii(36,0))*0.5;
 }
 
 AliTPCPointCorrection::AliTPCPointCorrection(const Text_t *name, const Text_t *title):
@@ -66,11 +83,30 @@ AliTPCPointCorrection::AliTPCPointCorrection(const Text_t *name, const Text_t *t
   fParamOutRVersion(0),
   fErrorsOutR(38),
   fErrorsOutZ(38),
-  fParamOutZVersion(0)
+  fParamOutZVersion(0),
+  //
+  // 
+  //
+  fXIO(0),
+  fXmiddle(0),
+  fXquadrant(0),
+  fArraySectorIntParam(36), // array of sector alignment parameters
+  fArraySectorIntCovar(36), // array of sector alignment covariances 
+  //
+  // Kalman filter for global alignment
+  //
+  fSectorParam(0),     // Kalman parameter   for A side
+  fSectorCovar(0)     // Kalman covariance  for A side 
+  
 {
   //
   // Non default constructor
   //
+  AliTPCROC * roc = AliTPCROC::Instance();
+  fXquadrant = roc->GetPadRowRadii(36,53);
+  fXmiddle   = ( roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
+  fXIO       = ( roc->GetPadRowRadii(0,roc->GetNRows(0)-1)+roc->GetPadRowRadii(36,0))*0.5;
+
 }
 
 AliTPCPointCorrection::~AliTPCPointCorrection(){
@@ -151,7 +187,8 @@ Double_t AliTPCPointCorrection::CorrectionOutR0(Bool_t isGlobal, Bool_t type,  D
   Double_t radius = TMath::Sqrt(cx*cx+cy*cy);  
   AliTPCROC * roc = AliTPCROC::Instance();
   Double_t xout = roc->GetPadRowRadiiUp(roc->GetNRows(37)-1);
-  Double_t dout = xout-radius;
+  Double_t dout = xout-radius;  
+  if (dout<0) return 0;
   //drift
   Double_t dr   = 0.5 - TMath::Abs(cz/250.);
   //
@@ -171,9 +208,69 @@ Double_t AliTPCPointCorrection::CorrectionOutR0(Bool_t isGlobal, Bool_t type,  D
   result+=(*vec)[5]*eoutp*dr*dr;
   result+=(*vec)[6]*eoutm*dr*dr;
   return result;
+}
+
+Double_t AliTPCPointCorrection::RPhiCOGCorrection(Int_t sector, Int_t padrow, Float_t pad, Float_t cy, Float_t y, Float_t z, Float_t ky,Float_t qMax, Float_t threshold){
+  //
+  // Calculates COG corection in RPHI direction
+  // cluster and track position  y is supposed to be corrected before for other effects
+  // (e.g ExB and alignemnt)
+  // Rphi distortion dependeds on the distance to the edge-pad, distance to the wire edge and
+  // relative distance to the center of the pad. Therefore the y position is trnsfromed to the 
+  // pad coordiante frame (correction offset (ExB alignemnt) substracted). 
+  //   
+  // Input parameters:
+  //
+  // sector - sector number - 0-71  - cl.GetDetector()
+  // padrow - padrow number -0-63 -IROC 0-95 OROC - cl->GetRow()
+  // pad    - mean pad number  - cl->GetPad()
+  // cy     - cluster y        - cl->GetY()
+  //  y     - track -or cluster y
+  //  qMax  - cluster max charge - cl-.GetMax()
+  //  threshold - clusterer threshold
+  //
+  AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
+  if (!clparam) {
+    AliFatal("TPC OCDB not initialized"); 
+    return 0;
+  }
+  Int_t padtype=0;
+  if (sector>=36) padtype = (padrow>64)?2:1;
+  Double_t padwidth=(padtype==0)? 0.4:0.6;
 
+  Double_t sigma = clparam->GetRMS0(0,padtype,250-TMath::Abs(z),ky)/padwidth;
+  //
+  Int_t   npads   =  AliTPCROC::Instance()->GetNPads(sector,padrow);
+  Float_t yshift  =  TMath::Abs(cy)-TMath::Abs((-npads*0.5+pad)*padwidth);   // the clusters can be corrected before
+                                            // shift to undo correction
+  
+  y -= yshift*((y>0)?1.:-1.);                             // y in the sector coordinate
+  Double_t y0     = npads*0.5*padwidth;
+  Double_t dy     = (TMath::Abs(y0)-TMath::Abs(y))/padwidth-0.5;  // distance to  the center of pad0   
+  //
+  Double_t padcenter = TMath::Nint(dy);
+  Double_t sumw=0;
+  Double_t sumwi=0;
+  for (Double_t ip=padcenter-2.5; ip<=padcenter+2.5;ip++){
+    Double_t weightGaus = qMax*TMath::Exp(-(dy-ip)*(dy-ip)/(2*sigma*sigma));
+    Double_t ypad       = (ip-npads*0.5)*padwidth;
+    Double_t weightGain = GetEdgeQ0(sector,padrow,ypad);
+    //
+    Double_t weight = TMath::Nint(weightGaus*weightGain);
+    if (ip<0 &&weight> threshold) weight = threshold;  // this is done in cl finder
+    if (weight<0) continue;
+    sumw+=weight;
+    sumwi+=weight*(ip);    
+  }
+  Double_t result =0;
+  if (sumw>0) result = sumwi/sumw;
+  result = (result-dy)*padwidth;
+  return result;
 }
 
+
+
+
 Double_t AliTPCPointCorrection::CorrectionOutZ0(Bool_t isGlobal, Bool_t type,  Double_t cx, Double_t cy, Double_t cz, Int_t sector){
   //
   // return dR correction - for correction version 0 
@@ -195,6 +292,7 @@ Double_t AliTPCPointCorrection::CorrectionOutZ0(Bool_t isGlobal, Bool_t type,  D
   AliTPCROC * roc = AliTPCROC::Instance();
   Double_t xout = roc->GetPadRowRadiiUp(roc->GetNRows(37)-1);
   Double_t dout = xout-radius;
+  if (dout<0) return 0;
   //drift
   Double_t dr   = 0.5 - TMath::Abs(cz/250.);
   //
@@ -217,5 +315,162 @@ Double_t AliTPCPointCorrection::CorrectionOutZ0(Bool_t isGlobal, Bool_t type,  D
 
 }
 
+Double_t  AliTPCPointCorrection::GetEdgeQ0(Int_t sector, Int_t padrow, Float_t y){
+  //
+  /* TF1 fexp("fexp","1-exp(-[0]*(x-[1]))",0,20)
+          | param [0] | param [1]
+     IROC | 4.71413   | 1.39558
+     OROC1| 2.11437   | 1.52643
+     OROC2| 2.15082   | 1.53537 
+  */
+
+  Double_t params[2]={100,0};
+  if (sector<36){
+    params[0]=4.71413; params[1]=1.39558;
+  }
+  if (sector>=36){
+    if (padrow<64) {  params[0]=2.114; params[1]=1.526;}
+    if (padrow>=64){  params[0]=2.15; params[1]=1.535;}
+  }
+  Double_t result= 1;
+  Double_t xrow  = AliTPCROC::Instance()->GetPadRowRadii(sector,padrow);
+  Double_t ymax  = TMath::Tan(TMath::Pi()/18.)*xrow;
+  Double_t dedge = ymax-TMath::Abs(y);
+  if (dedge-params[1]<0)             return 0;
+  if (dedge>10.*params[1]) return 1;
+  result = 1.-TMath::Exp(-params[0]*(dedge-params[1]));
+  return result;
+}
+
+Double_t AliTPCPointCorrection::SRPhiCOGCorrection(Int_t sector, Int_t padrow, Float_t pad, Float_t cy, Float_t y, Float_t z, Float_t ky,Float_t qMax, Float_t threshold){
+  return fgInstance->RPhiCOGCorrection(sector, padrow, pad, cy, y,z, ky, qMax, threshold);
+}
+
+Double_t AliTPCPointCorrection::SGetEdgeQ0(Int_t sector, Int_t padrow, Float_t y){
+  //
+  //
+  return fgInstance->GetEdgeQ0(sector, padrow, y);
+}
+
+void     AliTPCPointCorrection::AddCorrectionSector(TObjArray & sideAPar, TObjArray &sideCPar, TObjArray & sideACov, TObjArray &sideCCov, Bool_t reset){
+  //
+  //
+  //
+  for (Int_t isec=0;isec<36;isec++){
+    TMatrixD * mat1 = (TMatrixD*)(sideAPar.At(isec));
+    TMatrixD * mat1Covar = (TMatrixD*)(sideACov.At(isec));
+    if (!mat1) continue;
+    TMatrixD * mat0 = (TMatrixD*)(fArraySectorIntParam.At(isec));
+    TMatrixD * mat0Covar = (TMatrixD*)(fArraySectorIntCovar.At(isec));
+    if (!mat0) {
+      fArraySectorIntParam.AddAt(mat1->Clone(),isec); 
+      fArraySectorIntCovar.AddAt(mat1Covar->Clone(),isec); 
+      continue;
+    }
+    (*mat0Covar)=(*mat1Covar);      
+    if (reset)   (*mat0)=(*mat1);
+    if (!reset) (*mat0)+=(*mat1);
+  }
+
+  for (Int_t isec=0;isec<36;isec++){
+    TMatrixD * mat1 = (TMatrixD*)(sideCPar.At(isec));
+    TMatrixD * mat1Covar = (TMatrixD*)(sideCCov.At(isec));
+    if (!mat1) continue;
+    TMatrixD * mat0 = (TMatrixD*)(fArraySectorIntParam.At(isec));
+    TMatrixD * mat0Covar = (TMatrixD*)(fArraySectorIntCovar.At(isec));
+    if (!mat0) {
+      fArraySectorIntParam.AddAt(mat1->Clone(),isec); 
+      fArraySectorIntCovar.AddAt(mat1Covar->Clone(),isec); 
+      continue;
+    }
+    (*mat0Covar)=(*mat1Covar);      
+    if (reset)   (*mat0)=(*mat1);
+    if (!reset) (*mat0)+=(*mat1);
+  }
+}
+
+
+Double_t AliTPCPointCorrection::GetCorrectionSector(Int_t coord, Int_t sector, Double_t lx, Double_t ly, Double_t /*lz*/, Int_t quadrant){
+  //
+  // Get position correction for given sector
+  //
+
+  TMatrixD * param = (TMatrixD*)fArraySectorIntParam.At(sector%36);
+  if (!param) return 0;
+  if (quadrant<0){   //recalc quadrant if not specified
+    if (lx<fXIO) quadrant=0;
+    if(lx>fXIO){
+      if (lx<fXquadrant) {
+       if (ly<0) quadrant=1;
+       if (ly>0) quadrant=2;
+      }
+      if (lx>=fXquadrant) {
+       if (ly<0) quadrant=3;
+       if (ly>0) quadrant=4;
+      }
+    }    
+  }
+  Double_t a10 = (*param)(quadrant*6+0,0);
+  Double_t a20 = (*param)(quadrant*6+1,0);
+  Double_t a21 = (*param)(quadrant*6+2,0);
+  Double_t dx  = (*param)(quadrant*6+3,0);
+  Double_t dy  = (*param)(quadrant*6+4,0);
+  Double_t dz  = (*param)(quadrant*6+5,0);
+  Double_t deltaX = lx-fXmiddle;
+  if (coord==0) return dx;
+  if (coord==1) return dy+deltaX*a10;
+  if (coord==2) return dz+deltaX*a20+ly*a21;
+  if (coord==3) return a10;
+  if (coord==4) return a20;
+  if (coord==5) return a21;
+  //
+  return 0;
+}
+
+Double_t AliTPCPointCorrection::SGetCorrectionSector(Int_t coord, Int_t sector, Double_t lx, Double_t ly, Double_t lz, Int_t quadrant){
+  //
+  //
+  //
+  if (!Instance()) return 0;
+  return Instance()->GetCorrectionSector(coord,sector,lx,ly,lz, quadrant);
+}
+
+
+
+Double_t AliTPCPointCorrection::GetCorrection(Int_t coord, Int_t sector, Double_t lx, Double_t ly, Double_t /*lz*/){
+  //
+  // Get position correction for given sector
+  //
+  if (!fSectorParam) return 0;
+  
+  Double_t a10 = (*fSectorParam)(sector*6+0,0);
+  Double_t a20 = (*fSectorParam)(sector*6+1,0);
+  Double_t a21 = (*fSectorParam)(sector*6+2,0);
+  Double_t dx  = (*fSectorParam)(sector*6+3,0);
+  Double_t dy  = (*fSectorParam)(sector*6+4,0);
+  Double_t dz  = (*fSectorParam)(sector*6+5,0);
+  Double_t deltaX = lx-fXIO;
+  //
+  if (coord==0) return dx;
+  if (coord==1) return dy+deltaX*a10;
+  if (coord==2) return dz+deltaX*a20+ly*a21;
+  if (coord==3) return a10;
+  if (coord==4) return a20;
+  if (coord==5) return a21;
+  return 0;
+}
+
+Double_t AliTPCPointCorrection::SGetCorrection(Int_t coord, Int_t sector, Double_t lx, Double_t ly, Double_t lz){
+  //
+  //
+  //
+  if (!Instance()) return 0;
+  return Instance()->GetCorrection(coord,sector,lx,ly,lz);
+}
+
+
+
+
+