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