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