ATO-98 1.) distinguish between the refX radius and extrapolation radius 2.) Adding...
authormivanov <marian.ivanov@cern.ch>
Sat, 8 Nov 2014 16:46:26 +0000 (17:46 +0100)
committermivanov <marian.ivanov@cern.ch>
Sat, 22 Nov 2014 21:51:00 +0000 (22:51 +0100)
TPC/Calib/AliTPCCorrectionFit.cxx
TPC/Calib/AliTPCCorrectionFit.h
TPC/CalibMacros/RegisterCorrection.C

index d00e62e..8e81c39 100644 (file)
@@ -106,11 +106,14 @@ AliTPCCorrectionFit::~AliTPCCorrectionFit() {
   //
 }
 
-
-Double_t AliTPCCorrectionFit::EvalAt(Double_t phi, Double_t refX, Double_t theta, Int_t corr, Int_t ptype, Float_t wt, Float_t t1, Float_t t2){
+/*
+Double_t AliTPCCorrectionFit::EvalAt(Double_t phi, Double_t refX, Double_t evalX, Double_t theta, Int_t corr, Int_t ptype, Float_t wt, Float_t t1, Float_t t2){
   //
-  // Evalution at point using the lienar approximation
+  // Evaluation distortion at point using the linear approximation
   //
+  //  refX  -  reference X where phi and snp is calculated
+  //  evalX  - ref X where the distortion is evaluated
+
   if (wt<50){
     AliTPCCorrection* pcorr = AliTPCCorrection::GetVisualCorrection(corr);     
     pcorr->SetOmegaTauT1T2(wt,t1,t2);
@@ -119,15 +122,15 @@ Double_t AliTPCCorrectionFit::EvalAt(Double_t phi, Double_t refX, Double_t theta
   if (sector<0) sector+=18;
   Double_t y85=AliTPCCorrection::GetCorrSector(sector,85,theta,1,corr);
   Double_t y245=AliTPCCorrection::GetCorrSector(sector,245,theta,1,corr);
-  if (ptype==0) return y85+(y245-y85)*(refX-85.)/(245.-85.);
+  if (ptype==0) return y85+(y245-y85)*(evalX-85.)/(245.-85.);
   if (ptype==2) return (y245-y85)/(245.-85.);
   return 0;
 }
+*/
 
 
 
-
-Double_t AliTPCCorrectionFit::EvalAtPar(Double_t phi0, Double_t snp, Double_t refX, Double_t theta, Int_t corr, Int_t ptype, Int_t nsteps,Float_t wt, Float_t t1, Float_t t2){
+Double_t AliTPCCorrectionFit::EvalAtPar(Double_t phi0, Double_t snp, Double_t refX, Double_t evalX, Double_t theta, Int_t corr, Int_t ptype, Int_t nsteps,Float_t wt, Float_t t1, Float_t t2){
   //
   // Fit the distortion along the line with the parabolic model
   // We assume that the track are primaries  - where the vertex is at (0,0,0)
@@ -135,7 +138,8 @@ Double_t AliTPCCorrectionFit::EvalAtPar(Double_t phi0, Double_t snp, Double_t re
   // Parameters:
   //  phi0  -  phi at the entrance of the TPC
   //  snp   -  local inclination angle at the entrance of the TPC
-  //  refX  -  ref X where the distortion is evanluated
+  //  refX  -  reference X where phi and snp is calculated
+  //  evalX  - ref X where the distortion is evaluated
   //  theta -  inclination angle
   //  corr  -  internal number of the correction and distortion 
   //  ptype -  0 - local y distortion
@@ -145,23 +149,28 @@ Double_t AliTPCCorrectionFit::EvalAtPar(Double_t phi0, Double_t snp, Double_t re
   //           4 - distortion in the curvature
   //  nsteps - number of fit points
   //
-  // return value -  distortion at point refX with type ptype
+  // return value -  distortion at point evalX with type ptype
   //
   static TLinearFitter fitter(3,"pol2"); 
   fitter.ClearPoints();
   if (nsteps<3) nsteps=3;
   AliTPCCorrection* pcorr = AliTPCCorrection::GetVisualCorrection(corr);       
-  if pcorr==0) return 0;
+  if (pcorr==0) return 0;
   if (wt<50){
     pcorr->SetOmegaTauT1T2(wt,t1,t2);
   }
   Double_t deltaX=(245-85)/(nsteps);
+  Double_t curv=2.*snp/refX;
+  Double_t dPhi0=TMath::ASin(snp);
+
   for (Int_t istep=0; istep<(nsteps+1); istep++){
     //
     Double_t localX =85.+deltaX*istep;
-    Double_t localPhi=phi0+deltaX*snp*istep;
+    Double_t dPhi =TMath::ASin(0.5*localX*curv)-dPhi0;
+    Double_t localPhi=phi0+dPhi;
     Double_t sector = 9*localPhi/TMath::Pi();
     if (sector<0) sector+=18;
+    if (sector>18) sector-=18;
     Double_t dy=AliTPCCorrection::GetCorrSector(sector,localX,theta,1,corr);
     Double_t dlocalX=AliTPCCorrection::GetCorrSector(sector,localX,theta,0,corr);
     Double_t x[1]={localX-dlocalX};
@@ -172,30 +181,37 @@ Double_t AliTPCCorrectionFit::EvalAtPar(Double_t phi0, Double_t snp, Double_t re
   par[0]=fitter.GetParameter(0);
   par[1]=fitter.GetParameter(1);
   par[2]=fitter.GetParameter(2);
-
+  
   Double_t schi2= TMath::Sqrt(fitter.GetChisquare()/nsteps);
   if (schi2> fgMaxChi2HelixAt){
-    ::Error("AliTPCCorrectionFit::EvalAtHelix",TString::Format("%s\tbad chi2:\t%2.2f",corr->GetName(),schi2).Data());
+    ::Error("AliTPCCorrectionFit::EvalAtHelix",TString::Format("%s\tbad chi2:\t%2.2f",pcorr->GetName(),schi2).Data());
     return 0;
   }
-  if (ptype==0) return par[0]+par[1]*refX+par[2]*refX*refX;
-  if (ptype==2) return par[1]+2*par[2]*refX;
+  if (ptype==0) return par[0]+par[1]*evalX+par[2]*evalX*evalX;
+  if (ptype==2) return par[1]+2*par[2]*evalX;
   if (ptype==4) return par[2];
   if (ptype==5) return schi2;
   return 0;
 }
 
 
-Double_t AliTPCCorrectionFit::EvalAtHelix(Double_t phi0, Double_t snp, Double_t refX, Double_t theta, Int_t corr, Int_t ptype, Int_t nsteps, Float_t wt, Float_t t1, Float_t t2){
+Double_t AliTPCCorrectionFit::EvalAtHelix(Double_t phi0, Double_t snp, Double_t refX, Double_t evalX, Double_t theta, Int_t corr, Int_t ptype, Int_t nsteps, Float_t wt, Float_t t1, Float_t t2){
   //
   // Fit the distortion along the line with the helix model
   // FIXME - original trajectory to be changed - AliHelix to be used
   // We assume that the track are primaries  - where the vertex is at (0,0,0)
   //
+  // This procedure should emulate distortions as calibrated in the AliTPCcalibTime class
+  //     a.) AliTPCcalibTime::FillResHistoTPCTOF                     -  residuals, phi0 and snp  evaluated at reference refX  TOF clusters (refX=evalX)
+  //     b.) AliTPCcalibTime::FillResHistoTPCTRD                     -  residuals, phi0 and snp  evaluated at reference refX  TPCout radius (refX=evalX)
+  //     c.) AliTPCcalibTime::FillResHistoTPCITS                     -  residuals, phi0 and snp  evaluated at reference refX  TPCin radius (refX=evalX)
+  //     d.) AliTPCcalibTime::FillResHistoTPC  (vertex)              -  phi0 and snp  evaluated at reference refX  TPCin radius, residuals at vertex (refX!=evalX)
+  //         *         
   // Parameters:
   //  phi0  -  phi at the entrance of the TPC
   //  snp   -  local inclination angle at the entrance of the TPC
-  //  refX  -  ref X where the distortion is evanluated
+  //  refX  -  reference X where phi and snp is calculated
+  //  evalX  - ref X where the distortion is evaluated
   //  theta -  inclination angle
   //  corr  -  internal number of the correction and distortion 
   //  ptype -  0 - local y distortion
@@ -209,7 +225,7 @@ Double_t AliTPCCorrectionFit::EvalAtHelix(Double_t phi0, Double_t snp, Double_t
   //
   static TLinearFitter fitter(3,"pol2"); 
   fitter.ClearPoints();
-
+  //
   if (nsteps<4) nsteps=4;
   Double_t deltaX=(245-85)/(nsteps);
   AliRieman rieman(nsteps);
@@ -217,12 +233,16 @@ Double_t AliTPCCorrectionFit::EvalAtHelix(Double_t phi0, Double_t snp, Double_t
     AliTPCCorrection* pcorr = AliTPCCorrection::GetVisualCorrection(corr);     
     pcorr->SetOmegaTauT1T2(wt,t1,t2);
   }
+  Double_t curv=2.*snp/refX;
+  Double_t dPhi0=TMath::ASin(snp);
   for (Int_t istep=0; istep<(nsteps+1); istep++){
     //
     Double_t localX =85.+deltaX*istep;
-    Double_t localPhi=phi0+deltaX*snp*istep;
+    Double_t dPhi =TMath::ASin(0.5*localX*curv)-dPhi0;
+    Double_t localPhi=phi0+dPhi;
     Double_t sector = 9*localPhi/TMath::Pi();
     if (sector<0) sector+=18;
+    if (sector>18) sector-=18;
     Double_t dy=AliTPCCorrection::GetCorrSector(sector,localX,theta,1,corr);
     Double_t dlocalX=AliTPCCorrection::GetCorrSector(sector,localX,theta,0,corr);
     Double_t x[1]={localX-dlocalX};
@@ -238,8 +258,8 @@ Double_t AliTPCCorrectionFit::EvalAtHelix(Double_t phi0, Double_t snp, Double_t
     ::Error("AliTPCCorrectionFit::EvalAtHelix",TString::Format("Bad chi2\t%2.2f",schi2).Data());
     return 0;
   }
-  if (ptype==0) return rieman.GetYat(refX);
-  if (ptype==2) return rieman.GetDYat(refX);
+  if (ptype==0) return rieman.GetYat(evalX);
+  if (ptype==2) return rieman.GetDYat(evalX);
   if (ptype==4) return rieman.GetC();
   return 0;
 }
index b3ca34e..b52b135 100644 (file)
@@ -20,9 +20,9 @@ public:
   //
   //
   //
-  static Double_t EvalAt(Double_t phi, Double_t refX, Double_t theta, Int_t corr, Int_t ptype, Float_t wt=100, Float_t t1=1, Float_t t2=1);
-  static Double_t EvalAtPar(Double_t phi, Double_t snp, Double_t refX, Double_t theta, Int_t corr, Int_t ptype, Int_t nstep, Float_t wt=100,  Float_t t1=1, Float_t t2=1);
-  static Double_t EvalAtHelix(Double_t phi, Double_t snp, Double_t refX, Double_t theta, Int_t corr, Int_t ptype, Int_t nstep,  Float_t wt=100, Float_t t1=1, Float_t t2=1);
+  //static Double_t EvalAt(Double_t phi, Double_t refX, Double_t theta,Double_t evalX,  Int_t corr, Int_t ptype, Float_t wt=100, Float_t t1=1, Float_t t2=1);
+  static Double_t EvalAtPar(Double_t phi, Double_t snp, Double_t refX, Double_t evalX, Double_t theta, Int_t corr, Int_t ptype, Int_t nstep, Float_t wt=100,  Float_t t1=1, Float_t t2=1);
+  static Double_t EvalAtHelix(Double_t phi, Double_t snp, Double_t refX, Double_t evalX, Double_t theta, Int_t corr, Int_t ptype, Int_t nstep,  Float_t wt=100, Float_t t1=1, Float_t t2=1);
   //
   // Make distortion maps
   //
index 0287bef..8590c37 100644 (file)
@@ -78,6 +78,8 @@
 #include  "TGeoGlobalMagField.h"
 #include  "TROOT.h"
 #include "TLinearFitter.h"
+#include "TStopwatch.h"
+#include "AliTPCCorrectionFit.h"
 #endif
 
 
@@ -1436,15 +1438,40 @@ AliTPCComposedCorrection * GetCorrectionFromFile(){
 void TestParExample(){
   //
   // dz shift example: AliTPCCorrection::AddVisualCorrection(rocDzIROCA,705); 
-  // => parabolic fit not appropriate - helix fit differnt results
+  // => parabolic fit and helix fit agrees - once significant ammount of points used
+  //    160 point - agreement  ~2%; 
+  //     80 points - agreement ~5%
   AliTPCCorrection* corr = AliTPCCorrection::GetVisualCorrection(705);  
   corr->SetOmegaTauT1T2(0.33,1,1);
-  TF1 f705Par("f705","AliTPCCorrectionFit::EvalAtPar(0,0,x,0.1,705,0,20)",0,360);
+  TF1 f705Par("f705","AliTPCCorrectionFit::EvalAtPar(0,0,85, x,0.1,705,0,80)",0,360);
   f705Par.SetLineColor(2);f705Par.SetNpx(500);
-  TF1 f705Helix("f705Helix","AliTPCCorrectionFit::EvalAtHelix(0,0,x,0.1,705,0,20)",0,360);
+  TF1 f705Helix("f705Helix","AliTPCCorrectionFit::EvalAtHelix(0,0,85,x,0.1,705,0,80)",0,360);
   f705Helix.SetLineColor(4);f705Helix.SetNpx(500);
   f705Helix.Draw();
   f705Par.Draw("same");
+
 }
 
 
+
+void TestFitSpeed(Int_t nEvals){
+  //
+  //  test speed of helix fit/ resp. parabolic fir
+  //
+  TStopwatch timerh; 
+  ::Info("TestFitSpeed","Helix fit");
+  for (Int_t i=0;i<nEvals; i++) AliTPCCorrectionFit::EvalAtPar(0,0,85,gRandom->Rndm(),0.1,705,0,80); 
+  timerh.Print();
+  TStopwatch timerP; 
+  ::Info("TestFitSpeed","Parabolicfit");
+  for (Int_t i=0;i<nEvals; i++) AliTPCCorrectionFit::EvalAtPar(0,0,85,gRandom->Rndm(),0.1,705,0,80); 
+  timerP.Print();
+  /*
+    For the test system CPU comsumption is the same:
+    I-TestFitSpeed: Helix fit
+    Real time 0:00:03, CP time 3.310
+    I-TestFitSpeed: Parabolicfit
+    Real time 0:00:03, CP time 3.280
+  */
+}
+