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