Added protection, fixed dynamic cast
[u/mrichter/AliRoot.git] / TPC / AliTPCFCVoltError3D.cxx
CommitLineData
c9cbd2f2 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// AliTPCFCVoltError3D class //
19// The class calculates the space point distortions due to residual voltage //
20// errors on the Field Cage boundaries (rods) of the TPC in 3D. //
21// //
22// The class allows "effective Omega Tau" corrections. //
23// //
24// NOTE: This class is capable of calculating z distortions due to //
25// drift velocity changes in dependence of the electric field!!! //
26// //
27// date: 08/08/2010 //
28// Authors: Jim Thomas, Stefan Rossegger //
29// //
30// Example usage : //
31// AliTPCFCVoltError3D fcerror; //
32//////////////////////////////////////////////////////////////////////////////
33
34#include "AliMagF.h"
35#include "TGeoGlobalMagField.h"
36#include "AliTPCcalibDB.h"
37#include "AliTPCParam.h"
38#include "AliLog.h"
39#include "TMatrixD.h"
40
41#include "TMath.h"
42#include "AliTPCROC.h"
43#include "AliTPCFCVoltError3D.h"
44
45ClassImp(AliTPCFCVoltError3D)
46
47AliTPCFCVoltError3D::AliTPCFCVoltError3D()
48 : AliTPCCorrection("FieldCageVoltErrors","FieldCage (Rods) Voltage Errors"),
49 fC0(0.),fC1(0.),
50 fInitLookUp(kFALSE)
51{
52 //
53 // default constructor
54 //
55
56 // flags for filled 'basic lookup tables'
57 for (Int_t i=0; i<5; i++){
58 fInitLookUpBasic[i]= kFALSE;
59 }
60
61 // Boundary settings
62 for (Int_t i=0; i<36; i++){
63 fRodVoltShiftA[i] = 0;
64 fRodVoltShiftC[i] = 0;
65 }
66 for (Int_t i=0; i<2; i++){
67 fRotatedClipVoltA[i] = 0;
68 fRotatedClipVoltC[i] = 0;
69 }
70 //
71 for (Int_t i=0; i<18; i++){
72 fOFCRodShiftA[i] = 0;
73 fOFCRodShiftC[i] = 0;
74 }
75
76 // Array which will contain the solution according to the setted boundary conditions
77 // it represents a sum up of the 4 basic look up tables (created individually)
78 // see InitFCVoltError3D() function
79 for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
80 fLookUpErOverEz[k] = new TMatrixD(kNR,kNZ);
81 fLookUpEphiOverEz[k] = new TMatrixD(kNR,kNZ);
82 fLookUpDeltaEz[k] = new TMatrixD(kNR,kNZ);
83 }
84
85 for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
86 fLookUpBasic1ErOverEz[k] = 0;
87 fLookUpBasic1EphiOverEz[k] = 0;
88 fLookUpBasic1DeltaEz[k] = 0;
89
90 fLookUpBasic2ErOverEz[k] = 0;
91 fLookUpBasic2EphiOverEz[k] = 0;
92 fLookUpBasic2DeltaEz[k] = 0;
93
94 fLookUpBasic3ErOverEz[k] = 0;
95 fLookUpBasic3EphiOverEz[k] = 0;
96 fLookUpBasic3DeltaEz[k] = 0;
97
98 fLookUpBasic4ErOverEz[k] = 0;
99 fLookUpBasic4EphiOverEz[k] = 0;
100 fLookUpBasic4DeltaEz[k] = 0;
101
102 fLookUpBasic5ErOverEz[k] = 0;
103 fLookUpBasic5EphiOverEz[k] = 0;
104 fLookUpBasic5DeltaEz[k] = 0;
105 }
106
107}
108
109AliTPCFCVoltError3D::~AliTPCFCVoltError3D() {
110 //
111 // destructor
112 //
113
114 for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
115 delete fLookUpErOverEz[k];
116 delete fLookUpEphiOverEz[k];
117 delete fLookUpDeltaEz[k];
118 }
119
120 for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
121 delete fLookUpBasic1ErOverEz[k]; // does nothing if pointer is zero!
122 delete fLookUpBasic1EphiOverEz[k];
123 delete fLookUpBasic1DeltaEz[k];
124
125 delete fLookUpBasic2ErOverEz[k]; // does nothing if pointer is zero!
126 delete fLookUpBasic2EphiOverEz[k];
127 delete fLookUpBasic2DeltaEz[k];
128
129 delete fLookUpBasic3ErOverEz[k]; // does nothing if pointer is zero!
130 delete fLookUpBasic3EphiOverEz[k];
131 delete fLookUpBasic3DeltaEz[k];
132
133 delete fLookUpBasic4ErOverEz[k]; // does nothing if pointer is zero!
134 delete fLookUpBasic4EphiOverEz[k];
135 delete fLookUpBasic4DeltaEz[k];
136
137 delete fLookUpBasic5ErOverEz[k]; // does nothing if pointer is zero!
138 delete fLookUpBasic5EphiOverEz[k];
139 delete fLookUpBasic5DeltaEz[k];
140 }
141}
142
143void AliTPCFCVoltError3D::Init() {
144 //
145 // Initialization funtion
146 //
147
148 AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
149 if (!magF) AliError("Magneticd field - not initialized");
150 Double_t bzField = magF->SolenoidField()/10.; //field in T
151 AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
152 if (!param) AliError("Parameters - not initialized");
153 Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us] // From dataBase: to be updated: per second (ideally)
154 Double_t ezField = 400; // [V/cm] // to be updated: never (hopefully)
155 Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
156 // Correction Terms for effective omegaTau; obtained by a laser calibration run
157 SetOmegaTauT1T2(wt,fT1,fT2);
158
159 InitFCVoltError3D();
160}
161
162void AliTPCFCVoltError3D::Update(const TTimeStamp &/*timeStamp*/) {
163 //
164 // Update function
165 //
166 AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
167 if (!magF) AliError("Magneticd field - not initialized");
168 Double_t bzField = magF->SolenoidField()/10.; //field in T
169 AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
170 if (!param) AliError("Parameters - not initialized");
171 Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us] // From dataBase: to be updated: per second (ideally)
172 Double_t ezField = 400; // [V/cm] // to be updated: never (hopefully)
173 Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
174 // Correction Terms for effective omegaTau; obtained by a laser calibration run
175 SetOmegaTauT1T2(wt,fT1,fT2);
176
177
178}
179
180
181
182void AliTPCFCVoltError3D::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {
183 //
184 // Calculates the correction due e.g. residual voltage errors on the TPC boundaries
185 //
186
187 if (!fInitLookUp) {
188 AliInfo("Lookup table was not initialized! Perform the inizialisation now ...");
189 InitFCVoltError3D();
190 return;
191 }
192
193 Int_t order = 1 ; // FIXME: hardcoded? Linear interpolation = 1, Quadratic = 2
194 // note that the poisson solution was linearly mirroed on this grid!
195 Double_t intEr, intEphi, intDeltaEz;
196 Double_t r, phi, z ;
197 Int_t sign;
198
199 r = TMath::Sqrt( x[0]*x[0] + x[1]*x[1] ) ;
200 phi = TMath::ATan2(x[1],x[0]) ;
201 if ( phi < 0 ) phi += TMath::TwoPi() ; // Table uses phi from 0 to 2*Pi
202 z = x[2] ; // Create temporary copy of x[2]
203
204 if ( (roc%36) < 18 ) {
205 sign = 1; // (TPC A side)
206 } else {
207 sign = -1; // (TPC C side)
208 }
209
210 if ( sign==1 && z < fgkZOffSet ) z = fgkZOffSet; // Protect against discontinuity at CE
211 if ( sign==-1 && z > -fgkZOffSet ) z = -fgkZOffSet; // Protect against discontinuity at CE
212
213
214 if ( (sign==1 && z<0) || (sign==-1 && z>0) ) // just a consistency check
215 AliError("ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!");
216
217 // Get the Er and Ephi field integrals plus the integral over DeltaEz
218 intEr = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
219 fgkRList, fgkZList, fgkPhiList, fLookUpErOverEz );
220 intEphi = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
221 fgkRList, fgkZList, fgkPhiList, fLookUpEphiOverEz );
222 intDeltaEz = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
223 fgkRList, fgkZList, fgkPhiList, fLookUpDeltaEz );
224
225 // printf("%lf %lf %lf\n",intEr,intEphi,intDeltaEz);
226
227 // Calculate distorted position
228 if ( r > 0.0 ) {
229 phi = phi + ( fC0*intEphi - fC1*intEr ) / r;
230 r = r + ( fC0*intEr + fC1*intEphi );
231 }
232
233 // Calculate correction in cartesian coordinates
234 dx[0] = r * TMath::Cos(phi) - x[0];
235 dx[1] = r * TMath::Sin(phi) - x[1];
236 dx[2] = intDeltaEz; // z distortion - (internally scaled with driftvelocity dependency
237 // on the Ez field plus the actual ROC misalignment (if set TRUE)
238
239}
240
241void AliTPCFCVoltError3D::InitFCVoltError3D() {
242 //
243 // Initialization of the Lookup table which contains the solutions of the
244 // Dirichlet boundary problem
245 // Calculation of the single 3D-Poisson solver is done just if needed
246 // (see basic lookup tables in header file)
247 //
248
249 const Int_t order = 1 ; // Linear interpolation = 1, Quadratic = 2
250 const Float_t gridSizeR = (fgkOFCRadius-fgkIFCRadius) / (kRows-1) ;
251 const Float_t gridSizeZ = fgkTPCZ0 / (kColumns-1) ;
252 const Float_t gridSizePhi = TMath::TwoPi() / ( 18.0 * kPhiSlicesPerSector);
253
254 // temporary arrays to create the boundary conditions
255 TMatrixD *arrayofArrayV[kPhiSlices], *arrayofCharge[kPhiSlices] ;
256 for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
257 arrayofArrayV[k] = new TMatrixD(kRows,kColumns) ;
258 arrayofCharge[k] = new TMatrixD(kRows,kColumns) ;
259 }
260 Double_t innerList[kPhiSlices], outerList[kPhiSlices] ;
261
262 // list of point as used in the poisson relation and the interpolation (during sum up)
263 Double_t rlist[kRows], zedlist[kColumns] , philist[kPhiSlices];
264 for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
265 philist[k] = gridSizePhi * k;
266 for ( Int_t i = 0 ; i < kRows ; i++ ) {
267 rlist[i] = fgkIFCRadius + i*gridSizeR ;
268 for ( Int_t j = 0 ; j < kColumns ; j++ ) { // Fill Vmatrix with Boundary Conditions
269 zedlist[j] = j * gridSizeZ ;
270 }
271 }
272 }
273
274 // ==========================================================================
275 // Solve Poisson's equation in 3D cylindrical coordinates by relaxation technique
276 // Allow for different size grid spacing in R and Z directions
277
278 const Int_t symmetry[5] = {1,1,-1,-1,1}; // shifted rod (2x), rotated clip (2x), only rod shift on OFC (1x)
279
280 // check which lookup tables are needed ---------------------------------
281
282 Bool_t needTable[5] ={kFALSE,kFALSE,kFALSE,kFALSE,kFALSE};
283
284 // Shifted rods & strips
285 for ( Int_t rod = 0 ; rod < 18 ; rod++ ) {
286 if (fRodVoltShiftA[rod]!=0) needTable[0]=kTRUE;
287 if (fRodVoltShiftC[rod]!=0) needTable[0]=kTRUE;
288 }
289 for ( Int_t rod = 18 ; rod < 36 ; rod++ ) {
290 if (fRodVoltShiftA[rod]!=0) needTable[1]=kTRUE;
291 if (fRodVoltShiftC[rod]!=0) needTable[1]=kTRUE;
292 }
293 // Rotated clips
294 if (fRotatedClipVoltA[0]!=0 || fRotatedClipVoltC[0]!=0) needTable[2]=kTRUE;
295 if (fRotatedClipVoltA[1]!=0 || fRotatedClipVoltC[1]!=0) needTable[3]=kTRUE;
296
297 // shifted Copper rods on OFC
298 for ( Int_t rod = 0 ; rod < 18 ; rod++ ) {
299 if (fOFCRodShiftA[rod]!=0) needTable[4]=kTRUE;
300 if (fOFCRodShiftC[rod]!=0) needTable[4]=kTRUE;
301 }
302
303
304 // reserve the arrays for the basic lookup tables ----------------------
305 if (needTable[0] && !fInitLookUpBasic[0]) {
306 for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) { // Possibly make an extra table to be used for 0 == 360
307 fLookUpBasic1ErOverEz[k] = new TMatrixD(kRows,kColumns);
308 fLookUpBasic1EphiOverEz[k] = new TMatrixD(kRows,kColumns);
309 fLookUpBasic1DeltaEz[k] = new TMatrixD(kRows,kColumns);
310 // will be deleted in destructor
311 }
312 }
313 if (needTable[1] && !fInitLookUpBasic[1]) {
314 for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) { // Possibly make an extra table to be used for 0 == 360
315 fLookUpBasic2ErOverEz[k] = new TMatrixD(kRows,kColumns);
316 fLookUpBasic2EphiOverEz[k] = new TMatrixD(kRows,kColumns);
317 fLookUpBasic2DeltaEz[k] = new TMatrixD(kRows,kColumns);
318 // will be deleted in destructor
319 }
320 }
321 if (needTable[2] && !fInitLookUpBasic[2]) {
322 for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) { // Possibly make an extra table to be used for 0 == 360
323 fLookUpBasic3ErOverEz[k] = new TMatrixD(kRows,kColumns);
324 fLookUpBasic3EphiOverEz[k] = new TMatrixD(kRows,kColumns);
325 fLookUpBasic3DeltaEz[k] = new TMatrixD(kRows,kColumns);
326 // will be deleted in destructor
327 }
328 }
329 if (needTable[3] && !fInitLookUpBasic[3]) {
330 for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) { // Possibly make an extra table to be used for 0 == 360
331 fLookUpBasic4ErOverEz[k] = new TMatrixD(kRows,kColumns);
332 fLookUpBasic4EphiOverEz[k] = new TMatrixD(kRows,kColumns);
333 fLookUpBasic4DeltaEz[k] = new TMatrixD(kRows,kColumns);
334 // will be deleted in destructor
335 }
336 }
337 if (needTable[4] && !fInitLookUpBasic[4]) {
338 for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) { // Possibly make an extra table to be used for 0 == 360
339 fLookUpBasic5ErOverEz[k] = new TMatrixD(kRows,kColumns);
340 fLookUpBasic5EphiOverEz[k] = new TMatrixD(kRows,kColumns);
341 fLookUpBasic5DeltaEz[k] = new TMatrixD(kRows,kColumns);
342 // will be deleted in destructor
343 }
344 }
345
346 // Set bondaries and solve Poisson's equation --------------------------
347
348 for (Int_t look=0; look<5; look++) {
349
350 Float_t inner18[18] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } ;
351 Float_t outer18[18] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } ;
352
353 if (needTable[look] && !fInitLookUpBasic[look]) {
354
355 // Specify which rod is the reference; other distortions calculated by summing over multiple rotations of refrence
356 // Note: the interpolation later on depends on it!! Do not change if not really needed!
357 if (look==0) {
358 AliInfo(Form("IFC - ROD&Strip shift : Filling table (~ %d sec)",3*(int)(kPhiSlices/5)));
359 inner18[0] = 1;
360 } else if (look==1) {
361 AliInfo(Form("OFC - ROD&Strip shift : Filling table (~ %d sec)",3*(int)(kPhiSlices/5)));
362 outer18[0] = 1;
363 } else if (look==2) {
364 AliInfo(Form("IFC - Clip rot. : Filling table (~ %d sec)",3*(int)(kPhiSlices/5)));
365 inner18[0] = 1;
366 } else if (look==3) {
367 AliInfo(Form("OFC - Clip rot. : Filling table (~ %d sec)",3*(int)(kPhiSlices/5)));
368 outer18[0] = 1;
369 } else if (look==4) {
370 AliInfo(Form("OFC - CopperRod shift : Filling table (~ %d sec)",3*(int)(kPhiSlices/5)));
371 outer18[0] = 1;
372 }
373 // Build potentials on boundary stripes for a rotated clip (SYMMETRY==-1) or a shifted rod (SYMMETRY==1)
374 if (look!=4) {
375 // in these cases, the strips follow the actual rod shift (linear interpolation between the rods)
376 for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
377 Int_t slices = kPhiSlicesPerSector ;
378 Int_t ipoint = k/slices ;
379 innerList[k] = (((ipoint+1)*slices-k)*-1*inner18[ipoint]+(k-ipoint*slices)*inner18[ipoint+1])/slices ;
380 outerList[k] = (((ipoint+1)*slices-k)*-1*outer18[ipoint]+(k-ipoint*slices)*outer18[ipoint+1])/slices ;
381 if ( (k%slices) == 0 && symmetry[look] == -1 ) { innerList[k] = 0.0 ; outerList[k] = 0.0 ; }
382 // above, force Zero if Anti-Sym
383 }
384 } else {
385 // in this case, we assume that the strips stay at the correct position, but the copper plated OFC-rods move
386 // the distortion is then much more localized around the rod (indicated by real data)
387 for ( Int_t k = 0 ; k < kPhiSlices ; k++ )
388 innerList[k] = outerList[k] = 0;
389 outerList[0] = outer18[0]; // point at rod
390 outerList[1] = outer18[0]/4; // point close to rod (educated-`guessed` scaling)
391 }
392
393 // Fill arrays with initial conditions. V on the boundary and Charge in the volume.
394 for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
395 TMatrixD &arrayV = *arrayofArrayV[k] ;
396 TMatrixD &charge = *arrayofCharge[k] ;
397 for ( Int_t i = 0 ; i < kRows ; i++ ) {
398 for ( Int_t j = 0 ; j < kColumns ; j++ ) { // Fill Vmatrix with Boundary Conditions
399 arrayV(i,j) = 0.0 ;
400 charge(i,j) = 0.0 ;
401 if ( i == 0 ) arrayV(i,j) = innerList[k] ;
402 if ( i == kRows-1 ) arrayV(i,j) = outerList[k] ;
403 }
404 }
405 // no charge in the volume
406 for ( Int_t i = 1 ; i < kRows-1 ; i++ ) {
407 for ( Int_t j = 1 ; j < kColumns-1 ; j++ ) {
408 charge(i,j) = 0.0 ;
409 }
410 }
411 }
412
413 // Solve Poisson's equation in 3D cylindrical coordinates by relaxation technique
414 // Allow for different size grid spacing in R and Z directions
415 if (look==0) {
416 PoissonRelaxation3D( arrayofArrayV, arrayofCharge,
417 fLookUpBasic1ErOverEz, fLookUpBasic1EphiOverEz, fLookUpBasic1DeltaEz,
418 kRows, kColumns, kPhiSlices, gridSizePhi, kIterations, symmetry[0]) ;
419 AliInfo("IFC - ROD&Strip shift : done ");
420 } else if (look==1) {
421 PoissonRelaxation3D( arrayofArrayV, arrayofCharge,
422 fLookUpBasic2ErOverEz, fLookUpBasic2EphiOverEz, fLookUpBasic2DeltaEz,
423 kRows, kColumns, kPhiSlices, gridSizePhi, kIterations, symmetry[1]) ;
424
425 AliInfo("OFC - ROD&Strip shift : done ");
426 } else if (look==2) {
427 PoissonRelaxation3D( arrayofArrayV, arrayofCharge,
428 fLookUpBasic3ErOverEz, fLookUpBasic3EphiOverEz, fLookUpBasic3DeltaEz,
429 kRows, kColumns, kPhiSlices, gridSizePhi, kIterations, symmetry[2]) ;
430 AliInfo("IFC - Clip rot. : done ");
431 } else if (look==3) {
432 PoissonRelaxation3D( arrayofArrayV, arrayofCharge,
433 fLookUpBasic4ErOverEz, fLookUpBasic4EphiOverEz, fLookUpBasic4DeltaEz,
434 kRows, kColumns, kPhiSlices, gridSizePhi, kIterations, symmetry[3]) ;
435 AliInfo("OFC - Clip rot. : done ");
436 } else if (look==4) {
437 PoissonRelaxation3D( arrayofArrayV, arrayofCharge,
438 fLookUpBasic5ErOverEz, fLookUpBasic5EphiOverEz, fLookUpBasic5DeltaEz,
439 kRows, kColumns, kPhiSlices, gridSizePhi, kIterations, symmetry[4]) ;
440 AliInfo("OFC - CopperRod shift : done ");
441 }
442
443 fInitLookUpBasic[look] = kTRUE;
444 }
445 }
446
447
448 // =============================================================================
449 // Create the final lookup table with corresponding ROD shifts or clip rotations
450
451 Float_t erIntegralSum = 0.0 ;
452 Float_t ephiIntegralSum = 0.0 ;
453 Float_t ezIntegralSum = 0.0 ;
454
455 Double_t phiPrime = 0. ;
456 Double_t erIntegral = 0. ;
457 Double_t ephiIntegral = 0. ;
458 Double_t ezIntegral = 0. ;
459
460
461 AliInfo("Calculation of combined Look-up Table");
462
463 // Interpolate and sum the Basic lookup tables onto the standard grid
464 Double_t r, phi, z ;
465 for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
466 phi = fgkPhiList[k] ;
467
468 TMatrixD &erOverEz = *fLookUpErOverEz[k] ;
469 TMatrixD &ephiOverEz = *fLookUpEphiOverEz[k];
470 TMatrixD &deltaEz = *fLookUpDeltaEz[k] ;
471
472 for ( Int_t i = 0 ; i < kNZ ; i++ ) {
473 z = TMath::Abs(fgkZList[i]) ; // Symmetric solution in Z that depends only on ABS(Z)
474 for ( Int_t j = 0 ; j < kNR ; j++ ) {
475 r = fgkRList[j] ;
476 // Interpolate basicLookup tables; once for each rod, then sum the results
477
478 erIntegralSum = 0.0 ;
479 ephiIntegralSum = 0.0 ;
480 ezIntegralSum = 0.0 ;
481
482 // SHIFTED RODS =========================================================
483
484 // IFC ROD SHIFTS +++++++++++++++++++++++++++++
485 for ( Int_t rod = 0 ; rod < 18 ; rod++ ) {
486
487 erIntegral = ephiIntegral = ezIntegral = 0.0 ;
488
489 if ( fRodVoltShiftA[rod] == 0 && fgkZList[i] > 0) continue ;
490 if ( fRodVoltShiftC[rod] == 0 && fgkZList[i] < 0) continue ;
491
492 // Rotate to a position where we have maps and prepare to sum
493 phiPrime = phi - rod*kPhiSlicesPerSector*gridSizePhi ;
494
495 if ( phiPrime < 0 ) phiPrime += TMath::TwoPi() ; // Stay in range from 0 to TwoPi
496
497 if ( (phiPrime >= 0) && (phiPrime <= kPhiSlices*gridSizePhi) ) {
498
499 erIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
500 rlist, zedlist, philist, fLookUpBasic1ErOverEz );
501 ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
502 rlist, zedlist, philist, fLookUpBasic1EphiOverEz);
503 ezIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
504 rlist, zedlist, philist, fLookUpBasic1DeltaEz );
505
506 } else if ( (phiPrime <= TMath::TwoPi()) && (phiPrime >= (TMath::TwoPi()-kPhiSlices*gridSizePhi)) ){
507
508 phiPrime = TMath::TwoPi() - phiPrime ;
509
510 erIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
511 rlist, zedlist, philist, fLookUpBasic1ErOverEz );
512 ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
513 rlist, zedlist, philist, fLookUpBasic1EphiOverEz);
514 ezIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
515 rlist, zedlist, philist, fLookUpBasic1DeltaEz );
516
517 // Flip symmetry of phi integral for symmetric boundary conditions
518 if ( symmetry[0] == 1 ) ephiIntegral *= -1 ;
519 // Flip symmetry of r integral if anti-symmetric boundary conditions
520 if ( symmetry[0] == -1 ) erIntegral *= -1 ;
521
522 }
523
524 if ( fgkZList[i] > 0 ) {
525 erIntegralSum += fRodVoltShiftA[rod]*erIntegral ;
526 ephiIntegralSum += fRodVoltShiftA[rod]*ephiIntegral ;
527 ezIntegralSum += fRodVoltShiftA[rod]*ezIntegral;
528 } else if ( fgkZList[i] < 0 ) {
529 erIntegralSum += fRodVoltShiftC[rod]*erIntegral ;
530 ephiIntegralSum += fRodVoltShiftC[rod]*ephiIntegral ;
531 ezIntegralSum -= fRodVoltShiftC[rod]*ezIntegral;
532 }
533 }
534
535 // OFC ROD SHIFTS +++++++++++++++++++++++++++++
536 for ( Int_t rod = 18 ; rod < 36 ; rod++ ) {
537
538 erIntegral = ephiIntegral = ezIntegral = 0.0 ;
539
540 if ( fRodVoltShiftA[rod] == 0 && fgkZList[i] > 0) continue ;
541 if ( fRodVoltShiftC[rod] == 0 && fgkZList[i] < 0) continue ;
542
543 // Rotate to a position where we have maps and prepare to sum
544 phiPrime = phi - (rod-18)*kPhiSlicesPerSector*gridSizePhi ;
545
546
547 if ( phiPrime < 0 ) phiPrime += TMath::TwoPi() ; // Stay in range from 0 to TwoPi
548
549 if ( (phiPrime >= 0) && (phiPrime <= kPhiSlices*gridSizePhi) ) {
550
551 erIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
552 rlist, zedlist, philist, fLookUpBasic2ErOverEz );
553 ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
554 rlist, zedlist, philist, fLookUpBasic2EphiOverEz);
555 ezIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
556 rlist, zedlist, philist, fLookUpBasic2DeltaEz );
557
558 } else if ( (phiPrime <= TMath::TwoPi()) && (phiPrime >= (TMath::TwoPi()-kPhiSlices*gridSizePhi)) ){
559
560 phiPrime = TMath::TwoPi() - phiPrime ;
561
562 erIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
563 rlist, zedlist, philist, fLookUpBasic2ErOverEz );
564 ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
565 rlist, zedlist, philist, fLookUpBasic2EphiOverEz);
566 ezIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
567 rlist, zedlist, philist, fLookUpBasic2DeltaEz );
568
569 // Flip symmetry of phi integral for symmetric boundary conditions
570 if ( symmetry[1] == 1 ) ephiIntegral *= -1 ;
571 // Flip symmetry of r integral if anti-symmetric boundary conditions
572 if ( symmetry[1] == -1 ) erIntegral *= -1 ;
573
574 }
575
576 if ( fgkZList[i] > 0 ) {
577 erIntegralSum += fRodVoltShiftA[rod]*erIntegral ;
578 ephiIntegralSum += fRodVoltShiftA[rod]*ephiIntegral ;
579 ezIntegralSum += fRodVoltShiftA[rod]*ezIntegral;
580 } else if ( fgkZList[i] < 0 ) {
581 erIntegralSum += fRodVoltShiftC[rod]*erIntegral ;
582 ephiIntegralSum += fRodVoltShiftC[rod]*ephiIntegral ;
583 ezIntegralSum -= fRodVoltShiftC[rod]*ezIntegral;
584 }
585
586 } // rod loop - shited ROD
587
588
589 // Rotated clips =========================================================
590
591 Int_t rodIFC = 11; // resistor rod positions, IFC
592 Int_t rodOFC = 3; // resistor rod positions, OFC
593 // just one rod on IFC and OFC
594
595 // IFC ROTATED CLIP +++++++++++++++++++++++++++++
596 for ( Int_t rod = rodIFC ; rod < rodIFC+1 ; rod++ ) { // loop over 1 to keep the "ignore"-functionality
597
598 erIntegral = ephiIntegral = ezIntegral = 0.0 ;
599 if ( fRotatedClipVoltA[0] == 0 && fgkZList[i] > 0) continue ;
600 if ( fRotatedClipVoltC[0] == 0 && fgkZList[i] < 0) continue ;
601
602 // Rotate to a position where we have maps and prepare to sum
603 phiPrime = phi - rod*kPhiSlicesPerSector*gridSizePhi ;
604
605 if ( phiPrime < 0 ) phiPrime += TMath::TwoPi() ; // Stay in range from 0 to TwoPi
606
607 if ( (phiPrime >= 0) && (phiPrime <= kPhiSlices*gridSizePhi) ) {
608
609 erIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
610 rlist, zedlist, philist, fLookUpBasic3ErOverEz );
611 ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
612 rlist, zedlist, philist, fLookUpBasic3EphiOverEz);
613 ezIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
614 rlist, zedlist, philist, fLookUpBasic3DeltaEz );
615
616 } else if ( (phiPrime <= TMath::TwoPi()) && (phiPrime >= (TMath::TwoPi()-kPhiSlices*gridSizePhi)) ){
617
618 phiPrime = TMath::TwoPi() - phiPrime ;
619
620 erIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
621 rlist, zedlist, philist, fLookUpBasic3ErOverEz );
622 ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
623 rlist, zedlist, philist, fLookUpBasic3EphiOverEz);
624 ezIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
625 rlist, zedlist, philist, fLookUpBasic3DeltaEz );
626
627 // Flip symmetry of phi integral for symmetric boundary conditions
628 if ( symmetry[2] == 1 ) ephiIntegral *= -1 ;
629 // Flip symmetry of r integral if anti-symmetric boundary conditions
630 if ( symmetry[2] == -1 ) erIntegral *= -1 ;
631
632 }
633
634 if ( fgkZList[i] > 0 ) {
635 erIntegralSum += fRotatedClipVoltA[0]*erIntegral ;
636 ephiIntegralSum += fRotatedClipVoltA[0]*ephiIntegral ;
637 ezIntegralSum += fRotatedClipVoltA[0]*ezIntegral;
638 } else if ( fgkZList[i] < 0 ) {
639 erIntegralSum += fRotatedClipVoltC[0]*erIntegral ;
640 ephiIntegralSum += fRotatedClipVoltC[0]*ephiIntegral ;
641 ezIntegralSum -= fRotatedClipVoltC[0]*ezIntegral;
642 }
643 }
644
645 // OFC: ROTATED CLIP +++++++++++++++++++++++++++++
646 for ( Int_t rod = rodOFC ; rod < rodOFC+1 ; rod++ ) { // loop over 1 to keep the "ignore"-functionality
647
648 erIntegral = ephiIntegral = ezIntegral = 0.0 ;
649
650 if ( fRotatedClipVoltA[1] == 0 && fgkZList[i] > 0) continue ;
651 if ( fRotatedClipVoltC[1] == 0 && fgkZList[i] < 0) continue ;
652
653 // Rotate to a position where we have maps and prepare to sum
654 phiPrime = phi - rod*kPhiSlicesPerSector*gridSizePhi ;
655
656
657 if ( phiPrime < 0 ) phiPrime += TMath::TwoPi() ; // Stay in range from 0 to TwoPi
658
659 if ( (phiPrime >= 0) && (phiPrime <= kPhiSlices*gridSizePhi) ) {
660
661 erIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
662 rlist, zedlist, philist, fLookUpBasic4ErOverEz );
663 ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
664 rlist, zedlist, philist, fLookUpBasic4EphiOverEz);
665 ezIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
666 rlist, zedlist, philist, fLookUpBasic4DeltaEz );
667
668 } else if ( (phiPrime <= TMath::TwoPi()) && (phiPrime >= (TMath::TwoPi()-kPhiSlices*gridSizePhi)) ){
669
670 phiPrime = TMath::TwoPi() - phiPrime ;
671
672 erIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
673 rlist, zedlist, philist, fLookUpBasic4ErOverEz );
674 ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
675 rlist, zedlist, philist, fLookUpBasic4EphiOverEz);
676 ezIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
677 rlist, zedlist, philist, fLookUpBasic4DeltaEz );
678
679 // Flip symmetry of phi integral for symmetric boundary conditions
680 if ( symmetry[3] == 1 ) ephiIntegral *= -1 ;
681 // Flip symmetry of r integral if anti-symmetric boundary conditions
682 if ( symmetry[3] == -1 ) erIntegral *= -1 ;
683
684 }
685
686 if ( fgkZList[i] > 0 ) {
687 erIntegralSum += fRotatedClipVoltA[1]*erIntegral ;
688 ephiIntegralSum += fRotatedClipVoltA[1]*ephiIntegral ;
689 ezIntegralSum += fRotatedClipVoltA[1]*ezIntegral;
690 } else if ( fgkZList[i] < 0 ) {
691 erIntegralSum += fRotatedClipVoltC[1]*erIntegral ;
692 ephiIntegralSum += fRotatedClipVoltC[1]*ephiIntegral ;
693 ezIntegralSum -= fRotatedClipVoltC[1]*ezIntegral;
694 }
695 }
696
697 // OFC Cooper rod shift +++++++++++++++++++++++++++++
698 for ( Int_t rod = 0 ; rod < 18 ; rod++ ) {
699
700 erIntegral = ephiIntegral = ezIntegral = 0.0 ;
701
702 if ( fOFCRodShiftA[rod] == 0 && fgkZList[i] > 0) continue ;
703 if ( fOFCRodShiftC[rod] == 0 && fgkZList[i] < 0) continue ;
704
705 // Rotate to a position where we have maps and prepare to sum
706 phiPrime = phi - rod*kPhiSlicesPerSector*gridSizePhi ;
707
708 if ( phiPrime < 0 ) phiPrime += TMath::TwoPi() ; // Stay in range from 0 to TwoPi
709
710 if ( (phiPrime >= 0) && (phiPrime <= kPhiSlices*gridSizePhi) ) {
711
712 erIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
713 rlist, zedlist, philist, fLookUpBasic5ErOverEz );
714 ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
715 rlist, zedlist, philist, fLookUpBasic5EphiOverEz);
716 ezIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
717 rlist, zedlist, philist, fLookUpBasic5DeltaEz );
718
719 } else if ( (phiPrime <= TMath::TwoPi()) && (phiPrime >= (TMath::TwoPi()-kPhiSlices*gridSizePhi)) ){
720
721 phiPrime = TMath::TwoPi() - phiPrime ;
722
723 erIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
724 rlist, zedlist, philist, fLookUpBasic5ErOverEz );
725 ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
726 rlist, zedlist, philist, fLookUpBasic5EphiOverEz);
727 ezIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
728 rlist, zedlist, philist, fLookUpBasic5DeltaEz );
729
730 // Flip symmetry of phi integral for symmetric boundary conditions
731 if ( symmetry[0] == 1 ) ephiIntegral *= -1 ;
732 // Flip symmetry of r integral if anti-symmetric boundary conditions
733 if ( symmetry[0] == -1 ) erIntegral *= -1 ;
734
735 }
736
737 if ( fgkZList[i] > 0 ) {
738 erIntegralSum += fOFCRodShiftA[rod]*erIntegral ;
739 ephiIntegralSum += fOFCRodShiftA[rod]*ephiIntegral ;
740 ezIntegralSum += fOFCRodShiftA[rod]*ezIntegral;
741 } else if ( fgkZList[i] < 0 ) {
742 erIntegralSum += fOFCRodShiftC[rod]*erIntegral ;
743 ephiIntegralSum += fOFCRodShiftC[rod]*ephiIntegral ;
744 ezIntegralSum -= fOFCRodShiftC[rod]*ezIntegral;
745 }
746 }
747
748 // put the sum into the final lookup table
749 erOverEz(j,i) = erIntegralSum;
750 ephiOverEz(j,i) = ephiIntegralSum;
751 deltaEz(j,i) = ezIntegralSum;
752
753 // if (j==1) printf("%lf %lf %lf: %lf %lf %lf\n",r, phi, z, erIntegralSum,ephiIntegralSum,ezIntegralSum);
754
755 } // endo r loop
756 } // end of z loop
757 } // end of phi loop
758
759
760 // clear the temporary arrays lists
761 for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
762 delete arrayofArrayV[k];
763 delete arrayofCharge[k];
764 }
765
766 AliInfo(" done");
767 fInitLookUp = kTRUE;
768
769}
770
771void AliTPCFCVoltError3D::Print(const Option_t* option) const {
772 //
773 // Print function to check the settings of the Rod shifts and the rotated clips
774 // option=="a" prints the C0 and C1 coefficents for calibration purposes
775 //
776
777 TString opt = option; opt.ToLower();
778 printf("%s\n",GetTitle());
779 printf(" - ROD shifts (residual voltage settings): 40V correspond to 1mm of shift\n");
780 Int_t count = 0;
781 printf(" : A-side - IFC (Rod & Strips) \n ");
782 for (Int_t i=0; i<18;i++) {
783 if (fRodVoltShiftA[i]!=0) {
784 printf("Rod %2d: %3.1f V \t",i,fRodVoltShiftA[i]);
785 count++;
786 }
787 if (count==0&&i==17)
788 printf("-> all at 0 V\n");
789 else if (i==17){
790 printf("\n");
791 count=0;
792 }
793 }
794 printf(" : C-side - IFC (Rod & Strips) \n ");
795 for (Int_t i=0; i<18;i++) {
796 if (fRodVoltShiftC[i]!=0) {
797 printf("Rod %2d: %3.1f V \t",i,fRodVoltShiftC[i]);
798 count++;
799 }
800 if (count==0&&i==17)
801 printf("-> all at 0 V\n");
802 else if (i==17){
803 printf("\n");
804 count=0;
805 }
806 }
807 printf(" : A-side - OFC (only Rod) \n ");
808 for (Int_t i=18; i<36;i++) {
809 if (fRodVoltShiftA[i]!=0) {
810 printf("Rod %2d: %3.1f V \t",i,fRodVoltShiftA[i]);
811 count++;
812 }
813 if (count==0&&i==35)
814 printf("-> all at 0 V\n");
815 else if (i==35) {
816 printf("\n");
817 count=0;
818 }
819 }
820 printf(" : C-side - OFC (only Rod) \n ");
821 for (Int_t i=18; i<36;i++) {
822 if (fRodVoltShiftC[i]!=0) {
823 printf("Rod %2d: %3.1f V \t",i,fRodVoltShiftC[i]);
824 count++;
825 }
826 if (count==0&&i==35)
827 printf("-> all at 0 V\n");
828 else if (i==35){
829 printf("\n");
830 count=0;
831 }
832 }
833
834 printf(" - Rotated clips (residual voltage settings): 40V correspond to 1mm of 'offset'\n");
835 if (fRotatedClipVoltA[0]!=0) { printf(" A side (IFC): HV Rod: %3.1f V \n",fRotatedClipVoltA[0]); count++; }
836 if (fRotatedClipVoltA[1]!=0) { printf(" A side (OFC): HV Rod: %3.1f V \n",fRotatedClipVoltA[1]); count++; }
837 if (fRotatedClipVoltC[0]!=0) { printf(" C side (IFC): HV Rod: %3.1f V \n",fRotatedClipVoltC[0]); count++; }
838 if (fRotatedClipVoltC[1]!=0) { printf(" C side (OFC): HV Rod: %3.1f V \n",fRotatedClipVoltC[1]); count++; }
839 if (count==0)
840 printf(" -> no rotated clips \n");
841
842 count=0;
843 printf(" - Copper ROD shifts (without strips):\n");
844 printf(" : A-side - OFC (Copper Rod shift) \n ");
845 for (Int_t i=0; i<18;i++) {
846 if (fOFCRodShiftA[i]!=0) {
847 printf("Rod %2d: %3.1f V \t",i,fOFCRodShiftA[i]);
848 count++;
849 }
850 if (count==0&&i==17)
851 printf("-> all at 0 V\n");
852 else if (i==17){
853 printf("\n");
854 count=0;
855 }
856 }
857 printf(" : C-side - OFC (Copper Rod shift) \n ");
858 for (Int_t i=0; i<18;i++) {
859 if (fOFCRodShiftC[i]!=0) {
860 printf("Rod %2d: %3.1f V \t",i,fOFCRodShiftC[i]);
861 count++;
862 }
863 if (count==0&&i==17)
864 printf("-> all at 0 V\n");
865 else if (i==17){
866 printf("\n");
867 count=0;
868 }
869 }
870
871 if (opt.Contains("a")) { // Print all details
872 printf(" - T1: %1.4f, T2: %1.4f \n",fT1,fT2);
873 printf(" - C1: %1.4f, C0: %1.4f \n",fC1,fC0);
874 }
875
876 if (!fInitLookUp) AliError("Lookup table was not initialized! You should do InitFCVoltError3D() ...");
877
878}