BField.cxx

Go to the documentation of this file.

// $Id: BField.cxx,v 1.39 2009/09/16 16:55:44 rhatcher Exp $
//
// BField
//
// Implementation code for BField
//
// Author:  R. Hatcher 2000.06.20
//
#include "BField/BField.h"

#include "BField/BfldLoanPool.h"
#include "BField/BfldHandlerRect2d.h"
#include "BField/BfldHandlerVoronoi.h"
#include "BField/BfldCache.h"
#include "BField/BfldMap.h"
#include "DcsUser/CoilTools.h"

#include "TMath.h"

#include "MessageService/MsgService.h"
CVSID("$Id: BField.cxx,v 1.39 2009/09/16 16:55:44 rhatcher Exp $");

ClassImp(BField)

#include <cassert>
#include <iomanip>

const Int_t kNoMapUsed = -9999;
static const TVector3 zeroBfld(0.,0.,0.);

#if 0
static std::ostream& operator<<(std::ostream& os, const TVector3& tv3)
{
  int w=9, p=4; // 1mm precision
  int prec = os.precision();
  os << "[" << setw(w)   << setprecision(p) << tv3.x()
     << "," << setw(w)   << setprecision(p) << tv3.y()
     << "," << setw(w+1) << setprecision(p) << tv3.z()
     << "]";
   os << resetiosflags(ios::floatfield);  // reset to "%g" format
   os << setprecision(prec);  // restore precision
  return os;
}
#endif

//_____________________________________________________________________________
BField::BField()
  : fLoanPool(0), fHandler(0), fCache(0), 
    fCurrentPlaneMap(0), fOwnedPlaneMap(0)

{ 
   // Default constructor -- needed if classes are to inherit this interface
   Init(); 
   InitFlagsFromLoanPool();
   fNoField        = true;
}

//_____________________________________________________________________________
BField::BField(VldContext vldc, Int_t coarseness, Int_t useEverywhere)
  : fNoField(false), fVldContext(vldc), fLoanPool(0), fHandler(0), fCache(0),
    fPowerOff(false), fDegaussed(false), fCoilCurrent(0),
    fCurrentPlaneMap(0), fOwnedPlaneMap(0)
{
   // constructor
   //
   // "useEverywhere" is a temporary ruse for bypassing the need
   // for geometry during development, it assumes that a old 2d rectangular
   // grid is relevant *everywhere* and this is the ID.
   //
   // coarseness should be set to -1 to use a rectangular grid, in which
   // case useEverywhere should specify which of the old IDs to use.

   Init();
   InitFlagsFromLoanPool();
   fUseEverywhere = useEverywhere;

   // set VldRange
   VldTimeStamp start = VldTimeStamp::GetBOT();
   VldTimeStamp end   = VldTimeStamp::GetEOT();
   fVldRange = VldRange(vldc.GetDetector(),
                        vldc.GetSimFlag(),
                        start,
                        end,
                        "fake");

   if (vldc.GetDetector() == Detector::kCalDet) {
     fNoField = true;
     return;
   }

   // set the grid and attach appropriate handler
   // do this before we decide about UseEverywhere
   // (which might screw up the VldContext)
   SetGridHandler(vldc.GetDetector(),coarseness);   

   // have we been told by the package config to always
   // use a particular map?
   int forceEverywhere = fLoanPool->GetForceUseEverywhere();
   if (forceEverywhere) {
     if ( forceEverywhere == 1 ) {
       // special case 1 is not a particular map but to use
       // the code's choice of "best for this VldContext".
       forceEverywhere = BfldCache::GetDefaultMapVariant(fVldContext);
       MAXMSG("Bfld",Msg::kInfo,5)
         << "Bfield::ctor BfldLoanPool forceUseEverywhere=1"
         << " yields map " << forceEverywhere << " for "
         << endl << "    " << vldc.AsString("c") << "."
         << endl;
     }
     MAXMSG("Bfld",Msg::kInfo,5)
       << "BField::ctor useEverywhere was " << fUseEverywhere
       << ", forced by BfldLoanPool config to use "
       << forceEverywhere << "."
       << endl;
     fUseEverywhere = forceEverywhere;
   }

   // if we've been told to use a particular map (either from the
   // ctor args or from the package config), have we been
   // told by the package config to ignore this?
   if (fUseEverywhere) {
     int ignore = fLoanPool->GetIgnoreUseEverywhere();
     MAXMSG("Bfld",Msg::kInfo,5) 
       << "BField::ctor useEverywhere=" << fUseEverywhere
       << ((ignore)?" (now ignored)":" (active)") << "."
       << endl;
     if (ignore) fUseEverywhere = 0;
     else {
       // use fOwnedPlaneMap for keeping "everwhere" info
       fOwnedPlaneMap->SetMapVariant(0,fUseEverywhere);
       fOwnedPlaneMap->SetScale(0,1.0);
       fVldContext = VldContext();  // !NOTE! mucking up VldContext!
       MAXMSG("Bfld",Msg::kInfo,5) 
         << "BField::ctor setting fake VldContext to avoid "
         << "pulling geometry tables from DBI."
         << endl;
     }
   }

   // attach an appropriate cache we get from a loan pool
   // also look up coil current info
   ResetVldContext(vldc);

}

//_____________________________________________________________________________
BField::BField(const BField& that)
  : TObject(that)
{
  // (deep) copy ctor accounts for keeping ref counts straight
  *this = that; 
}

//_____________________________________________________________________________
BField& BField::operator=(const BField& that)
{
   // assignment operator (deep)

   if ( this == &that ) return *this;  // do nothing if asked to self-assign
   TObject::operator=(that);

   // delete owned handler, decrement cache reference
   fLoanPool = 0;
   if (fHandler) { delete fHandler;        fHandler = 0; }
   if (fCache)   { fCache->DecrementRef(); fCache = 0; } // not owned, can't delete it

   Init();

   // copy data
   fNoField             = that.fNoField;
   fVldContext          = that.fVldContext;
   fVldRange            = that.fVldRange;
   fLastMapVariant      = that.fLastMapVariant;

   fDoBHCorrection      = that.fDoBHCorrection;
   fDoSlotCorrection    = that.fDoSlotCorrection;
   fDoInterPlaneField   = that.fDoInterPlaneField;
   fDoSMGapAndEndField  = that.fDoSMGapAndEndField;

   fUseDCSCoilDir       = that.fUseDCSCoilDir;
   fUseDCSCurrent       = that.fUseDCSCurrent;
   fGapLineSrc          = that.fGapLineSrc;
   fApplyBdotScale      = that.fApplyBdotScale;

   fUseEverywhere       = that.fUseEverywhere; 
   fNoFieldBeyondZ      = that.fNoFieldBeyondZ;

   // attach the same cache as the copy-ee
   if ( that.fCache ) {
     fCache = that.fCache;
     fCache->IncrementRef();
   }
   else {
     // hey! this shouldn't happen
     assert(0);
   } 

   // translate chosen grid back into coarseness
   Int_t coarseness = BfldGrid::GetCoarseness(that.fGrid);

   // set the grid and attach appropriate handler
   SetGridHandler(fVldContext.GetDetector(),coarseness);

   fHandler->SetCache(fCache);
   fHandler->SetLastPolygAsSeed();

   return *this;
}

//_____________________________________________________________________________
void BField::ResetVldContext(const VldContext& vldc)
{
   // for any given VldContext we *should* look up coil status + current
 
   fVldContext = vldc;
   Detector::Detector_t det = vldc.GetDetector();

   VldTimeStamp start = VldTimeStamp::GetBOT();
   VldTimeStamp end   = VldTimeStamp::GetEOT();
   fVldRange = VldRange(det,vldc.GetSimFlag(),start,end,"");

   if (!fHandler) SetGridHandler(vldc.GetDetector(),-1);

   if (fCache) fCache->DecrementRef(); // unreference previous cache, if any
   
   // attach a cache we get from a loan pool
   fCache = fLoanPool->GetCache(fVldContext);
   fCache->IncrementRef();  // let cache know we have a pointer to it
   fVldRange.TrimTo(fCache->GetVldRange());

   fHandler->SetCache(fCache);
   fHandler->SetLastPolygAsSeed();

   fNoField   = false;
   fPowerOff  = false;
   fDegaussed = false;

   // here we determine the coil current and whether the detector
   // is degaussed (if I=0).  Do it early as special cases 
   // (eg. collar/coil or interplane ) might need it.

   // rwh: for now we just say get ready to multiply by 1.0
   switch (fVldContext.GetDetector()) {
   // nominal values in Amp-turns
   case Detector::kNear: 
     fCoilCurrent = 40.0 * 1000.;  // nominal "on" value
     if ( fUseDCSCurrent ) {
       std::pair<Float_t,Float_t> currents = CoilTools::CoilCurrent(vldc);
       fCoilCurrent = 8.0 * TMath::Abs(currents.first); // 8 planks * ~ -4973A
     }
     break;
   case Detector::kFar:
     fCoilCurrent = 15.2 * 1000.;  // nominal "on" value
     if ( fUseDCSCurrent ) {
       std::pair<Float_t,Float_t> currents = CoilTools::CoilCurrent(vldc);
       Float_t avg_cur = 0.5 * (currents.first+currents.second);
       fCoilCurrent = 190 * TMath::Abs(avg_cur); // 190 turns * ~ 80A
     }
     break;
   default:
     fCoilCurrent = 0.0;
     fNoField = true;
   }
   // should scale fCoilCurrent by ratio of readback to nominal
   // from DCS readback for non-saturated part of field ...

   if (fUseDCSCoilDir == 1 || fUseDCSCoilDir == 2 ) {
     if ( ! fUseDCSCurrent && ! CoilTools::IsOK(vldc)) fCoilCurrent = 0.0;
     else if (CoilTools::IsReverse(vldc)) fCoilCurrent = -fCoilCurrent;
   }
}

//_____________________________________________________________________________
void BField::Init()
{
   // initialize member data
   fNoField        = false;
   fVldRange       = VldRange();
   fGrid           = BfldGrid::kUndefined;
   fLastMapVariant = kNoMapUsed;
   fUseEverywhere  = 0; 

   // get access to the BfldLoanPool
   if (!fLoanPool) fLoanPool = BfldLoanPool::Instance();

   if (fOwnedPlaneMap) { delete fOwnedPlaneMap;  fOwnedPlaneMap = 0; }
   fOwnedPlaneMap = new BfldDbiPlaneMap;
}

//_____________________________________________________________________________
void BField::InitFlagsFromLoanPool()
{
   // initialize member data from loan pool

   // get access to the BfldLoanPool
   if (!fLoanPool) fLoanPool = BfldLoanPool::Instance();

   fDoBHCorrection      = fLoanPool->GetDefaultDoBHCorrection();
   fDoSlotCorrection    = fLoanPool->GetDefaultDoSlotCorrection();
   fDoInterPlaneField   = fLoanPool->GetDefaultDoInterPlaneField();
   fDoSMGapAndEndField  = fLoanPool->GetDefaultDoSMGapAndEndField();

   fUseDCSCoilDir       = fLoanPool->GetDefaultUseDCSCoilDir();
   fUseDCSCurrent       = fLoanPool->GetDefaultUseDCSCurrent();
   fGapLineSrc          = fLoanPool->GetDefaultGapLineSrc();
   fApplyBdotScale      = fLoanPool->GetDefaultApplyBdotScale();

   fNoFieldBeyondZ      = fLoanPool->GetDefaultNoFieldBeyondZ();
}

//_____________________________________________________________________________
BField::~BField()
{
   // destructor

   if (fHandler)       { delete fHandler;        fHandler = 0; }
   // fCache is not owned, can't delete it just change ref count
   if (fCache)         { fCache->DecrementRef(); fCache = 0; } 
   if (fOwnedPlaneMap) { delete fOwnedPlaneMap;  fOwnedPlaneMap = 0; }

}

//_____________________________________________________________________________
void BField::SetInterpMethod(BfldInterpMethod::InterpMethod_t method)
{
   // pass on interpolation method message to BFLHandler
   if (fHandler) fHandler->SetInterpMethod(method);
}

//_____________________________________________________________________________
void BField::SetGridHandler(Detector::Detector_t detector, Int_t coarseness)
{
   // translate Detector and "coarseness" into Grid enum
   fGrid = BfldGrid::GetGrid(detector,coarseness);
   
   // attach a BfldHandler (delete old if there is one)
   if (fHandler) { delete fHandler; fHandler = 0; }

   // check for special case of "rect2d" grid request
   if ( BfldGrid::kRect2dGrid == fGrid ) {
      // attach a special handler
      fHandler = new BfldHandlerRect2d;
   } else {
      fHandler = new BfldHandlerVoronoi;
   }

   // set the default interpolation method
   // BField's method simply passes on message to BFLHandler
   SetInterpMethod(BfldInterpMethod::kDefault);

   fLastMapVariant = kNoMapUsed;

}

//_____________________________________________________________________________
TVector3 BField::GetBField(TVector3& posGlobal, Bool_t isUVZ)
{
   // Return the magnetic field for any point in space
   // for the currently configured detector
   // if isUVZ position is in 

   if (fNoField) return zeroBfld;

   if (posGlobal.Z() > fNoFieldBeyondZ ) return zeroBfld;

   // stash this way for special cases ( coil, SM gap, interplane )
   fPositionIsUVZ = isUVZ;
   
   bool isPosSteelUVZ = isUVZ;
   int doLocalTransform = 0;

   if ( fDegaussed && fCoilCurrent == 0 ) return zeroBfld;

   bool was_in_steel = true;

   // calculate which plate (if any) this "z" corresponds to
   if (fUseEverywhere) {

     MAXMSG("Bfld",Msg::kWarning,4) 
       << "BField::GetBField no coil region handling or " << endl
       << "map selection based on plane with useEverywhere=" 
       << fUseEverywhere << "."
       << endl;
     fCurrentPlaneMap = fOwnedPlaneMap;

   } else {

      // check for special cases of geometry (via cache)
      //  0 =: don't check, outside steel use either 
      //       none (no InterPlaneField) or interplane gap map
      // <0 =: check position, but use zero when up/down/gap
      // >0 =: check position, try best to model field
      if ( fDoSMGapAndEndField ) {
        Ugli::SMRegion_t ismregion = fCache->InSMRegion(posGlobal,isUVZ);
        if ( ismregion != Ugli::kInSM1 && ismregion != Ugli::kInSM2 ) {
          // in SM gap or beyond SM ends
          if ( fDoSMGapAndEndField < 0 ) return zeroBfld; // no attempt to model
          return SMGapAndEndField(posGlobal,ismregion);
        }
      }

      // use geometry (via cache) to determine plane
      fCurrentPlaneMap = fCache->FindPlaneMap(posGlobal);
      was_in_steel = fCache->FindWasInSteel(true);
      // check that there was data

      if ( ! fCurrentPlaneMap ) {
        // nope?  probably because lacking BfldDbiPlaneMap table
        // fake something up to mimic old behaviour
        Int_t    imap  = BfldCache::GetDefaultMapVariant(fVldContext);
        Double_t scale = BfldCache::GetDefaultScale(fVldContext);
        MAXMSG("Bfld",Msg::kWarning,5)
          << "GetDefault(MapVariant,Scale) returned (" 
          << imap << "," << scale << ")"
          << ", called because missing BfldDbiPlaneMap?"
          << endl;
        fOwnedPlaneMap->SetMapVariant(0,imap);
        fOwnedPlaneMap->SetScale(0,scale);   
        fCurrentPlaneMap = fOwnedPlaneMap;
      }
   }

   // by default use the full plane map pair
   Int_t indxpair = BfldDbiPlaneMap::kFullSteelA;

   TVector3 posSteel(posGlobal);
   doLocalTransform = fCache->GetDoLocalTransform();
   if (doLocalTransform > 0) {
     posSteel = fCache->GetPositionInSteel();
     //if (doLocalTransform > 1) isPosSteelUVZ = false;
   }
   bool indetail = fCurrentPlaneMap->IsInDetail(posSteel.x(),posSteel.y());

   if ( ! was_in_steel ) {
     //
     // what to do if not (considered) in the steel
     //
     // if asked to do nothing ... return zero field
     if ( 0 == fDoInterPlaneField ) return zeroBfld;
     indxpair = BfldDbiPlaneMap::kDetailGapA;  // default to gap detail map
     if ( ! indetail ) {
       // not in detail region?
       if ( 3 == fDoInterPlaneField ) return zeroBfld;
       indxpair = BfldDbiPlaneMap::kFullGapA;
       // fall back to line current if no full gap map exists
       if (fCurrentPlaneMap->IsPairNull(indxpair) || 
           2 == fDoInterPlaneField) {
         TVector3 blinesrc(0,0,0);
         if ( fGapLineSrc & 1 ) {
           TVector3 posCoil = posSteel; // for now!!
           blinesrc += BFromLineSource(posCoil,fCoilCurrent);
         }
         if ( fGapLineSrc & 2 ) {
           TVector3 posRelRtnCoil = posGlobal;
           // account for position of return coil
           switch ( fVldContext.GetDetector() ) {
           case Detector::kNear:
             posRelRtnCoil -= TVector3(-1.9,-1.9,0);
             break;
           case Detector::kFar:
             posRelRtnCoil -= TVector3( 0.0,-4.8,0);
             break;
           default:
             posRelRtnCoil = TVector3(999.,999.,999.);
             break;
           }
           // note negative sign for fCoilCurrent for return leg
           blinesrc += BFromLineSource(posRelRtnCoil,-fCoilCurrent);
         }
         return blinesrc;
       }
     }
   } else {
     //
     // what to do if (considered) in the steel
     // set the right map pair for the conditions and position
     //
     if ( fPowerOff ) indxpair = BfldDbiPlaneMap::kPowerOffA;
     else if ( indetail ) {
       indxpair = BfldDbiPlaneMap::kDetailSteelA;
       if (fCurrentPlaneMap->IsPairNull(indxpair)) 
         indxpair = BfldDbiPlaneMap::kFullSteelA; // no detail map, revert
     }
   }

   Bool_t isnull[2];  // check if scale or mapVariant is 0
   isnull[0] = fCurrentPlaneMap->IsNull(indxpair);
   isnull[1] = fCurrentPlaneMap->IsNull(indxpair+1);

   Double_t coilCurrMap[2] = { 0, 0 };

   if ( isnull[0] && isnull[1] ) return zeroBfld;  // I've got nut'n

   TVector3 bfieldSum(zeroBfld);  // start with zero field

   for (UInt_t iv = 0; iv<2; ++iv) {  // loop over 2 possible contributions

     if ( isnull[iv] ) continue; // no contribution
     
     // variant refers to (in part) the steel chemistry
     // end-plate info to get right "variant"
     Int_t mapVariant  = fCurrentPlaneMap->GetMapVariant(indxpair+iv);
     Double_t scaleMap = fCurrentPlaneMap->GetScale(indxpair+iv);
     // if using DCS for sign, then ignore sign associated with scale
     if ( (fUseDCSCoilDir&1) ) scaleMap = TMath::Abs(scaleMap);

     // this should never happen
     if (mapVariant == 0) {
       MAXMSG("Bfld",Msg::kInfo,5)
         << "BField mapVariant=0 indxpair=" << indxpair << " iv=" << iv
         << " isnull=" << (isnull[iv]?"true":"false")
         << " scaleMap=" << scaleMap
         << endl;
       MAXMSG("Bfld",Msg::kInfo,5) << *fCurrentPlaneMap << endl;
       continue;
     }

     BfldMap* bmap = SetupHandlerForMap(mapVariant);
 
     coilCurrMap[iv] = bmap->GetGeneratedCoilCurrent();
     Double_t scaleCoil = fCoilCurrent/coilCurrMap[iv]; 
     // probably this effect isn't linear when field is saturated...
     Double_t scale = scaleCoil * scaleMap;
     //MAXMSG("Bfld",Msg::kDebug,5)
     //  << "Bfld scaleCoil " << scaleCoil << " CoilCurrent " << fCoilCurrent
     //  << " MapCurrent " << coilCurrMap[iv] << endl;

     // The handler is now configured ... interogate it for the field
     bfieldSum += (scale*fHandler->GetBField(posSteel,isPosSteelUVZ));

   } 

   if ( was_in_steel && doLocalTransform > 1 ) {
     // get the steel, to tranform B_local to B_global
     // note Ugli uses globalInXYZ, while BField uses isUVZ
     // always pass true here ... don't double rotate!
     UgliSteelPlnHandle usph = fCache->GetCurrentSteelPlnHandle();
     bfieldSum = usph.LocalToGlobalVect(bfieldSum,true);
   }

   // should these be done by the handler on the pre-(scaled&summed)
   // version before summing? i.e. which is more correct:
   //  (a)  bhcorr(slotcor(map0))*scale0 + bhcorr(slotcor(map1))*scale1
   //  (b)  bhcorr(slotcor(map0*scale0 + map1*scale1))
   // this choice is (b) otherwise we should move these to the 
   // BfldHandler and make it aware of the flags and initialize
   // it with the fCurrentPlaneMap so it can pick up the parameters

   if      ( fDoSlotCorrection   ) ApplySlotFactorCorr(bfieldSum);
   if      ( fDoBHCorrection < 0 ) ApplyBHFactorCorr(bfieldSum);
   else if ( fDoBHCorrection > 0 ) ApplyBHCurveCorr(bfieldSum);

   if ( fApplyBdotScale ) {
     // don't apply to the gap cases
     if ( indxpair < BfldDbiPlaneMap::kDetailGapA ) {
       Double_t bdot = fCurrentPlaneMap->GetBdotScale();
       // don't apply if for some reason bdot value is == 0 (bad DB value)
       if (bdot == 0.0) {
         MAXMSG("Bfld",Msg::kWarning,5)
           << " BdotScale was zero for plane " << fCurrentPlaneMap->GetPlaneId()
           << "." << endl;
       } else {
         bfieldSum *= bdot;
       }
     }
   }

   return bfieldSum;
}
                                
//_____________________________________________________________________________
BfldMap* BField::SetupHandlerForMap(Int_t mapVariant)
{
  // Ensure that the BfldHandler is properly configured
  // 
  if ( fLastMapVariant != mapVariant ) {
    //
    // map and/or grid has changed ... reinitialize handler
    //

    // Get the right base field map
    BfldMap  *bmap = fLoanPool->GetMap(fGrid,mapVariant);
    fHandler->SetMap(bmap);
    
    // Get the right Mesh diagram from the LoanPool
    // Give our BfldHandler a pointer to the v-diagram
    BfldMesh *mesh  = fLoanPool->GetMesh(fGrid,mapVariant);
    fHandler->SetMesh(mesh);
    
    
    fLastMapVariant = mapVariant;
    
  }  // map variant configured
  
  return fHandler->GetMap();
}

//_____________________________________________________________________________
void BField::ApplySlotFactorCorr(TVector3& b)
{
  //
  // Slot (gap) correction formula:
  //
  //
  //   B'  =   Bo                                 for Bo > Bcut
  //           Bo [ 1 - D * ( 1 - Bo/Bcut )^2 ]   
  //

#ifdef DEBUG
  MAXMSG("Bfld",Msg::kInfo,20) 
    << "ApplySlotFactorCorr " 
    << " SlotFactor " << fCurrentPlaneMap->GetSlotFactor()
    << " SlotCutoff " << fCurrentPlaneMap->GetSlotCutoff()
    << " mag " << b.Mag()
    << endl;
#endif

  Double_t slotfactor = fCurrentPlaneMap->GetSlotFactor();
  if ( slotfactor == 0 ) return;
  
  Double_t slotcutoff = fCurrentPlaneMap->GetSlotCutoff();
  if ( slotcutoff == 0 ) return;
     
  Double_t magbfld    = b.Mag();
  if ( magbfld > slotcutoff ) return;
     
  Double_t inner = 1.0 - magbfld/slotcutoff;
  Double_t slotscale = 1.0 - slotfactor * inner*inner;

#ifdef DEBUG
  MAXMSG("Bfld",Msg::kInfo,20) 
    << "ApplySlotFactorCorr " << magbfld << " " << slotscale << endl;
#endif  

  // apply correction
  b *= slotscale;
}

//_____________________________________________________________________________
void BField::ApplyBHFactorCorr(TVector3& b)
{
  // 
  // Steel chemistry correction (BH)
  //
  //           Bo                     for Bo > Bcut
  //   B'  =   Bo [ 1 + G*(1-Bo/C) ]  for Bcut/3 < Bo < Bcut
  //           Bo [ 1 + 2/3*G*Bo/C ]  for Bo < Bcut/3
  //

#ifdef DEBUG
  MAXMSG("Bfld",Msg::kInfo,20) 
    << "ApplyBHFactorCorr " 
    << " BHFactor " << fCurrentPlaneMap->GetBHFactor()
    << " BHCutoff " << fCurrentPlaneMap->GetBHCutoff()
    << " mag " << b.Mag()
    << endl;
#endif

  Double_t bhfactor = fCurrentPlaneMap->GetBHFactor();
  if ( bhfactor == 0 ) return;

  Double_t bhcutoff = fCurrentPlaneMap->GetBHCutoff();
  if ( bhcutoff == 0 ) return;

  Double_t magbfld  = b.Mag();
  if ( magbfld > bhcutoff) return;

  Double_t bhscale = 1.0;
  if ( magbfld > bhcutoff/3.0 ) {
    bhscale = 1.0 + bhfactor * ( 1.0 - magbfld/bhcutoff );
  }
  else {
    bhscale = 1.0 + (2.0/3.0)*bhfactor*magbfld/bhcutoff;
  }

#ifdef DEBUG
  MAXMSG("Bfld",Msg::kInfo,20) 
    << "ApplyBhFactorCorr " << bhscale << endl;
#endif 

  // apply correction
  b *= bhscale;
}

//_____________________________________________________________________________
void BField::ApplyBHCurveCorr(TVector3& /* b */ )
{
  // Use B-H curve correction
  MAXMSG("Bfld",Msg::kWarning,20)
    << "BField::ApplyBHCurveCorr not yet implemented!"
    << endl;
  // do nothing!
}

//_____________________________________________________________________________
TVector3 BField::BFromLineSource(TVector3& posRelCoil, Double_t current)
{
  // Field in air from a line source along z axis

  static double lastcurrent = 9999;
  if (lastcurrent != current) {
    MAXMSG("Bfld",Msg::kInfo,5)
      << "BFromLineSource current is " << current << endl;
    lastcurrent = current;
  }

  // B = u0*I/2pi*r = 0.003/r
  const Double_t u0by2pi = 12.566370614e-7 / (2.0 * TMath::Pi());
  
  Double_t r = posRelCoil.Perp();
  // protect against divide by zero
  r = TMath::Max(r,Munits::cm);  // no line src smaller than 1 cm
  Double_t brbyr = ( u0by2pi * current / r )/r;
  Double_t bx    =  posRelCoil.Y() * brbyr;
  Double_t by    = -posRelCoil.X() * brbyr;
  return TVector3(bx,by,0.0);

}

//_____________________________________________________________________________
TVector3 BField::SMGapAndEndFieldNear(const TVector3& position,
                                      Ugli::SMRegion_t iregion)
{
  // Return the field in the region up-/down- stream of detector
  // don't forget to check fIsUVZ

  int w=9, p=4; // 1mm precision
  int prec = cout.precision(); // hopefully this is the right output stream 
  MAXMSG("Bfld",Msg::kWarning,20)
    << "BField::SMGapAndEndFieldNear "
    << endl
    << "not yet implemented for SMRegion " << Ugli::AsString(iregion) << " @ "
    << "[" << setw(w)   << setprecision(p) << position.x()
    << "," << setw(w)   << setprecision(p) << position.y()
    << "," << setw(w+1) << setprecision(p) << position.z()
    << "]"
    << resetiosflags(ios::floatfield)  // reset to "%g" format
    << setprecision(prec)              // restore precision
    << endl;
  return zeroBfld;

}

//_____________________________________________________________________________
TVector3 BField::SMGapAndEndFieldFar(const TVector3& position,
                                     Ugli::SMRegion_t iregion)
{
  // Return the field in the region up-/down- stream of detector
  // or in the gap between the supermodules
  // don't forget to check fIsUVZ

  int w=9, p=4; // 1mm precision
  int prec = cout.precision(); // hopefully this is the right output stream 
  MAXMSG("Bfld",Msg::kWarning,20)
    << "BField::SMGapAndEndFieldFar "
    << endl
    << "not yet implemented for SMRegion " << Ugli::AsString(iregion) << " @ "
    << "[" << setw(w)   << setprecision(p) << position.x()
    << "," << setw(w)   << setprecision(p) << position.y()
    << "," << setw(w+1) << setprecision(p) << position.z()
    << "]"
    << resetiosflags(ios::floatfield)  // reset to "%g" format
    << setprecision(prec)              // restore precision
    << endl;
  return zeroBfld;

}

//_____________________________________________________________________________

Int_t BField::GetDoLocalTransform() const
{  return fCache->GetDoLocalTransform(); }

Int_t BField::GetRequireInZTest() const
{  return fCache->GetRequireInZTest(); }

Double_t BField::GetZTolerance() const
{ return fCache->GetZTolerance(); }

void BField::SetDoLocalTransform(Int_t iflg)
{ 
  // Just pass along the value
  // NOTE: this can affect other BField objects that share the same cache!
  fCache->SetDoLocalTransform(iflg); 
} 

void BField::SetRequireInZTest(Int_t ival)
{ 
  // Just pass along the value
  // NOTE: this can affect other BField objects that share the same cache!
  fCache->SetRequireInZTest(ival); 
} 

void  BField::SetZTolerance(Double_t zeps)
{
  // Just pass along the value
  // NOTE: this can affect other BField objects that share the same cache!
  fCache->SetZTolerance(zeps);
}


//_____________________________________________________________________________
void BField::Print(Option_t * /* option */) const
{
   // output info about this BField object

   MSG("Bfld",Msg::kInfo) 
     << "BField: created from VldContext " << fVldContext
     << endl << " VldRange: " << fVldRange << endl;
   MSG("Bfld",Msg::kInfo)
     << " Co-ords: " << (fPositionIsUVZ?"UVZ":"XYZ")
     << ", NoField: " << (fNoField?"yes":"no")
     << endl
     << " Power: " << (fPowerOff?"off":"on")
     << ", Degaussed: " << (fDegaussed?"yes":"no")
     << ", Current: " << fCoilCurrent
     << endl
     << " UseDCSCoilDir: " << fUseDCSCoilDir
     << ", UseDCSCurrent: " << fUseDCSCurrent
     << ", GapLineSrc: " << fGapLineSrc
     << ", InterPlaneField: " << fDoInterPlaneField
     << ", SMGapAndEndField: " << fDoSMGapAndEndField
     << ", ApplyBdotScale: " << (fApplyBdotScale?"yes":"no")
     << endl
     << " BHCorr: " << fDoBHCorrection
     << ", SlotCorr: " << fDoSlotCorrection
     << endl;
   if (fGrid) 
     MSG("Bfld",Msg::kInfo) 
       << " Grid: " << BfldGrid::AsString(fGrid) << endl;
   if (fCache)   {
     MSG("Bfld",Msg::kInfo) << " BfldCache: ";
     fCache->Print();
   }
   if (fHandler) {
     MSG("Bfld",Msg::kInfo) << " Handler: ";
     fHandler->Print();
   }
}

//_____________________________________________________________________________

---