STEER/AliFieldMap.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 // Class to read the magnetic field map and to calculate
  20 // the magnetic field in a given point. For the moment uses a simple
  21 // linear interpolation between the nodes of a rectangular grid
  22 // Author: Andreas Morsch <andreas.morsch@cern.ch>
  23 //-----------------------------------------------------------------------
  24
  25 #include <TSystem.h>
  26 #include <TVector.h>
  27
  28 #include "AliFieldMap.h"
  29
  30 ClassImp(AliFieldMap)
  31
  32 //_______________________________________________________________________
  33 AliFieldMap::AliFieldMap():
  34   fXbeg(0),
  35   fYbeg(0),
  36   fZbeg(0),
  37   fXend(0),
  38   fYend(0),
  39   fZend(0),
  40   fXdel(0),
  41   fYdel(0),
  42   fZdel(0),
  43   fXdeli(0),
  44   fYdeli(0),
  45   fZdeli(0),
  46   fXn(0),
  47   fYn(0),
  48   fZn(0),
  49   fWriteEnable(0),
  50   fB(0)
  51 {
  52   //
  53   // Standard constructor
  54   //
  55   SetWriteEnable();
  56 }
  57
  58 //_______________________________________________________________________
  59 AliFieldMap::AliFieldMap(const char *name, const char *title):
  60   TNamed(name,title),
  61   fXbeg(0),
  62   fYbeg(0),
  63   fZbeg(0),
  64   fXend(0),
  65   fYend(0),
  66   fZend(0),
  67   fXdel(0),
  68   fYdel(0),
  69   fZdel(0),
  70   fXdeli(0),
  71   fYdeli(0),
  72   fZdeli(0),
  73   fXn(0),
  74   fYn(0),
  75   fZn(0),
  76   fWriteEnable(0),
  77   fB(0)
  78 {
  79   //
  80   // Standard constructor
  81   //
  82   ReadField();
  83   SetWriteEnable();
  84 }
  85
  86 //_______________________________________________________________________
  87 AliFieldMap::~AliFieldMap()
  88 {
  89   //
  90   // Destructor
  91   //
  92   delete fB;
  93 }
  94
  95 //_______________________________________________________________________
  96 AliFieldMap::AliFieldMap(const AliFieldMap &map):
  97   TNamed(map),
  98   fXbeg(0),
  99   fYbeg(0),
 100   fZbeg(0),
 101   fXend(0),
 102   fYend(0),
 103   fZend(0),
 104   fXdel(0),
 105   fYdel(0),
 106   fZdel(0),
 107   fXdeli(0),
 108   fYdeli(0),
 109   fZdeli(0),
 110   fXn(0),
 111   fYn(0),
 112   fZn(0),
 113   fWriteEnable(0),
 114   fB(0)
 115 {
 116   //
 117   // Copy constructor
 118   //
 119   map.Copy(*this);
 120 }
 121
 122 //_______________________________________________________________________
 123 void AliFieldMap::ReadField()
 124 {
 125   //
 126   // Method to read the magnetic field map from file
 127   //
 128   FILE* magfile;
 129   //  FILE* endf = fopen("end.table", "r");
 130   //  FILE* out  = fopen("out", "w");
 131
 132   Int_t   ix, iy, iz, ipx, ipy, ipz;
 133   Float_t bx, by, bz;
 134   char *fname = 0;
 135   printf("%s: Reading Magnetic Field Map %s from file %s\n",
 136          ClassName(),fName.Data(),fTitle.Data());
 137
 138   fname   = gSystem->ExpandPathName(fTitle.Data());
 139   magfile = fopen(fname,"r");
 140   delete [] fname;
 141
 142   fscanf(magfile,"%d %d %d %f %f %f %f %f %f",
 143          &fXn, &fYn, &fZn, &fXbeg, &fYbeg, &fZbeg, &fXdel, &fYdel, &fZdel);
 144
 145   fXdeli = 1./fXdel;
 146   fYdeli = 1./fYdel;
 147   fZdeli = 1./fZdel;
 148   fXend  = fXbeg + (fXn-1)*fXdel;
 149   fYend  = fYbeg + (fYn-1)*fYdel;
 150   fZend  = fZbeg + (fZn-1)*fZdel;
 151
 152   Int_t nDim   = fXn*fYn*fZn;
 153
 154   //  Float_t x,y,z,b;
 155
 156   fB = new TVector(3*nDim);
 157   if (magfile) {
 158       for (ix = 0; ix < fXn; ix++) {
 159           ipx=ix*3*(fZn*fYn);
 160           for (iy = 0; iy < fYn; iy++) {
 161               ipy=ipx+iy*3*fZn;
 162               for (iz = 0; iz < fZn; iz++) {
 163                   ipz=ipy+iz*3;
 164
 165                   if (iz == -1) {
 166 //                    fscanf(endf,"%f %f %f", &bx,&by,&bz);
 167                   } else if (iz > -1) {
 168                       fscanf(magfile," %f %f %f", &bx, &by, &bz);
 169                   } else {
 170                       continue;
 171                   }
 172
 173 //                fscanf(magfile,"%f %f %f %f %f %f %f ",
 174 //                       &x, &y, &z, &bx,&by,&bz, &b);
 175 //                fprintf(out, "%15.8e %15.8e %15.8e \n", bx, by, bz);
 176
 177                   (*fB)(ipz+2) = 10.*bz;
 178                   (*fB)(ipz+1) = 10.*by;
 179                   (*fB)(ipz  ) = 10.*bx;
 180               } //iz
 181           } // iy
 182       } // ix
 183 /*
 184 //
 185 // this part for interpolation between z = 700 and 720 cm to get
 186 // z = 710 cm
 187 //
 188       for (ix = 0; ix < fXn; ix++) {
 189           ipx=ix*3*(fZn*fYn);
 190           for (iy = 0; iy < fYn; iy++) {
 191               ipy=ipx+iy*3*fZn;
 192               Float_t bxx = (Bx(ix,iy,0) + Bx(ix,iy,2))/2.;
 193               Float_t byy = (By(ix,iy,0) + By(ix,iy,2))/2.;
 194               Float_t bzz = (Bz(ix,iy,0) + Bz(ix,iy,2))/2.;
 195               ipz=ipy+3;
 196               (*fB)(ipz+2) = bzz;
 197               (*fB)(ipz+1) = byy;
 198               (*fB)(ipz  ) = bxx;
 199           } // iy
 200       } // ix
 201 */
 202   } else {
 203       printf("%s: File %s not found !\n",ClassName(),fTitle.Data());
 204       exit(1);
 205   } // if mafile
 206 }
 207
 208 //_______________________________________________________________________
 209 // void AliFieldMap::Field(Float_t *x, Float_t *b)
 210 // {
 211 //   //
 212 //   // Use simple interpolation to obtain field at point x
 213 //   //
 214 //     Double_t ratx, raty, ratz, hix, hiy, hiz, ratx1, raty1, ratz1,
 215 //      bhyhz, bhylz, blyhz, blylz, bhz, blz, xl[3];
 216 //     const Double_t kone=1;
 217 //     Int_t ix, iy, iz;
 218 //     b[0]=b[1]=b[2]=0;
 219 //     //
 220
 221 //     xl[0]=TMath::Abs(x[0])-fXbeg;
 222 //     xl[1]=TMath::Abs(x[1])-fYbeg;
 223 //     xl[2]=x[2]-fZbeg;
 224
 225 //     hix=xl[0]*fXdeli;
 226 //     ratx=hix-int(hix);
 227 //     ix=int(hix);
 228
 229 //     hiy=xl[1]*fYdeli;
 230 //     raty=hiy-int(hiy);
 231 //     iy=int(hiy);
 232
 233 //     hiz=xl[2]*fZdeli;
 234 //     ratz=hiz-int(hiz);
 235 //     iz=int(hiz);
 236
 237 //     ratx1=kone-ratx;
 238 //     raty1=kone-raty;
 239 //     ratz1=kone-ratz;
 240
 241 //     bhyhz = Bx(ix  ,iy+1,iz+1)*ratx1+Bx(ix+1,iy+1,iz+1)*ratx;
 242 //     bhylz = Bx(ix  ,iy+1,iz  )*ratx1+Bx(ix+1,iy+1,iz  )*ratx;
 243 //     blyhz = Bx(ix  ,iy  ,iz+1)*ratx1+Bx(ix+1,iy  ,iz+1)*ratx;
 244 //     blylz = Bx(ix  ,iy  ,iz  )*ratx1+Bx(ix+1,iy  ,iz  )*ratx;
 245 //     bhz   = blyhz             *raty1+bhyhz             *raty;
 246 //     blz   = blylz             *raty1+bhylz             *raty;
 247 //     b[0]  = blz               *ratz1+bhz               *ratz;
 248 //     //
 249 //     bhyhz = By(ix  ,iy+1,iz+1)*ratx1+By(ix+1,iy+1,iz+1)*ratx;
 250 //     bhylz = By(ix  ,iy+1,iz  )*ratx1+By(ix+1,iy+1,iz  )*ratx;
 251 //     blyhz = By(ix  ,iy  ,iz+1)*ratx1+By(ix+1,iy  ,iz+1)*ratx;
 252 //     blylz = By(ix  ,iy  ,iz  )*ratx1+By(ix+1,iy  ,iz  )*ratx;
 253 //     bhz   = blyhz             *raty1+bhyhz             *raty;
 254 //     blz   = blylz             *raty1+bhylz             *raty;
 255 //     b[1]  = blz               *ratz1+bhz               *ratz;
 256 //     //
 257 //     bhyhz = Bz(ix  ,iy+1,iz+1)*ratx1+Bz(ix+1,iy+1,iz+1)*ratx;
 258 //     bhylz = Bz(ix  ,iy+1,iz  )*ratx1+Bz(ix+1,iy+1,iz  )*ratx;
 259 //     blyhz = Bz(ix  ,iy  ,iz+1)*ratx1+Bz(ix+1,iy  ,iz+1)*ratx;
 260 //     blylz = Bz(ix  ,iy  ,iz  )*ratx1+Bz(ix+1,iy  ,iz  )*ratx;
 261 //     bhz   = blyhz             *raty1+bhyhz             *raty;
 262 //     blz   = blylz             *raty1+bhylz             *raty;
 263 //     b[2]  = blz               *ratz1+bhz               *ratz;
 264 // }
 265 //_______________________________________________________________________
 266 void AliFieldMap::Field(Float_t *x, Float_t *b)
 267 {
 268   //
 269   // Use simple interpolation to obtain field at point x
 270   // Faster version which uses direct access to the TVector
 271   // instead of Bx(...), By(...), Bz(...) interfaces.
 272   // The interpolation is done first in Z (4 points),
 273   // then in Y (2 points), and finally in X direction (one point).
 274
 275     Double_t h=(x[2]-fZbeg)*fZdeli;
 276     Int_t i=int(h); // ix
 277     register Int_t index=i;
 278     Double_t ratz=h-i;
 279     Double_t ratz1=1-ratz;
 280
 281     h=(TMath::Abs(x[1])-fYbeg)*fYdeli;
 282     i=int(h);       // iy
 283     index+=(i*fZn);
 284     Double_t raty=h-i;
 285     Double_t raty1=1-raty;
 286
 287     h=(TMath::Abs(x[0])-fXbeg)*fXdeli;
 288     i=int(h);       // iz
 289     index+=(i*fZn*fYn);
 290     index*=3; //index now points to the beginning of the field data (ix,iy,iz)
 291     Double_t ratx=h-i;
 292     Double_t ratx1=1-ratx;
 293
 294     // The way we retrive data from fB should be optimized
 295     // with respect to the CPU cache
 296
 297     Float_t bx00=(*fB)[index++]*ratz1;
 298     Float_t by00=(*fB)[index++]*ratz1;
 299     Float_t bz00=(*fB)[index++]*ratz1;
 300     bx00+=(*fB)[index++]*ratz;
 301     by00+=(*fB)[index++]*ratz;
 302     bz00+=(*fB)[index]*ratz;
 303
 304     index+=(3*fZn-5);
 305
 306     Float_t bx01=(*fB)[index++]*ratz1;
 307     Float_t by01=(*fB)[index++]*ratz1;
 308     Float_t bz01=(*fB)[index++]*ratz1;
 309     bx01+=(*fB)[index++]*ratz;
 310     by01+=(*fB)[index++]*ratz;
 311     bz01+=(*fB)[index]*ratz;
 312
 313     index+=(3*fYn*fZn-5);
 314
 315     Float_t bx11=(*fB)[index++]*ratz1;
 316     Float_t by11=(*fB)[index++]*ratz1;
 317     Float_t bz11=(*fB)[index++]*ratz1;
 318     bx11+=(*fB)[index++]*ratz;
 319     by11+=(*fB)[index++]*ratz;
 320     bz11+=(*fB)[index]*ratz;
 321
 322     index-=(3*fZn+5);
 323
 324     Float_t bx10=(*fB)[index++]*ratz1;
 325     Float_t by10=(*fB)[index++]*ratz1;
 326     Float_t bz10=(*fB)[index++]*ratz1;
 327     bx10+=(*fB)[index++]*ratz;
 328     by10+=(*fB)[index++]*ratz;
 329     bz10+=(*fB)[index]*ratz;
 330
 331     b[0]= (bx00*raty1+bx01*raty)*ratx1 +(bx10*raty1+bx11*raty)*ratx;
 332     b[1]= (by00*raty1+by01*raty)*ratx1 +(by10*raty1+by11*raty)*ratx;
 333     b[2]= (bz00*raty1+bz01*raty)*ratx1 +(bz10*raty1+bz11*raty)*ratx;
 334 }
 335
 336 //_______________________________________________________________________
 337 void AliFieldMap::Copy(AliFieldMap & /* magf */) const
 338 {
 339   //
 340   // Copy *this onto magf -- Not implemented
 341   //
 342   Fatal("Copy","Not implemented!\n");
 343 }
 344
 345 //_______________________________________________________________________
 346 AliFieldMap & AliFieldMap::operator =(const AliFieldMap &magf)
 347 {
 348   magf.Copy(*this);
 349   return *this;
 350 }
 351
 352 //_______________________________________________________________________
 353 void AliFieldMap::Streamer(TBuffer &R__b)
 354 {
 355    // Stream an object of class AliFieldMap.
 356     TVector* save = 0;
 357
 358     if (R__b.IsReading()) {
 359         AliFieldMap::Class()->ReadBuffer(R__b, this);
 360     } else {
 361         if (!fWriteEnable) {
 362             save = fB;
 363             fB = 0;
 364         }
 365         AliFieldMap::Class()->WriteBuffer(R__b, this);
 366         if (!fWriteEnable) fB = save;
 367     }
 368 }