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