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

// _________________________________________________________________
//
// Begin_Html
//   <h2> AliTPCSpaceCharge class  </h2>    
//   The class calculates the space point distortions due to a rotational 
//   symmetric space charge distribution with the TPC drift volume. 
//   <p>
//   The class uses the PoissonRelaxation2D to calculate the resulting 
//   electrical field inhomogeneities in the (r,z)-plane. Then, the 
//   Langevin-integral formalism is used to calculate the space point distortions. 
//   <p>
//   The class assumes, that the distortions scales linearly with the magnitude 
//   of the space charge distribution $\rho(r,z)$. The in here assumed distribution is 
//   $$\rho(r,z) = \frac{(A-B\,z)}{r^2} $$ wherein the factors A and B scale with the
//   event multiplicity and the interaction rate.
//   <p>
//   The scaling factor can be set via the function SetCorrectionFactor. An example of 
//   the shape of the distortions is given below.
//
//   MI modification - 22.05.2013
//   As an optional input the Space charge histogram RZ is used in case it is provided 
//     - using the SetInputSpaceCharge function
//
// End_Html
//
// Begin_Macro(source)
//   {
//   gROOT->SetStyle("Plain"); gStyle->SetPalette(1);
//   TCanvas *c2 = new TCanvas("cAliTPCSpaceCharge","cAliTPCSpaceCharge",500,300); 
//   AliTPCSpaceCharge sc;
//   sc.SetOmegaTauT1T2(-0.32,1,1); // B=0.5 Tesla
//   sc.SetCorrectionFactor(0.0015);
//   sc.CreateHistoDRinZR(0.)->Draw("surf2");
//   return c2;
//   } 
// End_Macro
//
// Begin_Html
//   <p>
//   Date: 23/08/2010 <br>
// Authors: Jim Thomas, Stefan Rossegger 
// End_Html 
// _________________________________________________________________


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

#include "TMath.h"
#include "AliTPCROC.h"
#include "AliTPCSpaceCharge.h"

ClassImp(AliTPCSpaceCharge)

AliTPCSpaceCharge::AliTPCSpaceCharge()
  : AliTPCCorrection("SpaceCharge2D","Space Charge 2D"),
    fC0(0.),fC1(0.),fCorrectionFactor(0.001),fSpaceChargeHistogram(0),
    fInitLookUp(kFALSE)
{
  //
  // default constructor
  //
 
}

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


void AliTPCSpaceCharge::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);

  InitSpaceChargeDistortion(); // fill the look up table
}

void AliTPCSpaceCharge::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);

  //  SetCorrectionFactor(1.); // should come from some database

}


void AliTPCSpaceCharge::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {
  //
  // Calculates the correction due the Space Charge effect within the TPC drift volume
  //   

  if (!fInitLookUp) {
    AliInfo("Lookup table was not initialized! Perform the inizialisation now ...");
    InitSpaceChargeDistortion();
  }
  Int_t   order     = 1 ;    // FIXME: hardcoded? Linear interpolation = 1, Quadratic = 2         
                        
  Double_t intEr, intEphi, intdEz;
  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!");

  // Efield is symmetric in phi - 2D calculation
  intEphi = 0.0; 
  // Get the E field integrals
  Interpolate2DEdistortion( order, r, z, fLookUpErOverEz, intEr );
  // Get DeltaEz field integral
  Interpolate2DEdistortion( order, r, z, fLookUpDeltaEz, intdEz );
  
 
  // Calculate distorted position
  if ( r > 0.0 ) {
    phi =  phi + fCorrectionFactor *( fC0*intEphi - fC1*intEr ) / r;      
    r   =  r   + fCorrectionFactor *( fC0*intEr   + fC1*intEphi );  
  }
  Double_t dz = intdEz*fCorrectionFactor;
 
  // Calculate correction in cartesian coordinates
  dx[0] = - (r * TMath::Cos(phi) - x[0]);
  dx[1] = - (r * TMath::Sin(phi) - x[1]); 
  dx[2] = - dz;  // z distortion - (internally scaled with driftvelocity dependency 
                 // on the Ez field 

}

void AliTPCSpaceCharge::InitSpaceChargeDistortion() {
  //
  // Initialization of the Lookup table which contains the solutions of the 
  // poisson problem
  //

  const Float_t  gridSizeR   =  (fgkOFCRadius-fgkIFCRadius) / (kRows-1) ;
  const Float_t  gridSizeZ   =  fgkTPCZ0 / (kColumns-1) ;

  TMatrixD voltArray(kRows,kColumns);        // dummy boundary vectors
  TMatrixD chargeDensity(kRows,kColumns);    // charge
  TMatrixD arrayErOverEz(kRows,kColumns);    // solution in Er
  TMatrixD arrayDeltaEz(kRows,kColumns);    // solution in Ez

  Double_t  rList[kRows], zedList[kColumns] ;
  
  // Fill arrays with initial conditions.  V on the boundary and ChargeDensity in the volume.      
  for ( Int_t j = 0 ; j < kColumns ; j++ ) {
    Double_t zed = j*gridSizeZ ;
    zedList[j] = zed ;
    for ( Int_t i = 0 ; i < kRows ; i++ )  {
      Double_t radius = fgkIFCRadius + i*gridSizeR ;
      rList[i]           = radius ;
      voltArray(i,j)        = 0;  // Initialize voltArray to zero - not used in this class
      chargeDensity(i,j)     = 0;  // Initialize ChargeDensity to zero
    }
  }      

  // Fill the initial conditions
  for ( Int_t j = 1 ; j < kColumns-1 ; j++ ) {
    Double_t zed = j*gridSizeZ ;
    for ( Int_t i = 1 ; i < kRows-1 ; i++ ) { 
      Double_t radius = fgkIFCRadius + i*gridSizeR ;

      Double_t zterm = (fgkTPCZ0-zed) * (fgkOFCRadius*fgkOFCRadius - fgkIFCRadius*fgkIFCRadius) / fgkTPCZ0 ;
      // for 1/R**2 charge density in the TPC; then integrated in Z due to drifting ions
      chargeDensity(i,j) = zterm / ( TMath::Log(fgkOFCRadius/fgkIFCRadius) * ( radius*radius ) ) ;              
    }
  }
  // Fill the initial space charge in case histogram exist
  if (fSpaceChargeHistogram){
    for ( Int_t j = 1 ; j < kColumns-1 ; j++ ) {
      Double_t zed = j*gridSizeZ ;
      for ( Int_t i = 1 ; i < kRows-1 ; i++ ) { 
	Double_t radius = fgkIFCRadius + i*gridSizeR ;
	
	Double_t zterm = (fgkTPCZ0-zed) * (fgkOFCRadius*fgkOFCRadius - fgkIFCRadius*fgkIFCRadius) / fgkTPCZ0 ;
	// for 1/R**2 charge density in the TPC; then integrated in Z due to drifting ions
	chargeDensity(i,j) = fSpaceChargeHistogram->Interpolate(radius,zed);
      }
    }
  }


  // Solve the electrosatic problem in 2D 

  PoissonRelaxation2D( voltArray, chargeDensity, arrayErOverEz, arrayDeltaEz, kRows, kColumns, kIterations ) ;
  
  //Interpolate results onto standard grid for Electric Fields
  Int_t ilow=0, jlow=0 ;
  Double_t z,r;
  Float_t saveEr[2], saveEz[2] ;	      
  for ( Int_t i = 0 ; i < kNZ ; ++i )  {
    z = TMath::Abs( fgkZList[i] ) ; // assume symmetric behaviour on A and C side
    for ( Int_t j = 0 ; j < kNR ; ++j ) {

      // Linear interpolation !!
      r = fgkRList[j] ;
      Search( kRows,   rList, r, ilow ) ;          // Note switch - R in rows and Z in columns
      Search( kColumns, zedList, z, jlow ) ;
      if ( ilow < 0 ) ilow = 0 ;                   // check if out of range
      if ( jlow < 0 ) jlow = 0 ;   
      if ( ilow + 1  >=  kRows - 1 ) ilow =  kRows - 2 ;	      
      if ( jlow + 1  >=  kColumns - 1 ) jlow =  kColumns - 2 ; 
    
      saveEr[0] = arrayErOverEz(ilow,jlow) + 
	(arrayErOverEz(ilow,jlow+1)-arrayErOverEz(ilow,jlow))*(z-zedList[jlow])/gridSizeZ ;
      saveEr[1] = arrayErOverEz(ilow+1,jlow) + 
	(arrayErOverEz(ilow+1,jlow+1)-arrayErOverEz(ilow+1,jlow))*(z-zedList[jlow])/gridSizeZ ;
      saveEz[0] = arrayDeltaEz(ilow,jlow) + 
	(arrayDeltaEz(ilow,jlow+1)-arrayDeltaEz(ilow,jlow))*(z-zedList[jlow])/gridSizeZ ;
      saveEz[1] = arrayDeltaEz(ilow+1,jlow) + 
	(arrayDeltaEz(ilow+1,jlow+1)-arrayDeltaEz(ilow+1,jlow))*(z-zedList[jlow])/gridSizeZ ;

      
      fLookUpErOverEz[i][j] = saveEr[0] + (saveEr[1]-saveEr[0])*(r-rList[ilow])/gridSizeR ;
      fLookUpDeltaEz[i][j]  = saveEz[0] + (saveEz[1]-saveEz[0])*(r-rList[ilow])/gridSizeR ;

      if (fgkZList[i]<0)  fLookUpDeltaEz[i][j] *= -1; // C side is negative z
    }
  }
  
  fInitLookUp = kTRUE;

}

void AliTPCSpaceCharge::Print(const Option_t* option) const {
  //
  // Print function to check the settings of the boundary vectors
  // option=="a" prints the C0 and C1 coefficents for calibration purposes
  //

  TString opt = option; opt.ToLower();
  printf("%s\n",GetTitle());
  printf(" - Space Charge effects assuming a radial symmetric z over r^2 SC-distribution.\n");
  printf("   SC correction factor: %f \n",fCorrectionFactor);

  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 InitSpaceChargeDistortion() ...");

}
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
b4caed64	16	// _________________________________________________________________
	17	//
	18	// Begin_Html
	19	// <h2> AliTPCSpaceCharge class </h2>
	20	// The class calculates the space point distortions due to a rotational
	21	// symmetric space charge distribution with the TPC drift volume.
	22	// <p>
	23	// The class uses the PoissonRelaxation2D to calculate the resulting
	24	// electrical field inhomogeneities in the (r,z)-plane. Then, the
	25	// Langevin-integral formalism is used to calculate the space point distortions.
	26	// <p>
	27	// The class assumes, that the distortions scales linearly with the magnitude
	28	// of the space charge distribution $\rho(r,z)$. The in here assumed distribution is
	29	// $$\rho(r,z) = \frac{(A-B\,z)}{r^2} $$ wherein the factors A and B scale with the
	30	// event multiplicity and the interaction rate.
	31	// <p>
	32	// The scaling factor can be set via the function SetCorrectionFactor. An example of
	33	// the shape of the distortions is given below.
276c23d3	34	//
	35	// MI modification - 22.05.2013
	36	// As an optional input the Space charge histogram RZ is used in case it is provided
	37	// - using the SetInputSpaceCharge function
	38	//
b4caed64	39	// End_Html
b4caed64	40	//
6a1caa6b	41	// Begin_Macro(source)
b4caed64	42	// {
b4caed64	43	// gROOT->SetStyle("Plain"); gStyle->SetPalette(1);
6a1caa6b	44	// TCanvas *c2 = new TCanvas("cAliTPCSpaceCharge","cAliTPCSpaceCharge",500,300);
b4caed64	45	// AliTPCSpaceCharge sc;
	46	// sc.SetOmegaTauT1T2(-0.32,1,1); // B=0.5 Tesla
	47	// sc.SetCorrectionFactor(0.0015);
	48	// sc.CreateHistoDRinZR(0.)->Draw("surf2");
	49	// return c2;
	50	// }
	51	// End_Macro
	52	//
	53	// Begin_Html
	54	// <p>
	55	// Date: 23/08/2010 <br>
	56	// Authors: Jim Thomas, Stefan Rossegger
	57	// End_Html
	58	// _________________________________________________________________
	59
	60
c9cbd2f2	61
	62	#include "AliMagF.h"
	63	#include "TGeoGlobalMagField.h"
	64	#include "AliTPCcalibDB.h"
	65	#include "AliTPCParam.h"
	66	#include "AliLog.h"
	67	#include "TMatrixD.h"
276c23d3	68	#include "TH2.h"
c9cbd2f2	69
	70	#include "TMath.h"
	71	#include "AliTPCROC.h"
	72	#include "AliTPCSpaceCharge.h"
	73
	74	ClassImp(AliTPCSpaceCharge)
	75
	76	AliTPCSpaceCharge::AliTPCSpaceCharge()
	77	: AliTPCCorrection("SpaceCharge2D","Space Charge 2D"),
276c23d3	78	fC0(0.),fC1(0.),fCorrectionFactor(0.001),fSpaceChargeHistogram(0),
c9cbd2f2	79	fInitLookUp(kFALSE)
	80	{
	81	//
	82	// default constructor
	83	//
	84
	85	}
	86
	87	AliTPCSpaceCharge::~AliTPCSpaceCharge() {
	88	//
	89	// default destructor
	90	//
	91	}
	92
	93
	94
	95	void AliTPCSpaceCharge::Init() {
	96	//
	97	// Initialization funtion
	98	//
	99
	100	AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
	101	if (!magF) AliError("Magneticd field - not initialized");
	102	Double_t bzField = magF->SolenoidField()/10.; //field in T
	103	AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
	104	if (!param) AliError("Parameters - not initialized");
	105	Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us] // From dataBase: to be updated: per second (ideally)
	106	Double_t ezField = 400; // [V/cm] // to be updated: never (hopefully)
	107	Double_t wt = -10.0 * (bzField10) vdrift / ezField ;
	108	// Correction Terms for effective omegaTau; obtained by a laser calibration run
	109	SetOmegaTauT1T2(wt,fT1,fT2);
	110
	111	InitSpaceChargeDistortion(); // fill the look up table
	112	}
	113
	114	void AliTPCSpaceCharge::Update(const TTimeStamp &/timeStamp/) {
	115	//
	116	// Update function
	117	//
	118	AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
	119	if (!magF) AliError("Magneticd field - not initialized");
	120	Double_t bzField = magF->SolenoidField()/10.; //field in T
	121	AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
	122	if (!param) AliError("Parameters - not initialized");
	123	Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us] // From dataBase: to be updated: per second (ideally)
	124	Double_t ezField = 400; // [V/cm] // to be updated: never (hopefully)
	125	Double_t wt = -10.0 * (bzField10) vdrift / ezField ;
	126	// Correction Terms for effective omegaTau; obtained by a laser calibration run
	127	SetOmegaTauT1T2(wt,fT1,fT2);
	128
	129	// SetCorrectionFactor(1.); // should come from some database
	130
	131	}
	132
	133
	134
	135	void AliTPCSpaceCharge::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {
	136	//
	137	// Calculates the correction due the Space Charge effect within the TPC drift volume
	138	//
	139
	140	if (!fInitLookUp) {
	141	AliInfo("Lookup table was not initialized! Perform the inizialisation now ...");
	142	InitSpaceChargeDistortion();
143	}
144	Int_t order = 1 ; // FIXME: hardcoded? Linear interpolation = 1, Quadratic = 2
145
146	Double_t intEr, intEphi, intdEz;
147	Double_t r, phi, z ;
148	Int_t sign;
149
150	r = TMath::Sqrt( x[0]x[0] + x[1]x[1] ) ;
151	phi = TMath::ATan2(x[1],x[0]) ;
152	if ( phi < 0 ) phi += TMath::TwoPi() ; // Table uses phi from 0 to 2*Pi
153	z = x[2] ; // Create temporary copy of x[2]
154
155	if ( (roc%36) < 18 ) {
156	sign = 1; // (TPC A side)
157	} else {
158	sign = -1; // (TPC C side)
159	}
160
161	if ( sign==1 && z < fgkZOffSet ) z = fgkZOffSet; // Protect against discontinuity at CE
162	if ( sign==-1 && z > -fgkZOffSet ) z = -fgkZOffSet; // Protect against discontinuity at CE
163
164
165	if ( (sign==1 && z<0) \|\| (sign==-1 && z>0) ) // just a consistency check
166	AliError("ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!");
167
168	// Efield is symmetric in phi - 2D calculation
169	intEphi = 0.0;
170	// Get the E field integrals
171	Interpolate2DEdistortion( order, r, z, fLookUpErOverEz, intEr );
172	// Get DeltaEz field integral
173	Interpolate2DEdistortion( order, r, z, fLookUpDeltaEz, intdEz );
174
175
176	// Calculate distorted position
177	if ( r > 0.0 ) {
178	phi = phi + fCorrectionFactor ( fC0intEphi - fC1*intEr ) / r;
179	r = r + fCorrectionFactor ( fC0intEr + fC1*intEphi );
180	}
181	Double_t dz = intdEz*fCorrectionFactor;
182
183	// Calculate correction in cartesian coordinates
752b0cc7	184	dx[0] = - (r * TMath::Cos(phi) - x[0]);
	185	dx[1] = - (r * TMath::Sin(phi) - x[1]);
	186	dx[2] = - dz; // z distortion - (internally scaled with driftvelocity dependency
	187	// on the Ez field
c9cbd2f2	188
	189	}
	190
	191	void AliTPCSpaceCharge::InitSpaceChargeDistortion() {
	192	//
	193	// Initialization of the Lookup table which contains the solutions of the
	194	// poisson problem
	195	//
	196
	197	const Float_t gridSizeR = (fgkOFCRadius-fgkIFCRadius) / (kRows-1) ;
	198	const Float_t gridSizeZ = fgkTPCZ0 / (kColumns-1) ;
	199
	200	TMatrixD voltArray(kRows,kColumns); // dummy boundary vectors
	201	TMatrixD chargeDensity(kRows,kColumns); // charge
	202	TMatrixD arrayErOverEz(kRows,kColumns); // solution in Er
	203	TMatrixD arrayDeltaEz(kRows,kColumns); // solution in Ez
	204
	205	Double_t rList[kRows], zedList[kColumns] ;
	206
	207	// Fill arrays with initial conditions. V on the boundary and ChargeDensity in the volume.
	208	for ( Int_t j = 0 ; j < kColumns ; j++ ) {
	209	Double_t zed = j*gridSizeZ ;
	210	zedList[j] = zed ;
	211	for ( Int_t i = 0 ; i < kRows ; i++ ) {
	212	Double_t radius = fgkIFCRadius + i*gridSizeR ;
	213	rList[i] = radius ;
	214	voltArray(i,j) = 0; // Initialize voltArray to zero - not used in this class
	215	chargeDensity(i,j) = 0; // Initialize ChargeDensity to zero
	216	}
	217	}
	218
	219	// Fill the initial conditions
	220	for ( Int_t j = 1 ; j < kColumns-1 ; j++ ) {
	221	Double_t zed = j*gridSizeZ ;
	222	for ( Int_t i = 1 ; i < kRows-1 ; i++ ) {
	223	Double_t radius = fgkIFCRadius + i*gridSizeR ;
	224
	225	Double_t zterm = (fgkTPCZ0-zed) * (fgkOFCRadiusfgkOFCRadius - fgkIFCRadiusfgkIFCRadius) / fgkTPCZ0 ;
	226	// for 1/R**2 charge density in the TPC; then integrated in Z due to drifting ions
	227	chargeDensity(i,j) = zterm / ( TMath::Log(fgkOFCRadius/fgkIFCRadius) * ( radius*radius ) ) ;
	228	}
	229	}
276c23d3	230	// Fill the initial space charge in case histogram exist
	231	if (fSpaceChargeHistogram){
	232	for ( Int_t j = 1 ; j < kColumns-1 ; j++ ) {
	233	Double_t zed = j*gridSizeZ ;
	234	for ( Int_t i = 1 ; i < kRows-1 ; i++ ) {
	235	Double_t radius = fgkIFCRadius + i*gridSizeR ;
	236
	237	Double_t zterm = (fgkTPCZ0-zed) * (fgkOFCRadiusfgkOFCRadius - fgkIFCRadiusfgkIFCRadius) / fgkTPCZ0 ;
	238	// for 1/R**2 charge density in the TPC; then integrated in Z due to drifting ions
	239	chargeDensity(i,j) = fSpaceChargeHistogram->Interpolate(radius,zed);
	240	}
	241	}
	242	}
c9cbd2f2	243
	244
	245	// Solve the electrosatic problem in 2D
	246
	247	PoissonRelaxation2D( voltArray, chargeDensity, arrayErOverEz, arrayDeltaEz, kRows, kColumns, kIterations ) ;
	248
	249	//Interpolate results onto standard grid for Electric Fields
	250	Int_t ilow=0, jlow=0 ;
	251	Double_t z,r;
	252	Float_t saveEr[2], saveEz[2] ;
	253	for ( Int_t i = 0 ; i < kNZ ; ++i ) {
	254	z = TMath::Abs( fgkZList[i] ) ; // assume symmetric behaviour on A and C side
	255	for ( Int_t j = 0 ; j < kNR ; ++j ) {
	256
	257	// Linear interpolation !!
	258	r = fgkRList[j] ;
	259	Search( kRows, rList, r, ilow ) ; // Note switch - R in rows and Z in columns
	260	Search( kColumns, zedList, z, jlow ) ;
	261	if ( ilow < 0 ) ilow = 0 ; // check if out of range
	262	if ( jlow < 0 ) jlow = 0 ;
	263	if ( ilow + 1 >= kRows - 1 ) ilow = kRows - 2 ;
	264	if ( jlow + 1 >= kColumns - 1 ) jlow = kColumns - 2 ;
	265
	266	saveEr[0] = arrayErOverEz(ilow,jlow) +
	267	(arrayErOverEz(ilow,jlow+1)-arrayErOverEz(ilow,jlow))*(z-zedList[jlow])/gridSizeZ ;
	268	saveEr[1] = arrayErOverEz(ilow+1,jlow) +
	269	(arrayErOverEz(ilow+1,jlow+1)-arrayErOverEz(ilow+1,jlow))*(z-zedList[jlow])/gridSizeZ ;
	270	saveEz[0] = arrayDeltaEz(ilow,jlow) +
	271	(arrayDeltaEz(ilow,jlow+1)-arrayDeltaEz(ilow,jlow))*(z-zedList[jlow])/gridSizeZ ;
	272	saveEz[1] = arrayDeltaEz(ilow+1,jlow) +
	273	(arrayDeltaEz(ilow+1,jlow+1)-arrayDeltaEz(ilow+1,jlow))*(z-zedList[jlow])/gridSizeZ ;
	274
	275
	276	fLookUpErOverEz[i][j] = saveEr[0] + (saveEr[1]-saveEr[0])*(r-rList[ilow])/gridSizeR ;
	277	fLookUpDeltaEz[i][j] = saveEz[0] + (saveEz[1]-saveEz[0])*(r-rList[ilow])/gridSizeR ;
15687d71	278
15687d71	279	if (fgkZList[i]<0) fLookUpDeltaEz[i][j] *= -1; // C side is negative z
c9cbd2f2	280	}
	281	}
	282
	283	fInitLookUp = kTRUE;
	284
	285	}
	286
	287	void AliTPCSpaceCharge::Print(const Option_t* option) const {
	288	//
	289	// Print function to check the settings of the boundary vectors
	290	// option=="a" prints the C0 and C1 coefficents for calibration purposes
	291	//
	292
	293	TString opt = option; opt.ToLower();
	294	printf("%s\n",GetTitle());
	295	printf(" - Space Charge effects assuming a radial symmetric z over r^2 SC-distribution.\n");
	296	printf(" SC correction factor: %f \n",fCorrectionFactor);
	297
	298	if (opt.Contains("a")) { // Print all details
	299	printf(" - T1: %1.4f, T2: %1.4f \n",fT1,fT2);
	300	printf(" - C1: %1.4f, C0: %1.4f \n",fC1,fC0);
	301	}
	302
	303	if (!fInitLookUp) AliError("Lookup table was not initialized! You should do InitSpaceChargeDistortion() ...");
	304
	305	}