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