TPC/AliTPCSpaceCharge3D.cxx

   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 // AliTPCSpaceCharge3D class                                                //
  19 // The class calculates the space point distortions due to a space charge   //
  20 // effect ....                                                              //
  21 // Method of calculation:                                                   //
  22 // The analytical solution for the poisson problem in 3D (cylindrical coord)//
  23 // is used in form of look up tables. PieceOfCake (POC) voxel were pre-     //
  24 // calculated and can be sumed up (weighted) according to the what is needed//
  25 //                                                                          //
  26 // The class allows "effective Omega Tau" corrections.                      //
  27 //                                                                          //
  28 // NOTE: This class is not  capable of calculation z distortions due to     //
  29 //       drift velocity change in dependence of the electric field!!!       //
  30 //                                                                          //
  31 // date: 19/06/2010                                                         //
  32 // Authors: Stefan Rossegger                                                //
  33 //                                                                          //
  34 // Example usage:                                                           //
  35 //////////////////////////////////////////////////////////////////////////////
  36
  37 #include "AliMagF.h"
  38 #include "TGeoGlobalMagField.h"
  39 #include "AliTPCcalibDB.h"
  40 #include "AliTPCParam.h"
  41 #include "AliLog.h"
  42 #include "TH2F.h"
  43 #include "TH3F.h"
  44 #include "TFile.h"
  45 #include "TVector.h"
  46 #include "TMatrix.h"
  47 #include "TMatrixD.h"
  48
  49 #include "TMath.h"
  50 #include "AliTPCROC.h"
  51 #include "AliTPCSpaceCharge3D.h"
  52
  53 ClassImp(AliTPCSpaceCharge3D)
  54
  55 AliTPCSpaceCharge3D::AliTPCSpaceCharge3D()
  56   : AliTPCCorrection("SpaceCharge3D","Space Charge - 3D"),
  57     fC0(0.),fC1(0.),
  58     fCorrectionFactor(1.),
  59     fInitLookUp(kFALSE),
  60     fSCDataFileName(""),
  61     fSCLookUpPOCsFileName3D(""),
  62     fSCLookUpPOCsFileNameRZ(""),
  63     fSCLookUpPOCsFileNameRPhi(""),
  64     fSCdensityInRZ(0),
  65     fSCdensityInRPhiA(0),
  66     fSCdensityInRPhiC(0)
  67 {
  68   //
  69   // default constructor
  70   //
  71
  72   // Array which will contain the solution according to the setted charge density distribution
  73   // see InitSpaceCharge3DDistortion() function
  74   for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
  75     fLookUpErOverEz[k]   =  new TMatrixD(kNR,kNZ);
  76     fLookUpEphiOverEz[k] =  new TMatrixD(kNR,kNZ);
  77     fLookUpDeltaEz[k]    =  new TMatrixD(kNR,kNZ);
  78     fSCdensityDistribution[k] = new TMatrixD(kNR,kNZ);
  79   }
  80   fSCdensityInRZ   = new TMatrixD(kNR,kNZ);
  81   fSCdensityInRPhiA = new TMatrixD(kNR,kNPhi);
  82   fSCdensityInRPhiC = new TMatrixD(kNR,kNPhi);
  83
  84   // location of the precalculated look up tables
  85
  86   fSCLookUpPOCsFileName3D="$(ALICE_ROOT)/TPC/Calib/maps/sc_3D_raw_18-18-26_17p-18p-25p-MN30.root"; // rough estimate
  87   fSCLookUpPOCsFileNameRZ="$(ALICE_ROOT)/TPC/Calib/maps/sc_radSym_35-01-51_34p-01p-50p_MN60.root";
  88   fSCLookUpPOCsFileNameRPhi="$(ALICE_ROOT)/TPC/Calib/maps/sc_cChInZ_35-144-26_34p-18p-01p-MN30.root";
  89   //  fSCLookUpPOCsFileNameRPhi="$(ALICE_ROOT)/TPC/Calib/maps/sc_cChInZ_35-36-26_34p-18p-01p-MN40.root";
  90
  91
  92
  93   // standard location of the space charge distibution ... can be changes
  94   fSCDataFileName="$(ALICE_ROOT)/TPC/Calib/maps/sc_3D_distribution_Sim.root";
  95
  96   //  SetSCDataFileName(fSCDataFileName.Data()); // should be done by the user
  97
  98
  99 }
 100
 101 AliTPCSpaceCharge3D::~AliTPCSpaceCharge3D() {
 102   //
 103   // default destructor
 104   //
 105
 106   for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
 107     delete fLookUpErOverEz[k];
 108     delete fLookUpEphiOverEz[k];
 109     delete fLookUpDeltaEz[k];
 110     delete fSCdensityDistribution[k];
 111   }
 112   delete fSCdensityInRZ;
 113   delete fSCdensityInRPhiA;
 114   delete fSCdensityInRPhiC;
 115
 116 }
 117
 118
 119 void AliTPCSpaceCharge3D::Init() {
 120   //
 121   // Initialization funtion
 122   //
 123
 124   AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
 125   if (!magF) AliError("Magneticd field - not initialized");
 126   Double_t bzField = magF->SolenoidField()/10.; //field in T
 127   AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
 128   if (!param) AliError("Parameters - not initialized");
 129   Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us]   // From dataBase: to be updated: per second (ideally)
 130   Double_t ezField = 400; // [V/cm]   // to be updated: never (hopefully)
 131   Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
 132   // Correction Terms for effective omegaTau; obtained by a laser calibration run
 133   SetOmegaTauT1T2(wt,fT1,fT2);
 134
 135   InitSpaceCharge3DDistortion(); // fill the look up table
 136 }
 137
 138 void AliTPCSpaceCharge3D::Update(const TTimeStamp &/*timeStamp*/) {
 139   //
 140   // Update function
 141   //
 142   AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
 143   if (!magF) AliError("Magneticd field - not initialized");
 144   Double_t bzField = magF->SolenoidField()/10.; //field in T
 145   AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
 146   if (!param) AliError("Parameters - not initialized");
 147   Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us]  // From dataBase: to be updated: per second (ideally)
 148   Double_t ezField = 400; // [V/cm]   // to be updated: never (hopefully)
 149   Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
 150   // Correction Terms for effective omegaTau; obtained by a laser calibration run
 151   SetOmegaTauT1T2(wt,fT1,fT2);
 152
 153   //  SetCorrectionFactor(1.); // should come from some database
 154
 155 }
 156
 157
 158 void AliTPCSpaceCharge3D::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {
 159   //
 160   // Calculates the correction due the Space Charge effect within the TPC drift volume
 161   //
 162
 163   if (!fInitLookUp) {
 164     AliInfo("Lookup table was not initialized! Performing the inizialisation now ...");
 165     InitSpaceCharge3DDistortion();
 166   }
 167
 168   Int_t   order     = 1 ;    // FIXME: hardcoded? Linear interpolation = 1, Quadratic = 2
 169
 170   Double_t intEr, intEphi, intdEz ;
 171   Double_t r, phi, z ;
 172   Int_t    sign;
 173
 174   r      =  TMath::Sqrt( x[0]*x[0] + x[1]*x[1] ) ;
 175   phi    =  TMath::ATan2(x[1],x[0]) ;
 176   if ( phi < 0 ) phi += TMath::TwoPi() ;                   // Table uses phi from 0 to 2*Pi
 177   z      =  x[2] ;                                         // Create temporary copy of x[2]
 178
 179   if ( (roc%36) < 18 ) {
 180     sign =  1;       // (TPC A side)
 181   } else {
 182     sign = -1;       // (TPC C side)
 183   }
 184
 185   if ( sign==1  && z <  fgkZOffSet ) z =  fgkZOffSet;    // Protect against discontinuity at CE
 186   if ( sign==-1 && z > -fgkZOffSet ) z = -fgkZOffSet;    // Protect against discontinuity at CE
 187
 188
 189   if ( (sign==1 && z<0) || (sign==-1 && z>0) ) // just a consistency check
 190     AliError("ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!");
 191
 192   // Get the Er and Ephi field integrals plus the integral over DeltaEz
 193   intEr      = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
 194                                   fgkRList, fgkZList, fgkPhiList, fLookUpErOverEz  );
 195   intEphi    = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
 196                                   fgkRList, fgkZList, fgkPhiList, fLookUpEphiOverEz);
 197   intdEz = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
 198                                   fgkRList, fgkZList, fgkPhiList, fLookUpDeltaEz   );
 199
 200   // Calculate distorted position
 201   if ( r > 0.0 ) {
 202     phi =  phi + fCorrectionFactor *( fC0*intEphi - fC1*intEr ) / r;
 203     r   =  r   + fCorrectionFactor *( fC0*intEr   + fC1*intEphi );
 204   }
 205   Double_t dz = intdEz * fCorrectionFactor * fgkdvdE;
 206
 207   // Calculate correction in cartesian coordinates
 208   dx[0] = - (r * TMath::Cos(phi) - x[0]);
 209   dx[1] = - (r * TMath::Sin(phi) - x[1]);
 210   dx[2] = - dz;  // z distortion - (scaled with driftvelocity dependency on the Ez field and the overall scaling factor)
 211
 212 }
 213
 214 void AliTPCSpaceCharge3D::InitSpaceCharge3DDistortion() {
 215  //
 216   // Initialization of the Lookup table which contains the solutions of the
 217   // "space charge" (poisson) problem - Faster and more accureate
 218   //
 219   // Method: Weighted sum-up of the different fields within the look up table
 220   // but using two lookup tables with higher granularity in the (r,z) and the (rphi)- plane to emulate
 221   // more realistic space charges. (r,z) from primary ionisation. (rphi) for possible Gating leaks
 222
 223   if (fInitLookUp) {
 224     AliInfo("Lookup table was already initialized!  Doing it again anyway ...");
 225     //    return;
 226   }
 227
 228   // ------------------------------------------------------------------------------------------------------
 229   // step 1: lookup table in rz, fine grid, radial symetric, to emulate primary ionization
 230
 231   AliInfo("Step 1: Preparation of the weighted look-up tables.");
 232
 233   // lookup table in rz, fine grid
 234
 235   TFile *fZR = new TFile(fSCLookUpPOCsFileNameRZ.Data(),"READ");
 236   if ( !fZR ) {
 237     AliError("Precalculated POC-looup-table in ZR could not be found");
 238     return;
 239   }
 240
 241   // units are in [m]
 242   TVector *gridf1  = (TVector*) fZR->Get("constants");
 243   TVector &grid1 = *gridf1;
 244   TMatrix *coordf1  = (TMatrix*) fZR->Get("coordinates");
 245   TMatrix &coord1 = *coordf1;
 246   TMatrix *coordPOCf1  = (TMatrix*) fZR->Get("POCcoord");
 247   TMatrix &coordPOC1 = *coordPOCf1;
 248
 249   Int_t rows      = (Int_t)grid1(0);   // number of points in r direction  - from RZ or RPhi table
 250   Int_t phiSlices = (Int_t)grid1(1);   // number of points in phi          - from RPhi table
 251   Int_t columns   = (Int_t)grid1(2);   // number of points in z direction  - from RZ table
 252
 253   Float_t gridSizeR   =  (fgkOFCRadius-fgkIFCRadius)/(rows-1);   // unit in [cm]
 254   Float_t gridSizeZ   =  fgkTPCZ0/(columns-1);                  // unit in [cm]
 255
 256   // temporary matrices needed for the calculation // for rotational symmetric RZ table, phislices is 1
 257
 258   TMatrixD *arrayofErA[phiSlices], *arrayofdEzA[phiSlices];
 259   TMatrixD *arrayofErC[phiSlices], *arrayofdEzC[phiSlices];
 260
 261   TMatrixD *arrayofEroverEzA[phiSlices], *arrayofDeltaEzA[phiSlices];
 262   TMatrixD *arrayofEroverEzC[phiSlices], *arrayofDeltaEzC[phiSlices];
 263
 264
 265   for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
 266
 267     arrayofErA[k]   =   new TMatrixD(rows,columns) ;
 268     arrayofdEzA[k]  =   new TMatrixD(rows,columns) ;
 269     arrayofErC[k]   =   new TMatrixD(rows,columns) ;
 270     arrayofdEzC[k]  =   new TMatrixD(rows,columns) ;
 271
 272     arrayofEroverEzA[k]   =   new TMatrixD(rows,columns) ;
 273     arrayofDeltaEzA[k]    =   new TMatrixD(rows,columns) ;
 274     arrayofEroverEzC[k]   =   new TMatrixD(rows,columns) ;
 275     arrayofDeltaEzC[k]    =   new TMatrixD(rows,columns) ;
 276
 277     // zero initialization not necessary, it is done in the constructor of TMatrix
 278
 279   }
 280
 281   // list of points as used during sum up
 282   Double_t  rlist1[rows], zedlist1[columns];// , philist1[phiSlices];
 283   for ( Int_t i = 0 ; i < rows ; i++ )    {
 284     rlist1[i] = fgkIFCRadius + i*gridSizeR ;
 285     for ( Int_t j = 0 ; j < columns ; j++ ) {
 286       zedlist1[j]  = j * gridSizeZ ;
 287     }
 288   }
 289
 290   TTree *treePOC = (TTree*)fZR->Get("POCall");
 291
 292   TVector *bEr  = 0;   //TVector *bEphi= 0;
 293   TVector *bEz  = 0;
 294
 295   treePOC->SetBranchAddress("Er",&bEr);
 296   treePOC->SetBranchAddress("Ez",&bEz);
 297
 298
 299   // Read the complete tree and do a weighted sum-up over the POC configurations
 300   // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 301
 302   Int_t treeNumPOC = (Int_t)treePOC->GetEntries(); // Number of POC conf. in the look-up table
 303   Int_t ipC = 0; // POC Conf. counter (note: different to the POC number in the tree!)
 304
 305   for (Int_t itreepC=0; itreepC<treeNumPOC; itreepC++) { // ------------- loop over POC configurations in tree
 306
 307     treePOC->GetEntry(itreepC);
 308
 309     // center of the POC voxel in [meter]
 310     Double_t r0 = coordPOC1(ipC,0);
 311     Double_t phi0 = coordPOC1(ipC,1);
 312     Double_t z0 = coordPOC1(ipC,2);
 313
 314     ipC++; // POC configuration counter
 315
 316     // weights (charge density) at POC position on the A and C side (in C/m^3/e0)
 317     // note: coordinates are in [cm]
 318     Double_t weightA = GetSpaceChargeDensity(r0*100,phi0, z0*100, 1);  // partial load in r,z
 319     Double_t weightC = GetSpaceChargeDensity(r0*100,phi0,-z0*100, 1);  // partial load in r,z
 320
 321     // Summing up the vector components according to their weight
 322
 323     Int_t ip = 0;
 324     for ( Int_t j = 0 ; j < columns ; j++ ) {
 325       for ( Int_t i = 0 ; i < rows ; i++ )    {
 326         for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
 327
 328           // check wether the coordinates were screwed
 329           if (TMath::Abs((coord1(0,ip)*100-rlist1[i]))>1 ||
 330               TMath::Abs((coord1(2,ip)*100-zedlist1[j])>1)) {
 331             AliError("internal error: coordinate system was screwed during the sum-up");
 332             printf("sum-up: (r,z)=(%lf,%lf)\n",rlist1[i],zedlist1[j]);
 333             printf("lookup: (r,z)=(%lf,%lf)\n",coord1(0,ip)*100,coord1(2,ip)*100);
 334             AliError("Don't trust the results of the space charge calculation!");
 335           }
 336
 337           // unfortunately, the lookup tables were produced to be faster for phi symmetric charges
 338           // This will be the most frequent usage (hopefully)
 339           // That's why we have to do this here ...
 340
 341           TMatrixD &erA   =  *arrayofErA[k]  ;
 342           TMatrixD &dEzA  =  *arrayofdEzA[k]   ;
 343
 344           TMatrixD &erC   =  *arrayofErC[k]  ;
 345           TMatrixD &dEzC  =  *arrayofdEzC[k]   ;
 346
 347           // Sum up - Efield values in [V/m] -> transition to [V/cm]
 348           erA(i,j) += ((*bEr)(ip)) * weightA /100;
 349           erC(i,j) += ((*bEr)(ip)) * weightC /100;
 350           dEzA(i,j)  += ((*bEz)(ip)) * weightA /100;
 351           dEzC(i,j)  += ((*bEz)(ip)) * weightC /100;
 352
 353           // increase the counter
 354           ip++;
 355         }
 356       }
 357     } // end coordinate loop
 358   } // end POC loop
 359
 360
 361   // -------------------------------------------------------------------------------
 362   // Division by the Ez (drift) field and integration along z
 363
 364   //  AliInfo("Step 1: Division and integration");
 365
 366   Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = Electric Field (V/cm) Magnitude ~ -400 V/cm;
 367
 368   for ( Int_t k = 0 ; k < phiSlices ; k++ ) { // phi loop
 369
 370     // matrices holding the solution - summation of POC charges // see above
 371     TMatrixD &erA   =  *arrayofErA[k]  ;
 372     TMatrixD &dezA  =  *arrayofdEzA[k]   ;
 373     TMatrixD &erC   =  *arrayofErC[k]  ;
 374     TMatrixD &dezC  =  *arrayofdEzC[k]   ;
 375
 376     // matrices which will contain the integrated fields (divided by the drift field)
 377     TMatrixD &erOverEzA   =  *arrayofEroverEzA[k]  ;
 378     TMatrixD &deltaEzA    =  *arrayofDeltaEzA[k];
 379     TMatrixD &erOverEzC   =  *arrayofEroverEzC[k]  ;
 380     TMatrixD &deltaEzC    =  *arrayofDeltaEzC[k];
 381
 382     for ( Int_t i = 0 ; i < rows ; i++ )    { // r loop
 383       for ( Int_t j = columns-1 ; j >= 0 ; j-- ) {// z loop
 384         // Count backwards to facilitate integration over Z
 385
 386         Int_t index = 1 ; // Simpsons rule if N=odd.If N!=odd then add extra point
 387                           // by trapezoidal rule.
 388
 389         erOverEzA(i,j) = 0; //ephiOverEzA(i,j) = 0;
 390         deltaEzA(i,j) = 0;
 391         erOverEzC(i,j) = 0; //ephiOverEzC(i,j) = 0;
 392         deltaEzC(i,j) = 0;
 393
 394         for ( Int_t m = j ; m < columns ; m++ ) { // integration
 395
 396           erOverEzA(i,j)   += index*(gridSizeZ/3.0)*erA(i,m)/(-1*ezField) ;
 397           erOverEzC(i,j)   += index*(gridSizeZ/3.0)*erC(i,m)/(-1*ezField)  ;
 398           deltaEzA(i,j)    += index*(gridSizeZ/3.0)*dezA(i,m)/(-1) ;
 399           deltaEzC(i,j)    += index*(gridSizeZ/3.0)*dezC(i,m)/(-1) ;
 400
 401           if ( index != 4 )  index = 4; else index = 2 ;
 402
 403         }
 404
 405         if ( index == 4 ) {
 406           erOverEzA(i,j)   -= (gridSizeZ/3.0)*erA(i,columns-1)/(-1*ezField) ;
 407           erOverEzC(i,j)   -= (gridSizeZ/3.0)*erC(i,columns-1)/(-1*ezField) ;
 408           deltaEzA(i,j)    -= (gridSizeZ/3.0)*dezA(i,columns-1)/(-1) ;
 409           deltaEzC(i,j)    -= (gridSizeZ/3.0)*dezC(i,columns-1)/(-1) ;
 410         }
 411         if ( index == 2 ) {
 412           erOverEzA(i,j)   += (gridSizeZ/3.0)*(0.5*erA(i,columns-2)-2.5*erA(i,columns-1))/(-1*ezField) ;
 413           erOverEzC(i,j)   += (gridSizeZ/3.0)*(0.5*erC(i,columns-2)-2.5*erC(i,columns-1))/(-1*ezField) ;
 414           deltaEzA(i,j)    += (gridSizeZ/3.0)*(0.5*dezA(i,columns-2)-2.5*dezA(i,columns-1))/(-1) ;
 415           deltaEzC(i,j)    += (gridSizeZ/3.0)*(0.5*dezC(i,columns-2)-2.5*dezC(i,columns-1))/(-1) ;
 416         }
 417         if ( j == columns-2 ) {
 418           erOverEzA(i,j)   = (gridSizeZ/3.0)*(1.5*erA(i,columns-2)+1.5*erA(i,columns-1))/(-1*ezField) ;
 419           erOverEzC(i,j)   = (gridSizeZ/3.0)*(1.5*erC(i,columns-2)+1.5*erC(i,columns-1))/(-1*ezField) ;
 420           deltaEzA(i,j)    = (gridSizeZ/3.0)*(1.5*dezA(i,columns-2)+1.5*dezA(i,columns-1))/(-1) ;
 421           deltaEzC(i,j)    = (gridSizeZ/3.0)*(1.5*dezC(i,columns-2)+1.5*dezC(i,columns-1))/(-1) ;
 422         }
 423         if ( j == columns-1 ) {
 424           erOverEzA(i,j)   = 0;
 425           erOverEzC(i,j)   = 0;
 426           deltaEzA(i,j)    = 0;
 427           deltaEzC(i,j)    = 0;
 428         }
 429       }
 430     }
 431
 432   }
 433
 434   //  AliInfo("Step 1: Interpolation to Standard grid");
 435
 436   // -------------------------------------------------------------------------------
 437   // Interpolate results onto the standard grid which is used for all AliTPCCorrections classes
 438
 439   const Int_t order  = 1  ;  // Linear interpolation = 1, Quadratic = 2
 440
 441   Double_t  r, z;//phi, z ;
 442   for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
 443     //    phi = fgkPhiList[k] ;
 444
 445     // final lookup table
 446     TMatrixD &erOverEzFinal   =  *fLookUpErOverEz[k]  ;
 447     TMatrixD &deltaEzFinal    =  *fLookUpDeltaEz[k]   ;
 448
 449     // calculated and integrated tables - just one phi slice
 450     TMatrixD &erOverEzA   =  *arrayofEroverEzA[0]  ;
 451     TMatrixD &deltaEzA    =  *arrayofDeltaEzA[0];
 452     TMatrixD &erOverEzC   =  *arrayofEroverEzC[0]  ;
 453     TMatrixD &deltaEzC    =  *arrayofDeltaEzC[0];
 454
 455
 456     for ( Int_t j = 0 ; j < kNZ ; j++ ) {
 457
 458       z = TMath::Abs(fgkZList[j]) ;  // z position is symmetric
 459
 460       for ( Int_t i = 0 ; i < kNR ; i++ ) {
 461         r = fgkRList[i] ;
 462
 463         // Interpolate Lookup tables onto standard grid
 464         if (fgkZList[j]>0) {
 465           erOverEzFinal(i,j)   = Interpolate2DTable(order, r, z, rows, columns, rlist1, zedlist1, erOverEzA  );
 466           deltaEzFinal(i,j)    = Interpolate2DTable(order, r, z, rows, columns, rlist1, zedlist1, deltaEzA   );
 467         } else {
 468           erOverEzFinal(i,j)   = Interpolate2DTable(order, r, z, rows, columns, rlist1, zedlist1, erOverEzC  );
 469           deltaEzFinal(i,j)    = - Interpolate2DTable(order, r, z, rows, columns, rlist1, zedlist1, deltaEzC   );
 470           // negative coordinate system on C side
 471         }
 472
 473       } // end r loop
 474     } // end z loop
 475   } // end phi loop
 476
 477
 478   // clear the temporary arrays lists
 479   for ( Int_t k = 0 ; k < phiSlices ; k++ )  {
 480
 481     delete arrayofErA[k];
 482     delete arrayofdEzA[k];
 483     delete arrayofErC[k];
 484     delete arrayofdEzC[k];
 485
 486     delete arrayofEroverEzA[k];
 487     delete arrayofDeltaEzA[k];
 488     delete arrayofEroverEzC[k];
 489     delete arrayofDeltaEzC[k];
 490
 491   }
 492
 493   fZR->Close();
 494
 495   // ------------------------------------------------------------------------------------------------------
 496   // Step 2: Load and sum up lookup table in rphi, fine grid, to emulate for example a GG leak
 497
 498   //  AliInfo("Step 2: Preparation of the weighted look-up table");
 499
 500   TFile *fRPhi = new TFile(fSCLookUpPOCsFileNameRPhi.Data(),"READ");
 501   if ( !fRPhi ) {
 502     AliError("Precalculated POC-looup-table in RPhi could not be found");
 503     return;
 504   }
 505
 506   // units are in [m]
 507   TVector *gridf2  = (TVector*) fRPhi->Get("constants");
 508   TVector &grid2 = *gridf2;
 509   TMatrix *coordf2  = (TMatrix*) fRPhi->Get("coordinates");
 510   TMatrix &coord2 = *coordf2;
 511   TMatrix *coordPOCf2  = (TMatrix*) fRPhi->Get("POCcoord");
 512   TMatrix &coordPOC2 = *coordPOCf2;
 513
 514   rows      = (Int_t)grid2(0);   // number of points in r direction
 515   phiSlices = (Int_t)grid2(1);   // number of points in phi
 516   columns   = (Int_t)grid2(2);   // number of points in z direction
 517
 518   gridSizeR   =  (fgkOFCRadius-fgkIFCRadius)/(rows-1);   // unit in [cm]
 519   Float_t gridSizePhi =  TMath::TwoPi()/phiSlices;         // unit in [rad]
 520   gridSizeZ   =  fgkTPCZ0/(columns-1);                  // unit in [cm]
 521
 522   // list of points as used during sum up
 523   Double_t  rlist2[rows], philist2[phiSlices], zedlist2[columns];
 524   for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
 525     philist2[k] =  gridSizePhi * k;
 526     for ( Int_t i = 0 ; i < rows ; i++ )    {
 527       rlist2[i] = fgkIFCRadius + i*gridSizeR ;
 528       for ( Int_t j = 0 ; j < columns ; j++ ) {
 529         zedlist2[j]  = j * gridSizeZ ;
 530       }
 531     }
 532   } // only done once
 533
 534   // temporary matrices needed for the calculation
 535
 536   TMatrixD *arrayofErA2[phiSlices], *arrayofEphiA2[phiSlices], *arrayofdEzA2[phiSlices];
 537   TMatrixD *arrayofErC2[phiSlices], *arrayofEphiC2[phiSlices], *arrayofdEzC2[phiSlices];
 538
 539   TMatrixD *arrayofEroverEzA2[phiSlices], *arrayofEphioverEzA2[phiSlices], *arrayofDeltaEzA2[phiSlices];
 540   TMatrixD *arrayofEroverEzC2[phiSlices], *arrayofEphioverEzC2[phiSlices], *arrayofDeltaEzC2[phiSlices];
 541
 542
 543   for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
 544
 545     arrayofErA2[k]   =   new TMatrixD(rows,columns) ;
 546     arrayofEphiA2[k] =   new TMatrixD(rows,columns) ;
 547     arrayofdEzA2[k]  =   new TMatrixD(rows,columns) ;
 548     arrayofErC2[k]   =   new TMatrixD(rows,columns) ;
 549     arrayofEphiC2[k] =   new TMatrixD(rows,columns) ;
 550     arrayofdEzC2[k]  =   new TMatrixD(rows,columns) ;
 551
 552     arrayofEroverEzA2[k]   =   new TMatrixD(rows,columns) ;
 553     arrayofEphioverEzA2[k] =   new TMatrixD(rows,columns) ;
 554     arrayofDeltaEzA2[k]    =   new TMatrixD(rows,columns) ;
 555     arrayofEroverEzC2[k]   =   new TMatrixD(rows,columns) ;
 556     arrayofEphioverEzC2[k] =   new TMatrixD(rows,columns) ;
 557     arrayofDeltaEzC2[k]    =   new TMatrixD(rows,columns) ;
 558
 559     // zero initialization not necessary, it is done in the constructor of TMatrix
 560
 561   }
 562
 563
 564   treePOC = (TTree*)fRPhi->Get("POCall");
 565
 566   //  TVector *bEr  = 0;   // done above
 567   TVector *bEphi= 0;
 568   //  TVector *bEz  = 0;   // done above
 569
 570   treePOC->SetBranchAddress("Er",&bEr);
 571   treePOC->SetBranchAddress("Ephi",&bEphi);
 572   treePOC->SetBranchAddress("Ez",&bEz);
 573
 574   // Read the complete tree and do a weighted sum-up over the POC configurations
 575   // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 576
 577   treeNumPOC = (Int_t)treePOC->GetEntries(); // Number of POC conf. in the look-up table
 578   ipC = 0; // POC Conf. counter (note: different to the POC number in the tree!)
 579
 580   for (Int_t itreepC=0; itreepC<treeNumPOC; itreepC++) { // ------------- loop over POC configurations in tree
 581
 582     treePOC->GetEntry(itreepC);
 583
 584     // center of the POC voxel in [meter]
 585     Double_t r0 = coordPOC2(ipC,0);
 586     Double_t phi0 = coordPOC2(ipC,1);
 587     //    Double_t z0 = coordPOC2(ipC,2);
 588
 589      // weights (charge density) at POC position on the A and C side (in C/m^3/e0)
 590     // note: coordinates are in [cm]
 591     Double_t weightA = GetSpaceChargeDensity(r0*100,phi0, 0.499, 2);  // partial load in r,phi
 592     Double_t weightC = GetSpaceChargeDensity(r0*100,phi0,-0.499, 2);  // partial load in r,phi
 593
 594     //    printf("-----\n%lf %lf : %e %e\n",r0,phi0,weightA,weightC);
 595
 596     // Summing up the vector components according to their weight
 597
 598     Int_t ip = 0;
 599     for ( Int_t j = 0 ; j < columns ; j++ ) {
 600       for ( Int_t i = 0 ; i < rows ; i++ )    {
 601         for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
 602
 603           // check wether the coordinates were screwed
 604           if (TMath::Abs((coord2(0,ip)*100-rlist2[i]))>1 ||
 605               TMath::Abs((coord2(1,ip)-philist2[k]))>1 ||
 606               TMath::Abs((coord2(2,ip)*100-zedlist2[j]))>1) {
 607             AliError("internal error: coordinate system was screwed during the sum-up");
 608             printf("lookup: (r,phi,z)=(%lf,%lf,%lf)\n",coord2(0,ip)*100,coord2(1,ip),coord2(2,ip)*100);
 609             printf("sum-up: (r,phi,z)=(%lf,%lf,%lf)\n",rlist2[i],philist2[k],zedlist2[j]);
 610             AliError("Don't trust the results of the space charge calculation!");
 611           }
 612
 613           // unfortunately, the lookup tables were produced to be faster for phi symmetric charges
 614           // This will be the most frequent usage (hopefully)
 615           // That's why we have to do this here ...
 616
 617           TMatrixD &erA   =  *arrayofErA2[k]  ;
 618           TMatrixD &ephiA =  *arrayofEphiA2[k];
 619           TMatrixD &dEzA  =  *arrayofdEzA2[k] ;
 620
 621           TMatrixD &erC   =  *arrayofErC2[k]  ;
 622           TMatrixD &ephiC =  *arrayofEphiC2[k];
 623           TMatrixD &dEzC  =  *arrayofdEzC2[k]   ;
 624
 625           // Sum up - Efield values in [V/m] -> transition to [V/cm]
 626           erA(i,j) += ((*bEr)(ip)) * weightA /100;
 627           erC(i,j) += ((*bEr)(ip)) * weightC /100;
 628           ephiA(i,j) += ((*bEphi)(ip)) * weightA/100; // [V/rad]
 629           ephiC(i,j) += ((*bEphi)(ip)) * weightC/100; // [V/rad]
 630           dEzA(i,j)  += ((*bEz)(ip)) * weightA /100;
 631           dEzC(i,j)  += ((*bEz)(ip)) * weightC /100;
 632
 633           // increase the counter
 634           ip++;
 635         }
 636       }
 637     } // end coordinate loop
 638
 639
 640     // Rotation and summation in the rest of the dPhiSteps
 641     // which were not stored in the this tree due to storage & symmetry reasons
 642
 643
 644     Int_t phiPoints = (Int_t) grid2(1);
 645     Int_t phiPOC    = (Int_t) grid2(4);
 646
 647     //   printf("%d %d\n",phiPOC,flagRadSym);
 648
 649     for (Int_t phiiC = 1; phiiC<phiPOC; phiiC++) { // just used for non-radial symetric table
 650
 651       Double_t phi0R = phiiC*phi0*2 + phi0; // rotate further
 652
 653       // weights (charge density) at POC position on the A and C side (in C/m^3/e0)
 654       // note: coordinates are in [cm] // ecxept z
 655       weightA = GetSpaceChargeDensity(r0*100,phi0R, 0.499, 2);  // partial load in r,phi
 656       weightC = GetSpaceChargeDensity(r0*100,phi0R,-0.499, 2);  // partial load in r,phi
 657
 658       // printf("%lf %lf : %e %e\n",r0,phi0R,weightA,weightC);
 659
 660       // Summing up the vector components according to their weight
 661       ip = 0;
 662       for ( Int_t j = 0 ; j < columns ; j++ ) {
 663         for ( Int_t i = 0 ; i < rows ; i++ )    {
 664           for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
 665
 666             // Note: rotating the coordinated during the sum up
 667
 668             Int_t rotVal = (phiPoints/phiPOC)*phiiC;
 669             Int_t ipR = -1;
 670
 671             if ((ip%phiPoints)>=rotVal) {
 672               ipR = ip-rotVal;
 673             } else {
 674               ipR = ip+(phiPoints-rotVal);
 675             }
 676
 677             // unfortunately, the lookup tables were produced to be faster for phi symmetric charges
 678             // This will be the most frequent usage
 679             // That's why we have to do this here and not outside the loop ...
 680
 681             TMatrixD &erA   =  *arrayofErA2[k]  ;
 682             TMatrixD &ephiA =  *arrayofEphiA2[k];
 683             TMatrixD &dEzA  =  *arrayofdEzA2[k]   ;
 684
 685             TMatrixD &erC   =  *arrayofErC2[k]  ;
 686             TMatrixD &ephiC =  *arrayofEphiC2[k];
 687             TMatrixD &dEzC  =  *arrayofdEzC2[k]   ;
 688
 689             // Sum up - Efield values in [V/m] -> transition to [V/cm]
 690             erA(i,j) += ((*bEr)(ipR)) * weightA /100;
 691             erC(i,j) += ((*bEr)(ipR)) * weightC /100;
 692             ephiA(i,j) += ((*bEphi)(ipR)) * weightA/100; // [V/rad]
 693             ephiC(i,j) += ((*bEphi)(ipR)) * weightC/100; // [V/rad]
 694             dEzA(i,j)  += ((*bEz)(ipR)) * weightA /100;
 695             dEzC(i,j)  += ((*bEz)(ipR)) * weightC /100;
 696
 697             // increase the counter
 698             ip++;
 699           }
 700         }
 701       } // end coordinate loop
 702
 703     } // end phi-POC summation (phiiC)
 704
 705     ipC++; // POC configuration counter
 706
 707     //   printf("POC: (r,phi,z) = (%lf %lf %lf) | weight(A,C): %03.1lf %03.1lf\n",r0,phi0,z0, weightA, weightC);
 708
 709   }
 710
 711
 712
 713
 714   // -------------------------------------------------------------------------------
 715   // Division by the Ez (drift) field and integration along z
 716
 717   //  AliInfo("Step 2: Division and integration");
 718
 719
 720   for ( Int_t k = 0 ; k < phiSlices ; k++ ) { // phi loop
 721
 722     // matrices holding the solution - summation of POC charges // see above
 723     TMatrixD &erA   =  *arrayofErA2[k]  ;
 724     TMatrixD &ephiA =  *arrayofEphiA2[k];
 725     TMatrixD &dezA  =  *arrayofdEzA2[k]   ;
 726     TMatrixD &erC   =  *arrayofErC2[k]  ;
 727     TMatrixD &ephiC =  *arrayofEphiC2[k];
 728     TMatrixD &dezC  =  *arrayofdEzC2[k]   ;
 729
 730     // matrices which will contain the integrated fields (divided by the drift field)
 731     TMatrixD &erOverEzA   =  *arrayofEroverEzA2[k]  ;
 732     TMatrixD &ephiOverEzA =  *arrayofEphioverEzA2[k];
 733     TMatrixD &deltaEzA    =  *arrayofDeltaEzA2[k];
 734     TMatrixD &erOverEzC   =  *arrayofEroverEzC2[k]  ;
 735     TMatrixD &ephiOverEzC =  *arrayofEphioverEzC2[k];
 736     TMatrixD &deltaEzC    =  *arrayofDeltaEzC2[k];
 737
 738     for ( Int_t i = 0 ; i < rows ; i++ )    { // r loop
 739       for ( Int_t j = columns-1 ; j >= 0 ; j-- ) {// z loop
 740         // Count backwards to facilitate integration over Z
 741
 742         Int_t index = 1 ; // Simpsons rule if N=odd.If N!=odd then add extra point by trapezoidal rule.
 743
 744         erOverEzA(i,j) = 0;
 745         ephiOverEzA(i,j) = 0;
 746         deltaEzA(i,j) = 0;
 747         erOverEzC(i,j) = 0;
 748         ephiOverEzC(i,j) = 0;
 749         deltaEzC(i,j) = 0;
 750
 751         for ( Int_t m = j ; m < columns ; m++ ) { // integration
 752
 753           erOverEzA(i,j)   += index*(gridSizeZ/3.0)*erA(i,m)/(-1*ezField) ;
 754           erOverEzC(i,j)   += index*(gridSizeZ/3.0)*erC(i,m)/(-1*ezField)  ;
 755           ephiOverEzA(i,j) += index*(gridSizeZ/3.0)*ephiA(i,m)/(-1*ezField)  ;
 756           ephiOverEzC(i,j) += index*(gridSizeZ/3.0)*ephiC(i,m)/(-1*ezField)  ;
 757           deltaEzA(i,j)    += index*(gridSizeZ/3.0)*dezA(i,m)/(-1) ;
 758           deltaEzC(i,j)    += index*(gridSizeZ/3.0)*dezC(i,m)/(-1) ;
 759
 760           if ( index != 4 )  index = 4; else index = 2 ;
 761
 762         }
 763
 764         if ( index == 4 ) {
 765           erOverEzA(i,j)   -= (gridSizeZ/3.0)*erA(i,columns-1)/(-1*ezField) ;
 766           erOverEzC(i,j)   -= (gridSizeZ/3.0)*erC(i,columns-1)/(-1*ezField) ;
 767           ephiOverEzA(i,j) -= (gridSizeZ/3.0)*ephiA(i,columns-1)/(-1*ezField) ;
 768           ephiOverEzC(i,j) -= (gridSizeZ/3.0)*ephiC(i,columns-1)/(-1*ezField) ;
 769           deltaEzA(i,j)    -= (gridSizeZ/3.0)*dezA(i,columns-1)/(-1) ;
 770           deltaEzC(i,j)    -= (gridSizeZ/3.0)*dezC(i,columns-1)/(-1) ;
 771         }
 772         if ( index == 2 ) {
 773           erOverEzA(i,j)   += (gridSizeZ/3.0)*(0.5*erA(i,columns-2)-2.5*erA(i,columns-1))/(-1*ezField) ;
 774           erOverEzC(i,j)   += (gridSizeZ/3.0)*(0.5*erC(i,columns-2)-2.5*erC(i,columns-1))/(-1*ezField) ;
 775           ephiOverEzA(i,j) += (gridSizeZ/3.0)*(0.5*ephiA(i,columns-2)-2.5*ephiA(i,columns-1))/(-1*ezField) ;
 776           ephiOverEzC(i,j) += (gridSizeZ/3.0)*(0.5*ephiC(i,columns-2)-2.5*ephiC(i,columns-1))/(-1*ezField) ;
 777           deltaEzA(i,j)    += (gridSizeZ/3.0)*(0.5*dezA(i,columns-2)-2.5*dezA(i,columns-1))/(-1) ;
 778           deltaEzC(i,j)    += (gridSizeZ/3.0)*(0.5*dezC(i,columns-2)-2.5*dezC(i,columns-1))/(-1) ;
 779         }
 780         if ( j == columns-2 ) {
 781           erOverEzA(i,j)   = (gridSizeZ/3.0)*(1.5*erA(i,columns-2)+1.5*erA(i,columns-1))/(-1*ezField) ;
 782           erOverEzC(i,j)   = (gridSizeZ/3.0)*(1.5*erC(i,columns-2)+1.5*erC(i,columns-1))/(-1*ezField) ;
 783           ephiOverEzA(i,j) = (gridSizeZ/3.0)*(1.5*ephiA(i,columns-2)+1.5*ephiA(i,columns-1))/(-1*ezField) ;
 784           ephiOverEzC(i,j) = (gridSizeZ/3.0)*(1.5*ephiC(i,columns-2)+1.5*ephiC(i,columns-1))/(-1*ezField) ;
 785           deltaEzA(i,j)    = (gridSizeZ/3.0)*(1.5*dezA(i,columns-2)+1.5*dezA(i,columns-1))/(-1) ;
 786           deltaEzC(i,j)    = (gridSizeZ/3.0)*(1.5*dezC(i,columns-2)+1.5*dezC(i,columns-1))/(-1) ;
 787         }
 788         if ( j == columns-1 ) {
 789           erOverEzA(i,j)   = 0;
 790           erOverEzC(i,j)   = 0;
 791           ephiOverEzA(i,j) = 0;
 792           ephiOverEzC(i,j) = 0;
 793           deltaEzA(i,j)    = 0;
 794           deltaEzC(i,j)    = 0;
 795         }
 796       }
 797     }
 798
 799   }
 800
 801   AliInfo("Step 2: Interpolation to Standard grid");
 802
 803   // -------------------------------------------------------------------------------
 804   // Interpolate results onto the standard grid which is used for all AliTPCCorrections classes
 805
 806
 807   for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
 808     Double_t phi = fgkPhiList[k] ;
 809
 810     // final lookup table
 811     TMatrixD &erOverEzFinal   =  *fLookUpErOverEz[k]  ;
 812     TMatrixD &ephiOverEzFinal =  *fLookUpEphiOverEz[k];
 813     TMatrixD &deltaEzFinal    =  *fLookUpDeltaEz[k]   ;
 814
 815     for ( Int_t j = 0 ; j < kNZ ; j++ ) {
 816
 817       z = TMath::Abs(fgkZList[j]) ;  // z position is symmetric
 818
 819       for ( Int_t i = 0 ; i < kNR ; i++ ) {
 820         r = fgkRList[i] ;
 821
 822         // Interpolate Lookup tables onto standard grid
 823         if (fgkZList[j]>0) {
 824           erOverEzFinal(i,j)   += Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
 825                                                rlist2, zedlist2, philist2, arrayofEroverEzA2  );
 826           ephiOverEzFinal(i,j) += Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
 827                                                rlist2, zedlist2, philist2, arrayofEphioverEzA2);
 828           deltaEzFinal(i,j)    += Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
 829                                                rlist2, zedlist2, philist2, arrayofDeltaEzA2   );
 830         } else {
 831           erOverEzFinal(i,j)   += Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
 832                                                rlist2, zedlist2, philist2, arrayofEroverEzC2  );
 833           ephiOverEzFinal(i,j) += Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
 834                                                rlist2, zedlist2, philist2, arrayofEphioverEzC2);
 835           deltaEzFinal(i,j)  -=  Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
 836                                                rlist2, zedlist2, philist2, arrayofDeltaEzC2   );
 837         }
 838
 839       } // end r loop
 840     } // end z loop
 841   } // end phi loop
 842
 843
 844   // clear the temporary arrays lists
 845   for ( Int_t k = 0 ; k < phiSlices ; k++ )  {
 846
 847     delete arrayofErA2[k];
 848     delete arrayofEphiA2[k];
 849     delete arrayofdEzA2[k];
 850     delete arrayofErC2[k];
 851     delete arrayofEphiC2[k];
 852     delete arrayofdEzC2[k];
 853
 854     delete arrayofEroverEzA2[k];
 855     delete arrayofEphioverEzA2[k];
 856     delete arrayofDeltaEzA2[k];
 857     delete arrayofEroverEzC2[k];
 858     delete arrayofEphioverEzC2[k];
 859     delete arrayofDeltaEzC2[k];
 860
 861   }
 862
 863   fRPhi->Close();
 864
 865   // FINISHED
 866
 867   fInitLookUp = kTRUE;
 868
 869 }
 870
 871 void AliTPCSpaceCharge3D::InitSpaceCharge3DDistortionCourse() {
 872   //
 873   // Initialization of the Lookup table which contains the solutions of the
 874   // "space charge" (poisson) problem
 875   //
 876   // The sum-up uses a look-up table which contains different discretized Space charge fields
 877   // in order to calculate the corresponding field deviations due to a given (discretized)
 878   // space charge distribution ....
 879   //
 880   // Method of calculation: Weighted sum-up of the different fields within the look up table
 881   // Note: Full 3d version: Course and slow ...
 882
 883   if (fInitLookUp) {
 884     AliInfo("Lookup table was already initialized!");
 885     //    return;
 886   }
 887
 888   AliInfo("Preparation of the weighted look-up table");
 889
 890   TFile *f = new TFile(fSCLookUpPOCsFileName3D.Data(),"READ");
 891   if ( !f ) {
 892     AliError("Precalculated POC-looup-table could not be found");
 893     return;
 894   }
 895
 896   // units are in [m]
 897   TVector *gridf  = (TVector*) f->Get("constants");
 898   TVector &grid = *gridf;
 899   TMatrix *coordf  = (TMatrix*) f->Get("coordinates");
 900   TMatrix &coord = *coordf;
 901   TMatrix *coordPOCf  = (TMatrix*) f->Get("POCcoord");
 902   TMatrix &coordPOC = *coordPOCf;
 903
 904   Bool_t flagRadSym = 0;
 905   if (grid(1)==1 && grid(4)==1) {
 906     //   AliInfo("LOOK UP TABLE IS RADIAL SYMETTRIC - Field in Phi is ZERO");
 907     flagRadSym=1;
 908   }
 909
 910   Int_t rows      = (Int_t)grid(0);   // number of points in r direction
 911   Int_t phiSlices = (Int_t)grid(1);   // number of points in phi
 912   Int_t columns   = (Int_t)grid(2);   // number of points in z direction
 913
 914   const Float_t gridSizeR   =  (fgkOFCRadius-fgkIFCRadius)/(rows-1);   // unit in [cm]
 915   const Float_t gridSizePhi =  TMath::TwoPi()/phiSlices;         // unit in [rad]
 916   const Float_t gridSizeZ   =  fgkTPCZ0/(columns-1);                  // unit in [cm]
 917
 918   // temporary matrices needed for the calculation
 919   TMatrixD *arrayofErA[phiSlices], *arrayofEphiA[phiSlices], *arrayofdEzA[phiSlices];
 920   TMatrixD *arrayofErC[phiSlices], *arrayofEphiC[phiSlices], *arrayofdEzC[phiSlices];
 921
 922   TMatrixD *arrayofEroverEzA[phiSlices], *arrayofEphioverEzA[phiSlices], *arrayofDeltaEzA[phiSlices];
 923   TMatrixD *arrayofEroverEzC[phiSlices], *arrayofEphioverEzC[phiSlices], *arrayofDeltaEzC[phiSlices];
 924
 925
 926   for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
 927
 928     arrayofErA[k]   =   new TMatrixD(rows,columns) ;
 929     arrayofEphiA[k] =   new TMatrixD(rows,columns) ; // zero if radial symmetric
 930     arrayofdEzA[k]  =   new TMatrixD(rows,columns) ;
 931     arrayofErC[k]   =   new TMatrixD(rows,columns) ;
 932     arrayofEphiC[k] =   new TMatrixD(rows,columns) ; // zero if radial symmetric
 933     arrayofdEzC[k]  =   new TMatrixD(rows,columns) ;
 934
 935     arrayofEroverEzA[k]   =   new TMatrixD(rows,columns) ;
 936     arrayofEphioverEzA[k] =   new TMatrixD(rows,columns) ; // zero if radial symmetric
 937     arrayofDeltaEzA[k]    =   new TMatrixD(rows,columns) ;
 938     arrayofEroverEzC[k]   =   new TMatrixD(rows,columns) ;
 939     arrayofEphioverEzC[k] =   new TMatrixD(rows,columns) ; // zero if radial symmetric
 940     arrayofDeltaEzC[k]    =   new TMatrixD(rows,columns) ;
 941
 942     // Set the values to zero the lookup tables
 943     // not necessary, it is done in the constructor of TMatrix - code deleted
 944
 945   }
 946
 947   // list of points as used in the interpolation (during sum up)
 948   Double_t  rlist[rows], zedlist[columns] , philist[phiSlices];
 949   for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
 950     philist[k] =  gridSizePhi * k;
 951     for ( Int_t i = 0 ; i < rows ; i++ )    {
 952       rlist[i] = fgkIFCRadius + i*gridSizeR ;
 953       for ( Int_t j = 0 ; j < columns ; j++ ) {
 954         zedlist[j]  = j * gridSizeZ ;
 955       }
 956     }
 957   } // only done once
 958
 959
 960   TTree *treePOC = (TTree*)f->Get("POCall");
 961
 962   TVector *bEr  = 0;   TVector *bEphi= 0;   TVector *bEz  = 0;
 963
 964   treePOC->SetBranchAddress("Er",&bEr);
 965   if (!flagRadSym) treePOC->SetBranchAddress("Ephi",&bEphi);
 966   treePOC->SetBranchAddress("Ez",&bEz);
 967
 968
 969   // Read the complete tree and do a weighted sum-up over the POC configurations
 970   // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 971
 972   Int_t treeNumPOC = (Int_t)treePOC->GetEntries(); // Number of POC conf. in the look-up table
 973   Int_t ipC = 0; // POC Conf. counter (note: different to the POC number in the tree!)
 974
 975   for (Int_t itreepC=0; itreepC<treeNumPOC; itreepC++) { // ------------- loop over POC configurations in tree
 976
 977     treePOC->GetEntry(itreepC);
 978
 979     // center of the POC voxel in [meter]
 980     Double_t r0 = coordPOC(ipC,0);
 981     Double_t phi0 = coordPOC(ipC,1);
 982     Double_t z0 = coordPOC(ipC,2);
 983
 984     ipC++; // POC configuration counter
 985
 986     // weights (charge density) at POC position on the A and C side (in C/m^3/e0)
 987     // note: coordinates are in [cm]
 988     Double_t weightA = GetSpaceChargeDensity(r0*100,phi0, z0*100);
 989     Double_t weightC = GetSpaceChargeDensity(r0*100,phi0,-z0*100);
 990
 991     // Summing up the vector components according to their weight
 992
 993     Int_t ip = 0;
 994     for ( Int_t j = 0 ; j < columns ; j++ ) {
 995       for ( Int_t i = 0 ; i < rows ; i++ )    {
 996         for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
 997
 998           // check wether the coordinates were screwed
 999           if (TMath::Abs((coord(0,ip)*100-rlist[i]))>1 ||
1000               TMath::Abs((coord(1,ip)-philist[k]))>1 ||
1001               TMath::Abs((coord(2,ip)*100-zedlist[j]))>1) {
1002             AliError("internal error: coordinate system was screwed during the sum-up");
1003             printf("lookup: (r,phi,z)=(%lf,%lf,%lf)\n",coord(0,ip)*100,coord(1,ip),coord(2,ip)*100);
1004             printf("sum-up: (r,phi,z)=(%lf,%lf,%lf)\n",rlist[i],philist[k],zedlist[j]);
1005             AliError("Don't trust the results of the space charge calculation!");
1006           }
1007
1008           // unfortunately, the lookup tables were produced to be faster for phi symmetric charges
1009           // This will be the most frequent usage (hopefully)
1010           // That's why we have to do this here ...
1011
1012           TMatrixD &erA   =  *arrayofErA[k]  ;
1013           TMatrixD &ephiA =  *arrayofEphiA[k];
1014           TMatrixD &dEzA  =  *arrayofdEzA[k]   ;
1015
1016           TMatrixD &erC   =  *arrayofErC[k]  ;
1017           TMatrixD &ephiC =  *arrayofEphiC[k];
1018           TMatrixD &dEzC  =  *arrayofdEzC[k]   ;
1019
1020           // Sum up - Efield values in [V/m] -> transition to [V/cm]
1021           erA(i,j) += ((*bEr)(ip)) * weightA /100;
1022           erC(i,j) += ((*bEr)(ip)) * weightC /100;
1023           if (!flagRadSym) {
1024             ephiA(i,j) += ((*bEphi)(ip)) * weightA/100; // [V/rad]
1025             ephiC(i,j) += ((*bEphi)(ip)) * weightC/100; // [V/rad]
1026           }
1027           dEzA(i,j)  += ((*bEz)(ip)) * weightA /100;
1028           dEzC(i,j)  += ((*bEz)(ip)) * weightC /100;
1029
1030           // increase the counter
1031           ip++;
1032         }
1033       }
1034     } // end coordinate loop
1035
1036
1037     // Rotation and summation in the rest of the dPhiSteps
1038     // which were not stored in the this tree due to storage & symmetry reasons
1039
1040     Int_t phiPoints = (Int_t) grid(1);
1041     Int_t phiPOC    = (Int_t) grid(4);
1042
1043     //   printf("%d %d\n",phiPOC,flagRadSym);
1044
1045     for (Int_t phiiC = 1; phiiC<phiPOC; phiiC++) { // just used for non-radial symetric table
1046
1047       r0 = coordPOC(ipC,0);
1048       phi0 = coordPOC(ipC,1);
1049       z0 = coordPOC(ipC,2);
1050
1051       ipC++; // POC conf. counter
1052
1053       // weights (charge density) at POC position on the A and C side (in C/m^3/e0)
1054       // note: coordinates are in [cm]
1055       weightA = GetSpaceChargeDensity(r0*100,phi0, z0*100);
1056       weightC = GetSpaceChargeDensity(r0*100,phi0,-z0*100);
1057
1058       //     printf("%lf %lf %lf: %e %e\n",r0,phi0,z0,weightA,weightC);
1059
1060       // Summing up the vector components according to their weight
1061       ip = 0;
1062       for ( Int_t j = 0 ; j < columns ; j++ ) {
1063         for ( Int_t i = 0 ; i < rows ; i++ )    {
1064           for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
1065
1066             // Note: rotating the coordinated during the sum up
1067
1068             Int_t rotVal = (phiPoints/phiPOC)*phiiC;
1069             Int_t ipR = -1;
1070
1071             if ((ip%phiPoints)>=rotVal) {
1072               ipR = ip-rotVal;
1073             } else {
1074               ipR = ip+(phiPoints-rotVal);
1075             }
1076
1077             // unfortunately, the lookup tables were produced to be faster for phi symmetric charges
1078             // This will be the most frequent usage
1079             // That's why we have to do this here and not outside the loop ...
1080
1081             TMatrixD &erA   =  *arrayofErA[k]  ;
1082             TMatrixD &ephiA =  *arrayofEphiA[k];
1083             TMatrixD &dEzA  =  *arrayofdEzA[k]   ;
1084
1085             TMatrixD &erC   =  *arrayofErC[k]  ;
1086             TMatrixD &ephiC =  *arrayofEphiC[k];
1087             TMatrixD &dEzC  =  *arrayofdEzC[k]   ;
1088
1089             // Sum up - Efield values in [V/m] -> transition to [V/cm]
1090             erA(i,j) += ((*bEr)(ipR)) * weightA /100;
1091             erC(i,j) += ((*bEr)(ipR)) * weightC /100;
1092             if (!flagRadSym) {
1093               ephiA(i,j) += ((*bEphi)(ipR)) * weightA/100; // [V/rad]
1094               ephiC(i,j) += ((*bEphi)(ipR)) * weightC/100; // [V/rad]
1095             }
1096             dEzA(i,j)  += ((*bEz)(ipR)) * weightA /100;
1097             dEzC(i,j)  += ((*bEz)(ipR)) * weightC /100;
1098
1099             // increase the counter
1100             ip++;
1101           }
1102         }
1103       } // end coordinate loop
1104
1105     } // end phi-POC summation (phiiC)
1106
1107
1108     // printf("POC: (r,phi,z) = (%lf %lf %lf) | weight(A,C): %03.1lf %03.1lf\n",r0,phi0,z0, weightA, weightC);
1109
1110   }
1111
1112
1113
1114   // -------------------------------------------------------------------------------
1115   // Division by the Ez (drift) field and integration along z
1116
1117   AliInfo("Division and integration");
1118
1119   Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = Electric Field (V/cm) Magnitude ~ -400 V/cm;
1120
1121   for ( Int_t k = 0 ; k < phiSlices ; k++ ) { // phi loop
1122
1123     // matrices holding the solution - summation of POC charges // see above
1124     TMatrixD &erA   =  *arrayofErA[k]  ;
1125     TMatrixD &ephiA =  *arrayofEphiA[k];
1126     TMatrixD &dezA  =  *arrayofdEzA[k]   ;
1127     TMatrixD &erC   =  *arrayofErC[k]  ;
1128     TMatrixD &ephiC =  *arrayofEphiC[k];
1129     TMatrixD &dezC  =  *arrayofdEzC[k]   ;
1130
1131     // matrices which will contain the integrated fields (divided by the drift field)
1132     TMatrixD &erOverEzA   =  *arrayofEroverEzA[k]  ;
1133     TMatrixD &ephiOverEzA =  *arrayofEphioverEzA[k];
1134     TMatrixD &deltaEzA    =  *arrayofDeltaEzA[k];
1135     TMatrixD &erOverEzC   =  *arrayofEroverEzC[k]  ;
1136     TMatrixD &ephiOverEzC =  *arrayofEphioverEzC[k];
1137     TMatrixD &deltaEzC    =  *arrayofDeltaEzC[k];
1138
1139     for ( Int_t i = 0 ; i < rows ; i++ )    { // r loop
1140       for ( Int_t j = columns-1 ; j >= 0 ; j-- ) {// z loop
1141         // Count backwards to facilitate integration over Z
1142
1143         Int_t index = 1 ; // Simpsons rule if N=odd.If N!=odd then add extra point by trapezoidal rule.
1144
1145         erOverEzA(i,j) = 0; ephiOverEzA(i,j) = 0;  deltaEzA(i,j) = 0;
1146         erOverEzC(i,j) = 0; ephiOverEzC(i,j) = 0;  deltaEzC(i,j) = 0;
1147
1148         for ( Int_t m = j ; m < columns ; m++ ) { // integration
1149
1150           erOverEzA(i,j)   += index*(gridSizeZ/3.0)*erA(i,m)/(-1*ezField) ;
1151           erOverEzC(i,j)   += index*(gridSizeZ/3.0)*erC(i,m)/(-1*ezField)  ;
1152           if (!flagRadSym) {
1153             ephiOverEzA(i,j) += index*(gridSizeZ/3.0)*ephiA(i,m)/(-1*ezField)  ;
1154             ephiOverEzC(i,j) += index*(gridSizeZ/3.0)*ephiC(i,m)/(-1*ezField)  ;
1155           }
1156           deltaEzA(i,j)    += index*(gridSizeZ/3.0)*dezA(i,m)/(-1) ;
1157           deltaEzC(i,j)    += index*(gridSizeZ/3.0)*dezC(i,m)/(-1) ;
1158
1159           if ( index != 4 )  index = 4; else index = 2 ;
1160
1161         }
1162
1163         if ( index == 4 ) {
1164           erOverEzA(i,j)   -= (gridSizeZ/3.0)*erA(i,columns-1)/(-1*ezField) ;
1165           erOverEzC(i,j)   -= (gridSizeZ/3.0)*erC(i,columns-1)/(-1*ezField) ;
1166           if (!flagRadSym) {
1167             ephiOverEzA(i,j) -= (gridSizeZ/3.0)*ephiA(i,columns-1)/(-1*ezField) ;
1168             ephiOverEzC(i,j) -= (gridSizeZ/3.0)*ephiC(i,columns-1)/(-1*ezField) ;
1169           }
1170           deltaEzA(i,j)    -= (gridSizeZ/3.0)*dezA(i,columns-1)/(-1) ;
1171           deltaEzC(i,j)    -= (gridSizeZ/3.0)*dezC(i,columns-1)/(-1) ;
1172         }
1173         if ( index == 2 ) {
1174           erOverEzA(i,j)   += (gridSizeZ/3.0)*(0.5*erA(i,columns-2)-2.5*erA(i,columns-1))/(-1*ezField) ;
1175           erOverEzC(i,j)   += (gridSizeZ/3.0)*(0.5*erC(i,columns-2)-2.5*erC(i,columns-1))/(-1*ezField) ;
1176           if (!flagRadSym) {
1177             ephiOverEzA(i,j) += (gridSizeZ/3.0)*(0.5*ephiA(i,columns-2)-2.5*ephiA(i,columns-1))/(-1*ezField) ;
1178             ephiOverEzC(i,j) += (gridSizeZ/3.0)*(0.5*ephiC(i,columns-2)-2.5*ephiC(i,columns-1))/(-1*ezField) ;
1179           }
1180           deltaEzA(i,j)    += (gridSizeZ/3.0)*(0.5*dezA(i,columns-2)-2.5*dezA(i,columns-1))/(-1) ;
1181           deltaEzC(i,j)    += (gridSizeZ/3.0)*(0.5*dezC(i,columns-2)-2.5*dezC(i,columns-1))/(-1) ;
1182         }
1183         if ( j == columns-2 ) {
1184           erOverEzA(i,j)   = (gridSizeZ/3.0)*(1.5*erA(i,columns-2)+1.5*erA(i,columns-1))/(-1*ezField) ;
1185           erOverEzC(i,j)   = (gridSizeZ/3.0)*(1.5*erC(i,columns-2)+1.5*erC(i,columns-1))/(-1*ezField) ;
1186           if (!flagRadSym) {
1187             ephiOverEzA(i,j) = (gridSizeZ/3.0)*(1.5*ephiA(i,columns-2)+1.5*ephiA(i,columns-1))/(-1*ezField) ;
1188             ephiOverEzC(i,j) = (gridSizeZ/3.0)*(1.5*ephiC(i,columns-2)+1.5*ephiC(i,columns-1))/(-1*ezField) ;
1189           }
1190           deltaEzA(i,j)    = (gridSizeZ/3.0)*(1.5*dezA(i,columns-2)+1.5*dezA(i,columns-1))/(-1) ;
1191           deltaEzC(i,j)    = (gridSizeZ/3.0)*(1.5*dezC(i,columns-2)+1.5*dezC(i,columns-1))/(-1) ;
1192         }
1193         if ( j == columns-1 ) {
1194           erOverEzA(i,j)   = 0;
1195           erOverEzC(i,j)   = 0;
1196           if (!flagRadSym) {
1197             ephiOverEzA(i,j) = 0;
1198             ephiOverEzC(i,j) = 0;
1199           }
1200           deltaEzA(i,j)    = 0;
1201           deltaEzC(i,j)    = 0;
1202         }
1203       }
1204     }
1205
1206   }
1207
1208
1209
1210   AliInfo("Interpolation to Standard grid");
1211
1212   // -------------------------------------------------------------------------------
1213   // Interpolate results onto the standard grid which is used for all AliTPCCorrections classes
1214
1215   const Int_t order  = 1  ;  // Linear interpolation = 1, Quadratic = 2
1216
1217   Double_t  r, phi, z ;
1218   for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
1219     phi = fgkPhiList[k] ;
1220
1221     TMatrixD &erOverEz   =  *fLookUpErOverEz[k]  ;
1222     TMatrixD &ephiOverEz =  *fLookUpEphiOverEz[k];
1223     TMatrixD &deltaEz    =  *fLookUpDeltaEz[k]   ;
1224
1225     for ( Int_t j = 0 ; j < kNZ ; j++ ) {
1226
1227       z = TMath::Abs(fgkZList[j]) ;  // z position is symmetric
1228
1229       for ( Int_t i = 0 ; i < kNR ; i++ ) {
1230         r = fgkRList[i] ;
1231
1232         // Interpolate Lookup tables onto standard grid
1233         if (fgkZList[j]>0) {
1234           erOverEz(i,j)   = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
1235                                                rlist, zedlist, philist, arrayofEroverEzA  );
1236           ephiOverEz(i,j) = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
1237                                                rlist, zedlist, philist, arrayofEphioverEzA);
1238           deltaEz(i,j)    = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
1239                                                rlist, zedlist, philist, arrayofDeltaEzA   );
1240         } else {
1241           erOverEz(i,j)   = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
1242                                                rlist, zedlist, philist, arrayofEroverEzC  );
1243           ephiOverEz(i,j) = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
1244                                                rlist, zedlist, philist, arrayofEphioverEzC);
1245           deltaEz(i,j)  = - Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
1246                                                rlist, zedlist, philist, arrayofDeltaEzC   );
1247           // negative coordinate system on C side
1248         }
1249
1250       } // end r loop
1251     } // end z loop
1252   } // end phi loop
1253
1254
1255   // clear the temporary arrays lists
1256   for ( Int_t k = 0 ; k < phiSlices ; k++ )  {
1257
1258     delete arrayofErA[k];
1259     delete arrayofEphiA[k];
1260     delete arrayofdEzA[k];
1261     delete arrayofErC[k];
1262     delete arrayofEphiC[k];
1263     delete arrayofdEzC[k];
1264
1265     delete arrayofEroverEzA[k];
1266     delete arrayofEphioverEzA[k];
1267     delete arrayofDeltaEzA[k];
1268     delete arrayofEroverEzC[k];
1269     delete arrayofEphioverEzC[k];
1270     delete arrayofDeltaEzC[k];
1271
1272   }
1273
1274   fInitLookUp = kTRUE;
1275
1276 }
1277
1278
1279 void AliTPCSpaceCharge3D::SetSCDataFileName(TString fname) {
1280   //
1281   // Set & load the Space charge density distribution from a file
1282   // (linear interpolation onto a standard grid)
1283   //
1284
1285
1286   fSCDataFileName = fname;
1287
1288   TFile *f = new TFile(fSCDataFileName.Data(),"READ");
1289   if (!f) {
1290     AliError(Form("File %s, which should contain the space charge distribution, could not be found",
1291                   fSCDataFileName.Data()));
1292     return;
1293   }
1294
1295   TH2F *densityRZ = (TH2F*) f->Get("SpaceChargeInRZ");
1296   if (!densityRZ) {
1297     AliError(Form("The indicated file (%s) does not contain a histogram called %s",
1298                   fSCDataFileName.Data(),"SpaceChargeInRZ"));
1299     return;
1300   }
1301
1302   TH3F *densityRPhi = (TH3F*) f->Get("SpaceChargeInRPhi");
1303   if (!densityRPhi) {
1304     AliError(Form("The indicated file (%s) does not contain a histogram called %s",
1305                   fSCDataFileName.Data(),"SpaceChargeInRPhi"));
1306     return;
1307   }
1308
1309
1310   Double_t  r, phi, z ;
1311
1312   TMatrixD &scDensityInRZ   =  *fSCdensityInRZ;
1313   TMatrixD &scDensityInRPhiA   =  *fSCdensityInRPhiA;
1314   TMatrixD &scDensityInRPhiC   =  *fSCdensityInRPhiC;
1315   for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
1316     phi = fgkPhiList[k] ;
1317     TMatrixD &scDensity   =  *fSCdensityDistribution[k]  ;
1318     for ( Int_t j = 0 ; j < kNZ ; j++ ) {
1319       z = fgkZList[j] ;
1320       for ( Int_t i = 0 ; i < kNR ; i++ ) {
1321         r = fgkRList[i] ;
1322
1323         // partial load in (r,z)
1324         if (k==0) // do just once
1325           scDensityInRZ(i,j) =  densityRZ->Interpolate(r,z);
1326
1327         // partial load in (r,phi)
1328         if ( j==0 || j == kNZ/2 ) {
1329           if (z>0)
1330             scDensityInRPhiA(i,k) =  densityRPhi->Interpolate(r,phi,0.499);  // A side
1331           else
1332             scDensityInRPhiC(i,k) =  densityRPhi->Interpolate(r,phi,-0.499); // C side
1333         }
1334
1335         // Full 3D configuration ...
1336         if (z>0)
1337            scDensity(i,j) = scDensityInRZ(i,j) + scDensityInRPhiA(i,k);
1338         else
1339            scDensity(i,j) = scDensityInRZ(i,j) + scDensityInRPhiC(i,k);
1340       }
1341     }
1342   }
1343
1344   f->Close();
1345
1346   fInitLookUp = kFALSE;
1347
1348
1349 }
1350
1351
1352 Float_t  AliTPCSpaceCharge3D::GetSpaceChargeDensity(Float_t r, Float_t phi, Float_t z, Int_t mode) {
1353   //
1354   // returns the (input) space charge density at a given point according
1355   // Note: input in [cm], output in [C/m^3/e0] !!
1356   //
1357
1358   if (!fSCdensityDistribution || !fSCdensityInRZ || !fSCdensityInRPhiA || !fSCdensityInRPhiC ) {
1359     printf("Irgend a schaaas is nuul - argg\n");
1360     return 0.;
1361   }
1362
1363   while (phi<0) phi += TMath::TwoPi();
1364   while (phi>TMath::TwoPi()) phi -= TMath::TwoPi();
1365
1366
1367   // Float_t sc =fSCdensityDistribution->Interpolate(r0,phi0,z0);
1368   Int_t order = 1; //
1369   Float_t sc = 0;
1370
1371   if (mode == 0) { // return full load
1372     sc = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
1373                             fgkRList, fgkZList, fgkPhiList, fSCdensityDistribution );
1374
1375   } else if (mode == 1) { // return partial load in (r,z)
1376     TMatrixD &scDensityInRZ   =  *fSCdensityInRZ;
1377     sc = Interpolate2DTable(order, r, z, kNR, kNZ, fgkRList, fgkZList, scDensityInRZ );
1378
1379   } else if (mode == 2) { // return partial load in (r,phi)
1380
1381     if (z>0) {
1382       TMatrixD &scDensityInRPhi   =  *fSCdensityInRPhiA;
1383       sc = Interpolate2DTable(order, r, phi, kNR, kNPhi, fgkRList, fgkPhiList, scDensityInRPhi );
1384     } else {
1385       TMatrixD &scDensityInRPhi   =  *fSCdensityInRPhiC;
1386       sc = Interpolate2DTable(order, r, phi, kNR, kNPhi, fgkRList, fgkPhiList, scDensityInRPhi );
1387     }
1388
1389   } else {
1390     // should i give a warning?
1391     sc = 0;
1392   }
1393
1394   //  printf("%lf %lf %lf: %lf\n",r,phi,z,sc);
1395
1396   return sc;
1397 }
1398
1399
1400 TH2F * AliTPCSpaceCharge3D::CreateHistoSCinXY(Float_t z, Int_t nx, Int_t ny, Int_t mode) {
1401   //
1402   // return a simple histogramm containing the space charge distribution (input for the calculation)
1403   //
1404
1405   TH2F *h=CreateTH2F("spaceCharge",GetTitle(),"x [cm]","y [cm]","#rho_{sc} [C/m^{3}/e_{0}]",
1406                      nx,-250.,250.,ny,-250.,250.);
1407
1408   for (Int_t iy=1;iy<=ny;++iy) {
1409     Double_t yp = h->GetYaxis()->GetBinCenter(iy);
1410     for (Int_t ix=1;ix<=nx;++ix) {
1411       Double_t xp = h->GetXaxis()->GetBinCenter(ix);
1412
1413       Float_t r = TMath::Sqrt(xp*xp+yp*yp);
1414       Float_t phi = TMath::ATan2(yp,xp);
1415
1416       if (85.<=r && r<=250.) {
1417         Float_t sc = GetSpaceChargeDensity(r,phi,z,mode)/fgke0; // in [C/m^3/e0]
1418         h->SetBinContent(ix,iy,sc);
1419       } else {
1420         h->SetBinContent(ix,iy,0.);
1421       }
1422     }
1423   }
1424
1425   return h;
1426 }
1427
1428 TH2F * AliTPCSpaceCharge3D::CreateHistoSCinZR(Float_t phi, Int_t nz, Int_t nr,Int_t mode ) {
1429   //
1430   // return a simple histogramm containing the space charge distribution (input for the calculation)
1431   //
1432
1433   TH2F *h=CreateTH2F("spaceCharge",GetTitle(),"z [cm]","r [cm]","#rho_{sc} [C/m^{3}/e_{0}]",
1434                      nz,-250.,250.,nr,85.,250.);
1435
1436   for (Int_t ir=1;ir<=nr;++ir) {
1437     Float_t r = h->GetYaxis()->GetBinCenter(ir);
1438     for (Int_t iz=1;iz<=nz;++iz) {
1439       Float_t z = h->GetXaxis()->GetBinCenter(iz);
1440       Float_t sc = GetSpaceChargeDensity(r,phi,z,mode)/fgke0; // in [C/m^3/e0]
1441       h->SetBinContent(iz,ir,sc);
1442     }
1443   }
1444
1445   return h;
1446 }
1447
1448 void AliTPCSpaceCharge3D::WriteChargeDistributionToFile(const char* fname) {
1449   //
1450   // Example on how to write a Space charge distribution into a File
1451   //  (see below: estimate from scaling STAR measurements to Alice)
1452   // Charge distribution is splitted into two (RZ and RPHI) in order to speed up
1453   // the needed calculation time
1454   //
1455
1456   TFile *f = new TFile(fname,"RECREATE");
1457
1458   // some grid, not too course
1459   Int_t nr = 350;
1460   Int_t nphi = 180;
1461   Int_t nz = 500;
1462
1463   Double_t dr = (fgkOFCRadius-fgkIFCRadius)/(nr+1);
1464   Double_t dphi = TMath::TwoPi()/(nphi+1);
1465   Double_t dz = 500./(nz+1);
1466   Double_t safty = 0.; // due to a root bug which does not interpolate the boundary (first and last bin) correctly
1467
1468
1469   // Charge distribution in ZR (rotational symmetric) ------------------
1470
1471   TH2F *histoZR = new TH2F("chargeZR","chargeZR",
1472                            nr,fgkIFCRadius-dr-safty,fgkOFCRadius+dr+safty,
1473                            nz,-250-dz-safty,250+dz+safty);
1474
1475   for (Int_t ir=1;ir<=nr;++ir) {
1476     Double_t rp = histoZR->GetXaxis()->GetBinCenter(ir);
1477     for (Int_t iz=1;iz<=nz;++iz) {
1478       Double_t zp = histoZR->GetYaxis()->GetBinCenter(iz);
1479
1480       // recalculation to meter
1481       Double_t lZ = 2.5; // approx. TPC drift length
1482       Double_t rpM = rp/100.; // in [m]
1483       Double_t zpM = TMath::Abs(zp/100.); // in [m]
1484
1485       // setting of mb multiplicity and Interaction rate
1486       Double_t multiplicity = 950;
1487       Double_t intRate = 7800;
1488
1489       // calculation of "scaled" parameters
1490       Double_t a = multiplicity*intRate/79175;
1491       Double_t b = a/lZ;
1492
1493       Double_t charge = (a - b*zpM)/(rpM*rpM); // charge in [C/m^3/e0]
1494
1495       charge = charge*fgke0; // [C/m^3]
1496
1497       if (zp<0) charge *= 0.9; // e.g. slightly less on C side due to front absorber
1498
1499       //  charge = 0; // for tests
1500       histoZR->SetBinContent(ir,iz,charge);
1501     }
1502   }
1503
1504   histoZR->Write("SpaceChargeInRZ");
1505
1506
1507   // Charge distribution in RPhi (e.g. Floating GG wire) ------------
1508
1509   TH3F *histoRPhi = new TH3F("chargeRPhi","chargeRPhi",
1510                              nr,fgkIFCRadius-dr-safty,fgkOFCRadius+dr+safty,
1511                              nphi,0-dphi-safty,TMath::TwoPi()+dphi+safty,
1512                              2,-1,1); // z part - to allow A and C side differences
1513
1514   // some 'arbitrary' GG leaks
1515   Int_t   nGGleaks = 5;
1516   Double_t secPosA[5]    = {3,6,6,11,13};         // sector
1517   Double_t radialPosA[5] = {125,100,160,200,230}; // radius in cm
1518   Double_t secPosC[5]    = {1,8,12,15,15};        // sector
1519   Double_t radialPosC[5] = {245,120,140,120,190}; // radius in cm
1520
1521   for (Int_t ir=1;ir<=nr;++ir) {
1522     Double_t rp = histoRPhi->GetXaxis()->GetBinCenter(ir);
1523     for (Int_t iphi=1;iphi<=nphi;++iphi) {
1524       Double_t phip = histoRPhi->GetYaxis()->GetBinCenter(iphi);
1525       for (Int_t iz=1;iz<=2;++iz) {
1526         Double_t zp = histoRPhi->GetZaxis()->GetBinCenter(iz);
1527
1528         Double_t charge = 0;
1529
1530         for (Int_t igg = 0; igg<nGGleaks; igg++) { // loop over GG leaks
1531
1532           // A side
1533           Double_t secPos = secPosA[igg];
1534           Double_t radialPos = radialPosA[igg];
1535
1536           if (zp<0) { // C side
1537             secPos = secPosC[igg];
1538             radialPos = radialPosC[igg];
1539           }
1540
1541           // some 'arbitrary' GG leaks
1542           if (  (phip<(TMath::Pi()/9*(secPos+1)) && phip>(TMath::Pi()/9*secPos) ) ) { // sector slice
1543             if ( rp>(radialPos-2.5) && rp<(radialPos+2.5))  // 5 cm slice
1544               charge = 300;
1545           }
1546
1547         }
1548
1549         charge = charge*fgke0; // [C/m^3]
1550
1551         histoRPhi->SetBinContent(ir,iphi,iz,charge);
1552       }
1553     }
1554   }
1555
1556   histoRPhi->Write("SpaceChargeInRPhi");
1557
1558   f->Close();
1559
1560 }
1561
1562
1563 void AliTPCSpaceCharge3D::Print(const Option_t* option) const {
1564   //
1565   // Print function to check the settings of the boundary vectors
1566   // option=="a" prints the C0 and C1 coefficents for calibration purposes
1567   //
1568
1569   TString opt = option; opt.ToLower();
1570   printf("%s\n",GetTitle());
1571   printf(" - Space Charge effect with arbitrary 3D charge density (as input).\n");
1572   printf("   SC correction factor: %f \n",fCorrectionFactor);
1573
1574   if (opt.Contains("a")) { // Print all details
1575     printf(" - T1: %1.4f, T2: %1.4f \n",fT1,fT2);
1576     printf(" - C1: %1.4f, C0: %1.4f \n",fC1,fC0);
1577   }
1578
1579   if (!fInitLookUp) AliError("Lookup table was not initialized! You should do InitSpaceCharge3DDistortion() ...");
1580
1581 }