o first try to create a LUT from the residual distortions

diff --git a/TPC/Base/AliTPCCorrectionLookupTable.cxx b/TPC/Base/AliTPCCorrectionLookupTable.cxx
index 3ab8da44359d2da42a1a68d59788851adbcf743c..177f4fea2c95d47c8f460eb525cd043ccd9047cf 100644 (file)
--- a/TPC/Base/AliTPCCorrectionLookupTable.cxx
+++ b/TPC/Base/AliTPCCorrectionLookupTable.cxx
@@ -264,6 +264,66 @@ void AliTPCCorrectionLookupTable::InitTables()
   }
 }
 
+//_________________________________________________________________________________________
+void AliTPCCorrectionLookupTable::CreateLookupTableFromResidualDistortion(THn &resDist)
+{
+  //
+  // create lookup table from residual distortions stored in a 3d histogram
+  // assume dimensions are r, phi, z
+  //
+  if (fNR==0) {
+    AliError("Limits are not set yet. Please use one of the Set..Limits functions first");
+    return;
+  }
+  
+  TStopwatch s;
+  
+  ResetTables();
+  InitTables();
+
+  Double_t x[3]={0.,0.,0.};
+  
+  for (Int_t iPhi=0; iPhi<fNPhi; ++iPhi){
+    const Double_t phi=fLimitsPhi(iPhi);
+    x[1]=phi;
+    //
+    TMatrixF &mDxDist   = *fLookUpDxDist[iPhi];
+    TMatrixF &mDyDist   = *fLookUpDyDist[iPhi];
+    TMatrixF &mDzDist   = *fLookUpDzDist[iPhi];
+    //
+    TMatrixF &mDxCorr   = *fLookUpDxCorr[iPhi];
+    TMatrixF &mDyCorr   = *fLookUpDyCorr[iPhi];
+    TMatrixF &mDzCorr   = *fLookUpDzCorr[iPhi];
+    
+    for (Int_t ir=0; ir<fNR; ++ir){
+      const Double_t r=fLimitsR(ir);
+      x[0]=r;
+      
+      for (Int_t iz=0; iz<fNZ; ++iz){
+        const Double_t z=fLimitsZ(iz);
+        x[2]=z;
+
+        const Double_t drphi = resDist.GetBinContent(resDist.GetBin(x));
+        Double_t dx[3]={0.,drphi,0.};
+        
+        // transform rphi distortions (local y, so dy') to a global distortion
+        // assume no radial distortion (dx' = 0)
+        // assume no residual distortion in z for the moment
+        Double_t cs=TMath::Cos(phi), sn=TMath::Sin(phi), lx=dx[0];
+        dx[0]=lx*cs - dx[1]*sn; dx[1]=lx*sn + dx[1]*cs;
+
+        mDxDist(ir,iz)=dx[0];
+        mDyDist(ir,iz)=dx[1];
+        mDzDist(ir,iz)=dx[2];
+
+        mDxCorr(ir,iz)=-dx[0];
+        mDyCorr(ir,iz)=-dx[1];
+        mDzCorr(ir,iz)=-dx[2];
+      }
+    }
+  }
+}
+  
 //_________________________________________________________________________________________
 void AliTPCCorrectionLookupTable::InitTablesPhiBin(Int_t iPhi)
 {

diff --git a/TPC/Base/AliTPCCorrectionLookupTable.h b/TPC/Base/AliTPCCorrectionLookupTable.h
index 9a9e6a8c55ed3193c8c0b5f7d1515c22b4282676..67f7e89f2fd455480a8bb31d7ef2fd616a5f7c7f 100644 (file)
--- a/TPC/Base/AliTPCCorrectionLookupTable.h
+++ b/TPC/Base/AliTPCCorrectionLookupTable.h
@@ -9,6 +9,7 @@
 #include "AliTPCCorrection.h"
 #include <TVectorD.h>
 #include <TMatrixFfwd.h>
+#include <THn.h>
 
 class AliTPCCorrectionLookupTable : public AliTPCCorrection {
 
@@ -23,6 +24,8 @@ public:
   void CreateLookupTable(AliTPCCorrection &tpcCorr, Float_t stepSize=5.);
   void CreateLookupTableSinglePhi(AliTPCCorrection &tpcCorr, Int_t iPhi, Float_t stepSize=5.);
 
+  void CreateLookupTableFromResidualDistortion(THn &resDist);
+  
   void MergePhiTables(const char* files);
 
   void   SetFillCorrection(Bool_t fill) { fFillCorrection=fill;   }

diff --git a/TPC/Upgrade/macros/AnaDelta.C b/TPC/Upgrade/macros/AnaDelta.C
index 791c87d7e926b32bf882932a993be4c4ed8ae51e..32d38bdf91a06f906eac9af85736e4292c89198a 100644 (file)
--- a/TPC/Upgrade/macros/AnaDelta.C
+++ b/TPC/Upgrade/macros/AnaDelta.C
@@ -96,7 +96,11 @@ void AnaDeltaTree(TString file, TString outFile="deltas_tree.root")
   hn->GetAxis(0)->Set(vR  ->GetNrows()-1, vR  ->GetMatrixArray());
   hn->GetAxis(1)->Set(vPhi->GetNrows()-1, vPhi->GetMatrixArray());
   hn->GetAxis(2)->Set(vZ  ->GetNrows()-1, vZ  ->GetMatrixArray());
-  
+
+  hn->GetAxis(0)->SetNameTitle("r","r (cm)");
+  hn->GetAxis(1)->SetNameTitle("phi","#varphi");
+  hn->GetAxis(2)->SetNameTitle("z","z (cm)");
+  hn->GetAxis(3)->SetNameTitle("drphi","#Delta(r#varphi)");
   
   for (Int_t iev=0; iev<t->GetEntries(); ++iev) {
     t->GetEntry(iev);
@@ -302,6 +306,17 @@ void DumpHn(THn *hn, TTreeSRedirector &stream)
   TAxis *ar   = hn->GetAxis(0);
   TAxis *aphi = hn->GetAxis(1);
   TAxis *az   = hn->GetAxis(2);
+
+  // output Hn
+  const Int_t nbins=3;
+  Int_t bins[nbins]    = {1,1,1};
+  Double_t xmin[nbins] = {0.,0.,0.};
+  Double_t xmax[nbins] = {1.,1.,1.};
+  THnF hnRes("hnRes", "hnRes", nbins, bins, xmin, xmax);
+
+  ar  ->Copy(*hnRes.GetAxis(0));
+  aphi->Copy(*hnRes.GetAxis(1));
+  az  ->Copy(*hnRes.GetAxis(2));
   
   
   for (Int_t iz=0; iz<az->GetNbins(); ++iz) {
@@ -325,40 +340,54 @@ void DumpHn(THn *hn, TTreeSRedirector &stream)
         Double_t cr   = ar->GetBinCenter(ir+1);
         Double_t cphi = aphi->GetBinCenter(iphi+1);
         Double_t cz   = az->GetBinCenter(iz+1);
-        hProj->SetNameTitle(Form("h_%02d_%02d_%02d",iz+1, iphi+1, ir+1),
-                            Form("z,phi,r: %02d,%02d,%02d (%.2f, %.2f, %.2f)",iz+1,iphi+1,ir+1, cz, cphi, cr )
+        hProj->SetNameTitle(Form("h_%02d_%02d_%02d",iz, iphi, ir),
+                            Form("z,phi,r: %02d,%02d,%02d (%.2f, %.2f, %.2f)",iz,iphi,ir, cz, cphi, cr )
         );
-        hProj->Fit(&fg,"QR");
+        hProj->Fit(&fg,"LMQR");
         arrFits.Add(hProj);
 
-        Double_t mean     = fg.GetParameter(1);
-        Double_t meanErr  = fg.GetParError(1);
-        Double_t sigma    = fg.GetParameter(2);
-        Double_t sigmaErr = fg.GetParError(2);
+        Float_t mean     = fg.GetParameter(1);
+        Float_t meanErr  = fg.GetParError(1);
+        Float_t sigma    = fg.GetParameter(2);
+        Float_t sigmaErr = fg.GetParError(2);
         Int_t    entries  = hProj->GetEntries();
-        Double_t chi2ndf  = fg.GetChisquare()/fg.GetNDF();
-
+        Float_t chi2ndf  = fg.GetChisquare()/fg.GetNDF();
+        Float_t mean2    = hProj->GetMean();
+        Float_t meanErr2 = hProj->GetMeanError();
+        Float_t rms2     = hProj->GetRMS();
+        Float_t rmsErr2  = hProj->GetRMSError();
+        
         stream << "d" <<
-        "ir="        << ir       <<
-        "iphi="      << iphi     <<
-        "iz="        << iz       <<
-        "cr="        << cr       <<
-        "cphi="      << cphi     <<
-        "cz="        << cz       <<
-        "mean="      << mean     <<
-        "meanErr="   << meanErr  <<
-        "sigma="     << sigma    <<
-        "sigmaErr="  << sigmaErr <<
-        "entries="   << entries  <<
-        "chi2ndf="   << chi2ndf  <<
+        "ir="          << ir       <<
+        "iphi="        << iphi     <<
+        "iz="          << iz       <<
+        "cr="          << cr       <<
+        "cphi="        << cphi     <<
+        "cz="          << cz       <<
+        "mean="        << mean     <<
+        "meanErr="     << meanErr  <<
+        "sigma="       << sigma    <<
+        "sigmaErr="    << sigmaErr <<
+        "histMean="    << mean2    <<
+        "histMeanErr=" << meanErr2 <<
+        "histRMS="     << rms2    <<
+        "histRMSErr="  << rmsErr2 <<
+        "entries="     << entries  <<
+        "chi2ndf="     << chi2ndf  <<
         "\n";
+
+        Double_t x[nbins]={cr, cphi, cz};
+        if (meanErr<0.3) hnRes.Fill(x,mean);
       }
     }
     stream.GetFile()->cd();
     arrFits.Write(0x0,TObject::kSingleKey);
     gROOT->cd();
   }
-  
+
+  stream.GetFile()->cd();
+  hnRes.Write();
+  gROOT->cd();
 }
 
 //______________________________________________________________________________