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