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