Fix needed to connect friends to trees and chains.
[u/mrichter/AliRoot.git] / STEER / AliMagFCheb.cxx
CommitLineData
0eea9d4d 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/* $Id$ */
17
18///////////////////////////////////////////////////////////////////////////////////
19// //
20// Wrapper for the set of mag.field parameterizations by Chebyshev polinomials //
21// To obtain the field in cartesian coordinates/components use //
22// Field(float* xyz, float* bxyz); //
23// For cylindrical coordinates/components: //
24// FieldCyl(float* rphiz, float* brphiz) //
25// //
26// For the moment only the solenoid part is parameterized in the volume defined //
27// by R<500, -550<Z<550 cm //
28// //
29// The region R<423 cm, -343.3<Z<481.3 for 30kA and -343.3<Z<481.3 for 12kA //
30// is parameterized using measured data while outside the Tosca calculation //
31// is used (matched to data on the boundary of the measurements) //
32// //
33// If the querried point is outside the validity region no the return values //
34// for the field components are set to 0. //
35// //
36///////////////////////////////////////////////////////////////////////////////////
37
0bc7b414 38#include "AliLog.h"
0eea9d4d 39#include "AliMagFCheb.h"
40
41ClassImp(AliMagFCheb)
42
43
44
45//__________________________________________________________________________________________
0bc7b414 46AliMagFCheb::AliMagFCheb() :
47 TNamed(),
48 fNParamsSol(0),
49 fNSegZSol(0),
50 fNParamsDip(0),
51 fSegZSol(0),
52 fSegRSol(0),
53 fNSegRSol(0),
54 fSegZIdSol(0),
55 fMinZSol(0.),
56 fMaxZSol(0.),
57 fMaxRSol(0.),
58 fParamsSol(0),
59 fParamsDip(0)
60{
61 Init0();
62}
63
0eea9d4d 64AliMagFCheb::AliMagFCheb(const char* inputFile) :
65 TNamed("Field Map", inputFile),
66 fNParamsSol(0),
67 fNSegZSol(0),
68 fNParamsDip(0),
69 fSegZSol(0),
70 fSegRSol(0),
71 fNSegRSol(0),
72 fSegZIdSol(0),
73 fMinZSol(0.),
74 fMaxZSol(0.),
0bc7b414 75 fMaxRSol(0.),
76 fParamsSol(0),
77 fParamsDip(0)
0eea9d4d 78{
79 Init0();
80 LoadData(inputFile);
81}
82
c437b1a5 83//__________________________________________________________________________________________
84AliMagFCheb::AliMagFCheb(const AliMagFCheb& src)
85 : TNamed(src),
86 fNParamsSol(src.fNParamsSol),
87 fNSegZSol(src.fNSegZSol),
88 fNParamsDip(src.fNParamsDip),
89 fSegZSol(0),
90 fSegRSol(0),
91 fNSegRSol(0),
92 fSegZIdSol(0),
93 fMinZSol(src.fMinZSol),
94 fMaxZSol(src.fMaxZSol),
95 fMaxRSol(src.fMaxRSol),
96 fParamsSol(0),
97 fParamsDip(0)
0bc7b414 98{
c437b1a5 99// Copy constructor
100 if (src.fSegZSol) {
101 fSegZSol = new Float_t[fNSegZSol];
102 for (int i=fNSegZSol;i--;) fSegZSol[i] = src.fSegZSol[i];
103 }
104 if (src.fSegRSol) {
105 fSegRSol = new Float_t[fNParamsSol];
106 for (int i=fNParamsSol;i--;) fSegRSol[i] = src.fSegRSol[i];
107 }
108 if (src.fNSegRSol) {
109 fNSegRSol = new Int_t[fNSegZSol];
110 for (int i=fNSegZSol;i--;) fNSegRSol[i] = src.fNSegRSol[i];
111 }
112 if (src.fSegZIdSol) {
113 fSegZIdSol = new Int_t[fNSegZSol];
114 for (int i=fNSegZSol;i--;) fSegZIdSol[i] = src.fSegZIdSol[i];
115 }
116 if (src.fParamsSol) {
117 fParamsSol = new TObjArray(fNParamsSol);
118 for (int i=0;i<fNParamsSol;i++) {
119 AliCheb3D* pr = src.GetParamSol(i);
120 if (pr) fParamsSol->AddAtAndExpand(new AliCheb3D(*pr),i);
121 }
122 }
123 if (src.fParamsDip) {
124 fParamsDip = new TObjArray(fNParamsDip);
125 for (int i=0;i<fNParamsDip;i++) {
126 AliCheb3D* pr = src.GetParamDip(i);
127 if (pr) fParamsDip->AddAtAndExpand(new AliCheb3D(*pr),i);
128 }
129 }
0bc7b414 130 //
0bc7b414 131}
132
c437b1a5 133
134AliMagFCheb& AliMagFCheb::operator=(const AliMagFCheb& rhs)
0bc7b414 135{
c437b1a5 136// Assignment operator
137 if (this != &rhs) {
138 Clear();
139 SetName(rhs.GetName());
140 SetTitle(rhs.GetTitle());
141 fNParamsSol = rhs.fNParamsSol;
142 fNSegZSol = rhs.fNSegZSol;
143 fMinZSol = rhs.fMinZSol;
144 fMaxZSol = rhs.fMaxZSol;
145 fMaxRSol = rhs.fMaxRSol;
146 fNParamsDip = rhs.fNParamsDip;
147 fSegZSol = fSegRSol = 0;
148 fNSegRSol = fSegZIdSol = 0;
149 fParamsSol = fParamsDip = 0;
150 //
151 if (rhs.fSegZSol) {
152 fSegZSol = new Float_t[fNSegZSol];
153 for (int i=fNSegZSol;i--;) fSegZSol[i] = rhs.fSegZSol[i];
154 }
155 if (rhs.fSegRSol) {
156 fSegRSol = new Float_t[fNParamsSol];
157 for (int i=fNParamsSol;i--;) fSegRSol[i] = rhs.fSegRSol[i];
158 }
159 if (rhs.fNSegRSol) {
160 fNSegRSol = new Int_t[fNSegZSol];
161 for (int i=fNSegZSol;i--;) fNSegRSol[i] = rhs.fNSegRSol[i];
162 }
163 if (rhs.fSegZIdSol) {
164 fSegZIdSol = new Int_t[fNSegZSol];
165 for (int i=fNSegZSol;i--;) fSegZIdSol[i] = rhs.fSegZIdSol[i];
166 }
167 if (rhs.fParamsSol) {
168 fParamsSol = new TObjArray(fNParamsSol);
169 for (int i=0;i<fNParamsSol;i++) {
170 AliCheb3D* pr = rhs.GetParamSol(i);
171 if (pr) fParamsSol->AddAtAndExpand(new AliCheb3D(*pr),i);
172 }
173 }
174 if (rhs.fParamsDip) {
175 fParamsDip = new TObjArray(fNParamsDip);
176 for (int i=0;i<fNParamsDip;i++) {
177 AliCheb3D* pr = rhs.GetParamDip(i);
178 if (pr) fParamsDip->AddAtAndExpand(new AliCheb3D(*pr),i);
179 }
180 }
181 }
182 return *this;
183 //
0bc7b414 184}
185
c437b1a5 186
0eea9d4d 187//__________________________________________________________________________________________
188AliMagFCheb::~AliMagFCheb()
189{
190 if (fNParamsSol) {
191 delete fParamsSol;
192 delete[] fSegZSol;
193 delete[] fSegRSol;
194 delete[] fNSegRSol;
195 delete[] fSegZIdSol;
196 }
197 //
198 // Dipole part ...
199 if (fNParamsDip) {
200 delete fParamsDip;
201 }
202}
203
204//__________________________________________________________________________________________
205void AliMagFCheb::Init0()
206{
207 // Solenoid part
208 fNParamsSol = 0;
209 fNSegZSol = 0;
210 //
211 fSegZSol = 0;
212 fSegRSol = 0;
213 //
214 fNSegRSol = 0;
215 fSegZIdSol = 0;
216 //
217 fMinZSol = fMaxZSol = fMaxRSol = fMaxRSol = 0;
218 fParamsSol = 0;
219 //
220 // Dipole part ...
221 fNParamsDip = 0;
222 fParamsDip = 0;
223 //
224}
225
226//__________________________________________________________________________________________
227void AliMagFCheb::AddParamSol(AliCheb3D* param)
228{
229 // adds new parameterization piece for Sol
230 // NOTE: pieces must be added strictly in increasing R then increasing Z order
231 //
232 if (!fParamsSol) fParamsSol = new TObjArray();
233 fParamsSol->Add(param);
234 fNParamsSol++;
235 //
236}
237
238//__________________________________________________________________________________________
239void AliMagFCheb::AddParamDip(AliCheb3D* param)
240{
241 // adds new parameterization piece for Dipole
242 //
243 if (!fParamsDip) fParamsDip = new TObjArray();
244 fParamsDip->Add(param);
245 fNParamsDip++;
246 //
247}
248
249//__________________________________________________________________________________________
250void AliMagFCheb::BuildTableSol()
251{
252 // build the indexes for each parameterization of Solenoid
253 //
254 const float kSafety=0.001;
255 //
256 fSegRSol = new Float_t[fNParamsSol];
257 float *tmpbufF = new float[fNParamsSol+1];
258 int *tmpbufI = new int[fNParamsSol+1];
259 int *tmpbufI1 = new int[fNParamsSol+1];
260 //
261 // count number of Z slices and number of R slices in each Z slice
262 for (int ip=0;ip<fNParamsSol;ip++) {
263 if (ip==0 || (GetParamSol(ip)->GetBoundMax(2)-GetParamSol(ip-1)->GetBoundMax(2))>kSafety) { // new Z slice
264 tmpbufF[fNSegZSol] = GetParamSol(ip)->GetBoundMax(2); //
265 tmpbufI[fNSegZSol] = 0;
266 tmpbufI1[fNSegZSol++] = ip;
267 }
268 fSegRSol[ip] = GetParamSol(ip)->GetBoundMax(0); // upper R
269 tmpbufI[fNSegZSol-1]++;
270 }
271 //
272 fSegZSol = new Float_t[fNSegZSol];
273 fSegZIdSol = new Int_t[fNSegZSol];
274 fNSegRSol = new Int_t[fNSegZSol];
275 for (int iz=0;iz<fNSegZSol;iz++) {
276 fSegZSol[iz] = tmpbufF[iz];
277 fNSegRSol[iz] = tmpbufI[iz];
278 fSegZIdSol[iz] = tmpbufI1[iz];
279 }
280 //
281 fMinZSol = GetParamSol(0)->GetBoundMin(2);
282 fMaxZSol = GetParamSol(fNParamsSol-1)->GetBoundMax(2);
283 fMaxRSol = GetParamSol(fNParamsSol-1)->GetBoundMax(0);
284 //
285 delete[] tmpbufF;
286 delete[] tmpbufI;
287 delete[] tmpbufI1;
288 //
289 //
290}
291
292//__________________________________________________________________________________________
293void AliMagFCheb::Field(Float_t *xyz, Float_t *b) const
294{
295 // compute field in cartesian coordinates
296 float rphiz[3];
297 if (xyz[2]>GetMaxZSol() || xyz[2]<GetMinZSol()) {for (int i=3;i--;) b[i]=0; return;}
298 //
299 // Sol region
300 // convert coordinates to cyl system
301 rphiz[0] = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
302 rphiz[1] = TMath::ATan2(xyz[1],xyz[0]);
303 rphiz[2] = xyz[2];
304 if (rphiz[0]>GetMaxRSol()) {for (int i=3;i--;) b[i]=0; return;}
305 //
306 FieldCylSol(rphiz,b);
307 //
308 // convert field to cartesian system
309 float btr = TMath::Sqrt(b[0]*b[0]+b[1]*b[1]);
310 float psiPLUSphi = rphiz[1] + TMath::ATan2(b[1],b[0]);
311 b[0] = btr*TMath::Cos(psiPLUSphi);
312 b[1] = btr*TMath::Sin(psiPLUSphi);
313 //
314}
315
316//__________________________________________________________________________________________
317void AliMagFCheb::FieldCylSol(Float_t *rphiz, Float_t *b) const
318{
319 // compute Solenoid field in Cylindircal coordinates
320 // note: the check for the point being inside the parameterized region is done outside
321 float &r = rphiz[0];
322 float &z = rphiz[2];
323 int SolZId = 0;
324 while (z>fSegZSol[SolZId]) ++SolZId; // find Z segment
325 int SolRId = fSegZIdSol[SolZId]; // first R segment for this Z
326 while (r>fSegRSol[SolRId]) ++SolRId; // find R segment
327 GetParamSol( SolRId )->Eval(rphiz,b);
328 //
329}
330
331//__________________________________________________________________________________________
332void AliMagFCheb::Print(Option_t *) const
333{
334 printf("Alice magnetic field parameterized by Chebyshev polynomials\n");
335 printf("Segmentation for Solenoid (%+.2f<Z<%+.2f cm | R<%.2f cm)\n",fMinZSol,fMaxZSol,fMaxRSol);
336 //
337 for (int iz=0;iz<fNSegZSol;iz++) {
338 AliCheb3D* param = GetParamSol( fSegZIdSol[iz] );
339 printf("*** Z Segment %2d (%+7.2f<Z<%+7.2f)\t***\n",iz,param->GetBoundMin(2),param->GetBoundMax(2));
340 for (int ir=0;ir<fNSegRSol[iz];ir++) {
341 param = GetParamSol( fSegZIdSol[iz]+ir );
342 printf(" R Segment %2d (%+7.2f<R<%+7.2f, Precision: %.1e) (ID=%2d)\n",ir, param->GetBoundMin(0),param->GetBoundMax(0),
343 param->GetPrecision(),fSegZIdSol[iz]+ir);
344 }
345 }
346}
347
348//_______________________________________________
349#ifdef _INC_CREATION_ALICHEB3D_
350void AliMagFCheb::SaveData(const char* outfile) const
351{
352 // writes coefficients data to output text file
353 TString strf = outfile;
354 gSystem->ExpandPathName(strf);
355 FILE* stream = fopen(strf,"w+");
356 //
357 // Sol part
358 fprintf(stream,"# Set of Chebyshev parameterizations for ALICE magnetic field\nSTART %s\n",GetName());
359 fprintf(stream,"START SOLENOID\n#Number of pieces\n%d\n",fNParamsSol);
360 for (int ip=0;ip<fNParamsSol;ip++) GetParamSol(ip)->SaveData(stream);
361 fprintf(stream,"#\nEND SOLENOID\n");
362 //
363 // Dip part
364 fprintf(stream,"START DIPOLE\n#Number of pieces\n%d\n",fNParamsDip);
365 for (int ip=0;ip<fNParamsDip;ip++) GetParamDip(ip)->SaveData(stream);
366 fprintf(stream,"#\nEND DIPOLE\n");
367 //
368 fprintf(stream,"#\nEND %s\n",GetName());
369 //
370 fclose(stream);
371 //
372}
373#endif
374
375//_______________________________________________
376void AliMagFCheb::LoadData(const char* inpfile)
377{
378 // read coefficients data from the text file
379 //
380 TString strf = inpfile;
381 gSystem->ExpandPathName(strf);
382 FILE* stream = fopen(strf,"r");
383 if (!stream) {
384 printf("Did not find input file %s\n",strf.Data());
385 return;
386 }
387 //
388 TString buffs;
389 AliCheb3DCalc::ReadLine(buffs,stream);
390 if (!buffs.BeginsWith("START")) {Error("LoadData","Expected: \"START <name>\", found \"%s\"\nStop\n",buffs.Data());exit(1);}
391 if (buffs.First(' ')>0) SetName(buffs.Data()+buffs.First(' ')+1);
392 //
393 // Solenoid part
394 AliCheb3DCalc::ReadLine(buffs,stream);
395 if (!buffs.BeginsWith("START SOLENOID")) {Error("LoadData","Expected: \"START SOLENOID\", found \"%s\"\nStop\n",buffs.Data());exit(1);}
396 AliCheb3DCalc::ReadLine(buffs,stream); // nparam
397 int nparSol = buffs.Atoi();
398 //
399 for (int ip=0;ip<nparSol;ip++) {
400 AliCheb3D* cheb = new AliCheb3D();
401 cheb->LoadData(stream);
402 AddParamSol(cheb);
403 }
404 //
405 AliCheb3DCalc::ReadLine(buffs,stream);
406 if (!buffs.BeginsWith("END SOLENOID")) {Error("LoadData","Expected \"END SOLENOID\", found \"%s\"\nStop\n",buffs.Data());exit(1);}
407 //
408 // Dipole part
409 AliCheb3DCalc::ReadLine(buffs,stream);
410 if (!buffs.BeginsWith("START DIPOLE")) {Error("LoadData","Expected: \"START DIPOLE\", found \"%s\"\nStop\n",buffs.Data());exit(1);}
411 AliCheb3DCalc::ReadLine(buffs,stream); // nparam
412 int nparDip = buffs.Atoi();
413 //
414 for (int ip=0;ip<nparDip;ip++) {
415 AliCheb3D* cheb = new AliCheb3D();
416 cheb->LoadData(stream);
417 AddParamDip(cheb);
418 }
419 //
420 AliCheb3DCalc::ReadLine(buffs,stream);
421 if (!buffs.BeginsWith("END DIPOLE")) {Error("LoadData","Expected \"END DIPOLE\", found \"%s\"\nStop\n",GetName(),buffs.Data());exit(1);}
422 //
423 AliCheb3DCalc::ReadLine(buffs,stream);
424 if (!buffs.BeginsWith("END") || !buffs.Contains(GetName())) {Error("LoadData","Expected: \"END %s\", found \"%s\"\nStop\n",GetName(),buffs.Data());exit(1);}
425 //
426 fclose(stream);
427 BuildTableSol();
428 // BuildDipTable();
429 printf("Loaded magnetic field \"%s\" from %s\n",GetName(),strf.Data());
430 //
431}