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