From 1a26e63a000f67a44b43881a9238255dc167bd3a Mon Sep 17 00:00:00 2001
From: shahoian <ruben.shahoyan@cern.ch>
Date: Mon, 11 Jul 2016 20:39:56 +0200
Subject: [PATCH] Method for inversion consistency check (write tree with
 closure test)

---
 TPC/TPCcalib/AliTPCDcalibRes.cxx | 68 +++++++++++++++++++++++++++++++-
 TPC/TPCcalib/AliTPCDcalibRes.h   | 13 +++++-
 TPC/TPCcalib/TPCcalibLinkDef.h   |  1 +
 3 files changed, 79 insertions(+), 3 deletions(-)

diff --git a/TPC/TPCcalib/AliTPCDcalibRes.cxx b/TPC/TPCcalib/AliTPCDcalibRes.cxx
index 2b3c89cb62c..966e1b4aeeb 100644
--- a/TPC/TPCcalib/AliTPCDcalibRes.cxx
+++ b/TPC/TPCcalib/AliTPCDcalibRes.cxx
@@ -564,8 +564,8 @@ void AliTPCDcalibRes::ReProcessFromResVoxTree(const char* resTreeFile, Bool_t ba
   if (!LoadDataFromResTree(resTreeFile)) return;
   ReProcessResiduals();
   //
-  delete fChebCorr; fChebCorr = 0;
-  delete fChebDist; fChebDist = 0;
+  //  delete fChebCorr; fChebCorr = 0;
+  //  delete fChebDist; fChebDist = 0;
   //
   CreateCorrectionObject();
   //
@@ -3823,6 +3823,7 @@ void AliTPCDcalibRes::CreateCorrectionObject()
     return;
   }
   //
+  delete fChebCorr; fChebCorr = 0;
   TString name = Form("run%d_%lld_%lld",fRun,fTMin,fTMax);
   fChebCorr = new AliTPCChebCorr(name.Data(),name.Data(),
 				 fChebPhiSlicePerSector,fChebZSlicePerSide,1.0f);
@@ -4350,6 +4351,7 @@ void AliTPCDcalibRes::CreateDistortionObject()
     return;
   }
   TString name = Form("run%d_%lld_%lld_InvDist",fRun,fTMin,fTMax);
+  delete fChebDist; fChebDist = 0;
   fChebDist = new AliTPCChebDist(name.Data(),name.Data(),fChebPhiSlicePerSector,fChebZSlicePerSide,1.0f);
   fChebDist->SetUseFloatPrec(kFALSE);
   fChebDist->SetRun(fRun);
@@ -4408,3 +4410,65 @@ void AliTPCDcalibRes::BringToActiveBoundary(int sect36, float xyz[3]) const
   if (xyz[kResZ]*side<0)    xyz[kResZ] = 0;
   //
 }
+
+//____________________________________________
+void AliTPCDcalibRes::WriteDistCorTestTree(int nx, int nphi2sect,int nzside)
+{
+  // write test tree with distortions and corrections
+  const float kMaxY2X = TMath::Tan(0.5f*kSecDPhi);
+  //
+  dcTst_t dc, *dcp=&dc;
+  //
+  if (!fChebDist) {
+    AliError("Cheb distortion is not ready yet");
+    return;
+  }
+  fChebDist->Init();
+  //
+  TFile* flOut = new TFile("DistCorrTest.root","recreate");
+  TTree* resTree = new TTree("dc","DistCorr test");
+  resTree->Branch("dc", &dcp);
+  //
+  float dz2x = 1./nzside;
+  float dy2x = 2.*kMaxY2X/nphi2sect;
+  float dx = (kMaxX-kMinX)/nx;
+  float dist[4],corr[4];
+  for (int is=0;is<kNSect2;is++) { 
+    dc.bsec = is;
+    for (int ix=0;ix<nx;ix++) {
+      float x = kMinX + dx*(ix+0.5);
+      dc.bvox[kVoxX] = ix;
+      dc.xyz[kResX] = x;
+      //
+      for (int ip=0;ip<nphi2sect;ip++) {
+	float y2x = -kMaxY2X + dy2x*(0.5+ip);
+	dc.bvox[kVoxF] = ip;
+	dc.xyz[kResY] = y2x*x;
+	//
+	for (int iz=0;iz<nzside;iz++) {
+	  float z2x = dz2x*(0.5+iz);
+	  dc.bvox[kVoxZ] = iz;
+	  if (is>=kNSect) z2x = -z2x;
+	  dc.xyz[kResZ] = z2x*x;
+	  //
+	  fChebDist->Eval(is,x,y2x,z2x,dist);
+	  float xd = x + dist[kResX];
+	  float yd = y2x*x+dist[kResY];
+	  float zd = z2x*x+dist[kResZ];
+	  GetSmoothEstimate(is,xd,yd/xd,zd/xd,(0x1<<4)-1,corr); 
+	  for (int v=kResDim;v--;) {
+	    dc.dxyz[v] = dist[v];
+	    dc.cxyz[v] = corr[v];
+	  }
+	  //
+	  resTree->Fill();
+	}
+	//
+      }
+    }
+  }
+  resTree->Write("", TObject::kOverwrite);
+  delete resTree;
+  flOut->Close();
+  delete flOut;
+}
diff --git a/TPC/TPCcalib/AliTPCDcalibRes.h b/TPC/TPCcalib/AliTPCDcalibRes.h
index f04be168d22..2a58c830e01 100644
--- a/TPC/TPCcalib/AliTPCDcalibRes.h
+++ b/TPC/TPCcalib/AliTPCDcalibRes.h
@@ -113,7 +113,7 @@ class AliTPCDcalibRes: public TNamed
     bres_t() {memset(this,0,sizeof(bres_t));}
   };
  
-  struct delta_t { // structure to organized the input from delta trees
+  struct delta_t { // structure to organize the input from delta trees
     TVectorF *vecDYTRD,*vecDZTRD,*vecDYITS,*vecDZITS,*vecDYTOF,*vecDZTOF,*vecZ,*vecR,*vecSec,*vecPhi;
     AliExternalTrackParam* param;
     Double_t tofBC;
@@ -124,6 +124,16 @@ class AliTPCDcalibRes: public TNamed
     //
     delta_t() {memset(this,0,sizeof(delta_t));}
   };
+  
+  struct dcTst_t {  // struct for distortion/correction check
+    Double32_t xyz[3];
+    Double32_t dxyz[4]; // [-20.,20.,15]
+    Double32_t cxyz[4]; // [-20.,20.,15]
+    UChar_t bsec;          // sector
+    UChar_t bvox[kVoxDim]; // voxel bin info: VoxF,kVoxX,kVoxZ
+    //
+    dcTst_t() {memset(this,0,sizeof(dcTst_t));}
+  };
 
  public:
 
@@ -182,6 +192,7 @@ class AliTPCDcalibRes: public TNamed
   void WriteStatHistos();
   void LoadStatHistos();
   void WriteResTree();
+  void WriteDistCorTestTree(int nx=250, int nphi2sect=30,int nzside=10);
   Bool_t LoadDataFromResTree(const char* resTreeFile);
 
   void  FixAlignmentBug(int sect, float q2pt, float bz, float& alp, float& x, float &z, float &deltaY, float &deltaZ);
diff --git a/TPC/TPCcalib/TPCcalibLinkDef.h b/TPC/TPCcalib/TPCcalibLinkDef.h
index c0e1a9ad085..382be86ef84 100644
--- a/TPC/TPCcalib/TPCcalibLinkDef.h
+++ b/TPC/TPCcalib/TPCcalibLinkDef.h
@@ -56,6 +56,7 @@
 #pragma link C++ class AliTPCDcalibRes::dts_t+;
 #pragma link C++ class AliTPCDcalibRes::dtc_t+;
 #pragma link C++ class AliTPCDcalibRes::dtv_t+;
+#pragma link C++ class AliTPCDcalibRes::dcTst_t+;
 #pragma link C++ class AliTPCDcalibRes::bres_t+;
 #pragma link C++ class AliTPCDcalibRes::delta_t+;