]> 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
- bug fix
[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),
60 fSCDataFileName((char*)"$(ALICE_ROOT)/TPC/Calib/maps/sc_3D_distribution_Sim.root"),
61 fSCLookUpPOCsFileName((char*)"$(ALICE_ROOT)/TPC/Calib/maps/sc_3D_raw_18-18-26_17p-18p-25p-MN30.root")
62{
63 //
64 // default constructor
65 //
66
67 // Array which will contain the solution according to the setted charge density distribution
68 // see InitSpaceCharge3DDistortion() function
69 for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
70 fLookUpErOverEz[k] = new TMatrixD(kNR,kNZ);
71 fLookUpEphiOverEz[k] = new TMatrixD(kNR,kNZ);
72 fLookUpDeltaEz[k] = new TMatrixD(kNR,kNZ);
73 fSCdensityDistribution[k] = new TMatrixD(kNR,kNZ);
74 }
75
76 printf("%s\n",fSCDataFileName);
77 printf("%s\n",fSCLookUpPOCsFileName);
78 SetSCDataFileName(fSCDataFileName);
79
80}
81
82AliTPCSpaceCharge3D::~AliTPCSpaceCharge3D() {
83 //
84 // default destructor
85 //
86
87 for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
88 delete fLookUpErOverEz[k];
89 delete fLookUpEphiOverEz[k];
90 delete fLookUpDeltaEz[k];
91 delete fSCdensityDistribution[k];
92 }
93
94 // delete fSCdensityDistribution;
95}
96
97
98
99void AliTPCSpaceCharge3D::Init() {
100 //
101 // Initialization funtion
102 //
103
104 AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
105 if (!magF) AliError("Magneticd field - not initialized");
106 Double_t bzField = magF->SolenoidField()/10.; //field in T
107 AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
108 if (!param) AliError("Parameters - not initialized");
109 Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us] // From dataBase: to be updated: per second (ideally)
110 Double_t ezField = 400; // [V/cm] // to be updated: never (hopefully)
111 Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
112 // Correction Terms for effective omegaTau; obtained by a laser calibration run
113 SetOmegaTauT1T2(wt,fT1,fT2);
114
115 InitSpaceCharge3DDistortion(); // fill the look up table
116}
117
118void AliTPCSpaceCharge3D::Update(const TTimeStamp &/*timeStamp*/) {
119 //
120 // Update function
121 //
122 AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
123 if (!magF) AliError("Magneticd field - not initialized");
124 Double_t bzField = magF->SolenoidField()/10.; //field in T
125 AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
126 if (!param) AliError("Parameters - not initialized");
127 Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us] // From dataBase: to be updated: per second (ideally)
128 Double_t ezField = 400; // [V/cm] // to be updated: never (hopefully)
129 Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
130 // Correction Terms for effective omegaTau; obtained by a laser calibration run
131 SetOmegaTauT1T2(wt,fT1,fT2);
132
133 // SetCorrectionFactor(1.); // should come from some database
134
135}
136
137
138
139void AliTPCSpaceCharge3D::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {
140 //
141 // Calculates the correction due the Space Charge effect within the TPC drift volume
142 //
143
144 if (!fInitLookUp) {
145 AliInfo("Lookup table was not initialized! Performing the inizialisation now ...");
146 InitSpaceCharge3DDistortion();
147 }
148
149 Int_t order = 1 ; // FIXME: hardcoded? Linear interpolation = 1, Quadratic = 2
150
151 Double_t intEr, intEphi, intdEz ;
152 Double_t r, phi, z ;
153 Int_t sign;
154
155 r = TMath::Sqrt( x[0]*x[0] + x[1]*x[1] ) ;
156 phi = TMath::ATan2(x[1],x[0]) ;
157 if ( phi < 0 ) phi += TMath::TwoPi() ; // Table uses phi from 0 to 2*Pi
158 z = x[2] ; // Create temporary copy of x[2]
159
160 if ( (roc%36) < 18 ) {
161 sign = 1; // (TPC A side)
162 } else {
163 sign = -1; // (TPC C side)
164 }
165
166 if ( sign==1 && z < fgkZOffSet ) z = fgkZOffSet; // Protect against discontinuity at CE
167 if ( sign==-1 && z > -fgkZOffSet ) z = -fgkZOffSet; // Protect against discontinuity at CE
168
169
170 if ( (sign==1 && z<0) || (sign==-1 && z>0) ) // just a consistency check
171 AliError("ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!");
172
173 // Get the Er and Ephi field integrals plus the integral over DeltaEz
174 intEr = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
175 fgkRList, fgkZList, fgkPhiList, fLookUpErOverEz );
176 intEphi = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
177 fgkRList, fgkZList, fgkPhiList, fLookUpEphiOverEz);
178 intdEz = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
179 fgkRList, fgkZList, fgkPhiList, fLookUpDeltaEz );
180
181 // Calculate distorted position
182 if ( r > 0.0 ) {
183 phi = phi + fCorrectionFactor *( fC0*intEphi - fC1*intEr ) / r;
184 r = r + fCorrectionFactor *( fC0*intEr + fC1*intEphi );
185 }
186 Double_t dz = intdEz * fCorrectionFactor * fgkdvdE;
187
188 // Calculate correction in cartesian coordinates
189 dx[0] = r * TMath::Cos(phi) - x[0];
190 dx[1] = r * TMath::Sin(phi) - x[1];
191 dx[2] = dz; // z distortion - (scaled with driftvelocity dependency on the Ez field and the overall scaling factor)
192
193}
194
195
196void AliTPCSpaceCharge3D::InitSpaceCharge3DDistortion() {
197 //
198 // Initialization of the Lookup table which contains the solutions of the
199 // "space charge" (poisson) problem
200 //
201 // The sum-up uses a look-up table which contains different discretized Space charge fields
202 // in order to calculate the corresponding field deviations due to a given (discretized)
203 // space charge distribution ....
204 //
205 // Method of calculation: Weighted sum up of the different fields within the look up table
206 //
207
208 if (fInitLookUp) {
209 AliInfo("Lookup table was already initialized!");
210 // return;
211 }
212
213 AliInfo("Preparation of the weighted look-up table");
214
215 TFile *f = new TFile(fSCLookUpPOCsFileName);
216 if (!f) {
217 AliError("Precalculated POC-looup-table could not be found");
218 return;
219 }
220
221 // units are in [m]
222 TVector *gridf = (TVector*) f->Get("constants");
223 TVector &grid = *gridf;
224 TMatrix *coordf = (TMatrix*) f->Get("coordinates");
225 TMatrix &coord = *coordf;
226 TMatrix *coordPOCf = (TMatrix*) f->Get("POCcoord");
227 TMatrix &coordPOC = *coordPOCf;
228
229 Bool_t flagRadSym = 0;
230 if (grid(1)==1 && grid(4)==1) {
231 AliInfo("LOOK UP TABLE IS RADIAL SYMETTRIC - Field in Phi is ZERO");
232 flagRadSym=1;
233 }
234
235 Int_t rows = (Int_t)grid(0); // number of points in r direction
236 Int_t phiSlices = (Int_t)grid(1); // number of points in phi
237 Int_t columns = (Int_t)grid(2); // number of points in z direction
238
239 const Float_t gridSizeR = grid(6)*100; // unit in [cm]
240 const Float_t gridSizePhi = grid(7); // unit in [rad]
241 const Float_t gridSizeZ = grid(8)*100; // unit in [cm]
242
243 // temporary matrices needed for the calculation
244 TMatrixD *arrayofErA[phiSlices], *arrayofEphiA[phiSlices], *arrayofdEzA[phiSlices];
245 TMatrixD *arrayofErC[phiSlices], *arrayofEphiC[phiSlices], *arrayofdEzC[phiSlices];
246
247 TMatrixD *arrayofEroverEzA[phiSlices], *arrayofEphioverEzA[phiSlices], *arrayofDeltaEzA[phiSlices];
248 TMatrixD *arrayofEroverEzC[phiSlices], *arrayofEphioverEzC[phiSlices], *arrayofDeltaEzC[phiSlices];
249
250
251 for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
252
253 arrayofErA[k] = new TMatrixD(rows,columns) ;
254 arrayofEphiA[k] = new TMatrixD(rows,columns) ; // zero if radial symmetric
255 arrayofdEzA[k] = new TMatrixD(rows,columns) ;
256 arrayofErC[k] = new TMatrixD(rows,columns) ;
257 arrayofEphiC[k] = new TMatrixD(rows,columns) ; // zero if radial symmetric
258 arrayofdEzC[k] = new TMatrixD(rows,columns) ;
259
260 arrayofEroverEzA[k] = new TMatrixD(rows,columns) ;
261 arrayofEphioverEzA[k] = new TMatrixD(rows,columns) ; // zero if radial symmetric
262 arrayofDeltaEzA[k] = new TMatrixD(rows,columns) ;
263 arrayofEroverEzC[k] = new TMatrixD(rows,columns) ;
264 arrayofEphioverEzC[k] = new TMatrixD(rows,columns) ; // zero if radial symmetric
265 arrayofDeltaEzC[k] = new TMatrixD(rows,columns) ;
266
267 // Set the values to zero the lookup tables
268 // not necessary, it is done in the constructor of TMatrix - code deleted
269
270 }
271
272 // list of points as used in the interpolation (during sum up)
273 Double_t rlist[rows], zedlist[columns] , philist[phiSlices];
274 for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
275 philist[k] = gridSizePhi * k;
276 for ( Int_t i = 0 ; i < rows ; i++ ) {
277 rlist[i] = fgkIFCRadius + i*gridSizeR ;
278 for ( Int_t j = 0 ; j < columns ; j++ ) {
279 zedlist[j] = j * gridSizeZ ;
280 }
281 }
282 } // only done once
283
284
285 // Prepare summation Vectors
286 // Int_t numPosInd = (Int_t)(grid(0)*grid(1)*grid(2)); // Number of fulcrums in the look-up table
287 // Int_t numPOC = (Int_t)(grid(3)*grid(4)*grid(5)); // Number of POC conf. in the look-up table
288
289 TTree *treePOC = (TTree*)f->Get("POCall");
290
291 TVector *bEr = 0; TVector *bEphi= 0; TVector *bEz = 0;
292
293 treePOC->SetBranchAddress("Er",&bEr);
294 if (!flagRadSym) treePOC->SetBranchAddress("Ephi",&bEphi);
295 treePOC->SetBranchAddress("Ez",&bEz);
296
297 /* TBranch *b2Er = treePOC->GetBranch("Er");
298 TBranch *b2Ephi =0;
299 if (!flagRadSym) b2Ephi = treePOC->GetBranch("Ephi");
300 TBranch *b2Ez = treePOC->GetBranch("Ez");
301 */
302
303 // Read the complete tree and do a weighted sum-up over the POC configurations
304 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
305
306 Int_t treeNumPOC = (Int_t)treePOC->GetEntries(); // Number of POC conf. in the look-up table
307 Int_t ipC = 0; // POC Conf. counter (note: different to the POC number in the tree!)
308
309 for (Int_t itreepC=0; itreepC<treeNumPOC; itreepC++) { // loop over POC configurations in tree
310
311 // if (!flagRadSym) cout<<ipC<<endl;
312
313 treePOC->GetEntry(itreepC);
314
315 /* b2Er->GetEntry(itreepC);
316 if (!flagRadSym) b2Ephi->GetEntry(itreepC);
317 b2Ez->GetEntry(itreepC);
318 */
319
320 // center of the POC voxel in [meter]
321 Double_t r0 = coordPOC(ipC,0);
322 Double_t phi0 = coordPOC(ipC,1);
323 Double_t z0 = coordPOC(ipC,2);
324
325 ipC++; // POC conf. counter
326
327 // weights (charge density) at POC position on the A and C side (in C/m^3/e0)
328 // note: coordinates are in [cm]
329 Double_t weightA = GetSpaceChargeDensity(r0*100,phi0, z0*100);
330 Double_t weightC = GetSpaceChargeDensity(r0*100,phi0,-z0*100);
331
332 // Summing up the vector components according to their weight
333
334 Int_t ip = 0;
335 for ( Int_t j = 0 ; j < columns ; j++ ) {
336 for ( Int_t i = 0 ; i < rows ; i++ ) {
337 for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
338
339 // check wether the coordinates were screwed
340 if ((coord(0,ip)*100-rlist[i])>0.01 ||
341 (coord(1,ip)-philist[k])>0.01 ||
342 (coord(2,ip)*100-zedlist[j])>0.01) {
343 AliError("internal error: coordinate system was screwed during the sum-up");
344 printf("lookup: (r,phi,z)=(%lf,%lf,%lf)\n",coord(0,ip)*100,coord(1,ip),coord(2,ip)*100);
345 printf("sum-up: (r,phi,z)=(%lf,%lf,%lf)\n",rlist[i],philist[k],zedlist[j]);
346 AliError("Don't trust the results of the space charge calibration!");
347 }
348
349 // unfortunately, the lookup tables were produced to be faster for phi symmetric charges
350 // This will be the most frequent usage (hopefully)
351 // That's why we have to do this here ...
352
353 TMatrixD &erA = *arrayofErA[k] ;
354 TMatrixD &ephiA = *arrayofEphiA[k];
355 TMatrixD &dEzA = *arrayofdEzA[k] ;
356
357 TMatrixD &erC = *arrayofErC[k] ;
358 TMatrixD &ephiC = *arrayofEphiC[k];
359 TMatrixD &dEzC = *arrayofdEzC[k] ;
360
361 // Sum up - Efield values in [V/m] -> transition to [V/cm]
362 erA(i,j) += ((*bEr)(ip)) * weightA /100;
363 erC(i,j) += ((*bEr)(ip)) * weightC /100;
364 if (!flagRadSym) {
365 ephiA(i,j) += ((*bEphi)(ip)) * weightA; // [V/rad]
366 ephiC(i,j) += ((*bEphi)(ip)) * weightC; // [V/rad]
367 }
368 dEzA(i,j) += ((*bEz)(ip)) * weightA /100;
369 dEzC(i,j) += ((*bEz)(ip)) * weightC /100;
370
371 // increase the counter
372 ip++;
373 }
374 }
375 }
376
377 // Rotation and summation in the rest of the dPhiSteps
378 // which were not stored in the tree due to storage & symmetry reasons
379
380 Int_t phiPoints = (Int_t) grid(1);
381 Int_t phiPOC = (Int_t) grid(4);
382
383 // printf("%d %d\n",phiPOC,flagRadSym);
384
385 for (Int_t phiiC = 1; phiiC<phiPOC; phiiC++) { // just used for non-radial symetric table
386
387 r0 = coordPOC(ipC,0);
388 phi0 = coordPOC(ipC,1);
389 z0 = coordPOC(ipC,2);
390
391 ipC++; // POC conf. counter
392
393 // weights (charge density) at POC position on the A and C side (in C/m^3/e0)
394 // note: coordinates are in [cm]
395 weightA = GetSpaceChargeDensity(r0*100,phi0, z0*100);
396 weightC = GetSpaceChargeDensity(r0*100,phi0,-z0*100);
397
398
399 // Summing up the vector components according to their weight
400 ip = 0;
401 for ( Int_t j = 0 ; j < columns ; j++ ) {
402 for ( Int_t i = 0 ; i < rows ; i++ ) {
403 for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
404
405 // Note: rotating the coordinated during the sum up
406
407 Int_t rotVal = (phiPoints/phiPOC)*phiiC;
408 Int_t ipR = -1;
409
410 if ((ip%phiPoints)>=rotVal) {
411 ipR = ip-rotVal;
412 } else {
413 ipR = ip+(phiPoints-rotVal);
414 }
415
416 // unfortunately, the lookup tables were produced to be faster for phi symmetric charges
417 // This will be the most frequent usage (hopefully)
418 // That's why we have to do this here and not outside the loop ...
419
420 TMatrixD &erA = *arrayofErA[k] ;
421 TMatrixD &ephiA = *arrayofEphiA[k];
422 TMatrixD &dEzA = *arrayofdEzA[k] ;
423
424 TMatrixD &erC = *arrayofErC[k] ;
425 TMatrixD &ephiC = *arrayofEphiC[k];
426 TMatrixD &dEzC = *arrayofdEzC[k] ;
427
428 // Sum up - Efield values in [V/m] -> transition to [V/cm]
429 erA(i,j) += ((*bEr)(ipR)) * weightA /100;
430 erC(i,j) += ((*bEr)(ipR)) * weightC /100;
431 if (!flagRadSym) {
432 ephiA(i,j) += ((*bEphi)(ipR)) * weightA; // [V/rad]
433 ephiC(i,j) += ((*bEphi)(ipR)) * weightC; // [V/rad]
434 }
435 dEzA(i,j) += ((*bEz)(ipR)) * weightA /100;
436 dEzC(i,j) += ((*bEz)(ipR)) * weightC /100;
437
438 // increase the counter
439 ip++;
440 }
441 }
442 } // end coordinate loop
443
444 } // end phi-POC summation (phiiC)
445
446
447 // printf("POC: (r,phi,z) = (%lf %lf %lf) | weight(A,C): %03.1lf %03.1lf\n",r0,phi0,z0, weightA, weightC);
448
449 } // end POC loop
450
451 AliInfo("Division and integration");
452
453 // -------------------------------------------------------------------------------
454 // Division by the Ez (drift) field and integration along z
455
456 Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = Electric Field (V/cm) Magnitude ~ -400 V/cm;
457
458 for ( Int_t k = 0 ; k < phiSlices ; k++ ) { // phi loop
459
460 // matrices holding the solution - summation of POC charges
461 TMatrixD &erA = *arrayofErA[k] ;
462 TMatrixD &ephiA = *arrayofEphiA[k];
463 TMatrixD &dezA = *arrayofdEzA[k] ;
464 TMatrixD &erC = *arrayofErC[k] ;
465 TMatrixD &ephiC = *arrayofEphiC[k];
466 TMatrixD &dezC = *arrayofdEzC[k] ;
467
468 // matrices which contain the integrated fields (divided by the drift field)
469 TMatrixD &erOverEzA = *arrayofEroverEzA[k] ;
470 TMatrixD &ephiOverEzA = *arrayofEphioverEzA[k];
471 TMatrixD &deltaEzA = *arrayofDeltaEzA[k];
472 TMatrixD &erOverEzC = *arrayofEroverEzC[k] ;
473 TMatrixD &ephiOverEzC = *arrayofEphioverEzC[k];
474 TMatrixD &deltaEzC = *arrayofDeltaEzC[k];
475
476 for ( Int_t i = 0 ; i < rows ; i++ ) { // r loop
477 for ( Int_t j = columns-1 ; j >= 0 ; j-- ) {// z loop
478 // Count backwards to facilitate integration over Z
479
480 Int_t index = 1 ; // Simpsons rule if N=odd.If N!=odd then add extra point by trapezoidal rule.
481
482 erOverEzA(i,j) = 0; ephiOverEzA(i,j) = 0; deltaEzA(i,j) = 0;
483 erOverEzC(i,j) = 0; ephiOverEzC(i,j) = 0; deltaEzC(i,j) = 0;
484
485 for ( Int_t m = j ; m < columns ; m++ ) { // integration
486
487 erOverEzA(i,j) += index*(gridSizeZ/3.0)*erA(i,m)/(-1*ezField) ;
488 erOverEzC(i,j) += index*(gridSizeZ/3.0)*erC(i,m)/(-1*ezField) ;
489 if (!flagRadSym) {
490 ephiOverEzA(i,j) += index*(gridSizeZ/3.0)*ephiA(i,m)/(-1*ezField) ;
491 ephiOverEzC(i,j) += index*(gridSizeZ/3.0)*ephiC(i,m)/(-1*ezField) ;
492 }
493 deltaEzA(i,j) += index*(gridSizeZ/3.0)*dezA(i,m) ;
494 deltaEzC(i,j) += index*(gridSizeZ/3.0)*dezC(i,m) ;
495
496 if ( index != 4 ) index = 4; else index = 2 ;
497
498 }
499
500 if ( index == 4 ) {
501 erOverEzA(i,j) -= (gridSizeZ/3.0)*erA(i,columns-1)/(-1*ezField) ;
502 erOverEzC(i,j) -= (gridSizeZ/3.0)*erC(i,columns-1)/(-1*ezField) ;
503 if (!flagRadSym) {
504 ephiOverEzA(i,j) -= (gridSizeZ/3.0)*ephiA(i,columns-1)/(-1*ezField) ;
505 ephiOverEzC(i,j) -= (gridSizeZ/3.0)*ephiC(i,columns-1)/(-1*ezField) ;
506 }
507 deltaEzA(i,j) -= (gridSizeZ/3.0)*dezA(i,columns-1) ;
508 deltaEzC(i,j) -= (gridSizeZ/3.0)*dezC(i,columns-1) ;
509 }
510 if ( index == 2 ) {
511 erOverEzA(i,j) += (gridSizeZ/3.0)*(0.5*erA(i,columns-2)-2.5*erA(i,columns-1))/(-1*ezField) ;
512 erOverEzC(i,j) += (gridSizeZ/3.0)*(0.5*erC(i,columns-2)-2.5*erC(i,columns-1))/(-1*ezField) ;
513 if (!flagRadSym) {
514 ephiOverEzA(i,j) += (gridSizeZ/3.0)*(0.5*ephiA(i,columns-2)-2.5*ephiA(i,columns-1))/(-1*ezField) ;
515 ephiOverEzC(i,j) += (gridSizeZ/3.0)*(0.5*ephiC(i,columns-2)-2.5*ephiC(i,columns-1))/(-1*ezField) ;
516 }
517 deltaEzA(i,j) += (gridSizeZ/3.0)*(0.5*dezA(i,columns-2)-2.5*dezA(i,columns-1)) ;
518 deltaEzC(i,j) += (gridSizeZ/3.0)*(0.5*dezC(i,columns-2)-2.5*dezC(i,columns-1)) ;
519 }
520 if ( j == columns-2 ) {
521 erOverEzA(i,j) = (gridSizeZ/3.0)*(1.5*erA(i,columns-2)+1.5*erA(i,columns-1))/(-1*ezField) ;
522 erOverEzC(i,j) = (gridSizeZ/3.0)*(1.5*erC(i,columns-2)+1.5*erC(i,columns-1))/(-1*ezField) ;
523 if (!flagRadSym) {
524 ephiOverEzA(i,j) = (gridSizeZ/3.0)*(1.5*ephiA(i,columns-2)+1.5*ephiA(i,columns-1))/(-1*ezField) ;
525 ephiOverEzC(i,j) = (gridSizeZ/3.0)*(1.5*ephiC(i,columns-2)+1.5*ephiC(i,columns-1))/(-1*ezField) ;
526 }
527 deltaEzA(i,j) = (gridSizeZ/3.0)*(1.5*dezA(i,columns-2)+1.5*dezA(i,columns-1)) ;
528 deltaEzC(i,j) = (gridSizeZ/3.0)*(1.5*dezC(i,columns-2)+1.5*dezC(i,columns-1)) ;
529 }
530 if ( j == columns-1 ) {
531 erOverEzA(i,j) = 0; erOverEzC(i,j) = 0;
532 if (!flagRadSym) {
533 ephiOverEzA(i,j) = 0; ephiOverEzC(i,j) = 0;
534 }
535 deltaEzA(i,j) = 0; deltaEzC(i,j) = 0;
536 }
537 }
538 }
539
540 }
541
542 AliInfo("Interpolation to Standard grid");
543
544 // -------------------------------------------------------------------------------
545 // Interpolate results onto the standard grid which is used for all AliTPCCorrections
546
547 const Int_t order = 1 ; // Linear interpolation = 1, Quadratic = 2
548
549 Double_t r, phi, z ;
550 for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
551 phi = fgkPhiList[k] ;
552
553 TMatrixD &erOverEz = *fLookUpErOverEz[k] ;
554 TMatrixD &ephiOverEz = *fLookUpEphiOverEz[k];
555 TMatrixD &deltaEz = *fLookUpDeltaEz[k] ;
556
557 for ( Int_t j = 0 ; j < kNZ ; j++ ) {
558
559 z = TMath::Abs(fgkZList[j]) ; // z position is symmetric
560
561 for ( Int_t i = 0 ; i < kNR ; i++ ) {
562 r = fgkRList[i] ;
563
564 // Interpolate Lookup tables onto standard grid
565 if (fgkZList[j]>0) {
566 erOverEz(i,j) = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
567 rlist, zedlist, philist, arrayofEroverEzA );
568 ephiOverEz(i,j) = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
569 rlist, zedlist, philist, arrayofEphioverEzA);
570 deltaEz(i,j) = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
571 rlist, zedlist, philist, arrayofDeltaEzA );
572 } else {
573 erOverEz(i,j) = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
574 rlist, zedlist, philist, arrayofEroverEzC );
575 ephiOverEz(i,j) = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
576 rlist, zedlist, philist, arrayofEphioverEzC);
577 deltaEz(i,j) = - Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
578 rlist, zedlist, philist, arrayofDeltaEzC );
579 // negative coordinate system on C side
580 }
581
582 } // end r loop
583 } // end z loop
584 } // end phi loop
585
586
587 // clear the temporary arrays lists
588 for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
589
590 delete arrayofErA[k];
591 delete arrayofEphiA[k];
592 delete arrayofdEzA[k];
593 delete arrayofErC[k];
594 delete arrayofEphiC[k];
595 delete arrayofdEzC[k];
596
597 delete arrayofEroverEzA[k];
598 delete arrayofEphioverEzA[k];
599 delete arrayofDeltaEzA[k];
600 delete arrayofEroverEzC[k];
601 delete arrayofEphioverEzC[k];
602 delete arrayofDeltaEzC[k];
603
604 }
605
606
607 fInitLookUp = kTRUE;
608
609}
610
611
612void AliTPCSpaceCharge3D::SetSCDataFileName(char *const fname) {
613 //
614 // Set & load the Space charge density distribution from a file
615 // (linear interpolation onto a standard grid)
616 //
617
618 fSCDataFileName = fname;
619
620 TFile *f = new TFile(fSCDataFileName,"READ");
621 if (!f) {
622 AliError(Form("File %s, which should contain the space charge distribution, could not be found",
623 fSCDataFileName));
624 return;
625 }
626
627 TH3F *density = (TH3F*) f->Get("SpaceChargeDistribution");
628 if (!density) {
629 AliError(Form("The indicated file (%s) does not contain a histogram called %s",
630 fSCDataFileName,"SpaceChargeDistribution"));
631 return;
632 }
633
634
635 Double_t r, phi, z ;
636 for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
637 phi = fgkPhiList[k] ;
638 TMatrixD &scdensity = *fSCdensityDistribution[k] ;
639 for ( Int_t j = 0 ; j < kNZ ; j++ ) {
640 z = TMath::Abs(fgkZList[j]) ; // z position is symmetric
641 for ( Int_t i = 0 ; i < kNR ; i++ ) {
642 r = fgkRList[i] ;
643 scdensity(i,j) = density->Interpolate(r,phi,z); // quite slow
644 // printf("%lf %lf %lf: %lf\n",r,phi,z,scdensity(i,j));
645 }
646 }
647 }
648
649
650 f->Close();
651
652 fInitLookUp = kFALSE;
653
654
655}
656
657
658Float_t AliTPCSpaceCharge3D::GetSpaceChargeDensity(Float_t r, Float_t phi, Float_t z) {
659 //
660 // returns the (input) space charge density at a given point according
661 // Note: input in [cm], output in [C/m^3/e0] !!
662 //
663
664 if (!fSCdensityDistribution) {
665 printf("Scheiss is nuul - argg\n");
666 return 0.;
667 }
668
669 while (phi<0) phi += TMath::TwoPi();
670 while (phi>TMath::TwoPi()) phi -= TMath::TwoPi();
671
672
673 // Float_t sc =fSCdensityDistribution->Interpolate(r0,phi0,z0);
674 Int_t order = 1; //
675 Float_t sc = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
676 fgkRList, fgkZList, fgkPhiList, fSCdensityDistribution );
677
678 // printf("%lf %lf %lf: %lf\n",r,phi,z,sc);
679
680 return sc;
681}
682
683
684TH2F * AliTPCSpaceCharge3D::CreateHistoSCinXY(Float_t z, Int_t nx, Int_t ny) {
685 //
686 // return a simple histogramm containing the space charge distribution (input for the calculation)
687 //
688
689 TH2F *h=CreateTH2F("spaceCharge",GetTitle(),"x [cm]","y [cm]","#rho_{sc} [C/m^{3}/e_{0}]",
690 nx,-250.,250.,ny,-250.,250.);
691
692 for (Int_t iy=1;iy<=ny;++iy) {
693 Double_t yp = h->GetYaxis()->GetBinCenter(iy);
694 for (Int_t ix=1;ix<=nx;++ix) {
695 Double_t xp = h->GetXaxis()->GetBinCenter(ix);
696
697 Float_t r = TMath::Sqrt(xp*xp+yp*yp);
698 Float_t phi = TMath::ATan2(yp,xp);
699
700 if (85.<=r && r<=250.) {
701 Float_t sc = GetSpaceChargeDensity(r,phi,z)/fgke0; // in [C/m^3/e0]
702 h->SetBinContent(ix,iy,sc);
703 } else {
704 h->SetBinContent(ix,iy,0.);
705 }
706 }
707 }
708
709 return h;
710}
711
712TH2F * AliTPCSpaceCharge3D::CreateHistoSCinZR(Float_t phi, Int_t nz, Int_t nr) {
713 //
714 // return a simple histogramm containing the space charge distribution (input for the calculation)
715 //
716
717 TH2F *h=CreateTH2F("spaceCharge",GetTitle(),"z [cm]","r [cm]","#rho_{sc} [C/m^{3}/e_{0}]",
718 nz,-250.,250.,nr,85.,250.);
719
720 for (Int_t ir=1;ir<=nr;++ir) {
721 Float_t r = h->GetYaxis()->GetBinCenter(ir);
722 for (Int_t iz=1;iz<=nz;++iz) {
723 Float_t z = h->GetXaxis()->GetBinCenter(iz);
724 Float_t sc = GetSpaceChargeDensity(r,phi,z)/fgke0; // in [C/m^3/e0]
725 h->SetBinContent(iz,ir,sc);
726 }
727 }
728
729 return h;
730}
731
732void AliTPCSpaceCharge3D::WriteChargeDistributionToFile(const char* fname) {
733 //
734 // Example on how to write a Space charge distribution into a File
735 // (see below: estimate from scaling STAR measurements to Alice)
736 //
737
738 TFile *f = new TFile(fname,"RECREATE");
739
740 // some grid, not too course
741 Int_t nr = 72;
742 Int_t nphi = 180;
743 Int_t nz = 250;
744
745 Double_t dr = (fgkOFCRadius-fgkIFCRadius)/(nr+1);
746 Double_t dphi = TMath::TwoPi()/(nphi+1);
747 Double_t dz = 500./(nz+1);
748 Double_t safty = 0.; // due to a root bug which does not interpolate the boundary (first and last bin) correctly
749
750 TH3F *histo = new TH3F("charge","charge",
751 nr,fgkIFCRadius-dr-safty,fgkOFCRadius+dr+safty,
752 nphi,0-dphi-safty,TMath::TwoPi()+dphi+safty,
753 nz,-250-dz-safty,250+dz+safty);
754
755 for (Int_t ir=1;ir<=nr;++ir) {
756 Double_t rp = histo->GetXaxis()->GetBinCenter(ir);
757 for (Int_t iphi=1;iphi<=nphi;++iphi) {
758 Double_t phip = histo->GetYaxis()->GetBinCenter(iphi);
759 for (Int_t iz=1;iz<=nz;++iz) {
760 Double_t zp = histo->GetZaxis()->GetBinCenter(iz);
761
762
763 zp = TMath::Abs(zp); // symmetric on both sides of the TPC
764
765 // recalculation to meter
766 Double_t lZ = 2.5; // approx. TPC drift length
767 Double_t rpM = rp/100.; // in [m]
768 Double_t zpM = zp/100.; // in [m]
769
770 // setting of mb multiplicity and Interaction rate
771 Double_t multiplicity = 950;
772 Double_t intRate = 8000;
773
774 // calculation of "scaled" parameters
775 Double_t a = multiplicity*intRate/79175;
776 Double_t b = a/lZ;
777
778 Double_t charge = (a - b * zpM)/(rpM*rpM); // charge in [C/m^3/e0]
779
780 // some 'arbitrary' GG leaks
781 if ( (phip<(TMath::Pi()/9*3) && phip>(TMath::Pi()/9*2) ) ) { //||
782 // (phip<(TMath::Pi()/9*7) && phip>(TMath::Pi()/9*6) ) ||
783 // (phip<(TMath::Pi()/9*16) && phip>(TMath::Pi()/9*13) ) ){
784 if (rp>130&&rp<135)
785 charge = 300;
786 }
787
788 // printf("%lf %lf %lf: %lf\n",rp,phip,zp,charge);
789
790 charge = charge*fgke0; // [C/m^3]
791
792 histo->SetBinContent(ir,iphi,iz,charge);
793 }
794 }
795 }
796
797
798 histo->Write("SpaceChargeDistribution");
799 f->Close();
800
801}
802
803
804void AliTPCSpaceCharge3D::Print(const Option_t* option) const {
805 //
806 // Print function to check the settings of the boundary vectors
807 // option=="a" prints the C0 and C1 coefficents for calibration purposes
808 //
809
810 TString opt = option; opt.ToLower();
811 printf("%s\n",GetTitle());
812 printf(" - Space Charge effect with arbitrary 3D charge density (as input).\n");
813 printf(" SC correction factor: %f \n",fCorrectionFactor);
814
815 if (opt.Contains("a")) { // Print all details
816 printf(" - T1: %1.4f, T2: %1.4f \n",fT1,fT2);
817 printf(" - C1: %1.4f, C0: %1.4f \n",fC1,fC0);
818 }
819
820 if (!fInitLookUp) AliError("Lookup table was not initialized! You should do InitSpaceCharge3DDistortion() ...");
821
822}