BfldMapRect2d.cxx

Go to the documentation of this file.

// $Id: BfldMapRect2d.cxx,v 1.26 2006/06/23 19:27:07 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"
#include "MessageService/MsgFormat.h"
CVSID("$Id: BfldMapRect2d.cxx,v 1.26 2006/06/23 19:27:07 rhatcher Exp $");

#include <cassert>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <sstream>

#include <vector>

ClassImp(BfldMapRect2d)

//_____________________________________________________________________________
BfldMapRect2d::BfldMapRect2d() 
   : BfldMap()
{
  this->Clear();
}

//_____________________________________________________________________________
BfldMapRect2d::BfldMapRect2d(BfldGrid::Grid_t grid, Int_t variant) 
   : BfldMap(grid,variant)
{

   this->Clear();

   // read in a Map

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

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

   char string[255];

   // format and expand file name
   const char* search_path = ".:$BMAPPATH";
   sprintf(string,"bfld_%3d.dat",varabs);
   char* expanded = gSystem->Which(search_path,string,kReadPermission);
   if (!expanded) {
      MSG("Bfld",Msg::kError) << "can not find input file: '" 
                              << fDataSource <<"'"<< endl;
      return;
   }
   fDataSource = expanded;

   delete [] expanded;  // user responsible for deleting returned string

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

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

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

   fIsUVZ = ( fHeader.find("UVZ") != string::npos );
   Double_t dirScale  = 1;

   // search the title string for a coil current
   // generally a token like 15kAturns, 40kAturns
   string::size_type loc_aturn = fHeader.find("Aturn");
   // alas sometimes the header has the value string "40k A-turns"
   if (loc_aturn == string::npos) loc_aturn = fHeader.find("A-turn");

   if (loc_aturn != string::npos) {
     // find last space before 'Aturn'
     string::size_type loc_space = fHeader.find_last_of(" ",loc_aturn); 
     std::string sval = fHeader.substr(loc_space,loc_aturn-loc_space);
     MSG("Bfld",Msg::kDebug) << " found Aturns: '" << sval << "'" << endl;
     string::size_type loc_k = sval.find("k");
     if (loc_k != string::npos) {
       sval.replace(loc_k,1," ");
       MSG("Bfld",Msg::kDebug) << " found 'k' new string '" << sval << "'" << endl;
     }
     std::istringstream inputString(sval);
     inputString >> fGeneratedCoilCurrent;
     if (loc_k != string::npos) fGeneratedCoilCurrent *= 1000;

     if ( fHeader.find("CurrentForward") == string::npos ) {
       MSG("Bfld",Msg::kInfo) 
         << "BfldMapRect2d: maps were created focusing mu+;"
         << " correct this with an extra -1" << endl;
       dirScale = -1;
     }   

     MSG("Bfld",Msg::kInfo) 
       << "BfldMapRect2d: GeneratedCoilCurrent=" 
       << setprecision(4) << (fGeneratedCoilCurrent/1000.)
       << "kA*turns based on header info."
       << resetiosflags(ios::floatfield)  // reset to "%g" format
       << endl;
   }
   else {
     MSG("Bfld",Msg::kError) 
       << "BfldMapRect2d: could not find coil current in input file header "
       << 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) 
     << "BfldMapRect2d:"
     << setprecision(6)
     << " x0 " << fX0 << " dx " << fDx << " nx " << fNx 
     << " y0 " << fY0 << " dy " << fDy << " ny " << fNy 
     << resetiosflags(ios::floatfield)  // reset to "%g" format
     << endl;

   // next 3 elements are irrelevant (z0,dz,nz)
   Float_t z0, dz;
   Int_t   nz;
   from >> z0 >> dz >> nz;

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

   // replace "x" with " "
   quads = quads.ReplaceAll("x"," ");
   // quads should be of the form ' 00 00 00 00  FNAL2'
   MSG("Bfld",Msg::kDebug) << "quads string '" 
                           << quads << "'" << endl;

   fFileForm = 0;
   if (quads.Contains("FNAL"))  fFileForm = 1;
   if (quads.Contains("FNAL2")) fFileForm = 2;

   // 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::kSynopsis) 
     << "BfldMapRect2d: fileform " << fFileForm
     << "  quadflags (" 
     << fQuadFlags[0] << "," 
     << fQuadFlags[1] << "," 
     << fQuadFlags[2] << "," 
     << fQuadFlags[3] << ")  use them " 
     << (fQuadrant ? "'true'" : "'false'") 
     << endl;

   ResizeVectors(fNx*fNy);

   // build a temporary mesh so we can can use the functions
   // for finding the generator from x,y etc
   // but this isn't the one that will go into the loan pool
   BfldMeshRect2d tmpMesh(this); 
   
   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.;

   std::vector<bool> seen_gen(fNx*fNy);

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

      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;
         // avoid unitialized values
         z    = 0.0;
         bmag = TMath::Sqrt(bx*bx+by*by+bz*bz);
         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;
      }

      if (!from.good()) {
        if (ilines != fNx*fNy)  //  read one too far ...
          MSG("Bfld",Msg::kInfo) 
            << "BfldMapRect2d: early EOF on " << endl 
            << "  "
            << fDataSource 
            << endl
            << "  saw " << ilines << " entries, expected " 
            << fNx*fNy << "; could be due to suppressed (0,0) points."
            << endl;
        break;
      }

      // The above puts field values in kGauss
      // convert them to base units
      // Also make correction for reversed current direction when generated
      bx   *= Munits::kilogauss * dirScale;
      by   *= Munits::kilogauss * dirScale;
      bz   *= Munits::kilogauss * dirScale;
      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

      // when reading points, find the nearest generator
      // not the lower left one
      indx = tmpMesh.NearestGenerator(x,y);

      if ( seen_gen[indx] ) {
        MAXMSG("Bfld",Msg::kWarning,20)
          << "Duplicate entry for indx " << indx
          << " x= " << x << ", y= " << y << endl;
      }
      else {
        seen_gen[indx] = true;
      }

      // cross check! for bad data points
      Int_t jxx, jyy;
      tmpMesh.GeneratorToIndices(indx,jxx,jyy);
      Float_t xx,yy;
      xx = fX0 + jxx*fDx;
      yy = fY0 + jyy*fDy;
      const Float_t dxy_epsilon = 0.0005;
      static MsgFormat f9p6("f9.6"); 
      static MsgFormat f5p1("f5.1"); 
      if ( TMath::Abs(x-xx) > dxy_epsilon*fDx ||
           TMath::Abs(y-yy) > dxy_epsilon*fDy    ) {
        MAXMSG("Bfld",Msg::kWarning,20)
          << " [" << setw(5) << ilines << "] off grid? "
          << endl 
          <<   "  (x,y)=(" << f9p6(x)  << "," << f9p6(y)  << ") "
          << " expected (" << f9p6(xx) << "," << f9p6(yy) << ") "
          << endl
          << "   error is " << f5p1(100.*TMath::Abs(x-xx)/fDx) << "% of dx "
          <<           ", " << f5p1(100.*TMath::Abs(y-yy)/fDy) << "% of dy "
          << endl;
      }

      AddGenerator(indx,bx,by,bz);

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

   }

   MSG("Bfld",Msg::kDebug) 
     << "BfldMapRect2d: read " << ilines << " entries from "
     << fDataSource << endl;

   fFilledOkay = kTRUE;
}

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

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

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

}

//_____________________________________________________________________________
void BfldMapRect2d::Clear(const Option_t*)
{
   // Clear to defined state.

   fX0           = 0.;
   fDx           = 0.;
   fNx           = 0;
   fY0           = 0.;
   fDy           = 0.;
   fNy           = 0;
   fQuadrant     = false;
   fQuadFlags[0] = 0;
   fQuadFlags[1] = 0;
   fQuadFlags[2] = 0;
   fQuadFlags[3] = 0;
   fFileForm     = 0;

}

//_____________________________________________________________________________

---