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