]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/AliMagFCheb.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Moving AliDSCValue, AliDCSSensor And AliDCSSensorArray from CDB to STEER to avoud...
[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
38#include "AliMagFCheb.h"
39
40ClassImp(AliMagFCheb)
41
42
43
44//__________________________________________________________________________________________
45AliMagFCheb::AliMagFCheb(const char* inputFile) :
46 TNamed("Field Map", inputFile),
47 fNParamsSol(0),
48 fNSegZSol(0),
49 fNParamsDip(0),
50 fSegZSol(0),
51 fSegRSol(0),
52 fNSegRSol(0),
53 fSegZIdSol(0),
54 fMinZSol(0.),
55 fMaxZSol(0.),
56 fMaxRSol(0.)
57{
58 Init0();
59 LoadData(inputFile);
60}
61
62//__________________________________________________________________________________________
63AliMagFCheb::~AliMagFCheb()
64{
65 if (fNParamsSol) {
66 delete fParamsSol;
67 delete[] fSegZSol;
68 delete[] fSegRSol;
69 delete[] fNSegRSol;
70 delete[] fSegZIdSol;
71 }
72 //
73 // Dipole part ...
74 if (fNParamsDip) {
75 delete fParamsDip;
76 }
77}
78
79//__________________________________________________________________________________________
80void AliMagFCheb::Init0()
81{
82 // Solenoid part
83 fNParamsSol = 0;
84 fNSegZSol = 0;
85 //
86 fSegZSol = 0;
87 fSegRSol = 0;
88 //
89 fNSegRSol = 0;
90 fSegZIdSol = 0;
91 //
92 fMinZSol = fMaxZSol = fMaxRSol = fMaxRSol = 0;
93 fParamsSol = 0;
94 //
95 // Dipole part ...
96 fNParamsDip = 0;
97 fParamsDip = 0;
98 //
99}
100
101//__________________________________________________________________________________________
102void AliMagFCheb::AddParamSol(AliCheb3D* param)
103{
104 // adds new parameterization piece for Sol
105 // NOTE: pieces must be added strictly in increasing R then increasing Z order
106 //
107 if (!fParamsSol) fParamsSol = new TObjArray();
108 fParamsSol->Add(param);
109 fNParamsSol++;
110 //
111}
112
113//__________________________________________________________________________________________
114void AliMagFCheb::AddParamDip(AliCheb3D* param)
115{
116 // adds new parameterization piece for Dipole
117 //
118 if (!fParamsDip) fParamsDip = new TObjArray();
119 fParamsDip->Add(param);
120 fNParamsDip++;
121 //
122}
123
124//__________________________________________________________________________________________
125void AliMagFCheb::BuildTableSol()
126{
127 // build the indexes for each parameterization of Solenoid
128 //
129 const float kSafety=0.001;
130 //
131 fSegRSol = new Float_t[fNParamsSol];
132 float *tmpbufF = new float[fNParamsSol+1];
133 int *tmpbufI = new int[fNParamsSol+1];
134 int *tmpbufI1 = new int[fNParamsSol+1];
135 //
136 // count number of Z slices and number of R slices in each Z slice
137 for (int ip=0;ip<fNParamsSol;ip++) {
138 if (ip==0 || (GetParamSol(ip)->GetBoundMax(2)-GetParamSol(ip-1)->GetBoundMax(2))>kSafety) { // new Z slice
139 tmpbufF[fNSegZSol] = GetParamSol(ip)->GetBoundMax(2); //
140 tmpbufI[fNSegZSol] = 0;
141 tmpbufI1[fNSegZSol++] = ip;
142 }
143 fSegRSol[ip] = GetParamSol(ip)->GetBoundMax(0); // upper R
144 tmpbufI[fNSegZSol-1]++;
145 }
146 //
147 fSegZSol = new Float_t[fNSegZSol];
148 fSegZIdSol = new Int_t[fNSegZSol];
149 fNSegRSol = new Int_t[fNSegZSol];
150 for (int iz=0;iz<fNSegZSol;iz++) {
151 fSegZSol[iz] = tmpbufF[iz];
152 fNSegRSol[iz] = tmpbufI[iz];
153 fSegZIdSol[iz] = tmpbufI1[iz];
154 }
155 //
156 fMinZSol = GetParamSol(0)->GetBoundMin(2);
157 fMaxZSol = GetParamSol(fNParamsSol-1)->GetBoundMax(2);
158 fMaxRSol = GetParamSol(fNParamsSol-1)->GetBoundMax(0);
159 //
160 delete[] tmpbufF;
161 delete[] tmpbufI;
162 delete[] tmpbufI1;
163 //
164 //
165}
166
167//__________________________________________________________________________________________
168void AliMagFCheb::Field(Float_t *xyz, Float_t *b) const
169{
170 // compute field in cartesian coordinates
171 float rphiz[3];
172 if (xyz[2]>GetMaxZSol() || xyz[2]<GetMinZSol()) {for (int i=3;i--;) b[i]=0; return;}
173 //
174 // Sol region
175 // convert coordinates to cyl system
176 rphiz[0] = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
177 rphiz[1] = TMath::ATan2(xyz[1],xyz[0]);
178 rphiz[2] = xyz[2];
179 if (rphiz[0]>GetMaxRSol()) {for (int i=3;i--;) b[i]=0; return;}
180 //
181 FieldCylSol(rphiz,b);
182 //
183 // convert field to cartesian system
184 float btr = TMath::Sqrt(b[0]*b[0]+b[1]*b[1]);
185 float psiPLUSphi = rphiz[1] + TMath::ATan2(b[1],b[0]);
186 b[0] = btr*TMath::Cos(psiPLUSphi);
187 b[1] = btr*TMath::Sin(psiPLUSphi);
188 //
189}
190
191//__________________________________________________________________________________________
192void AliMagFCheb::FieldCylSol(Float_t *rphiz, Float_t *b) const
193{
194 // compute Solenoid field in Cylindircal coordinates
195 // note: the check for the point being inside the parameterized region is done outside
196 float &r = rphiz[0];
197 float &z = rphiz[2];
198 int SolZId = 0;
199 while (z>fSegZSol[SolZId]) ++SolZId; // find Z segment
200 int SolRId = fSegZIdSol[SolZId]; // first R segment for this Z
201 while (r>fSegRSol[SolRId]) ++SolRId; // find R segment
202 GetParamSol( SolRId )->Eval(rphiz,b);
203 //
204}
205
206//__________________________________________________________________________________________
207void AliMagFCheb::Print(Option_t *) const
208{
209 printf("Alice magnetic field parameterized by Chebyshev polynomials\n");
210 printf("Segmentation for Solenoid (%+.2f<Z<%+.2f cm | R<%.2f cm)\n",fMinZSol,fMaxZSol,fMaxRSol);
211 //
212 for (int iz=0;iz<fNSegZSol;iz++) {
213 AliCheb3D* param = GetParamSol( fSegZIdSol[iz] );
214 printf("*** Z Segment %2d (%+7.2f<Z<%+7.2f)\t***\n",iz,param->GetBoundMin(2),param->GetBoundMax(2));
215 for (int ir=0;ir<fNSegRSol[iz];ir++) {
216 param = GetParamSol( fSegZIdSol[iz]+ir );
217 printf(" R Segment %2d (%+7.2f<R<%+7.2f, Precision: %.1e) (ID=%2d)\n",ir, param->GetBoundMin(0),param->GetBoundMax(0),
218 param->GetPrecision(),fSegZIdSol[iz]+ir);
219 }
220 }
221}
222
223//_______________________________________________
224#ifdef _INC_CREATION_ALICHEB3D_
225void AliMagFCheb::SaveData(const char* outfile) const
226{
227 // writes coefficients data to output text file
228 TString strf = outfile;
229 gSystem->ExpandPathName(strf);
230 FILE* stream = fopen(strf,"w+");
231 //
232 // Sol part
233 fprintf(stream,"# Set of Chebyshev parameterizations for ALICE magnetic field\nSTART %s\n",GetName());
234 fprintf(stream,"START SOLENOID\n#Number of pieces\n%d\n",fNParamsSol);
235 for (int ip=0;ip<fNParamsSol;ip++) GetParamSol(ip)->SaveData(stream);
236 fprintf(stream,"#\nEND SOLENOID\n");
237 //
238 // Dip part
239 fprintf(stream,"START DIPOLE\n#Number of pieces\n%d\n",fNParamsDip);
240 for (int ip=0;ip<fNParamsDip;ip++) GetParamDip(ip)->SaveData(stream);
241 fprintf(stream,"#\nEND DIPOLE\n");
242 //
243 fprintf(stream,"#\nEND %s\n",GetName());
244 //
245 fclose(stream);
246 //
247}
248#endif
249
250//_______________________________________________
251void AliMagFCheb::LoadData(const char* inpfile)
252{
253 // read coefficients data from the text file
254 //
255 TString strf = inpfile;
256 gSystem->ExpandPathName(strf);
257 FILE* stream = fopen(strf,"r");
258 if (!stream) {
259 printf("Did not find input file %s\n",strf.Data());
260 return;
261 }
262 //
263 TString buffs;
264 AliCheb3DCalc::ReadLine(buffs,stream);
265 if (!buffs.BeginsWith("START")) {Error("LoadData","Expected: \"START <name>\", found \"%s\"\nStop\n",buffs.Data());exit(1);}
266 if (buffs.First(' ')>0) SetName(buffs.Data()+buffs.First(' ')+1);
267 //
268 // Solenoid part
269 AliCheb3DCalc::ReadLine(buffs,stream);
270 if (!buffs.BeginsWith("START SOLENOID")) {Error("LoadData","Expected: \"START SOLENOID\", found \"%s\"\nStop\n",buffs.Data());exit(1);}
271 AliCheb3DCalc::ReadLine(buffs,stream); // nparam
272 int nparSol = buffs.Atoi();
273 //
274 for (int ip=0;ip<nparSol;ip++) {
275 AliCheb3D* cheb = new AliCheb3D();
276 cheb->LoadData(stream);
277 AddParamSol(cheb);
278 }
279 //
280 AliCheb3DCalc::ReadLine(buffs,stream);
281 if (!buffs.BeginsWith("END SOLENOID")) {Error("LoadData","Expected \"END SOLENOID\", found \"%s\"\nStop\n",buffs.Data());exit(1);}
282 //
283 // Dipole part
284 AliCheb3DCalc::ReadLine(buffs,stream);
285 if (!buffs.BeginsWith("START DIPOLE")) {Error("LoadData","Expected: \"START DIPOLE\", found \"%s\"\nStop\n",buffs.Data());exit(1);}
286 AliCheb3DCalc::ReadLine(buffs,stream); // nparam
287 int nparDip = buffs.Atoi();
288 //
289 for (int ip=0;ip<nparDip;ip++) {
290 AliCheb3D* cheb = new AliCheb3D();
291 cheb->LoadData(stream);
292 AddParamDip(cheb);
293 }
294 //
295 AliCheb3DCalc::ReadLine(buffs,stream);
296 if (!buffs.BeginsWith("END DIPOLE")) {Error("LoadData","Expected \"END DIPOLE\", found \"%s\"\nStop\n",GetName(),buffs.Data());exit(1);}
297 //
298 AliCheb3DCalc::ReadLine(buffs,stream);
299 if (!buffs.BeginsWith("END") || !buffs.Contains(GetName())) {Error("LoadData","Expected: \"END %s\", found \"%s\"\nStop\n",GetName(),buffs.Data());exit(1);}
300 //
301 fclose(stream);
302 BuildTableSol();
303 // BuildDipTable();
304 printf("Loaded magnetic field \"%s\" from %s\n",GetName(),strf.Data());
305 //
306}