]>
Commit | Line | Data |
---|---|---|
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 | ||
41 | ClassImp(AliMagFCheb) | |
42 | ||
43 | ||
44 | ||
45 | //__________________________________________________________________________________________ | |
0bc7b414 | 46 | AliMagFCheb::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 | 64 | AliMagFCheb::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 | //__________________________________________________________________________________________ |
84 | AliMagFCheb::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 | |
134 | AliMagFCheb& 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 | //__________________________________________________________________________________________ |
188 | AliMagFCheb::~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 | //__________________________________________________________________________________________ | |
205 | void 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 | //__________________________________________________________________________________________ | |
227 | void 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 | //__________________________________________________________________________________________ | |
239 | void 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 | //__________________________________________________________________________________________ | |
250 | void 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 | //__________________________________________________________________________________________ | |
293 | void 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 | //__________________________________________________________________________________________ | |
317 | void 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 | //__________________________________________________________________________________________ | |
332 | void 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_ | |
350 | void 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 | //_______________________________________________ | |
376 | void 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 | } |