]>
Commit | Line | Data |
---|---|---|
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 | ||
53 | ClassImp(AliTPCSpaceCharge3D) | |
54 | ||
55 | AliTPCSpaceCharge3D::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 | ||
82 | AliTPCSpaceCharge3D::~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 | ||
99 | void 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 | ||
118 | void 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 | ||
139 | void 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 | ||
196 | void 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 | ||
612 | void 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 | ||
658 | Float_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 | ||
684 | TH2F * 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 | ||
712 | TH2F * 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 | ||
732 | void 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 | ||
804 | void 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 | } |