Line data Source code
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 : #include "AliMagWrapCheb.h"
17 : #include "AliLog.h"
18 : #include <TSystem.h>
19 : #include <TArrayF.h>
20 : #include <TArrayI.h>
21 :
22 176 : ClassImp(AliMagWrapCheb)
23 :
24 : //__________________________________________________________________________________________
25 5 : AliMagWrapCheb::AliMagWrapCheb() :
26 5 : fNParamsSol(0),fNZSegSol(0),fNPSegSol(0),fNRSegSol(0),
27 5 : fSegZSol(0),fSegPSol(0),fSegRSol(0),
28 5 : fBegSegPSol(0),fNSegPSol(0),fBegSegRSol(0),fNSegRSol(0),fSegIDSol(0),fMinZSol(1.e6),fMaxZSol(-1.e6),fParamsSol(0),fMaxRSol(0),
29 : //
30 5 : fNParamsTPC(0),fNZSegTPC(0),fNPSegTPC(0),fNRSegTPC(0),
31 5 : fSegZTPC(0),fSegPTPC(0),fSegRTPC(0),
32 5 : fBegSegPTPC(0),fNSegPTPC(0),fBegSegRTPC(0),fNSegRTPC(0),fSegIDTPC(0),fMinZTPC(1.e6),fMaxZTPC(-1.e6),fParamsTPC(0),fMaxRTPC(0),
33 : //
34 5 : fNParamsTPCRat(0),fNZSegTPCRat(0),fNPSegTPCRat(0),fNRSegTPCRat(0),
35 5 : fSegZTPCRat(0),fSegPTPCRat(0),fSegRTPCRat(0),
36 5 : fBegSegPTPCRat(0),fNSegPTPCRat(0),fBegSegRTPCRat(0),fNSegRTPCRat(0),fSegIDTPCRat(0),fMinZTPCRat(1.e6),fMaxZTPCRat(-1.e6),fParamsTPCRat(0),fMaxRTPCRat(0),
37 : //
38 5 : fNParamsDip(0),fNZSegDip(0),fNYSegDip(0),fNXSegDip(0),
39 5 : fSegZDip(0),fSegYDip(0),fSegXDip(0),
40 5 : fBegSegYDip(0),fNSegYDip(0),fBegSegXDip(0),fNSegXDip(0),fSegIDDip(0),fMinZDip(1.e6),fMaxZDip(-1.e6),fParamsDip(0)
41 : //
42 : #ifdef _MAGCHEB_CACHE_
43 5 : ,fCacheSol(0),fCacheDip(0),fCacheTPCInt(0),fCacheTPCRat(0)
44 : #endif
45 : //
46 25 : {
47 : // default constructor
48 10 : }
49 :
50 : //__________________________________________________________________________________________
51 : AliMagWrapCheb::AliMagWrapCheb(const AliMagWrapCheb& src) :
52 0 : TNamed(src),
53 0 : fNParamsSol(0),fNZSegSol(0),fNPSegSol(0),fNRSegSol(0),
54 0 : fSegZSol(0),fSegPSol(0),fSegRSol(0),
55 0 : fBegSegPSol(0),fNSegPSol(0),fBegSegRSol(0),fNSegRSol(0),fSegIDSol(0),fMinZSol(1.e6),fMaxZSol(-1.e6),fParamsSol(0),fMaxRSol(0),
56 : //
57 0 : fNParamsTPC(0),fNZSegTPC(0),fNPSegTPC(0),fNRSegTPC(0),
58 0 : fSegZTPC(0),fSegPTPC(0),fSegRTPC(0),
59 0 : fBegSegPTPC(0),fNSegPTPC(0),fBegSegRTPC(0),fNSegRTPC(0),fSegIDTPC(0),fMinZTPC(1.e6),fMaxZTPC(-1.e6),fParamsTPC(0),fMaxRTPC(0),
60 : //
61 0 : fNParamsTPCRat(0),fNZSegTPCRat(0),fNPSegTPCRat(0),fNRSegTPCRat(0),
62 0 : fSegZTPCRat(0),fSegPTPCRat(0),fSegRTPCRat(0),
63 0 : fBegSegPTPCRat(0),fNSegPTPCRat(0),fBegSegRTPCRat(0),fNSegRTPCRat(0),fSegIDTPCRat(0),fMinZTPCRat(1.e6),fMaxZTPCRat(-1.e6),fParamsTPCRat(0),fMaxRTPCRat(0),
64 : //
65 0 : fNParamsDip(0),fNZSegDip(0),fNYSegDip(0),fNXSegDip(0),
66 0 : fSegZDip(0),fSegYDip(0),fSegXDip(0),
67 0 : fBegSegYDip(0),fNSegYDip(0),fBegSegXDip(0),fNSegXDip(0),fSegIDDip(0),fMinZDip(1.e6),fMaxZDip(-1.e6),fParamsDip(0)
68 : //
69 : #ifdef _MAGCHEB_CACHE_
70 0 : ,fCacheSol(0),fCacheDip(0),fCacheTPCInt(0),fCacheTPCRat(0)
71 : #endif
72 0 : {
73 : // copy constructor
74 0 : CopyFrom(src);
75 : //
76 0 : }
77 :
78 : //__________________________________________________________________________________________
79 : void AliMagWrapCheb::CopyFrom(const AliMagWrapCheb& src)
80 : {
81 : // copy method
82 0 : Clear();
83 0 : SetName(src.GetName());
84 0 : SetTitle(src.GetTitle());
85 : //
86 0 : fNParamsSol = src.fNParamsSol;
87 0 : fNZSegSol = src.fNZSegSol;
88 0 : fNPSegSol = src.fNPSegSol;
89 0 : fNRSegSol = src.fNRSegSol;
90 0 : fMinZSol = src.fMinZSol;
91 0 : fMaxZSol = src.fMaxZSol;
92 0 : fMaxRSol = src.fMaxRSol;
93 0 : if (src.fNParamsSol) {
94 0 : memcpy(fSegZSol = new Float_t[fNZSegSol], src.fSegZSol, sizeof(Float_t)*fNZSegSol);
95 0 : memcpy(fSegPSol = new Float_t[fNPSegSol], src.fSegPSol, sizeof(Float_t)*fNPSegSol);
96 0 : memcpy(fSegRSol = new Float_t[fNRSegSol], src.fSegRSol, sizeof(Float_t)*fNRSegSol);
97 0 : memcpy(fBegSegPSol= new Int_t[fNZSegSol], src.fBegSegPSol, sizeof(Int_t)*fNZSegSol);
98 0 : memcpy(fNSegPSol = new Int_t[fNZSegSol], src.fNSegPSol, sizeof(Int_t)*fNZSegSol);
99 0 : memcpy(fBegSegRSol= new Int_t[fNPSegSol], src.fBegSegRSol, sizeof(Int_t)*fNPSegSol);
100 0 : memcpy(fNSegRSol = new Int_t[fNPSegSol], src.fNSegRSol, sizeof(Int_t)*fNPSegSol);
101 0 : memcpy(fSegIDSol = new Int_t[fNRSegSol], src.fSegIDSol, sizeof(Int_t)*fNRSegSol);
102 0 : fParamsSol = new TObjArray(fNParamsSol);
103 0 : for (int i=0;i<fNParamsSol;i++) fParamsSol->AddAtAndExpand(new AliCheb3D(*src.GetParamSol(i)),i);
104 0 : }
105 : //
106 0 : fNParamsTPC = src.fNParamsTPC;
107 0 : fNZSegTPC = src.fNZSegTPC;
108 0 : fNPSegTPC = src.fNPSegTPC;
109 0 : fNRSegTPC = src.fNRSegTPC;
110 0 : fMinZTPC = src.fMinZTPC;
111 0 : fMaxZTPC = src.fMaxZTPC;
112 0 : fMaxRTPC = src.fMaxRTPC;
113 0 : if (src.fNParamsTPC) {
114 0 : memcpy(fSegZTPC = new Float_t[fNZSegTPC], src.fSegZTPC, sizeof(Float_t)*fNZSegTPC);
115 0 : memcpy(fSegPTPC = new Float_t[fNPSegTPC], src.fSegPTPC, sizeof(Float_t)*fNPSegTPC);
116 0 : memcpy(fSegRTPC = new Float_t[fNRSegTPC], src.fSegRTPC, sizeof(Float_t)*fNRSegTPC);
117 0 : memcpy(fBegSegPTPC= new Int_t[fNZSegTPC], src.fBegSegPTPC, sizeof(Int_t)*fNZSegTPC);
118 0 : memcpy(fNSegPTPC = new Int_t[fNZSegTPC], src.fNSegPTPC, sizeof(Int_t)*fNZSegTPC);
119 0 : memcpy(fBegSegRTPC= new Int_t[fNPSegTPC], src.fBegSegRTPC, sizeof(Int_t)*fNPSegTPC);
120 0 : memcpy(fNSegRTPC = new Int_t[fNPSegTPC], src.fNSegRTPC, sizeof(Int_t)*fNPSegTPC);
121 0 : memcpy(fSegIDTPC = new Int_t[fNRSegTPC], src.fSegIDTPC, sizeof(Int_t)*fNRSegTPC);
122 0 : fParamsTPC = new TObjArray(fNParamsTPC);
123 0 : for (int i=0;i<fNParamsTPC;i++) fParamsTPC->AddAtAndExpand(new AliCheb3D(*src.GetParamTPCInt(i)),i);
124 0 : }
125 : //
126 0 : fNParamsTPCRat = src.fNParamsTPCRat;
127 0 : fNZSegTPCRat = src.fNZSegTPCRat;
128 0 : fNPSegTPCRat = src.fNPSegTPCRat;
129 0 : fNRSegTPCRat = src.fNRSegTPCRat;
130 0 : fMinZTPCRat = src.fMinZTPCRat;
131 0 : fMaxZTPCRat = src.fMaxZTPCRat;
132 0 : fMaxRTPCRat = src.fMaxRTPCRat;
133 0 : if (src.fNParamsTPCRat) {
134 0 : memcpy(fSegZTPCRat = new Float_t[fNZSegTPCRat], src.fSegZTPCRat, sizeof(Float_t)*fNZSegTPCRat);
135 0 : memcpy(fSegPTPCRat = new Float_t[fNPSegTPCRat], src.fSegPTPCRat, sizeof(Float_t)*fNPSegTPCRat);
136 0 : memcpy(fSegRTPCRat = new Float_t[fNRSegTPCRat], src.fSegRTPCRat, sizeof(Float_t)*fNRSegTPCRat);
137 0 : memcpy(fBegSegPTPCRat= new Int_t[fNZSegTPCRat], src.fBegSegPTPCRat, sizeof(Int_t)*fNZSegTPCRat);
138 0 : memcpy(fNSegPTPCRat = new Int_t[fNZSegTPCRat], src.fNSegPTPCRat, sizeof(Int_t)*fNZSegTPCRat);
139 0 : memcpy(fBegSegRTPCRat= new Int_t[fNPSegTPCRat], src.fBegSegRTPCRat, sizeof(Int_t)*fNPSegTPCRat);
140 0 : memcpy(fNSegRTPCRat = new Int_t[fNPSegTPCRat], src.fNSegRTPCRat, sizeof(Int_t)*fNPSegTPCRat);
141 0 : memcpy(fSegIDTPCRat = new Int_t[fNRSegTPCRat], src.fSegIDTPCRat, sizeof(Int_t)*fNRSegTPCRat);
142 0 : fParamsTPCRat = new TObjArray(fNParamsTPCRat);
143 0 : for (int i=0;i<fNParamsTPCRat;i++) fParamsTPCRat->AddAtAndExpand(new AliCheb3D(*src.GetParamTPCRatInt(i)),i);
144 0 : }
145 : //
146 0 : fNParamsDip = src.fNParamsDip;
147 0 : fNZSegDip = src.fNZSegDip;
148 0 : fNYSegDip = src.fNYSegDip;
149 0 : fNXSegDip = src.fNXSegDip;
150 0 : fMinZDip = src.fMinZDip;
151 0 : fMaxZDip = src.fMaxZDip;
152 0 : if (src.fNParamsDip) {
153 0 : memcpy(fSegZDip = new Float_t[fNZSegDip], src.fSegZDip, sizeof(Float_t)*fNZSegDip);
154 0 : memcpy(fSegYDip = new Float_t[fNYSegDip], src.fSegYDip, sizeof(Float_t)*fNYSegDip);
155 0 : memcpy(fSegXDip = new Float_t[fNXSegDip], src.fSegXDip, sizeof(Float_t)*fNXSegDip);
156 0 : memcpy(fBegSegYDip= new Int_t[fNZSegDip], src.fBegSegYDip, sizeof(Int_t)*fNZSegDip);
157 0 : memcpy(fNSegYDip = new Int_t[fNZSegDip], src.fNSegYDip, sizeof(Int_t)*fNZSegDip);
158 0 : memcpy(fBegSegXDip= new Int_t[fNYSegDip], src.fBegSegXDip, sizeof(Int_t)*fNYSegDip);
159 0 : memcpy(fNSegXDip = new Int_t[fNYSegDip], src.fNSegXDip, sizeof(Int_t)*fNYSegDip);
160 0 : memcpy(fSegIDDip = new Int_t[fNXSegDip], src.fSegIDDip, sizeof(Int_t)*fNXSegDip);
161 0 : fParamsDip = new TObjArray(fNParamsDip);
162 0 : for (int i=0;i<fNParamsDip;i++) fParamsDip->AddAtAndExpand(new AliCheb3D(*src.GetParamDip(i)),i);
163 0 : }
164 : //
165 0 : }
166 :
167 : //__________________________________________________________________________________________
168 : AliMagWrapCheb& AliMagWrapCheb::operator=(const AliMagWrapCheb& rhs)
169 : {
170 : // assignment
171 0 : if (this != &rhs) {
172 0 : Clear();
173 0 : CopyFrom(rhs);
174 0 : }
175 0 : return *this;
176 : //
177 : }
178 :
179 : //__________________________________________________________________________________________
180 : void AliMagWrapCheb::Clear(const Option_t *)
181 : {
182 : // clear all dynamic parts
183 10 : if (fNParamsSol) {
184 5 : fParamsSol->SetOwner(kTRUE);
185 15 : delete fParamsSol; fParamsSol = 0;
186 15 : delete[] fSegZSol; fSegZSol = 0;
187 15 : delete[] fSegPSol; fSegPSol = 0;
188 15 : delete[] fSegRSol; fSegRSol = 0;
189 15 : delete[] fBegSegPSol; fBegSegPSol = 0;
190 15 : delete[] fNSegPSol; fNSegPSol = 0;
191 15 : delete[] fBegSegRSol; fBegSegRSol = 0;
192 15 : delete[] fNSegRSol; fNSegRSol = 0;
193 15 : delete[] fSegIDSol; fSegIDSol = 0;
194 5 : }
195 5 : fNParamsSol = fNZSegSol = fNPSegSol = fNRSegSol = 0;
196 5 : fMinZSol = 1e6;
197 5 : fMaxZSol = -1e6;
198 5 : fMaxRSol = 0;
199 : //
200 5 : if (fNParamsTPC) {
201 5 : fParamsTPC->SetOwner(kTRUE);
202 15 : delete fParamsTPC; fParamsTPC = 0;
203 15 : delete[] fSegZTPC; fSegZTPC = 0;
204 15 : delete[] fSegPTPC; fSegPTPC = 0;
205 15 : delete[] fSegRTPC; fSegRTPC = 0;
206 15 : delete[] fBegSegPTPC; fBegSegPTPC = 0;
207 15 : delete[] fNSegPTPC; fNSegPTPC = 0;
208 15 : delete[] fBegSegRTPC; fBegSegRTPC = 0;
209 15 : delete[] fNSegRTPC; fNSegRTPC = 0;
210 15 : delete[] fSegIDTPC; fSegIDTPC = 0;
211 5 : }
212 5 : fNParamsTPC = fNZSegTPC = fNPSegTPC = fNRSegTPC = 0;
213 5 : fMinZTPC = 1e6;
214 5 : fMaxZTPC = -1e6;
215 5 : fMaxRTPC = 0;
216 : //
217 5 : if (fNParamsTPCRat) {
218 5 : fParamsTPCRat->SetOwner(kTRUE);
219 15 : delete fParamsTPCRat; fParamsTPCRat = 0;
220 15 : delete[] fSegZTPCRat; fSegZTPCRat = 0;
221 15 : delete[] fSegPTPCRat; fSegPTPCRat = 0;
222 15 : delete[] fSegRTPCRat; fSegRTPCRat = 0;
223 15 : delete[] fBegSegPTPCRat; fBegSegPTPCRat = 0;
224 15 : delete[] fNSegPTPCRat; fNSegPTPCRat = 0;
225 15 : delete[] fBegSegRTPCRat; fBegSegRTPCRat = 0;
226 15 : delete[] fNSegRTPCRat; fNSegRTPCRat = 0;
227 15 : delete[] fSegIDTPCRat; fSegIDTPCRat = 0;
228 5 : }
229 5 : fNParamsTPCRat = fNZSegTPCRat = fNPSegTPCRat = fNRSegTPCRat = 0;
230 5 : fMinZTPCRat = 1e6;
231 5 : fMaxZTPCRat = -1e6;
232 5 : fMaxRTPCRat = 0;
233 : //
234 5 : if (fNParamsDip) {
235 5 : fParamsDip->SetOwner(kTRUE);
236 15 : delete fParamsDip; fParamsDip = 0;
237 15 : delete[] fSegZDip; fSegZDip = 0;
238 15 : delete[] fSegYDip; fSegYDip = 0;
239 15 : delete[] fSegXDip; fSegXDip = 0;
240 15 : delete[] fBegSegYDip; fBegSegYDip = 0;
241 15 : delete[] fNSegYDip; fNSegYDip = 0;
242 15 : delete[] fBegSegXDip; fBegSegXDip = 0;
243 15 : delete[] fNSegXDip; fNSegXDip = 0;
244 15 : delete[] fSegIDDip; fSegIDDip = 0;
245 5 : }
246 5 : fNParamsDip = fNZSegDip = fNYSegDip = fNXSegDip = 0;
247 5 : fMinZDip = 1e6;
248 5 : fMaxZDip = -1e6;
249 : //
250 : #ifdef _MAGCHEB_CACHE_
251 5 : fCacheSol = 0;
252 5 : fCacheDip = 0;
253 5 : fCacheTPCInt = 0;
254 5 : fCacheTPCRat = 0;
255 : #endif
256 : //
257 5 : }
258 :
259 : //__________________________________________________________________________________________
260 : void AliMagWrapCheb::Field(const Double_t *xyz, Double_t *b) const
261 : {
262 : // compute field in cartesian coordinates. If point is outside of the parameterized region
263 : // get it at closest valid point
264 5752352 : Double_t rphiz[3];
265 : //
266 : #ifndef _BRING_TO_BOUNDARY_ // exact matching to fitted volume is requested
267 2876176 : b[0] = b[1] = b[2] = 0;
268 : #endif
269 : //
270 2876176 : if (xyz[2]>fMinZSol) {
271 2783475 : CartToCyl(xyz,rphiz);
272 : //
273 : #ifdef _MAGCHEB_CACHE_
274 5566950 : if (fCacheSol && fCacheSol->IsInside(rphiz))
275 2671237 : fCacheSol->Eval(rphiz,b);
276 : else
277 : #endif //_MAGCHEB_CACHE_
278 112238 : FieldCylSol(rphiz,b);
279 : // convert field to cartesian system
280 2783475 : CylToCartCylB(rphiz, b,b);
281 2783475 : return;
282 : }
283 : //
284 : #ifdef _MAGCHEB_CACHE_
285 185399 : if (fCacheDip && fCacheDip->IsInside(xyz)) {
286 71415 : fCacheDip->Eval(xyz,b); // check the cache first
287 71415 : return;
288 : }
289 : #else //_MAGCHEB_CACHE_
290 : AliCheb3D* fCacheDip = 0;
291 : #endif //_MAGCHEB_CACHE_
292 21286 : int iddip = FindDipSegment(xyz);
293 21286 : if (iddip>=0) {
294 21286 : fCacheDip = GetParamDip(iddip);
295 : //
296 : #ifndef _BRING_TO_BOUNDARY_
297 21995 : if (!fCacheDip->IsInside(xyz)) return;
298 : #endif //_BRING_TO_BOUNDARY_
299 : //
300 20577 : fCacheDip->Eval(xyz,b);
301 20577 : }
302 : //
303 2896753 : }
304 :
305 : //__________________________________________________________________________________________
306 : Double_t AliMagWrapCheb::GetBz(const Double_t *xyz) const
307 : {
308 : // compute Bz for the point in cartesian coordinates. If point is outside of the parameterized region
309 : // get it at closest valid point
310 594848 : Double_t rphiz[3];
311 : //
312 297424 : if (xyz[2]>fMinZSol) {
313 : //
314 297424 : CartToCyl(xyz,rphiz);
315 : //
316 : #ifdef _MAGCHEB_CACHE_
317 876200 : if (fCacheSol && fCacheSol->IsInside(rphiz)) return fCacheSol->Eval(rphiz,2);
318 : #endif //_MAGCHEB_CACHE_
319 16067 : return FieldCylSolBz(rphiz);
320 : }
321 : //
322 : #ifdef _MAGCHEB_CACHE_
323 0 : if (fCacheDip && fCacheDip->IsInside(xyz)) return fCacheDip->Eval(xyz,2); // check the cache first
324 : //
325 : #else //_MAGCHEB_CACHE_
326 : AliCheb3D* fCacheDip = 0;
327 : #endif //_MAGCHEB_CACHE_
328 : //
329 0 : int iddip = FindDipSegment(xyz);
330 0 : if (iddip>=0) {
331 0 : fCacheDip = GetParamDip(iddip);
332 : //
333 : #ifndef _BRING_TO_BOUNDARY_
334 0 : if (!fCacheDip->IsInside(xyz)) return 0.;
335 : #endif // _BRING_TO_BOUNDARY_
336 : //
337 0 : return fCacheDip->Eval(xyz,2);
338 : //
339 : }
340 : //
341 0 : return 0;
342 : //
343 297424 : }
344 :
345 :
346 : //__________________________________________________________________________________________
347 : void AliMagWrapCheb::Print(Option_t *) const
348 : {
349 : // print info
350 0 : printf("Alice magnetic field parameterized by Chebyshev polynomials\n");
351 0 : printf("Segmentation for Solenoid (%+.2f<Z<%+.2f cm | R<%.2f cm)\n",fMinZSol,fMaxZSol,fMaxRSol);
352 : //
353 0 : if (fParamsSol) {
354 0 : for (int i=0;i<fNParamsSol;i++) {
355 0 : printf("SOL%4d ",i);
356 0 : GetParamSol(i)->Print();
357 : }
358 0 : }
359 : //
360 0 : printf("Segmentation for TPC field integral (%+.2f<Z<%+.2f cm | R<%.2f cm)\n",fMinZTPC,fMaxZTPC,fMaxRTPC);
361 : //
362 0 : if (fParamsTPC) {
363 0 : for (int i=0;i<fNParamsTPC;i++) {
364 0 : printf("TPC%4d ",i);
365 0 : GetParamTPCInt(i)->Print();
366 : }
367 0 : }
368 : //
369 0 : printf("Segmentation for TPC field ratios integral (%+.2f<Z<%+.2f cm | R<%.2f cm)\n",fMinZTPCRat,fMaxZTPCRat,fMaxRTPCRat);
370 : //
371 0 : if (fParamsTPCRat) {
372 0 : for (int i=0;i<fNParamsTPCRat;i++) {
373 0 : printf("TPCRat%4d ",i);
374 0 : GetParamTPCRatInt(i)->Print();
375 : }
376 0 : }
377 : //
378 0 : printf("Segmentation for Dipole (%+.2f<Z<%+.2f cm)\n",fMinZDip,fMaxZDip);
379 0 : if (fParamsDip) {
380 0 : for (int i=0;i<fNParamsDip;i++) {
381 0 : printf("DIP%4d ",i);
382 0 : GetParamDip(i)->Print();
383 : }
384 0 : }
385 : //
386 0 : }
387 :
388 : //__________________________________________________________________________________________________
389 : Int_t AliMagWrapCheb::FindDipSegment(const Double_t *xyz) const
390 : {
391 : // find the segment containing point xyz. If it is outside find the closest segment
392 42572 : if (!fNParamsDip) return -1;
393 21286 : int xid,yid,zid = TMath::BinarySearch(fNZSegDip,fSegZDip,(Float_t)xyz[2]); // find zsegment
394 : //
395 : Bool_t reCheck = kFALSE;
396 21286 : while(1) {
397 21286 : int ysegBeg = fBegSegYDip[zid];
398 : //
399 362451 : for (yid=0;yid<fNSegYDip[zid];yid++) if (xyz[1]<fSegYDip[ysegBeg+yid]) break;
400 21286 : if ( --yid < 0 ) yid = 0;
401 21286 : yid += ysegBeg;
402 : //
403 21286 : int xsegBeg = fBegSegXDip[yid];
404 297798 : for (xid=0;xid<fNSegXDip[yid];xid++) if (xyz[0]<fSegXDip[xsegBeg+xid]) break;
405 : //
406 21286 : if ( --xid < 0) xid = 0;
407 21286 : xid += xsegBeg;
408 : //
409 : // to make sure that due to the precision problems we did not pick the next Zbin
410 42572 : if (!reCheck && (xyz[2] - fSegZDip[zid] < 3.e-5) && zid &&
411 0 : !GetParamDip(fSegIDDip[xid])->IsInside(xyz)) { // check the previous Z bin
412 0 : zid--;
413 : reCheck = kTRUE;
414 0 : continue;
415 : }
416 21286 : break;
417 : }
418 : // AliInfo(Form("%+.2f %+.2f %+.2f %d %d %d %4d",xyz[0],xyz[1],xyz[2],xid,yid,zid,fSegIDDip[xid]));
419 21286 : return fSegIDDip[xid];
420 21286 : }
421 :
422 : //__________________________________________________________________________________________________
423 : Int_t AliMagWrapCheb::FindSolSegment(const Double_t *rpz) const
424 : {
425 : // find the segment containing point xyz. If it is outside find the closest segment
426 256610 : if (!fNParamsSol) return -1;
427 128305 : int rid,pid,zid = TMath::BinarySearch(fNZSegSol,fSegZSol,(Float_t)rpz[2]); // find zsegment
428 : //
429 : Bool_t reCheck = kFALSE;
430 128305 : while(1) {
431 128306 : int psegBeg = fBegSegPSol[zid];
432 2478892 : for (pid=0;pid<fNSegPSol[zid];pid++) if (rpz[1]<fSegPSol[psegBeg+pid]) break;
433 128306 : if ( --pid < 0 ) pid = 0;
434 128306 : pid += psegBeg;
435 : //
436 128306 : int rsegBeg = fBegSegRSol[pid];
437 1694339 : for (rid=0;rid<fNSegRSol[pid];rid++) if (rpz[0]<fSegRSol[rsegBeg+rid]) break;
438 128306 : if ( --rid < 0) rid = 0;
439 128306 : rid += rsegBeg;
440 : //
441 : // to make sure that due to the precision problems we did not pick the next Zbin
442 264993 : if (!reCheck && (rpz[2] - fSegZSol[zid] < 3.e-5) && zid &&
443 8382 : !GetParamSol(fSegIDSol[rid])->IsInside(rpz)) { // check the previous Z bin
444 1 : zid--;
445 : reCheck = kTRUE;
446 1 : continue;
447 : }
448 128305 : break;
449 : }
450 : // AliInfo(Form("%+.2f %+.4f %+.2f %d %d %d %4d",rpz[0],rpz[1],rpz[2],rid,pid,zid,fSegIDSol[rid]));
451 128305 : return fSegIDSol[rid];
452 128305 : }
453 :
454 : //__________________________________________________________________________________________________
455 : Int_t AliMagWrapCheb::FindTPCSegment(const Double_t *rpz) const
456 : {
457 : // find the segment containing point xyz. If it is outside find the closest segment
458 0 : if (!fNParamsTPC) return -1;
459 0 : int rid,pid,zid = TMath::BinarySearch(fNZSegTPC,fSegZTPC,(Float_t)rpz[2]); // find zsegment
460 : //
461 : Bool_t reCheck = kFALSE;
462 0 : while(1) {
463 0 : int psegBeg = fBegSegPTPC[zid];
464 : //
465 0 : for (pid=0;pid<fNSegPTPC[zid];pid++) if (rpz[1]<fSegPTPC[psegBeg+pid]) break;
466 0 : if ( --pid < 0 ) pid = 0;
467 0 : pid += psegBeg;
468 : //
469 0 : int rsegBeg = fBegSegRTPC[pid];
470 0 : for (rid=0;rid<fNSegRTPC[pid];rid++) if (rpz[0]<fSegRTPC[rsegBeg+rid]) break;
471 0 : if ( --rid < 0) rid = 0;
472 0 : rid += rsegBeg;
473 : //
474 : // to make sure that due to the precision problems we did not pick the next Zbin
475 0 : if (!reCheck && (rpz[2] - fSegZTPC[zid] < 3.e-5) && zid &&
476 0 : !GetParamTPCInt(fSegIDTPC[rid])->IsInside(rpz)) { // check the previous Z bin
477 0 : zid--;
478 : reCheck = kTRUE;
479 0 : continue;
480 : }
481 0 : break;
482 : }
483 : // AliInfo(Form("%+.2f %+.4f %+.2f %d %d %d %4d",rpz[0],rpz[1],rpz[2],rid,pid,zid,fSegIDTPC[rid]));
484 0 : return fSegIDTPC[rid];
485 0 : }
486 :
487 : //__________________________________________________________________________________________________
488 : Int_t AliMagWrapCheb::FindTPCRatSegment(const Double_t *rpz) const
489 : {
490 : // find the segment containing point xyz. If it is outside find the closest segment
491 0 : if (!fNParamsTPCRat) return -1;
492 0 : int rid,pid,zid = TMath::BinarySearch(fNZSegTPCRat,fSegZTPCRat,(Float_t)rpz[2]); // find zsegment
493 : //
494 : Bool_t reCheck = kFALSE;
495 0 : while(1) {
496 0 : int psegBeg = fBegSegPTPCRat[zid];
497 : //
498 0 : for (pid=0;pid<fNSegPTPCRat[zid];pid++) if (rpz[1]<fSegPTPCRat[psegBeg+pid]) break;
499 0 : if ( --pid < 0 ) pid = 0;
500 0 : pid += psegBeg;
501 : //
502 0 : int rsegBeg = fBegSegRTPCRat[pid];
503 0 : for (rid=0;rid<fNSegRTPCRat[pid];rid++) if (rpz[0]<fSegRTPCRat[rsegBeg+rid]) break;
504 0 : if ( --rid < 0) rid = 0;
505 0 : rid += rsegBeg;
506 : //
507 : // to make sure that due to the precision problems we did not pick the next Zbin
508 0 : if (!reCheck && (rpz[2] - fSegZTPCRat[zid] < 3.e-5) && zid &&
509 0 : !GetParamTPCRatInt(fSegIDTPCRat[rid])->IsInside(rpz)) { // check the previous Z bin
510 0 : zid--;
511 : reCheck = kTRUE;
512 0 : continue;
513 : }
514 0 : break;
515 : }
516 : // AliInfo(Form("%+.2f %+.4f %+.2f %d %d %d %4d",rpz[0],rpz[1],rpz[2],rid,pid,zid,fSegIDTPCRat[rid]));
517 0 : return fSegIDTPCRat[rid];
518 0 : }
519 :
520 :
521 : //__________________________________________________________________________________________
522 : void AliMagWrapCheb::GetTPCInt(const Double_t *xyz, Double_t *b) const
523 : {
524 : // compute TPC region field integral in cartesian coordinates.
525 : // If point is outside of the parameterized region get it at closeset valid point
526 : static Double_t rphiz[3];
527 : //
528 : // TPCInt region
529 : // convert coordinates to cyl system
530 0 : CartToCyl(xyz,rphiz);
531 : #ifndef _BRING_TO_BOUNDARY_
532 0 : if ( (rphiz[2]>GetMaxZTPCInt()||rphiz[2]<GetMinZTPCInt()) ||
533 0 : rphiz[0]>GetMaxRTPCInt()) {for (int i=3;i--;) b[i]=0; return;}
534 : #endif
535 : //
536 0 : GetTPCIntCyl(rphiz,b);
537 : //
538 : // convert field to cartesian system
539 0 : CylToCartCylB(rphiz, b,b);
540 : //
541 0 : }
542 :
543 : //__________________________________________________________________________________________
544 : void AliMagWrapCheb::GetTPCRatInt(const Double_t *xyz, Double_t *b) const
545 : {
546 : // compute TPCRat region field integral in cartesian coordinates.
547 : // If point is outside of the parameterized region get it at closeset valid point
548 : static Double_t rphiz[3];
549 : //
550 : // TPCRatInt region
551 : // convert coordinates to cyl system
552 0 : CartToCyl(xyz,rphiz);
553 : #ifndef _BRING_TO_BOUNDARY_
554 0 : if ( (rphiz[2]>GetMaxZTPCRatInt()||rphiz[2]<GetMinZTPCRatInt()) ||
555 0 : rphiz[0]>GetMaxRTPCRatInt()) {for (int i=3;i--;) b[i]=0; return;}
556 : #endif
557 : //
558 0 : GetTPCRatIntCyl(rphiz,b);
559 : //
560 : // convert field to cartesian system
561 0 : CylToCartCylB(rphiz, b,b);
562 : //
563 0 : }
564 :
565 : //__________________________________________________________________________________________
566 : void AliMagWrapCheb::FieldCylSol(const Double_t *rphiz, Double_t *b) const
567 : {
568 : // compute Solenoid field in Cylindircal coordinates
569 : // note: if the point is outside the volume get the field in closest parameterized point
570 224476 : int id = FindSolSegment(rphiz);
571 112238 : if (id>=0) {
572 : #ifndef _MAGCHEB_CACHE_
573 : AliCheb3D* fCacheSol = 0;
574 : #endif
575 112238 : fCacheSol = GetParamSol(id);
576 : //
577 : #ifndef _BRING_TO_BOUNDARY_ // exact matching to fitted volume is requested
578 147227 : if (!fCacheSol->IsInside(rphiz)) return;
579 : #endif
580 77249 : fCacheSol->Eval(rphiz,b);
581 77249 : }
582 : //
583 189487 : }
584 :
585 : //__________________________________________________________________________________________
586 : Double_t AliMagWrapCheb::FieldCylSolBz(const Double_t *rphiz) const
587 : {
588 : // compute Solenoid field in Cylindircal coordinates
589 : // note: if the point is outside the volume get the field in closest parameterized point
590 32134 : int id = FindSolSegment(rphiz);
591 16067 : if (id<0) return 0.;
592 : //
593 : #ifndef _MAGCHEB_CACHE_
594 : AliCheb3D* fCacheSol = 0;
595 : #endif
596 : //
597 16067 : fCacheSol = GetParamSol(id);
598 : #ifndef _BRING_TO_BOUNDARY_
599 48201 : return fCacheSol->IsInside(rphiz) ? fCacheSol->Eval(rphiz,2) : 0;
600 : #else
601 : return fCacheSol->Eval(rphiz,2);
602 : #endif
603 : //
604 16067 : }
605 :
606 : //__________________________________________________________________________________________
607 : void AliMagWrapCheb::GetTPCIntCyl(const Double_t *rphiz, Double_t *b) const
608 : {
609 : // compute field integral in TPC region in Cylindircal coordinates
610 : // note: the check for the point being inside the parameterized region is done outside
611 : //
612 : #ifdef _MAGCHEB_CACHE_
613 : //
614 0 : if (fCacheTPCInt && fCacheTPCInt->IsInside(rphiz)) {
615 0 : fCacheTPCInt->Eval(rphiz,b);
616 0 : return;
617 : }
618 : #else //_MAGCHEB_CACHE_
619 : AliCheb3D* fCacheTPCInt = 0;
620 : #endif //_MAGCHEB_CACHE_
621 : //
622 0 : int id = FindTPCSegment(rphiz);
623 0 : if (id>=0) {
624 : // if (id>=fNParamsTPC) {
625 : // b[0] = b[1] = b[2] = 0;
626 : // AliError(Form("Wrong TPCParam segment %d",id));
627 : // b[0] = b[1] = b[2] = 0;
628 : // return;
629 : // }
630 0 : fCacheTPCInt = GetParamTPCInt(id);
631 0 : if (fCacheTPCInt->IsInside(rphiz)) {
632 0 : fCacheTPCInt->Eval(rphiz,b);
633 0 : return;
634 : }
635 : }
636 : //
637 0 : b[0] = b[1] = b[2] = 0;
638 : //
639 0 : }
640 :
641 : //__________________________________________________________________________________________
642 : void AliMagWrapCheb::GetTPCRatIntCyl(const Double_t *rphiz, Double_t *b) const
643 : {
644 : // compute field integral in TPCRat region in Cylindircal coordinates
645 : // note: the check for the point being inside the parameterized region is done outside
646 : //
647 : #ifdef _MAGCHEB_CACHE_
648 0 : if (fCacheTPCRat && fCacheTPCRat->IsInside(rphiz)) {
649 0 : fCacheTPCRat->Eval(rphiz,b);
650 0 : return;
651 : }
652 : #else
653 : AliCheb3D* fCacheTPCRat = 0;
654 : #endif //_MAGCHEB_CACHE_
655 : //
656 0 : int id = FindTPCRatSegment(rphiz);
657 0 : if (id>=0) {
658 : // if (id>=fNParamsTPCRat) {
659 : // AliError(Form("Wrong TPCRatParam segment %d",id));
660 : // b[0] = b[1] = b[2] = 0;
661 : // return;
662 : // }
663 0 : fCacheTPCRat = GetParamTPCRatInt(id);
664 0 : if (fCacheTPCRat->IsInside(rphiz)) {
665 0 : fCacheTPCRat->Eval(rphiz,b);
666 0 : return;
667 : }
668 : }
669 : //
670 0 : b[0] = b[1] = b[2] = 0;
671 : //
672 0 : }
673 :
674 :
675 : #ifdef _INC_CREATION_ALICHEB3D_
676 : //_______________________________________________
677 : void AliMagWrapCheb::LoadData(const char* inpfile)
678 : {
679 : // read coefficients data from the text file
680 : //
681 0 : TString strf = inpfile;
682 0 : gSystem->ExpandPathName(strf);
683 0 : FILE* stream = fopen(strf,"r");
684 0 : if (!stream) {
685 0 : printf("Did not find input file %s\n",strf.Data());
686 0 : return;
687 : }
688 : //
689 0 : TString buffs;
690 0 : AliCheb3DCalc::ReadLine(buffs,stream);
691 0 : if (!buffs.BeginsWith("START")) AliFatalF("Expected: \"START <name>\", found \"%s\"",buffs.Data());
692 0 : if (buffs.First(' ')>0) SetName(buffs.Data()+buffs.First(' ')+1);
693 : //
694 : // Solenoid part -----------------------------------------------------------
695 0 : AliCheb3DCalc::ReadLine(buffs,stream);
696 0 : if (!buffs.BeginsWith("START SOLENOID")) AliFatalF("Expected: \"START SOLENOID\", found \"%s\"",buffs.Data());
697 0 : AliCheb3DCalc::ReadLine(buffs,stream); // nparam
698 0 : int nparSol = buffs.Atoi();
699 : //
700 0 : for (int ip=0;ip<nparSol;ip++) {
701 0 : AliCheb3D* cheb = new AliCheb3D();
702 0 : cheb->LoadData(stream);
703 0 : AddParamSol(cheb);
704 : }
705 : //
706 0 : AliCheb3DCalc::ReadLine(buffs,stream);
707 0 : if (!buffs.BeginsWith("END SOLENOID")) AliFatalF("Expected \"END SOLENOID\", found \"%s\"",buffs.Data());
708 : //
709 : // TPCInt part -----------------------------------------------------------
710 0 : AliCheb3DCalc::ReadLine(buffs,stream);
711 0 : if (!buffs.BeginsWith("START TPCINT")) AliFatalF("Expected: \"START TPCINT\", found \"%s\"",buffs.Data());
712 : //
713 0 : AliCheb3DCalc::ReadLine(buffs,stream); // nparam
714 0 : int nparTPCInt = buffs.Atoi();
715 : //
716 0 : for (int ip=0;ip<nparTPCInt;ip++) {
717 0 : AliCheb3D* cheb = new AliCheb3D();
718 0 : cheb->LoadData(stream);
719 0 : AddParamTPCInt(cheb);
720 : }
721 : //
722 0 : AliCheb3DCalc::ReadLine(buffs,stream);
723 0 : if (!buffs.BeginsWith("END TPCINT")) AliFatalF("Expected \"END TPCINT\", found \"%s\"",buffs.Data());
724 : //
725 : // TPCRatInt part -----------------------------------------------------------
726 0 : AliCheb3DCalc::ReadLine(buffs,stream);
727 0 : if (!buffs.BeginsWith("START TPCRatINT")) AliFatalF("Expected: \"START TPCRatINT\", found \"%s\"",buffs.Data());
728 0 : AliCheb3DCalc::ReadLine(buffs,stream); // nparam
729 0 : int nparTPCRatInt = buffs.Atoi();
730 : //
731 0 : for (int ip=0;ip<nparTPCRatInt;ip++) {
732 0 : AliCheb3D* cheb = new AliCheb3D();
733 0 : cheb->LoadData(stream);
734 0 : AddParamTPCRatInt(cheb);
735 : }
736 : //
737 0 : AliCheb3DCalc::ReadLine(buffs,stream);
738 0 : if (!buffs.BeginsWith("END TPCRatINT")) AliFatalF("Expected \"END TPCRatINT\", found \"%s\"",buffs.Data());
739 : //
740 : // Dipole part -----------------------------------------------------------
741 0 : AliCheb3DCalc::ReadLine(buffs,stream);
742 0 : if (!buffs.BeginsWith("START DIPOLE")) AliFatalF("Expected: \"START DIPOLE\", found \"%s\"",buffs.Data());
743 : //
744 0 : AliCheb3DCalc::ReadLine(buffs,stream); // nparam
745 0 : int nparDip = buffs.Atoi();
746 : //
747 0 : for (int ip=0;ip<nparDip;ip++) {
748 0 : AliCheb3D* cheb = new AliCheb3D();
749 0 : cheb->LoadData(stream);
750 0 : AddParamDip(cheb);
751 : }
752 : //
753 0 : AliCheb3DCalc::ReadLine(buffs,stream);
754 0 : if (!buffs.BeginsWith("END DIPOLE")) AliFatalF("Expected \"END DIPOLE\", found \"%s\"",buffs.Data());
755 : //
756 0 : AliCheb3DCalc::ReadLine(buffs,stream);
757 0 : if (!buffs.BeginsWith("END ") && !buffs.Contains(GetName())) AliFatalF("Expected: \"END %s\", found \"%s\"",GetName(),buffs.Data());
758 : //
759 : // ---------------------------------------------------------------------------
760 0 : fclose(stream);
761 0 : BuildTableSol();
762 0 : BuildTableDip();
763 0 : BuildTableTPCInt();
764 0 : BuildTableTPCRatInt();
765 : //
766 0 : printf("Loaded magnetic field \"%s\" from %s\n",GetName(),strf.Data());
767 : //
768 0 : }
769 :
770 : //__________________________________________________________________________________________
771 : void AliMagWrapCheb::BuildTableSol()
772 : {
773 : // build lookup table
774 0 : BuildTable(fNParamsSol,fParamsSol,
775 0 : fNZSegSol,fNPSegSol,fNRSegSol,
776 0 : fMinZSol,fMaxZSol,
777 0 : &fSegZSol,&fSegPSol,&fSegRSol,
778 0 : &fBegSegPSol,&fNSegPSol,
779 0 : &fBegSegRSol,&fNSegRSol,
780 0 : &fSegIDSol);
781 0 : }
782 :
783 : //__________________________________________________________________________________________
784 : void AliMagWrapCheb::BuildTableDip()
785 : {
786 : // build lookup table
787 0 : BuildTable(fNParamsDip,fParamsDip,
788 0 : fNZSegDip,fNYSegDip,fNXSegDip,
789 0 : fMinZDip,fMaxZDip,
790 0 : &fSegZDip,&fSegYDip,&fSegXDip,
791 0 : &fBegSegYDip,&fNSegYDip,
792 0 : &fBegSegXDip,&fNSegXDip,
793 0 : &fSegIDDip);
794 0 : }
795 :
796 : //__________________________________________________________________________________________
797 : void AliMagWrapCheb::BuildTableTPCInt()
798 : {
799 : // build lookup table
800 0 : BuildTable(fNParamsTPC,fParamsTPC,
801 0 : fNZSegTPC,fNPSegTPC,fNRSegTPC,
802 0 : fMinZTPC,fMaxZTPC,
803 0 : &fSegZTPC,&fSegPTPC,&fSegRTPC,
804 0 : &fBegSegPTPC,&fNSegPTPC,
805 0 : &fBegSegRTPC,&fNSegRTPC,
806 0 : &fSegIDTPC);
807 0 : }
808 :
809 : //__________________________________________________________________________________________
810 : void AliMagWrapCheb::BuildTableTPCRatInt()
811 : {
812 : // build lookup table
813 0 : BuildTable(fNParamsTPCRat,fParamsTPCRat,
814 0 : fNZSegTPCRat,fNPSegTPCRat,fNRSegTPCRat,
815 0 : fMinZTPCRat,fMaxZTPCRat,
816 0 : &fSegZTPCRat,&fSegPTPCRat,&fSegRTPCRat,
817 0 : &fBegSegPTPCRat,&fNSegPTPCRat,
818 0 : &fBegSegRTPCRat,&fNSegRTPCRat,
819 0 : &fSegIDTPCRat);
820 0 : }
821 :
822 : #endif
823 :
824 : //_______________________________________________
825 : #ifdef _INC_CREATION_ALICHEB3D_
826 :
827 : //__________________________________________________________________________________________
828 0 : AliMagWrapCheb::AliMagWrapCheb(const char* inputFile) :
829 0 : fNParamsSol(0),fNZSegSol(0),fNPSegSol(0),fNRSegSol(0),
830 0 : fSegZSol(0),fSegPSol(0),fSegRSol(0),
831 0 : fBegSegPSol(0),fNSegPSol(0),fBegSegRSol(0),fNSegRSol(0),fSegIDSol(0),fMinZSol(1.e6),fMaxZSol(-1.e6),fParamsSol(0),fMaxRSol(0),
832 : //
833 0 : fNParamsTPC(0),fNZSegTPC(0),fNPSegTPC(0),fNRSegTPC(0),
834 0 : fSegZTPC(0),fSegPTPC(0),fSegRTPC(0),
835 0 : fBegSegPTPC(0),fNSegPTPC(0),fBegSegRTPC(0),fNSegRTPC(0),fSegIDTPC(0),fMinZTPC(1.e6),fMaxZTPC(-1.e6),fParamsTPC(0),fMaxRTPC(0),
836 : //
837 0 : fNParamsTPCRat(0),fNZSegTPCRat(0),fNPSegTPCRat(0),fNRSegTPCRat(0),
838 0 : fSegZTPCRat(0),fSegPTPCRat(0),fSegRTPCRat(0),
839 0 : fBegSegPTPCRat(0),fNSegPTPCRat(0),fBegSegRTPCRat(0),fNSegRTPCRat(0),fSegIDTPCRat(0),fMinZTPCRat(1.e6),fMaxZTPCRat(-1.e6),fParamsTPCRat(0),fMaxRTPCRat(0),
840 : //
841 0 : fNParamsDip(0),fNZSegDip(0),fNYSegDip(0),fNXSegDip(0),
842 0 : fSegZDip(0),fSegYDip(0),fSegXDip(0),
843 0 : fBegSegYDip(0),fNSegYDip(0),fBegSegXDip(0),fNSegXDip(0),fSegIDDip(0),fMinZDip(1.e6),fMaxZDip(-1.e6),fParamsDip(0)
844 : #ifdef _MAGCHEB_CACHE_
845 0 : ,fCacheSol(0),fCacheDip(0),fCacheTPCInt(0),fCacheTPCRat(0)
846 : #endif
847 : //
848 0 : {
849 : // construct from coeffs from the text file
850 0 : LoadData(inputFile);
851 0 : }
852 :
853 : //__________________________________________________________________________________________
854 : void AliMagWrapCheb::AddParamSol(const AliCheb3D* param)
855 : {
856 : // adds new parameterization piece for Sol
857 : // NOTE: pieces must be added strictly in increasing R then increasing Z order
858 : //
859 0 : if (!fParamsSol) fParamsSol = new TObjArray();
860 0 : fParamsSol->Add( (AliCheb3D*)param );
861 0 : fNParamsSol++;
862 0 : if (fMaxRSol<param->GetBoundMax(0)) fMaxRSol = param->GetBoundMax(0);
863 : //
864 0 : }
865 :
866 : //__________________________________________________________________________________________
867 : void AliMagWrapCheb::AddParamTPCInt(const AliCheb3D* param)
868 : {
869 : // adds new parameterization piece for TPCInt
870 : // NOTE: pieces must be added strictly in increasing R then increasing Z order
871 : //
872 0 : if (!fParamsTPC) fParamsTPC = new TObjArray();
873 0 : fParamsTPC->Add( (AliCheb3D*)param);
874 0 : fNParamsTPC++;
875 0 : if (fMaxRTPC<param->GetBoundMax(0)) fMaxRTPC = param->GetBoundMax(0);
876 : //
877 0 : }
878 :
879 : //__________________________________________________________________________________________
880 : void AliMagWrapCheb::AddParamTPCRatInt(const AliCheb3D* param)
881 : {
882 : // adds new parameterization piece for TPCRatInt
883 : // NOTE: pieces must be added strictly in increasing R then increasing Z order
884 : //
885 0 : if (!fParamsTPCRat) fParamsTPCRat = new TObjArray();
886 0 : fParamsTPCRat->Add( (AliCheb3D*)param);
887 0 : fNParamsTPCRat++;
888 0 : if (fMaxRTPCRat<param->GetBoundMax(0)) fMaxRTPCRat = param->GetBoundMax(0);
889 : //
890 0 : }
891 :
892 : //__________________________________________________________________________________________
893 : void AliMagWrapCheb::AddParamDip(const AliCheb3D* param)
894 : {
895 : // adds new parameterization piece for Dipole
896 : //
897 0 : if (!fParamsDip) fParamsDip = new TObjArray();
898 0 : fParamsDip->Add( (AliCheb3D*)param);
899 0 : fNParamsDip++;
900 : //
901 0 : }
902 :
903 : //__________________________________________________________________________________________
904 : void AliMagWrapCheb::ResetDip()
905 : {
906 : // clean Dip field (used for update)
907 0 : if (fNParamsDip) {
908 0 : delete fParamsDip; fParamsDip = 0;
909 0 : delete[] fSegZDip; fSegZDip = 0;
910 0 : delete[] fSegXDip; fSegXDip = 0;
911 0 : delete[] fSegYDip; fSegYDip = 0;
912 0 : delete[] fBegSegYDip; fBegSegYDip = 0;
913 0 : delete[] fNSegYDip; fNSegYDip = 0;
914 0 : delete[] fBegSegXDip; fBegSegXDip = 0;
915 0 : delete[] fNSegXDip; fNSegXDip = 0;
916 0 : delete[] fSegIDDip; fSegIDDip = 0;
917 0 : }
918 0 : fNParamsDip = fNZSegDip = fNXSegDip = fNYSegDip = 0;
919 0 : fMinZDip = 1e6;
920 0 : fMaxZDip = -1e6;
921 : //
922 0 : }
923 :
924 : //__________________________________________________________________________________________
925 : void AliMagWrapCheb::ResetSol()
926 : {
927 : // clean Sol field (used for update)
928 0 : if (fNParamsSol) {
929 0 : delete fParamsSol; fParamsSol = 0;
930 0 : delete[] fSegZSol; fSegZSol = 0;
931 0 : delete[] fSegPSol; fSegPSol = 0;
932 0 : delete[] fSegRSol; fSegRSol = 0;
933 0 : delete[] fBegSegPSol; fBegSegPSol = 0;
934 0 : delete[] fNSegPSol; fNSegPSol = 0;
935 0 : delete[] fBegSegRSol; fBegSegRSol = 0;
936 0 : delete[] fNSegRSol; fNSegRSol = 0;
937 0 : delete[] fSegIDSol; fSegIDSol = 0;
938 0 : }
939 0 : fNParamsSol = fNZSegSol = fNPSegSol = fNRSegSol = 0;
940 0 : fMinZSol = 1e6;
941 0 : fMaxZSol = -1e6;
942 0 : fMaxRSol = 0;
943 : //
944 0 : }
945 :
946 : //__________________________________________________________________________________________
947 : void AliMagWrapCheb::ResetTPCInt()
948 : {
949 : // clean TPC field integral (used for update)
950 0 : if (fNParamsTPC) {
951 0 : delete fParamsTPC; fParamsTPC = 0;
952 0 : delete[] fSegZTPC; fSegZTPC = 0;
953 0 : delete[] fSegPTPC; fSegPTPC = 0;
954 0 : delete[] fSegRTPC; fSegRTPC = 0;
955 0 : delete[] fBegSegPTPC; fBegSegPTPC = 0;
956 0 : delete[] fNSegPTPC; fNSegPTPC = 0;
957 0 : delete[] fBegSegRTPC; fBegSegRTPC = 0;
958 0 : delete[] fNSegRTPC; fNSegRTPC = 0;
959 0 : delete[] fSegIDTPC; fSegIDTPC = 0;
960 0 : }
961 0 : fNParamsTPC = fNZSegTPC = fNPSegTPC = fNRSegTPC = 0;
962 0 : fMinZTPC = 1e6;
963 0 : fMaxZTPC = -1e6;
964 0 : fMaxRTPC = 0;
965 : //
966 0 : }
967 :
968 : //__________________________________________________________________________________________
969 : void AliMagWrapCheb::ResetTPCRatInt()
970 : {
971 : // clean TPCRat field integral (used for update)
972 0 : if (fNParamsTPCRat) {
973 0 : delete fParamsTPCRat; fParamsTPCRat = 0;
974 0 : delete[] fSegZTPCRat; fSegZTPCRat = 0;
975 0 : delete[] fSegPTPCRat; fSegPTPCRat = 0;
976 0 : delete[] fSegRTPCRat; fSegRTPCRat = 0;
977 0 : delete[] fBegSegPTPCRat; fBegSegPTPCRat = 0;
978 0 : delete[] fNSegPTPCRat; fNSegPTPCRat = 0;
979 0 : delete[] fBegSegRTPCRat; fBegSegRTPCRat = 0;
980 0 : delete[] fNSegRTPCRat; fNSegRTPCRat = 0;
981 0 : delete[] fSegIDTPCRat; fSegIDTPCRat = 0;
982 0 : }
983 0 : fNParamsTPCRat = fNZSegTPCRat = fNPSegTPCRat = fNRSegTPCRat = 0;
984 0 : fMinZTPCRat = 1e6;
985 0 : fMaxZTPCRat = -1e6;
986 0 : fMaxRTPCRat = 0;
987 : //
988 0 : }
989 :
990 :
991 : //__________________________________________________
992 : void AliMagWrapCheb::BuildTable(Int_t npar,TObjArray *parArr, Int_t &nZSeg, Int_t &nYSeg, Int_t &nXSeg,
993 : Float_t &minZ,Float_t &maxZ,
994 : Float_t **segZ,Float_t **segY,Float_t **segX,
995 : Int_t **begSegY,Int_t **nSegY,
996 : Int_t **begSegX,Int_t **nSegX,
997 : Int_t **segID)
998 : {
999 : // build lookup table for dipole
1000 : //
1001 0 : if (npar<1) return;
1002 0 : TArrayF segYArr,segXArr;
1003 0 : TArrayI begSegYDipArr,begSegXDipArr;
1004 0 : TArrayI nSegYDipArr,nSegXDipArr;
1005 0 : TArrayI segIDArr;
1006 0 : float *tmpSegZ,*tmpSegY,*tmpSegX;
1007 : //
1008 : // create segmentation in Z
1009 0 : nZSeg = SegmentDimension(&tmpSegZ, parArr, npar, 2, 1,-1, 1,-1, 1,-1) - 1;
1010 0 : nYSeg = 0;
1011 0 : nXSeg = 0;
1012 : //
1013 : // for each Z slice create segmentation in Y
1014 0 : begSegYDipArr.Set(nZSeg);
1015 0 : nSegYDipArr.Set(nZSeg);
1016 0 : float xyz[3];
1017 0 : for (int iz=0;iz<nZSeg;iz++) {
1018 0 : printf("\nZSegment#%d %+e : %+e\n",iz,tmpSegZ[iz],tmpSegZ[iz+1]);
1019 0 : int ny = SegmentDimension(&tmpSegY, parArr, npar, 1,
1020 0 : 1,-1, 1,-1, tmpSegZ[iz],tmpSegZ[iz+1]) - 1;
1021 0 : segYArr.Set(ny + nYSeg);
1022 0 : for (int iy=0;iy<ny;iy++) segYArr[nYSeg+iy] = tmpSegY[iy];
1023 0 : begSegYDipArr[iz] = nYSeg;
1024 0 : nSegYDipArr[iz] = ny;
1025 0 : printf(" Found %d YSegments, to start from %d\n",ny, begSegYDipArr[iz]);
1026 : //
1027 : // for each slice in Z and Y create segmentation in X
1028 0 : begSegXDipArr.Set(nYSeg+ny);
1029 0 : nSegXDipArr.Set(nYSeg+ny);
1030 0 : xyz[2] = (tmpSegZ[iz]+tmpSegZ[iz+1])/2.; // mean Z of this segment
1031 : //
1032 0 : for (int iy=0;iy<ny;iy++) {
1033 0 : int isg = nYSeg+iy;
1034 0 : printf("\n YSegment#%d %+e : %+e\n",iy, tmpSegY[iy],tmpSegY[iy+1]);
1035 0 : int nx = SegmentDimension(&tmpSegX, parArr, npar, 0,
1036 0 : 1,-1, tmpSegY[iy],tmpSegY[iy+1], tmpSegZ[iz],tmpSegZ[iz+1]) - 1;
1037 : //
1038 0 : segXArr.Set(nx + nXSeg);
1039 0 : for (int ix=0;ix<nx;ix++) segXArr[nXSeg+ix] = tmpSegX[ix];
1040 0 : begSegXDipArr[isg] = nXSeg;
1041 0 : nSegXDipArr[isg] = nx;
1042 0 : printf(" Found %d XSegments, to start from %d\n",nx, begSegXDipArr[isg]);
1043 : //
1044 0 : segIDArr.Set(nXSeg+nx);
1045 : //
1046 : // find corresponding params
1047 0 : xyz[1] = (tmpSegY[iy]+tmpSegY[iy+1])/2.; // mean Y of this segment
1048 : //
1049 0 : for (int ix=0;ix<nx;ix++) {
1050 0 : xyz[0] = (tmpSegX[ix]+tmpSegX[ix+1])/2.; // mean X of this segment
1051 0 : for (int ipar=0;ipar<npar;ipar++) {
1052 0 : AliCheb3D* cheb = (AliCheb3D*) parArr->At(ipar);
1053 0 : if (!cheb->IsInside(xyz)) continue;
1054 0 : segIDArr[nXSeg+ix] = ipar;
1055 0 : break;
1056 : }
1057 : }
1058 0 : nXSeg += nx;
1059 : //
1060 0 : delete[] tmpSegX;
1061 : }
1062 0 : delete[] tmpSegY;
1063 0 : nYSeg += ny;
1064 : }
1065 : //
1066 0 : minZ = tmpSegZ[0];
1067 0 : maxZ = tmpSegZ[nZSeg];
1068 0 : (*segZ) = new Float_t[nZSeg];
1069 0 : for (int i=nZSeg;i--;) (*segZ)[i] = tmpSegZ[i];
1070 0 : delete[] tmpSegZ;
1071 : //
1072 0 : (*segY) = new Float_t[nYSeg];
1073 0 : (*segX) = new Float_t[nXSeg];
1074 0 : (*begSegY) = new Int_t[nZSeg];
1075 0 : (*nSegY) = new Int_t[nZSeg];
1076 0 : (*begSegX) = new Int_t[nYSeg];
1077 0 : (*nSegX) = new Int_t[nYSeg];
1078 0 : (*segID) = new Int_t[nXSeg];
1079 : //
1080 0 : for (int i=nYSeg;i--;) (*segY)[i] = segYArr[i];
1081 0 : for (int i=nXSeg;i--;) (*segX)[i] = segXArr[i];
1082 0 : for (int i=nZSeg;i--;) {(*begSegY)[i] = begSegYDipArr[i]; (*nSegY)[i] = nSegYDipArr[i];}
1083 0 : for (int i=nYSeg;i--;) {(*begSegX)[i] = begSegXDipArr[i]; (*nSegX)[i] = nSegXDipArr[i];}
1084 0 : for (int i=nXSeg;i--;) {(*segID)[i] = segIDArr[i];}
1085 : //
1086 0 : }
1087 :
1088 : /*
1089 : //__________________________________________________
1090 : void AliMagWrapCheb::BuildTableDip()
1091 : {
1092 : // build lookup table for dipole
1093 : //
1094 : if (fNParamsDip<1) return;
1095 : TArrayF segY,segX;
1096 : TArrayI begSegYDip,begSegXDip;
1097 : TArrayI nsegYDip,nsegXDip;
1098 : TArrayI segID;
1099 : float *tmpSegZ,*tmpSegY,*tmpSegX;
1100 : //
1101 : // create segmentation in Z
1102 : fNZSegDip = SegmentDimension(&tmpSegZ, fParamsDip, fNParamsDip, 2, 1,-1, 1,-1, 1,-1) - 1;
1103 : fNYSegDip = 0;
1104 : fNXSegDip = 0;
1105 : //
1106 : // for each Z slice create segmentation in Y
1107 : begSegYDip.Set(fNZSegDip);
1108 : nsegYDip.Set(fNZSegDip);
1109 : float xyz[3];
1110 : for (int iz=0;iz<fNZSegDip;iz++) {
1111 : printf("\nZSegment#%d %+e : %+e\n",iz,tmpSegZ[iz],tmpSegZ[iz+1]);
1112 : int ny = SegmentDimension(&tmpSegY, fParamsDip, fNParamsDip, 1,
1113 : 1,-1, 1,-1, tmpSegZ[iz],tmpSegZ[iz+1]) - 1;
1114 : segY.Set(ny + fNYSegDip);
1115 : for (int iy=0;iy<ny;iy++) segY[fNYSegDip+iy] = tmpSegY[iy];
1116 : begSegYDip[iz] = fNYSegDip;
1117 : nsegYDip[iz] = ny;
1118 : printf(" Found %d YSegments, to start from %d\n",ny, begSegYDip[iz]);
1119 : //
1120 : // for each slice in Z and Y create segmentation in X
1121 : begSegXDip.Set(fNYSegDip+ny);
1122 : nsegXDip.Set(fNYSegDip+ny);
1123 : xyz[2] = (tmpSegZ[iz]+tmpSegZ[iz+1])/2.; // mean Z of this segment
1124 : //
1125 : for (int iy=0;iy<ny;iy++) {
1126 : int isg = fNYSegDip+iy;
1127 : printf("\n YSegment#%d %+e : %+e\n",iy, tmpSegY[iy],tmpSegY[iy+1]);
1128 : int nx = SegmentDimension(&tmpSegX, fParamsDip, fNParamsDip, 0,
1129 : 1,-1, tmpSegY[iy],tmpSegY[iy+1], tmpSegZ[iz],tmpSegZ[iz+1]) - 1;
1130 : //
1131 : segX.Set(nx + fNXSegDip);
1132 : for (int ix=0;ix<nx;ix++) segX[fNXSegDip+ix] = tmpSegX[ix];
1133 : begSegXDip[isg] = fNXSegDip;
1134 : nsegXDip[isg] = nx;
1135 : printf(" Found %d XSegments, to start from %d\n",nx, begSegXDip[isg]);
1136 : //
1137 : segID.Set(fNXSegDip+nx);
1138 : //
1139 : // find corresponding params
1140 : xyz[1] = (tmpSegY[iy]+tmpSegY[iy+1])/2.; // mean Y of this segment
1141 : //
1142 : for (int ix=0;ix<nx;ix++) {
1143 : xyz[0] = (tmpSegX[ix]+tmpSegX[ix+1])/2.; // mean X of this segment
1144 : for (int ipar=0;ipar<fNParamsDip;ipar++) {
1145 : AliCheb3D* cheb = (AliCheb3D*) fParamsDip->At(ipar);
1146 : if (!cheb->IsInside(xyz)) continue;
1147 : segID[fNXSegDip+ix] = ipar;
1148 : break;
1149 : }
1150 : }
1151 : fNXSegDip += nx;
1152 : //
1153 : delete[] tmpSegX;
1154 : }
1155 : delete[] tmpSegY;
1156 : fNYSegDip += ny;
1157 : }
1158 : //
1159 : fMinZDip = tmpSegZ[0];
1160 : fMaxZDip = tmpSegZ[fNZSegDip];
1161 : fSegZDip = new Float_t[fNZSegDip];
1162 : for (int i=fNZSegDip;i--;) fSegZDip[i] = tmpSegZ[i];
1163 : delete[] tmpSegZ;
1164 : //
1165 : fSegYDip = new Float_t[fNYSegDip];
1166 : fSegXDip = new Float_t[fNXSegDip];
1167 : fBegSegYDip = new Int_t[fNZSegDip];
1168 : fNSegYDip = new Int_t[fNZSegDip];
1169 : fBegSegXDip = new Int_t[fNYSegDip];
1170 : fNSegXDip = new Int_t[fNYSegDip];
1171 : fSegIDDip = new Int_t[fNXSegDip];
1172 : //
1173 : for (int i=fNYSegDip;i--;) fSegYDip[i] = segY[i];
1174 : for (int i=fNXSegDip;i--;) fSegXDip[i] = segX[i];
1175 : for (int i=fNZSegDip;i--;) {fBegSegYDip[i] = begSegYDip[i]; fNSegYDip[i] = nsegYDip[i];}
1176 : for (int i=fNYSegDip;i--;) {fBegSegXDip[i] = begSegXDip[i]; fNSegXDip[i] = nsegXDip[i];}
1177 : for (int i=fNXSegDip;i--;) {fSegIDDip[i] = segID[i];}
1178 : //
1179 : }
1180 : */
1181 :
1182 : //________________________________________________________________
1183 : void AliMagWrapCheb::SaveData(const char* outfile) const
1184 : {
1185 : // writes coefficients data to output text file
1186 0 : TString strf = outfile;
1187 0 : gSystem->ExpandPathName(strf);
1188 0 : FILE* stream = fopen(strf,"w+");
1189 : //
1190 : // Sol part ---------------------------------------------------------
1191 0 : fprintf(stream,"# Set of Chebyshev parameterizations for ALICE magnetic field\nSTART %s\n",GetName());
1192 0 : fprintf(stream,"START SOLENOID\n#Number of pieces\n%d\n",fNParamsSol);
1193 0 : for (int ip=0;ip<fNParamsSol;ip++) GetParamSol(ip)->SaveData(stream);
1194 0 : fprintf(stream,"#\nEND SOLENOID\n");
1195 : //
1196 : // TPCInt part ---------------------------------------------------------
1197 : // fprintf(stream,"# Set of Chebyshev parameterizations for ALICE magnetic field\nSTART %s\n",GetName());
1198 0 : fprintf(stream,"START TPCINT\n#Number of pieces\n%d\n",fNParamsTPC);
1199 0 : for (int ip=0;ip<fNParamsTPC;ip++) GetParamTPCInt(ip)->SaveData(stream);
1200 0 : fprintf(stream,"#\nEND TPCINT\n");
1201 : //
1202 : // TPCRatInt part ---------------------------------------------------------
1203 : // fprintf(stream,"# Set of Chebyshev parameterizations for ALICE magnetic field\nSTART %s\n",GetName());
1204 0 : fprintf(stream,"START TPCRatINT\n#Number of pieces\n%d\n",fNParamsTPCRat);
1205 0 : for (int ip=0;ip<fNParamsTPCRat;ip++) GetParamTPCRatInt(ip)->SaveData(stream);
1206 0 : fprintf(stream,"#\nEND TPCRatINT\n");
1207 : //
1208 : // Dip part ---------------------------------------------------------
1209 0 : fprintf(stream,"START DIPOLE\n#Number of pieces\n%d\n",fNParamsDip);
1210 0 : for (int ip=0;ip<fNParamsDip;ip++) GetParamDip(ip)->SaveData(stream);
1211 0 : fprintf(stream,"#\nEND DIPOLE\n");
1212 : //
1213 0 : fprintf(stream,"#\nEND %s\n",GetName());
1214 : //
1215 0 : fclose(stream);
1216 : //
1217 0 : }
1218 :
1219 : Int_t AliMagWrapCheb::SegmentDimension(float** seg,const TObjArray* par,int npar, int dim,
1220 : float xmn,float xmx,float ymn,float ymx,float zmn,float zmx)
1221 : {
1222 : // find all boundaries in deimension dim for boxes in given region.
1223 : // if mn>mx for given projection the check is not done for it.
1224 0 : float *tmpC = new float[2*npar];
1225 0 : int *tmpInd = new int[2*npar];
1226 : int nseg0 = 0;
1227 0 : for (int ip=0;ip<npar;ip++) {
1228 0 : AliCheb3D* cheb = (AliCheb3D*) par->At(ip);
1229 0 : if (xmn<xmx && (cheb->GetBoundMin(0)>(xmx+xmn)/2 || cheb->GetBoundMax(0)<(xmn+xmx)/2)) continue;
1230 0 : if (ymn<ymx && (cheb->GetBoundMin(1)>(ymx+ymn)/2 || cheb->GetBoundMax(1)<(ymn+ymx)/2)) continue;
1231 0 : if (zmn<zmx && (cheb->GetBoundMin(2)>(zmx+zmn)/2 || cheb->GetBoundMax(2)<(zmn+zmx)/2)) continue;
1232 : //
1233 0 : tmpC[nseg0++] = cheb->GetBoundMin(dim);
1234 0 : tmpC[nseg0++] = cheb->GetBoundMax(dim);
1235 0 : }
1236 : // range Dim's boundaries in increasing order
1237 0 : TMath::Sort(nseg0,tmpC,tmpInd,kFALSE);
1238 : // count number of really different Z's
1239 : int nseg = 0;
1240 : float cprev = -1e6;
1241 0 : for (int ip=0;ip<nseg0;ip++) {
1242 0 : if (TMath::Abs(cprev-tmpC[ tmpInd[ip] ])>1e-4) {
1243 0 : cprev = tmpC[ tmpInd[ip] ];
1244 0 : nseg++;
1245 0 : }
1246 0 : else tmpInd[ip] = -1; // supress redundant Z
1247 : }
1248 : //
1249 0 : *seg = new float[nseg]; // create final Z segmenations
1250 : nseg = 0;
1251 0 : for (int ip=0;ip<nseg0;ip++) if (tmpInd[ip]>=0) (*seg)[nseg++] = tmpC[ tmpInd[ip] ];
1252 : //
1253 0 : delete[] tmpC;
1254 0 : delete[] tmpInd;
1255 0 : return nseg;
1256 : }
1257 :
1258 : #endif
1259 :
|
