STEER/AliMagFCheb.cxx

   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 "AliLog.h"
  39 #include "AliMagFCheb.h"
  40
  41 ClassImp(AliMagFCheb)
  42
  43
  44
  45 //__________________________________________________________________________________________
  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
  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.),
  75     fMaxRSol(0.),
  76     fParamsSol(0),
  77     fParamsDip(0)
  78 {
  79   Init0();
  80   LoadData(inputFile);
  81 }
  82
  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)
  98 {
  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   }
 130   //
 131 }
 132
 133
 134 AliMagFCheb& AliMagFCheb::operator=(const AliMagFCheb& rhs)
 135 {
 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     //
 184 }
 185
 186
 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 }