First version of the SDD DA calibration classes. AliITSOnlineSDDBase - for measuremen...
[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
0bc7b414 83AliMagFCheb::AliMagFCheb(const AliMagFCheb &cheb) :
84 TNamed(),
85 fNParamsSol(0),
86 fNSegZSol(0),
87 fNParamsDip(0),
88 fSegZSol(0),
89 fSegRSol(0),
90 fNSegRSol(0),
91 fSegZIdSol(0),
92 fMinZSol(0.),
93 fMaxZSol(0.),
94 fMaxRSol(0.),
95 fParamsSol(0),
96 fParamsDip(0)
97{
98 //
99 // Copy constructor for AliMC
100 //
101 cheb.Copy(*this);
102}
103
104//_______________________________________________________________________
105void AliMagFCheb::Copy(TObject &) const
106{
107 //dummy Copy function
108 AliFatal("Not implemented!");
109}
110
0eea9d4d 111//__________________________________________________________________________________________
112AliMagFCheb::~AliMagFCheb()
113{
114 if (fNParamsSol) {
115 delete fParamsSol;
116 delete[] fSegZSol;
117 delete[] fSegRSol;
118 delete[] fNSegRSol;
119 delete[] fSegZIdSol;
120 }
121 //
122 // Dipole part ...
123 if (fNParamsDip) {
124 delete fParamsDip;
125 }
126}
127
128//__________________________________________________________________________________________
129void AliMagFCheb::Init0()
130{
131 // Solenoid part
132 fNParamsSol = 0;
133 fNSegZSol = 0;
134 //
135 fSegZSol = 0;
136 fSegRSol = 0;
137 //
138 fNSegRSol = 0;
139 fSegZIdSol = 0;
140 //
141 fMinZSol = fMaxZSol = fMaxRSol = fMaxRSol = 0;
142 fParamsSol = 0;
143 //
144 // Dipole part ...
145 fNParamsDip = 0;
146 fParamsDip = 0;
147 //
148}
149
150//__________________________________________________________________________________________
151void AliMagFCheb::AddParamSol(AliCheb3D* param)
152{
153 // adds new parameterization piece for Sol
154 // NOTE: pieces must be added strictly in increasing R then increasing Z order
155 //
156 if (!fParamsSol) fParamsSol = new TObjArray();
157 fParamsSol->Add(param);
158 fNParamsSol++;
159 //
160}
161
162//__________________________________________________________________________________________
163void AliMagFCheb::AddParamDip(AliCheb3D* param)
164{
165 // adds new parameterization piece for Dipole
166 //
167 if (!fParamsDip) fParamsDip = new TObjArray();
168 fParamsDip->Add(param);
169 fNParamsDip++;
170 //
171}
172
173//__________________________________________________________________________________________
174void AliMagFCheb::BuildTableSol()
175{
176 // build the indexes for each parameterization of Solenoid
177 //
178 const float kSafety=0.001;
179 //
180 fSegRSol = new Float_t[fNParamsSol];
181 float *tmpbufF = new float[fNParamsSol+1];
182 int *tmpbufI = new int[fNParamsSol+1];
183 int *tmpbufI1 = new int[fNParamsSol+1];
184 //
185 // count number of Z slices and number of R slices in each Z slice
186 for (int ip=0;ip<fNParamsSol;ip++) {
187 if (ip==0 || (GetParamSol(ip)->GetBoundMax(2)-GetParamSol(ip-1)->GetBoundMax(2))>kSafety) { // new Z slice
188 tmpbufF[fNSegZSol] = GetParamSol(ip)->GetBoundMax(2); //
189 tmpbufI[fNSegZSol] = 0;
190 tmpbufI1[fNSegZSol++] = ip;
191 }
192 fSegRSol[ip] = GetParamSol(ip)->GetBoundMax(0); // upper R
193 tmpbufI[fNSegZSol-1]++;
194 }
195 //
196 fSegZSol = new Float_t[fNSegZSol];
197 fSegZIdSol = new Int_t[fNSegZSol];
198 fNSegRSol = new Int_t[fNSegZSol];
199 for (int iz=0;iz<fNSegZSol;iz++) {
200 fSegZSol[iz] = tmpbufF[iz];
201 fNSegRSol[iz] = tmpbufI[iz];
202 fSegZIdSol[iz] = tmpbufI1[iz];
203 }
204 //
205 fMinZSol = GetParamSol(0)->GetBoundMin(2);
206 fMaxZSol = GetParamSol(fNParamsSol-1)->GetBoundMax(2);
207 fMaxRSol = GetParamSol(fNParamsSol-1)->GetBoundMax(0);
208 //
209 delete[] tmpbufF;
210 delete[] tmpbufI;
211 delete[] tmpbufI1;
212 //
213 //
214}
215
216//__________________________________________________________________________________________
217void AliMagFCheb::Field(Float_t *xyz, Float_t *b) const
218{
219 // compute field in cartesian coordinates
220 float rphiz[3];
221 if (xyz[2]>GetMaxZSol() || xyz[2]<GetMinZSol()) {for (int i=3;i--;) b[i]=0; return;}
222 //
223 // Sol region
224 // convert coordinates to cyl system
225 rphiz[0] = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
226 rphiz[1] = TMath::ATan2(xyz[1],xyz[0]);
227 rphiz[2] = xyz[2];
228 if (rphiz[0]>GetMaxRSol()) {for (int i=3;i--;) b[i]=0; return;}
229 //
230 FieldCylSol(rphiz,b);
231 //
232 // convert field to cartesian system
233 float btr = TMath::Sqrt(b[0]*b[0]+b[1]*b[1]);
234 float psiPLUSphi = rphiz[1] + TMath::ATan2(b[1],b[0]);
235 b[0] = btr*TMath::Cos(psiPLUSphi);
236 b[1] = btr*TMath::Sin(psiPLUSphi);
237 //
238}
239
240//__________________________________________________________________________________________
241void AliMagFCheb::FieldCylSol(Float_t *rphiz, Float_t *b) const
242{
243 // compute Solenoid field in Cylindircal coordinates
244 // note: the check for the point being inside the parameterized region is done outside
245 float &r = rphiz[0];
246 float &z = rphiz[2];
247 int SolZId = 0;
248 while (z>fSegZSol[SolZId]) ++SolZId; // find Z segment
249 int SolRId = fSegZIdSol[SolZId]; // first R segment for this Z
250 while (r>fSegRSol[SolRId]) ++SolRId; // find R segment
251 GetParamSol( SolRId )->Eval(rphiz,b);
252 //
253}
254
255//__________________________________________________________________________________________
256void AliMagFCheb::Print(Option_t *) const
257{
258 printf("Alice magnetic field parameterized by Chebyshev polynomials\n");
259 printf("Segmentation for Solenoid (%+.2f<Z<%+.2f cm | R<%.2f cm)\n",fMinZSol,fMaxZSol,fMaxRSol);
260 //
261 for (int iz=0;iz<fNSegZSol;iz++) {
262 AliCheb3D* param = GetParamSol( fSegZIdSol[iz] );
263 printf("*** Z Segment %2d (%+7.2f<Z<%+7.2f)\t***\n",iz,param->GetBoundMin(2),param->GetBoundMax(2));
264 for (int ir=0;ir<fNSegRSol[iz];ir++) {
265 param = GetParamSol( fSegZIdSol[iz]+ir );
266 printf(" R Segment %2d (%+7.2f<R<%+7.2f, Precision: %.1e) (ID=%2d)\n",ir, param->GetBoundMin(0),param->GetBoundMax(0),
267 param->GetPrecision(),fSegZIdSol[iz]+ir);
268 }
269 }
270}
271
272//_______________________________________________
273#ifdef _INC_CREATION_ALICHEB3D_
274void AliMagFCheb::SaveData(const char* outfile) const
275{
276 // writes coefficients data to output text file
277 TString strf = outfile;
278 gSystem->ExpandPathName(strf);
279 FILE* stream = fopen(strf,"w+");
280 //
281 // Sol part
282 fprintf(stream,"# Set of Chebyshev parameterizations for ALICE magnetic field\nSTART %s\n",GetName());
283 fprintf(stream,"START SOLENOID\n#Number of pieces\n%d\n",fNParamsSol);
284 for (int ip=0;ip<fNParamsSol;ip++) GetParamSol(ip)->SaveData(stream);
285 fprintf(stream,"#\nEND SOLENOID\n");
286 //
287 // Dip part
288 fprintf(stream,"START DIPOLE\n#Number of pieces\n%d\n",fNParamsDip);
289 for (int ip=0;ip<fNParamsDip;ip++) GetParamDip(ip)->SaveData(stream);
290 fprintf(stream,"#\nEND DIPOLE\n");
291 //
292 fprintf(stream,"#\nEND %s\n",GetName());
293 //
294 fclose(stream);
295 //
296}
297#endif
298
299//_______________________________________________
300void AliMagFCheb::LoadData(const char* inpfile)
301{
302 // read coefficients data from the text file
303 //
304 TString strf = inpfile;
305 gSystem->ExpandPathName(strf);
306 FILE* stream = fopen(strf,"r");
307 if (!stream) {
308 printf("Did not find input file %s\n",strf.Data());
309 return;
310 }
311 //
312 TString buffs;
313 AliCheb3DCalc::ReadLine(buffs,stream);
314 if (!buffs.BeginsWith("START")) {Error("LoadData","Expected: \"START <name>\", found \"%s\"\nStop\n",buffs.Data());exit(1);}
315 if (buffs.First(' ')>0) SetName(buffs.Data()+buffs.First(' ')+1);
316 //
317 // Solenoid part
318 AliCheb3DCalc::ReadLine(buffs,stream);
319 if (!buffs.BeginsWith("START SOLENOID")) {Error("LoadData","Expected: \"START SOLENOID\", found \"%s\"\nStop\n",buffs.Data());exit(1);}
320 AliCheb3DCalc::ReadLine(buffs,stream); // nparam
321 int nparSol = buffs.Atoi();
322 //
323 for (int ip=0;ip<nparSol;ip++) {
324 AliCheb3D* cheb = new AliCheb3D();
325 cheb->LoadData(stream);
326 AddParamSol(cheb);
327 }
328 //
329 AliCheb3DCalc::ReadLine(buffs,stream);
330 if (!buffs.BeginsWith("END SOLENOID")) {Error("LoadData","Expected \"END SOLENOID\", found \"%s\"\nStop\n",buffs.Data());exit(1);}
331 //
332 // Dipole part
333 AliCheb3DCalc::ReadLine(buffs,stream);
334 if (!buffs.BeginsWith("START DIPOLE")) {Error("LoadData","Expected: \"START DIPOLE\", found \"%s\"\nStop\n",buffs.Data());exit(1);}
335 AliCheb3DCalc::ReadLine(buffs,stream); // nparam
336 int nparDip = buffs.Atoi();
337 //
338 for (int ip=0;ip<nparDip;ip++) {
339 AliCheb3D* cheb = new AliCheb3D();
340 cheb->LoadData(stream);
341 AddParamDip(cheb);
342 }
343 //
344 AliCheb3DCalc::ReadLine(buffs,stream);
345 if (!buffs.BeginsWith("END DIPOLE")) {Error("LoadData","Expected \"END DIPOLE\", found \"%s\"\nStop\n",GetName(),buffs.Data());exit(1);}
346 //
347 AliCheb3DCalc::ReadLine(buffs,stream);
348 if (!buffs.BeginsWith("END") || !buffs.Contains(GetName())) {Error("LoadData","Expected: \"END %s\", found \"%s\"\nStop\n",GetName(),buffs.Data());exit(1);}
349 //
350 fclose(stream);
351 BuildTableSol();
352 // BuildDipTable();
353 printf("Loaded magnetic field \"%s\" from %s\n",GetName(),strf.Data());
354 //
355}