[u/mrichter/AliRoot.git] / TPC / AliTPCROCVoltError3D.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

//////////////////////////////////////////////////////////////////////////////
//                                                                          //
// AliTPCROCVoltError3D class                                               //
// The class calculates the space point distortions due to residual voltage //
// errors on Read Out Chambers of the TPC in 3D.                            //
//                                                                          //
// The class allows "effective Omega Tau" corrections.                      // 
//                                                                          //
// NOTE: This class is capable of calculating z distortions due to          //
//       misalignment and the vd dependency on the residual drift field     //
//                                                                          //
// date: 08/08/2010                                                         //
// Authors: Jim Thomas, Stefan Rossegger                                    //
//                                                                          //
// Example usage :                                                          //
//  AliTPCROCVoltError3D ROCerror;                                            //
//////////////////////////////////////////////////////////////////////////////

#include "AliMagF.h"
#include "TGeoGlobalMagField.h"
#include "AliTPCcalibDB.h"
#include "AliTPCParam.h"
#include "AliLog.h"
#include "TMatrixD.h"
#include "TFile.h"

#include "TMath.h"
#include "AliTPCROC.h"
#include "AliTPCROCVoltError3D.h"

ClassImp(AliTPCROCVoltError3D)

AliTPCROCVoltError3D::AliTPCROCVoltError3D()
  : AliTPCCorrection("ROCVoltErrors","ROC z alignment Errors"),
    fC0(0.),fC1(0.),
    fROCdisplacement(kTRUE),
    fInitLookUp(kFALSE),
    fROCDataFileName(""),  
    fdzDataLinFit(0)
{
  //
  // default constructor
  //

  // Array which will contain the solution according to the setted boundary conditions
  // main input: z alignment of the Read Out chambers
  // see InitROCVoltError3D() function
  for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
    fLookUpErOverEz[k]   =  new TMatrixD(kNR,kNZ);  
    fLookUpEphiOverEz[k] =  new TMatrixD(kNR,kNZ);
    fLookUpDeltaEz[k]    =  new TMatrixD(kNR,kNZ);   
  }
  fROCDataFileName="$ALICE_ROOT/TPC/Calib/maps/TPCROCdzSurvey.root";
  SetROCDataFileName(fROCDataFileName.Data()); // initialization of fdzDataLinFit is included

}

AliTPCROCVoltError3D::~AliTPCROCVoltError3D() {
  //
  // destructor
  //
  
  for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
    delete fLookUpErOverEz[k];
    delete fLookUpEphiOverEz[k];
    delete fLookUpDeltaEz[k];
  }

  delete fdzDataLinFit;
}

void AliTPCROCVoltError3D::SetROCData(TMatrixD * matrix){
  //
  // Set a z alignment map of the chambers not via a file, but
  // directly via a TMatrix(72,3), where dz = p0 + p1*lx + p2*ly
  //
  if (!fdzDataLinFit) fdzDataLinFit=new TMatrixD(*matrix);
  else *fdzDataLinFit = *matrix;
}


void AliTPCROCVoltError3D::Init() {
  //
  // Initialization funtion
  //
  
  AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
  if (!magF) AliError("Magneticd field - not initialized");
  Double_t bzField = magF->SolenoidField()/10.; //field in T
  AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
  if (!param) AliError("Parameters - not initialized");
  Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us]   // From dataBase: to be updated: per second (ideally)
  Double_t ezField = 400; // [V/cm]   // to be updated: never (hopefully)
  Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ; 
  // Correction Terms for effective omegaTau; obtained by a laser calibration run
  SetOmegaTauT1T2(wt,fT1,fT2);

  if (!fInitLookUp) InitROCVoltError3D();
}

void AliTPCROCVoltError3D::Update(const TTimeStamp &/*timeStamp*/) {
  //
  // Update function 
  //
  AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
  if (!magF) AliError("Magneticd field - not initialized");
  Double_t bzField = magF->SolenoidField()/10.; //field in T
  AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
  if (!param) AliError("Parameters - not initialized");
  Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us]  // From dataBase: to be updated: per second (ideally)
  Double_t ezField = 400; // [V/cm]   // to be updated: never (hopefully)
  Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ; 
  // Correction Terms for effective omegaTau; obtained by a laser calibration run
  SetOmegaTauT1T2(wt,fT1,fT2);

}

void  AliTPCROCVoltError3D::SetROCDataFileName(const char * fname) {
  //
  // Set / load the ROC data (linear fit of ROC misalignments)
  //

  fROCDataFileName = fname;
  
  TFile f(fROCDataFileName.Data(),"READ");
  TMatrixD *m = (TMatrixD*) f.Get("dzSurveyLinFitData");
  TMatrixD &mf = *m;

  // prepare some space

  if (fdzDataLinFit) delete fdzDataLinFit;
  fdzDataLinFit = new TMatrixD(72,3);
  TMatrixD &dataIntern = *fdzDataLinFit;
  
  for (Int_t iroc=0;iroc<72;iroc++) {
    dataIntern(iroc,0) = mf(iroc,0);  // z0 offset
    dataIntern(iroc,1) = mf(iroc,1);  // slope in x
    dataIntern(iroc,2) = mf(iroc,2);  // slope in y
  }

  f.Close();

  fInitLookUp = kFALSE;

}

void AliTPCROCVoltError3D::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {
  //
  // Calculates the correction due e.g. residual voltage errors on the TPC boundaries
  //   

  if (!fInitLookUp) {
    AliInfo("Lookup table was not initialized! Perform the inizialisation now ...");
    InitROCVoltError3D();
    return;
  }

  Int_t   order     = 1 ;               // FIXME: hardcoded? Linear interpolation = 1, Quadratic = 2         

  Double_t intEr, intEphi, intDeltaEz;
  Double_t r, phi, z ;
  Int_t    sign;

  r      =  TMath::Sqrt( x[0]*x[0] + x[1]*x[1] ) ;
  phi    =  TMath::ATan2(x[1],x[0]) ;
  if ( phi < 0 ) phi += TMath::TwoPi() ;                   // Table uses phi from 0 to 2*Pi
  z      =  x[2] ;                                         // Create temporary copy of x[2]

  if ( (roc%36) < 18 ) {
    sign =  1;       // (TPC A side)
  } else {
    sign = -1;       // (TPC C side)
  }
  
  if ( sign==1  && z <  fgkZOffSet ) z =  fgkZOffSet;    // Protect against discontinuity at CE
  if ( sign==-1 && z > -fgkZOffSet ) z = -fgkZOffSet;    // Protect against discontinuity at CE
  

  if ( (sign==1 && z<0) || (sign==-1 && z>0) ) // just a consistency check
    AliError("ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!");

  // Get the Er and Ephi field integrals plus the integral over DeltaEz 
  intEr      = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi, 
				  fgkRList, fgkZList, fgkPhiList, fLookUpErOverEz  );
  intEphi    = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi, 
				  fgkRList, fgkZList, fgkPhiList, fLookUpEphiOverEz);
  intDeltaEz = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi, 
				  fgkRList, fgkZList, fgkPhiList, fLookUpDeltaEz   );

  //  printf("%lf %lf %lf\n",intEr,intEphi,intDeltaEz);

  // Calculate distorted position
  if ( r > 0.0 ) {
    phi =  phi + ( fC0*intEphi - fC1*intEr ) / r;      
    r   =  r   + ( fC0*intEr   + fC1*intEphi );  
  }
  
  // Calculate correction in cartesian coordinates
  dx[0] = r * TMath::Cos(phi) - x[0];
  dx[1] = r * TMath::Sin(phi) - x[1]; 
  dx[2] = intDeltaEz;  // z distortion - (internally scaled with driftvelocity dependency 
                       // on the Ez field plus the actual ROC misalignment (if set TRUE)

}

void AliTPCROCVoltError3D::InitROCVoltError3D() {
  //
  // Initialization of the Lookup table which contains the solutions of the 
  // Dirichlet boundary problem
  // Calculation of the single 3D-Poisson solver is done just if needed
  // (see basic lookup tables in header file)
  //

  const Int_t   order       =    1  ;  // Linear interpolation = 1, Quadratic = 2  
  const Float_t gridSizeR   =  (fgkOFCRadius-fgkIFCRadius) / (kRows-1) ;
  const Float_t gridSizeZ   =  fgkTPCZ0 / (kColumns-1) ;
  const Float_t gridSizePhi =  TMath::TwoPi() / ( 18.0 * kPhiSlicesPerSector);

  // temporary arrays to create the boundary conditions
  TMatrixD *arrayofArrayV[kPhiSlices], *arrayofCharge[kPhiSlices] ; 
  TMatrixD *arrayofEroverEz[kPhiSlices], *arrayofEphioverEz[kPhiSlices], *arrayofDeltaEz[kPhiSlices] ; 

  for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
    arrayofArrayV[k]     =   new TMatrixD(kRows,kColumns) ;
    arrayofCharge[k]     =   new TMatrixD(kRows,kColumns) ;
    arrayofEroverEz[k]   =   new TMatrixD(kRows,kColumns) ;
    arrayofEphioverEz[k] =   new TMatrixD(kRows,kColumns) ;
    arrayofDeltaEz[k]    =   new TMatrixD(kRows,kColumns) ;
  }
  
  // list of point as used in the poisson relation and the interpolation (during sum up)
  Double_t  rlist[kRows], zedlist[kColumns] , philist[kPhiSlices];
  for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
    philist[k] =  gridSizePhi * k;
    for ( Int_t i = 0 ; i < kRows ; i++ )    {
      rlist[i] = fgkIFCRadius + i*gridSizeR ;
      for ( Int_t j = 0 ; j < kColumns ; j++ ) { // Fill Vmatrix with Boundary Conditions
	zedlist[j]  = j * gridSizeZ ;
      }
    }
  }

  // ==========================================================================
  // Solve Poisson's equation in 3D cylindrical coordinates by relaxation technique
  // Allow for different size grid spacing in R and Z directions
  
  const Int_t   symmetry = 0;
 
  // Set bondaries and solve Poisson's equation --------------------------
  
  if ( !fInitLookUp ) {
    
    AliInfo(Form("Solving the poisson equation (~ %d sec)",2*10*(int)(kPhiSlices/10)));
    
    for ( Int_t side = 0 ; side < 2 ; side++ ) {  // Solve Poisson3D twice; once for +Z and once for -Z
      
      for ( Int_t k = 0 ; k < kPhiSlices ; k++ )  {
	TMatrixD &arrayV    =  *arrayofArrayV[k] ;
	TMatrixD &charge    =  *arrayofCharge[k] ;
	
	//Fill arrays with initial conditions.  V on the boundary and Charge in the volume.
	for ( Int_t i = 0 ; i < kRows ; i++ ) {
	  for ( Int_t j = 0 ; j < kColumns ; j++ ) {  // Fill Vmatrix with Boundary Conditions
	    arrayV(i,j) = 0.0 ; 
	    charge(i,j) = 0.0 ;

	    Float_t radius0 = rlist[i] ;
	    Float_t phi0    = gridSizePhi * k ;
	    
	    // To avoid problems at sector boundaries, use an average of +- 1 degree from actual phi location
	    if ( j == (kColumns-1) ) 
	      arrayV(i,j) = 0.5*  ( GetROCVoltOffset( side, radius0, phi0+0.02 ) + GetROCVoltOffset( side, radius0, phi0-0.02 ) ) ;

	  }
	}      
	
	for ( Int_t i = 1 ; i < kRows-1 ; i++ ) { 
	  for ( Int_t j = 1 ; j < kColumns-1 ; j++ ) {
	    charge(i,j)  =  0.0 ;
	  }
	}
      }      
      
      // Solve Poisson's equation in 3D cylindrical coordinates by relaxation technique
      // Allow for different size grid spacing in R and Z directions
      
      PoissonRelaxation3D( arrayofArrayV, arrayofCharge, 
			   arrayofEroverEz, arrayofEphioverEz, arrayofDeltaEz,
			   kRows, kColumns, kPhiSlices, gridSizePhi, kIterations, 
			   symmetry, fROCdisplacement) ;
      
      
      //Interpolate results onto a custom grid which is used just for these calculations.
      Double_t  r, phi, z ;
      for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
	phi = fgkPhiList[k] ;
	
	TMatrixD &erOverEz   =  *fLookUpErOverEz[k]  ;
	TMatrixD &ephiOverEz =  *fLookUpEphiOverEz[k];
	TMatrixD &deltaEz    =  *fLookUpDeltaEz[k]   ;
	
	for ( Int_t j = 0 ; j < kNZ ; j++ ) {

	  z = TMath::Abs(fgkZList[j]) ;  // Symmetric solution in Z that depends only on ABS(Z)
  
	  if ( side == 0 &&  fgkZList[j] < 0 ) continue; // Skip rest of this loop if on the wrong side
	  if ( side == 1 &&  fgkZList[j] > 0 ) continue; // Skip rest of this loop if on the wrong side
	  
	  for ( Int_t i = 0 ; i < kNR ; i++ ) { 
	    r = fgkRList[i] ;

	    // Interpolate basicLookup tables; once for each rod, then sum the results
	    erOverEz(i,j)   = Interpolate3DTable(order, r, z, phi, kRows, kColumns, kPhiSlices, 
						 rlist, zedlist, philist, arrayofEroverEz  );
	    ephiOverEz(i,j) = Interpolate3DTable(order, r, z, phi, kRows, kColumns, kPhiSlices,
						 rlist, zedlist, philist, arrayofEphioverEz);
	    deltaEz(i,j)    = Interpolate3DTable(order, r, z, phi, kRows, kColumns, kPhiSlices,
						 rlist, zedlist, philist, arrayofDeltaEz  );

	    if (side == 1)  deltaEz(i,j) = -  deltaEz(i,j); // negative coordinate system on C side

	  } // end r loop
	}// end z loop
      }// end phi loop

      if ( side == 0 ) AliInfo(" A side done");
      if ( side == 1 ) AliInfo(" C side done");
    } // end side loop
  }
  
  // clear the temporary arrays lists
  for ( Int_t k = 0 ; k < kPhiSlices ; k++ )  {
    delete arrayofArrayV[k];
    delete arrayofCharge[k];
    delete arrayofEroverEz[k];  
    delete arrayofEphioverEz[k];
    delete arrayofDeltaEz[k];
  }
 

  fInitLookUp = kTRUE;

}


Float_t AliTPCROCVoltError3D::GetROCVoltOffset(Int_t side, Float_t r0, Float_t phi0) {
  // 
  // Returns the dz alignment data (in voltage equivalents) at 
  // the given position
  //

  Float_t xp = r0*TMath::Cos(phi0);
  Float_t yp = r0*TMath::Sin(phi0);
  
  // phi0 should be between 0 and 2pi 
  if (phi0<0) phi0+=TMath::TwoPi();
  Int_t roc = (Int_t)TMath::Floor((TMath::RadToDeg()*phi0)/20);
  if (side==1) roc+=18; // C side
  if (r0>132) roc+=36;  // OROC 
  
  // linear-plane data:  z = z0 + kx*lx + ky*ly (rotation in local coordinates)
  TMatrixD &fitData = *fdzDataLinFit;

  // local coordinates                                                          
  Double_t secAlpha = TMath::DegToRad()*(10.+20.*(((Int_t)roc)%18));
  Float_t lx = xp*TMath::Cos(secAlpha)+yp*TMath::Sin(secAlpha);
  Float_t ly = yp*TMath::Cos(secAlpha)-xp*TMath::Sin(secAlpha);

  Float_t dz = fitData(roc,0)+fitData(roc,1)*lx + fitData(roc,2)*ly; // value in cm

  // aproximated Voltage-offset-aquivalent to the z misalignment
  // (linearly scaled with the z position)
  Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm; 
  Float_t voltOff = dz*ezField;            // values in "Volt equivalents"

  return voltOff;
}

TH2F * AliTPCROCVoltError3D::CreateHistoOfZSurvey(Int_t side, Int_t nx, Int_t ny) {
  //
  // return a simple histogramm containing the input to the poisson solver
  // (z positions of the Read-out chambers, linearly interpolated)

  char hname[100];
  if (side==0) sprintf(hname,"survey_dz_Aside");
  if (side==1) sprintf(hname,"survey_dz_Cside");

  TH2F *h = new TH2F(hname,hname,nx,-250.,250.,ny,-250.,250.);

  for (Int_t iy=1;iy<=ny;++iy) {
    Double_t yp = h->GetYaxis()->GetBinCenter(iy);
    for (Int_t ix=1;ix<=nx;++ix) {
      Double_t xp = h->GetXaxis()->GetBinCenter(ix);
    
      Float_t r0 = TMath::Sqrt(xp*xp+yp*yp);
      Float_t phi0 = TMath::ATan2(yp,xp); 
   
      Float_t dz = GetROCVoltOffset(side,r0,phi0); // in [volt]

      Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm; 
      dz = dz/ezField;    // in [cm]

      if (85.<=r0 && r0<=245.) {
	h->SetBinContent(ix,iy,dz); 
      } else {
	h->SetBinContent(ix,iy,0.);
      }
    }
  }
  
  h->GetXaxis()->SetTitle("x [cm]");
  h->GetYaxis()->SetTitle("y [cm]");
  h->GetZaxis()->SetTitle("dz [cm]");
  h->SetStats(0);
  //  h->DrawCopy("colz");

  return h;
} 

void AliTPCROCVoltError3D::Print(const Option_t* option) const {
  //
  // Print function to check the settings of the Rod shifts and the rotated clips
  // option=="a" prints the C0 and C1 coefficents for calibration purposes
  //

  TString opt = option; opt.ToLower();
  printf("%s\n",GetTitle());
  printf(" - Voltage settings on the TPC Read-Out chambers - linearly interpolated\n");
  printf("   info: Check the following data-file for more details: %s \n",fROCDataFileName.Data());

  if (opt.Contains("a")) { // Print all details
    printf(" - T1: %1.4f, T2: %1.4f \n",fT1,fT2);
    printf(" - C1: %1.4f, C0: %1.4f \n",fC1,fC0);
  }

  if (!fInitLookUp) AliError("Lookup table was not initialized! You should do InitROCVoltError3D() ...");

}
Commit	Line	Data
c9cbd2f2	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	//////////////////////////////////////////////////////////////////////////////
	17	// //
	18	// AliTPCROCVoltError3D class //
	19	// The class calculates the space point distortions due to residual voltage //
	20	// errors on Read Out Chambers of the TPC in 3D. //
	21	// //
	22	// The class allows "effective Omega Tau" corrections. //
	23	// //
	24	// NOTE: This class is capable of calculating z distortions due to //
	25	// misalignment and the vd dependency on the residual drift field //
	26	// //
	27	// date: 08/08/2010 //
	28	// Authors: Jim Thomas, Stefan Rossegger //
	29	// //
	30	// Example usage : //
	31	// AliTPCROCVoltError3D ROCerror; //
	32	//////////////////////////////////////////////////////////////////////////////
	33
	34	#include "AliMagF.h"
	35	#include "TGeoGlobalMagField.h"
	36	#include "AliTPCcalibDB.h"
	37	#include "AliTPCParam.h"
	38	#include "AliLog.h"
	39	#include "TMatrixD.h"
	40	#include "TFile.h"
	41
	42	#include "TMath.h"
	43	#include "AliTPCROC.h"
	44	#include "AliTPCROCVoltError3D.h"
	45
	46	ClassImp(AliTPCROCVoltError3D)
	47
	48	AliTPCROCVoltError3D::AliTPCROCVoltError3D()
	49	: AliTPCCorrection("ROCVoltErrors","ROC z alignment Errors"),
	50	fC0(0.),fC1(0.),
	51	fROCdisplacement(kTRUE),
	52	fInitLookUp(kFALSE),
acf5907b	53	fROCDataFileName(""),
c9cbd2f2	54	fdzDataLinFit(0)
	55	{
	56	//
	57	// default constructor
	58	//
	59
	60	// Array which will contain the solution according to the setted boundary conditions
	61	// main input: z alignment of the Read Out chambers
	62	// see InitROCVoltError3D() function
	63	for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
	64	fLookUpErOverEz[k] = new TMatrixD(kNR,kNZ);
	65	fLookUpEphiOverEz[k] = new TMatrixD(kNR,kNZ);
	66	fLookUpDeltaEz[k] = new TMatrixD(kNR,kNZ);
	67	}
acf5907b	68	fROCDataFileName="$ALICE_ROOT/TPC/Calib/maps/TPCROCdzSurvey.root";
acf5907b	69	SetROCDataFileName(fROCDataFileName.Data()); // initialization of fdzDataLinFit is included
c9cbd2f2	70
	71	}
	72
	73	AliTPCROCVoltError3D::~AliTPCROCVoltError3D() {
	74	//
	75	// destructor
	76	//
	77
	78	for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
	79	delete fLookUpErOverEz[k];
	80	delete fLookUpEphiOverEz[k];
	81	delete fLookUpDeltaEz[k];
	82	}
	83
	84	delete fdzDataLinFit;
	85	}
	86
acf5907b	87	void AliTPCROCVoltError3D::SetROCData(TMatrixD * matrix){
2bbac918	88	//
c756f562	89	// Set a z alignment map of the chambers not via a file, but
c756f562	90	// directly via a TMatrix(72,3), where dz = p0 + p1lx + p2ly
2bbac918	91	//
	92	if (!fdzDataLinFit) fdzDataLinFit=new TMatrixD(*matrix);
	93	else fdzDataLinFit = matrix;
	94	}
	95
	96
c9cbd2f2	97	void AliTPCROCVoltError3D::Init() {
	98	//
	99	// Initialization funtion
	100	//
	101
	102	AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
	103	if (!magF) AliError("Magneticd field - not initialized");
	104	Double_t bzField = magF->SolenoidField()/10.; //field in T
	105	AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
	106	if (!param) AliError("Parameters - not initialized");
	107	Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us] // From dataBase: to be updated: per second (ideally)
	108	Double_t ezField = 400; // [V/cm] // to be updated: never (hopefully)
	109	Double_t wt = -10.0 * (bzField10) vdrift / ezField ;
	110	// Correction Terms for effective omegaTau; obtained by a laser calibration run
	111	SetOmegaTauT1T2(wt,fT1,fT2);
	112
35ae345f	113	if (!fInitLookUp) InitROCVoltError3D();
c9cbd2f2	114	}
	115
	116	void AliTPCROCVoltError3D::Update(const TTimeStamp &/timeStamp/) {
	117	//
	118	// Update function
	119	//
	120	AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
	121	if (!magF) AliError("Magneticd field - not initialized");
	122	Double_t bzField = magF->SolenoidField()/10.; //field in T
	123	AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
	124	if (!param) AliError("Parameters - not initialized");
	125	Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us] // From dataBase: to be updated: per second (ideally)
	126	Double_t ezField = 400; // [V/cm] // to be updated: never (hopefully)
	127	Double_t wt = -10.0 * (bzField10) vdrift / ezField ;
	128	// Correction Terms for effective omegaTau; obtained by a laser calibration run
	129	SetOmegaTauT1T2(wt,fT1,fT2);
	130
	131	}
	132
acf5907b	133	void AliTPCROCVoltError3D::SetROCDataFileName(const char * fname) {
c9cbd2f2	134	//
	135	// Set / load the ROC data (linear fit of ROC misalignments)
	136	//
	137
	138	fROCDataFileName = fname;
	139
acf5907b	140	TFile f(fROCDataFileName.Data(),"READ");
c9cbd2f2	141	TMatrixD m = (TMatrixD) f.Get("dzSurveyLinFitData");
	142	TMatrixD &mf = *m;
	143
	144	// prepare some space
	145
	146	if (fdzDataLinFit) delete fdzDataLinFit;
	147	fdzDataLinFit = new TMatrixD(72,3);
	148	TMatrixD &dataIntern = *fdzDataLinFit;
	149
	150	for (Int_t iroc=0;iroc<72;iroc++) {
	151	dataIntern(iroc,0) = mf(iroc,0); // z0 offset
	152	dataIntern(iroc,1) = mf(iroc,1); // slope in x
	153	dataIntern(iroc,2) = mf(iroc,2); // slope in y
	154	}
	155
	156	f.Close();
	157
	158	fInitLookUp = kFALSE;
	159
	160	}
	161
	162	void AliTPCROCVoltError3D::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {
	163	//
	164	// Calculates the correction due e.g. residual voltage errors on the TPC boundaries
	165	//
	166
	167	if (!fInitLookUp) {
	168	AliInfo("Lookup table was not initialized! Perform the inizialisation now ...");
	169	InitROCVoltError3D();
	170	return;
	171	}
	172
	173	Int_t order = 1 ; // FIXME: hardcoded? Linear interpolation = 1, Quadratic = 2
	174
	175	Double_t intEr, intEphi, intDeltaEz;
	176	Double_t r, phi, z ;
	177	Int_t sign;
	178
	179	r = TMath::Sqrt( x[0]x[0] + x[1]x[1] ) ;
	180	phi = TMath::ATan2(x[1],x[0]) ;
	181	if ( phi < 0 ) phi += TMath::TwoPi() ; // Table uses phi from 0 to 2*Pi
	182	z = x[2] ; // Create temporary copy of x[2]
	183
	184	if ( (roc%36) < 18 ) {
	185	sign = 1; // (TPC A side)
	186	} else {
	187	sign = -1; // (TPC C side)
	188	}
	189
	190	if ( sign==1 && z < fgkZOffSet ) z = fgkZOffSet; // Protect against discontinuity at CE
	191	if ( sign==-1 && z > -fgkZOffSet ) z = -fgkZOffSet; // Protect against discontinuity at CE
	192
	193
	194	if ( (sign==1 && z<0) \|\| (sign==-1 && z>0) ) // just a consistency check
	195	AliError("ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!");
	196
	197	// Get the Er and Ephi field integrals plus the integral over DeltaEz
	198	intEr = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
	199	fgkRList, fgkZList, fgkPhiList, fLookUpErOverEz );
	200	intEphi = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
	201	fgkRList, fgkZList, fgkPhiList, fLookUpEphiOverEz);
	202	intDeltaEz = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
	203	fgkRList, fgkZList, fgkPhiList, fLookUpDeltaEz );
	204
205	// printf("%lf %lf %lf\n",intEr,intEphi,intDeltaEz);
206
207	// Calculate distorted position
208	if ( r > 0.0 ) {
209	phi = phi + ( fC0intEphi - fC1intEr ) / r;
210	r = r + ( fC0intEr + fC1intEphi );
211	}
212
213	// Calculate correction in cartesian coordinates
214	dx[0] = r * TMath::Cos(phi) - x[0];
215	dx[1] = r * TMath::Sin(phi) - x[1];
216	dx[2] = intDeltaEz; // z distortion - (internally scaled with driftvelocity dependency
217	// on the Ez field plus the actual ROC misalignment (if set TRUE)
218
219	}
220
221	void AliTPCROCVoltError3D::InitROCVoltError3D() {
222	//
223	// Initialization of the Lookup table which contains the solutions of the
224	// Dirichlet boundary problem
225	// Calculation of the single 3D-Poisson solver is done just if needed
226	// (see basic lookup tables in header file)
227	//
228
229	const Int_t order = 1 ; // Linear interpolation = 1, Quadratic = 2
230	const Float_t gridSizeR = (fgkOFCRadius-fgkIFCRadius) / (kRows-1) ;
231	const Float_t gridSizeZ = fgkTPCZ0 / (kColumns-1) ;
232	const Float_t gridSizePhi = TMath::TwoPi() / ( 18.0 * kPhiSlicesPerSector);
233
234	// temporary arrays to create the boundary conditions
235	TMatrixD arrayofArrayV[kPhiSlices], arrayofCharge[kPhiSlices] ;
236	TMatrixD arrayofEroverEz[kPhiSlices], arrayofEphioverEz[kPhiSlices], *arrayofDeltaEz[kPhiSlices] ;
237
238	for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
239	arrayofArrayV[k] = new TMatrixD(kRows,kColumns) ;
240	arrayofCharge[k] = new TMatrixD(kRows,kColumns) ;
241	arrayofEroverEz[k] = new TMatrixD(kRows,kColumns) ;
242	arrayofEphioverEz[k] = new TMatrixD(kRows,kColumns) ;
243	arrayofDeltaEz[k] = new TMatrixD(kRows,kColumns) ;
244	}
245
246	// list of point as used in the poisson relation and the interpolation (during sum up)
247	Double_t rlist[kRows], zedlist[kColumns] , philist[kPhiSlices];
248	for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
249	philist[k] = gridSizePhi * k;
250	for ( Int_t i = 0 ; i < kRows ; i++ ) {
251	rlist[i] = fgkIFCRadius + i*gridSizeR ;
252	for ( Int_t j = 0 ; j < kColumns ; j++ ) { // Fill Vmatrix with Boundary Conditions
253	zedlist[j] = j * gridSizeZ ;
254	}
255	}
256	}
257
258	// ==========================================================================
259	// Solve Poisson's equation in 3D cylindrical coordinates by relaxation technique
260	// Allow for different size grid spacing in R and Z directions
261
262	const Int_t symmetry = 0;
263
264	// Set bondaries and solve Poisson's equation --------------------------
265
266	if ( !fInitLookUp ) {
267
268	AliInfo(Form("Solving the poisson equation (~ %d sec)",210(int)(kPhiSlices/10)));
269
270	for ( Int_t side = 0 ; side < 2 ; side++ ) { // Solve Poisson3D twice; once for +Z and once for -Z
271
272	for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
273	TMatrixD &arrayV = *arrayofArrayV[k] ;
274	TMatrixD &charge = *arrayofCharge[k] ;
275
276	//Fill arrays with initial conditions. V on the boundary and Charge in the volume.
277	for ( Int_t i = 0 ; i < kRows ; i++ ) {
278	for ( Int_t j = 0 ; j < kColumns ; j++ ) { // Fill Vmatrix with Boundary Conditions
279	arrayV(i,j) = 0.0 ;
280	charge(i,j) = 0.0 ;
281
282	Float_t radius0 = rlist[i] ;
283	Float_t phi0 = gridSizePhi * k ;
284
285	// To avoid problems at sector boundaries, use an average of +- 1 degree from actual phi location
286	if ( j == (kColumns-1) )
287	arrayV(i,j) = 0.5* ( GetROCVoltOffset( side, radius0, phi0+0.02 ) + GetROCVoltOffset( side, radius0, phi0-0.02 ) ) ;
288
289	}
290	}
291
292	for ( Int_t i = 1 ; i < kRows-1 ; i++ ) {
293	for ( Int_t j = 1 ; j < kColumns-1 ; j++ ) {
294	charge(i,j) = 0.0 ;
295	}
296	}
297	}
298
299	// Solve Poisson's equation in 3D cylindrical coordinates by relaxation technique
300	// Allow for different size grid spacing in R and Z directions
301
302	PoissonRelaxation3D( arrayofArrayV, arrayofCharge,
303	arrayofEroverEz, arrayofEphioverEz, arrayofDeltaEz,
304	kRows, kColumns, kPhiSlices, gridSizePhi, kIterations,
305	symmetry, fROCdisplacement) ;
306
307
308	//Interpolate results onto a custom grid which is used just for these calculations.
309	Double_t r, phi, z ;
310	for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
311	phi = fgkPhiList[k] ;
312
313	TMatrixD &erOverEz = *fLookUpErOverEz[k] ;
314	TMatrixD &ephiOverEz = *fLookUpEphiOverEz[k];
315	TMatrixD &deltaEz = *fLookUpDeltaEz[k] ;
316
317	for ( Int_t j = 0 ; j < kNZ ; j++ ) {
318
319	z = TMath::Abs(fgkZList[j]) ; // Symmetric solution in Z that depends only on ABS(Z)
320
321	if ( side == 0 && fgkZList[j] < 0 ) continue; // Skip rest of this loop if on the wrong side
322	if ( side == 1 && fgkZList[j] > 0 ) continue; // Skip rest of this loop if on the wrong side
323
324	for ( Int_t i = 0 ; i < kNR ; i++ ) {
325	r = fgkRList[i] ;
326
327	// Interpolate basicLookup tables; once for each rod, then sum the results
328	erOverEz(i,j) = Interpolate3DTable(order, r, z, phi, kRows, kColumns, kPhiSlices,
329	rlist, zedlist, philist, arrayofEroverEz );
330	ephiOverEz(i,j) = Interpolate3DTable(order, r, z, phi, kRows, kColumns, kPhiSlices,
331	rlist, zedlist, philist, arrayofEphioverEz);
332	deltaEz(i,j) = Interpolate3DTable(order, r, z, phi, kRows, kColumns, kPhiSlices,
333	rlist, zedlist, philist, arrayofDeltaEz );
334
335	if (side == 1) deltaEz(i,j) = - deltaEz(i,j); // negative coordinate system on C side
336
337	} // end r loop
338	}// end z loop
339	}// end phi loop
340
341	if ( side == 0 ) AliInfo(" A side done");
342	if ( side == 1 ) AliInfo(" C side done");
343	} // end side loop
344	}
345
346	// clear the temporary arrays lists
347	for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
348	delete arrayofArrayV[k];
349	delete arrayofCharge[k];
350	delete arrayofEroverEz[k];
351	delete arrayofEphioverEz[k];
352	delete arrayofDeltaEz[k];
353	}
354
355
356	fInitLookUp = kTRUE;
357
358	}
359
360
361	Float_t AliTPCROCVoltError3D::GetROCVoltOffset(Int_t side, Float_t r0, Float_t phi0) {
362	//
363	// Returns the dz alignment data (in voltage equivalents) at
364	// the given position
365	//
366
367	Float_t xp = r0*TMath::Cos(phi0);
368	Float_t yp = r0*TMath::Sin(phi0);
369
370	// phi0 should be between 0 and 2pi
371	if (phi0<0) phi0+=TMath::TwoPi();
372	Int_t roc = (Int_t)TMath::Floor((TMath::RadToDeg()*phi0)/20);
373	if (side==1) roc+=18; // C side
374	if (r0>132) roc+=36; // OROC
375
c756f562	376	// linear-plane data: z = z0 + kxlx + kyly (rotation in local coordinates)
c9cbd2f2	377	TMatrixD &fitData = *fdzDataLinFit;
c756f562	378
	379	// local coordinates
	380	Double_t secAlpha = TMath::DegToRad()(10.+20.(((Int_t)roc)%18));
	381	Float_t lx = xpTMath::Cos(secAlpha)+ypTMath::Sin(secAlpha);
	382	Float_t ly = ypTMath::Cos(secAlpha)-xpTMath::Sin(secAlpha);
	383
	384	Float_t dz = fitData(roc,0)+fitData(roc,1)lx + fitData(roc,2)ly; // value in cm
c9cbd2f2	385
	386	// aproximated Voltage-offset-aquivalent to the z misalignment
	387	// (linearly scaled with the z position)
	388	Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm;
	389	Float_t voltOff = dz*ezField; // values in "Volt equivalents"
	390
	391	return voltOff;
	392	}
	393
	394	TH2F * AliTPCROCVoltError3D::CreateHistoOfZSurvey(Int_t side, Int_t nx, Int_t ny) {
	395	//
	396	// return a simple histogramm containing the input to the poisson solver
	397	// (z positions of the Read-out chambers, linearly interpolated)
	398
	399	char hname[100];
	400	if (side==0) sprintf(hname,"survey_dz_Aside");
	401	if (side==1) sprintf(hname,"survey_dz_Cside");
	402
	403	TH2F *h = new TH2F(hname,hname,nx,-250.,250.,ny,-250.,250.);
	404
	405	for (Int_t iy=1;iy<=ny;++iy) {
	406	Double_t yp = h->GetYaxis()->GetBinCenter(iy);
	407	for (Int_t ix=1;ix<=nx;++ix) {
	408	Double_t xp = h->GetXaxis()->GetBinCenter(ix);
	409
	410	Float_t r0 = TMath::Sqrt(xpxp+ypyp);
	411	Float_t phi0 = TMath::ATan2(yp,xp);
	412
	413	Float_t dz = GetROCVoltOffset(side,r0,phi0); // in [volt]
	414
	415	Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm;
	416	dz = dz/ezField; // in [cm]
	417
	418	if (85.<=r0 && r0<=245.) {
	419	h->SetBinContent(ix,iy,dz);
	420	} else {
	421	h->SetBinContent(ix,iy,0.);
	422	}
	423	}
	424	}
	425
	426	h->GetXaxis()->SetTitle("x [cm]");
	427	h->GetYaxis()->SetTitle("y [cm]");
	428	h->GetZaxis()->SetTitle("dz [cm]");
	429	h->SetStats(0);
	430	// h->DrawCopy("colz");
	431
	432	return h;
	433	}
	434
	435	void AliTPCROCVoltError3D::Print(const Option_t* option) const {
	436	//
	437	// Print function to check the settings of the Rod shifts and the rotated clips
	438	// option=="a" prints the C0 and C1 coefficents for calibration purposes
	439	//
	440
	441	TString opt = option; opt.ToLower();
	442	printf("%s\n",GetTitle());
	443	printf(" - Voltage settings on the TPC Read-Out chambers - linearly interpolated\n");
acf5907b	444	printf(" info: Check the following data-file for more details: %s \n",fROCDataFileName.Data());
c9cbd2f2	445
	446	if (opt.Contains("a")) { // Print all details
	447	printf(" - T1: %1.4f, T2: %1.4f \n",fT1,fT2);
	448	printf(" - C1: %1.4f, C0: %1.4f \n",fC1,fC0);
	449	}
	450
	451	if (!fInitLookUp) AliError("Lookup table was not initialized! You should do InitROCVoltError3D() ...");
	452
	453	}