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