1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 ///////////////////////////////////////////////////////////////////////////////////
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) //
26 // For the moment only the solenoid part is parameterized in the volume defined //
27 // by R<500, -550<Z<550 cm //
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) //
33 // If the querried point is outside the validity region no the return values //
34 // for the field components are set to 0. //
36 ///////////////////////////////////////////////////////////////////////////////////
39 #include "AliMagFCheb.h"
45 //__________________________________________________________________________________________
46 AliMagFCheb::AliMagFCheb() :
64 AliMagFCheb::AliMagFCheb(const char* inputFile) :
65 TNamed("Field Map", inputFile),
83 //__________________________________________________________________________________________
84 AliMagFCheb::AliMagFCheb(const AliMagFCheb& src)
86 fNParamsSol(src.fNParamsSol),
87 fNSegZSol(src.fNSegZSol),
88 fNParamsDip(src.fNParamsDip),
93 fMinZSol(src.fMinZSol),
94 fMaxZSol(src.fMaxZSol),
95 fMaxRSol(src.fMaxRSol),
101 fSegZSol = new Float_t[fNSegZSol];
102 for (int i=fNSegZSol;i--;) fSegZSol[i] = src.fSegZSol[i];
105 fSegRSol = new Float_t[fNParamsSol];
106 for (int i=fNParamsSol;i--;) fSegRSol[i] = src.fSegRSol[i];
109 fNSegRSol = new Int_t[fNSegZSol];
110 for (int i=fNSegZSol;i--;) fNSegRSol[i] = src.fNSegRSol[i];
112 if (src.fSegZIdSol) {
113 fSegZIdSol = new Int_t[fNSegZSol];
114 for (int i=fNSegZSol;i--;) fSegZIdSol[i] = src.fSegZIdSol[i];
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);
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);
134 AliMagFCheb& AliMagFCheb::operator=(const AliMagFCheb& rhs)
136 // Assignment operator
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;
152 fSegZSol = new Float_t[fNSegZSol];
153 for (int i=fNSegZSol;i--;) fSegZSol[i] = rhs.fSegZSol[i];
156 fSegRSol = new Float_t[fNParamsSol];
157 for (int i=fNParamsSol;i--;) fSegRSol[i] = rhs.fSegRSol[i];
160 fNSegRSol = new Int_t[fNSegZSol];
161 for (int i=fNSegZSol;i--;) fNSegRSol[i] = rhs.fNSegRSol[i];
163 if (rhs.fSegZIdSol) {
164 fSegZIdSol = new Int_t[fNSegZSol];
165 for (int i=fNSegZSol;i--;) fSegZIdSol[i] = rhs.fSegZIdSol[i];
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);
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);
187 //__________________________________________________________________________________________
188 AliMagFCheb::~AliMagFCheb()
204 //__________________________________________________________________________________________
205 void AliMagFCheb::Init0()
217 fMinZSol = fMaxZSol = fMaxRSol = fMaxRSol = 0;
226 //__________________________________________________________________________________________
227 void AliMagFCheb::AddParamSol(AliCheb3D* param)
229 // adds new parameterization piece for Sol
230 // NOTE: pieces must be added strictly in increasing R then increasing Z order
232 if (!fParamsSol) fParamsSol = new TObjArray();
233 fParamsSol->Add(param);
238 //__________________________________________________________________________________________
239 void AliMagFCheb::AddParamDip(AliCheb3D* param)
241 // adds new parameterization piece for Dipole
243 if (!fParamsDip) fParamsDip = new TObjArray();
244 fParamsDip->Add(param);
249 //__________________________________________________________________________________________
250 void AliMagFCheb::BuildTableSol()
252 // build the indexes for each parameterization of Solenoid
254 const float kSafety=0.001;
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];
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;
268 fSegRSol[ip] = GetParamSol(ip)->GetBoundMax(0); // upper R
269 tmpbufI[fNSegZSol-1]++;
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];
281 fMinZSol = GetParamSol(0)->GetBoundMin(2);
282 fMaxZSol = GetParamSol(fNParamsSol-1)->GetBoundMax(2);
283 fMaxRSol = GetParamSol(fNParamsSol-1)->GetBoundMax(0);
292 //__________________________________________________________________________________________
293 void AliMagFCheb::Field(Float_t *xyz, Float_t *b) const
295 // compute field in cartesian coordinates
297 if (xyz[2]>GetMaxZSol() || xyz[2]<GetMinZSol()) {for (int i=3;i--;) b[i]=0; return;}
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]);
304 if (rphiz[0]>GetMaxRSol()) {for (int i=3;i--;) b[i]=0; return;}
306 FieldCylSol(rphiz,b);
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);
316 //__________________________________________________________________________________________
317 void AliMagFCheb::FieldCylSol(Float_t *rphiz, Float_t *b) const
319 // compute Solenoid field in Cylindircal coordinates
320 // note: the check for the point being inside the parameterized region is done outside
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);
331 //__________________________________________________________________________________________
332 void AliMagFCheb::Print(Option_t *) const
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);
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);
348 //_______________________________________________
349 #ifdef _INC_CREATION_ALICHEB3D_
350 void AliMagFCheb::SaveData(const char* outfile) const
352 // writes coefficients data to output text file
353 TString strf = outfile;
354 gSystem->ExpandPathName(strf);
355 FILE* stream = fopen(strf,"w+");
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");
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");
368 fprintf(stream,"#\nEND %s\n",GetName());
375 //_______________________________________________
376 void AliMagFCheb::LoadData(const char* inpfile)
378 // read coefficients data from the text file
380 TString strf = inpfile;
381 gSystem->ExpandPathName(strf);
382 FILE* stream = fopen(strf,"r");
384 printf("Did not find input file %s\n",strf.Data());
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);
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();
399 for (int ip=0;ip<nparSol;ip++) {
400 AliCheb3D* cheb = new AliCheb3D();
401 cheb->LoadData(stream);
405 AliCheb3DCalc::ReadLine(buffs,stream);
406 if (!buffs.BeginsWith("END SOLENOID")) {Error("LoadData","Expected \"END SOLENOID\", found \"%s\"\nStop\n",buffs.Data());exit(1);}
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();
414 for (int ip=0;ip<nparDip;ip++) {
415 AliCheb3D* cheb = new AliCheb3D();
416 cheb->LoadData(stream);
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);}
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);}
429 printf("Loaded magnetic field \"%s\" from %s\n",GetName(),strf.Data());