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