---

BfldMapRect2d.cxx

Go to the documentation of this file.

// $Id: BfldMapRect2d.cxx,v 1.8 2002/06/13 15:02:10 rhatcher Exp $
// 
// BfldMapRect2d 
// 
// This structure provides services for holding the magnetic field components
// at a series of points (cf. BfldMesh).  It encapsulates knowledge 
// about the field values themselves, an possibly how to adjust them
// for external effects (chemistry, coil current, etc.) but doesn't
// explicitly hold any knowledge about the relative spatial placement
// of the field values; that it must glean from BfldMesh.
//
// This rectangular 2-d grid specialization actually does know info
// about the position info (as it was in the file that contained the
// field values), but it provides an interface so that BfldMeshRect2d
// can access it.
//
// Author:  R. Hatcher 2000.06.20
//
#include "BField/BfldMapRect2d.h"

#include "BField/BfldMeshRect2d.h"
#include "Conventions/Munits.h"

// necessary so that one can check "InheritsFrom"
#include "TClass.h"
// for expanding filename
#include "TSystem.h"
#include "TMath.h"

#include "MessageService/MsgService.h"
CVSID("$Id: BfldMapRect2d.cxx,v 1.8 2002/06/13 15:02:10 rhatcher Exp $");

#include <iostream>
#include <fstream>
#include <iomanip>

ClassImp(BfldMapRect2d)

//______________________________________________________________________________
BfldMapRect2d::BfldMapRect2d() 
   : BfldMap()
{

}

//______________________________________________________________________________
BfldMapRect2d::BfldMapRect2d(BfldGrid::Grid_t grid, Int_t variant) 
   : BfldMap(grid,variant), fReadOk(kFALSE)
{
   // read in a Map

   // strip off "-" sign if necessary
   Int_t varabs = TMath::Abs(variant);

   MSG("Bfld",Msg::kInfo) << "BfldMapRect2d ctor map " << varabs << endl;

   char string[255];

   // format and expand file name
   sprintf(string,"$BMAPPATH/bfld_%3d.dat",varabs);
   fFilename = string;
   gSystem->ExpandPathName(fFilename);

   // open the file for reading
   ifstream from;
   from.open(fFilename,ios::in);
   if (!from) {
      MSG("Bfld",Msg::kError) << "can not open input file: '" 
                              << fFilename <<"'"<< endl;
      return;
   }

   // read in the header line
   from.getline(string,sizeof(string),'\n');
   fHeader = string;

   
   MSG("Bfld",Msg::kInfo) << "read file: " << fFilename << endl;
   MSG("Bfld",Msg::kInfo) << "   header: " << fHeader << endl;

   from >> fX0 >> fDx >> fNx >> fY0 >> fDy >> fNy;
   // lengths read in from map files are in cm, convert to base units
   fX0 *= Munits::cm;
   fDx *= Munits::cm;
   fY0 *= Munits::cm;
   fDy *= Munits::cm;


   MSG("Bfld",Msg::kInfo) 
      << " x0 " << fX0 << " dx " << fDx << " nx " << fNx 
      << " y0 " << fY0 << " dy " << fDy << " ny " << fNy << endl;

   // read in remainder of the line for any quadflags
   from.getline(string,sizeof(string),'\n');
   TString quads = string;

   // first 19 characters are useless junk (z0,dz,nz)
   // replace "x" with " "
   quads = quads.Remove(0,19);
   quads = quads.ReplaceAll("x"," ");
   // quads should be of the form ' 00 00 00 00  FNAL2'

   fFileForm = 0;
   if (quads.Contains("FNAL"))  fFileForm = 1;
   if (quads.Contains("FNAL2")) fFileForm = 2;
   MSG("Bfld",Msg::kInfo) << "fileform " << fFileForm << endl;

   // remove any "FNAL" flag
   if (quads.Length() > 12) quads = quads.Remove(12);

   // scan the remaining part for the quad flags
   if (quads.Length() > 0) {
      sscanf(quads,"%3d%3d%3d%3d",
             fQuadFlags,fQuadFlags+1,fQuadFlags+2,fQuadFlags+3);
   } else {
      MSG("Bfld",Msg::kInfo) << "no quads info - use defaults" << endl;
      fQuadFlags[0] = 0;
      fQuadFlags[1] = 2;
      fQuadFlags[2] = 3;
      fQuadFlags[3] = 1;
   }

   // if x0=0,y0=0 assume map only has one quadrant 
   fQuadrant = ( 0 == fX0 && 0 == fY0 );

   MSG("Bfld",Msg::kInfo) << "quadflags (" 
                          << fQuadFlags[0] << "," 
                          << fQuadFlags[1] << "," 
                          << fQuadFlags[2] << "," 
                          << fQuadFlags[3] << ")  use them " 
                          << (fQuadrant ? "'true'" : "'false'") << endl;

   // allocate a container
   fFieldVectors = new TClonesArray("TVector3",fNx*fNy);

   Int_t   ii, indx;
   Float_t x,y,z;
   Float_t bx,by,bz,bmag;

   const Float_t GaussToKGauss = 0.001;
   const Float_t TeslaToKGauss = 10.;

   // start reading the file
   for (Int_t i=0; i < fNx*fNy; i++) {

      switch (fFileForm) {
      case 0:
         from >> x >> y >> z >> bx >> by >> bz >> bmag;
         // LLNL maps were in gauss ... need kGauss
         bx   = bx   * GaussToKGauss;
         by   = by   * GaussToKGauss;
         bz   = bz   * GaussToKGauss;
         bmag = bmag * GaussToKGauss;
         break;
      case 1:
         from >> ii >> x >> y >> bx >> by;
         // FNAL (version 1) maps were in Tesla
         bx   = bx   * TeslaToKGauss;
         by   = by   * TeslaToKGauss;
         bz   = 0.0;
         break;
      case 2:
         from >> x >> y >> z >> bx >> by >> bz >> bmag;
         // FNAL (version 2) maps were in Tesla
         bx   = bx   * TeslaToKGauss;
         by   = by   * TeslaToKGauss;
         bz   = bz   * TeslaToKGauss;
         bmag = bmag * TeslaToKGauss;
         break;
      default:
         MSG("Bfld",Msg::kFatal) << " bad FileForm " << fFileForm << endl;
         break;
      }
      // the above puts field values in kGauss
      // convert them to base units
      bx   *= Munits::kilogauss;
      by   *= Munits::kilogauss;
      bz   *= Munits::kilogauss;
      bmag *= Munits::kilogauss;
      // positions read in are given in cm
      x    *= Munits::cm;
      y    *= Munits::cm;
      z    *= Munits::cm;

      // all the maps aren't in consistent order (ascending/descending x,y)
      // so we have to actually interpret the (x,y) values to determine 
      // the array index.
      // we'll put them in the array in a way that we can lookup 
      // based on fX0,fDx,fNx,fY0,fDy,fNy

      indx = LowerLeftGenerator(x,y);

      // create new TVector3 and put it in the TClonesArray
      // what a crazy syntax!
      new((*fFieldVectors)[indx]) TVector3(bx,by,bz);

#ifdef BLAHBLAH
      // one *might* think of using [i] rather than [indx]
      if (indx != i &&  i < 2*fNx ) {
         MSG("Bfld",Msg::kWarning) << " i " << i << " indx " << indx
                                   << " (" << x << "," << y << ")" << endl;
      }
#endif


   }

   MSG("Bfld",Msg::kInfo) << "BfldMapRect2d file read in" << endl;

   fReadOk = kTRUE;
}

//______________________________________________________________________________
BfldMapRect2d::~BfldMapRect2d()
{
   // dtor deletes owned components

   MSG("Bfld",Msg::kDebug) << "~BfldMapRect2d variant " << fVariant << endl;

   fGrid    = BfldGrid::kUndefined;
   fVariant = -1;

   if (fFieldVectors) fFieldVectors->Delete();

   delete fFieldVectors;
   fFieldVectors = 0;

}

//______________________________________________________________________________
TVector3 BfldMapRect2d::GetBField(Int_t generator)
{
   TVector3 bvector = GetRawField(generator);
   // apply corrections for coil, chemistry, end plane-ness...
   return bvector;
}

//______________________________________________________________________________
TVector3 BfldMapRect2d::GetRawField(Int_t generator)
{
   if (fReadOk && generator > 0 && generator <= fFieldVectors->GetLast()) {
      TVector3* raw = dynamic_cast<TVector3*>(fFieldVectors->At(generator));
      if (raw) return *raw;
      MSG("Bfld",Msg::kError)
         << " no raw field for generator " << generator << endl;
   }
   return TVector3(0.,0.,0.);
}

//______________________________________________________________________________
void BfldMapRect2d::SetMesh(BfldMesh *mesh)
{
   fMesh = mesh;
   BfldMeshRect2d* rect2dMesh = dynamic_cast<BfldMeshRect2d*>(mesh);
   if ( rect2dMesh ) {
      rect2dMesh->SetMapRect2d(this);
   } else {
      MSG("Bfld",Msg::kWarning) << "BfldMapRect2d::SetMesh wasn't passed a Rect2 object" << endl;
   }
}

//______________________________________________________________________________
Int_t BfldMapRect2d::IndicesToGenerator(Int_t ix, Int_t iy)
{
   // return the generator (indx) based on x-y indices

   return ix*fNy + iy;
}

//______________________________________________________________________________
void BfldMapRect2d::GeneratorToIndices(Int_t gen, Int_t &ix, Int_t &iy)
{
   // return the x-y indices base on generator number

   // composed via ix*fNy + iy;

   ix = gen/fNy;
   iy = gen%fNy;

}

//______________________________________________________________________________
Int_t BfldMapRect2d::LowerLeftGenerator(Float_t x, Float_t y)
{
   // return the generator (indx) associated with the (x,y) position

   if ( fQuadrant ) {
      // symmetry case -- assume quadrant symmetry
      x = TMath::Abs(x);
      y = TMath::Abs(y);
   }

   // add epsilon to deal with roundoff errors introduced in
   // converting from cm in map file to Munits base length
   const Float_t epsilon = 1.0e-5;

   Float_t xx = (x - fX0) / fDx + epsilon;
   Float_t yy = (y - fY0) / fDy + epsilon;
   Int_t ix = (int)xx;
   Int_t iy = (int)yy;

   ix = TMath::Min(ix,fNx-1);
   iy = TMath::Min(iy,fNy-1);

   ix = TMath::Max(ix,0);
   iy = TMath::Max(iy,0);

   return IndicesToGenerator(ix,iy);
}

//______________________________________________________________________________
Int_t BfldMapRect2d::LowerLeftGenerator(TVector3 xyz)
{
   // return the generator (indx) associated with the (x,y) position

   return LowerLeftGenerator(xyz.X(),xyz.Y());
}

//______________________________________________________________________________
Int_t BfldMapRect2d::NeighborGenerator(Int_t gen, Int_t dnx, Int_t dny)
{
   // return the generator (indx) which is a neighbor to the input one
   // return -1 if no such index

   Int_t ix = 0;
   Int_t iy = 0;

   GeneratorToIndices(gen,ix,iy);

   ix = ix + dnx;
   iy = iy + dny;

   if (ix<0 || ix>=fNx) return -1;
   if (iy<0 || iy>=fNy) return -1;

   return IndicesToGenerator(ix,iy);

}

//______________________________________________________________________________

---