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