AliTPCSpaceCharge3D: fixed wrong calib maps paths
[u/mrichter/AliRoot.git] / TPC / TPCbase / AliTPCSpaceCharge3D.cxx
CommitLineData
cc3e558a 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
7d855b04 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
b8e57390 44/// ![Picture from ROOT macro](AliTPCSpaceCharge3D_cxx_2829f39.png)
45///
46/// \author Stefan Rossegger
47/// \date 19/06/2010
b4caed64 48
cc3e558a 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"
92a85338 65#include "AliSysInfo.h"
cc3e558a 66
7d855b04 67/// \cond CLASSIMP
cc3e558a 68ClassImp(AliTPCSpaceCharge3D)
7d855b04 69/// \endcond
cc3e558a 70
71AliTPCSpaceCharge3D::AliTPCSpaceCharge3D()
72 : AliTPCCorrection("SpaceCharge3D","Space Charge - 3D"),
73 fC0(0.),fC1(0.),
74 fCorrectionFactor(1.),
75 fInitLookUp(kFALSE),
15687d71 76 fSCDataFileName(""),
77 fSCLookUpPOCsFileName3D(""),
78 fSCLookUpPOCsFileNameRZ(""),
79 fSCLookUpPOCsFileNameRPhi(""),
80 fSCdensityInRZ(0),
7d855b04 81 fSCdensityInRPhiA(0),
92a85338 82 fSCdensityInRPhiC(0),
83 fSpaceChargeHistogram3D(0),
84 fSpaceChargeHistogramRPhi(0),
85 fSpaceChargeHistogramRZ(0)
cc3e558a 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++ ) {
7d855b04 94 fLookUpErOverEz[k] = new TMatrixF(kNR,kNZ);
2bf29b72 95 fLookUpEphiOverEz[k] = new TMatrixF(kNR,kNZ);
7d855b04 96 fLookUpDeltaEz[k] = new TMatrixF(kNR,kNZ);
2bf29b72 97 fSCdensityDistribution[k] = new TMatrixF(kNR,kNZ);
cc3e558a 98 }
15687d71 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
9bcafb19 105 fSCLookUpPOCsFileName3D="$(ALICE_ROOT)/TPC/TPCcalib/maps/sc_3D_raw_18-18-26_17p-18p-25p-MN30.root"; // rough estimate
106 fSCLookUpPOCsFileNameRZ="$(ALICE_ROOT)/TPC/TPCcalib/maps/sc_radSym_35-01-51_34p-01p-50p_MN60.root";
107 fSCLookUpPOCsFileNameRPhi="$(ALICE_ROOT)/TPC/TPCcalib/maps/sc_cChInZ_35-144-26_34p-18p-01p-MN30.root";
752b0cc7 108 // fSCLookUpPOCsFileNameRPhi="$(ALICE_ROOT)/TPC/Calib/maps/sc_cChInZ_35-36-26_34p-18p-01p-MN40.root";
7d855b04 109
15687d71 110
111
112 // standard location of the space charge distibution ... can be changes
9bcafb19 113 fSCDataFileName="$(ALICE_ROOT)/TPC/TPCcalib/maps/sc_3D_distribution_Sim.root";
15687d71 114
115 // SetSCDataFileName(fSCDataFileName.Data()); // should be done by the user
cc3e558a 116
cc3e558a 117
118}
119
120AliTPCSpaceCharge3D::~AliTPCSpaceCharge3D() {
7d855b04 121 /// default destructor
122
cc3e558a 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 }
15687d71 129 delete fSCdensityInRZ;
130 delete fSCdensityInRPhiA;
131 delete fSCdensityInRPhiC;
92a85338 132 delete fSpaceChargeHistogram3D;
133 delete fSpaceChargeHistogramRPhi;
134 delete fSpaceChargeHistogramRZ;
cc3e558a 135}
136
137
cc3e558a 138void AliTPCSpaceCharge3D::Init() {
7d855b04 139 /// Initialization funtion
140
cc3e558a 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)
7d855b04 148 Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
cc3e558a 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
155void AliTPCSpaceCharge3D::Update(const TTimeStamp &/*timeStamp*/) {
7d855b04 156 /// Update function
157
cc3e558a 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)
7d855b04 165 Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
cc3e558a 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
cc3e558a 174void AliTPCSpaceCharge3D::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {
7d855b04 175 /// Calculates the correction due the Space Charge effect within the TPC drift volume
cc3e558a 176
177 if (!fInitLookUp) {
178 AliInfo("Lookup table was not initialized! Performing the inizialisation now ...");
179 InitSpaceCharge3DDistortion();
180 }
181
7d855b04 182 Int_t order = 1 ; // FIXME: hardcoded? Linear interpolation = 1, Quadratic = 2
183
2bf29b72 184 Float_t intEr, intEphi, intdEz ;
cc3e558a 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 }
7d855b04 198
cc3e558a 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
7d855b04 201
cc3e558a 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
7d855b04 206 // Get the Er and Ephi field integrals plus the integral over DeltaEz
207 intEr = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
cc3e558a 208 fgkRList, fgkZList, fgkPhiList, fLookUpErOverEz );
7d855b04 209 intEphi = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
cc3e558a 210 fgkRList, fgkZList, fgkPhiList, fLookUpEphiOverEz);
7d855b04 211 intdEz = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
cc3e558a 212 fgkRList, fgkZList, fgkPhiList, fLookUpDeltaEz );
213
214 // Calculate distorted position
215 if ( r > 0.0 ) {
7d855b04 216 phi = phi + fCorrectionFactor *( fC0*intEphi - fC1*intEr ) / r;
217 r = r + fCorrectionFactor *( fC0*intEr + fC1*intEphi );
cc3e558a 218 }
219 Double_t dz = intdEz * fCorrectionFactor * fgkdvdE;
7d855b04 220
cc3e558a 221 // Calculate correction in cartesian coordinates
752b0cc7 222 dx[0] = - (r * TMath::Cos(phi) - x[0]);
7d855b04 223 dx[1] = - (r * TMath::Sin(phi) - x[1]);
752b0cc7 224 dx[2] = - dz; // z distortion - (scaled with driftvelocity dependency on the Ez field and the overall scaling factor)
cc3e558a 225
226}
227
cc3e558a 228void AliTPCSpaceCharge3D::InitSpaceCharge3DDistortion() {
7d855b04 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
15687d71 235
236 if (fInitLookUp) {
237 AliInfo("Lookup table was already initialized! Doing it again anyway ...");
db59c7fb 238 return;
15687d71 239 }
7d855b04 240
15687d71 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.");
7d855b04 245
15687d71 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;
7d855b04 252 }
15687d71 253
254 // units are in [m]
7d855b04 255 TVector *gridf1 = (TVector*) fZR->Get("constants");
15687d71 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;
7d855b04 261
15687d71 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
7d855b04 270
271 TMatrixD *arrayofErA[kNPhiSlices], *arrayofdEzA[kNPhiSlices];
272 TMatrixD *arrayofErC[kNPhiSlices], *arrayofdEzC[kNPhiSlices];
15687d71 273
9f3b99e2 274 TMatrixD *arrayofEroverEzA[kNPhiSlices], *arrayofDeltaEzA[kNPhiSlices];
275 TMatrixD *arrayofEroverEzC[kNPhiSlices], *arrayofDeltaEzC[kNPhiSlices];
15687d71 276
7d855b04 277
15687d71 278 for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
7d855b04 279
15687d71 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
7d855b04 290 // zero initialization not necessary, it is done in the constructor of TMatrix
15687d71 291
292 }
7d855b04 293
15687d71 294 // list of points as used during sum up
9f3b99e2 295 Double_t rlist1[kNRows], zedlist1[kNColumns];// , philist1[phiSlices];
15687d71 296 for ( Int_t i = 0 ; i < rows ; i++ ) {
297 rlist1[i] = fgkIFCRadius + i*gridSizeR ;
7d855b04 298 for ( Int_t j = 0 ; j < columns ; j++ ) {
15687d71 299 zedlist1[j] = j * gridSizeZ ;
300 }
301 }
7d855b04 302
15687d71 303 TTree *treePOC = (TTree*)fZR->Get("POCall");
304
305 TVector *bEr = 0; //TVector *bEphi= 0;
306 TVector *bEz = 0;
7d855b04 307
15687d71 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 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
7d855b04 314
15687d71 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
7d855b04 319
15687d71 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
7d855b04 333
15687d71 334 // Summing up the vector components according to their weight
335
336 Int_t ip = 0;
7d855b04 337 for ( Int_t j = 0 ; j < columns ; j++ ) {
15687d71 338 for ( Int_t i = 0 ; i < rows ; i++ ) {
339 for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
7d855b04 340
15687d71 341 // check wether the coordinates were screwed
7d855b04 342 if (TMath::Abs((coord1(0,ip)*100-rlist1[i]))>1 ||
343 TMath::Abs((coord1(2,ip)*100-zedlist1[j])>1)) {
15687d71 344 AliError("internal error: coordinate system was screwed during the sum-up");
9f3b99e2 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);
15687d71 347 AliError("Don't trust the results of the space charge calculation!");
348 }
7d855b04 349
15687d71 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] ;
7d855b04 356
15687d71 357 TMatrixD &erC = *arrayofErC[k] ;
358 TMatrixD &dEzC = *arrayofdEzC[k] ;
7d855b04 359
15687d71 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] ;
7d855b04 393 TMatrixD &deltaEzC = *arrayofDeltaEzC[k];
394
15687d71 395 for ( Int_t i = 0 ; i < rows ; i++ ) { // r loop
7d855b04 396 for ( Int_t j = columns-1 ; j >= 0 ; j-- ) {// z loop
15687d71 397 // Count backwards to facilitate integration over Z
398
7d855b04 399 Int_t index = 1 ; // Simpsons rule if N=odd.If N!=odd then add extra point
400 // by trapezoidal rule.
15687d71 401
402 erOverEzA(i,j) = 0; //ephiOverEzA(i,j) = 0;
403 deltaEzA(i,j) = 0;
7d855b04 404 erOverEzC(i,j) = 0; //ephiOverEzC(i,j) = 0;
15687d71 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 ) {
7d855b04 437 erOverEzA(i,j) = 0;
15687d71 438 erOverEzC(i,j) = 0;
7d855b04 439 deltaEzA(i,j) = 0;
15687d71 440 deltaEzC(i,j) = 0;
441 }
442 }
443 }
444
445 }
7d855b04 446
15687d71 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
7d855b04 452 const Int_t order = 1 ; // Linear interpolation = 1, Quadratic = 2
15687d71 453
454 Double_t r, z;//phi, z ;
455 for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
456 // phi = fgkPhiList[k] ;
7d855b04 457
15687d71 458 // final lookup table
2bf29b72 459 TMatrixF &erOverEzFinal = *fLookUpErOverEz[k] ;
460 TMatrixF &deltaEzFinal = *fLookUpDeltaEz[k] ;
7d855b04 461
15687d71 462 // calculated and integrated tables - just one phi slice
463 TMatrixD &erOverEzA = *arrayofEroverEzA[0] ;
464 TMatrixD &deltaEzA = *arrayofDeltaEzA[0];
465 TMatrixD &erOverEzC = *arrayofEroverEzC[0] ;
7d855b04 466 TMatrixD &deltaEzC = *arrayofDeltaEzC[0];
467
468
15687d71 469 for ( Int_t j = 0 ; j < kNZ ; j++ ) {
470
471 z = TMath::Abs(fgkZList[j]) ; // z position is symmetric
7d855b04 472
473 for ( Int_t i = 0 ; i < kNR ; i++ ) {
15687d71 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
7d855b04 490
15687d71 491 // clear the temporary arrays lists
492 for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
493
7d855b04 494 delete arrayofErA[k];
15687d71 495 delete arrayofdEzA[k];
7d855b04 496 delete arrayofErC[k];
15687d71 497 delete arrayofdEzC[k];
498
7d855b04 499 delete arrayofEroverEzA[k];
15687d71 500 delete arrayofDeltaEzA[k];
7d855b04 501 delete arrayofEroverEzC[k];
15687d71 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
7d855b04 510
15687d71 511 // AliInfo("Step 2: Preparation of the weighted look-up table");
7d855b04 512
15687d71 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;
7d855b04 517 }
15687d71 518
519 // units are in [m]
7d855b04 520 TVector *gridf2 = (TVector*) fRPhi->Get("constants");
15687d71 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
7d855b04 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
15687d71 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]
7d855b04 534
15687d71 535 // list of points as used during sum up
7d855b04 536 Double_t rlist2[kNRows], philist2[kNPhiSlices], zedlist2[kNColumns];
15687d71 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 ;
7d855b04 541 for ( Int_t j = 0 ; j < columns ; j++ ) {
15687d71 542 zedlist2[j] = j * gridSizeZ ;
543 }
544 }
545 } // only done once
7d855b04 546
547 // temporary matrices needed for the calculation
15687d71 548
9f3b99e2 549 TMatrixD *arrayofErA2[kNPhiSlices], *arrayofEphiA2[kNPhiSlices], *arrayofdEzA2[kNPhiSlices];
7d855b04 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];
15687d71 554
15687d71 555
15687d71 556 for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
7d855b04 557
15687d71 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) ;
7d855b04 566 arrayofEphioverEzA2[k] = new TMatrixD(rows,columns) ;
15687d71 567 arrayofDeltaEzA2[k] = new TMatrixD(rows,columns) ;
568 arrayofEroverEzC2[k] = new TMatrixD(rows,columns) ;
7d855b04 569 arrayofEphioverEzC2[k] = new TMatrixD(rows,columns) ;
15687d71 570 arrayofDeltaEzC2[k] = new TMatrixD(rows,columns) ;
571
7d855b04 572 // zero initialization not necessary, it is done in the constructor of TMatrix
15687d71 573
574 }
7d855b04 575
576
15687d71 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 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
7d855b04 589
15687d71 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
7d855b04 594
15687d71 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
7d855b04 606
9f3b99e2 607 // printf("-----\n%f %f : %e %e\n",r0,phi0,weightA,weightC);
15687d71 608
609 // Summing up the vector components according to their weight
610
611 Int_t ip = 0;
7d855b04 612 for ( Int_t j = 0 ; j < columns ; j++ ) {
15687d71 613 for ( Int_t i = 0 ; i < rows ; i++ ) {
614 for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
7d855b04 615
15687d71 616 // check wether the coordinates were screwed
7d855b04 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) {
15687d71 620 AliError("internal error: coordinate system was screwed during the sum-up");
9f3b99e2 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]);
15687d71 623 AliError("Don't trust the results of the space charge calculation!");
624 }
7d855b04 625
15687d71 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] ;
7d855b04 633
15687d71 634 TMatrixD &erC = *arrayofErC2[k] ;
635 TMatrixD &ephiC = *arrayofEphiC2[k];
636 TMatrixD &dEzC = *arrayofdEzC2[k] ;
7d855b04 637
15687d71 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
7d855b04 651
652
15687d71 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
7d855b04 656
15687d71 657 Int_t phiPoints = (Int_t) grid2(1);
658 Int_t phiPOC = (Int_t) grid2(4);
7d855b04 659
15687d71 660 // printf("%d %d\n",phiPOC,flagRadSym);
7d855b04 661
662 for (Int_t phiiC = 1; phiiC<phiPOC; phiiC++) { // just used for non-radial symetric table
663
15687d71 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
7d855b04 670
9f3b99e2 671 // printf("%f %f : %e %e\n",r0,phi0R,weightA,weightC);
7d855b04 672
15687d71 673 // Summing up the vector components according to their weight
674 ip = 0;
7d855b04 675 for ( Int_t j = 0 ; j < columns ; j++ ) {
15687d71 676 for ( Int_t i = 0 ; i < rows ; i++ ) {
677 for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
7d855b04 678
15687d71 679 // Note: rotating the coordinated during the sum up
7d855b04 680
15687d71 681 Int_t rotVal = (phiPoints/phiPOC)*phiiC;
682 Int_t ipR = -1;
7d855b04 683
15687d71 684 if ((ip%phiPoints)>=rotVal) {
685 ipR = ip-rotVal;
686 } else {
687 ipR = ip+(phiPoints-rotVal);
688 }
7d855b04 689
15687d71 690 // unfortunately, the lookup tables were produced to be faster for phi symmetric charges
7d855b04 691 // This will be the most frequent usage
15687d71 692 // That's why we have to do this here and not outside the loop ...
7d855b04 693
15687d71 694 TMatrixD &erA = *arrayofErA2[k] ;
695 TMatrixD &ephiA = *arrayofEphiA2[k];
696 TMatrixD &dEzA = *arrayofdEzA2[k] ;
7d855b04 697
15687d71 698 TMatrixD &erC = *arrayofErC2[k] ;
699 TMatrixD &ephiC = *arrayofEphiC2[k];
700 TMatrixD &dEzC = *arrayofdEzC2[k] ;
7d855b04 701
15687d71 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
9f3b99e2 720 // printf("POC: (r,phi,z) = (%f %f %f) | weight(A,C): %03.1lf %03.1lf\n",r0,phi0,z0, weightA, weightC);
7d855b04 721
15687d71 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];
7d855b04 749 TMatrixD &deltaEzC = *arrayofDeltaEzC2[k];
750
15687d71 751 for ( Int_t i = 0 ; i < rows ; i++ ) { // r loop
7d855b04 752 for ( Int_t j = columns-1 ; j >= 0 ; j-- ) {// z loop
15687d71 753 // Count backwards to facilitate integration over Z
754
7d855b04 755 Int_t index = 1 ; // Simpsons rule if N=odd.If N!=odd then add extra point by trapezoidal rule.
15687d71 756
7d855b04 757 erOverEzA(i,j) = 0;
15687d71 758 ephiOverEzA(i,j) = 0;
759 deltaEzA(i,j) = 0;
7d855b04 760 erOverEzC(i,j) = 0;
761 ephiOverEzC(i,j) = 0;
15687d71 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 ) {
7d855b04 802 erOverEzA(i,j) = 0;
15687d71 803 erOverEzC(i,j) = 0;
7d855b04 804 ephiOverEzA(i,j) = 0;
15687d71 805 ephiOverEzC(i,j) = 0;
7d855b04 806 deltaEzA(i,j) = 0;
15687d71 807 deltaEzC(i,j) = 0;
808 }
809 }
810 }
811
812 }
7d855b04 813
15687d71 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] ;
7d855b04 822
15687d71 823 // final lookup table
2bf29b72 824 TMatrixF &erOverEzFinal = *fLookUpErOverEz[k] ;
825 TMatrixF &ephiOverEzFinal = *fLookUpEphiOverEz[k];
826 TMatrixF &deltaEzFinal = *fLookUpDeltaEz[k] ;
7d855b04 827
15687d71 828 for ( Int_t j = 0 ; j < kNZ ; j++ ) {
829
830 z = TMath::Abs(fgkZList[j]) ; // z position is symmetric
7d855b04 831
832 for ( Int_t i = 0 ; i < kNR ; i++ ) {
15687d71 833 r = fgkRList[i] ;
834
835 // Interpolate Lookup tables onto standard grid
836 if (fgkZList[j]>0) {
7d855b04 837 erOverEzFinal(i,j) += Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
15687d71 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 {
7d855b04 844 erOverEzFinal(i,j) += Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
15687d71 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
7d855b04 855
856
15687d71 857 // clear the temporary arrays lists
858 for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
859
7d855b04 860 delete arrayofErA2[k];
15687d71 861 delete arrayofEphiA2[k];
862 delete arrayofdEzA2[k];
7d855b04 863 delete arrayofErC2[k];
15687d71 864 delete arrayofEphiC2[k];
865 delete arrayofdEzC2[k];
866
7d855b04 867 delete arrayofEroverEzA2[k];
15687d71 868 delete arrayofEphioverEzA2[k];
869 delete arrayofDeltaEzA2[k];
7d855b04 870 delete arrayofEroverEzC2[k];
15687d71 871 delete arrayofEphioverEzC2[k];
872 delete arrayofDeltaEzC2[k];
873
874 }
7d855b04 875
15687d71 876 fRPhi->Close();
7d855b04 877
15687d71 878 // FINISHED
879
880 fInitLookUp = kTRUE;
881
882}
883
884void AliTPCSpaceCharge3D::InitSpaceCharge3DDistortionCourse() {
7d855b04 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 ...
cc3e558a 894
895 if (fInitLookUp) {
896 AliInfo("Lookup table was already initialized!");
897 // return;
898 }
899
900 AliInfo("Preparation of the weighted look-up table");
7d855b04 901
15687d71 902 TFile *f = new TFile(fSCLookUpPOCsFileName3D.Data(),"READ");
903 if ( !f ) {
cc3e558a 904 AliError("Precalculated POC-looup-table could not be found");
905 return;
906 }
907
908 // units are in [m]
7d855b04 909 TVector *gridf = (TVector*) f->Get("constants");
cc3e558a 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;
7d855b04 915
cc3e558a 916 Bool_t flagRadSym = 0;
917 if (grid(1)==1 && grid(4)==1) {
15687d71 918 // AliInfo("LOOK UP TABLE IS RADIAL SYMETTRIC - Field in Phi is ZERO");
cc3e558a 919 flagRadSym=1;
920 }
921
7d855b04 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
cc3e558a 925
15687d71 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]
7d855b04 929
cc3e558a 930 // temporary matrices needed for the calculation
9f3b99e2 931 TMatrixD *arrayofErA[kNPhiSlices], *arrayofEphiA[kNPhiSlices], *arrayofdEzA[kNPhiSlices];
932 TMatrixD *arrayofErC[kNPhiSlices], *arrayofEphiC[kNPhiSlices], *arrayofdEzC[kNPhiSlices];
cc3e558a 933
9f3b99e2 934 TMatrixD *arrayofEroverEzA[kNPhiSlices], *arrayofEphioverEzA[kNPhiSlices], *arrayofDeltaEzA[kNPhiSlices];
935 TMatrixD *arrayofEroverEzC[kNPhiSlices], *arrayofEphioverEzC[kNPhiSlices], *arrayofDeltaEzC[kNPhiSlices];
cc3e558a 936
7d855b04 937
cc3e558a 938 for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
7d855b04 939
cc3e558a 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 }
7d855b04 958
cc3e558a 959 // list of points as used in the interpolation (during sum up)
9f3b99e2 960 Double_t rlist[kNRows], zedlist[kNColumns] , philist[kNPhiSlices];
cc3e558a 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 ;
7d855b04 965 for ( Int_t j = 0 ; j < columns ; j++ ) {
cc3e558a 966 zedlist[j] = j * gridSizeZ ;
967 }
968 }
969 } // only done once
7d855b04 970
971
cc3e558a 972 TTree *treePOC = (TTree*)f->Get("POCall");
973
974 TVector *bEr = 0; TVector *bEphi= 0; TVector *bEz = 0;
7d855b04 975
cc3e558a 976 treePOC->SetBranchAddress("Er",&bEr);
977 if (!flagRadSym) treePOC->SetBranchAddress("Ephi",&bEphi);
978 treePOC->SetBranchAddress("Ez",&bEz);
979
cc3e558a 980
981 // Read the complete tree and do a weighted sum-up over the POC configurations
982 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
7d855b04 983
cc3e558a 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
15687d71 987 for (Int_t itreepC=0; itreepC<treeNumPOC; itreepC++) { // ------------- loop over POC configurations in tree
7d855b04 988
cc3e558a 989 treePOC->GetEntry(itreepC);
990
cc3e558a 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
15687d71 996 ipC++; // POC configuration counter
cc3e558a 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]
92a85338 1000 Double_t weightA = GetSpaceChargeDensity(r0*100,phi0, z0*100,0);
1001 Double_t weightC = GetSpaceChargeDensity(r0*100,phi0,-z0*100,0);
7d855b04 1002
cc3e558a 1003 // Summing up the vector components according to their weight
1004
1005 Int_t ip = 0;
7d855b04 1006 for ( Int_t j = 0 ; j < columns ; j++ ) {
cc3e558a 1007 for ( Int_t i = 0 ; i < rows ; i++ ) {
1008 for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
7d855b04 1009
cc3e558a 1010 // check wether the coordinates were screwed
7d855b04 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) {
cc3e558a 1014 AliError("internal error: coordinate system was screwed during the sum-up");
9f3b99e2 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]);
15687d71 1017 AliError("Don't trust the results of the space charge calculation!");
cc3e558a 1018 }
7d855b04 1019
cc3e558a 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];
15687d71 1026 TMatrixD &dEzA = *arrayofdEzA[k] ;
7d855b04 1027
cc3e558a 1028 TMatrixD &erC = *arrayofErC[k] ;
1029 TMatrixD &ephiC = *arrayofEphiC[k];
15687d71 1030 TMatrixD &dEzC = *arrayofdEzC[k] ;
7d855b04 1031
cc3e558a 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) {
15687d71 1036 ephiA(i,j) += ((*bEphi)(ip)) * weightA/100; // [V/rad]
1037 ephiC(i,j) += ((*bEphi)(ip)) * weightC/100; // [V/rad]
cc3e558a 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 }
15687d71 1046 } // end coordinate loop
7d855b04 1047
1048
cc3e558a 1049 // Rotation and summation in the rest of the dPhiSteps
15687d71 1050 // which were not stored in the this tree due to storage & symmetry reasons
cc3e558a 1051
1052 Int_t phiPoints = (Int_t) grid(1);
1053 Int_t phiPOC = (Int_t) grid(4);
7d855b04 1054
15687d71 1055 // printf("%d %d\n",phiPOC,flagRadSym);
7d855b04 1056
1057 for (Int_t phiiC = 1; phiiC<phiPOC; phiiC++) { // just used for non-radial symetric table
1058
cc3e558a 1059 r0 = coordPOC(ipC,0);
1060 phi0 = coordPOC(ipC,1);
1061 z0 = coordPOC(ipC,2);
7d855b04 1062
cc3e558a 1063 ipC++; // POC conf. counter
7d855b04 1064
cc3e558a 1065 // weights (charge density) at POC position on the A and C side (in C/m^3/e0)
1066 // note: coordinates are in [cm]
7d855b04 1067 weightA = GetSpaceChargeDensity(r0*100,phi0, z0*100,0);
92a85338 1068 weightC = GetSpaceChargeDensity(r0*100,phi0,-z0*100,0);
7d855b04 1069
9f3b99e2 1070 // printf("%f %f %f: %e %e\n",r0,phi0,z0,weightA,weightC);
7d855b04 1071
cc3e558a 1072 // Summing up the vector components according to their weight
1073 ip = 0;
7d855b04 1074 for ( Int_t j = 0 ; j < columns ; j++ ) {
cc3e558a 1075 for ( Int_t i = 0 ; i < rows ; i++ ) {
1076 for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
7d855b04 1077
cc3e558a 1078 // Note: rotating the coordinated during the sum up
7d855b04 1079
cc3e558a 1080 Int_t rotVal = (phiPoints/phiPOC)*phiiC;
1081 Int_t ipR = -1;
7d855b04 1082
cc3e558a 1083 if ((ip%phiPoints)>=rotVal) {
1084 ipR = ip-rotVal;
1085 } else {
1086 ipR = ip+(phiPoints-rotVal);
1087 }
7d855b04 1088
cc3e558a 1089 // unfortunately, the lookup tables were produced to be faster for phi symmetric charges
7d855b04 1090 // This will be the most frequent usage
cc3e558a 1091 // That's why we have to do this here and not outside the loop ...
7d855b04 1092
cc3e558a 1093 TMatrixD &erA = *arrayofErA[k] ;
1094 TMatrixD &ephiA = *arrayofEphiA[k];
1095 TMatrixD &dEzA = *arrayofdEzA[k] ;
7d855b04 1096
cc3e558a 1097 TMatrixD &erC = *arrayofErC[k] ;
1098 TMatrixD &ephiC = *arrayofEphiC[k];
1099 TMatrixD &dEzC = *arrayofdEzC[k] ;
7d855b04 1100
cc3e558a 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) {
15687d71 1105 ephiA(i,j) += ((*bEphi)(ipR)) * weightA/100; // [V/rad]
1106 ephiC(i,j) += ((*bEphi)(ipR)) * weightC/100; // [V/rad]
cc3e558a 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)
7d855b04 1118
cc3e558a 1119
9f3b99e2 1120 // printf("POC: (r,phi,z) = (%f %f %f) | weight(A,C): %03.1lf %03.1lf\n",r0,phi0,z0, weightA, weightC);
7d855b04 1121
15687d71 1122 }
1123
cc3e558a 1124
cc3e558a 1125
1126 // -------------------------------------------------------------------------------
1127 // Division by the Ez (drift) field and integration along z
1128
15687d71 1129 AliInfo("Division and integration");
1130
cc3e558a 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
15687d71 1135 // matrices holding the solution - summation of POC charges // see above
cc3e558a 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
15687d71 1143 // matrices which will contain the integrated fields (divided by the drift field)
cc3e558a 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];
7d855b04 1149 TMatrixD &deltaEzC = *arrayofDeltaEzC[k];
1150
cc3e558a 1151 for ( Int_t i = 0 ; i < rows ; i++ ) { // r loop
7d855b04 1152 for ( Int_t j = columns-1 ; j >= 0 ; j-- ) {// z loop
cc3e558a 1153 // Count backwards to facilitate integration over Z
1154
7d855b04 1155 Int_t index = 1 ; // Simpsons rule if N=odd.If N!=odd then add extra point by trapezoidal rule.
cc3e558a 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 }
15687d71 1168 deltaEzA(i,j) += index*(gridSizeZ/3.0)*dezA(i,m)/(-1) ;
1169 deltaEzC(i,j) += index*(gridSizeZ/3.0)*dezC(i,m)/(-1) ;
cc3e558a 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 }
15687d71 1182 deltaEzA(i,j) -= (gridSizeZ/3.0)*dezA(i,columns-1)/(-1) ;
1183 deltaEzC(i,j) -= (gridSizeZ/3.0)*dezC(i,columns-1)/(-1) ;
cc3e558a 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 }
15687d71 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) ;
cc3e558a 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 }
15687d71 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) ;
cc3e558a 1204 }
1205 if ( j == columns-1 ) {
7d855b04 1206 erOverEzA(i,j) = 0;
15687d71 1207 erOverEzC(i,j) = 0;
cc3e558a 1208 if (!flagRadSym) {
7d855b04 1209 ephiOverEzA(i,j) = 0;
15687d71 1210 ephiOverEzC(i,j) = 0;
cc3e558a 1211 }
7d855b04 1212 deltaEzA(i,j) = 0;
15687d71 1213 deltaEzC(i,j) = 0;
cc3e558a 1214 }
1215 }
1216 }
1217
1218 }
15687d71 1219
7d855b04 1220
1221
cc3e558a 1222 AliInfo("Interpolation to Standard grid");
1223
1224 // -------------------------------------------------------------------------------
15687d71 1225 // Interpolate results onto the standard grid which is used for all AliTPCCorrections classes
cc3e558a 1226
7d855b04 1227 const Int_t order = 1 ; // Linear interpolation = 1, Quadratic = 2
cc3e558a 1228
1229 Double_t r, phi, z ;
1230 for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
1231 phi = fgkPhiList[k] ;
7d855b04 1232
2bf29b72 1233 TMatrixF &erOverEz = *fLookUpErOverEz[k] ;
1234 TMatrixF &ephiOverEz = *fLookUpEphiOverEz[k];
1235 TMatrixF &deltaEz = *fLookUpDeltaEz[k] ;
7d855b04 1236
cc3e558a 1237 for ( Int_t j = 0 ; j < kNZ ; j++ ) {
1238
1239 z = TMath::Abs(fgkZList[j]) ; // z position is symmetric
7d855b04 1240
1241 for ( Int_t i = 0 ; i < kNR ; i++ ) {
cc3e558a 1242 r = fgkRList[i] ;
1243
1244 // Interpolate Lookup tables onto standard grid
1245 if (fgkZList[j]>0) {
7d855b04 1246 erOverEz(i,j) = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
cc3e558a 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 {
7d855b04 1253 erOverEz(i,j) = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
cc3e558a 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
7d855b04 1266
cc3e558a 1267 // clear the temporary arrays lists
1268 for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
1269
7d855b04 1270 delete arrayofErA[k];
cc3e558a 1271 delete arrayofEphiA[k];
1272 delete arrayofdEzA[k];
7d855b04 1273 delete arrayofErC[k];
cc3e558a 1274 delete arrayofEphiC[k];
1275 delete arrayofdEzC[k];
1276
7d855b04 1277 delete arrayofEroverEzA[k];
cc3e558a 1278 delete arrayofEphioverEzA[k];
1279 delete arrayofDeltaEzA[k];
7d855b04 1280 delete arrayofEroverEzC[k];
cc3e558a 1281 delete arrayofEphioverEzC[k];
1282 delete arrayofDeltaEzC[k];
1283
1284 }
1285
cc3e558a 1286 fInitLookUp = kTRUE;
1287
1288}
1289
1290
15687d71 1291void AliTPCSpaceCharge3D::SetSCDataFileName(TString fname) {
7d855b04 1292 /// Set & load the Space charge density distribution from a file
1293 /// (linear interpolation onto a standard grid)
1294
cc3e558a 1295
1296 fSCDataFileName = fname;
1297
15687d71 1298 TFile *f = new TFile(fSCDataFileName.Data(),"READ");
7d855b04 1299 if (!f) {
cc3e558a 1300 AliError(Form("File %s, which should contain the space charge distribution, could not be found",
15687d71 1301 fSCDataFileName.Data()));
cc3e558a 1302 return;
1303 }
7d855b04 1304
15687d71 1305 TH2F *densityRZ = (TH2F*) f->Get("SpaceChargeInRZ");
7d855b04 1306 if (!densityRZ) {
cc3e558a 1307 AliError(Form("The indicated file (%s) does not contain a histogram called %s",
15687d71 1308 fSCDataFileName.Data(),"SpaceChargeInRZ"));
1309 return;
1310 }
1311
1312 TH3F *densityRPhi = (TH3F*) f->Get("SpaceChargeInRPhi");
7d855b04 1313 if (!densityRPhi) {
15687d71 1314 AliError(Form("The indicated file (%s) does not contain a histogram called %s",
1315 fSCDataFileName.Data(),"SpaceChargeInRPhi"));
cc3e558a 1316 return;
1317 }
7d855b04 1318
cc3e558a 1319
1320 Double_t r, phi, z ;
15687d71 1321
1322 TMatrixD &scDensityInRZ = *fSCdensityInRZ;
1323 TMatrixD &scDensityInRPhiA = *fSCdensityInRPhiA;
1324 TMatrixD &scDensityInRPhiC = *fSCdensityInRPhiC;
cc3e558a 1325 for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
1326 phi = fgkPhiList[k] ;
2bf29b72 1327 TMatrixF &scDensity = *fSCdensityDistribution[k] ;
cc3e558a 1328 for ( Int_t j = 0 ; j < kNZ ; j++ ) {
7d855b04 1329 z = fgkZList[j] ;
1330 for ( Int_t i = 0 ; i < kNR ; i++ ) {
cc3e558a 1331 r = fgkRList[i] ;
15687d71 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 ) {
7d855b04 1339 if (z>0)
15687d71 1340 scDensityInRPhiA(i,k) = densityRPhi->Interpolate(r,phi,0.499); // A side
7d855b04 1341 else
15687d71 1342 scDensityInRPhiC(i,k) = densityRPhi->Interpolate(r,phi,-0.499); // C side
1343 }
1344
1345 // Full 3D configuration ...
7d855b04 1346 if (z>0)
1347 scDensity(i,j) = scDensityInRZ(i,j) + scDensityInRPhiA(i,k);
15687d71 1348 else
7d855b04 1349 scDensity(i,j) = scDensityInRZ(i,j) + scDensityInRPhiC(i,k);
cc3e558a 1350 }
1351 }
1352 }
1353
cc3e558a 1354 f->Close();
1355
1356 fInitLookUp = kFALSE;
1357
7d855b04 1358
cc3e558a 1359}
1360
92a85338 1361void AliTPCSpaceCharge3D::SetInputSpaceCharge(TH3 * hisSpaceCharge3D, TH2 * hisRPhi, TH2* hisRZ, Double_t norm){
7d855b04 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
92a85338 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++ ) {
7d855b04 1386 z = fgkZList[j] ;
1387 for ( Int_t i = 0 ; i < kNR ; i++ ) {
92a85338 1388 // Full 3D configuration ...
1389 r = fgkRList[i] ;
7d855b04 1390 if (r>rmin && r<rmax && z>zmin && z< zmax){
92a85338 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) {
7d855b04 1397 if (z>0)
92a85338 1398 scDensityInRPhiA(i,k) = norm*fSpaceChargeHistogramRPhi->Interpolate(phi,r); // A side
7d855b04 1399 else
92a85338 1400 scDensityInRPhiC(i,k) = norm*fSpaceChargeHistogramRPhi->Interpolate(phi+TMath::TwoPi(),r); // C side
1401 }
7d855b04 1402
1403 if (z>0)
1404 scDensity(i,j) = norm*fSpaceChargeHistogram3D->Interpolate(phi,r,z);
92a85338 1405 else
1406 scDensity(i,j) = norm*fSpaceChargeHistogram3D->Interpolate(phi,r,z);
1407 }
1408 }
1409 }
1410 }
7d855b04 1411
92a85338 1412 fInitLookUp = kFALSE;
1413
1414}
1415
cc3e558a 1416
15687d71 1417Float_t AliTPCSpaceCharge3D::GetSpaceChargeDensity(Float_t r, Float_t phi, Float_t z, Int_t mode) {
7d855b04 1418 /// returns the (input) space charge density at a given point according
1419 /// Note: input in [cm], output in [C/m^3/e0] !!
cc3e558a 1420
cc3e558a 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; //
15687d71 1427 Float_t sc = 0;
cc3e558a 1428
15687d71 1429 if (mode == 0) { // return full load
7d855b04 1430 sc = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
15687d71 1431 fgkRList, fgkZList, fgkPhiList, fSCdensityDistribution );
7d855b04 1432
15687d71 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 );
7d855b04 1436
15687d71 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 }
7d855b04 1451
9f3b99e2 1452 // printf("%f %f %f: %f\n",r,phi,z,sc);
7d855b04 1453
cc3e558a 1454 return sc;
1455}
1456
1457
15687d71 1458TH2F * AliTPCSpaceCharge3D::CreateHistoSCinXY(Float_t z, Int_t nx, Int_t ny, Int_t mode) {
7d855b04 1459 /// return a simple histogramm containing the space charge distribution (input for the calculation)
cc3e558a 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);
7d855b04 1468
cc3e558a 1469 Float_t r = TMath::Sqrt(xp*xp+yp*yp);
7d855b04 1470 Float_t phi = TMath::ATan2(yp,xp);
1471
cc3e558a 1472 if (85.<=r && r<=250.) {
15687d71 1473 Float_t sc = GetSpaceChargeDensity(r,phi,z,mode)/fgke0; // in [C/m^3/e0]
7d855b04 1474 h->SetBinContent(ix,iy,sc);
cc3e558a 1475 } else {
1476 h->SetBinContent(ix,iy,0.);
1477 }
1478 }
1479 }
7d855b04 1480
cc3e558a 1481 return h;
7d855b04 1482}
cc3e558a 1483
15687d71 1484TH2F * AliTPCSpaceCharge3D::CreateHistoSCinZR(Float_t phi, Int_t nz, Int_t nr,Int_t mode ) {
7d855b04 1485 /// return a simple histogramm containing the space charge distribution (input for the calculation)
cc3e558a 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);
15687d71 1494 Float_t sc = GetSpaceChargeDensity(r,phi,z,mode)/fgke0; // in [C/m^3/e0]
cc3e558a 1495 h->SetBinContent(iz,ir,sc);
1496 }
1497 }
1498
1499 return h;
7d855b04 1500}
cc3e558a 1501
1502void AliTPCSpaceCharge3D::WriteChargeDistributionToFile(const char* fname) {
7d855b04 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
cc3e558a 1507
1508 TFile *f = new TFile(fname,"RECREATE");
7d855b04 1509
cc3e558a 1510 // some grid, not too course
15687d71 1511 Int_t nr = 350;
cc3e558a 1512 Int_t nphi = 180;
15687d71 1513 Int_t nz = 500;
cc3e558a 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
15687d71 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);
7d855b04 1526
cc3e558a 1527 for (Int_t ir=1;ir<=nr;++ir) {
15687d71 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);
7d855b04 1531
15687d71 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]
7d855b04 1536
15687d71 1537 // setting of mb multiplicity and Interaction rate
1538 Double_t multiplicity = 950;
1539 Double_t intRate = 7800;
cc3e558a 1540
15687d71 1541 // calculation of "scaled" parameters
1542 Double_t a = multiplicity*intRate/79175;
1543 Double_t b = a/lZ;
7d855b04 1544
15687d71 1545 Double_t charge = (a - b*zpM)/(rpM*rpM); // charge in [C/m^3/e0]
7d855b04 1546
15687d71 1547 charge = charge*fgke0; // [C/m^3]
7d855b04 1548
15687d71 1549 if (zp<0) charge *= 0.9; // e.g. slightly less on C side due to front absorber
cc3e558a 1550
15687d71 1551 // charge = 0; // for tests
7d855b04 1552 histoZR->SetBinContent(ir,iz,charge);
15687d71 1553 }
1554 }
7d855b04 1555
15687d71 1556 histoZR->Write("SpaceChargeInRZ");
7d855b04 1557
15687d71 1558
1559 // Charge distribution in RPhi (e.g. Floating GG wire) ------------
7d855b04 1560
15687d71 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
7d855b04 1565
15687d71 1566 // some 'arbitrary' GG leaks
1567 Int_t nGGleaks = 5;
1568 Double_t secPosA[5] = {3,6,6,11,13}; // sector
752b0cc7 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
15687d71 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);
7d855b04 1579
15687d71 1580 Double_t charge = 0;
7d855b04 1581
15687d71 1582 for (Int_t igg = 0; igg<nGGleaks; igg++) { // loop over GG leaks
7d855b04 1583
15687d71 1584 // A side
7d855b04 1585 Double_t secPos = secPosA[igg];
15687d71 1586 Double_t radialPos = radialPosA[igg];
1587
1588 if (zp<0) { // C side
7d855b04 1589 secPos = secPosC[igg];
15687d71 1590 radialPos = radialPosC[igg];
7d855b04 1591 }
15687d71 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
7d855b04 1596 charge = 300;
15687d71 1597 }
7d855b04 1598
1599 }
1600
cc3e558a 1601 charge = charge*fgke0; // [C/m^3]
1602
7d855b04 1603 histoRPhi->SetBinContent(ir,iphi,iz,charge);
cc3e558a 1604 }
1605 }
1606 }
1607
15687d71 1608 histoRPhi->Write("SpaceChargeInRPhi");
cc3e558a 1609
cc3e558a 1610 f->Close();
7d855b04 1611
cc3e558a 1612}
1613
1614
1615void AliTPCSpaceCharge3D::Print(const Option_t* option) const {
7d855b04 1616 /// Print function to check the settings of the boundary vectors
1617 /// option=="a" prints the C0 and C1 coefficents for calibration purposes
cc3e558a 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);
7d855b04 1627 }
1628
cc3e558a 1629 if (!fInitLookUp) AliError("Lookup table was not initialized! You should do InitSpaceCharge3DDistortion() ...");
1630
7d855b04 1631}
92a85338 1632
1633
1634
1635void AliTPCSpaceCharge3D::InitSpaceCharge3DPoisson(Int_t kRows, Int_t kColumns, Int_t kPhiSlices, Int_t kIterations){
7d855b04 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;
92a85338 1644 //
7d855b04 1645 const Int_t order = 1 ; // Linear interpolation = 1, Quadratic = 2
92a85338 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
7d855b04 1651 TMatrixD *arrayofArrayV[kPhiSlices], *arrayofCharge[kPhiSlices] ;
1652 TMatrixD *arrayofEroverEz[kPhiSlices], *arrayofEphioverEz[kPhiSlices], *arrayofDeltaEz[kPhiSlices] ;
92a85338 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 }
7d855b04 1661
92a85338 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
7d855b04 1677
92a85338 1678 const Int_t symmetry = 0;
7d855b04 1679
92a85338 1680 // Set bondaries and solve Poisson's equation --------------------------
7d855b04 1681
92a85338 1682 if ( !fInitLookUp ) {
7d855b04 1683
92a85338 1684 AliInfo(Form("Solving the poisson equation (~ %d sec)",2*10*(int)(kPhiSlices/10)));
7d855b04 1685
92a85338 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] ;
7d855b04 1691
92a85338 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
7d855b04 1695// arrayV(i,j) = 0.0 ;
92a85338 1696// charge(i,j) = 0.0 ;
1697
1698// // Float_t radius0 = rlist[i] ;
1699// // Float_t phi0 = gridSizePhi * k ;
7d855b04 1700
92a85338 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// }
7d855b04 1709// }
1710
1711 for ( Int_t i = 1 ; i < kRows-1 ; i++ ) {
1712 for ( Int_t j = 1 ; j < kColumns-1 ; j++ ) {
92a85338 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]);
7d855b04 1717 arrayV(i,j) = 0.0 ;
92a85338 1718 charge(i,j) = fSpaceChargeHistogram3D->Interpolate(phi0,radius0,z0);
1719 }
1720 }
7d855b04 1721 }
92a85338 1722 AliSysInfo::AddStamp("RunPoisson", 2,side,0);
7d855b04 1723
92a85338 1724 // Solve Poisson's equation in 3D cylindrical coordinates by relaxation technique
1725 // Allow for different size grid spacing in R and Z directions
7d855b04 1726
1727 // PoissonRelaxation3D( arrayofArrayV, arrayofCharge,
92a85338 1728 // arrayofEroverEz, arrayofEphioverEz, arrayofDeltaEz,
7d855b04 1729 // kRows, kColumns, kPhiSlices, gridSizePhi, kIterations,
92a85338 1730 // symmetry , fROCdisplacement) ;
7d855b04 1731 PoissonRelaxation3D( arrayofArrayV, arrayofCharge,
92a85338 1732 arrayofEroverEz, arrayofEphioverEz, arrayofDeltaEz,
7d855b04 1733 kRows, kColumns, kPhiSlices, gridSizePhi, kIterations,
92a85338 1734 symmetry ) ;
7d855b04 1735
92a85338 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] ;
7d855b04 1740
92a85338 1741 TMatrixF &erOverEz = *fLookUpErOverEz[k] ;
1742 TMatrixF &ephiOverEz = *fLookUpEphiOverEz[k];
1743 TMatrixF &deltaEz = *fLookUpDeltaEz[k] ;
7d855b04 1744
92a85338 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)
7d855b04 1748
92a85338 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
7d855b04 1751
1752 for ( Int_t i = 0 ; i < kNR ; i++ ) {
92a85338 1753 r = fgkRList[i] ;
1754
1755 // Interpolate basicLookup tables; once for each rod, then sum the results
7d855b04 1756 erOverEz(i,j) = Interpolate3DTable(order, r, z, phi, kRows, kColumns, kPhiSlices,
92a85338 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
7d855b04 1768 AliSysInfo::AddStamp("Interpolate Poisson", 3,side,0);
92a85338 1769 if ( side == 0 ) AliInfo(" A side done");
1770 if ( side == 1 ) AliInfo(" C side done");
1771 } // end side loop
1772 }
7d855b04 1773
92a85338 1774 // clear the temporary arrays lists
1775 for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
1776 delete arrayofArrayV[k];
1777 delete arrayofCharge[k];
7d855b04 1778 delete arrayofEroverEz[k];
92a85338 1779 delete arrayofEphioverEz[k];
1780 delete arrayofDeltaEz[k];
1781 }
7d855b04 1782
92a85338 1783
1784 fInitLookUp = kTRUE;
1785
cc3e558a 1786}
92a85338 1787
8a94851d 1788
1789
1790AliTPCSpaceCharge3D * AliTPCSpaceCharge3D::MakeCorrection22D(const char* fileName, Double_t multiplicity, Double_t intRate, Double_t epsIROC, Double_t epsOROC,
7d855b04 1791 Double_t gasfactor,
8a94851d 1792 Double_t radialScaling){
7d855b04 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");
8a94851d 1809 // some grid, not too coarse
1810 const Int_t nr = 350;
1811 const Int_t nphi = 180;
1812 const Int_t nz = 500;
7d855b04 1813 const Double_t kROROC=134; // hardwired OROC radius
8a94851d 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]
7d855b04 1831 Double_t integralRadialExponent2 = TMath::Power(radiusOuter,radialExponent+1) * 1./(radialExponent+1)
8a94851d 1832 - TMath::Power(radiusInner,radialExponent+1) * 1./(radialExponent+1);
7d855b04 1833
1834 radialExponent = -radialScaling; // user set
8a94851d 1835 Double_t integralRadialExponentUser = 0.;
1836 if(radialScaling > 1 + 0.000001 || radialScaling < 1 - 0.000001 ) // to avoid n = -1
7d855b04 1837 integralRadialExponentUser = TMath::Power(radiusOuter,radialExponent+1) * 1./(radialExponent+1)
8a94851d 1838 - TMath::Power(radiusInner,radialExponent+1) * 1./(radialExponent+1);
1839 else
1840 integralRadialExponentUser = TMath::Log(radiusOuter) - TMath::Log(radiusInner);
7d855b04 1841
8a94851d 1842 Double_t normRadialExponent = integralRadialExponent2 / integralRadialExponentUser;
7d855b04 1843
8a94851d 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);
7d855b04 1848 Double_t eps = (rp <kROROC) ? epsIROC:epsOROC;
8a94851d 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]
7d855b04 1853
8a94851d 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
7d855b04 1858
8a94851d 1859 charge = charge*fgke0; // [C/m^3]
1860 if (zp<0) charge *= 0.9; // Slightly less on C side due to front absorber
1861
7d855b04 1862 histoZR->SetBinContent(ir, iz, charge);
8a94851d 1863 }
1864 }
7d855b04 1865
8a94851d 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,
7d855b04 1873 2, -1, 1); // z part - to allow A and C side differences
8a94851d 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