M AliTPCcalibCalib.cxx - start refit form TPC out instead of track...
[u/mrichter/AliRoot.git] / TPC / AliTPCGGVoltError.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 // AliTPCGGVoltError class                                                //
19 // The class calculates the electric field and space point distortions    //
20 // due a Gating Grid (GG) Error voltage. It uses the exact calculation    //
21 // technique based on bessel functions. (original code from STAR)         //
22 // The class allows "effective Omega Tau" corrections.                    // 
23 //                                                                        //
24 // date: 27/04/2010                                                       //
25 // Authors: Jim Thomas, Stefan Rossegger, Magnus Mager                    //
26 //                                                                        //
27 // Example usage:                                                         //
28 //  AliTPCGGVoltError GGerror;                                            //
29 //  GGerror.SetOmegaTauT1T2(0.32,1.,1.); // values ideally from OCDB      //
30 //  GGerror.SetDeltaVGGA(50.);           // voltage offset A-side         //
31 //  GGerror.SetDeltaVGGC(50.);           // voltage offset C-side         //
32 //  GGerror.InitGGVoltErrorDistortion(); // initialization of the look up //
33 //  // plot dRPhi distortions ...                                         //
34 //  GGerror.CreateHistoDRPhiinZR(1.,100,100)->Draw("surf2");              //
35 ////////////////////////////////////////////////////////////////////////////
36
37
38
39 #include "AliTPCGGVoltError.h"
40 #include <TMath.h>
41
42 AliTPCGGVoltError::AliTPCGGVoltError()
43   : AliTPCCorrection("GGVoltError","GatingGrid (GG) Voltage Error"),
44     fC0(0.),fC1(0.),
45     fDeltaVGGA(0.),fDeltaVGGC(0.)
46 {
47   //
48   // default constructor
49   //
50 }
51
52 AliTPCGGVoltError::~AliTPCGGVoltError() {
53   //
54   // default destructor
55   //
56 }
57
58 void AliTPCGGVoltError::Init() {
59   //
60   // Initialization funtion (not used at the moment)
61   //
62   
63   // Set default parameters
64   // FIXME: Ask the database for these entries
65   
66   Double_t vdrift = 2.6; // [cm/us]   // From dataBase: to be updated: per second (ideally)
67   Double_t bzField = -0.5; // [Tesla] // From dataBase: to be updated: per run
68
69   Double_t ezField = 400; // [V/cm]   // to be updated: never (hopefully)
70   Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ; 
71
72   // Correction Terms for effective omegaTau; obtained by a laser calibration run
73   Double_t t1 = 0.9;   // ideally from database
74   Double_t t2 = 1.5;   // ideally from database
75
76   SetOmegaTauT1T2(wt,t1,t2);
77
78   SetDeltaVGGA(0.0);//  ideally from the database
79   SetDeltaVGGC(0.0);//  ideally from the database
80 }
81
82 void AliTPCGGVoltError::Update(const TTimeStamp &/*timeStamp*/) {
83   //
84   // Update function 
85   //
86
87   Double_t vdrift = 2.6; // [cm/us]   // From dataBase: to be updated: per second (ideally)
88   Double_t bzField = -0.5; // [Tesla] // From dataBase: to be updated: per run
89
90   Double_t ezField = 400; // [V/cm]   // to be updated: never (hopefully)
91   Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ; 
92
93   // Correction Terms for effective omegaTau; obtained by a laser calibration run
94   Double_t t1 = 0.9;   // ideally from database
95   Double_t t2 = 1.5;   // ideally from database
96
97   SetOmegaTauT1T2(wt,t1,t2);
98 }
99
100
101
102 void AliTPCGGVoltError::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {
103
104   //
105   // Gated Grid Voltage Error
106   //
107   // Calculates the effect of having an incorrect voltage on the A or C end plate Gated Grids.
108   //
109   // Electrostatic Equations from StarNote SN0253 by Howard Wieman.
110   //
111   
112   Int_t   order     = 1 ;               // FIXME: hardcoded? Linear interpolation = 1, Quadratic = 2         
113  
114   Double_t intEr, intEphi ;
115   Double_t r, phi, z ;
116   Int_t    sign ;
117
118   Double_t deltaVGG;
119   
120   r   = TMath::Sqrt( x[0]*x[0] + x[1]*x[1] );
121   phi = TMath::ATan2(x[1],x[0]);
122   if ( phi < 0 ) phi += TMath::TwoPi();                   // Table uses phi from 0 to 2*Pi
123   z   = x[2] ;
124
125   if ( (roc%36) < 18 ) {
126     sign =  1; 
127     deltaVGG = fDeltaVGGA;           // (TPC End A)
128   } else {
129     sign = -1;                       // (TPC End C)
130     deltaVGG = fDeltaVGGC; 
131   }
132
133   if ( sign==1  && z <  fgkZOffSet ) z =  fgkZOffSet;    // Protect against discontinuity at CE
134   if ( sign==-1 && z > -fgkZOffSet ) z = -fgkZOffSet;    // Protect against discontinuity at CE
135
136   Interpolate2DEdistortion( order, r, z, fGGVoltErrorER, intEr );
137   intEphi = 0.0;  // Efield is symmetric in phi
138
139   // Calculate distorted position
140   if ( r > 0.0 ) {
141     phi =  phi + deltaVGG*( fC0*intEphi - fC1*intEr ) / r;      
142     r   =  r   + deltaVGG*( fC0*intEr   + fC1*intEphi );  
143   }
144   
145   // Calculate correction in cartesian coordinates
146   dx[0] = r * TMath::Cos(phi) - x[0];
147   dx[1] = r * TMath::Sin(phi) - x[1]; 
148   dx[2] = 0.; // z distortion not implemented (1st order distortions)
149
150 }
151
152
153 Float_t AliTPCGGVoltError::GetIntErOverEz(const Float_t x[],const Short_t roc) {
154   //
155   // This function is purely for calibration purposes
156   // Calculates the integral (int Er/Ez dz) for the setted GG voltage offset 
157   // 
158   
159   Int_t   order     = 1 ;     // FIXME: so far hardcoded? Linear interpolation = 1, Quadratic = 2         
160   
161   Double_t intEr;
162   Double_t r, phi, z ;
163   Int_t    sign ;
164   
165   Double_t deltaVGG;
166   
167   r   = TMath::Sqrt( x[0]*x[0] + x[1]*x[1] );
168   phi = TMath::ATan2(x[1],x[0]);
169   if ( phi < 0 ) phi += TMath::TwoPi();        // Table uses phi from 0 to 2*Pi
170   z   = x[2] ;
171
172   if ( (roc%36) < 18 ) {
173     sign =  1; 
174     deltaVGG = fDeltaVGGA;           // (TPC End A)
175   } else {
176     sign = -1;                       // (TPC End C)
177     deltaVGG = fDeltaVGGC; 
178   }
179
180   if ( sign==1  && z <  fgkZOffSet ) z =  fgkZOffSet;    // Protect against discontinuity at CE
181   if ( sign==-1 && z > -fgkZOffSet ) z = -fgkZOffSet;    // Protect against discontinuity at CE
182
183   Interpolate2DEdistortion(order, r, z, fGGVoltErrorER, intEr );
184
185   return (intEr*deltaVGG);
186
187 }
188
189 void AliTPCGGVoltError::InitGGVoltErrorDistortion() {
190   //
191   // Initialization of the Lookup table which contains the solutions of the GG Error problem
192   //
193
194   Double_t r,z;
195   Int_t nterms = 100 ;
196   for ( Int_t i = 0 ; i < kNZ ; ++i ) {
197     z = fgkZList[i] ;
198     for ( Int_t j = 0 ; j < kNR ; ++j ) {
199       r = fgkRList[j] ;
200       fGGVoltErrorER[i][j] = 0.0 ;          
201       Double_t intz = 0.0 ;
202       for ( Int_t n = 1 ; n < nterms ; ++n ) {
203         Double_t k    =  n * TMath::Pi() / fgkTPC_Z0 ;
204         Double_t ein  =  0 ;                    // Error potential on the IFC
205         Double_t eout =  0 ;                    // Error potential on the OFC
206         if ( z < 0 ) {
207           ein   =  -2.0 / ( k * (fgkCathodeV - fgkGG) ) ;       
208           eout  =  -2.0 / ( k * (fgkCathodeV - fgkGG) ) ;       
209         }
210         if ( z == 0 ) continue ;
211         if ( z > 0 ) {
212           ein   =  -2.0 / ( k * (fgkCathodeV - fgkGG) ) ;       
213           eout  =  -2.0 / ( k * (fgkCathodeV - fgkGG) ) ;       
214         }
215         Double_t an   =  ein  * TMath::BesselK0( k*fgkOFCRadius ) - eout * TMath::BesselK0( k*fgkIFCRadius ) ;
216         Double_t bn   =  eout * TMath::BesselI0( k*fgkIFCRadius ) - ein  * TMath::BesselI0( k*fgkOFCRadius ) ;
217         Double_t numerator =
218           an * TMath::BesselI1( k*r ) - bn * TMath::BesselK1( k*r ) ;
219         Double_t denominator =
220           TMath::BesselK0( k*fgkOFCRadius ) * TMath::BesselI0( k*fgkIFCRadius ) -
221           TMath::BesselK0( k*fgkIFCRadius ) * TMath::BesselI0( k*fgkOFCRadius ) ;
222         Double_t zterm = TMath::Cos( k*(fgkTPC_Z0-TMath::Abs(z)) ) - 1 ;
223         intz += zterm * numerator / denominator ;
224         // Assume series converges, break if small terms
225         if ( n>10 && TMath::Abs(intz)*1.e-10 > TMath::Abs(numerator/denominator) ) break;   
226       }
227       fGGVoltErrorER[i][j] = (Double_t) intz ;
228
229     }
230   }
231 }
232
233
234
235 void AliTPCGGVoltError::Print(Option_t* option) const {
236   //
237   // Print function to check the settings (e.g. voltage offsets)
238   // option=="a" prints the C0 and C1 coefficents for calibration purposes
239   //
240
241   TString opt = option; opt.ToLower();
242   printf("%s\n",GetTitle());
243   printf(" - GG Voltage offset: A-side: %3.1f V, C-side: %3.1f V \n",fDeltaVGGA,fDeltaVGGC);  
244   if (opt.Contains("a")) { // Print all details
245     printf(" - C1: %1.4f, C0: %1.4f \n",fC1,fC0);
246   }    
247
248
249   
250 }