New more precise light yields from Raphael. It is time to put them as an object in...
[u/mrichter/AliRoot.git] / TPC / AliTPCROCVoltError3D.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// AliTPCROCVoltError3D class //
19// The class calculates the space point distortions due to residual voltage //
20// errors on Read Out Chambers 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// misalignment and the vd dependency on the residual drift field //
26// //
27// date: 08/08/2010 //
28// Authors: Jim Thomas, Stefan Rossegger //
29// //
30// Example usage : //
31// AliTPCROCVoltError3D ROCerror; //
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#include "TFile.h"
41
42#include "TMath.h"
43#include "AliTPCROC.h"
44#include "AliTPCROCVoltError3D.h"
45
46ClassImp(AliTPCROCVoltError3D)
47
48AliTPCROCVoltError3D::AliTPCROCVoltError3D()
49 : AliTPCCorrection("ROCVoltErrors","ROC z alignment Errors"),
50 fC0(0.),fC1(0.),
51 fROCdisplacement(kTRUE),
52 fInitLookUp(kFALSE),
2b68ab9c 53 fROCDataFileName((char*)"$(ALICE_ROOT)/TPC/Calib/maps/TPCROCdzSurvey.root"),
c9cbd2f2 54 fdzDataLinFit(0)
55{
56 //
57 // default constructor
58 //
59
60 // Array which will contain the solution according to the setted boundary conditions
61 // main input: z alignment of the Read Out chambers
62 // see InitROCVoltError3D() function
63 for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
64 fLookUpErOverEz[k] = new TMatrixD(kNR,kNZ);
65 fLookUpEphiOverEz[k] = new TMatrixD(kNR,kNZ);
66 fLookUpDeltaEz[k] = new TMatrixD(kNR,kNZ);
67 }
68
69 SetROCDataFileName(fROCDataFileName); // initialization of fdzDataLinFit is included
70
71}
72
73AliTPCROCVoltError3D::~AliTPCROCVoltError3D() {
74 //
75 // destructor
76 //
77
78 for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
79 delete fLookUpErOverEz[k];
80 delete fLookUpEphiOverEz[k];
81 delete fLookUpDeltaEz[k];
82 }
83
84 delete fdzDataLinFit;
85}
86
2bbac918 87void AliTPCROCVoltError3D::SetDZMap(TMatrixD * matrix){
88 //
89 //
90 //
91 if (!fdzDataLinFit) fdzDataLinFit=new TMatrixD(*matrix);
92 else *fdzDataLinFit = *matrix;
93}
94
95
c9cbd2f2 96void AliTPCROCVoltError3D::Init() {
97 //
98 // Initialization funtion
99 //
100
101 AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
102 if (!magF) AliError("Magneticd field - not initialized");
103 Double_t bzField = magF->SolenoidField()/10.; //field in T
104 AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
105 if (!param) AliError("Parameters - not initialized");
106 Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us] // From dataBase: to be updated: per second (ideally)
107 Double_t ezField = 400; // [V/cm] // to be updated: never (hopefully)
108 Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
109 // Correction Terms for effective omegaTau; obtained by a laser calibration run
110 SetOmegaTauT1T2(wt,fT1,fT2);
111
35ae345f 112 if (!fInitLookUp) InitROCVoltError3D();
c9cbd2f2 113}
114
115void AliTPCROCVoltError3D::Update(const TTimeStamp &/*timeStamp*/) {
116 //
117 // Update function
118 //
119 AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
120 if (!magF) AliError("Magneticd field - not initialized");
121 Double_t bzField = magF->SolenoidField()/10.; //field in T
122 AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
123 if (!param) AliError("Parameters - not initialized");
124 Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us] // From dataBase: to be updated: per second (ideally)
125 Double_t ezField = 400; // [V/cm] // to be updated: never (hopefully)
126 Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
127 // Correction Terms for effective omegaTau; obtained by a laser calibration run
128 SetOmegaTauT1T2(wt,fT1,fT2);
129
130}
131
132void AliTPCROCVoltError3D::SetROCDataFileName(char *const fname) {
133 //
134 // Set / load the ROC data (linear fit of ROC misalignments)
135 //
136
137 fROCDataFileName = fname;
138
139 TFile f(fROCDataFileName,"READ");
140 TMatrixD *m = (TMatrixD*) f.Get("dzSurveyLinFitData");
141 TMatrixD &mf = *m;
142
143 // prepare some space
144
145 if (fdzDataLinFit) delete fdzDataLinFit;
146 fdzDataLinFit = new TMatrixD(72,3);
147 TMatrixD &dataIntern = *fdzDataLinFit;
148
149 for (Int_t iroc=0;iroc<72;iroc++) {
150 dataIntern(iroc,0) = mf(iroc,0); // z0 offset
151 dataIntern(iroc,1) = mf(iroc,1); // slope in x
152 dataIntern(iroc,2) = mf(iroc,2); // slope in y
153 }
154
155 f.Close();
156
157 fInitLookUp = kFALSE;
158
159}
160
161void AliTPCROCVoltError3D::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {
162 //
163 // Calculates the correction due e.g. residual voltage errors on the TPC boundaries
164 //
165
166 if (!fInitLookUp) {
167 AliInfo("Lookup table was not initialized! Perform the inizialisation now ...");
168 InitROCVoltError3D();
169 return;
170 }
171
172 Int_t order = 1 ; // FIXME: hardcoded? Linear interpolation = 1, Quadratic = 2
173
174 Double_t intEr, intEphi, intDeltaEz;
175 Double_t r, phi, z ;
176 Int_t sign;
177
178 r = TMath::Sqrt( x[0]*x[0] + x[1]*x[1] ) ;
179 phi = TMath::ATan2(x[1],x[0]) ;
180 if ( phi < 0 ) phi += TMath::TwoPi() ; // Table uses phi from 0 to 2*Pi
181 z = x[2] ; // Create temporary copy of x[2]
182
183 if ( (roc%36) < 18 ) {
184 sign = 1; // (TPC A side)
185 } else {
186 sign = -1; // (TPC C side)
187 }
188
189 if ( sign==1 && z < fgkZOffSet ) z = fgkZOffSet; // Protect against discontinuity at CE
190 if ( sign==-1 && z > -fgkZOffSet ) z = -fgkZOffSet; // Protect against discontinuity at CE
191
192
193 if ( (sign==1 && z<0) || (sign==-1 && z>0) ) // just a consistency check
194 AliError("ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!");
195
196 // Get the Er and Ephi field integrals plus the integral over DeltaEz
197 intEr = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
198 fgkRList, fgkZList, fgkPhiList, fLookUpErOverEz );
199 intEphi = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
200 fgkRList, fgkZList, fgkPhiList, fLookUpEphiOverEz);
201 intDeltaEz = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
202 fgkRList, fgkZList, fgkPhiList, fLookUpDeltaEz );
203
204 // printf("%lf %lf %lf\n",intEr,intEphi,intDeltaEz);
205
206 // Calculate distorted position
207 if ( r > 0.0 ) {
208 phi = phi + ( fC0*intEphi - fC1*intEr ) / r;
209 r = r + ( fC0*intEr + fC1*intEphi );
210 }
211
212 // Calculate correction in cartesian coordinates
213 dx[0] = r * TMath::Cos(phi) - x[0];
214 dx[1] = r * TMath::Sin(phi) - x[1];
215 dx[2] = intDeltaEz; // z distortion - (internally scaled with driftvelocity dependency
216 // on the Ez field plus the actual ROC misalignment (if set TRUE)
217
218}
219
220void AliTPCROCVoltError3D::InitROCVoltError3D() {
221 //
222 // Initialization of the Lookup table which contains the solutions of the
223 // Dirichlet boundary problem
224 // Calculation of the single 3D-Poisson solver is done just if needed
225 // (see basic lookup tables in header file)
226 //
227
228 const Int_t order = 1 ; // Linear interpolation = 1, Quadratic = 2
229 const Float_t gridSizeR = (fgkOFCRadius-fgkIFCRadius) / (kRows-1) ;
230 const Float_t gridSizeZ = fgkTPCZ0 / (kColumns-1) ;
231 const Float_t gridSizePhi = TMath::TwoPi() / ( 18.0 * kPhiSlicesPerSector);
232
233 // temporary arrays to create the boundary conditions
234 TMatrixD *arrayofArrayV[kPhiSlices], *arrayofCharge[kPhiSlices] ;
235 TMatrixD *arrayofEroverEz[kPhiSlices], *arrayofEphioverEz[kPhiSlices], *arrayofDeltaEz[kPhiSlices] ;
236
237 for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
238 arrayofArrayV[k] = new TMatrixD(kRows,kColumns) ;
239 arrayofCharge[k] = new TMatrixD(kRows,kColumns) ;
240 arrayofEroverEz[k] = new TMatrixD(kRows,kColumns) ;
241 arrayofEphioverEz[k] = new TMatrixD(kRows,kColumns) ;
242 arrayofDeltaEz[k] = new TMatrixD(kRows,kColumns) ;
243 }
244
245 // list of point as used in the poisson relation and the interpolation (during sum up)
246 Double_t rlist[kRows], zedlist[kColumns] , philist[kPhiSlices];
247 for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
248 philist[k] = gridSizePhi * k;
249 for ( Int_t i = 0 ; i < kRows ; i++ ) {
250 rlist[i] = fgkIFCRadius + i*gridSizeR ;
251 for ( Int_t j = 0 ; j < kColumns ; j++ ) { // Fill Vmatrix with Boundary Conditions
252 zedlist[j] = j * gridSizeZ ;
253 }
254 }
255 }
256
257 // ==========================================================================
258 // Solve Poisson's equation in 3D cylindrical coordinates by relaxation technique
259 // Allow for different size grid spacing in R and Z directions
260
261 const Int_t symmetry = 0;
262
263 // Set bondaries and solve Poisson's equation --------------------------
264
265 if ( !fInitLookUp ) {
266
267 AliInfo(Form("Solving the poisson equation (~ %d sec)",2*10*(int)(kPhiSlices/10)));
268
269 for ( Int_t side = 0 ; side < 2 ; side++ ) { // Solve Poisson3D twice; once for +Z and once for -Z
270
271 for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
272 TMatrixD &arrayV = *arrayofArrayV[k] ;
273 TMatrixD &charge = *arrayofCharge[k] ;
274
275 //Fill arrays with initial conditions. V on the boundary and Charge in the volume.
276 for ( Int_t i = 0 ; i < kRows ; i++ ) {
277 for ( Int_t j = 0 ; j < kColumns ; j++ ) { // Fill Vmatrix with Boundary Conditions
278 arrayV(i,j) = 0.0 ;
279 charge(i,j) = 0.0 ;
280
281 Float_t radius0 = rlist[i] ;
282 Float_t phi0 = gridSizePhi * k ;
283
284 // To avoid problems at sector boundaries, use an average of +- 1 degree from actual phi location
285 if ( j == (kColumns-1) )
286 arrayV(i,j) = 0.5* ( GetROCVoltOffset( side, radius0, phi0+0.02 ) + GetROCVoltOffset( side, radius0, phi0-0.02 ) ) ;
287
288 }
289 }
290
291 for ( Int_t i = 1 ; i < kRows-1 ; i++ ) {
292 for ( Int_t j = 1 ; j < kColumns-1 ; j++ ) {
293 charge(i,j) = 0.0 ;
294 }
295 }
296 }
297
298 // Solve Poisson's equation in 3D cylindrical coordinates by relaxation technique
299 // Allow for different size grid spacing in R and Z directions
300
301 PoissonRelaxation3D( arrayofArrayV, arrayofCharge,
302 arrayofEroverEz, arrayofEphioverEz, arrayofDeltaEz,
303 kRows, kColumns, kPhiSlices, gridSizePhi, kIterations,
304 symmetry, fROCdisplacement) ;
305
306
307 //Interpolate results onto a custom grid which is used just for these calculations.
308 Double_t r, phi, z ;
309 for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
310 phi = fgkPhiList[k] ;
311
312 TMatrixD &erOverEz = *fLookUpErOverEz[k] ;
313 TMatrixD &ephiOverEz = *fLookUpEphiOverEz[k];
314 TMatrixD &deltaEz = *fLookUpDeltaEz[k] ;
315
316 for ( Int_t j = 0 ; j < kNZ ; j++ ) {
317
318 z = TMath::Abs(fgkZList[j]) ; // Symmetric solution in Z that depends only on ABS(Z)
319
320 if ( side == 0 && fgkZList[j] < 0 ) continue; // Skip rest of this loop if on the wrong side
321 if ( side == 1 && fgkZList[j] > 0 ) continue; // Skip rest of this loop if on the wrong side
322
323 for ( Int_t i = 0 ; i < kNR ; i++ ) {
324 r = fgkRList[i] ;
325
326 // Interpolate basicLookup tables; once for each rod, then sum the results
327 erOverEz(i,j) = Interpolate3DTable(order, r, z, phi, kRows, kColumns, kPhiSlices,
328 rlist, zedlist, philist, arrayofEroverEz );
329 ephiOverEz(i,j) = Interpolate3DTable(order, r, z, phi, kRows, kColumns, kPhiSlices,
330 rlist, zedlist, philist, arrayofEphioverEz);
331 deltaEz(i,j) = Interpolate3DTable(order, r, z, phi, kRows, kColumns, kPhiSlices,
332 rlist, zedlist, philist, arrayofDeltaEz );
333
334 if (side == 1) deltaEz(i,j) = - deltaEz(i,j); // negative coordinate system on C side
335
336 } // end r loop
337 }// end z loop
338 }// end phi loop
339
340 if ( side == 0 ) AliInfo(" A side done");
341 if ( side == 1 ) AliInfo(" C side done");
342 } // end side loop
343 }
344
345 // clear the temporary arrays lists
346 for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
347 delete arrayofArrayV[k];
348 delete arrayofCharge[k];
349 delete arrayofEroverEz[k];
350 delete arrayofEphioverEz[k];
351 delete arrayofDeltaEz[k];
352 }
353
354
355 fInitLookUp = kTRUE;
356
357}
358
359
360Float_t AliTPCROCVoltError3D::GetROCVoltOffset(Int_t side, Float_t r0, Float_t phi0) {
361 //
362 // Returns the dz alignment data (in voltage equivalents) at
363 // the given position
364 //
365
366 Float_t xp = r0*TMath::Cos(phi0);
367 Float_t yp = r0*TMath::Sin(phi0);
368
369 // phi0 should be between 0 and 2pi
370 if (phi0<0) phi0+=TMath::TwoPi();
371 Int_t roc = (Int_t)TMath::Floor((TMath::RadToDeg()*phi0)/20);
372 if (side==1) roc+=18; // C side
373 if (r0>132) roc+=36; // OROC
374
375 // linear-plane data: z = z0 + kx*x + ky*y
376 TMatrixD &fitData = *fdzDataLinFit;
377 Float_t dz = fitData(roc,0)+fitData(roc,1)*xp + fitData(roc,2)*yp; // value in cm
378
379 // aproximated Voltage-offset-aquivalent to the z misalignment
380 // (linearly scaled with the z position)
381 Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm;
382 Float_t voltOff = dz*ezField; // values in "Volt equivalents"
383
384 return voltOff;
385}
386
387TH2F * AliTPCROCVoltError3D::CreateHistoOfZSurvey(Int_t side, Int_t nx, Int_t ny) {
388 //
389 // return a simple histogramm containing the input to the poisson solver
390 // (z positions of the Read-out chambers, linearly interpolated)
391
392 char hname[100];
393 if (side==0) sprintf(hname,"survey_dz_Aside");
394 if (side==1) sprintf(hname,"survey_dz_Cside");
395
396 TH2F *h = new TH2F(hname,hname,nx,-250.,250.,ny,-250.,250.);
397
398 for (Int_t iy=1;iy<=ny;++iy) {
399 Double_t yp = h->GetYaxis()->GetBinCenter(iy);
400 for (Int_t ix=1;ix<=nx;++ix) {
401 Double_t xp = h->GetXaxis()->GetBinCenter(ix);
402
403 Float_t r0 = TMath::Sqrt(xp*xp+yp*yp);
404 Float_t phi0 = TMath::ATan2(yp,xp);
405
406 Float_t dz = GetROCVoltOffset(side,r0,phi0); // in [volt]
407
408 Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm;
409 dz = dz/ezField; // in [cm]
410
411 if (85.<=r0 && r0<=245.) {
412 h->SetBinContent(ix,iy,dz);
413 } else {
414 h->SetBinContent(ix,iy,0.);
415 }
416 }
417 }
418
419 h->GetXaxis()->SetTitle("x [cm]");
420 h->GetYaxis()->SetTitle("y [cm]");
421 h->GetZaxis()->SetTitle("dz [cm]");
422 h->SetStats(0);
423 // h->DrawCopy("colz");
424
425 return h;
426}
427
428void AliTPCROCVoltError3D::Print(const Option_t* option) const {
429 //
430 // Print function to check the settings of the Rod shifts and the rotated clips
431 // option=="a" prints the C0 and C1 coefficents for calibration purposes
432 //
433
434 TString opt = option; opt.ToLower();
435 printf("%s\n",GetTitle());
436 printf(" - Voltage settings on the TPC Read-Out chambers - linearly interpolated\n");
437 printf(" info: Check the following data-file for more details: %s \n",fROCDataFileName);
438
439 if (opt.Contains("a")) { // Print all details
440 printf(" - T1: %1.4f, T2: %1.4f \n",fT1,fT2);
441 printf(" - C1: %1.4f, C0: %1.4f \n",fC1,fC0);
442 }
443
444 if (!fInitLookUp) AliError("Lookup table was not initialized! You should do InitROCVoltError3D() ...");
445
446}