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