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