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