1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 // _________________________________________________________________
19 // <h2>AliTPCFCVoltError3D class </h2>
20 // The class calculates the space point distortions due to residual voltage errors
21 // on the Field Cage (FC) boundaries (rods and strips) of the TPC in 3D. It uses
22 // the Poisson relaxation technique in three dimension as implemented in the parent class.
24 // Although the calculation is performed in 3D, the calculation time was significantly
25 // reduced by using certain symmetry conditions within the calculation.
27 // The input parameters can be set via the functions (e.g.) SetRodVoltShift(rod,dV) where
28 // rod is the number of the rod on which the voltage offset dV is set. The corresponding
29 // shift in z direction would be $dz=dV/40$ with an opposite sign for the C side. The
30 // rods are numbered in anti-clock-wise direction starting at $\phi=0$. Rods in the IFC
31 // are from 0 to 17, rods on the OFC go from 18 to 35. <br>
32 // This convention is following the offline numbering scheme of the ROCs.
34 // Different misalignment scenarios can be modeled:
36 // <li> A rotated clip scenario is only possible at two positions (Rod 11 at IFC, rod 3(+18) at OFC)
37 // and can be set via SetRotatedClipVolt. The key feature is that at the mentioned positions,
38 // the strip-ends were combined. At these positions, a systematic offset of one strip-end in
39 // respect to the other is possible.
40 // <li> A normal rod offset, where the strips follow the movement of the rod, can be set via
41 // SetRodVoltShift. It has a anti-mirrored signature in terms of distortions when compared
42 // to the rotated clip. This misalignment is possible at each single rod of the FC.
43 // <li> A simple rod offset, where the strips do not follow the shift, results in an even more
44 // localized distortion close to the rod. The difference to a rod shift, where the strips follow,
45 // is more dominant on the OFC. This effect can be set via the function SetCopperRodShift.
49 // Begin_Macro(source)
51 // gROOT->SetStyle("Plain"); gStyle->SetPalette(1);
52 // TCanvas *c2 = new TCanvas("c2","c2",500,450);
53 // AliTPCFCVoltError3D fc;
54 // fc.SetOmegaTauT1T2(0,1,1);
55 // fc.SetRotatedClipVoltA(0,40);
56 // fc.SetRodVoltShiftA(3,40);
57 // fc.SetCopperRodShiftA(7+18,40);
58 // fc.SetRodVoltShiftA(15+18,40);
59 // fc.CreateHistoDRPhiinXY(10)->Draw("cont4z");
66 // Date: 08/08/2010 <br>
67 // Authors: Jim Thomas, Stefan Rossegger
69 // _________________________________________________________________
73 #include "TGeoGlobalMagField.h"
74 #include "AliTPCcalibDB.h"
75 #include "AliTPCParam.h"
80 #include "AliTPCROC.h"
81 #include "AliTPCFCVoltError3D.h"
83 ClassImp(AliTPCFCVoltError3D)
85 AliTPCFCVoltError3D::AliTPCFCVoltError3D()
86 : AliTPCCorrection("FieldCageVoltErrors","FieldCage (Rods) Voltage Errors"),
91 // default constructor
94 // flags for filled 'basic lookup tables'
95 for (Int_t i=0; i<6; i++){
96 fInitLookUpBasic[i]= kFALSE;
100 for (Int_t i=0; i<36; i++){
101 fRodVoltShiftA[i] = 0;
102 fRodVoltShiftC[i] = 0;
104 for (Int_t i=0; i<2; i++){
105 fRotatedClipVoltA[i] = 0;
106 fRotatedClipVoltC[i] = 0;
109 for (Int_t i=0; i<36; i++){
110 fCopperRodShiftA[i] = 0;
111 fCopperRodShiftC[i] = 0;
114 // Array which will contain the solution according to the setted boundary conditions
115 // it represents a sum up of the 4 basic look up tables (created individually)
116 // see InitFCVoltError3D() function
117 for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
118 fLookUpErOverEz[k] = new TMatrixD(kNR,kNZ);
119 fLookUpEphiOverEz[k] = new TMatrixD(kNR,kNZ);
120 fLookUpDeltaEz[k] = new TMatrixD(kNR,kNZ);
123 for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
124 fLookUpBasic1ErOverEz[k] = 0;
125 fLookUpBasic1EphiOverEz[k] = 0;
126 fLookUpBasic1DeltaEz[k] = 0;
128 fLookUpBasic2ErOverEz[k] = 0;
129 fLookUpBasic2EphiOverEz[k] = 0;
130 fLookUpBasic2DeltaEz[k] = 0;
132 fLookUpBasic3ErOverEz[k] = 0;
133 fLookUpBasic3EphiOverEz[k] = 0;
134 fLookUpBasic3DeltaEz[k] = 0;
136 fLookUpBasic4ErOverEz[k] = 0;
137 fLookUpBasic4EphiOverEz[k] = 0;
138 fLookUpBasic4DeltaEz[k] = 0;
140 fLookUpBasic5ErOverEz[k] = 0;
141 fLookUpBasic5EphiOverEz[k] = 0;
142 fLookUpBasic5DeltaEz[k] = 0;
144 fLookUpBasic6ErOverEz[k] = 0;
145 fLookUpBasic6EphiOverEz[k] = 0;
146 fLookUpBasic6DeltaEz[k] = 0;
151 AliTPCFCVoltError3D::~AliTPCFCVoltError3D() {
156 for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
157 delete fLookUpErOverEz[k];
158 delete fLookUpEphiOverEz[k];
159 delete fLookUpDeltaEz[k];
162 for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
163 delete fLookUpBasic1ErOverEz[k]; // does nothing if pointer is zero!
164 delete fLookUpBasic1EphiOverEz[k];
165 delete fLookUpBasic1DeltaEz[k];
167 delete fLookUpBasic2ErOverEz[k]; // does nothing if pointer is zero!
168 delete fLookUpBasic2EphiOverEz[k];
169 delete fLookUpBasic2DeltaEz[k];
171 delete fLookUpBasic3ErOverEz[k]; // does nothing if pointer is zero!
172 delete fLookUpBasic3EphiOverEz[k];
173 delete fLookUpBasic3DeltaEz[k];
175 delete fLookUpBasic4ErOverEz[k]; // does nothing if pointer is zero!
176 delete fLookUpBasic4EphiOverEz[k];
177 delete fLookUpBasic4DeltaEz[k];
179 delete fLookUpBasic5ErOverEz[k]; // does nothing if pointer is zero!
180 delete fLookUpBasic5EphiOverEz[k];
181 delete fLookUpBasic5DeltaEz[k];
183 delete fLookUpBasic6ErOverEz[k]; // does nothing if pointer is zero!
184 delete fLookUpBasic6EphiOverEz[k];
185 delete fLookUpBasic6DeltaEz[k];
190 void AliTPCFCVoltError3D::Init() {
192 // Initialization funtion
195 AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
196 if (!magF) AliError("Magneticd field - not initialized");
197 Double_t bzField = magF->SolenoidField()/10.; //field in T
198 AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
199 if (!param) AliError("Parameters - not initialized");
200 Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us] // From dataBase: to be updated: per second (ideally)
201 Double_t ezField = 400; // [V/cm] // to be updated: never (hopefully)
202 Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
203 // Correction Terms for effective omegaTau; obtained by a laser calibration run
204 SetOmegaTauT1T2(wt,fT1,fT2);
206 if (!fInitLookUp) InitFCVoltError3D();
209 void AliTPCFCVoltError3D::Update(const TTimeStamp &/*timeStamp*/) {
213 AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
214 if (!magF) AliError("Magneticd field - not initialized");
215 Double_t bzField = magF->SolenoidField()/10.; //field in T
216 AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
217 if (!param) AliError("Parameters - not initialized");
218 Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us] // From dataBase: to be updated: per second (ideally)
219 Double_t ezField = 400; // [V/cm] // to be updated: never (hopefully)
220 Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
221 // Correction Terms for effective omegaTau; obtained by a laser calibration run
222 SetOmegaTauT1T2(wt,fT1,fT2);
229 void AliTPCFCVoltError3D::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {
231 // Calculates the correction due e.g. residual voltage errors on the TPC boundaries
235 AliInfo("Lookup table was not initialized! Perform the inizialisation now ...");
239 Int_t order = 1 ; // FIXME: hardcoded? Linear interpolation = 1, Quadratic = 2
240 // note that the poisson solution was linearly mirroed on this grid!
241 Double_t intEr, intEphi, intDeltaEz;
245 r = TMath::Sqrt( x[0]*x[0] + x[1]*x[1] ) ;
246 phi = TMath::ATan2(x[1],x[0]) ;
247 if ( phi < 0 ) phi += TMath::TwoPi() ; // Table uses phi from 0 to 2*Pi
248 z = x[2] ; // Create temporary copy of x[2]
250 if ( (roc%36) < 18 ) {
251 sign = 1; // (TPC A side)
253 sign = -1; // (TPC C side)
256 if ( sign==1 && z < fgkZOffSet ) z = fgkZOffSet; // Protect against discontinuity at CE
257 if ( sign==-1 && z > -fgkZOffSet ) z = -fgkZOffSet; // Protect against discontinuity at CE
260 if ( (sign==1 && z<0) || (sign==-1 && z>0) ) // just a consistency check
261 AliError("ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!");
263 // Get the Er and Ephi field integrals plus the integral over DeltaEz
264 intEr = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
265 fgkRList, fgkZList, fgkPhiList, fLookUpErOverEz );
266 intEphi = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
267 fgkRList, fgkZList, fgkPhiList, fLookUpEphiOverEz );
268 intDeltaEz = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
269 fgkRList, fgkZList, fgkPhiList, fLookUpDeltaEz );
271 // printf("%lf %lf %lf\n",intEr,intEphi,intDeltaEz);
273 // Calculate distorted position
275 phi = phi + ( fC0*intEphi - fC1*intEr ) / r;
276 r = r + ( fC0*intEr + fC1*intEphi );
279 // Calculate correction in cartesian coordinates
280 dx[0] = r * TMath::Cos(phi) - x[0];
281 dx[1] = r * TMath::Sin(phi) - x[1];
282 dx[2] = intDeltaEz; // z distortion - (internally scaled with driftvelocity dependency
283 // on the Ez field plus the actual ROC misalignment (if set TRUE)
287 void AliTPCFCVoltError3D::InitFCVoltError3D() {
289 // Initialization of the Lookup table which contains the solutions of the
290 // Dirichlet boundary problem
291 // Calculation of the single 3D-Poisson solver is done just if needed
292 // (see basic lookup tables in header file)
295 const Int_t order = 1 ; // Linear interpolation = 1, Quadratic = 2
296 const Float_t gridSizeR = (fgkOFCRadius-fgkIFCRadius) / (kRows-1) ;
297 const Float_t gridSizeZ = fgkTPCZ0 / (kColumns-1) ;
298 const Float_t gridSizePhi = TMath::TwoPi() / ( 18.0 * kPhiSlicesPerSector);
300 // temporary arrays to create the boundary conditions
301 TMatrixD *arrayofArrayV[kPhiSlices], *arrayofCharge[kPhiSlices] ;
302 for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
303 arrayofArrayV[k] = new TMatrixD(kRows,kColumns) ;
304 arrayofCharge[k] = new TMatrixD(kRows,kColumns) ;
306 Double_t innerList[kPhiSlices], outerList[kPhiSlices] ;
308 // list of point as used in the poisson relation and the interpolation (during sum up)
309 Double_t rlist[kRows], zedlist[kColumns] , philist[kPhiSlices];
310 for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
311 philist[k] = gridSizePhi * k;
312 for ( Int_t i = 0 ; i < kRows ; i++ ) {
313 rlist[i] = fgkIFCRadius + i*gridSizeR ;
314 for ( Int_t j = 0 ; j < kColumns ; j++ ) { // Fill Vmatrix with Boundary Conditions
315 zedlist[j] = j * gridSizeZ ;
320 // ==========================================================================
321 // Solve Poisson's equation in 3D cylindrical coordinates by relaxation technique
322 // Allow for different size grid spacing in R and Z directions
324 const Int_t symmetry[6] = {1,1,-1,-1,1,1}; // shifted rod (2x), rotated clip (2x), only rod shift on OFC (1x)
326 // check which lookup tables are needed ---------------------------------
328 Bool_t needTable[6] ={kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE};
330 // Shifted rods & strips
331 for ( Int_t rod = 0 ; rod < 18 ; rod++ ) {
332 if (fRodVoltShiftA[rod]!=0) needTable[0]=kTRUE;
333 if (fRodVoltShiftC[rod]!=0) needTable[0]=kTRUE;
335 for ( Int_t rod = 18 ; rod < 36 ; rod++ ) {
336 if (fRodVoltShiftA[rod]!=0) needTable[1]=kTRUE;
337 if (fRodVoltShiftC[rod]!=0) needTable[1]=kTRUE;
340 if (fRotatedClipVoltA[0]!=0 || fRotatedClipVoltC[0]!=0) needTable[2]=kTRUE;
341 if (fRotatedClipVoltA[1]!=0 || fRotatedClipVoltC[1]!=0) needTable[3]=kTRUE;
343 // shifted Copper rods
344 for ( Int_t rod = 0 ; rod < 18 ; rod++ ) {
345 if (fCopperRodShiftA[rod]!=0) needTable[4]=kTRUE;
346 if (fCopperRodShiftC[rod]!=0) needTable[4]=kTRUE;
348 // shifted Copper rods
349 for ( Int_t rod = 18; rod < 36 ; rod++ ) {
350 if (fCopperRodShiftA[rod]!=0) needTable[5]=kTRUE;
351 if (fCopperRodShiftC[rod]!=0) needTable[5]=kTRUE;
355 // reserve the arrays for the basic lookup tables ----------------------
356 if (needTable[0] && !fInitLookUpBasic[0]) {
357 for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) { // Possibly make an extra table to be used for 0 == 360
358 fLookUpBasic1ErOverEz[k] = new TMatrixD(kRows,kColumns);
359 fLookUpBasic1EphiOverEz[k] = new TMatrixD(kRows,kColumns);
360 fLookUpBasic1DeltaEz[k] = new TMatrixD(kRows,kColumns);
361 // will be deleted in destructor
364 if (needTable[1] && !fInitLookUpBasic[1]) {
365 for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) { // Possibly make an extra table to be used for 0 == 360
366 fLookUpBasic2ErOverEz[k] = new TMatrixD(kRows,kColumns);
367 fLookUpBasic2EphiOverEz[k] = new TMatrixD(kRows,kColumns);
368 fLookUpBasic2DeltaEz[k] = new TMatrixD(kRows,kColumns);
369 // will be deleted in destructor
372 if (needTable[2] && !fInitLookUpBasic[2]) {
373 for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) { // Possibly make an extra table to be used for 0 == 360
374 fLookUpBasic3ErOverEz[k] = new TMatrixD(kRows,kColumns);
375 fLookUpBasic3EphiOverEz[k] = new TMatrixD(kRows,kColumns);
376 fLookUpBasic3DeltaEz[k] = new TMatrixD(kRows,kColumns);
377 // will be deleted in destructor
380 if (needTable[3] && !fInitLookUpBasic[3]) {
381 for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) { // Possibly make an extra table to be used for 0 == 360
382 fLookUpBasic4ErOverEz[k] = new TMatrixD(kRows,kColumns);
383 fLookUpBasic4EphiOverEz[k] = new TMatrixD(kRows,kColumns);
384 fLookUpBasic4DeltaEz[k] = new TMatrixD(kRows,kColumns);
385 // will be deleted in destructor
388 if (needTable[4] && !fInitLookUpBasic[4]) {
389 for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) { // Possibly make an extra table to be used for 0 == 360
390 fLookUpBasic5ErOverEz[k] = new TMatrixD(kRows,kColumns);
391 fLookUpBasic5EphiOverEz[k] = new TMatrixD(kRows,kColumns);
392 fLookUpBasic5DeltaEz[k] = new TMatrixD(kRows,kColumns);
393 // will be deleted in destructor
396 if (needTable[5] && !fInitLookUpBasic[5]) {
397 for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) { // Possibly make an extra table to be used for 0 == 360
398 fLookUpBasic6ErOverEz[k] = new TMatrixD(kRows,kColumns);
399 fLookUpBasic6EphiOverEz[k] = new TMatrixD(kRows,kColumns);
400 fLookUpBasic6DeltaEz[k] = new TMatrixD(kRows,kColumns);
401 // will be deleted in destructor
405 // Set bondaries and solve Poisson's equation --------------------------
407 for (Int_t look=0; look<6; look++) {
409 Float_t inner18[18] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } ;
410 Float_t outer18[18] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } ;
412 if (needTable[look] && !fInitLookUpBasic[look]) {
414 // Specify which rod is the reference; other distortions calculated by summing over multiple rotations of refrence
415 // Note: the interpolation later on depends on it!! Do not change if not really needed!
417 AliInfo(Form("IFC - ROD&Strip shift : Filling table (~ %d sec)",3*(int)(kPhiSlices/5)));
419 } else if (look==1) {
420 AliInfo(Form("OFC - ROD&Strip shift : Filling table (~ %d sec)",3*(int)(kPhiSlices/5)));
422 } else if (look==2) {
423 AliInfo(Form("IFC - Clip rot. : Filling table (~ %d sec)",3*(int)(kPhiSlices/5)));
425 } else if (look==3) {
426 AliInfo(Form("OFC - Clip rot. : Filling table (~ %d sec)",3*(int)(kPhiSlices/5)));
428 } else if (look==4) {
429 AliInfo(Form("IFC - CopperRod shift : Filling table (~ %d sec)",3*(int)(kPhiSlices/5)));
431 } else if (look==5) {
432 AliInfo(Form("OFC - CopperRod shift : Filling table (~ %d sec)",3*(int)(kPhiSlices/5)));
435 // Build potentials on boundary stripes for a rotated clip (SYMMETRY==-1) or a shifted rod (SYMMETRY==1)
437 // in these cases, the strips follow the actual rod shift (linear interpolation between the rods)
438 for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
439 Int_t slices = kPhiSlicesPerSector ;
440 Int_t ipoint = k/slices ;
441 innerList[k] = (((ipoint+1)*slices-k)*inner18[ipoint]-(k-ipoint*slices)*inner18[ipoint+1])/slices ;
442 outerList[k] = (((ipoint+1)*slices-k)*outer18[ipoint]-(k-ipoint*slices)*outer18[ipoint+1])/slices ;
443 if ( (k%slices) == 0 && symmetry[look] == -1 ) { innerList[k] = 0.0 ; outerList[k] = 0.0 ; }
444 // above, force Zero if Anti-Sym
447 // in this case, we assume that the strips stay at the correct position, but the rods move
448 // the distortion is then much more localized around the rod (indicated by real data)
449 for ( Int_t k = 0 ; k < kPhiSlices ; k++ )
450 innerList[k] = outerList[k] = 0;
451 innerList[0] = inner18[0]; // point at rod
452 innerList[0] = inner18[0]/4*3; // point close to rod (educated guess)
454 // in this case, we assume that the strips stay at the correct position, but the copper plated OFC-rods move
455 // the distortion is then much more localized around the rod (indicated by real data)
456 for ( Int_t k = 0 ; k < kPhiSlices ; k++ )
457 innerList[k] = outerList[k] = 0;
458 outerList[0] = outer18[0]; // point at rod
459 outerList[1] = outer18[0]/4; // point close to rod (educated-`guessed` scaling)
462 // Fill arrays with initial conditions. V on the boundary and Charge in the volume.
463 for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
464 TMatrixD &arrayV = *arrayofArrayV[k] ;
465 TMatrixD &charge = *arrayofCharge[k] ;
466 for ( Int_t i = 0 ; i < kRows ; i++ ) {
467 for ( Int_t j = 0 ; j < kColumns ; j++ ) { // Fill Vmatrix with Boundary Conditions
470 if ( i == 0 ) arrayV(i,j) = innerList[k] ;
471 if ( i == kRows-1 ) arrayV(i,j) = outerList[k] ;
474 // no charge in the volume
475 for ( Int_t i = 1 ; i < kRows-1 ; i++ ) {
476 for ( Int_t j = 1 ; j < kColumns-1 ; j++ ) {
482 // Solve Poisson's equation in 3D cylindrical coordinates by relaxation technique
483 // Allow for different size grid spacing in R and Z directions
485 PoissonRelaxation3D( arrayofArrayV, arrayofCharge,
486 fLookUpBasic1ErOverEz, fLookUpBasic1EphiOverEz, fLookUpBasic1DeltaEz,
487 kRows, kColumns, kPhiSlices, gridSizePhi, kIterations, symmetry[0]) ;
488 AliInfo("IFC - ROD&Strip shift : done ");
489 } else if (look==1) {
490 PoissonRelaxation3D( arrayofArrayV, arrayofCharge,
491 fLookUpBasic2ErOverEz, fLookUpBasic2EphiOverEz, fLookUpBasic2DeltaEz,
492 kRows, kColumns, kPhiSlices, gridSizePhi, kIterations, symmetry[1]) ;
494 AliInfo("OFC - ROD&Strip shift : done ");
495 } else if (look==2) {
496 PoissonRelaxation3D( arrayofArrayV, arrayofCharge,
497 fLookUpBasic3ErOverEz, fLookUpBasic3EphiOverEz, fLookUpBasic3DeltaEz,
498 kRows, kColumns, kPhiSlices, gridSizePhi, kIterations, symmetry[2]) ;
499 AliInfo("IFC - Clip rot. : done ");
500 } else if (look==3) {
501 PoissonRelaxation3D( arrayofArrayV, arrayofCharge,
502 fLookUpBasic4ErOverEz, fLookUpBasic4EphiOverEz, fLookUpBasic4DeltaEz,
503 kRows, kColumns, kPhiSlices, gridSizePhi, kIterations, symmetry[3]) ;
504 AliInfo("OFC - Clip rot. : done ");
505 } else if (look==4) {
506 PoissonRelaxation3D( arrayofArrayV, arrayofCharge,
507 fLookUpBasic5ErOverEz, fLookUpBasic5EphiOverEz, fLookUpBasic5DeltaEz,
508 kRows, kColumns, kPhiSlices, gridSizePhi, kIterations, symmetry[4]) ;
509 AliInfo("IFC - CopperRod shift : done ");
510 } else if (look==5) {
511 PoissonRelaxation3D( arrayofArrayV, arrayofCharge,
512 fLookUpBasic6ErOverEz, fLookUpBasic6EphiOverEz, fLookUpBasic6DeltaEz,
513 kRows, kColumns, kPhiSlices, gridSizePhi, kIterations, symmetry[5]) ;
514 AliInfo("OFC - CopperRod shift : done ");
517 fInitLookUpBasic[look] = kTRUE;
522 // =============================================================================
523 // Create the final lookup table with corresponding ROD shifts or clip rotations
525 Float_t erIntegralSum = 0.0 ;
526 Float_t ephiIntegralSum = 0.0 ;
527 Float_t ezIntegralSum = 0.0 ;
529 Double_t phiPrime = 0. ;
530 Double_t erIntegral = 0. ;
531 Double_t ephiIntegral = 0. ;
532 Double_t ezIntegral = 0. ;
535 AliInfo("Calculation of combined Look-up Table");
537 // Interpolate and sum the Basic lookup tables onto the standard grid
539 for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
540 phi = fgkPhiList[k] ;
542 TMatrixD &erOverEz = *fLookUpErOverEz[k] ;
543 TMatrixD &ephiOverEz = *fLookUpEphiOverEz[k];
544 TMatrixD &deltaEz = *fLookUpDeltaEz[k] ;
546 for ( Int_t i = 0 ; i < kNZ ; i++ ) {
547 z = TMath::Abs(fgkZList[i]) ; // Symmetric solution in Z that depends only on ABS(Z)
548 for ( Int_t j = 0 ; j < kNR ; j++ ) {
550 // Interpolate basicLookup tables; once for each rod, then sum the results
552 erIntegralSum = 0.0 ;
553 ephiIntegralSum = 0.0 ;
554 ezIntegralSum = 0.0 ;
556 // SHIFTED RODS =========================================================
558 // IFC ROD SHIFTS +++++++++++++++++++++++++++++
559 for ( Int_t rod = 0 ; rod < 18 ; rod++ ) {
561 erIntegral = ephiIntegral = ezIntegral = 0.0 ;
563 if ( fRodVoltShiftA[rod] == 0 && fgkZList[i] > 0) continue ;
564 if ( fRodVoltShiftC[rod] == 0 && fgkZList[i] < 0) continue ;
566 // Rotate to a position where we have maps and prepare to sum
567 phiPrime = phi - rod*kPhiSlicesPerSector*gridSizePhi ;
569 if ( phiPrime < 0 ) phiPrime += TMath::TwoPi() ; // Stay in range from 0 to TwoPi
571 if ( (phiPrime >= 0) && (phiPrime <= kPhiSlices*gridSizePhi) ) {
573 erIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
574 rlist, zedlist, philist, fLookUpBasic1ErOverEz );
575 ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
576 rlist, zedlist, philist, fLookUpBasic1EphiOverEz);
577 ezIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
578 rlist, zedlist, philist, fLookUpBasic1DeltaEz );
580 } else if ( (phiPrime <= TMath::TwoPi()) && (phiPrime >= (TMath::TwoPi()-kPhiSlices*gridSizePhi)) ){
582 phiPrime = TMath::TwoPi() - phiPrime ;
584 erIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
585 rlist, zedlist, philist, fLookUpBasic1ErOverEz );
586 ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
587 rlist, zedlist, philist, fLookUpBasic1EphiOverEz);
588 ezIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
589 rlist, zedlist, philist, fLookUpBasic1DeltaEz );
591 // Flip symmetry of phi integral for symmetric boundary conditions
592 if ( symmetry[0] == 1 ) ephiIntegral *= -1 ;
593 // Flip symmetry of r integral if anti-symmetric boundary conditions
594 if ( symmetry[0] == -1 ) erIntegral *= -1 ;
598 if ( fgkZList[i] > 0 ) {
599 erIntegralSum += fRodVoltShiftA[rod]*erIntegral ;
600 ephiIntegralSum += fRodVoltShiftA[rod]*ephiIntegral ;
601 ezIntegralSum += fRodVoltShiftA[rod]*ezIntegral;
602 } else if ( fgkZList[i] < 0 ) {
603 erIntegralSum += fRodVoltShiftC[rod]*erIntegral ;
604 ephiIntegralSum += fRodVoltShiftC[rod]*ephiIntegral ;
605 ezIntegralSum -= fRodVoltShiftC[rod]*ezIntegral;
609 // OFC ROD SHIFTS +++++++++++++++++++++++++++++
610 for ( Int_t rod = 18 ; rod < 36 ; rod++ ) {
612 erIntegral = ephiIntegral = ezIntegral = 0.0 ;
614 if ( fRodVoltShiftA[rod] == 0 && fgkZList[i] > 0) continue ;
615 if ( fRodVoltShiftC[rod] == 0 && fgkZList[i] < 0) continue ;
617 // Rotate to a position where we have maps and prepare to sum
618 phiPrime = phi - (rod-18)*kPhiSlicesPerSector*gridSizePhi ;
620 if ( phiPrime < 0 ) phiPrime += TMath::TwoPi() ; // Stay in range from 0 to TwoPi
622 if ( (phiPrime >= 0) && (phiPrime <= kPhiSlices*gridSizePhi) ) {
624 erIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
625 rlist, zedlist, philist, fLookUpBasic2ErOverEz );
626 ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
627 rlist, zedlist, philist, fLookUpBasic2EphiOverEz);
628 ezIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
629 rlist, zedlist, philist, fLookUpBasic2DeltaEz );
631 } else if ( (phiPrime <= TMath::TwoPi()) && (phiPrime >= (TMath::TwoPi()-kPhiSlices*gridSizePhi)) ){
633 phiPrime = TMath::TwoPi() - phiPrime ;
635 erIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
636 rlist, zedlist, philist, fLookUpBasic2ErOverEz );
637 ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
638 rlist, zedlist, philist, fLookUpBasic2EphiOverEz);
639 ezIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
640 rlist, zedlist, philist, fLookUpBasic2DeltaEz );
642 // Flip symmetry of phi integral for symmetric boundary conditions
643 if ( symmetry[1] == 1 ) ephiIntegral *= -1 ;
644 // Flip symmetry of r integral if anti-symmetric boundary conditions
645 if ( symmetry[1] == -1 ) erIntegral *= -1 ;
649 if ( fgkZList[i] > 0 ) {
650 erIntegralSum += fRodVoltShiftA[rod]*erIntegral ;
651 ephiIntegralSum += fRodVoltShiftA[rod]*ephiIntegral ;
652 ezIntegralSum += fRodVoltShiftA[rod]*ezIntegral;
653 } else if ( fgkZList[i] < 0 ) {
654 erIntegralSum += fRodVoltShiftC[rod]*erIntegral ;
655 ephiIntegralSum += fRodVoltShiftC[rod]*ephiIntegral ;
656 ezIntegralSum -= fRodVoltShiftC[rod]*ezIntegral;
659 } // rod loop - shited ROD
662 // Rotated clips =========================================================
664 Int_t rodIFC = 11; // resistor rod positions, IFC
665 Int_t rodOFC = 3; // resistor rod positions, OFC
666 // just one rod on IFC and OFC
668 // IFC ROTATED CLIP +++++++++++++++++++++++++++++
669 for ( Int_t rod = rodIFC ; rod < rodIFC+1 ; rod++ ) { // loop over 1 to keep the "ignore"-functionality
671 erIntegral = ephiIntegral = ezIntegral = 0.0 ;
672 if ( fRotatedClipVoltA[0] == 0 && fgkZList[i] > 0) continue ;
673 if ( fRotatedClipVoltC[0] == 0 && fgkZList[i] < 0) continue ;
675 // Rotate to a position where we have maps and prepare to sum
676 phiPrime = phi - rod*kPhiSlicesPerSector*gridSizePhi ;
678 if ( phiPrime < 0 ) phiPrime += TMath::TwoPi() ; // Stay in range from 0 to TwoPi
680 if ( (phiPrime >= 0) && (phiPrime <= kPhiSlices*gridSizePhi) ) {
682 erIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
683 rlist, zedlist, philist, fLookUpBasic3ErOverEz );
684 ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
685 rlist, zedlist, philist, fLookUpBasic3EphiOverEz);
686 ezIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
687 rlist, zedlist, philist, fLookUpBasic3DeltaEz );
689 } else if ( (phiPrime <= TMath::TwoPi()) && (phiPrime >= (TMath::TwoPi()-kPhiSlices*gridSizePhi)) ){
691 phiPrime = TMath::TwoPi() - phiPrime ;
693 erIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
694 rlist, zedlist, philist, fLookUpBasic3ErOverEz );
695 ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
696 rlist, zedlist, philist, fLookUpBasic3EphiOverEz);
697 ezIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
698 rlist, zedlist, philist, fLookUpBasic3DeltaEz );
700 // Flip symmetry of phi integral for symmetric boundary conditions
701 if ( symmetry[2] == 1 ) ephiIntegral *= -1 ;
702 // Flip symmetry of r integral if anti-symmetric boundary conditions
703 if ( symmetry[2] == -1 ) erIntegral *= -1 ;
707 if ( fgkZList[i] > 0 ) {
708 erIntegralSum += fRotatedClipVoltA[0]*erIntegral ;
709 ephiIntegralSum += fRotatedClipVoltA[0]*ephiIntegral ;
710 ezIntegralSum += fRotatedClipVoltA[0]*ezIntegral;
711 } else if ( fgkZList[i] < 0 ) {
712 erIntegralSum += fRotatedClipVoltC[0]*erIntegral ;
713 ephiIntegralSum += fRotatedClipVoltC[0]*ephiIntegral ;
714 ezIntegralSum -= fRotatedClipVoltC[0]*ezIntegral;
718 // OFC: ROTATED CLIP +++++++++++++++++++++++++++++
719 for ( Int_t rod = rodOFC ; rod < rodOFC+1 ; rod++ ) { // loop over 1 to keep the "ignore"-functionality
721 erIntegral = ephiIntegral = ezIntegral = 0.0 ;
723 if ( fRotatedClipVoltA[1] == 0 && fgkZList[i] > 0) continue ;
724 if ( fRotatedClipVoltC[1] == 0 && fgkZList[i] < 0) continue ;
726 // Rotate to a position where we have maps and prepare to sum
727 phiPrime = phi - rod*kPhiSlicesPerSector*gridSizePhi ;
730 if ( phiPrime < 0 ) phiPrime += TMath::TwoPi() ; // Stay in range from 0 to TwoPi
732 if ( (phiPrime >= 0) && (phiPrime <= kPhiSlices*gridSizePhi) ) {
734 erIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
735 rlist, zedlist, philist, fLookUpBasic4ErOverEz );
736 ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
737 rlist, zedlist, philist, fLookUpBasic4EphiOverEz);
738 ezIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
739 rlist, zedlist, philist, fLookUpBasic4DeltaEz );
741 } else if ( (phiPrime <= TMath::TwoPi()) && (phiPrime >= (TMath::TwoPi()-kPhiSlices*gridSizePhi)) ){
743 phiPrime = TMath::TwoPi() - phiPrime ;
745 erIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
746 rlist, zedlist, philist, fLookUpBasic4ErOverEz );
747 ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
748 rlist, zedlist, philist, fLookUpBasic4EphiOverEz);
749 ezIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
750 rlist, zedlist, philist, fLookUpBasic4DeltaEz );
752 // Flip symmetry of phi integral for symmetric boundary conditions
753 if ( symmetry[3] == 1 ) ephiIntegral *= -1 ;
754 // Flip symmetry of r integral if anti-symmetric boundary conditions
755 if ( symmetry[3] == -1 ) erIntegral *= -1 ;
759 if ( fgkZList[i] > 0 ) {
760 erIntegralSum += fRotatedClipVoltA[1]*erIntegral ;
761 ephiIntegralSum += fRotatedClipVoltA[1]*ephiIntegral ;
762 ezIntegralSum += fRotatedClipVoltA[1]*ezIntegral;
763 } else if ( fgkZList[i] < 0 ) {
764 erIntegralSum += fRotatedClipVoltC[1]*erIntegral ;
765 ephiIntegralSum += fRotatedClipVoltC[1]*ephiIntegral ;
766 ezIntegralSum -= fRotatedClipVoltC[1]*ezIntegral;
770 // IFC Cooper rod shift +++++++++++++++++++++++++++++
771 for ( Int_t rod = 0 ; rod < 18 ; rod++ ) {
773 erIntegral = ephiIntegral = ezIntegral = 0.0 ;
775 if ( fCopperRodShiftA[rod] == 0 && fgkZList[i] > 0) continue ;
776 if ( fCopperRodShiftC[rod] == 0 && fgkZList[i] < 0) continue ;
778 // Rotate to a position where we have maps and prepare to sum
779 phiPrime = phi - rod*kPhiSlicesPerSector*gridSizePhi ;
781 if ( phiPrime < 0 ) phiPrime += TMath::TwoPi() ; // Stay in range from 0 to TwoPi
783 if ( (phiPrime >= 0) && (phiPrime <= kPhiSlices*gridSizePhi) ) {
785 erIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
786 rlist, zedlist, philist, fLookUpBasic5ErOverEz );
787 ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
788 rlist, zedlist, philist, fLookUpBasic5EphiOverEz);
789 ezIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
790 rlist, zedlist, philist, fLookUpBasic5DeltaEz );
792 } else if ( (phiPrime <= TMath::TwoPi()) && (phiPrime >= (TMath::TwoPi()-kPhiSlices*gridSizePhi)) ){
794 phiPrime = TMath::TwoPi() - phiPrime ;
796 erIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
797 rlist, zedlist, philist, fLookUpBasic5ErOverEz );
798 ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
799 rlist, zedlist, philist, fLookUpBasic5EphiOverEz);
800 ezIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
801 rlist, zedlist, philist, fLookUpBasic5DeltaEz );
803 // Flip symmetry of phi integral for symmetric boundary conditions
804 if ( symmetry[4] == 1 ) ephiIntegral *= -1 ;
805 // Flip symmetry of r integral if anti-symmetric boundary conditions
806 if ( symmetry[4] == -1 ) erIntegral *= -1 ;
810 if ( fgkZList[i] > 0 ) {
811 erIntegralSum += fCopperRodShiftA[rod]*erIntegral ;
812 ephiIntegralSum += fCopperRodShiftA[rod]*ephiIntegral ;
813 ezIntegralSum += fCopperRodShiftA[rod]*ezIntegral;
814 } else if ( fgkZList[i] < 0 ) {
815 erIntegralSum += fCopperRodShiftC[rod]*erIntegral ;
816 ephiIntegralSum += fCopperRodShiftC[rod]*ephiIntegral ;
817 ezIntegralSum -= fCopperRodShiftC[rod]*ezIntegral;
821 // OFC Cooper rod shift +++++++++++++++++++++++++++++
822 for ( Int_t rod = 18 ; rod < 36 ; rod++ ) {
824 erIntegral = ephiIntegral = ezIntegral = 0.0 ;
826 if ( fCopperRodShiftA[rod] == 0 && fgkZList[i] > 0) continue ;
827 if ( fCopperRodShiftC[rod] == 0 && fgkZList[i] < 0) continue ;
829 // Rotate to a position where we have maps and prepare to sum
830 phiPrime = phi - (rod-18)*kPhiSlicesPerSector*gridSizePhi ;
832 if ( phiPrime < 0 ) phiPrime += TMath::TwoPi() ; // Stay in range from 0 to TwoPi
834 if ( (phiPrime >= 0) && (phiPrime <= kPhiSlices*gridSizePhi) ) {
836 erIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
837 rlist, zedlist, philist, fLookUpBasic6ErOverEz );
838 ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
839 rlist, zedlist, philist, fLookUpBasic6EphiOverEz);
840 ezIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
841 rlist, zedlist, philist, fLookUpBasic6DeltaEz );
843 } else if ( (phiPrime <= TMath::TwoPi()) && (phiPrime >= (TMath::TwoPi()-kPhiSlices*gridSizePhi)) ){
845 phiPrime = TMath::TwoPi() - phiPrime ;
847 erIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
848 rlist, zedlist, philist, fLookUpBasic6ErOverEz );
849 ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
850 rlist, zedlist, philist, fLookUpBasic6EphiOverEz);
851 ezIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
852 rlist, zedlist, philist, fLookUpBasic6DeltaEz );
854 // Flip symmetry of phi integral for symmetric boundary conditions
855 if ( symmetry[5] == 1 ) ephiIntegral *= -1 ;
856 // Flip symmetry of r integral if anti-symmetric boundary conditions
857 if ( symmetry[5] == -1 ) erIntegral *= -1 ;
861 if ( fgkZList[i] > 0 ) {
862 erIntegralSum += fCopperRodShiftA[rod]*erIntegral ;
863 ephiIntegralSum += fCopperRodShiftA[rod]*ephiIntegral ;
864 ezIntegralSum += fCopperRodShiftA[rod]*ezIntegral;
865 } else if ( fgkZList[i] < 0 ) {
866 erIntegralSum += fCopperRodShiftC[rod]*erIntegral ;
867 ephiIntegralSum += fCopperRodShiftC[rod]*ephiIntegral ;
868 ezIntegralSum -= fCopperRodShiftC[rod]*ezIntegral;
872 // put the sum into the final lookup table
873 erOverEz(j,i) = erIntegralSum;
874 ephiOverEz(j,i) = ephiIntegralSum;
875 deltaEz(j,i) = ezIntegralSum;
877 // if (j==1) printf("%lf %lf %lf: %lf %lf %lf\n",r, phi, z, erIntegralSum,ephiIntegralSum,ezIntegralSum);
884 // clear the temporary arrays lists
885 for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
886 delete arrayofArrayV[k];
887 delete arrayofCharge[k];
895 void AliTPCFCVoltError3D::Print(const Option_t* option) const {
897 // Print function to check the settings of the Rod shifts and the rotated clips
898 // option=="a" prints the C0 and C1 coefficents for calibration purposes
901 TString opt = option; opt.ToLower();
902 printf("%s\n",GetTitle());
903 printf(" - ROD shifts (residual voltage settings): 40V correspond to 1mm of shift\n");
905 printf(" : A-side - (Rod & Strips) \n ");
906 for (Int_t i=0; i<36;i++) {
907 if (fRodVoltShiftA[i]!=0) {
908 printf("Rod%2d:%3.1fV ",i,fRodVoltShiftA[i]);
912 printf("-> all at 0 V\n");
918 printf(" : C-side - (Rod & Strips) \n ");
919 for (Int_t i=0; i<36;i++) {
920 if (fRodVoltShiftC[i]!=0) {
921 printf("Rod%2d:%3.1fV ",i,fRodVoltShiftC[i]);
925 printf("-> all at 0 V\n");
932 printf(" - Rotated clips (residual voltage settings): 40V correspond to 1mm of 'offset'\n");
933 if (fRotatedClipVoltA[0]!=0) { printf(" A side (IFC): HV Rod: %3.1f V \n",fRotatedClipVoltA[0]); count++; }
934 if (fRotatedClipVoltA[1]!=0) { printf(" A side (OFC): HV Rod: %3.1f V \n",fRotatedClipVoltA[1]); count++; }
935 if (fRotatedClipVoltC[0]!=0) { printf(" C side (IFC): HV Rod: %3.1f V \n",fRotatedClipVoltC[0]); count++; }
936 if (fRotatedClipVoltC[1]!=0) { printf(" C side (OFC): HV Rod: %3.1f V \n",fRotatedClipVoltC[1]); count++; }
938 printf(" -> no rotated clips \n");
941 printf(" - Copper ROD shifts (without strips):\n");
942 printf(" : A-side - (Copper Rod shift) \n ");
943 for (Int_t i=0; i<36;i++) {
944 if (fCopperRodShiftA[i]!=0) {
945 printf("Rod%2d:%3.1fV ",i,fCopperRodShiftA[i]);
949 printf("-> all at 0 V\n");
955 printf(" : C-side - (Copper Rod shift) \n ");
956 for (Int_t i=0; i<36;i++) {
957 if (fCopperRodShiftC[i]!=0) {
958 printf("Rod%2d:%3.1fV ",i,fCopperRodShiftC[i]);
962 printf("-> all at 0 V\n");
969 if (opt.Contains("a")) { // Print all details
970 printf(" - T1: %1.4f, T2: %1.4f \n",fT1,fT2);
971 printf(" - C1: %1.4f, C0: %1.4f \n",fC1,fC0);
974 if (!fInitLookUp) AliError("Lookup table was not initialized! You should do InitFCVoltError3D() ...");