]>
Commit | Line | Data |
---|---|---|
43e55a68 | 1 | // |
2 | // | |
3 | // This class is based on Jim Thomas's implementation of the ExB effects. | |
4 | // It was shortened and made to comply with AliROOT coding conventions, | |
5 | // but not changed in functionallity. | |
6 | // | |
7 | // | |
8 | ||
9 | ||
10 | /*********************************************************************** | |
11 | * | |
12 | * Author: Jim Thomas 11/4/2009 | |
13 | * | |
14 | *********************************************************************** | |
15 | * | |
16 | * Description: | |
17 | * Utilities for handling SpaceCharge and a few other distortions | |
18 | * that may occur in the ALICE TPC | |
19 | * | |
20 | ***********************************************************************/ | |
21 | ||
22 | #include "AliTPCDistortions.h" | |
23 | #include <AliMagF.h> | |
24 | #include <TMath.h> | |
25 | ||
26 | // | |
27 | // Standard maps for E and B Field Distortions | |
28 | // | |
29 | // Note the careful steps in Z around the Central Electrode due to possible discontinuities at CE. | |
30 | // Needed for interpolation tools on the grid. Also note that whenever we interpolate this grid, | |
31 | // we explicitly do not allow Z to get closer to CE than 0.2 cm. This gives three points for | |
32 | // quadratic interpolation in all cases (if you need it). | |
33 | ||
34 | const Double_t AliTPCDistortions::fgkRList[AliTPCDistortions::kNR] = { 84.0, 84.5, 85.0, 85.5, 86.0, 87.0, 88.0, | |
35 | 90.0, 92.0, 94.0, 96.0, 98.0, 100.0, 102.0, 104.0, 106.0, 108.0, | |
36 | 110.0, 112.0, 114.0, 116.0, 118.0, 120.0, 122.0, 124.0, 126.0, 128.0, | |
37 | 130.0, 132.0, 134.0, 136.0, 138.0, 140.0, 142.0, 144.0, 146.0, 148.0, | |
38 | 150.0, 152.0, 154.0, 156.0, 158.0, 160.0, 162.0, 164.0, 166.0, 168.0, | |
39 | 170.0, 172.0, 174.0, 176.0, 178.0, 180.0, 182.0, 184.0, 186.0, 188.0, | |
40 | 190.0, 192.0, 194.0, 196.0, 198.0, 200.0, 202.0, 204.0, 206.0, 208.0, | |
41 | 210.0, 212.0, 214.0, 216.0, 218.0, 220.0, 222.0, 224.0, 226.0, 228.0, | |
42 | 230.0, 232.0, 234.0, 236.0, 238.0, 240.0, 242.0, 244.0, 246.0, 248.0, | |
43 | 249.0, 249.5, 250.0, 251.5, 252.0 } ; | |
44 | ||
45 | const Double_t AliTPCDistortions::fgkPhiList[AliTPCDistortions::kNPhi] = { 0.0000, | |
46 | 2.0*TMath::Pi()* 1.0/18, 2.0*TMath::Pi()* 2.0/18, 2.0*TMath::Pi()* 3.0/18, | |
47 | 2.0*TMath::Pi()* 4.0/18, 2.0*TMath::Pi()* 5.0/18, 2.0*TMath::Pi()* 6.0/18, | |
48 | 2.0*TMath::Pi()* 7.0/18, 2.0*TMath::Pi()* 8.0/18, 2.0*TMath::Pi()* 9.0/18, | |
49 | 2.0*TMath::Pi()*10.0/18, 2.0*TMath::Pi()*11.0/18, 2.0*TMath::Pi()*12.0/18, | |
50 | 2.0*TMath::Pi()*13.0/18, 2.0*TMath::Pi()*14.0/18, 2.0*TMath::Pi()*15.0/18, | |
51 | 2.0*TMath::Pi()*16.0/18, 2.0*TMath::Pi()*17.0/18, 2.0*TMath::Pi() } ; | |
52 | // 19 planes of phi - 18+1 so can wrap around | |
53 | ||
54 | const Double_t AliTPCDistortions::fgkZList[AliTPCDistortions::kNZ] = | |
55 | { -249.5, -249.0, -248.5, -248.0, -247.0, -246.0, -245.0, -243.0, -242.0, -241.0, | |
56 | -240.0, -238.0, -236.0, -234.0, -232.0, -230.0, -228.0, -226.0, -224.0, -222.0, | |
57 | -220.0, -218.0, -216.0, -214.0, -212.0, -210.0, -208.0, -206.0, -204.0, -202.0, | |
58 | -200.0, -198.0, -196.0, -194.0, -192.0, -190.0, -188.0, -186.0, -184.0, -182.0, | |
59 | -180.0, -178.0, -176.0, -174.0, -172.0, -170.0, -168.0, -166.0, -164.0, -162.0, | |
60 | -160.0, -158.0, -156.0, -154.0, -152.0, -150.0, -148.0, -146.0, -144.0, -142.0, | |
61 | -140.0, -138.0, -136.0, -134.0, -132.0, -130.0, -128.0, -126.0, -124.0, -122.0, | |
62 | -120.0, -118.0, -116.0, -114.0, -112.0, -110.0, -108.0, -106.0, -104.0, -102.0, | |
63 | -100.0, -98.0, -96.0, -94.0, -92.0, -90.0, -88.0, -86.0, -84.0, -82.0, | |
64 | -80.0, -78.0, -76.0, -74.0, -72.0, -70.0, -68.0, -66.0, -64.0, -62.0, | |
65 | -60.0, -58.0, -56.0, -54.0, -52.0, -50.0, -48.0, -46.0, -44.0, -42.0, | |
66 | -40.0, -38.0, -36.0, -34.0, -32.0, -30.0, -28.0, -26.0, -24.0, -22.0, | |
67 | -20.0, -18.0, -16.0, -14.0, -12.0, -10.0, -8.0, -6.0, -4.0, -2.0, | |
68 | -1.0, -0.5, -0.2, -0.1, -0.05, 0.05, 0.1, 0.2, 0.5, 1.0, | |
69 | 2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0, | |
70 | 22.0, 24.0, 26.0, 28.0, 30.0, 32.0, 34.0, 36.0, 38.0, 40.0, | |
71 | 42.0, 44.0, 46.0, 48.0, 50.0, 52.0, 54.0, 56.0, 58.0, 60.0, | |
72 | 62.0, 64.0, 66.0, 68.0, 70.0, 72.0, 74.0, 76.0, 78.0, 80.0, | |
73 | 82.0, 84.0, 86.0, 88.0, 90.0, 92.0, 94.0, 96.0, 98.0, 100.0, | |
74 | 102.0, 104.0, 106.0, 108.0, 110.0, 112.0, 114.0, 116.0, 118.0, 120.0, | |
75 | 122.0, 124.0, 126.0, 128.0, 130.0, 132.0, 134.0, 136.0, 138.0, 140.0, | |
76 | 142.0, 144.0, 146.0, 148.0, 150.0, 152.0, 154.0, 156.0, 158.0, 160.0, | |
77 | 162.0, 164.0, 166.0, 168.0, 170.0, 172.0, 174.0, 176.0, 178.0, 180.0, | |
78 | 182.0, 184.0, 186.0, 188.0, 190.0, 192.0, 194.0, 196.0, 198.0, 200.0, | |
79 | 202.0, 204.0, 206.0, 208.0, 210.0, 212.0, 214.0, 216.0, 218.0, 220.0, | |
80 | 222.0, 224.0, 226.0, 228.0, 230.0, 232.0, 234.0, 236.0, 238.0, 240.0, | |
81 | 242.0, 243.0, 244.0, 245.0, 246.0, 247.0, 248.0, 248.5, 249.0, 249.5 } ; | |
82 | ||
83 | const Double_t AliTPCDistortions::fgkIFCRadius= 83.06; // Mean Radius of the Inner Field Cage ( 82.43 min, 83.70 max) (cm) | |
84 | const Double_t AliTPCDistortions::fgkOFCRadius=254.5; // Mean Radius of the Outer Field Cage (252.55 min, 256.45 max) (cm) | |
85 | const Double_t AliTPCDistortions::fgkTPC_Z0=249.7; // Z location of TPC Gated Grid (cm) | |
86 | const Double_t AliTPCDistortions::fgkZOffSet= 0.2; // Offset from CE: calculate all distortions closer to CE as if at this point | |
87 | const Double_t AliTPCDistortions::fgkCathodeV =-100000.0; // Cathode Voltage (volts) | |
88 | const Double_t AliTPCDistortions::fgkGG =-70.0; // Gating Grid voltage (volts) | |
89 | const Double_t AliTPCDistortions::fgkAliceDriftV= 2.73; // Drift Velocity (cm/microSec) Magnitude | |
90 | ||
91 | ||
92 | ||
93 | ||
94 | AliTPCDistortions* AliTPCDistortions::fgInstance = 0; | |
95 | ||
96 | ||
97 | ||
98 | ||
99 | //_ singleton implementation __________________________________________________ | |
100 | AliTPCDistortions* AliTPCDistortions::Instance() | |
101 | { | |
102 | // | |
103 | // Singleton implementation | |
104 | // Returns an instance of this class, it is created if neccessary | |
105 | // | |
106 | if (fgInstance == 0){ | |
107 | fgInstance = new AliTPCDistortions(); | |
108 | //fgInstance->Init(); | |
109 | } | |
110 | return fgInstance; | |
111 | } | |
112 | ||
113 | ||
43e55a68 | 114 | AliTPCDistortions::AliTPCDistortions() |
7f550ed6 | 115 | :fOmegaTau(0.),fT1(1.),fT2(1.),fC0(1.),fC1(0.),fC2(0.), |
43e55a68 | 116 | fBField(0),fXTwist(0.),fYTwist(0.),fIFCShift(0.),fDeltaVGGA(0.),fDeltaVGGC(0.), |
117 | fJLow(0),fKLow(0) | |
118 | { | |
119 | // | |
120 | // default constructor (all corrections are set to have zero effect) | |
121 | // | |
122 | ||
123 | InitIFCShiftDistortion(); | |
124 | InitGGVoltErrorDistortion(); | |
125 | } | |
126 | ||
127 | void AliTPCDistortions::RecalculateCs() { | |
128 | // | |
129 | // helper function to recalculate the coefficients C0-C2 | |
130 | // after a change of omega-tau or the tensor terms | |
131 | // | |
132 | fC0=1. / ( 1. + fT2*fT2*fOmegaTau*fOmegaTau ) ; | |
133 | fC1=fT1*fOmegaTau / ( 1. + fT1*fT1*fOmegaTau*fOmegaTau ) ; | |
134 | fC2=fT2*fT2*fOmegaTau*fOmegaTau / ( 1. + fT2*fT2*fOmegaTau*fOmegaTau ) ; | |
135 | } | |
136 | ||
137 | void AliTPCDistortions::UndoTwistDistortion(const Double_t x[],Double_t xprime[],Int_t roc) | |
138 | { | |
139 | // | |
140 | // Twist distortion | |
141 | // | |
142 | // Remove the effects of a simple "twist" of the TPC in the magnet. If there is | |
143 | // an angle between the E and B fields, there will be a distortion in the recorded | |
144 | // tracks. This routine takes out that distortion. | |
145 | // | |
146 | // The parameters are: | |
147 | // twist[0]: twist angle in x-z | |
148 | // twist[1]: twist angle in y-z | |
149 | // | |
150 | ||
151 | Double_t zdrift ; | |
152 | Int_t sign ; | |
153 | ||
154 | Double_t z = x[2] ; // Creat temporary copy of x[2] | |
155 | ||
156 | if ( roc < 36 ) sign = 1 ; // (TPC End A) | |
157 | else sign = -1 ; // (TPC End C) | |
158 | if ( roc <= 35 && z < fgkZOffSet ) z = fgkZOffSet ; // Protect against discontinuity at CE | |
159 | if ( roc >= 36 && z > -fgkZOffSet ) z = -fgkZOffSet ; // Protect against discontinuity at CE | |
160 | ||
161 | zdrift = sign * ( fgkTPC_Z0 - TMath::Abs(z) ) ; | |
162 | ||
163 | xprime[0] = x[0] - ( fC1 * fYTwist - fC2 * fXTwist ) * zdrift/1000 ; | |
164 | xprime[1] = x[1] - ( -1* fC1 * fXTwist - fC2 * fYTwist ) * zdrift/1000 ; | |
165 | xprime[2] = x[2] ; // Subtract to undo the distortion | |
166 | } | |
167 | ||
168 | ||
169 | void AliTPCDistortions::InitIFCShiftDistortion() { | |
170 | // cout << "AliTPCDistortions::IFCShift Please wait for the tables to fill ... ~5 seconds" << endl ; | |
171 | Int_t nterms = 100 ; | |
172 | Double_t r,z; | |
173 | for ( Int_t i = 0 ; i < kNZ ; ++i ) { | |
174 | z = TMath::Abs( fgkZList[i] ) ; | |
175 | for ( Int_t j = 0 ; j < kNR ; ++j ) { | |
176 | r = fgkRList[j] ; | |
177 | fShiftER[i][j] = 0.0 ; | |
178 | Double_t intz = 0.0 ; | |
179 | for ( Int_t n = 1 ; n < nterms ; ++n ) { | |
180 | Double_t k = (2*n-1) * TMath::Pi() / fgkTPC_Z0 ; | |
181 | Double_t Cn = -4.0 / ( k * fgkTPC_Z0 ) ; | |
182 | Double_t numerator = | |
183 | TMath::BesselK0( k*fgkOFCRadius ) * TMath::BesselI1( k*r ) + | |
184 | TMath::BesselK1( k*r ) * TMath::BesselI0( k*fgkOFCRadius ) ; | |
185 | Double_t denominator = | |
186 | TMath::BesselK0( k*fgkOFCRadius ) * TMath::BesselI0( k*fgkIFCRadius ) - | |
187 | TMath::BesselK0( k*fgkIFCRadius ) * TMath::BesselI0( k*fgkOFCRadius ) ; | |
188 | Double_t zterm = 1 + TMath::Cos( k*z ) ; | |
189 | Double_t qwe = numerator / denominator ; | |
190 | intz += Cn * zterm * qwe ; | |
191 | if ( n>10 && fabs(intz)*1.e-10 > fabs(qwe) ) break; | |
192 | } | |
193 | if ( fgkZList[i] < 0 ) intz = -1 * intz ; // Force AntiSymmetry of solutions in Z | |
194 | fShiftER[i][j] = intz ; | |
195 | } | |
196 | } | |
197 | } | |
198 | ||
199 | void AliTPCDistortions::UndoIFCShiftDistortion(const Double_t x[],Double_t xprime[],Int_t roc) | |
200 | { | |
201 | // | |
202 | // IFC Shift Distortion | |
203 | // | |
204 | // The Inner field cage of the TPC may not be perfectly aligned with the outer field cage | |
205 | // of the TPC. They can be shifted along the Z axis by up to a 1 mm. This causes a tilting | |
206 | // of the equi-potential lines inside the TPC and therefore a DCA error at the vertex. | |
207 | // The distortion is anti-symmetric in Z. | |
208 | // | |
209 | ||
210 | Int_t order = 1 ; // Linear interpolation = 1, Quadratic = 2 | |
211 | ||
212 | Double_t intEr, intEphi ; | |
213 | Double_t r, phi, z ; | |
214 | Int_t sign ; | |
215 | ||
216 | r = TMath::Sqrt( x[0]*x[0] + x[1]*x[1] ) ; | |
217 | phi = TMath::ATan2(x[1],x[0]) ; | |
218 | if ( phi < 0 ) phi += TMath::TwoPi() ; // Table uses phi from 0 to 2*Pi | |
219 | z = x[2] ; // Create temporary copy of x[2] | |
220 | ||
221 | if ( roc < 36 ) sign = 1 ; // (TPC End A) | |
222 | else sign = -1 ; // (TPC End C) | |
223 | if ( roc <= 35 && z < fgkZOffSet ) z = fgkZOffSet ; // Protect against discontinuity at CE | |
224 | if ( roc >= 36 && z > -fgkZOffSet ) z = -fgkZOffSet ; // Protect against discontinuity at CE | |
225 | ||
226 | Interpolate2DEdistortion( order, r, z, fShiftER, intEr ) ; | |
227 | intEphi = 0.0 ; // Efield is symmetric in phi | |
228 | ||
229 | // Subtract to Undo the distortions | |
230 | if ( r > 0.0 ) | |
231 | { | |
232 | phi = phi - fIFCShift*( fC0*intEphi - fC1*intEr ) / r ; | |
233 | r = r - fIFCShift*( fC0*intEr + fC1*intEphi ) ; | |
234 | } | |
235 | ||
236 | xprime[0] = r * TMath::Cos(phi) ; | |
237 | xprime[1] = r * TMath::Sin(phi) ; | |
238 | xprime[2] = x[2] ; | |
239 | ||
240 | } | |
241 | ||
242 | void AliTPCDistortions::InitGGVoltErrorDistortion() { | |
243 | // cout << "AliTPCDistortions::UndoGG VE Please wait for the tables to fill ... ~5 seconds" << endl ; | |
244 | Double_t r,z; | |
245 | Int_t nterms = 100 ; | |
246 | for ( Int_t i = 0 ; i < kNZ ; ++i ) { | |
247 | z = fgkZList[i] ; | |
248 | for ( Int_t j = 0 ; j < kNR ; ++j ) { | |
249 | r = fgkRList[j] ; | |
250 | fGGVoltErrorER[i][j] = 0.0 ; | |
251 | Double_t intz = 0.0 ; | |
252 | for ( Int_t n = 1 ; n < nterms ; ++n ) { | |
253 | Double_t k = n * TMath::Pi() / fgkTPC_Z0 ; | |
254 | Double_t ein = 0 ; // Error potential on the IFC | |
255 | Double_t eout = 0 ; // Error potential on the OFC | |
256 | if ( z < 0 ) { | |
257 | ein = 2.0 * ((z-fgkTPC_Z0)/fgkTPC_Z0) / ( k * (fgkCathodeV - fgkGG) ) ; | |
258 | eout = 2.0 * ((z-fgkTPC_Z0)/fgkTPC_Z0) / ( k * (fgkCathodeV - fgkGG) ) ; | |
259 | } | |
260 | if ( z == 0 ) continue ; | |
261 | if ( z > 0 ) { | |
262 | ein = 2.0 * ((fgkTPC_Z0-z)/fgkTPC_Z0) / ( k * (fgkCathodeV - fgkGG) ) ; | |
263 | eout = 2.0 * ((fgkTPC_Z0-z)/fgkTPC_Z0) / ( k * (fgkCathodeV - fgkGG) ) ; | |
264 | } | |
265 | Double_t an = ein * TMath::BesselK0( k*fgkOFCRadius ) - eout * TMath::BesselK0( k*fgkIFCRadius ) ; | |
266 | Double_t bn = eout * TMath::BesselI0( k*fgkIFCRadius ) - ein * TMath::BesselI0( k*fgkOFCRadius ) ; | |
267 | Double_t numerator = | |
268 | an * TMath::BesselI1( k*r ) - bn * TMath::BesselK1( k*r ) ; | |
269 | Double_t denominator = | |
270 | TMath::BesselK0( k*fgkOFCRadius ) * TMath::BesselI0( k*fgkIFCRadius ) - | |
271 | TMath::BesselK0( k*fgkIFCRadius ) * TMath::BesselI0( k*fgkOFCRadius ) ; | |
272 | Double_t zterm = TMath::Cos( k*(fgkTPC_Z0-TMath::Abs(z)) ) - 1 ; | |
273 | intz += zterm * numerator / denominator ; | |
274 | // Assume series converges, break if small terms | |
275 | if ( n>10 && fabs(intz)*1.e-10 > fabs(numerator/denominator) ) break; | |
276 | } | |
277 | fGGVoltErrorER[i][j] = intz ; | |
278 | } | |
279 | } | |
280 | } | |
281 | ||
282 | ||
283 | ||
284 | void AliTPCDistortions::UndoGGVoltErrorDistortion(const Double_t x[], Double_t xprime[],Int_t roc) | |
285 | { | |
286 | // | |
287 | // Gated Grid Voltage Error | |
288 | // | |
289 | // Calculate the effect of having an incorrect voltage on the A or C end plate Gated Grids. | |
290 | // | |
291 | // Electrostatic Equations from StarNote SN0253 by Howard Wieman. | |
292 | // Note that we use a funny coordinate system where Z==0 at the GG. | |
293 | // | |
294 | ||
295 | Int_t order = 1 ; // Linear interpolation = 1, Quadratic = 2 | |
296 | ||
297 | Double_t intEr, intEphi ; | |
298 | Double_t r, phi, z ; | |
299 | Int_t sign ; | |
300 | ||
301 | Double_t deltaVGG; | |
302 | ||
303 | r = TMath::Sqrt( x[0]*x[0] + x[1]*x[1] ) ; | |
304 | phi = TMath::ATan2(x[1],x[0]) ; | |
305 | if ( phi < 0 ) phi += TMath::TwoPi() ; // Table uses phi from 0 to 2*Pi | |
306 | z = x[2] ; | |
307 | ||
308 | if ( roc < 36 ) sign = 1 ; // (TPC End A) | |
309 | else sign = -1 ; // (TPC End C) | |
310 | if ( roc <= 35 && z < fgkZOffSet ) z = fgkZOffSet ; // Protect against discontinuity at CE | |
311 | if ( roc >= 36 && z > -fgkZOffSet ) z = -fgkZOffSet ; // Protect against discontinuity at CE | |
312 | ||
313 | if ( roc < 36 ) deltaVGG = fDeltaVGGA ; // (TPC End A) | |
314 | else deltaVGG = fDeltaVGGC ; // (TPC End C) | |
315 | ||
316 | Interpolate2DEdistortion( order, r, z, fGGVoltErrorER, intEr ) ; | |
317 | intEphi = 0.0 ; // Efield is symmetric in phi | |
318 | ||
319 | // Subtract to Undo the distortions | |
320 | if ( r > 0.0 ) | |
321 | { | |
322 | phi = phi - deltaVGG*( fC0*intEphi - fC1*intEr ) / r ; | |
323 | r = r - deltaVGG*( fC0*intEr + fC1*intEphi ) ; | |
324 | } | |
325 | ||
326 | xprime[0] = r * TMath::Cos(phi) ; | |
327 | xprime[1] = r * TMath::Sin(phi) ; | |
328 | xprime[2] = x[2] ; | |
329 | ||
330 | } | |
331 | ||
332 | void AliTPCDistortions::UndoExBShapeDistortion( const Double_t x[], Double_t xprime[], Int_t roc) | |
333 | { | |
334 | // | |
335 | // Example of using Ruben's integrals to do ExB Shape distortions | |
336 | // | |
337 | // | |
338 | ||
339 | xprime[0] = x[0] ; xprime[1] = x[1] ; xprime[2] = x[2] ; | |
340 | ||
341 | if (!fBField) return; | |
342 | ||
343 | Double_t intBStart[3], intBEnd[3] , xStart[3], xEnd[3] ; | |
344 | Double_t intBxOverBz, intByOverBz; | |
345 | //Double_t denominator ; | |
346 | Int_t sign; | |
347 | if ( roc < 36 ) sign = 1 ; // (TPC End A) | |
348 | else sign = -1 ; // (TPC End C) | |
349 | ||
350 | xStart[0] = x[0] ; xStart[1] = x[1] ; xStart[2] = x[2] ; | |
351 | xEnd [0] = x[0] ; xEnd [1] = x[1] ; xEnd [2] = sign * fgkTPC_Z0 ; | |
352 | ||
353 | fBField -> GetTPCInt(xStart, intBStart) ; | |
354 | fBField -> GetTPCInt(xEnd , intBEnd ) ; | |
355 | ||
356 | if ( TMath::Abs(intBStart[2] - intBEnd[2]) < 0.1 ) return ; // Protect against divide by zero, below | |
357 | ||
358 | intBxOverBz = (intBStart[0] - intBEnd[0]) ; | |
359 | intByOverBz = (intBStart[1] - intBEnd[1]) ; | |
360 | ||
361 | xprime[0] += ( fC2*intBxOverBz - fC1*intByOverBz ) ; | |
362 | xprime[1] += ( fC2*intByOverBz + fC1*intBxOverBz ) ; | |
363 | xprime[2] += 0.0 ; | |
364 | ||
365 | } | |
366 | ||
367 | ||
368 | void AliTPCDistortions::Interpolate2DEdistortion( const Int_t order, const Double_t r, const Double_t z, | |
369 | const Double_t er[kNZ][kNR], Double_t &er_value ) | |
370 | { | |
371 | // | |
372 | // Interpolate table - 2D interpolation | |
373 | // | |
374 | Double_t save_er[10] ; | |
375 | ||
376 | Search( kNZ, fgkZList, z, fJLow ) ; | |
377 | Search( kNR, fgkRList, r, fKLow ) ; | |
378 | if ( fJLow < 0 ) fJLow = 0 ; // check if out of range | |
379 | if ( fKLow < 0 ) fKLow = 0 ; | |
380 | if ( fJLow + order >= kNZ - 1 ) fJLow = kNZ - 1 - order ; | |
381 | if ( fKLow + order >= kNR - 1 ) fKLow = kNR - 1 - order ; | |
382 | ||
383 | for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) | |
384 | { | |
385 | save_er[j-fJLow] = Interpolate( &fgkRList[fKLow], &er[j][fKLow], order, r ) ; | |
386 | } | |
387 | er_value = Interpolate( &fgkZList[fJLow], save_er, order, z ) ; | |
388 | ||
389 | } | |
390 | ||
391 | ||
392 | Double_t AliTPCDistortions::Interpolate( const Double_t xArray[], const Double_t yArray[], | |
393 | const Int_t order, const Double_t x ) | |
394 | { | |
395 | // | |
396 | // Interpolate function Y(x) using linear (order=1) or quadratic (order=2) interpolation. | |
397 | // | |
398 | ||
399 | Double_t y ; | |
400 | ||
401 | if ( order == 2 ) // Quadratic Interpolation = 2 | |
402 | ||
403 | { | |
404 | y = (x-xArray[1]) * (x-xArray[2]) * yArray[0] / ( (xArray[0]-xArray[1]) * (xArray[0]-xArray[2]) ) ; | |
405 | y += (x-xArray[2]) * (x-xArray[0]) * yArray[1] / ( (xArray[1]-xArray[2]) * (xArray[1]-xArray[0]) ) ; | |
406 | y += (x-xArray[0]) * (x-xArray[1]) * yArray[2] / ( (xArray[2]-xArray[0]) * (xArray[2]-xArray[1]) ) ; | |
407 | ||
408 | } | |
409 | ||
410 | else // Linear Interpolation = 1 | |
411 | ||
412 | { | |
413 | y = yArray[0] + ( yArray[1]-yArray[0] ) * ( x-xArray[0] ) / ( xArray[1] - xArray[0] ) ; | |
414 | } | |
415 | ||
416 | return (y) ; | |
417 | ||
418 | } | |
419 | ||
420 | ||
421 | void AliTPCDistortions::Search( const Int_t n, const Double_t xArray[], const Double_t x, Int_t &low ) | |
422 | { | |
423 | // | |
424 | // Search an ordered table by starting at the most recently used point | |
425 | // | |
426 | ||
427 | Long_t middle, high ; | |
428 | Int_t ascend = 0, increment = 1 ; | |
429 | ||
430 | if ( xArray[n-1] >= xArray[0] ) ascend = 1 ; // Ascending ordered table if true | |
431 | ||
432 | if ( low < 0 || low > n-1 ) { low = -1 ; high = n ; } | |
433 | ||
434 | else // Ordered Search phase | |
435 | { | |
436 | if ( (Int_t)( x >= xArray[low] ) == ascend ) | |
437 | { | |
438 | if ( low == n-1 ) return ; | |
439 | high = low + 1 ; | |
440 | while ( (Int_t)( x >= xArray[high] ) == ascend ) | |
441 | { | |
442 | low = high ; | |
443 | increment *= 2 ; | |
444 | high = low + increment ; | |
445 | if ( high > n-1 ) { high = n ; break ; } | |
446 | } | |
447 | } | |
448 | else | |
449 | { | |
450 | if ( low == 0 ) { low = -1 ; return ; } | |
451 | high = low - 1 ; | |
452 | while ( (Int_t)( x < xArray[low] ) == ascend ) | |
453 | { | |
454 | high = low ; | |
455 | increment *= 2 ; | |
456 | if ( increment >= high ) { low = -1 ; break ; } | |
457 | else low = high - increment ; | |
458 | } | |
459 | } | |
460 | } | |
461 | ||
462 | while ( (high-low) != 1 ) // Binary Search Phase | |
463 | { | |
464 | middle = ( high + low ) / 2 ; | |
465 | if ( (Int_t)( x >= xArray[middle] ) == ascend ) | |
466 | low = middle ; | |
467 | else | |
468 | high = middle ; | |
469 | } | |
470 | ||
471 | if ( x == xArray[n-1] ) low = n-2 ; | |
472 | if ( x == xArray[0] ) low = 0 ; | |
473 | ||
474 | } | |
475 | ||
476 | ClassImp(AliTPCDistortions); |