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