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