First version of a class to correct for 3D Space Charges
authorsrossegg <srossegg@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 2 Sep 2010 17:43:47 +0000 (17:43 +0000)
committersrossegg <srossegg@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 2 Sep 2010 17:43:47 +0000 (17:43 +0000)
TPC/AliTPCSpaceCharge3D.cxx [new file with mode: 0644]
TPC/AliTPCSpaceCharge3D.h [new file with mode: 0644]

diff --git a/TPC/AliTPCSpaceCharge3D.cxx b/TPC/AliTPCSpaceCharge3D.cxx
new file mode 100644 (file)
index 0000000..56ff58b
--- /dev/null
@@ -0,0 +1,822 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+//////////////////////////////////////////////////////////////////////////////
+//                                                                          //
+// AliTPCSpaceCharge3D class                                                //
+// The class calculates the space point distortions due to a space charge   //
+// effect ....                                                              //
+// Method of calculation:                                                   //
+// The analytical solution for the poisson problem in 3D (cylindrical coord)//
+// is used in form of look up tables. PieceOfCake (POC) voxel were pre-     //
+// calculated and can be sumed up (weighted) according to the what is needed// 
+//                                                                          //
+// The class allows "effective Omega Tau" corrections.                      // 
+//                                                                          //
+// NOTE: This class is not  capable of calculation z distortions due to     //
+//       drift velocity change in dependence of the electric field!!!       //
+//                                                                          //
+// date: 19/06/2010                                                         //
+// Authors: Stefan Rossegger                                                //
+//                                                                          //
+// Example usage:                                                           //
+//////////////////////////////////////////////////////////////////////////////
+
+#include "AliMagF.h"
+#include "TGeoGlobalMagField.h"
+#include "AliTPCcalibDB.h"
+#include "AliTPCParam.h"
+#include "AliLog.h"
+#include "TH2F.h"
+#include "TH3F.h"
+#include "TFile.h"
+#include "TVector.h"
+#include "TMatrix.h"
+#include "TMatrixD.h"
+
+#include "TMath.h"
+#include "AliTPCROC.h"
+#include "AliTPCSpaceCharge3D.h"
+
+ClassImp(AliTPCSpaceCharge3D)
+
+AliTPCSpaceCharge3D::AliTPCSpaceCharge3D()
+  : AliTPCCorrection("SpaceCharge3D","Space Charge - 3D"),
+    fC0(0.),fC1(0.),
+    fCorrectionFactor(1.),
+    fInitLookUp(kFALSE),
+    fSCDataFileName((char*)"$(ALICE_ROOT)/TPC/Calib/maps/sc_3D_distribution_Sim.root"),
+    fSCLookUpPOCsFileName((char*)"$(ALICE_ROOT)/TPC/Calib/maps/sc_3D_raw_18-18-26_17p-18p-25p-MN30.root")
+{
+  //
+  // default constructor
+  //
+
+  // Array which will contain the solution according to the setted charge density distribution
+  // see InitSpaceCharge3DDistortion() 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); 
+    fSCdensityDistribution[k] = new TMatrixD(kNR,kNZ);
+  }
+
+  printf("%s\n",fSCDataFileName);
+  printf("%s\n",fSCLookUpPOCsFileName);
+  SetSCDataFileName(fSCDataFileName);
+
+}
+
+AliTPCSpaceCharge3D::~AliTPCSpaceCharge3D() {
+  //
+  // default destructor
+  //
+  for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
+    delete fLookUpErOverEz[k];
+    delete fLookUpEphiOverEz[k];
+    delete fLookUpDeltaEz[k];
+    delete fSCdensityDistribution[k];
+  }
+
+  //  delete fSCdensityDistribution;
+}
+
+
+
+void AliTPCSpaceCharge3D::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);
+
+  InitSpaceCharge3DDistortion(); // fill the look up table
+}
+
+void AliTPCSpaceCharge3D::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 AliTPCSpaceCharge3D::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! Performing the inizialisation now ...");
+    InitSpaceCharge3DDistortion();
+  }
+
+  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!");
+
+  // 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);
+  intdEz = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi, 
+                                 fgkRList, fgkZList, fgkPhiList, fLookUpDeltaEz   );
+
+  // 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 * fgkdvdE;
+  // 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 - (scaled with driftvelocity dependency on the Ez field and the overall scaling factor)
+
+}
+
+
+void AliTPCSpaceCharge3D::InitSpaceCharge3DDistortion() {
+  //
+  // Initialization of the Lookup table which contains the solutions of the 
+  // "space charge" (poisson) problem 
+  //
+  // The sum-up uses a look-up table which contains different discretized Space charge fields 
+  // in order to calculate the corresponding field deviations due to a given (discretized)
+  // space charge distribution ....
+  //
+  // Method of calculation: Weighted sum up of the different fields within the look up table
+  //
+
+  if (fInitLookUp) {
+    AliInfo("Lookup table was already initialized!");
+    //    return;
+  }
+
+  AliInfo("Preparation of the weighted look-up table");
+   
+  TFile *f = new TFile(fSCLookUpPOCsFileName);
+  if (!f) { 
+    AliError("Precalculated POC-looup-table could not be found");
+    return;
+  }
+
+  // units are in [m]
+  TVector *gridf  = (TVector*) f->Get("constants"); 
+  TVector &grid = *gridf;
+  TMatrix *coordf  = (TMatrix*) f->Get("coordinates");
+  TMatrix &coord = *coordf;
+  TMatrix *coordPOCf  = (TMatrix*) f->Get("POCcoord");
+  TMatrix &coordPOC = *coordPOCf;
+  
+  Bool_t flagRadSym = 0;
+  if (grid(1)==1 && grid(4)==1) {
+    AliInfo("LOOK UP TABLE IS RADIAL SYMETTRIC - Field in Phi is ZERO");
+    flagRadSym=1;
+  }
+
+  Int_t rows      = (Int_t)grid(0);   // number of points in r direction 
+  Int_t phiSlices = (Int_t)grid(1);   // number of points in phi         
+  Int_t columns   = (Int_t)grid(2);   // number of points in z direction 
+
+  const Float_t gridSizeR   =  grid(6)*100;    // unit in [cm]
+  const Float_t gridSizePhi   =  grid(7);  // unit in [rad]
+  const Float_t gridSizeZ =  grid(8)*100;      // unit in [cm]
+  // temporary matrices needed for the calculation
+  TMatrixD *arrayofErA[phiSlices], *arrayofEphiA[phiSlices], *arrayofdEzA[phiSlices];
+  TMatrixD *arrayofErC[phiSlices], *arrayofEphiC[phiSlices], *arrayofdEzC[phiSlices];
+
+  TMatrixD *arrayofEroverEzA[phiSlices], *arrayofEphioverEzA[phiSlices], *arrayofDeltaEzA[phiSlices];
+  TMatrixD *arrayofEroverEzC[phiSlices], *arrayofEphioverEzC[phiSlices], *arrayofDeltaEzC[phiSlices];
+
+  for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
+   
+    arrayofErA[k]   =   new TMatrixD(rows,columns) ;
+    arrayofEphiA[k] =   new TMatrixD(rows,columns) ; // zero if radial symmetric
+    arrayofdEzA[k]  =   new TMatrixD(rows,columns) ;
+    arrayofErC[k]   =   new TMatrixD(rows,columns) ;
+    arrayofEphiC[k] =   new TMatrixD(rows,columns) ; // zero if radial symmetric
+    arrayofdEzC[k]  =   new TMatrixD(rows,columns) ;
+
+    arrayofEroverEzA[k]   =   new TMatrixD(rows,columns) ;
+    arrayofEphioverEzA[k] =   new TMatrixD(rows,columns) ; // zero if radial symmetric
+    arrayofDeltaEzA[k]    =   new TMatrixD(rows,columns) ;
+    arrayofEroverEzC[k]   =   new TMatrixD(rows,columns) ;
+    arrayofEphioverEzC[k] =   new TMatrixD(rows,columns) ; // zero if radial symmetric
+    arrayofDeltaEzC[k]    =   new TMatrixD(rows,columns) ;
+
+    // Set the values to zero the lookup tables
+    // not necessary, it is done in the constructor of TMatrix - code deleted
+
+  }
+  // list of points as used in the interpolation (during sum up)
+  Double_t  rlist[rows], zedlist[columns] , philist[phiSlices];
+  for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
+    philist[k] =  gridSizePhi * k;
+    for ( Int_t i = 0 ; i < rows ; i++ )    {
+      rlist[i] = fgkIFCRadius + i*gridSizeR ;
+      for ( Int_t j = 0 ; j < columns ; j++ ) { 
+       zedlist[j]  = j * gridSizeZ ;
+      }
+    }
+  } // only done once
+  
+  
+  // Prepare summation Vectors
+  //  Int_t numPosInd = (Int_t)(grid(0)*grid(1)*grid(2)); // Number of fulcrums in the look-up table
+  //  Int_t numPOC = (Int_t)(grid(3)*grid(4)*grid(5));    // Number of POC conf. in the look-up table
+  
+  TTree *treePOC = (TTree*)f->Get("POCall");
+
+  TVector *bEr  = 0;   TVector *bEphi= 0;   TVector *bEz  = 0;
+  
+  treePOC->SetBranchAddress("Er",&bEr);
+  if (!flagRadSym) treePOC->SetBranchAddress("Ephi",&bEphi);
+  treePOC->SetBranchAddress("Ez",&bEz);
+
+  /*  TBranch *b2Er = treePOC->GetBranch("Er");
+      TBranch *b2Ephi =0;
+      if (!flagRadSym) b2Ephi = treePOC->GetBranch("Ephi");
+      TBranch *b2Ez = treePOC->GetBranch("Ez");
+  */
+
+  // Read the complete tree and do a weighted sum-up over the POC configurations
+  // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+  
+  Int_t treeNumPOC = (Int_t)treePOC->GetEntries(); // Number of POC conf. in the look-up table
+  Int_t ipC = 0; // POC Conf. counter (note: different to the POC number in the tree!)
+
+  for (Int_t itreepC=0; itreepC<treeNumPOC; itreepC++) { // loop over POC configurations in tree
+  
+    // if (!flagRadSym) cout<<ipC<<endl;
+
+    treePOC->GetEntry(itreepC);
+
+    /*    b2Er->GetEntry(itreepC); 
+         if (!flagRadSym) b2Ephi->GetEntry(itreepC); 
+         b2Ez->GetEntry(itreepC); 
+    */
+
+    // center of the POC voxel in [meter]
+    Double_t r0 = coordPOC(ipC,0);
+    Double_t phi0 = coordPOC(ipC,1);
+    Double_t z0 = coordPOC(ipC,2);
+
+    ipC++; // POC conf. counter
+
+    // weights (charge density) at POC position on the A and C side (in C/m^3/e0)
+    // note: coordinates are in [cm]
+    Double_t weightA = GetSpaceChargeDensity(r0*100,phi0, z0*100); 
+    Double_t weightC = GetSpaceChargeDensity(r0*100,phi0,-z0*100);
+  
+    // Summing up the vector components according to their weight
+
+    Int_t ip = 0;
+    for ( Int_t j = 0 ; j < columns ; j++ ) { 
+      for ( Int_t i = 0 ; i < rows ; i++ )    {
+       for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
+                 
+         // check wether the coordinates were screwed
+         if ((coord(0,ip)*100-rlist[i])>0.01 || 
+             (coord(1,ip)-philist[k])>0.01 || 
+             (coord(2,ip)*100-zedlist[j])>0.01) {
+           AliError("internal error: coordinate system was screwed during the sum-up");
+           printf("lookup: (r,phi,z)=(%lf,%lf,%lf)\n",coord(0,ip)*100,coord(1,ip),coord(2,ip)*100);
+           printf("sum-up: (r,phi,z)=(%lf,%lf,%lf)\n",rlist[i],philist[k],zedlist[j]);
+           AliError("Don't trust the results of the space charge calibration!");
+         }
+         
+         // unfortunately, the lookup tables were produced to be faster for phi symmetric charges
+         // This will be the most frequent usage (hopefully)
+         // That's why we have to do this here ...
+
+         TMatrixD &erA   =  *arrayofErA[k]  ;
+         TMatrixD &ephiA =  *arrayofEphiA[k];
+         TMatrixD &dEzA    =  *arrayofdEzA[k]   ;
+   
+         TMatrixD &erC   =  *arrayofErC[k]  ;
+         TMatrixD &ephiC =  *arrayofEphiC[k];
+         TMatrixD &dEzC    =  *arrayofdEzC[k]   ;
+   
+         // Sum up - Efield values in [V/m] -> transition to [V/cm]
+         erA(i,j) += ((*bEr)(ip)) * weightA /100;
+         erC(i,j) += ((*bEr)(ip)) * weightC /100;
+         if (!flagRadSym) {
+           ephiA(i,j) += ((*bEphi)(ip)) * weightA; // [V/rad]
+           ephiC(i,j) += ((*bEphi)(ip)) * weightC; // [V/rad]
+         }
+         dEzA(i,j)  += ((*bEz)(ip)) * weightA /100;
+         dEzC(i,j)  += ((*bEz)(ip)) * weightC /100;
+
+         // increase the counter
+         ip++;
+       }
+      }
+    }
+    
+    // Rotation and summation in the rest of the dPhiSteps
+    // which were not stored in the tree due to storage & symmetry reasons
+
+    Int_t phiPoints = (Int_t) grid(1);
+    Int_t phiPOC    = (Int_t) grid(4);
+    
+    // printf("%d %d\n",phiPOC,flagRadSym);
+    
+    for (Int_t phiiC = 1; phiiC<phiPOC; phiiC++) { // just used for non-radial symetric table
+      
+      r0 = coordPOC(ipC,0);
+      phi0 = coordPOC(ipC,1);
+      z0 = coordPOC(ipC,2);
+      
+      ipC++; // POC conf. counter
+      
+      // weights (charge density) at POC position on the A and C side (in C/m^3/e0)
+      // note: coordinates are in [cm]
+      weightA = GetSpaceChargeDensity(r0*100,phi0, z0*100); 
+      weightC = GetSpaceChargeDensity(r0*100,phi0,-z0*100);
+      
+      
+      // Summing up the vector components according to their weight
+      ip = 0;
+      for ( Int_t j = 0 ; j < columns ; j++ ) { 
+       for ( Int_t i = 0 ; i < rows ; i++ )    {
+         for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
+           
+           // Note: rotating the coordinated during the sum up
+           
+           Int_t rotVal = (phiPoints/phiPOC)*phiiC;
+           Int_t ipR = -1;
+           
+           if ((ip%phiPoints)>=rotVal) {
+             ipR = ip-rotVal;
+           } else {
+             ipR = ip+(phiPoints-rotVal);
+           }
+           
+           // unfortunately, the lookup tables were produced to be faster for phi symmetric charges
+           // This will be the most frequent usage (hopefully)
+           // That's why we have to do this here and not outside the loop ...
+           
+           TMatrixD &erA   =  *arrayofErA[k]  ;
+           TMatrixD &ephiA =  *arrayofEphiA[k];
+           TMatrixD &dEzA  =  *arrayofdEzA[k]   ;
+           
+           TMatrixD &erC   =  *arrayofErC[k]  ;
+           TMatrixD &ephiC =  *arrayofEphiC[k];
+           TMatrixD &dEzC  =  *arrayofdEzC[k]   ;
+       
+           // Sum up - Efield values in [V/m] -> transition to [V/cm]
+           erA(i,j) += ((*bEr)(ipR)) * weightA /100;
+           erC(i,j) += ((*bEr)(ipR)) * weightC /100;
+           if (!flagRadSym) {
+             ephiA(i,j) += ((*bEphi)(ipR)) * weightA; // [V/rad]
+             ephiC(i,j) += ((*bEphi)(ipR)) * weightC; // [V/rad]
+           }
+           dEzA(i,j)  += ((*bEz)(ipR)) * weightA /100;
+           dEzC(i,j)  += ((*bEz)(ipR)) * weightC /100;
+
+           // increase the counter
+           ip++;
+         }
+       }
+      } // end coordinate loop
+
+    } // end phi-POC summation (phiiC)
+   
+
+    // printf("POC: (r,phi,z) = (%lf %lf %lf) | weight(A,C): %03.1lf %03.1lf\n",r0,phi0,z0, weightA, weightC);
+
+  } // end POC loop
+  AliInfo("Division and integration");
+
+  // -------------------------------------------------------------------------------
+  // Division by the Ez (drift) field and integration along z
+
+  Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = Electric Field (V/cm) Magnitude ~ -400 V/cm;
+
+  for ( Int_t k = 0 ; k < phiSlices ; k++ ) { // phi loop
+
+    // matrices holding the solution - summation of POC charges
+    TMatrixD &erA   =  *arrayofErA[k]  ;
+    TMatrixD &ephiA =  *arrayofEphiA[k];
+    TMatrixD &dezA  =  *arrayofdEzA[k]   ;
+    TMatrixD &erC   =  *arrayofErC[k]  ;
+    TMatrixD &ephiC =  *arrayofEphiC[k];
+    TMatrixD &dezC  =  *arrayofdEzC[k]   ;
+
+    // matrices which contain the integrated fields (divided by the drift field)
+    TMatrixD &erOverEzA   =  *arrayofEroverEzA[k]  ;
+    TMatrixD &ephiOverEzA =  *arrayofEphioverEzA[k];
+    TMatrixD &deltaEzA    =  *arrayofDeltaEzA[k];
+    TMatrixD &erOverEzC   =  *arrayofEroverEzC[k]  ;
+    TMatrixD &ephiOverEzC =  *arrayofEphioverEzC[k];
+    TMatrixD &deltaEzC    =  *arrayofDeltaEzC[k];    
+    
+    for ( Int_t i = 0 ; i < rows ; i++ )    { // r loop
+      for ( Int_t j = columns-1 ; j >= 0 ; j-- ) {// z loop 
+       // Count backwards to facilitate integration over Z
+
+       Int_t index = 1 ; // Simpsons rule if N=odd.If N!=odd then add extra point by trapezoidal rule.  
+
+       erOverEzA(i,j) = 0; ephiOverEzA(i,j) = 0;  deltaEzA(i,j) = 0;
+       erOverEzC(i,j) = 0; ephiOverEzC(i,j) = 0;  deltaEzC(i,j) = 0;
+
+       for ( Int_t m = j ; m < columns ; m++ ) { // integration
+
+         erOverEzA(i,j)   += index*(gridSizeZ/3.0)*erA(i,m)/(-1*ezField) ;
+         erOverEzC(i,j)   += index*(gridSizeZ/3.0)*erC(i,m)/(-1*ezField)  ;
+         if (!flagRadSym) {
+           ephiOverEzA(i,j) += index*(gridSizeZ/3.0)*ephiA(i,m)/(-1*ezField)  ;
+           ephiOverEzC(i,j) += index*(gridSizeZ/3.0)*ephiC(i,m)/(-1*ezField)  ;
+         }
+         deltaEzA(i,j)    += index*(gridSizeZ/3.0)*dezA(i,m) ;
+         deltaEzC(i,j)    += index*(gridSizeZ/3.0)*dezC(i,m) ;
+
+         if ( index != 4 )  index = 4; else index = 2 ;
+
+       }
+
+       if ( index == 4 ) {
+         erOverEzA(i,j)   -= (gridSizeZ/3.0)*erA(i,columns-1)/(-1*ezField) ;
+         erOverEzC(i,j)   -= (gridSizeZ/3.0)*erC(i,columns-1)/(-1*ezField) ;
+         if (!flagRadSym) {
+           ephiOverEzA(i,j) -= (gridSizeZ/3.0)*ephiA(i,columns-1)/(-1*ezField) ;
+           ephiOverEzC(i,j) -= (gridSizeZ/3.0)*ephiC(i,columns-1)/(-1*ezField) ;
+         }
+         deltaEzA(i,j)    -= (gridSizeZ/3.0)*dezA(i,columns-1) ;
+         deltaEzC(i,j)    -= (gridSizeZ/3.0)*dezC(i,columns-1) ;
+       }
+       if ( index == 2 ) {
+         erOverEzA(i,j)   += (gridSizeZ/3.0)*(0.5*erA(i,columns-2)-2.5*erA(i,columns-1))/(-1*ezField) ;
+         erOverEzC(i,j)   += (gridSizeZ/3.0)*(0.5*erC(i,columns-2)-2.5*erC(i,columns-1))/(-1*ezField) ;
+         if (!flagRadSym) {
+           ephiOverEzA(i,j) += (gridSizeZ/3.0)*(0.5*ephiA(i,columns-2)-2.5*ephiA(i,columns-1))/(-1*ezField) ;
+           ephiOverEzC(i,j) += (gridSizeZ/3.0)*(0.5*ephiC(i,columns-2)-2.5*ephiC(i,columns-1))/(-1*ezField) ;
+         }
+         deltaEzA(i,j)    += (gridSizeZ/3.0)*(0.5*dezA(i,columns-2)-2.5*dezA(i,columns-1)) ;
+         deltaEzC(i,j)    += (gridSizeZ/3.0)*(0.5*dezC(i,columns-2)-2.5*dezC(i,columns-1)) ;
+       }
+       if ( j == columns-2 ) {
+         erOverEzA(i,j)   = (gridSizeZ/3.0)*(1.5*erA(i,columns-2)+1.5*erA(i,columns-1))/(-1*ezField) ;
+         erOverEzC(i,j)   = (gridSizeZ/3.0)*(1.5*erC(i,columns-2)+1.5*erC(i,columns-1))/(-1*ezField) ;
+         if (!flagRadSym) {
+           ephiOverEzA(i,j) = (gridSizeZ/3.0)*(1.5*ephiA(i,columns-2)+1.5*ephiA(i,columns-1))/(-1*ezField) ;
+           ephiOverEzC(i,j) = (gridSizeZ/3.0)*(1.5*ephiC(i,columns-2)+1.5*ephiC(i,columns-1))/(-1*ezField) ;
+         }
+         deltaEzA(i,j)    = (gridSizeZ/3.0)*(1.5*dezA(i,columns-2)+1.5*dezA(i,columns-1)) ;
+         deltaEzC(i,j)    = (gridSizeZ/3.0)*(1.5*dezC(i,columns-2)+1.5*dezC(i,columns-1)) ;
+       }
+       if ( j == columns-1 ) {
+         erOverEzA(i,j)   = 0;    erOverEzC(i,j)   = 0;
+         if (!flagRadSym) {
+           ephiOverEzA(i,j) = 0;  ephiOverEzC(i,j) = 0;
+         }
+         deltaEzA(i,j)    = 0;    deltaEzC(i,j)    = 0;
+       }
+      }
+    }
+
+  }
+  
+  AliInfo("Interpolation to Standard grid");
+
+  // -------------------------------------------------------------------------------
+  // Interpolate results onto the standard grid which is used for all AliTPCCorrections
+
+  const Int_t order  = 1  ;  // Linear interpolation = 1, Quadratic = 2  
+
+  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]) ;  // z position is symmetric
+    
+      for ( Int_t i = 0 ; i < kNR ; i++ ) { 
+       r = fgkRList[i] ;
+
+       // Interpolate Lookup tables onto standard grid
+       if (fgkZList[j]>0) {
+         erOverEz(i,j)   = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices, 
+                                              rlist, zedlist, philist, arrayofEroverEzA  );
+         ephiOverEz(i,j) = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
+                                              rlist, zedlist, philist, arrayofEphioverEzA);
+         deltaEz(i,j)    = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
+                                              rlist, zedlist, philist, arrayofDeltaEzA   );
+       } else {
+         erOverEz(i,j)   = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices, 
+                                              rlist, zedlist, philist, arrayofEroverEzC  );
+         ephiOverEz(i,j) = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
+                                              rlist, zedlist, philist, arrayofEphioverEzC);
+         deltaEz(i,j)  = - Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
+                                              rlist, zedlist, philist, arrayofDeltaEzC   );
+         // negative coordinate system on C side
+       }
+
+      } // end r loop
+    } // end z loop
+  } // end phi loop
+
+  
+  // clear the temporary arrays lists
+  for ( Int_t k = 0 ; k < phiSlices ; k++ )  {
+
+    delete arrayofErA[k];  
+    delete arrayofEphiA[k];
+    delete arrayofdEzA[k];
+    delete arrayofErC[k];  
+    delete arrayofEphiC[k];
+    delete arrayofdEzC[k];
+
+    delete arrayofEroverEzA[k];  
+    delete arrayofEphioverEzA[k];
+    delete arrayofDeltaEzA[k];
+    delete arrayofEroverEzC[k];  
+    delete arrayofEphioverEzC[k];
+    delete arrayofDeltaEzC[k];
+
+  }
+
+
+  fInitLookUp = kTRUE;
+
+}
+
+
+void AliTPCSpaceCharge3D::SetSCDataFileName(char *const fname) {
+  //
+  // Set & load the Space charge density distribution from a file 
+  // (linear interpolation onto a standard grid)
+  //
+
+  fSCDataFileName = fname;
+
+  TFile *f = new TFile(fSCDataFileName,"READ");
+  if (!f) { 
+    AliError(Form("File %s, which should contain the space charge distribution, could not be found",
+                 fSCDataFileName));
+    return;
+  }
+  TH3F *density = (TH3F*) f->Get("SpaceChargeDistribution");
+  if (!density) { 
+    AliError(Form("The indicated file (%s) does not contain a histogram called %s",
+                 fSCDataFileName,"SpaceChargeDistribution"));
+    return;
+  }
+
+  Double_t  r, phi, z ;
+  for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
+    phi = fgkPhiList[k] ;
+    TMatrixD &scdensity   =  *fSCdensityDistribution[k]  ;
+    for ( Int_t j = 0 ; j < kNZ ; j++ ) {
+      z = TMath::Abs(fgkZList[j]) ;  // z position is symmetric
+      for ( Int_t i = 0 ; i < kNR ; i++ ) { 
+       r = fgkRList[i] ;
+       scdensity(i,j) = density->Interpolate(r,phi,z); // quite slow
+       //      printf("%lf %lf %lf: %lf\n",r,phi,z,scdensity(i,j));
+      }
+    }
+  }
+
+  f->Close();
+
+  fInitLookUp = kFALSE;
+
+  
+}
+
+
+Float_t  AliTPCSpaceCharge3D::GetSpaceChargeDensity(Float_t r, Float_t phi, Float_t z) {
+  //
+  // returns the (input) space charge density at a given point according 
+  // Note: input in [cm], output in [C/m^3/e0] !!
+  //
+
+  if (!fSCdensityDistribution) {
+    printf("Scheiss is nuul - argg\n");
+    return 0.;
+  }
+
+  while (phi<0) phi += TMath::TwoPi();
+  while (phi>TMath::TwoPi()) phi -= TMath::TwoPi();
+
+
+  // Float_t sc =fSCdensityDistribution->Interpolate(r0,phi0,z0);
+  Int_t order = 1; //
+  Float_t sc = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi, 
+                                 fgkRList, fgkZList, fgkPhiList, fSCdensityDistribution );
+
+  //  printf("%lf %lf %lf: %lf\n",r,phi,z,sc);
+  
+  return sc;
+}
+
+
+TH2F * AliTPCSpaceCharge3D::CreateHistoSCinXY(Float_t z, Int_t nx, Int_t ny) {
+  //
+  // return a simple histogramm containing the space charge distribution (input for the calculation)
+  //
+
+  TH2F *h=CreateTH2F("spaceCharge",GetTitle(),"x [cm]","y [cm]","#rho_{sc} [C/m^{3}/e_{0}]",
+                    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 r = TMath::Sqrt(xp*xp+yp*yp);
+      Float_t phi = TMath::ATan2(yp,xp); 
+   
+      if (85.<=r && r<=250.) {
+       Float_t sc = GetSpaceChargeDensity(r,phi,z)/fgke0; // in [C/m^3/e0]
+       h->SetBinContent(ix,iy,sc); 
+      } else {
+       h->SetBinContent(ix,iy,0.);
+      }
+    }
+  }
+  
+  return h;
+} 
+
+TH2F * AliTPCSpaceCharge3D::CreateHistoSCinZR(Float_t phi, Int_t nz, Int_t nr) {
+  //
+  // return a simple histogramm containing the space charge distribution (input for the calculation)
+  //
+
+  TH2F *h=CreateTH2F("spaceCharge",GetTitle(),"z [cm]","r [cm]","#rho_{sc} [C/m^{3}/e_{0}]",
+                    nz,-250.,250.,nr,85.,250.);
+
+  for (Int_t ir=1;ir<=nr;++ir) {
+    Float_t r = h->GetYaxis()->GetBinCenter(ir);
+    for (Int_t iz=1;iz<=nz;++iz) {
+      Float_t z = h->GetXaxis()->GetBinCenter(iz);
+      Float_t sc = GetSpaceChargeDensity(r,phi,z)/fgke0; // in [C/m^3/e0]
+      h->SetBinContent(iz,ir,sc);
+    }
+  }
+
+  return h;
+} 
+
+void AliTPCSpaceCharge3D::WriteChargeDistributionToFile(const char* fname) {
+  //
+  // Example on how to write a Space charge distribution into a File
+  //  (see below: estimate from scaling STAR measurements to Alice)
+  //
+
+  TFile *f = new TFile(fname,"RECREATE");
+
+  // some grid, not too course
+  Int_t nr = 72;
+  Int_t nphi = 180;
+  Int_t nz = 250;
+
+  Double_t dr = (fgkOFCRadius-fgkIFCRadius)/(nr+1);
+  Double_t dphi = TMath::TwoPi()/(nphi+1);
+  Double_t dz = 500./(nz+1);
+  Double_t safty = 0.; // due to a root bug which does not interpolate the boundary (first and last bin) correctly
+
+  TH3F *histo = new TH3F("charge","charge",
+                        nr,fgkIFCRadius-dr-safty,fgkOFCRadius+dr+safty,
+                        nphi,0-dphi-safty,TMath::TwoPi()+dphi+safty,
+                        nz,-250-dz-safty,250+dz+safty);
+  for (Int_t ir=1;ir<=nr;++ir) {
+    Double_t rp = histo->GetXaxis()->GetBinCenter(ir);
+    for (Int_t iphi=1;iphi<=nphi;++iphi) {
+      Double_t phip = histo->GetYaxis()->GetBinCenter(iphi);
+      for (Int_t iz=1;iz<=nz;++iz) {
+       Double_t zp = histo->GetZaxis()->GetBinCenter(iz);
+
+
+       zp = TMath::Abs(zp); // symmetric on both sides of the TPC
+       
+       // recalculation to meter
+       Double_t lZ = 2.5; // approx. TPC drift length
+       Double_t rpM = rp/100.; // in [m]
+       Double_t zpM = zp/100.; // in [m]
+       
+       // setting of mb multiplicity and Interaction rate
+       Double_t multiplicity = 950;
+       Double_t intRate = 8000;
+
+       // calculation of "scaled" parameters
+       Double_t a = multiplicity*intRate/79175;
+       Double_t b = a/lZ;
+       
+       Double_t charge = (a - b * zpM)/(rpM*rpM); // charge in [C/m^3/e0]
+       
+       // some 'arbitrary' GG leaks
+       if (  (phip<(TMath::Pi()/9*3) && phip>(TMath::Pi()/9*2) ) ) { //||
+         //          (phip<(TMath::Pi()/9*7) && phip>(TMath::Pi()/9*6) ) ||
+         //          (phip<(TMath::Pi()/9*16) && phip>(TMath::Pi()/9*13) ) ){
+         if (rp>130&&rp<135) 
+           charge = 300; 
+       }
+
+       //      printf("%lf %lf %lf: %lf\n",rp,phip,zp,charge);
+       charge = charge*fgke0; // [C/m^3]
+
+       histo->SetBinContent(ir,iphi,iz,charge); 
+      }
+    }
+  }
+
+
+  histo->Write("SpaceChargeDistribution");
+  f->Close();
+  
+}
+
+
+void AliTPCSpaceCharge3D::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 effect with arbitrary 3D charge density (as input).\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 InitSpaceCharge3DDistortion() ...");
+
+}
diff --git a/TPC/AliTPCSpaceCharge3D.h b/TPC/AliTPCSpaceCharge3D.h
new file mode 100644 (file)
index 0000000..9de2b10
--- /dev/null
@@ -0,0 +1,82 @@
+#ifndef ALITPCSPACECHARGE3D_H
+#define ALITPCSPACECHARGE3D_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+////////////////////////////////////////////////////////////////////////////
+//                                                                        //
+// AliTPCSpaceCharge3D class                                           //
+// date: 19/06/2010                                                       //
+// Authors: Stefan Rossegger                                              //
+////////////////////////////////////////////////////////////////////////////
+
+#include "AliTPCCorrection.h"
+class TH3F;
+
+class AliTPCSpaceCharge3D : public AliTPCCorrection {
+public:
+  AliTPCSpaceCharge3D();
+  virtual ~AliTPCSpaceCharge3D();
+
+  // initialization and update functions
+  virtual void Init();
+  virtual void Update(const TTimeStamp &timeStamp);
+
+  // common setters and getters for tangled ExB effect
+  virtual void SetOmegaTauT1T2(Float_t omegaTau,Float_t t1,Float_t t2) {
+    fT1=t1; fT2=t2;
+    const Double_t wt0=t2*omegaTau;     fC0=1./(1.+wt0*wt0);
+    const Double_t wt1=t1*omegaTau;     fC1=wt1/(1.+wt1*wt1);
+  };
+  void SetC0C1(Float_t c0,Float_t c1) {fC0=c0;fC1=c1;} // CAUTION: USE WITH CARE
+  Float_t GetC0() const {return fC0;}
+  Float_t GetC1() const {return fC1;}
+
+  // setters and getters
+  void SetCorrectionFactor(Float_t correctionFactor) {fCorrectionFactor=correctionFactor;}
+  Float_t GetCorrectionFactor() const {return fCorrectionFactor;}
+
+  void  SetSCDataFileName(char *const fname);
+  char* GetSCDataFileName() const { return fSCDataFileName; }
+
+  void InitSpaceCharge3DDistortion();
+
+  Float_t GetSpaceChargeDensity(Float_t r, Float_t phi, Float_t z);
+  TH2F* CreateHistoSCinXY(Float_t z, Int_t nx=100, Int_t ny=100);
+  TH2F* CreateHistoSCinZR(Float_t phi, Int_t nz=100, Int_t nr=100);
+
+
+  void WriteChargeDistributionToFile(const char* fname = "SC-Alice.root");
+
+  virtual void Print(const Option_t* option="") const;
+
+protected:
+  virtual void GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]);
+
+private:
+
+  AliTPCSpaceCharge3D(const AliTPCSpaceCharge3D &);               // not implemented
+  AliTPCSpaceCharge3D &operator=(const AliTPCSpaceCharge3D &);    // not implemented
+
+  Float_t fC0; // coefficient C0                 (compare Jim Thomas's notes for definitions)
+  Float_t fC1; // coefficient C1                 (compare Jim Thomas's notes for definitions)
+  Float_t fCorrectionFactor;       // Space Charge Correction factor in comparison to initialized
+                                   // look up table which was created for M_mb = 900 and IR = 3000
+                                   // compare Internal Note Nr: ???
+
+  Bool_t fInitLookUp;                 // flag to check it the Look Up table was created
+
+  TMatrixD *fLookUpErOverEz[kNPhi];   // Array to store electric field integral (int Er/Ez)
+  TMatrixD *fLookUpEphiOverEz[kNPhi]; // Array to store electric field integral (int Er/Ez)
+  TMatrixD *fLookUpDeltaEz[kNPhi];    // Array to store electric field integral (int Er/Ez)
+
+  char *fSCDataFileName;         // file which contains the space charge distribution
+  char *fSCLookUpPOCsFileName;   // filename of the precalculated lookup tables (for individual voxels)
+
+  TMatrixD *fSCdensityDistribution[kNPhi]; // space charge distribution
+  ClassDef(AliTPCSpaceCharge3D,1); 
+};
+
+#endif