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