]>
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 | ||
38 | #include "AliMagFCheb.h" | |
39 | ||
40 | ClassImp(AliMagFCheb) | |
41 | ||
42 | ||
43 | ||
44 | //__________________________________________________________________________________________ | |
45 | AliMagFCheb::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 | //__________________________________________________________________________________________ | |
63 | AliMagFCheb::~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 | //__________________________________________________________________________________________ | |
80 | void 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 | //__________________________________________________________________________________________ | |
102 | void 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 | //__________________________________________________________________________________________ | |
114 | void 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 | //__________________________________________________________________________________________ | |
125 | void 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 | //__________________________________________________________________________________________ | |
168 | void 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 | //__________________________________________________________________________________________ | |
192 | void 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 | //__________________________________________________________________________________________ | |
207 | void 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_ | |
225 | void 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 | //_______________________________________________ | |
251 | void 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 | } |