TPC/TPCbase/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 /// \class AliTPCSpaceCharge3D
  17 /// \brief The class calculates the space point distortions due to an arbitrary space charge distribution in 3D.
  18 ///
  19 /// The method of calculation is based on the analytical solution for the Poisson
  20 /// problem in 3D (cylindrical coordinates). The solution is used in form of
  21 /// look up tables, where the pre calculated solutions for different voxel
  22 /// positions are stored. These voxel solutions can be summed up according
  23 /// to the weight of the position of the applied space charge distribution.
  24 /// Further details can be found in \cite[chap.5]{PhD-thesis_S.Rossegger}.
  25 ///
  26 /// The class also allows a simple scaling of the resulting distortions
  27 /// via the function SetCorrectionFactor. This speeds up the calculation
  28 /// time under the assumption, that the distortions scales linearly with the
  29 /// magnitude of the space charge distribution $\rho(r,z)$ and the shape stays
  30 /// the same at higher luminosities.
  31 ///
  32 /// In contrast to the implementation in 2D (see the class AliTPCSpaceChargeabove),
  33 /// the input charge distribution can be of arbitrary character. An example on how
  34 /// to produce a corresponding charge distribution can be found in the function
  35 /// WriteChargeDistributionToFile. In there, a $\rho(r,z) = (A-B\,z)/r^2$,
  36 /// with slightly different magnitude on the A and C side (due to the muon absorber),
  37 /// is superpositioned with a few leaking wires at arbitrary positions.
  38 ///
  39 /// Marian Ivanov change: 26.06.2013
  40 /// Usage of the realy 3D space charge map as an optional input
  41 /// SetInputSpaceCharge map.
  42 /// In case given map is used 2 2D maps are ignored and  scaling functions  $\rho(r,z) = (A-B\,z)/r^2$,
  43 /// will not work
  44
  45
  46 // End_Html
  47 //
  48 // Begin_Macro(source)
  49 //   {
  50 //   gROOT->SetStyle("Plain"); gStyle->SetPalette(1);
  51 //   TCanvas *c2 = new TCanvas("cAliTPCSpaceCharge3D","cAliTPCSpaceCharge3D",500,400);
  52 //   AliTPCSpaceCharge3D sc;
  53 //   sc.WriteChargeDistributionToFile("SC_zr2_GGleaks.root");
  54 //   sc.SetSCDataFileName("SC_zr2_GGleaks.root");
  55 //   sc.SetOmegaTauT1T2(0,1,1); // B=0
  56 //   sc.InitSpaceCharge3DDistortion();
  57 //   sc.CreateHistoDRinXY(15,300,300)->Draw("colz");
  58 //   return c2;
  59 //   }
  60 // End_Macro
  61 //
  62 // Begin_Html
  63 //   <p>
  64 //   Date: 19/06/2010  <br>
  65 //   Authors: Stefan Rossegger
  66 // End_Html
  67 // _________________________________________________________________
  68
  69
  70 #include "AliMagF.h"
  71 #include "TGeoGlobalMagField.h"
  72 #include "AliTPCcalibDB.h"
  73 #include "AliTPCParam.h"
  74 #include "AliLog.h"
  75 #include "TH2F.h"
  76 #include "TH3F.h"
  77 #include "TFile.h"
  78 #include "TVector.h"
  79 #include "TMatrix.h"
  80 #include "TMatrixD.h"
  81
  82 #include "TMath.h"
  83 #include "AliTPCROC.h"
  84 #include "AliTPCSpaceCharge3D.h"
  85 #include "AliSysInfo.h"
  86
  87 /// \cond CLASSIMP
  88 ClassImp(AliTPCSpaceCharge3D)
  89 /// \endcond
  90
  91 AliTPCSpaceCharge3D::AliTPCSpaceCharge3D()
  92   : AliTPCCorrection("SpaceCharge3D","Space Charge - 3D"),
  93     fC0(0.),fC1(0.),
  94     fCorrectionFactor(1.),
  95     fInitLookUp(kFALSE),
  96     fSCDataFileName(""),
  97     fSCLookUpPOCsFileName3D(""),
  98     fSCLookUpPOCsFileNameRZ(""),
  99     fSCLookUpPOCsFileNameRPhi(""),
 100     fSCdensityInRZ(0),
 101     fSCdensityInRPhiA(0),
 102     fSCdensityInRPhiC(0),
 103     fSpaceChargeHistogram3D(0),
 104     fSpaceChargeHistogramRPhi(0),
 105     fSpaceChargeHistogramRZ(0)
 106 {
 107   //
 108   // default constructor
 109   //
 110
 111   // Array which will contain the solution according to the setted charge density distribution
 112   // see InitSpaceCharge3DDistortion() function
 113   for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
 114     fLookUpErOverEz[k]   =  new TMatrixF(kNR,kNZ);
 115     fLookUpEphiOverEz[k] =  new TMatrixF(kNR,kNZ);
 116     fLookUpDeltaEz[k]    =  new TMatrixF(kNR,kNZ);
 117     fSCdensityDistribution[k] = new TMatrixF(kNR,kNZ);
 118   }
 119   fSCdensityInRZ   = new TMatrixD(kNR,kNZ);
 120   fSCdensityInRPhiA = new TMatrixD(kNR,kNPhi);
 121   fSCdensityInRPhiC = new TMatrixD(kNR,kNPhi);
 122
 123   // location of the precalculated look up tables
 124
 125   fSCLookUpPOCsFileName3D="$(ALICE_ROOT)/TPC/Calib/maps/sc_3D_raw_18-18-26_17p-18p-25p-MN30.root"; // rough estimate
 126   fSCLookUpPOCsFileNameRZ="$(ALICE_ROOT)/TPC/Calib/maps/sc_radSym_35-01-51_34p-01p-50p_MN60.root";
 127   fSCLookUpPOCsFileNameRPhi="$(ALICE_ROOT)/TPC/Calib/maps/sc_cChInZ_35-144-26_34p-18p-01p-MN30.root";
 128   //  fSCLookUpPOCsFileNameRPhi="$(ALICE_ROOT)/TPC/Calib/maps/sc_cChInZ_35-36-26_34p-18p-01p-MN40.root";
 129
 130
 131
 132   // standard location of the space charge distibution ... can be changes
 133   fSCDataFileName="$(ALICE_ROOT)/TPC/Calib/maps/sc_3D_distribution_Sim.root";
 134
 135   //  SetSCDataFileName(fSCDataFileName.Data()); // should be done by the user
 136
 137
 138 }
 139
 140 AliTPCSpaceCharge3D::~AliTPCSpaceCharge3D() {
 141   /// default destructor
 142
 143   for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
 144     delete fLookUpErOverEz[k];
 145     delete fLookUpEphiOverEz[k];
 146     delete fLookUpDeltaEz[k];
 147     delete fSCdensityDistribution[k];
 148   }
 149   delete fSCdensityInRZ;
 150   delete fSCdensityInRPhiA;
 151   delete fSCdensityInRPhiC;
 152   delete fSpaceChargeHistogram3D;
 153   delete fSpaceChargeHistogramRPhi;
 154   delete fSpaceChargeHistogramRZ;
 155 }
 156
 157
 158 void AliTPCSpaceCharge3D::Init() {
 159   /// Initialization funtion
 160
 161   AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
 162   if (!magF) AliError("Magneticd field - not initialized");
 163   Double_t bzField = magF->SolenoidField()/10.; //field in T
 164   AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
 165   if (!param) AliError("Parameters - not initialized");
 166   Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us]   // From dataBase: to be updated: per second (ideally)
 167   Double_t ezField = 400; // [V/cm]   // to be updated: never (hopefully)
 168   Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
 169   // Correction Terms for effective omegaTau; obtained by a laser calibration run
 170   SetOmegaTauT1T2(wt,fT1,fT2);
 171
 172   InitSpaceCharge3DDistortion(); // fill the look up table
 173 }
 174
 175 void AliTPCSpaceCharge3D::Update(const TTimeStamp &/*timeStamp*/) {
 176   /// Update function
 177
 178   AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
 179   if (!magF) AliError("Magneticd field - not initialized");
 180   Double_t bzField = magF->SolenoidField()/10.; //field in T
 181   AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
 182   if (!param) AliError("Parameters - not initialized");
 183   Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us]  // From dataBase: to be updated: per second (ideally)
 184   Double_t ezField = 400; // [V/cm]   // to be updated: never (hopefully)
 185   Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
 186   // Correction Terms for effective omegaTau; obtained by a laser calibration run
 187   SetOmegaTauT1T2(wt,fT1,fT2);
 188
 189   //  SetCorrectionFactor(1.); // should come from some database
 190
 191 }
 192
 193
 194 void AliTPCSpaceCharge3D::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {
 195   /// Calculates the correction due the Space Charge effect within the TPC drift volume
 196
 197   if (!fInitLookUp) {
 198     AliInfo("Lookup table was not initialized! Performing the inizialisation now ...");
 199     InitSpaceCharge3DDistortion();
 200   }
 201
 202   Int_t   order     = 1 ;    // FIXME: hardcoded? Linear interpolation = 1, Quadratic = 2
 203
 204   Float_t intEr, intEphi, intdEz ;
 205   Double_t r, phi, z ;
 206   Int_t    sign;
 207
 208   r      =  TMath::Sqrt( x[0]*x[0] + x[1]*x[1] ) ;
 209   phi    =  TMath::ATan2(x[1],x[0]) ;
 210   if ( phi < 0 ) phi += TMath::TwoPi() ;                   // Table uses phi from 0 to 2*Pi
 211   z      =  x[2] ;                                         // Create temporary copy of x[2]
 212
 213   if ( (roc%36) < 18 ) {
 214     sign =  1;       // (TPC A side)
 215   } else {
 216     sign = -1;       // (TPC C side)
 217   }
 218
 219   if ( sign==1  && z <  fgkZOffSet ) z =  fgkZOffSet;    // Protect against discontinuity at CE
 220   if ( sign==-1 && z > -fgkZOffSet ) z = -fgkZOffSet;    // Protect against discontinuity at CE
 221
 222
 223   if ( (sign==1 && z<0) || (sign==-1 && z>0) ) // just a consistency check
 224     AliError("ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!");
 225
 226   // Get the Er and Ephi field integrals plus the integral over DeltaEz
 227   intEr      = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
 228                                   fgkRList, fgkZList, fgkPhiList, fLookUpErOverEz  );
 229   intEphi    = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
 230                                   fgkRList, fgkZList, fgkPhiList, fLookUpEphiOverEz);
 231   intdEz = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
 232                                   fgkRList, fgkZList, fgkPhiList, fLookUpDeltaEz   );
 233
 234   // Calculate distorted position
 235   if ( r > 0.0 ) {
 236     phi =  phi + fCorrectionFactor *( fC0*intEphi - fC1*intEr ) / r;
 237     r   =  r   + fCorrectionFactor *( fC0*intEr   + fC1*intEphi );
 238   }
 239   Double_t dz = intdEz * fCorrectionFactor * fgkdvdE;
 240
 241   // Calculate correction in cartesian coordinates
 242   dx[0] = - (r * TMath::Cos(phi) - x[0]);
 243   dx[1] = - (r * TMath::Sin(phi) - x[1]);
 244   dx[2] = - dz;  // z distortion - (scaled with driftvelocity dependency on the Ez field and the overall scaling factor)
 245
 246 }
 247
 248 void AliTPCSpaceCharge3D::InitSpaceCharge3DDistortion() {
 249  /// Initialization of the Lookup table which contains the solutions of the
 250  /// "space charge" (poisson) problem - Faster and more accureate
 251  ///
 252  /// Method: Weighted sum-up of the different fields within the look up table
 253  /// but using two lookup tables with higher granularity in the (r,z) and the (rphi)- plane to emulate
 254  /// more realistic space charges. (r,z) from primary ionisation. (rphi) for possible Gating leaks
 255
 256   if (fInitLookUp) {
 257     AliInfo("Lookup table was already initialized!  Doing it again anyway ...");
 258     return;
 259   }
 260
 261   // ------------------------------------------------------------------------------------------------------
 262   // step 1: lookup table in rz, fine grid, radial symetric, to emulate primary ionization
 263
 264   AliInfo("Step 1: Preparation of the weighted look-up tables.");
 265
 266   // lookup table in rz, fine grid
 267
 268   TFile *fZR = new TFile(fSCLookUpPOCsFileNameRZ.Data(),"READ");
 269   if ( !fZR ) {
 270     AliError("Precalculated POC-looup-table in ZR could not be found");
 271     return;
 272   }
 273
 274   // units are in [m]
 275   TVector *gridf1  = (TVector*) fZR->Get("constants");
 276   TVector &grid1 = *gridf1;
 277   TMatrix *coordf1  = (TMatrix*) fZR->Get("coordinates");
 278   TMatrix &coord1 = *coordf1;
 279   TMatrix *coordPOCf1  = (TMatrix*) fZR->Get("POCcoord");
 280   TMatrix &coordPOC1 = *coordPOCf1;
 281
 282   Int_t rows      = (Int_t)grid1(0);   // number of points in r direction  - from RZ or RPhi table
 283   Int_t phiSlices = (Int_t)grid1(1);   // number of points in phi          - from RPhi table
 284   Int_t columns   = (Int_t)grid1(2);   // number of points in z direction  - from RZ table
 285
 286   Float_t gridSizeR   =  (fgkOFCRadius-fgkIFCRadius)/(rows-1);   // unit in [cm]
 287   Float_t gridSizeZ   =  fgkTPCZ0/(columns-1);                  // unit in [cm]
 288
 289   // temporary matrices needed for the calculation // for rotational symmetric RZ table, phislices is 1
 290
 291   TMatrixD *arrayofErA[kNPhiSlices], *arrayofdEzA[kNPhiSlices];
 292   TMatrixD *arrayofErC[kNPhiSlices], *arrayofdEzC[kNPhiSlices];
 293
 294   TMatrixD *arrayofEroverEzA[kNPhiSlices], *arrayofDeltaEzA[kNPhiSlices];
 295   TMatrixD *arrayofEroverEzC[kNPhiSlices], *arrayofDeltaEzC[kNPhiSlices];
 296
 297
 298   for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
 299
 300     arrayofErA[k]   =   new TMatrixD(rows,columns) ;
 301     arrayofdEzA[k]  =   new TMatrixD(rows,columns) ;
 302     arrayofErC[k]   =   new TMatrixD(rows,columns) ;
 303     arrayofdEzC[k]  =   new TMatrixD(rows,columns) ;
 304
 305     arrayofEroverEzA[k]   =   new TMatrixD(rows,columns) ;
 306     arrayofDeltaEzA[k]    =   new TMatrixD(rows,columns) ;
 307     arrayofEroverEzC[k]   =   new TMatrixD(rows,columns) ;
 308     arrayofDeltaEzC[k]    =   new TMatrixD(rows,columns) ;
 309
 310     // zero initialization not necessary, it is done in the constructor of TMatrix
 311
 312   }
 313
 314   // list of points as used during sum up
 315   Double_t  rlist1[kNRows], zedlist1[kNColumns];// , philist1[phiSlices];
 316   for ( Int_t i = 0 ; i < rows ; i++ )    {
 317     rlist1[i] = fgkIFCRadius + i*gridSizeR ;
 318     for ( Int_t j = 0 ; j < columns ; j++ ) {
 319       zedlist1[j]  = j * gridSizeZ ;
 320     }
 321   }
 322
 323   TTree *treePOC = (TTree*)fZR->Get("POCall");
 324
 325   TVector *bEr  = 0;   //TVector *bEphi= 0;
 326   TVector *bEz  = 0;
 327
 328   treePOC->SetBranchAddress("Er",&bEr);
 329   treePOC->SetBranchAddress("Ez",&bEz);
 330
 331
 332   // Read the complete tree and do a weighted sum-up over the POC configurations
 333   // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 334
 335   Int_t treeNumPOC = (Int_t)treePOC->GetEntries(); // Number of POC conf. in the look-up table
 336   Int_t ipC = 0; // POC Conf. counter (note: different to the POC number in the tree!)
 337
 338   for (Int_t itreepC=0; itreepC<treeNumPOC; itreepC++) { // ------------- loop over POC configurations in tree
 339
 340     treePOC->GetEntry(itreepC);
 341
 342     // center of the POC voxel in [meter]
 343     Double_t r0 = coordPOC1(ipC,0);
 344     Double_t phi0 = coordPOC1(ipC,1);
 345     Double_t z0 = coordPOC1(ipC,2);
 346
 347     ipC++; // POC configuration counter
 348
 349     // weights (charge density) at POC position on the A and C side (in C/m^3/e0)
 350     // note: coordinates are in [cm]
 351     Double_t weightA = GetSpaceChargeDensity(r0*100,phi0, z0*100, 1);  // partial load in r,z
 352     Double_t weightC = GetSpaceChargeDensity(r0*100,phi0,-z0*100, 1);  // partial load in r,z
 353
 354     // Summing up the vector components according to their weight
 355
 356     Int_t ip = 0;
 357     for ( Int_t j = 0 ; j < columns ; j++ ) {
 358       for ( Int_t i = 0 ; i < rows ; i++ )    {
 359         for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
 360
 361           // check wether the coordinates were screwed
 362           if (TMath::Abs((coord1(0,ip)*100-rlist1[i]))>1 ||
 363               TMath::Abs((coord1(2,ip)*100-zedlist1[j])>1)) {
 364             AliError("internal error: coordinate system was screwed during the sum-up");
 365             printf("sum-up: (r,z)=(%f,%f)\n",rlist1[i],zedlist1[j]);
 366             printf("lookup: (r,z)=(%f,%f)\n",coord1(0,ip)*100,coord1(2,ip)*100);
 367             AliError("Don't trust the results of the space charge calculation!");
 368           }
 369
 370           // unfortunately, the lookup tables were produced to be faster for phi symmetric charges
 371           // This will be the most frequent usage (hopefully)
 372           // That's why we have to do this here ...
 373
 374           TMatrixD &erA   =  *arrayofErA[k]  ;
 375           TMatrixD &dEzA  =  *arrayofdEzA[k]   ;
 376
 377           TMatrixD &erC   =  *arrayofErC[k]  ;
 378           TMatrixD &dEzC  =  *arrayofdEzC[k]   ;
 379
 380           // Sum up - Efield values in [V/m] -> transition to [V/cm]
 381           erA(i,j) += ((*bEr)(ip)) * weightA /100;
 382           erC(i,j) += ((*bEr)(ip)) * weightC /100;
 383           dEzA(i,j)  += ((*bEz)(ip)) * weightA /100;
 384           dEzC(i,j)  += ((*bEz)(ip)) * weightC /100;
 385
 386           // increase the counter
 387           ip++;
 388         }
 389       }
 390     } // end coordinate loop
 391   } // end POC loop
 392
 393
 394   // -------------------------------------------------------------------------------
 395   // Division by the Ez (drift) field and integration along z
 396
 397   //  AliInfo("Step 1: Division and integration");
 398
 399   Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = Electric Field (V/cm) Magnitude ~ -400 V/cm;
 400
 401   for ( Int_t k = 0 ; k < phiSlices ; k++ ) { // phi loop
 402
 403     // matrices holding the solution - summation of POC charges // see above
 404     TMatrixD &erA   =  *arrayofErA[k]  ;
 405     TMatrixD &dezA  =  *arrayofdEzA[k]   ;
 406     TMatrixD &erC   =  *arrayofErC[k]  ;
 407     TMatrixD &dezC  =  *arrayofdEzC[k]   ;
 408
 409     // matrices which will contain the integrated fields (divided by the drift field)
 410     TMatrixD &erOverEzA   =  *arrayofEroverEzA[k]  ;
 411     TMatrixD &deltaEzA    =  *arrayofDeltaEzA[k];
 412     TMatrixD &erOverEzC   =  *arrayofEroverEzC[k]  ;
 413     TMatrixD &deltaEzC    =  *arrayofDeltaEzC[k];
 414
 415     for ( Int_t i = 0 ; i < rows ; i++ )    { // r loop
 416       for ( Int_t j = columns-1 ; j >= 0 ; j-- ) {// z loop
 417         // Count backwards to facilitate integration over Z
 418
 419         Int_t index = 1 ; // Simpsons rule if N=odd.If N!=odd then add extra point
 420                           // by trapezoidal rule.
 421
 422         erOverEzA(i,j) = 0; //ephiOverEzA(i,j) = 0;
 423         deltaEzA(i,j) = 0;
 424         erOverEzC(i,j) = 0; //ephiOverEzC(i,j) = 0;
 425         deltaEzC(i,j) = 0;
 426
 427         for ( Int_t m = j ; m < columns ; m++ ) { // integration
 428
 429           erOverEzA(i,j)   += index*(gridSizeZ/3.0)*erA(i,m)/(-1*ezField) ;
 430           erOverEzC(i,j)   += index*(gridSizeZ/3.0)*erC(i,m)/(-1*ezField)  ;
 431           deltaEzA(i,j)    += index*(gridSizeZ/3.0)*dezA(i,m)/(-1) ;
 432           deltaEzC(i,j)    += index*(gridSizeZ/3.0)*dezC(i,m)/(-1) ;
 433
 434           if ( index != 4 )  index = 4; else index = 2 ;
 435
 436         }
 437
 438         if ( index == 4 ) {
 439           erOverEzA(i,j)   -= (gridSizeZ/3.0)*erA(i,columns-1)/(-1*ezField) ;
 440           erOverEzC(i,j)   -= (gridSizeZ/3.0)*erC(i,columns-1)/(-1*ezField) ;
 441           deltaEzA(i,j)    -= (gridSizeZ/3.0)*dezA(i,columns-1)/(-1) ;
 442           deltaEzC(i,j)    -= (gridSizeZ/3.0)*dezC(i,columns-1)/(-1) ;
 443         }
 444         if ( index == 2 ) {
 445           erOverEzA(i,j)   += (gridSizeZ/3.0)*(0.5*erA(i,columns-2)-2.5*erA(i,columns-1))/(-1*ezField) ;
 446           erOverEzC(i,j)   += (gridSizeZ/3.0)*(0.5*erC(i,columns-2)-2.5*erC(i,columns-1))/(-1*ezField) ;
 447           deltaEzA(i,j)    += (gridSizeZ/3.0)*(0.5*dezA(i,columns-2)-2.5*dezA(i,columns-1))/(-1) ;
 448           deltaEzC(i,j)    += (gridSizeZ/3.0)*(0.5*dezC(i,columns-2)-2.5*dezC(i,columns-1))/(-1) ;
 449         }
 450         if ( j == columns-2 ) {
 451           erOverEzA(i,j)   = (gridSizeZ/3.0)*(1.5*erA(i,columns-2)+1.5*erA(i,columns-1))/(-1*ezField) ;
 452           erOverEzC(i,j)   = (gridSizeZ/3.0)*(1.5*erC(i,columns-2)+1.5*erC(i,columns-1))/(-1*ezField) ;
 453           deltaEzA(i,j)    = (gridSizeZ/3.0)*(1.5*dezA(i,columns-2)+1.5*dezA(i,columns-1))/(-1) ;
 454           deltaEzC(i,j)    = (gridSizeZ/3.0)*(1.5*dezC(i,columns-2)+1.5*dezC(i,columns-1))/(-1) ;
 455         }
 456         if ( j == columns-1 ) {
 457           erOverEzA(i,j)   = 0;
 458           erOverEzC(i,j)   = 0;
 459           deltaEzA(i,j)    = 0;
 460           deltaEzC(i,j)    = 0;
 461         }
 462       }
 463     }
 464
 465   }
 466
 467   //  AliInfo("Step 1: Interpolation to Standard grid");
 468
 469   // -------------------------------------------------------------------------------
 470   // Interpolate results onto the standard grid which is used for all AliTPCCorrections classes
 471
 472   const Int_t order  = 1  ;  // Linear interpolation = 1, Quadratic = 2
 473
 474   Double_t  r, z;//phi, z ;
 475   for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
 476     //    phi = fgkPhiList[k] ;
 477
 478     // final lookup table
 479     TMatrixF &erOverEzFinal   =  *fLookUpErOverEz[k]  ;
 480     TMatrixF &deltaEzFinal    =  *fLookUpDeltaEz[k]   ;
 481
 482     // calculated and integrated tables - just one phi slice
 483     TMatrixD &erOverEzA   =  *arrayofEroverEzA[0]  ;
 484     TMatrixD &deltaEzA    =  *arrayofDeltaEzA[0];
 485     TMatrixD &erOverEzC   =  *arrayofEroverEzC[0]  ;
 486     TMatrixD &deltaEzC    =  *arrayofDeltaEzC[0];
 487
 488
 489     for ( Int_t j = 0 ; j < kNZ ; j++ ) {
 490
 491       z = TMath::Abs(fgkZList[j]) ;  // z position is symmetric
 492
 493       for ( Int_t i = 0 ; i < kNR ; i++ ) {
 494         r = fgkRList[i] ;
 495
 496         // Interpolate Lookup tables onto standard grid
 497         if (fgkZList[j]>0) {
 498           erOverEzFinal(i,j)   = Interpolate2DTable(order, r, z, rows, columns, rlist1, zedlist1, erOverEzA  );
 499           deltaEzFinal(i,j)    = Interpolate2DTable(order, r, z, rows, columns, rlist1, zedlist1, deltaEzA   );
 500         } else {
 501           erOverEzFinal(i,j)   = Interpolate2DTable(order, r, z, rows, columns, rlist1, zedlist1, erOverEzC  );
 502           deltaEzFinal(i,j)    = - Interpolate2DTable(order, r, z, rows, columns, rlist1, zedlist1, deltaEzC   );
 503           // negative coordinate system on C side
 504         }
 505
 506       } // end r loop
 507     } // end z loop
 508   } // end phi loop
 509
 510
 511   // clear the temporary arrays lists
 512   for ( Int_t k = 0 ; k < phiSlices ; k++ )  {
 513
 514     delete arrayofErA[k];
 515     delete arrayofdEzA[k];
 516     delete arrayofErC[k];
 517     delete arrayofdEzC[k];
 518
 519     delete arrayofEroverEzA[k];
 520     delete arrayofDeltaEzA[k];
 521     delete arrayofEroverEzC[k];
 522     delete arrayofDeltaEzC[k];
 523
 524   }
 525
 526   fZR->Close();
 527
 528   // ------------------------------------------------------------------------------------------------------
 529   // Step 2: Load and sum up lookup table in rphi, fine grid, to emulate for example a GG leak
 530
 531   //  AliInfo("Step 2: Preparation of the weighted look-up table");
 532
 533   TFile *fRPhi = new TFile(fSCLookUpPOCsFileNameRPhi.Data(),"READ");
 534   if ( !fRPhi ) {
 535     AliError("Precalculated POC-looup-table in RPhi could not be found");
 536     return;
 537   }
 538
 539   // units are in [m]
 540   TVector *gridf2  = (TVector*) fRPhi->Get("constants");
 541   TVector &grid2 = *gridf2;
 542   TMatrix *coordf2  = (TMatrix*) fRPhi->Get("coordinates");
 543   TMatrix &coord2 = *coordf2;
 544   TMatrix *coordPOCf2  = (TMatrix*) fRPhi->Get("POCcoord");
 545   TMatrix &coordPOC2 = *coordPOCf2;
 546
 547   rows      = (Int_t)grid2(0);   // number of points in r direction
 548   phiSlices = (Int_t)grid2(1);   // number of points in phi
 549   columns   = (Int_t)grid2(2);   // number of points in z direction
 550
 551   gridSizeR   =  (fgkOFCRadius-fgkIFCRadius)/(rows-1);   // unit in [cm]
 552   Float_t gridSizePhi =  TMath::TwoPi()/phiSlices;         // unit in [rad]
 553   gridSizeZ   =  fgkTPCZ0/(columns-1);                  // unit in [cm]
 554
 555   // list of points as used during sum up
 556   Double_t  rlist2[kNRows], philist2[kNPhiSlices], zedlist2[kNColumns];
 557   for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
 558     philist2[k] =  gridSizePhi * k;
 559     for ( Int_t i = 0 ; i < rows ; i++ )    {
 560       rlist2[i] = fgkIFCRadius + i*gridSizeR ;
 561       for ( Int_t j = 0 ; j < columns ; j++ ) {
 562         zedlist2[j]  = j * gridSizeZ ;
 563       }
 564     }
 565   } // only done once
 566
 567   // temporary matrices needed for the calculation
 568
 569   TMatrixD *arrayofErA2[kNPhiSlices], *arrayofEphiA2[kNPhiSlices], *arrayofdEzA2[kNPhiSlices];
 570   TMatrixD *arrayofErC2[kNPhiSlices], *arrayofEphiC2[kNPhiSlices], *arrayofdEzC2[kNPhiSlices];
 571
 572   TMatrixD *arrayofEroverEzA2[kNPhiSlices], *arrayofEphioverEzA2[kNPhiSlices], *arrayofDeltaEzA2[kNPhiSlices];
 573   TMatrixD *arrayofEroverEzC2[kNPhiSlices], *arrayofEphioverEzC2[kNPhiSlices], *arrayofDeltaEzC2[kNPhiSlices];
 574
 575
 576   for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
 577
 578     arrayofErA2[k]   =   new TMatrixD(rows,columns) ;
 579     arrayofEphiA2[k] =   new TMatrixD(rows,columns) ;
 580     arrayofdEzA2[k]  =   new TMatrixD(rows,columns) ;
 581     arrayofErC2[k]   =   new TMatrixD(rows,columns) ;
 582     arrayofEphiC2[k] =   new TMatrixD(rows,columns) ;
 583     arrayofdEzC2[k]  =   new TMatrixD(rows,columns) ;
 584
 585     arrayofEroverEzA2[k]   =   new TMatrixD(rows,columns) ;
 586     arrayofEphioverEzA2[k] =   new TMatrixD(rows,columns) ;
 587     arrayofDeltaEzA2[k]    =   new TMatrixD(rows,columns) ;
 588     arrayofEroverEzC2[k]   =   new TMatrixD(rows,columns) ;
 589     arrayofEphioverEzC2[k] =   new TMatrixD(rows,columns) ;
 590     arrayofDeltaEzC2[k]    =   new TMatrixD(rows,columns) ;
 591
 592     // zero initialization not necessary, it is done in the constructor of TMatrix
 593
 594   }
 595
 596
 597   treePOC = (TTree*)fRPhi->Get("POCall");
 598
 599   //  TVector *bEr  = 0;   // done above
 600   TVector *bEphi= 0;
 601   //  TVector *bEz  = 0;   // done above
 602
 603   treePOC->SetBranchAddress("Er",&bEr);
 604   treePOC->SetBranchAddress("Ephi",&bEphi);
 605   treePOC->SetBranchAddress("Ez",&bEz);
 606
 607   // Read the complete tree and do a weighted sum-up over the POC configurations
 608   // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 609
 610   treeNumPOC = (Int_t)treePOC->GetEntries(); // Number of POC conf. in the look-up table
 611   ipC = 0; // POC Conf. counter (note: different to the POC number in the tree!)
 612
 613   for (Int_t itreepC=0; itreepC<treeNumPOC; itreepC++) { // ------------- loop over POC configurations in tree
 614
 615     treePOC->GetEntry(itreepC);
 616
 617     // center of the POC voxel in [meter]
 618     Double_t r0 = coordPOC2(ipC,0);
 619     Double_t phi0 = coordPOC2(ipC,1);
 620     //    Double_t z0 = coordPOC2(ipC,2);
 621
 622      // weights (charge density) at POC position on the A and C side (in C/m^3/e0)
 623     // note: coordinates are in [cm]
 624     Double_t weightA = GetSpaceChargeDensity(r0*100,phi0, 0.499, 2);  // partial load in r,phi
 625     Double_t weightC = GetSpaceChargeDensity(r0*100,phi0,-0.499, 2);  // partial load in r,phi
 626
 627     //    printf("-----\n%f %f : %e %e\n",r0,phi0,weightA,weightC);
 628
 629     // Summing up the vector components according to their weight
 630
 631     Int_t ip = 0;
 632     for ( Int_t j = 0 ; j < columns ; j++ ) {
 633       for ( Int_t i = 0 ; i < rows ; i++ )    {
 634         for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
 635
 636           // check wether the coordinates were screwed
 637           if (TMath::Abs((coord2(0,ip)*100-rlist2[i]))>1 ||
 638               TMath::Abs((coord2(1,ip)-philist2[k]))>1 ||
 639               TMath::Abs((coord2(2,ip)*100-zedlist2[j]))>1) {
 640             AliError("internal error: coordinate system was screwed during the sum-up");
 641             printf("lookup: (r,phi,z)=(%f,%f,%f)\n",coord2(0,ip)*100,coord2(1,ip),coord2(2,ip)*100);
 642             printf("sum-up: (r,phi,z)=(%f,%f,%f)\n",rlist2[i],philist2[k],zedlist2[j]);
 643             AliError("Don't trust the results of the space charge calculation!");
 644           }
 645
 646           // unfortunately, the lookup tables were produced to be faster for phi symmetric charges
 647           // This will be the most frequent usage (hopefully)
 648           // That's why we have to do this here ...
 649
 650           TMatrixD &erA   =  *arrayofErA2[k]  ;
 651           TMatrixD &ephiA =  *arrayofEphiA2[k];
 652           TMatrixD &dEzA  =  *arrayofdEzA2[k] ;
 653
 654           TMatrixD &erC   =  *arrayofErC2[k]  ;
 655           TMatrixD &ephiC =  *arrayofEphiC2[k];
 656           TMatrixD &dEzC  =  *arrayofdEzC2[k]   ;
 657
 658           // Sum up - Efield values in [V/m] -> transition to [V/cm]
 659           erA(i,j) += ((*bEr)(ip)) * weightA /100;
 660           erC(i,j) += ((*bEr)(ip)) * weightC /100;
 661           ephiA(i,j) += ((*bEphi)(ip)) * weightA/100; // [V/rad]
 662           ephiC(i,j) += ((*bEphi)(ip)) * weightC/100; // [V/rad]
 663           dEzA(i,j)  += ((*bEz)(ip)) * weightA /100;
 664           dEzC(i,j)  += ((*bEz)(ip)) * weightC /100;
 665
 666           // increase the counter
 667           ip++;
 668         }
 669       }
 670     } // end coordinate loop
 671
 672
 673     // Rotation and summation in the rest of the dPhiSteps
 674     // which were not stored in the this tree due to storage & symmetry reasons
 675
 676
 677     Int_t phiPoints = (Int_t) grid2(1);
 678     Int_t phiPOC    = (Int_t) grid2(4);
 679
 680     //   printf("%d %d\n",phiPOC,flagRadSym);
 681
 682     for (Int_t phiiC = 1; phiiC<phiPOC; phiiC++) { // just used for non-radial symetric table
 683
 684       Double_t phi0R = phiiC*phi0*2 + phi0; // rotate further
 685
 686       // weights (charge density) at POC position on the A and C side (in C/m^3/e0)
 687       // note: coordinates are in [cm] // ecxept z
 688       weightA = GetSpaceChargeDensity(r0*100,phi0R, 0.499, 2);  // partial load in r,phi
 689       weightC = GetSpaceChargeDensity(r0*100,phi0R,-0.499, 2);  // partial load in r,phi
 690
 691       // printf("%f %f : %e %e\n",r0,phi0R,weightA,weightC);
 692
 693       // Summing up the vector components according to their weight
 694       ip = 0;
 695       for ( Int_t j = 0 ; j < columns ; j++ ) {
 696         for ( Int_t i = 0 ; i < rows ; i++ )    {
 697           for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
 698
 699             // Note: rotating the coordinated during the sum up
 700
 701             Int_t rotVal = (phiPoints/phiPOC)*phiiC;
 702             Int_t ipR = -1;
 703
 704             if ((ip%phiPoints)>=rotVal) {
 705               ipR = ip-rotVal;
 706             } else {
 707               ipR = ip+(phiPoints-rotVal);
 708             }
 709
 710             // unfortunately, the lookup tables were produced to be faster for phi symmetric charges
 711             // This will be the most frequent usage
 712             // That's why we have to do this here and not outside the loop ...
 713
 714             TMatrixD &erA   =  *arrayofErA2[k]  ;
 715             TMatrixD &ephiA =  *arrayofEphiA2[k];
 716             TMatrixD &dEzA  =  *arrayofdEzA2[k]   ;
 717
 718             TMatrixD &erC   =  *arrayofErC2[k]  ;
 719             TMatrixD &ephiC =  *arrayofEphiC2[k];
 720             TMatrixD &dEzC  =  *arrayofdEzC2[k]   ;
 721
 722             // Sum up - Efield values in [V/m] -> transition to [V/cm]
 723             erA(i,j) += ((*bEr)(ipR)) * weightA /100;
 724             erC(i,j) += ((*bEr)(ipR)) * weightC /100;
 725             ephiA(i,j) += ((*bEphi)(ipR)) * weightA/100; // [V/rad]
 726             ephiC(i,j) += ((*bEphi)(ipR)) * weightC/100; // [V/rad]
 727             dEzA(i,j)  += ((*bEz)(ipR)) * weightA /100;
 728             dEzC(i,j)  += ((*bEz)(ipR)) * weightC /100;
 729
 730             // increase the counter
 731             ip++;
 732           }
 733         }
 734       } // end coordinate loop
 735
 736     } // end phi-POC summation (phiiC)
 737
 738     ipC++; // POC configuration counter
 739
 740     //   printf("POC: (r,phi,z) = (%f %f %f) | weight(A,C): %03.1lf %03.1lf\n",r0,phi0,z0, weightA, weightC);
 741
 742   }
 743
 744
 745
 746
 747   // -------------------------------------------------------------------------------
 748   // Division by the Ez (drift) field and integration along z
 749
 750   //  AliInfo("Step 2: Division and integration");
 751
 752
 753   for ( Int_t k = 0 ; k < phiSlices ; k++ ) { // phi loop
 754
 755     // matrices holding the solution - summation of POC charges // see above
 756     TMatrixD &erA   =  *arrayofErA2[k]  ;
 757     TMatrixD &ephiA =  *arrayofEphiA2[k];
 758     TMatrixD &dezA  =  *arrayofdEzA2[k]   ;
 759     TMatrixD &erC   =  *arrayofErC2[k]  ;
 760     TMatrixD &ephiC =  *arrayofEphiC2[k];
 761     TMatrixD &dezC  =  *arrayofdEzC2[k]   ;
 762
 763     // matrices which will contain the integrated fields (divided by the drift field)
 764     TMatrixD &erOverEzA   =  *arrayofEroverEzA2[k]  ;
 765     TMatrixD &ephiOverEzA =  *arrayofEphioverEzA2[k];
 766     TMatrixD &deltaEzA    =  *arrayofDeltaEzA2[k];
 767     TMatrixD &erOverEzC   =  *arrayofEroverEzC2[k]  ;
 768     TMatrixD &ephiOverEzC =  *arrayofEphioverEzC2[k];
 769     TMatrixD &deltaEzC    =  *arrayofDeltaEzC2[k];
 770
 771     for ( Int_t i = 0 ; i < rows ; i++ )    { // r loop
 772       for ( Int_t j = columns-1 ; j >= 0 ; j-- ) {// z loop
 773         // Count backwards to facilitate integration over Z
 774
 775         Int_t index = 1 ; // Simpsons rule if N=odd.If N!=odd then add extra point by trapezoidal rule.
 776
 777         erOverEzA(i,j) = 0;
 778         ephiOverEzA(i,j) = 0;
 779         deltaEzA(i,j) = 0;
 780         erOverEzC(i,j) = 0;
 781         ephiOverEzC(i,j) = 0;
 782         deltaEzC(i,j) = 0;
 783
 784         for ( Int_t m = j ; m < columns ; m++ ) { // integration
 785
 786           erOverEzA(i,j)   += index*(gridSizeZ/3.0)*erA(i,m)/(-1*ezField) ;
 787           erOverEzC(i,j)   += index*(gridSizeZ/3.0)*erC(i,m)/(-1*ezField)  ;
 788           ephiOverEzA(i,j) += index*(gridSizeZ/3.0)*ephiA(i,m)/(-1*ezField)  ;
 789           ephiOverEzC(i,j) += index*(gridSizeZ/3.0)*ephiC(i,m)/(-1*ezField)  ;
 790           deltaEzA(i,j)    += index*(gridSizeZ/3.0)*dezA(i,m)/(-1) ;
 791           deltaEzC(i,j)    += index*(gridSizeZ/3.0)*dezC(i,m)/(-1) ;
 792
 793           if ( index != 4 )  index = 4; else index = 2 ;
 794
 795         }
 796
 797         if ( index == 4 ) {
 798           erOverEzA(i,j)   -= (gridSizeZ/3.0)*erA(i,columns-1)/(-1*ezField) ;
 799           erOverEzC(i,j)   -= (gridSizeZ/3.0)*erC(i,columns-1)/(-1*ezField) ;
 800           ephiOverEzA(i,j) -= (gridSizeZ/3.0)*ephiA(i,columns-1)/(-1*ezField) ;
 801           ephiOverEzC(i,j) -= (gridSizeZ/3.0)*ephiC(i,columns-1)/(-1*ezField) ;
 802           deltaEzA(i,j)    -= (gridSizeZ/3.0)*dezA(i,columns-1)/(-1) ;
 803           deltaEzC(i,j)    -= (gridSizeZ/3.0)*dezC(i,columns-1)/(-1) ;
 804         }
 805         if ( index == 2 ) {
 806           erOverEzA(i,j)   += (gridSizeZ/3.0)*(0.5*erA(i,columns-2)-2.5*erA(i,columns-1))/(-1*ezField) ;
 807           erOverEzC(i,j)   += (gridSizeZ/3.0)*(0.5*erC(i,columns-2)-2.5*erC(i,columns-1))/(-1*ezField) ;
 808           ephiOverEzA(i,j) += (gridSizeZ/3.0)*(0.5*ephiA(i,columns-2)-2.5*ephiA(i,columns-1))/(-1*ezField) ;
 809           ephiOverEzC(i,j) += (gridSizeZ/3.0)*(0.5*ephiC(i,columns-2)-2.5*ephiC(i,columns-1))/(-1*ezField) ;
 810           deltaEzA(i,j)    += (gridSizeZ/3.0)*(0.5*dezA(i,columns-2)-2.5*dezA(i,columns-1))/(-1) ;
 811           deltaEzC(i,j)    += (gridSizeZ/3.0)*(0.5*dezC(i,columns-2)-2.5*dezC(i,columns-1))/(-1) ;
 812         }
 813         if ( j == columns-2 ) {
 814           erOverEzA(i,j)   = (gridSizeZ/3.0)*(1.5*erA(i,columns-2)+1.5*erA(i,columns-1))/(-1*ezField) ;
 815           erOverEzC(i,j)   = (gridSizeZ/3.0)*(1.5*erC(i,columns-2)+1.5*erC(i,columns-1))/(-1*ezField) ;
 816           ephiOverEzA(i,j) = (gridSizeZ/3.0)*(1.5*ephiA(i,columns-2)+1.5*ephiA(i,columns-1))/(-1*ezField) ;
 817           ephiOverEzC(i,j) = (gridSizeZ/3.0)*(1.5*ephiC(i,columns-2)+1.5*ephiC(i,columns-1))/(-1*ezField) ;
 818           deltaEzA(i,j)    = (gridSizeZ/3.0)*(1.5*dezA(i,columns-2)+1.5*dezA(i,columns-1))/(-1) ;
 819           deltaEzC(i,j)    = (gridSizeZ/3.0)*(1.5*dezC(i,columns-2)+1.5*dezC(i,columns-1))/(-1) ;
 820         }
 821         if ( j == columns-1 ) {
 822           erOverEzA(i,j)   = 0;
 823           erOverEzC(i,j)   = 0;
 824           ephiOverEzA(i,j) = 0;
 825           ephiOverEzC(i,j) = 0;
 826           deltaEzA(i,j)    = 0;
 827           deltaEzC(i,j)    = 0;
 828         }
 829       }
 830     }
 831
 832   }
 833
 834   AliInfo("Step 2: Interpolation to Standard grid");
 835
 836   // -------------------------------------------------------------------------------
 837   // Interpolate results onto the standard grid which is used for all AliTPCCorrections classes
 838
 839
 840   for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
 841     Double_t phi = fgkPhiList[k] ;
 842
 843     // final lookup table
 844     TMatrixF &erOverEzFinal   =  *fLookUpErOverEz[k]  ;
 845     TMatrixF &ephiOverEzFinal =  *fLookUpEphiOverEz[k];
 846     TMatrixF &deltaEzFinal    =  *fLookUpDeltaEz[k]   ;
 847
 848     for ( Int_t j = 0 ; j < kNZ ; j++ ) {
 849
 850       z = TMath::Abs(fgkZList[j]) ;  // z position is symmetric
 851
 852       for ( Int_t i = 0 ; i < kNR ; i++ ) {
 853         r = fgkRList[i] ;
 854
 855         // Interpolate Lookup tables onto standard grid
 856         if (fgkZList[j]>0) {
 857           erOverEzFinal(i,j)   += Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
 858                                                rlist2, zedlist2, philist2, arrayofEroverEzA2  );
 859           ephiOverEzFinal(i,j) += Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
 860                                                rlist2, zedlist2, philist2, arrayofEphioverEzA2);
 861           deltaEzFinal(i,j)    += Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
 862                                                rlist2, zedlist2, philist2, arrayofDeltaEzA2   );
 863         } else {
 864           erOverEzFinal(i,j)   += Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
 865                                                rlist2, zedlist2, philist2, arrayofEroverEzC2  );
 866           ephiOverEzFinal(i,j) += Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
 867                                                rlist2, zedlist2, philist2, arrayofEphioverEzC2);
 868           deltaEzFinal(i,j)  -=  Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
 869                                                rlist2, zedlist2, philist2, arrayofDeltaEzC2   );
 870         }
 871
 872       } // end r loop
 873     } // end z loop
 874   } // end phi loop
 875
 876
 877   // clear the temporary arrays lists
 878   for ( Int_t k = 0 ; k < phiSlices ; k++ )  {
 879
 880     delete arrayofErA2[k];
 881     delete arrayofEphiA2[k];
 882     delete arrayofdEzA2[k];
 883     delete arrayofErC2[k];
 884     delete arrayofEphiC2[k];
 885     delete arrayofdEzC2[k];
 886
 887     delete arrayofEroverEzA2[k];
 888     delete arrayofEphioverEzA2[k];
 889     delete arrayofDeltaEzA2[k];
 890     delete arrayofEroverEzC2[k];
 891     delete arrayofEphioverEzC2[k];
 892     delete arrayofDeltaEzC2[k];
 893
 894   }
 895
 896   fRPhi->Close();
 897
 898   // FINISHED
 899
 900   fInitLookUp = kTRUE;
 901
 902 }
 903
 904 void AliTPCSpaceCharge3D::InitSpaceCharge3DDistortionCourse() {
 905   /// Initialization of the Lookup table which contains the solutions of the
 906   /// "space charge" (poisson) problem
 907   ///
 908   /// The sum-up uses a look-up table which contains different discretized Space charge fields
 909   /// in order to calculate the corresponding field deviations due to a given (discretized)
 910   /// space charge distribution ....
 911   ///
 912   /// Method of calculation: Weighted sum-up of the different fields within the look up table
 913   /// Note: Full 3d version: Course and slow ...
 914
 915   if (fInitLookUp) {
 916     AliInfo("Lookup table was already initialized!");
 917     //    return;
 918   }
 919
 920   AliInfo("Preparation of the weighted look-up table");
 921
 922   TFile *f = new TFile(fSCLookUpPOCsFileName3D.Data(),"READ");
 923   if ( !f ) {
 924     AliError("Precalculated POC-looup-table could not be found");
 925     return;
 926   }
 927
 928   // units are in [m]
 929   TVector *gridf  = (TVector*) f->Get("constants");
 930   TVector &grid = *gridf;
 931   TMatrix *coordf  = (TMatrix*) f->Get("coordinates");
 932   TMatrix &coord = *coordf;
 933   TMatrix *coordPOCf  = (TMatrix*) f->Get("POCcoord");
 934   TMatrix &coordPOC = *coordPOCf;
 935
 936   Bool_t flagRadSym = 0;
 937   if (grid(1)==1 && grid(4)==1) {
 938     //   AliInfo("LOOK UP TABLE IS RADIAL SYMETTRIC - Field in Phi is ZERO");
 939     flagRadSym=1;
 940   }
 941
 942   Int_t rows      = (Int_t)grid(0);   // number of points in r direction
 943   Int_t phiSlices = (Int_t)grid(1);   // number of points in phi
 944   Int_t columns   = (Int_t)grid(2);   // number of points in z direction
 945
 946   const Float_t gridSizeR   =  (fgkOFCRadius-fgkIFCRadius)/(rows-1);   // unit in [cm]
 947   const Float_t gridSizePhi =  TMath::TwoPi()/phiSlices;         // unit in [rad]
 948   const Float_t gridSizeZ   =  fgkTPCZ0/(columns-1);                  // unit in [cm]
 949
 950   // temporary matrices needed for the calculation
 951   TMatrixD *arrayofErA[kNPhiSlices], *arrayofEphiA[kNPhiSlices], *arrayofdEzA[kNPhiSlices];
 952   TMatrixD *arrayofErC[kNPhiSlices], *arrayofEphiC[kNPhiSlices], *arrayofdEzC[kNPhiSlices];
 953
 954   TMatrixD *arrayofEroverEzA[kNPhiSlices], *arrayofEphioverEzA[kNPhiSlices], *arrayofDeltaEzA[kNPhiSlices];
 955   TMatrixD *arrayofEroverEzC[kNPhiSlices], *arrayofEphioverEzC[kNPhiSlices], *arrayofDeltaEzC[kNPhiSlices];
 956
 957
 958   for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
 959
 960     arrayofErA[k]   =   new TMatrixD(rows,columns) ;
 961     arrayofEphiA[k] =   new TMatrixD(rows,columns) ; // zero if radial symmetric
 962     arrayofdEzA[k]  =   new TMatrixD(rows,columns) ;
 963     arrayofErC[k]   =   new TMatrixD(rows,columns) ;
 964     arrayofEphiC[k] =   new TMatrixD(rows,columns) ; // zero if radial symmetric
 965     arrayofdEzC[k]  =   new TMatrixD(rows,columns) ;
 966
 967     arrayofEroverEzA[k]   =   new TMatrixD(rows,columns) ;
 968     arrayofEphioverEzA[k] =   new TMatrixD(rows,columns) ; // zero if radial symmetric
 969     arrayofDeltaEzA[k]    =   new TMatrixD(rows,columns) ;
 970     arrayofEroverEzC[k]   =   new TMatrixD(rows,columns) ;
 971     arrayofEphioverEzC[k] =   new TMatrixD(rows,columns) ; // zero if radial symmetric
 972     arrayofDeltaEzC[k]    =   new TMatrixD(rows,columns) ;
 973
 974     // Set the values to zero the lookup tables
 975     // not necessary, it is done in the constructor of TMatrix - code deleted
 976
 977   }
 978
 979   // list of points as used in the interpolation (during sum up)
 980   Double_t  rlist[kNRows], zedlist[kNColumns] , philist[kNPhiSlices];
 981   for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
 982     philist[k] =  gridSizePhi * k;
 983     for ( Int_t i = 0 ; i < rows ; i++ )    {
 984       rlist[i] = fgkIFCRadius + i*gridSizeR ;
 985       for ( Int_t j = 0 ; j < columns ; j++ ) {
 986         zedlist[j]  = j * gridSizeZ ;
 987       }
 988     }
 989   } // only done once
 990
 991
 992   TTree *treePOC = (TTree*)f->Get("POCall");
 993
 994   TVector *bEr  = 0;   TVector *bEphi= 0;   TVector *bEz  = 0;
 995
 996   treePOC->SetBranchAddress("Er",&bEr);
 997   if (!flagRadSym) treePOC->SetBranchAddress("Ephi",&bEphi);
 998   treePOC->SetBranchAddress("Ez",&bEz);
 999
1000
1001   // Read the complete tree and do a weighted sum-up over the POC configurations
1002   // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1003
1004   Int_t treeNumPOC = (Int_t)treePOC->GetEntries(); // Number of POC conf. in the look-up table
1005   Int_t ipC = 0; // POC Conf. counter (note: different to the POC number in the tree!)
1006
1007   for (Int_t itreepC=0; itreepC<treeNumPOC; itreepC++) { // ------------- loop over POC configurations in tree
1008
1009     treePOC->GetEntry(itreepC);
1010
1011     // center of the POC voxel in [meter]
1012     Double_t r0 = coordPOC(ipC,0);
1013     Double_t phi0 = coordPOC(ipC,1);
1014     Double_t z0 = coordPOC(ipC,2);
1015
1016     ipC++; // POC configuration counter
1017
1018     // weights (charge density) at POC position on the A and C side (in C/m^3/e0)
1019     // note: coordinates are in [cm]
1020     Double_t weightA = GetSpaceChargeDensity(r0*100,phi0, z0*100,0);
1021     Double_t weightC = GetSpaceChargeDensity(r0*100,phi0,-z0*100,0);
1022
1023     // Summing up the vector components according to their weight
1024
1025     Int_t ip = 0;
1026     for ( Int_t j = 0 ; j < columns ; j++ ) {
1027       for ( Int_t i = 0 ; i < rows ; i++ )    {
1028         for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
1029
1030           // check wether the coordinates were screwed
1031           if (TMath::Abs((coord(0,ip)*100-rlist[i]))>1 ||
1032               TMath::Abs((coord(1,ip)-philist[k]))>1 ||
1033               TMath::Abs((coord(2,ip)*100-zedlist[j]))>1) {
1034             AliError("internal error: coordinate system was screwed during the sum-up");
1035             printf("lookup: (r,phi,z)=(%f,%f,%f)\n",coord(0,ip)*100,coord(1,ip),coord(2,ip)*100);
1036             printf("sum-up: (r,phi,z)=(%f,%f,%f)\n",rlist[i],philist[k],zedlist[j]);
1037             AliError("Don't trust the results of the space charge calculation!");
1038           }
1039
1040           // unfortunately, the lookup tables were produced to be faster for phi symmetric charges
1041           // This will be the most frequent usage (hopefully)
1042           // That's why we have to do this here ...
1043
1044           TMatrixD &erA   =  *arrayofErA[k]  ;
1045           TMatrixD &ephiA =  *arrayofEphiA[k];
1046           TMatrixD &dEzA  =  *arrayofdEzA[k]   ;
1047
1048           TMatrixD &erC   =  *arrayofErC[k]  ;
1049           TMatrixD &ephiC =  *arrayofEphiC[k];
1050           TMatrixD &dEzC  =  *arrayofdEzC[k]   ;
1051
1052           // Sum up - Efield values in [V/m] -> transition to [V/cm]
1053           erA(i,j) += ((*bEr)(ip)) * weightA /100;
1054           erC(i,j) += ((*bEr)(ip)) * weightC /100;
1055           if (!flagRadSym) {
1056             ephiA(i,j) += ((*bEphi)(ip)) * weightA/100; // [V/rad]
1057             ephiC(i,j) += ((*bEphi)(ip)) * weightC/100; // [V/rad]
1058           }
1059           dEzA(i,j)  += ((*bEz)(ip)) * weightA /100;
1060           dEzC(i,j)  += ((*bEz)(ip)) * weightC /100;
1061
1062           // increase the counter
1063           ip++;
1064         }
1065       }
1066     } // end coordinate loop
1067
1068
1069     // Rotation and summation in the rest of the dPhiSteps
1070     // which were not stored in the this tree due to storage & symmetry reasons
1071
1072     Int_t phiPoints = (Int_t) grid(1);
1073     Int_t phiPOC    = (Int_t) grid(4);
1074
1075     //   printf("%d %d\n",phiPOC,flagRadSym);
1076
1077     for (Int_t phiiC = 1; phiiC<phiPOC; phiiC++) { // just used for non-radial symetric table
1078
1079       r0 = coordPOC(ipC,0);
1080       phi0 = coordPOC(ipC,1);
1081       z0 = coordPOC(ipC,2);
1082
1083       ipC++; // POC conf. counter
1084
1085       // weights (charge density) at POC position on the A and C side (in C/m^3/e0)
1086       // note: coordinates are in [cm]
1087       weightA = GetSpaceChargeDensity(r0*100,phi0, z0*100,0);
1088       weightC = GetSpaceChargeDensity(r0*100,phi0,-z0*100,0);
1089
1090       //     printf("%f %f %f: %e %e\n",r0,phi0,z0,weightA,weightC);
1091
1092       // Summing up the vector components according to their weight
1093       ip = 0;
1094       for ( Int_t j = 0 ; j < columns ; j++ ) {
1095         for ( Int_t i = 0 ; i < rows ; i++ )    {
1096           for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
1097
1098             // Note: rotating the coordinated during the sum up
1099
1100             Int_t rotVal = (phiPoints/phiPOC)*phiiC;
1101             Int_t ipR = -1;
1102
1103             if ((ip%phiPoints)>=rotVal) {
1104               ipR = ip-rotVal;
1105             } else {
1106               ipR = ip+(phiPoints-rotVal);
1107             }
1108
1109             // unfortunately, the lookup tables were produced to be faster for phi symmetric charges
1110             // This will be the most frequent usage
1111             // That's why we have to do this here and not outside the loop ...
1112
1113             TMatrixD &erA   =  *arrayofErA[k]  ;
1114             TMatrixD &ephiA =  *arrayofEphiA[k];
1115             TMatrixD &dEzA  =  *arrayofdEzA[k]   ;
1116
1117             TMatrixD &erC   =  *arrayofErC[k]  ;
1118             TMatrixD &ephiC =  *arrayofEphiC[k];
1119             TMatrixD &dEzC  =  *arrayofdEzC[k]   ;
1120
1121             // Sum up - Efield values in [V/m] -> transition to [V/cm]
1122             erA(i,j) += ((*bEr)(ipR)) * weightA /100;
1123             erC(i,j) += ((*bEr)(ipR)) * weightC /100;
1124             if (!flagRadSym) {
1125               ephiA(i,j) += ((*bEphi)(ipR)) * weightA/100; // [V/rad]
1126               ephiC(i,j) += ((*bEphi)(ipR)) * weightC/100; // [V/rad]
1127             }
1128             dEzA(i,j)  += ((*bEz)(ipR)) * weightA /100;
1129             dEzC(i,j)  += ((*bEz)(ipR)) * weightC /100;
1130
1131             // increase the counter
1132             ip++;
1133           }
1134         }
1135       } // end coordinate loop
1136
1137     } // end phi-POC summation (phiiC)
1138
1139
1140     // printf("POC: (r,phi,z) = (%f %f %f) | weight(A,C): %03.1lf %03.1lf\n",r0,phi0,z0, weightA, weightC);
1141
1142   }
1143
1144
1145
1146   // -------------------------------------------------------------------------------
1147   // Division by the Ez (drift) field and integration along z
1148
1149   AliInfo("Division and integration");
1150
1151   Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = Electric Field (V/cm) Magnitude ~ -400 V/cm;
1152
1153   for ( Int_t k = 0 ; k < phiSlices ; k++ ) { // phi loop
1154
1155     // matrices holding the solution - summation of POC charges // see above
1156     TMatrixD &erA   =  *arrayofErA[k]  ;
1157     TMatrixD &ephiA =  *arrayofEphiA[k];
1158     TMatrixD &dezA  =  *arrayofdEzA[k]   ;
1159     TMatrixD &erC   =  *arrayofErC[k]  ;
1160     TMatrixD &ephiC =  *arrayofEphiC[k];
1161     TMatrixD &dezC  =  *arrayofdEzC[k]   ;
1162
1163     // matrices which will contain the integrated fields (divided by the drift field)
1164     TMatrixD &erOverEzA   =  *arrayofEroverEzA[k]  ;
1165     TMatrixD &ephiOverEzA =  *arrayofEphioverEzA[k];
1166     TMatrixD &deltaEzA    =  *arrayofDeltaEzA[k];
1167     TMatrixD &erOverEzC   =  *arrayofEroverEzC[k]  ;
1168     TMatrixD &ephiOverEzC =  *arrayofEphioverEzC[k];
1169     TMatrixD &deltaEzC    =  *arrayofDeltaEzC[k];
1170
1171     for ( Int_t i = 0 ; i < rows ; i++ )    { // r loop
1172       for ( Int_t j = columns-1 ; j >= 0 ; j-- ) {// z loop
1173         // Count backwards to facilitate integration over Z
1174
1175         Int_t index = 1 ; // Simpsons rule if N=odd.If N!=odd then add extra point by trapezoidal rule.
1176
1177         erOverEzA(i,j) = 0; ephiOverEzA(i,j) = 0;  deltaEzA(i,j) = 0;
1178         erOverEzC(i,j) = 0; ephiOverEzC(i,j) = 0;  deltaEzC(i,j) = 0;
1179
1180         for ( Int_t m = j ; m < columns ; m++ ) { // integration
1181
1182           erOverEzA(i,j)   += index*(gridSizeZ/3.0)*erA(i,m)/(-1*ezField) ;
1183           erOverEzC(i,j)   += index*(gridSizeZ/3.0)*erC(i,m)/(-1*ezField)  ;
1184           if (!flagRadSym) {
1185             ephiOverEzA(i,j) += index*(gridSizeZ/3.0)*ephiA(i,m)/(-1*ezField)  ;
1186             ephiOverEzC(i,j) += index*(gridSizeZ/3.0)*ephiC(i,m)/(-1*ezField)  ;
1187           }
1188           deltaEzA(i,j)    += index*(gridSizeZ/3.0)*dezA(i,m)/(-1) ;
1189           deltaEzC(i,j)    += index*(gridSizeZ/3.0)*dezC(i,m)/(-1) ;
1190
1191           if ( index != 4 )  index = 4; else index = 2 ;
1192
1193         }
1194
1195         if ( index == 4 ) {
1196           erOverEzA(i,j)   -= (gridSizeZ/3.0)*erA(i,columns-1)/(-1*ezField) ;
1197           erOverEzC(i,j)   -= (gridSizeZ/3.0)*erC(i,columns-1)/(-1*ezField) ;
1198           if (!flagRadSym) {
1199             ephiOverEzA(i,j) -= (gridSizeZ/3.0)*ephiA(i,columns-1)/(-1*ezField) ;
1200             ephiOverEzC(i,j) -= (gridSizeZ/3.0)*ephiC(i,columns-1)/(-1*ezField) ;
1201           }
1202           deltaEzA(i,j)    -= (gridSizeZ/3.0)*dezA(i,columns-1)/(-1) ;
1203           deltaEzC(i,j)    -= (gridSizeZ/3.0)*dezC(i,columns-1)/(-1) ;
1204         }
1205         if ( index == 2 ) {
1206           erOverEzA(i,j)   += (gridSizeZ/3.0)*(0.5*erA(i,columns-2)-2.5*erA(i,columns-1))/(-1*ezField) ;
1207           erOverEzC(i,j)   += (gridSizeZ/3.0)*(0.5*erC(i,columns-2)-2.5*erC(i,columns-1))/(-1*ezField) ;
1208           if (!flagRadSym) {
1209             ephiOverEzA(i,j) += (gridSizeZ/3.0)*(0.5*ephiA(i,columns-2)-2.5*ephiA(i,columns-1))/(-1*ezField) ;
1210             ephiOverEzC(i,j) += (gridSizeZ/3.0)*(0.5*ephiC(i,columns-2)-2.5*ephiC(i,columns-1))/(-1*ezField) ;
1211           }
1212           deltaEzA(i,j)    += (gridSizeZ/3.0)*(0.5*dezA(i,columns-2)-2.5*dezA(i,columns-1))/(-1) ;
1213           deltaEzC(i,j)    += (gridSizeZ/3.0)*(0.5*dezC(i,columns-2)-2.5*dezC(i,columns-1))/(-1) ;
1214         }
1215         if ( j == columns-2 ) {
1216           erOverEzA(i,j)   = (gridSizeZ/3.0)*(1.5*erA(i,columns-2)+1.5*erA(i,columns-1))/(-1*ezField) ;
1217           erOverEzC(i,j)   = (gridSizeZ/3.0)*(1.5*erC(i,columns-2)+1.5*erC(i,columns-1))/(-1*ezField) ;
1218           if (!flagRadSym) {
1219             ephiOverEzA(i,j) = (gridSizeZ/3.0)*(1.5*ephiA(i,columns-2)+1.5*ephiA(i,columns-1))/(-1*ezField) ;
1220             ephiOverEzC(i,j) = (gridSizeZ/3.0)*(1.5*ephiC(i,columns-2)+1.5*ephiC(i,columns-1))/(-1*ezField) ;
1221           }
1222           deltaEzA(i,j)    = (gridSizeZ/3.0)*(1.5*dezA(i,columns-2)+1.5*dezA(i,columns-1))/(-1) ;
1223           deltaEzC(i,j)    = (gridSizeZ/3.0)*(1.5*dezC(i,columns-2)+1.5*dezC(i,columns-1))/(-1) ;
1224         }
1225         if ( j == columns-1 ) {
1226           erOverEzA(i,j)   = 0;
1227           erOverEzC(i,j)   = 0;
1228           if (!flagRadSym) {
1229             ephiOverEzA(i,j) = 0;
1230             ephiOverEzC(i,j) = 0;
1231           }
1232           deltaEzA(i,j)    = 0;
1233           deltaEzC(i,j)    = 0;
1234         }
1235       }
1236     }
1237
1238   }
1239
1240
1241
1242   AliInfo("Interpolation to Standard grid");
1243
1244   // -------------------------------------------------------------------------------
1245   // Interpolate results onto the standard grid which is used for all AliTPCCorrections classes
1246
1247   const Int_t order  = 1  ;  // Linear interpolation = 1, Quadratic = 2
1248
1249   Double_t  r, phi, z ;
1250   for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
1251     phi = fgkPhiList[k] ;
1252
1253     TMatrixF &erOverEz   =  *fLookUpErOverEz[k]  ;
1254     TMatrixF &ephiOverEz =  *fLookUpEphiOverEz[k];
1255     TMatrixF &deltaEz    =  *fLookUpDeltaEz[k]   ;
1256
1257     for ( Int_t j = 0 ; j < kNZ ; j++ ) {
1258
1259       z = TMath::Abs(fgkZList[j]) ;  // z position is symmetric
1260
1261       for ( Int_t i = 0 ; i < kNR ; i++ ) {
1262         r = fgkRList[i] ;
1263
1264         // Interpolate Lookup tables onto standard grid
1265         if (fgkZList[j]>0) {
1266           erOverEz(i,j)   = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
1267                                                rlist, zedlist, philist, arrayofEroverEzA  );
1268           ephiOverEz(i,j) = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
1269                                                rlist, zedlist, philist, arrayofEphioverEzA);
1270           deltaEz(i,j)    = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
1271                                                rlist, zedlist, philist, arrayofDeltaEzA   );
1272         } else {
1273           erOverEz(i,j)   = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
1274                                                rlist, zedlist, philist, arrayofEroverEzC  );
1275           ephiOverEz(i,j) = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
1276                                                rlist, zedlist, philist, arrayofEphioverEzC);
1277           deltaEz(i,j)  = - Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
1278                                                rlist, zedlist, philist, arrayofDeltaEzC   );
1279           // negative coordinate system on C side
1280         }
1281
1282       } // end r loop
1283     } // end z loop
1284   } // end phi loop
1285
1286
1287   // clear the temporary arrays lists
1288   for ( Int_t k = 0 ; k < phiSlices ; k++ )  {
1289
1290     delete arrayofErA[k];
1291     delete arrayofEphiA[k];
1292     delete arrayofdEzA[k];
1293     delete arrayofErC[k];
1294     delete arrayofEphiC[k];
1295     delete arrayofdEzC[k];
1296
1297     delete arrayofEroverEzA[k];
1298     delete arrayofEphioverEzA[k];
1299     delete arrayofDeltaEzA[k];
1300     delete arrayofEroverEzC[k];
1301     delete arrayofEphioverEzC[k];
1302     delete arrayofDeltaEzC[k];
1303
1304   }
1305
1306   fInitLookUp = kTRUE;
1307
1308 }
1309
1310
1311 void AliTPCSpaceCharge3D::SetSCDataFileName(TString fname) {
1312   /// Set & load the Space charge density distribution from a file
1313   /// (linear interpolation onto a standard grid)
1314
1315
1316   fSCDataFileName = fname;
1317
1318   TFile *f = new TFile(fSCDataFileName.Data(),"READ");
1319   if (!f) {
1320     AliError(Form("File %s, which should contain the space charge distribution, could not be found",
1321                   fSCDataFileName.Data()));
1322     return;
1323   }
1324
1325   TH2F *densityRZ = (TH2F*) f->Get("SpaceChargeInRZ");
1326   if (!densityRZ) {
1327     AliError(Form("The indicated file (%s) does not contain a histogram called %s",
1328                   fSCDataFileName.Data(),"SpaceChargeInRZ"));
1329     return;
1330   }
1331
1332   TH3F *densityRPhi = (TH3F*) f->Get("SpaceChargeInRPhi");
1333   if (!densityRPhi) {
1334     AliError(Form("The indicated file (%s) does not contain a histogram called %s",
1335                   fSCDataFileName.Data(),"SpaceChargeInRPhi"));
1336     return;
1337   }
1338
1339
1340   Double_t  r, phi, z ;
1341
1342   TMatrixD &scDensityInRZ   =  *fSCdensityInRZ;
1343   TMatrixD &scDensityInRPhiA   =  *fSCdensityInRPhiA;
1344   TMatrixD &scDensityInRPhiC   =  *fSCdensityInRPhiC;
1345   for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
1346     phi = fgkPhiList[k] ;
1347     TMatrixF &scDensity   =  *fSCdensityDistribution[k]  ;
1348     for ( Int_t j = 0 ; j < kNZ ; j++ ) {
1349       z = fgkZList[j] ;
1350       for ( Int_t i = 0 ; i < kNR ; i++ ) {
1351         r = fgkRList[i] ;
1352
1353         // partial load in (r,z)
1354         if (k==0) // do just once
1355           scDensityInRZ(i,j) =  densityRZ->Interpolate(r,z);
1356
1357         // partial load in (r,phi)
1358         if ( j==0 || j == kNZ/2 ) {
1359           if (z>0)
1360             scDensityInRPhiA(i,k) =  densityRPhi->Interpolate(r,phi,0.499);  // A side
1361           else
1362             scDensityInRPhiC(i,k) =  densityRPhi->Interpolate(r,phi,-0.499); // C side
1363         }
1364
1365         // Full 3D configuration ...
1366         if (z>0)
1367            scDensity(i,j) = scDensityInRZ(i,j) + scDensityInRPhiA(i,k);
1368         else
1369            scDensity(i,j) = scDensityInRZ(i,j) + scDensityInRPhiC(i,k);
1370       }
1371     }
1372   }
1373
1374   f->Close();
1375
1376   fInitLookUp = kFALSE;
1377
1378
1379 }
1380
1381 void     AliTPCSpaceCharge3D::SetInputSpaceCharge(TH3 * hisSpaceCharge3D,  TH2 * hisRPhi, TH2* hisRZ, Double_t norm){
1382   /// Use 3D space charge map as an optional input
1383   /// The layout of the input histogram is assumed to be: (phi,r,z)
1384   /// Density histogram is expreseed is expected to bin in  C/m^3
1385   ///
1386   /// Standard histogram interpolation is used in order to use the density at center of voxel
1387
1388   fSpaceChargeHistogram3D = hisSpaceCharge3D;
1389   fSpaceChargeHistogramRPhi = hisRPhi;
1390   fSpaceChargeHistogramRZ = hisRZ;
1391
1392   Double_t  r, phi, z ;
1393   TMatrixD &scDensityInRZ   =  *fSCdensityInRZ;
1394   TMatrixD &scDensityInRPhiA   =  *fSCdensityInRPhiA;
1395   TMatrixD &scDensityInRPhiC   =  *fSCdensityInRPhiC;
1396   //
1397   Double_t rmin=hisSpaceCharge3D->GetYaxis()->GetBinCenter(0);
1398   Double_t rmax=hisSpaceCharge3D->GetYaxis()->GetBinUpEdge(hisSpaceCharge3D->GetYaxis()->GetNbins());
1399   Double_t zmin=hisSpaceCharge3D->GetZaxis()->GetBinCenter(0);
1400   Double_t zmax=hisSpaceCharge3D->GetZaxis()->GetBinCenter(hisSpaceCharge3D->GetZaxis()->GetNbins());
1401
1402   for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
1403     phi = fgkPhiList[k] ;
1404     TMatrixF &scDensity   =  *fSCdensityDistribution[k]  ;
1405     for ( Int_t j = 0 ; j < kNZ ; j++ ) {
1406       z = fgkZList[j] ;
1407       for ( Int_t i = 0 ; i < kNR ; i++ ) {
1408         // Full 3D configuration ...
1409         r = fgkRList[i] ;
1410         if (r>rmin && r<rmax && z>zmin && z< zmax){
1411           // partial load in (r,z)
1412           if (k==0) {
1413             if (fSpaceChargeHistogramRZ) scDensityInRZ(i,j) = norm*fSpaceChargeHistogramRZ->Interpolate(r,z) ;
1414           }
1415           // partial load in (r,phi)
1416           if ( (j==0 || j == kNZ/2) && fSpaceChargeHistogramRPhi) {
1417             if (z>0)
1418               scDensityInRPhiA(i,k) = norm*fSpaceChargeHistogramRPhi->Interpolate(phi,r);  // A side
1419             else
1420               scDensityInRPhiC(i,k) = norm*fSpaceChargeHistogramRPhi->Interpolate(phi+TMath::TwoPi(),r); // C side
1421           }
1422
1423           if (z>0)
1424             scDensity(i,j) = norm*fSpaceChargeHistogram3D->Interpolate(phi,r,z);
1425           else
1426             scDensity(i,j) = norm*fSpaceChargeHistogram3D->Interpolate(phi,r,z);
1427         }
1428       }
1429     }
1430   }
1431
1432   fInitLookUp = kFALSE;
1433
1434 }
1435
1436
1437 Float_t  AliTPCSpaceCharge3D::GetSpaceChargeDensity(Float_t r, Float_t phi, Float_t z, Int_t mode) {
1438   /// returns the (input) space charge density at a given point according
1439   /// Note: input in [cm], output in [C/m^3/e0] !!
1440
1441   while (phi<0) phi += TMath::TwoPi();
1442   while (phi>TMath::TwoPi()) phi -= TMath::TwoPi();
1443
1444
1445   // Float_t sc =fSCdensityDistribution->Interpolate(r0,phi0,z0);
1446   Int_t order = 1; //
1447   Float_t sc = 0;
1448
1449   if (mode == 0) { // return full load
1450     sc = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
1451                             fgkRList, fgkZList, fgkPhiList, fSCdensityDistribution );
1452
1453   } else if (mode == 1) { // return partial load in (r,z)
1454     TMatrixD &scDensityInRZ   =  *fSCdensityInRZ;
1455     sc = Interpolate2DTable(order, r, z, kNR, kNZ, fgkRList, fgkZList, scDensityInRZ );
1456
1457   } else if (mode == 2) { // return partial load in (r,phi)
1458
1459     if (z>0) {
1460       TMatrixD &scDensityInRPhi   =  *fSCdensityInRPhiA;
1461       sc = Interpolate2DTable(order, r, phi, kNR, kNPhi, fgkRList, fgkPhiList, scDensityInRPhi );
1462     } else {
1463       TMatrixD &scDensityInRPhi   =  *fSCdensityInRPhiC;
1464       sc = Interpolate2DTable(order, r, phi, kNR, kNPhi, fgkRList, fgkPhiList, scDensityInRPhi );
1465     }
1466
1467   } else {
1468     // should i give a warning?
1469     sc = 0;
1470   }
1471
1472   //  printf("%f %f %f: %f\n",r,phi,z,sc);
1473
1474   return sc;
1475 }
1476
1477
1478 TH2F * AliTPCSpaceCharge3D::CreateHistoSCinXY(Float_t z, Int_t nx, Int_t ny, Int_t mode) {
1479   /// return a simple histogramm containing the space charge distribution (input for the calculation)
1480
1481   TH2F *h=CreateTH2F("spaceCharge",GetTitle(),"x [cm]","y [cm]","#rho_{sc} [C/m^{3}/e_{0}]",
1482                      nx,-250.,250.,ny,-250.,250.);
1483
1484   for (Int_t iy=1;iy<=ny;++iy) {
1485     Double_t yp = h->GetYaxis()->GetBinCenter(iy);
1486     for (Int_t ix=1;ix<=nx;++ix) {
1487       Double_t xp = h->GetXaxis()->GetBinCenter(ix);
1488
1489       Float_t r = TMath::Sqrt(xp*xp+yp*yp);
1490       Float_t phi = TMath::ATan2(yp,xp);
1491
1492       if (85.<=r && r<=250.) {
1493         Float_t sc = GetSpaceChargeDensity(r,phi,z,mode)/fgke0; // in [C/m^3/e0]
1494         h->SetBinContent(ix,iy,sc);
1495       } else {
1496         h->SetBinContent(ix,iy,0.);
1497       }
1498     }
1499   }
1500
1501   return h;
1502 }
1503
1504 TH2F * AliTPCSpaceCharge3D::CreateHistoSCinZR(Float_t phi, Int_t nz, Int_t nr,Int_t mode ) {
1505   /// return a simple histogramm containing the space charge distribution (input for the calculation)
1506
1507   TH2F *h=CreateTH2F("spaceCharge",GetTitle(),"z [cm]","r [cm]","#rho_{sc} [C/m^{3}/e_{0}]",
1508                      nz,-250.,250.,nr,85.,250.);
1509
1510   for (Int_t ir=1;ir<=nr;++ir) {
1511     Float_t r = h->GetYaxis()->GetBinCenter(ir);
1512     for (Int_t iz=1;iz<=nz;++iz) {
1513       Float_t z = h->GetXaxis()->GetBinCenter(iz);
1514       Float_t sc = GetSpaceChargeDensity(r,phi,z,mode)/fgke0; // in [C/m^3/e0]
1515       h->SetBinContent(iz,ir,sc);
1516     }
1517   }
1518
1519   return h;
1520 }
1521
1522 void AliTPCSpaceCharge3D::WriteChargeDistributionToFile(const char* fname) {
1523   /// Example on how to write a Space charge distribution into a File
1524   ///  (see below: estimate from scaling STAR measurements to Alice)
1525   /// Charge distribution is splitted into two (RZ and RPHI) in order to speed up
1526   /// the needed calculation time
1527
1528   TFile *f = new TFile(fname,"RECREATE");
1529
1530   // some grid, not too course
1531   Int_t nr = 350;
1532   Int_t nphi = 180;
1533   Int_t nz = 500;
1534
1535   Double_t dr = (fgkOFCRadius-fgkIFCRadius)/(nr+1);
1536   Double_t dphi = TMath::TwoPi()/(nphi+1);
1537   Double_t dz = 500./(nz+1);
1538   Double_t safty = 0.; // due to a root bug which does not interpolate the boundary (first and last bin) correctly
1539
1540
1541   // Charge distribution in ZR (rotational symmetric) ------------------
1542
1543   TH2F *histoZR = new TH2F("chargeZR","chargeZR",
1544                            nr,fgkIFCRadius-dr-safty,fgkOFCRadius+dr+safty,
1545                            nz,-250-dz-safty,250+dz+safty);
1546
1547   for (Int_t ir=1;ir<=nr;++ir) {
1548     Double_t rp = histoZR->GetXaxis()->GetBinCenter(ir);
1549     for (Int_t iz=1;iz<=nz;++iz) {
1550       Double_t zp = histoZR->GetYaxis()->GetBinCenter(iz);
1551
1552       // recalculation to meter
1553       Double_t lZ = 2.5; // approx. TPC drift length
1554       Double_t rpM = rp/100.; // in [m]
1555       Double_t zpM = TMath::Abs(zp/100.); // in [m]
1556
1557       // setting of mb multiplicity and Interaction rate
1558       Double_t multiplicity = 950;
1559       Double_t intRate = 7800;
1560
1561       // calculation of "scaled" parameters
1562       Double_t a = multiplicity*intRate/79175;
1563       Double_t b = a/lZ;
1564
1565       Double_t charge = (a - b*zpM)/(rpM*rpM); // charge in [C/m^3/e0]
1566
1567       charge = charge*fgke0; // [C/m^3]
1568
1569       if (zp<0) charge *= 0.9; // e.g. slightly less on C side due to front absorber
1570
1571       //  charge = 0; // for tests
1572       histoZR->SetBinContent(ir,iz,charge);
1573     }
1574   }
1575
1576   histoZR->Write("SpaceChargeInRZ");
1577
1578
1579   // Charge distribution in RPhi (e.g. Floating GG wire) ------------
1580
1581   TH3F *histoRPhi = new TH3F("chargeRPhi","chargeRPhi",
1582                              nr,fgkIFCRadius-dr-safty,fgkOFCRadius+dr+safty,
1583                              nphi,0-dphi-safty,TMath::TwoPi()+dphi+safty,
1584                              2,-1,1); // z part - to allow A and C side differences
1585
1586   // some 'arbitrary' GG leaks
1587   Int_t   nGGleaks = 5;
1588   Double_t secPosA[5]    = {3,6,6,11,13};         // sector
1589   Double_t radialPosA[5] = {125,100,160,200,230}; // radius in cm
1590   Double_t secPosC[5]    = {1,8,12,15,15};        // sector
1591   Double_t radialPosC[5] = {245,120,140,120,190}; // radius in cm
1592
1593   for (Int_t ir=1;ir<=nr;++ir) {
1594     Double_t rp = histoRPhi->GetXaxis()->GetBinCenter(ir);
1595     for (Int_t iphi=1;iphi<=nphi;++iphi) {
1596       Double_t phip = histoRPhi->GetYaxis()->GetBinCenter(iphi);
1597       for (Int_t iz=1;iz<=2;++iz) {
1598         Double_t zp = histoRPhi->GetZaxis()->GetBinCenter(iz);
1599
1600         Double_t charge = 0;
1601
1602         for (Int_t igg = 0; igg<nGGleaks; igg++) { // loop over GG leaks
1603
1604           // A side
1605           Double_t secPos = secPosA[igg];
1606           Double_t radialPos = radialPosA[igg];
1607
1608           if (zp<0) { // C side
1609             secPos = secPosC[igg];
1610             radialPos = radialPosC[igg];
1611           }
1612
1613           // some 'arbitrary' GG leaks
1614           if (  (phip<(TMath::Pi()/9*(secPos+1)) && phip>(TMath::Pi()/9*secPos) ) ) { // sector slice
1615             if ( rp>(radialPos-2.5) && rp<(radialPos+2.5))  // 5 cm slice
1616               charge = 300;
1617           }
1618
1619         }
1620
1621         charge = charge*fgke0; // [C/m^3]
1622
1623         histoRPhi->SetBinContent(ir,iphi,iz,charge);
1624       }
1625     }
1626   }
1627
1628   histoRPhi->Write("SpaceChargeInRPhi");
1629
1630   f->Close();
1631
1632 }
1633
1634
1635 void AliTPCSpaceCharge3D::Print(const Option_t* option) const {
1636   /// Print function to check the settings of the boundary vectors
1637   /// option=="a" prints the C0 and C1 coefficents for calibration purposes
1638
1639   TString opt = option; opt.ToLower();
1640   printf("%s\n",GetTitle());
1641   printf(" - Space Charge effect with arbitrary 3D charge density (as input).\n");
1642   printf("   SC correction factor: %f \n",fCorrectionFactor);
1643
1644   if (opt.Contains("a")) { // Print all details
1645     printf(" - T1: %1.4f, T2: %1.4f \n",fT1,fT2);
1646     printf(" - C1: %1.4f, C0: %1.4f \n",fC1,fC0);
1647   }
1648
1649   if (!fInitLookUp) AliError("Lookup table was not initialized! You should do InitSpaceCharge3DDistortion() ...");
1650
1651 }
1652
1653
1654
1655 void AliTPCSpaceCharge3D::InitSpaceCharge3DPoisson(Int_t kRows, Int_t kColumns, Int_t kPhiSlices, Int_t kIterations){
1656   /// MI extension  - calculate E field
1657   /// - inspired by  AliTPCROCVoltError3D::InitROCVoltError3D()
1658   /// Initialization of the Lookup table which contains the solutions of the
1659   /// Dirichlet boundary problem
1660   /// Calculation of the single 3D-Poisson solver is done just if needed
1661   /// (see basic lookup tables in header file)
1662
1663   Int_t kPhiSlicesPerSector = kPhiSlices/18;
1664   //
1665   const Int_t   order       =    1  ;  // Linear interpolation = 1, Quadratic = 2
1666   const Float_t gridSizeR   =  (fgkOFCRadius-fgkIFCRadius) / (kRows-1) ;
1667   const Float_t gridSizeZ   =  fgkTPCZ0 / (kColumns-1) ;
1668   const Float_t gridSizePhi =  TMath::TwoPi() / ( 18.0 * kPhiSlicesPerSector);
1669
1670   // temporary arrays to create the boundary conditions
1671   TMatrixD *arrayofArrayV[kPhiSlices], *arrayofCharge[kPhiSlices] ;
1672   TMatrixD *arrayofEroverEz[kPhiSlices], *arrayofEphioverEz[kPhiSlices], *arrayofDeltaEz[kPhiSlices] ;
1673
1674   for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
1675     arrayofArrayV[k]     =   new TMatrixD(kRows,kColumns) ;
1676     arrayofCharge[k]     =   new TMatrixD(kRows,kColumns) ;
1677     arrayofEroverEz[k]   =   new TMatrixD(kRows,kColumns) ;
1678     arrayofEphioverEz[k] =   new TMatrixD(kRows,kColumns) ;
1679     arrayofDeltaEz[k]    =   new TMatrixD(kRows,kColumns) ;
1680   }
1681
1682   // list of point as used in the poisson relation and the interpolation (during sum up)
1683   Double_t  rlist[kRows], zedlist[kColumns] , philist[kPhiSlices];
1684   for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
1685     philist[k] =  gridSizePhi * k;
1686     for ( Int_t i = 0 ; i < kRows ; i++ )    {
1687       rlist[i] = fgkIFCRadius + i*gridSizeR ;
1688       for ( Int_t j = 0 ; j < kColumns ; j++ ) { // Fill Vmatrix with Boundary Conditions
1689         zedlist[j]  = j * gridSizeZ ;
1690       }
1691     }
1692   }
1693
1694   // ==========================================================================
1695   // Solve Poisson's equation in 3D cylindrical coordinates by relaxation technique
1696   // Allow for different size grid spacing in R and Z directions
1697
1698   const Int_t   symmetry = 0;
1699
1700   // Set bondaries and solve Poisson's equation --------------------------
1701
1702   if ( !fInitLookUp ) {
1703
1704     AliInfo(Form("Solving the poisson equation (~ %d sec)",2*10*(int)(kPhiSlices/10)));
1705
1706     for ( Int_t side = 0 ; side < 2 ; side++ ) {  // Solve Poisson3D twice; once for +Z and once for -Z
1707       AliSysInfo::AddStamp("RunSide", 1,side,0);
1708       for ( Int_t k = 0 ; k < kPhiSlices ; k++ )  {
1709         TMatrixD &arrayV    =  *arrayofArrayV[k] ;
1710         TMatrixD &charge    =  *arrayofCharge[k] ;
1711
1712         //Fill arrays with initial conditions.  V on the boundary and Charge in the volume.
1713 //      for ( Int_t i = 0 ; i < kRows ; i++ ) {
1714 //        for ( Int_t j = 0 ; j < kColumns ; j++ ) {  // Fill Vmatrix with Boundary Conditions
1715 //          arrayV(i,j) = 0.0 ;
1716 //          charge(i,j) = 0.0 ;
1717
1718 // //       Float_t radius0 = rlist[i] ;
1719 // //       Float_t phi0    = gridSizePhi * k ;
1720
1721 //          // To avoid problems at sector boundaries, use an average of +- 1 degree from actual phi location
1722 // //       if ( j == (kColumns-1) ) {
1723 // //         arrayV(i,j) = 0.5*  ( GetROCVoltOffset( side, radius0, phi0+0.02 ) + GetROCVoltOffset( side, radius0, phi0-0.02 ) ) ;
1724
1725 // //         if (side==1) // C side
1726 // //           arrayV(i,j) = -arrayV(i,j); // minus sign on the C side to allow a consistent usage of global z when setting the boundaries
1727 // //       }
1728 //        }
1729 //      }
1730
1731         for ( Int_t i = 1 ; i < kRows-1 ; i++ ) {
1732           for ( Int_t j = 1 ; j < kColumns-1 ; j++ ) {
1733             Float_t radius0 = rlist[i] ;
1734             Float_t phi0    = gridSizePhi * k ;
1735             Double_t z0 = zedlist[j];
1736             if (side==1) z0= -TMath::Abs(zedlist[j]);
1737             arrayV(i,j) = 0.0 ;
1738             charge(i,j)  =  fSpaceChargeHistogram3D->Interpolate(phi0,radius0,z0);
1739           }
1740         }
1741       }
1742       AliSysInfo::AddStamp("RunPoisson", 2,side,0);
1743
1744       // Solve Poisson's equation in 3D cylindrical coordinates by relaxation technique
1745       // Allow for different size grid spacing in R and Z directions
1746
1747       //      PoissonRelaxation3D( arrayofArrayV, arrayofCharge,
1748       //                           arrayofEroverEz, arrayofEphioverEz, arrayofDeltaEz,
1749       //                           kRows, kColumns, kPhiSlices, gridSizePhi, kIterations,
1750       //                           symmetry , fROCdisplacement) ;
1751       PoissonRelaxation3D( arrayofArrayV, arrayofCharge,
1752                            arrayofEroverEz, arrayofEphioverEz, arrayofDeltaEz,
1753                            kRows, kColumns, kPhiSlices, gridSizePhi, kIterations,
1754                            symmetry ) ;
1755
1756       //Interpolate results onto a custom grid which is used just for these calculations.
1757       Double_t  r, phi, z ;
1758       for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
1759         phi = fgkPhiList[k] ;
1760
1761         TMatrixF &erOverEz   =  *fLookUpErOverEz[k]  ;
1762         TMatrixF &ephiOverEz =  *fLookUpEphiOverEz[k];
1763         TMatrixF &deltaEz    =  *fLookUpDeltaEz[k]   ;
1764
1765         for ( Int_t j = 0 ; j < kNZ ; j++ ) {
1766
1767           z = TMath::Abs(fgkZList[j]) ;  // Symmetric solution in Z that depends only on ABS(Z)
1768
1769           if ( side == 0 &&  fgkZList[j] < 0 ) continue; // Skip rest of this loop if on the wrong side
1770           if ( side == 1 &&  fgkZList[j] > 0 ) continue; // Skip rest of this loop if on the wrong side
1771
1772           for ( Int_t i = 0 ; i < kNR ; i++ ) {
1773             r = fgkRList[i] ;
1774
1775             // Interpolate basicLookup tables; once for each rod, then sum the results
1776             erOverEz(i,j)   = Interpolate3DTable(order, r, z, phi, kRows, kColumns, kPhiSlices,
1777                                                  rlist, zedlist, philist, arrayofEroverEz  );
1778             ephiOverEz(i,j) = Interpolate3DTable(order, r, z, phi, kRows, kColumns, kPhiSlices,
1779                                                  rlist, zedlist, philist, arrayofEphioverEz);
1780             deltaEz(i,j)    = Interpolate3DTable(order, r, z, phi, kRows, kColumns, kPhiSlices,
1781                                                  rlist, zedlist, philist, arrayofDeltaEz  );
1782
1783             if (side == 1)  deltaEz(i,j) = -  deltaEz(i,j); // negative coordinate system on C side
1784
1785           } // end r loop
1786         }// end z loop
1787       }// end phi loop
1788       AliSysInfo::AddStamp("Interpolate Poisson", 3,side,0);
1789       if ( side == 0 ) AliInfo(" A side done");
1790       if ( side == 1 ) AliInfo(" C side done");
1791     } // end side loop
1792   }
1793
1794   // clear the temporary arrays lists
1795   for ( Int_t k = 0 ; k < kPhiSlices ; k++ )  {
1796     delete arrayofArrayV[k];
1797     delete arrayofCharge[k];
1798     delete arrayofEroverEz[k];
1799     delete arrayofEphioverEz[k];
1800     delete arrayofDeltaEz[k];
1801   }
1802
1803
1804   fInitLookUp = kTRUE;
1805
1806 }
1807
1808
1809
1810 AliTPCSpaceCharge3D * AliTPCSpaceCharge3D::MakeCorrection22D(const char* fileName, Double_t multiplicity, Double_t intRate, Double_t epsIROC, Double_t epsOROC,
1811                                                              Double_t gasfactor,
1812                                                              Double_t radialScaling){
1813   /// Origin: Christian Lippmann, CERN, Christian.Lippmann@cern.ch based on the internal note (xxx ...)
1814   /// adopted by Marian Ivanov (different epsilon in IROC and OROC)
1815   ///
1816   /// Charge distribution is splitted into two (RZ and RPHI) in order to speed up
1817   /// the needed calculation time.
1818   ///
1819   /// Explanation of variables:
1820   /// 1) multiplicity: charghed particle dn/deta for top 80% centrality (660 for 2011,
1821   ///    expect 950 for full energy)
1822   /// 2) intRate: Total interaction rate (e.g. 50kHz for the upgrade)
1823   /// 3) eps: Number of backdrifting ions per primary electron (0 for MWPC, e.g.10 for GEM)
1824   /// 4) gasfactor: Use different gas. E.g. Ar/CO2 has twice the primary ionization, ion drift
1825   ///    velocity factor 2.5 slower, so  gasfactor = 5.
1826
1827
1828   TFile *f = new TFile(fileName, "RECREATE");
1829   // some grid, not too coarse
1830   const Int_t nr   = 350;
1831   const Int_t nphi = 180;
1832   const Int_t nz   = 500;
1833   const Double_t kROROC=134;  // hardwired OROC radius
1834
1835   Double_t dr = (fgkOFCRadius-fgkIFCRadius)/(nr+1);
1836   Double_t dphi = TMath::TwoPi()/(nphi+1);
1837   Double_t dz = 500./(nz+1);
1838   Double_t safty = 0.; // due to a root bug which does not interpolate the boundary ..
1839   // .. (first and last bin) correctly
1840
1841   // Charge distribution in ZR (rotational symmetric) ------------------
1842
1843   TH2F *histoZR = new TH2F("chargeZR", "chargeZR",
1844                            nr, fgkIFCRadius-dr-safty, fgkOFCRadius+dr+safty,
1845                            nz, -250-dz-safty, 250+dz+safty);
1846
1847   // For the normalization to same integral as radial exponent = 2
1848   Double_t radialExponent             = -2.; // reference = 2
1849   Double_t radiusInner                = histoZR->GetXaxis()->GetBinCenter(1) / 100.;//in [m]
1850   Double_t radiusOuter                = histoZR->GetXaxis()->GetBinCenter(nr) / 100.;//in [m]
1851   Double_t integralRadialExponent2    = TMath::Power(radiusOuter,radialExponent+1) * 1./(radialExponent+1)
1852     - TMath::Power(radiusInner,radialExponent+1) * 1./(radialExponent+1);
1853
1854   radialExponent                      = -radialScaling; // user set
1855   Double_t integralRadialExponentUser = 0.;
1856   if(radialScaling > 1 + 0.000001 || radialScaling < 1 - 0.000001 ) // to avoid n = -1
1857     integralRadialExponentUser = TMath::Power(radiusOuter,radialExponent+1) * 1./(radialExponent+1)
1858       - TMath::Power(radiusInner,radialExponent+1) * 1./(radialExponent+1);
1859   else
1860     integralRadialExponentUser = TMath::Log(radiusOuter) - TMath::Log(radiusInner);
1861
1862   Double_t normRadialExponent         = integralRadialExponent2 / integralRadialExponentUser;
1863
1864   for (Int_t ir=1;ir<=nr;++ir) {
1865     Double_t rp = histoZR->GetXaxis()->GetBinCenter(ir);
1866     for (Int_t iz=1;iz<=nz;++iz) {
1867       Double_t zp = histoZR->GetYaxis()->GetBinCenter(iz);
1868       Double_t eps = (rp <kROROC) ? epsIROC:epsOROC;
1869       // recalculation to meter
1870       Double_t lZ = 2.5; // approx. TPC drift length
1871       Double_t rpM = rp/100.; // in [m]
1872       Double_t zpM = TMath::Abs(zp/100.); // in [m]
1873
1874       // calculation of "scaled" parameters
1875       Double_t a = multiplicity*intRate/76628;
1876       //Double_t charge = gasfactor * ( a / (rpM*rpM) * (1 - zpM/lZ) ); // charge in [C/m^3/e0], no IBF
1877       Double_t charge = normRadialExponent * gasfactor * ( a / (TMath::Power(rpM,radialScaling)) * (1 - zpM/lZ + eps) ); // charge in [C/m^3/e0], with IBF
1878
1879       charge = charge*fgke0;          // [C/m^3]
1880       if (zp<0) charge *= 0.9; // Slightly less on C side due to front absorber
1881
1882       histoZR->SetBinContent(ir, iz, charge);
1883     }
1884   }
1885
1886   histoZR->Write("SpaceChargeInRZ");
1887   //
1888   // Charge distribution in RPhi (e.g. Floating GG wire) ------------
1889   //
1890   TH3F *histoRPhi = new TH3F("chargeRPhi", "chargeRPhi",
1891                              nr, fgkIFCRadius-dr-safty, fgkOFCRadius+dr+safty,
1892                              nphi, 0-dphi-safty, TMath::TwoPi()+dphi+safty,
1893                              2, -1, 1); // z part - to allow A and C side differences
1894   histoRPhi->Write("SpaceChargeInRPhi");
1895   f->Close();
1896   //
1897   //
1898   //
1899   AliTPCSpaceCharge3D *spaceCharge = new AliTPCSpaceCharge3D;
1900   spaceCharge->SetSCDataFileName(fileName);
1901   spaceCharge->SetOmegaTauT1T2(0.325,1,1); // Ne CO2
1902   spaceCharge->InitSpaceCharge3DDistortion();
1903   spaceCharge->AddVisualCorrection(spaceCharge,1);
1904   //
1905   f = new TFile(fileName, "update");
1906   spaceCharge->Write("spaceCharge");
1907   f->Close();
1908   return spaceCharge;
1909 }
1910