[u/mrichter/AliRoot.git] / TPC / AliTPCGGVoltError.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.                  *
 **************************************************************************/

////////////////////////////////////////////////////////////////////////////
//                                                                        //
// AliTPCGGVoltError class                                                //
// The class calculates the electric field and space point distortions    //
// due a Gating Grid (GG) Error voltage. It uses the exact calculation    //
// technique based on bessel functions. (original code from STAR)         //
// The class allows "effective Omega Tau" corrections.                    // 
//                                                                        //
// date: 27/04/2010                                                       //
// Authors: Jim Thomas, Stefan Rossegger, Magnus Mager                    //
//                                                                        //
// Example usage:                                                         //
//  AliTPCGGVoltError GGerror;                                            //
//  GGerror.SetOmegaTauT1T2(0.32,1.,1.); // values ideally from OCDB      //
//  GGerror.SetDeltaVGGA(50.);           // voltage offset A-side         //
//  GGerror.SetDeltaVGGC(50.);           // voltage offset C-side         //
//  GGerror.InitGGVoltErrorDistortion(); // initialization of the look up //
//  // plot dRPhi distortions ...                                         //
//  GGerror.CreateHistoDRPhiinZR(1.,100,100)->Draw("surf2");              //
////////////////////////////////////////////////////////////////////////////


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

#include "AliTPCGGVoltError.h"
#include <TMath.h>

AliTPCGGVoltError::AliTPCGGVoltError()
  : AliTPCCorrection("GGVoltError","GatingGrid (GG) Voltage Error"),
    fC0(0.),fC1(0.),
    fDeltaVGGA(0.),fDeltaVGGC(0.),
    fInitLookUp(kFALSE)
{
  //
  // default constructor
  //
}

AliTPCGGVoltError::~AliTPCGGVoltError() {
  //
  // default destructor
  //
}

void AliTPCGGVoltError::Init() {
  //
  // Init 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 ; 
  //
  SetOmegaTauT1T2(wt,fT1,fT2);
  InitGGVoltErrorDistortion();
  //SetDeltaVGGA(0.0);//  ideally from the database
  //SetDeltaVGGC(0.0);//  ideally from the database
}

void AliTPCGGVoltError::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 ; 

  SetOmegaTauT1T2(wt,fT1,fT2);
  //  InitGGVoltErrorDistortion(); // not necessary in here since the Voltage should not change!
}


void AliTPCGGVoltError::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {

  //
  // Gated Grid Voltage Error
  //
  // Calculates the effect of having an incorrect voltage on the A or C end plate Gated Grids.
  //
  // Electrostatic Equations from StarNote SN0253 by Howard Wieman.
  //
  
  if (!fInitLookUp) AliError("Lookup table was not initialized! You should do InitGGVoltErrorDistortion() ...");
  
  Int_t   order     = 1 ;               // FIXME: hardcoded? Linear interpolation = 1, Quadratic = 2         
 
  Double_t intEr, intEphi ;
  Double_t r, phi, z ;
  Int_t    sign ;

  Double_t deltaVGG;
  
  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] ;

  if ( (roc%36) < 18 ) {
    sign =  1; 
    deltaVGG = fDeltaVGGA;           // (TPC End A)
  } else {
    sign = -1;                       // (TPC End C)
    deltaVGG = fDeltaVGGC; 
  }

  if ( sign==1  && z <  fgkZOffSet ) z =  fgkZOffSet;    // Protect against discontinuity at CE
  if ( sign==-1 && z > -fgkZOffSet ) z = -fgkZOffSet;    // Protect against discontinuity at CE

  Interpolate2DEdistortion( order, r, z, fGGVoltErrorER, intEr );
  intEphi = 0.0;  // Efield is symmetric in phi

  // Calculate distorted position
  if ( r > 0.0 ) {
    phi =  phi + deltaVGG*( fC0*intEphi - fC1*intEr ) / r;      
    r   =  r   + deltaVGG*( 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] = 0.; // z distortion not implemented (1st order distortions) - see e.g. AliTPCBoundaryVoltError-class


}


Float_t AliTPCGGVoltError::GetIntErOverEz(const Float_t x[],const Short_t roc) {
  //
  // This function is purely for calibration purposes
  // Calculates the integral (int Er/Ez dz) for the setted GG voltage offset 
  // 

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

  Int_t   order     = 1 ;     // FIXME: so far hardcoded? Linear interpolation = 1, Quadratic = 2         
  
  Double_t intEr;
  Double_t r, phi, z ;
  Int_t    sign ;
  
  Double_t deltaVGG;
  
  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] ;

  if ( (roc%36) < 18 ) {
    sign =  1; 
    deltaVGG = fDeltaVGGA;           // (TPC End A)
  } else {
    sign = -1;                       // (TPC End C)
    deltaVGG = fDeltaVGGC; 
  }

  if ( sign==1  && z <  fgkZOffSet ) z =  fgkZOffSet;    // Protect against discontinuity at CE
  if ( sign==-1 && z > -fgkZOffSet ) z = -fgkZOffSet;    // Protect against discontinuity at CE

  Interpolate2DEdistortion(order, r, z, fGGVoltErrorER, intEr );

  return (intEr*deltaVGG);

}

void AliTPCGGVoltError::InitGGVoltErrorDistortion() {
  //
  // Initialization of the Lookup table which contains the solutions of the GG Error problem
  //

  Double_t r,z;
  Int_t nterms = 100 ;
  for ( Int_t i = 0 ; i < kNZ ; ++i ) {
    z = fgkZList[i] ;
    for ( Int_t j = 0 ; j < kNR ; ++j ) {
      r = fgkRList[j] ;
      fGGVoltErrorER[i][j] = 0.0 ; 	    
      Double_t intz = 0.0 ;
      for ( Int_t n = 1 ; n < nterms ; ++n ) {
	Double_t k    =  n * TMath::Pi() / fgkTPCZ0 ;
	Double_t ein  =  0 ;                    // Error potential on the IFC
	Double_t eout =  0 ;                    // Error potential on the OFC
	if ( z < 0 ) {
	  ein   =  -2.0 / ( k * (fgkCathodeV - fgkGG) ) ;       
	  eout  =  -2.0 / ( k * (fgkCathodeV - fgkGG) ) ;       
	}
	if ( z == 0 ) continue ;
	if ( z > 0 ) {
	  ein   =  -2.0 / ( k * (fgkCathodeV - fgkGG) ) ;       
	  eout  =  -2.0 / ( k * (fgkCathodeV - fgkGG) ) ;       
	}
	Double_t an   =  ein  * TMath::BesselK0( k*fgkOFCRadius ) - eout * TMath::BesselK0( k*fgkIFCRadius ) ;
	Double_t bn   =  eout * TMath::BesselI0( k*fgkIFCRadius ) - ein  * TMath::BesselI0( k*fgkOFCRadius ) ;
	Double_t numerator =
	  an * TMath::BesselI1( k*r ) - bn * TMath::BesselK1( k*r ) ;
	Double_t denominator =
	  TMath::BesselK0( k*fgkOFCRadius ) * TMath::BesselI0( k*fgkIFCRadius ) -
	  TMath::BesselK0( k*fgkIFCRadius ) * TMath::BesselI0( k*fgkOFCRadius ) ;
	Double_t zterm = TMath::Cos( k*(fgkTPCZ0-TMath::Abs(z)) ) - 1 ;
	intz += zterm * numerator / denominator ;
	// Assume series converges, break if small terms
	if ( n>10 && TMath::Abs(intz)*1.e-10 > TMath::Abs(numerator/denominator) ) break;   
      }
      fGGVoltErrorER[i][j] = (Double_t) intz ;

    }
  }
  
  fInitLookUp = kTRUE;
}


void AliTPCGGVoltError::Print(const Option_t* option) const {
  //
  // Print function to check the settings (e.g. voltage offsets)
  // option=="a" prints the C0 and C1 coefficents for calibration purposes
  //

  TString opt = option; opt.ToLower();
  printf("%s\n",GetTitle());
  printf(" - GG Voltage offset: A-side: %3.1f V, C-side: %3.1f V \n",fDeltaVGGA,fDeltaVGGC);  
  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 InitGGVoltErrorDistortion() ...");

  
}
Commit	Line	Data
0116859c	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	// AliTPCGGVoltError class //
	19	// The class calculates the electric field and space point distortions //
	20	// due a Gating Grid (GG) Error voltage. It uses the exact calculation //
	21	// technique based on bessel functions. (original code from STAR) //
	22	// The class allows "effective Omega Tau" corrections. //
	23	// //
	24	// date: 27/04/2010 //
	25	// Authors: Jim Thomas, Stefan Rossegger, Magnus Mager //
	26	// //
	27	// Example usage: //
	28	// AliTPCGGVoltError GGerror; //
	29	// GGerror.SetOmegaTauT1T2(0.32,1.,1.); // values ideally from OCDB //
	30	// GGerror.SetDeltaVGGA(50.); // voltage offset A-side //
	31	// GGerror.SetDeltaVGGC(50.); // voltage offset C-side //
	32	// GGerror.InitGGVoltErrorDistortion(); // initialization of the look up //
	33	// // plot dRPhi distortions ... //
	34	// GGerror.CreateHistoDRPhiinZR(1.,100,100)->Draw("surf2"); //
	35	////////////////////////////////////////////////////////////////////////////
	36
	37
534fd34a	38	#include "AliMagF.h"
	39	#include "TGeoGlobalMagField.h"
	40	#include "AliTPCcalibDB.h"
	41	#include "AliTPCParam.h"
	42	#include "AliLog.h"
0116859c	43
	44	#include "AliTPCGGVoltError.h"
	45	#include <TMath.h>
	46
	47	AliTPCGGVoltError::AliTPCGGVoltError()
	48	: AliTPCCorrection("GGVoltError","GatingGrid (GG) Voltage Error"),
	49	fC0(0.),fC1(0.),
1b923461	50	fDeltaVGGA(0.),fDeltaVGGC(0.),
1b923461	51	fInitLookUp(kFALSE)
0116859c	52	{
	53	//
	54	// default constructor
	55	//
	56	}
	57
	58	AliTPCGGVoltError::~AliTPCGGVoltError() {
	59	//
	60	// default destructor
	61	//
	62	}
	63
e527a1b9	64	void AliTPCGGVoltError::Init() {
e527a1b9	65	//
b1f0a2a5	66	// Init function
534fd34a	67	//
	68	AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
	69	if (!magF) AliError("Magneticd field - not initialized");
	70	Double_t bzField = magF->SolenoidField()/10.; //field in T
	71	AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
	72	if (!param) AliError("Parameters - not initialized");
	73	Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us] // From dataBase: to be updated: per second (ideally)
e527a1b9	74	Double_t ezField = 400; // [V/cm] // to be updated: never (hopefully)
e527a1b9	75	Double_t wt = -10.0 * (bzField10) vdrift / ezField ;
534fd34a	76	//
534fd34a	77	SetOmegaTauT1T2(wt,fT1,fT2);
8b63d99c	78	InitGGVoltErrorDistortion();
534fd34a	79	//SetDeltaVGGA(0.0);// ideally from the database
534fd34a	80	//SetDeltaVGGC(0.0);// ideally from the database
e527a1b9	81	}
	82
	83	void AliTPCGGVoltError::Update(const TTimeStamp &/timeStamp/) {
	84	//
534fd34a	85	// Update function
e527a1b9	86	//
534fd34a	87	AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
	88	if (!magF) AliError("Magneticd field - not initialized");
	89	Double_t bzField = magF->SolenoidField()/10.; //field in T
	90	AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
	91	if (!param) AliError("Parameters - not initialized");
	92	Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us] // From dataBase: to be updated: per second (ideally)
e527a1b9	93	Double_t ezField = 400; // [V/cm] // to be updated: never (hopefully)
	94	Double_t wt = -10.0 * (bzField10) vdrift / ezField ;
	95
534fd34a	96	SetOmegaTauT1T2(wt,fT1,fT2);
1b923461	97	// InitGGVoltErrorDistortion(); // not necessary in here since the Voltage should not change!
e527a1b9	98	}
	99
	100
	101
0116859c	102	void AliTPCGGVoltError::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {
	103
	104	//
	105	// Gated Grid Voltage Error
	106	//
	107	// Calculates the effect of having an incorrect voltage on the A or C end plate Gated Grids.
	108	//
	109	// Electrostatic Equations from StarNote SN0253 by Howard Wieman.
	110	//
	111
1b923461	112	if (!fInitLookUp) AliError("Lookup table was not initialized! You should do InitGGVoltErrorDistortion() ...");
1b923461	113
0116859c	114	Int_t order = 1 ; // FIXME: hardcoded? Linear interpolation = 1, Quadratic = 2
	115
	116	Double_t intEr, intEphi ;
	117	Double_t r, phi, z ;
	118	Int_t sign ;
	119
	120	Double_t deltaVGG;
	121
	122	r = TMath::Sqrt( x[0]x[0] + x[1]x[1] );
	123	phi = TMath::ATan2(x[1],x[0]);
	124	if ( phi < 0 ) phi += TMath::TwoPi(); // Table uses phi from 0 to 2*Pi
	125	z = x[2] ;
	126
	127	if ( (roc%36) < 18 ) {
	128	sign = 1;
	129	deltaVGG = fDeltaVGGA; // (TPC End A)
	130	} else {
	131	sign = -1; // (TPC End C)
	132	deltaVGG = fDeltaVGGC;
	133	}
	134
	135	if ( sign==1 && z < fgkZOffSet ) z = fgkZOffSet; // Protect against discontinuity at CE
	136	if ( sign==-1 && z > -fgkZOffSet ) z = -fgkZOffSet; // Protect against discontinuity at CE
	137
	138	Interpolate2DEdistortion( order, r, z, fGGVoltErrorER, intEr );
	139	intEphi = 0.0; // Efield is symmetric in phi
	140
	141	// Calculate distorted position
	142	if ( r > 0.0 ) {
	143	phi = phi + deltaVGG( fC0intEphi - fC1*intEr ) / r;
	144	r = r + deltaVGG( fC0intEr + fC1*intEphi );
	145	}
	146
	147	// Calculate correction in cartesian coordinates
	148	dx[0] = r * TMath::Cos(phi) - x[0];
	149	dx[1] = r * TMath::Sin(phi) - x[1];
c9cbd2f2	150	dx[2] = 0.; // z distortion not implemented (1st order distortions) - see e.g. AliTPCBoundaryVoltError-class
0116859c	151
1b923461	152
1b923461	153
0116859c	154	}
	155
	156
	157	Float_t AliTPCGGVoltError::GetIntErOverEz(const Float_t x[],const Short_t roc) {
	158	//
	159	// This function is purely for calibration purposes
	160	// Calculates the integral (int Er/Ez dz) for the setted GG voltage offset
	161	//
1b923461	162
	163	if (!fInitLookUp) AliError("Lookup table was not initialized! You should do InitGGVoltErrorDistortion() ...");
	164
0116859c	165	Int_t order = 1 ; // FIXME: so far hardcoded? Linear interpolation = 1, Quadratic = 2
	166
	167	Double_t intEr;
	168	Double_t r, phi, z ;
	169	Int_t sign ;
	170
	171	Double_t deltaVGG;
	172
	173	r = TMath::Sqrt( x[0]x[0] + x[1]x[1] );
	174	phi = TMath::ATan2(x[1],x[0]);
	175	if ( phi < 0 ) phi += TMath::TwoPi(); // Table uses phi from 0 to 2*Pi
	176	z = x[2] ;
	177
	178	if ( (roc%36) < 18 ) {
	179	sign = 1;
	180	deltaVGG = fDeltaVGGA; // (TPC End A)
	181	} else {
	182	sign = -1; // (TPC End C)
	183	deltaVGG = fDeltaVGGC;
	184	}
	185
	186	if ( sign==1 && z < fgkZOffSet ) z = fgkZOffSet; // Protect against discontinuity at CE
	187	if ( sign==-1 && z > -fgkZOffSet ) z = -fgkZOffSet; // Protect against discontinuity at CE
	188
	189	Interpolate2DEdistortion(order, r, z, fGGVoltErrorER, intEr );
	190
	191	return (intEr*deltaVGG);
	192
	193	}
	194
	195	void AliTPCGGVoltError::InitGGVoltErrorDistortion() {
	196	//
	197	// Initialization of the Lookup table which contains the solutions of the GG Error problem
	198	//
	199
	200	Double_t r,z;
	201	Int_t nterms = 100 ;
	202	for ( Int_t i = 0 ; i < kNZ ; ++i ) {
	203	z = fgkZList[i] ;
	204	for ( Int_t j = 0 ; j < kNR ; ++j ) {
	205	r = fgkRList[j] ;
	206	fGGVoltErrorER[i][j] = 0.0 ;
	207	Double_t intz = 0.0 ;
	208	for ( Int_t n = 1 ; n < nterms ; ++n ) {
b1f0a2a5	209	Double_t k = n * TMath::Pi() / fgkTPCZ0 ;
0116859c	210	Double_t ein = 0 ; // Error potential on the IFC
	211	Double_t eout = 0 ; // Error potential on the OFC
	212	if ( z < 0 ) {
	213	ein = -2.0 / ( k * (fgkCathodeV - fgkGG) ) ;
	214	eout = -2.0 / ( k * (fgkCathodeV - fgkGG) ) ;
	215	}
	216	if ( z == 0 ) continue ;
	217	if ( z > 0 ) {
	218	ein = -2.0 / ( k * (fgkCathodeV - fgkGG) ) ;
	219	eout = -2.0 / ( k * (fgkCathodeV - fgkGG) ) ;
	220	}
	221	Double_t an = ein * TMath::BesselK0( kfgkOFCRadius ) - eout TMath::BesselK0( k*fgkIFCRadius ) ;
	222	Double_t bn = eout * TMath::BesselI0( kfgkIFCRadius ) - ein TMath::BesselI0( k*fgkOFCRadius ) ;
	223	Double_t numerator =
	224	an * TMath::BesselI1( kr ) - bn TMath::BesselK1( k*r ) ;
	225	Double_t denominator =
	226	TMath::BesselK0( kfgkOFCRadius ) TMath::BesselI0( k*fgkIFCRadius ) -
	227	TMath::BesselK0( kfgkIFCRadius ) TMath::BesselI0( k*fgkOFCRadius ) ;
b1f0a2a5	228	Double_t zterm = TMath::Cos( k*(fgkTPCZ0-TMath::Abs(z)) ) - 1 ;
0116859c	229	intz += zterm * numerator / denominator ;
0116859c	230	// Assume series converges, break if small terms
b9f518ba	231	if ( n>10 && TMath::Abs(intz)*1.e-10 > TMath::Abs(numerator/denominator) ) break;
0116859c	232	}
	233	fGGVoltErrorER[i][j] = (Double_t) intz ;
	234
	235	}
	236	}
1b923461	237
1b923461	238	fInitLookUp = kTRUE;
0116859c	239	}
	240
	241
	242
b1f0a2a5	243	void AliTPCGGVoltError::Print(const Option_t* option) const {
0116859c	244	//
	245	// Print function to check the settings (e.g. voltage offsets)
	246	// option=="a" prints the C0 and C1 coefficents for calibration purposes
	247	//
	248
	249	TString opt = option; opt.ToLower();
	250	printf("%s\n",GetTitle());
	251	printf(" - GG Voltage offset: A-side: %3.1f V, C-side: %3.1f V \n",fDeltaVGGA,fDeltaVGGC);
	252	if (opt.Contains("a")) { // Print all details
534fd34a	253	printf(" - T1: %1.4f, T2: %1.4f \n",fT1,fT2);
0116859c	254	printf(" - C1: %1.4f, C0: %1.4f \n",fC1,fC0);
	255	}
	256
1b923461	257	if (!fInitLookUp) AliError("Lookup table was not initialized! You should do InitGGVoltErrorDistortion() ...");
0116859c	258
	259
	260	}