]> git.uio.no Git - u/mrichter/AliRoot.git/blob - TPC/AliTPCROCVoltError3D.cxx
Global fit of alignment with E field distortion map (440 parameters)
[u/mrichter/AliRoot.git] / TPC / AliTPCROCVoltError3D.cxx
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
46 ClassImp(AliTPCROCVoltError3D)
47
48 AliTPCROCVoltError3D::AliTPCROCVoltError3D()
49   : AliTPCCorrection("ROCVoltErrors","ROC z alignment Errors"),
50     fC0(0.),fC1(0.),
51     fROCdisplacement(kTRUE),
52     fElectronArrivalCorrection(kTRUE),
53     fInitLookUp(kFALSE),
54     fROCDataFileName(""),  
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   }
69   fROCDataFileName="$ALICE_ROOT/TPC/Calib/maps/TPCROCdzSurvey.root";
70   SetROCDataFileName(fROCDataFileName.Data()); // initialization of fdzDataLinFit is included
71
72 }
73
74 AliTPCROCVoltError3D::~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
88 void AliTPCROCVoltError3D::SetROCData(TMatrixD * matrix){
89   //
90   // Set a z alignment map of the chambers not via a file, but
91   // directly via a TMatrix(72,3), where dz = p0 + p1*(lx-133.4) + p2*ly (all in cm)
92   //
93   if (!fdzDataLinFit) fdzDataLinFit=new TMatrixD(*matrix);
94   else *fdzDataLinFit = *matrix;
95 }
96
97
98 void 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
114   if (!fInitLookUp) InitROCVoltError3D();
115 }
116
117 void 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
134 void  AliTPCROCVoltError3D::SetROCDataFileName(const char * fname) {
135   //
136   // Set / load the ROC data (linear fit of ROC misalignments)
137   //
138
139   fROCDataFileName = fname;
140   
141   TFile f(fROCDataFileName.Data(),"READ");
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
163 void 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
174   Int_t   order     = 1 ;    // FIXME: hardcoded? Linear interpolation = 1, Quadratic = 2         
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
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
240 }
241
242 void 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
307             if ( j == (kColumns-1) ) 
308               arrayV(i,j) = 0.5*  ( GetROCVoltOffset( side, radius0, phi0+0.02 ) + GetROCVoltOffset( side, radius0, phi0-0.02 ) ) ;
309
310           }
311         }      
312         
313         for ( Int_t i = 1 ; i < kRows-1 ; i++ ) { 
314           for ( Int_t j = 1 ; j < kColumns-1 ; j++ ) {
315             charge(i,j)  =  0.0 ;
316           }
317         }
318       }      
319       
320       // Solve Poisson's equation in 3D cylindrical coordinates by relaxation technique
321       // Allow for different size grid spacing in R and Z directions
322       
323       PoissonRelaxation3D( arrayofArrayV, arrayofCharge, 
324                            arrayofEroverEz, arrayofEphioverEz, arrayofDeltaEz,
325                            kRows, kColumns, kPhiSlices, gridSizePhi, kIterations, 
326                            symmetry, fROCdisplacement) ;
327       
328       
329       //Interpolate results onto a custom grid which is used just for these calculations.
330       Double_t  r, phi, z ;
331       for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
332         phi = fgkPhiList[k] ;
333         
334         TMatrixD &erOverEz   =  *fLookUpErOverEz[k]  ;
335         TMatrixD &ephiOverEz =  *fLookUpEphiOverEz[k];
336         TMatrixD &deltaEz    =  *fLookUpDeltaEz[k]   ;
337         
338         for ( Int_t j = 0 ; j < kNZ ; j++ ) {
339
340           z = TMath::Abs(fgkZList[j]) ;  // Symmetric solution in Z that depends only on ABS(Z)
341   
342           if ( side == 0 &&  fgkZList[j] < 0 ) continue; // Skip rest of this loop if on the wrong side
343           if ( side == 1 &&  fgkZList[j] > 0 ) continue; // Skip rest of this loop if on the wrong side
344           
345           for ( Int_t i = 0 ; i < kNR ; i++ ) { 
346             r = fgkRList[i] ;
347
348             // Interpolate basicLookup tables; once for each rod, then sum the results
349             erOverEz(i,j)   = Interpolate3DTable(order, r, z, phi, kRows, kColumns, kPhiSlices, 
350                                                  rlist, zedlist, philist, arrayofEroverEz  );
351             ephiOverEz(i,j) = Interpolate3DTable(order, r, z, phi, kRows, kColumns, kPhiSlices,
352                                                  rlist, zedlist, philist, arrayofEphioverEz);
353             deltaEz(i,j)    = Interpolate3DTable(order, r, z, phi, kRows, kColumns, kPhiSlices,
354                                                  rlist, zedlist, philist, arrayofDeltaEz  );
355
356             if (side == 1)  deltaEz(i,j) = -  deltaEz(i,j); // negative coordinate system on C side
357
358           } // end r loop
359         }// end z loop
360       }// end phi loop
361
362       if ( side == 0 ) AliInfo(" A side done");
363       if ( side == 1 ) AliInfo(" C side done");
364     } // end side loop
365   }
366   
367   // clear the temporary arrays lists
368   for ( Int_t k = 0 ; k < kPhiSlices ; k++ )  {
369     delete arrayofArrayV[k];
370     delete arrayofCharge[k];
371     delete arrayofEroverEz[k];  
372     delete arrayofEphioverEz[k];
373     delete arrayofDeltaEz[k];
374   }
375  
376
377   fInitLookUp = kTRUE;
378
379 }
380
381
382 Float_t AliTPCROCVoltError3D::GetROCVoltOffset(Int_t side, Float_t r0, Float_t phi0) {
383   // 
384   // Returns the dz alignment data (in voltage equivalents) at 
385   // the given position
386   //
387
388   Float_t xp = r0*TMath::Cos(phi0);
389   Float_t yp = r0*TMath::Sin(phi0);
390   
391   // phi0 should be between 0 and 2pi 
392   if (phi0<0) phi0+=TMath::TwoPi();
393   Int_t roc = (Int_t)TMath::Floor((TMath::RadToDeg()*phi0)/20);
394   if (side==1) roc+=18; // C side
395   if (r0>132) roc+=36;  // OROC 
396   
397   // linear-plane data:  z = z0 + kx*lx + ky*ly (rotation in local coordinates)
398   TMatrixD &fitData = *fdzDataLinFit;
399
400   // local coordinates                                                          
401   Double_t secAlpha = TMath::DegToRad()*(10.+20.*(((Int_t)roc)%18));
402   Float_t lx = xp*TMath::Cos(secAlpha)+yp*TMath::Sin(secAlpha);
403   Float_t ly = yp*TMath::Cos(secAlpha)-xp*TMath::Sin(secAlpha);
404
405   // reference of rotation in lx is at the intersection between OROC and IROC
406   // necessary, since the Fitprozedure is otherwise useless
407   
408   AliTPCROC * rocInfo = AliTPCROC::Instance();
409   Double_t lxRef  = (rocInfo->GetPadRowRadii(0,62)+rocInfo->GetPadRowRadii(36,0))/2;
410   
411   Float_t dz = fitData(roc,0)+fitData(roc,1)*(lx-lxRef) + fitData(roc,2)*ly; // value in cm
412
413   // aproximated Voltage-offset-aquivalent to the z misalignment
414   // (linearly scaled with the z position)
415   Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm; 
416   Float_t voltOff = dz*ezField;            // values in "Volt equivalents"
417
418   return voltOff;
419 }
420
421 TH2F * AliTPCROCVoltError3D::CreateHistoOfZAlignment(Int_t side, Int_t nx, Int_t ny) {
422   //
423   // return a simple histogramm containing the input to the poisson solver
424   // (z positions of the Read-out chambers, linearly interpolated)
425
426   char hname[100];
427   if (side==0) sprintf(hname,"survey_dz_Aside");
428   if (side==1) sprintf(hname,"survey_dz_Cside");
429
430   TH2F *h = new TH2F(hname,hname,nx,-250.,250.,ny,-250.,250.);
431
432   for (Int_t iy=1;iy<=ny;++iy) {
433     Double_t yp = h->GetYaxis()->GetBinCenter(iy);
434     for (Int_t ix=1;ix<=nx;++ix) {
435       Double_t xp = h->GetXaxis()->GetBinCenter(ix);
436     
437       Float_t r0 = TMath::Sqrt(xp*xp+yp*yp);
438       Float_t phi0 = TMath::ATan2(yp,xp); 
439    
440       Float_t dz = GetROCVoltOffset(side,r0,phi0); // in [volt]
441
442       Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm; 
443       dz = dz/ezField;    // in [cm]
444
445       if (85.<=r0 && r0<=245.) {
446         h->SetBinContent(ix,iy,dz); 
447       } else {
448         h->SetBinContent(ix,iy,0.);
449       }
450     }
451   }
452   
453   h->GetXaxis()->SetTitle("x [cm]");
454   h->GetYaxis()->SetTitle("y [cm]");
455   h->GetZaxis()->SetTitle("dz [cm]");
456   h->SetStats(0);
457   //  h->DrawCopy("colz");
458
459   return h;
460
461
462 void AliTPCROCVoltError3D::Print(const Option_t* option) const {
463   //
464   // Print function to check the settings of the Rod shifts and the rotated clips
465   // option=="a" prints the C0 and C1 coefficents for calibration purposes
466   //
467
468   TString opt = option; opt.ToLower();
469   printf("%s\n",GetTitle());
470   printf(" - z aligmnet of the TPC Read-Out chambers \n");
471   printf("   (linear interpolation within the chamber:  dz = z0 + kx*(lx-133) + ky*ly [cm] ) \n");
472   printf("   Info: Check the following data-file for more details: %s \n",fROCDataFileName.Data());
473
474   if (opt.Contains("a")) { // Print all details
475     TMatrixD &fitData = *fdzDataLinFit;
476     printf(" A side:  IROC   ROCX=(z0,kx,ky): \n");
477     for (Int_t roc = 0; roc<18; roc++) 
478       printf("ROC%d:(%.2e,%.2e,%.2e) ",roc,fitData(roc,0),fitData(roc,1),fitData(roc,2));
479     printf("\n A side:  OROC   ROCX=(z0,kx,ky): \n");
480     for (Int_t roc = 36; roc<54; roc++) 
481       printf("ROC%d:(%.2e,%.2e,%.2e) ",roc,fitData(roc,0),fitData(roc,1),fitData(roc,2));
482     printf("\n C side:  IROC   ROCX=(z0,kx,ky): \n");
483     for (Int_t roc = 18; roc<36; roc++) 
484       printf("ROC%d:(%.2e,%.2e,%.2e) ",roc,fitData(roc,0),fitData(roc,1),fitData(roc,2));
485     printf("\n C side:  OROC   ROCX=(z0,kx,ky): \n");
486     for (Int_t roc = 54; roc<72; roc++) 
487       printf("ROC%d:(%.2e,%.2e,%.2e) ",roc,fitData(roc,0),fitData(roc,1),fitData(roc,2));
488     printf("\n\n");
489     printf(" - T1: %1.4f, T2: %1.4f \n",fT1,fT2);
490     printf(" - C1: %1.4f, C0: %1.4f \n",fC1,fC0);
491   }
492
493   if (!fInitLookUp) AliError("Lookup table was not initialized! You should do InitROCVoltError3D() ...");
494
495 }