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