BfldCache.cxx

Go to the documentation of this file.

// $Id: BfldCache.cxx,v 1.30 2008/03/19 04:29:31 rhatcher Exp $
// 
// BfldCache 
// 
// BfldCache caches info for the BfldHandler
//
// Author:  R. Hatcher 2000.06.20
//

#include "BField/BfldCache.h"
#include "BField/BfldLoanPool.h"

#include "Plex/PlexPlaneId.h"
#include "Conventions/Munits.h"

#include "MessageService/MsgService.h"
CVSID("$Id: BfldCache.cxx,v 1.30 2008/03/19 04:29:31 rhatcher Exp $");

ClassImp(BfldCache)

#include <cassert>
#include <float.h>  // FLT_EPSILON for floating point comparison
#include <iomanip>

//_____________________________________________________________________________
// Utility functions

bool BFieldWithinZ(Double_t z, Double_t zmin, Double_t zmax)
{ 
  return ( z >= zmin && z <= zmax ); 
}
bool BFieldOutsideZ(Double_t z, Double_t zmin, Double_t zmax)
{ 
  return ( z < zmin || z > zmax ); 
}

#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

//_____________________________________________________________________________
BfldCache::BfldCache()
  : fRef(0), fSteelCache(), fSteelZmin(+9999.), fSteelZmax(-9999.),
    fDbiMapEntries(0), fPlaneMapCache(0)
{
   MSG("Bfld",Msg::kDebug) << "BfldCache - default ctor" << endl;
   InitFlags();
}

//_____________________________________________________________________________
BfldCache::BfldCache(const VldContext& vldc)
  : fDetector(vldc.GetDetector()), fRef(0), fUgliGeomHandle(vldc),
    fSteelCache(), fSteelZmin(+9999.), fSteelZmax(-9999.), 
    fDoLocalTransform(0), fInZTest(0), fZTolerance(0),
    fPlaneMaps(vldc,0,Dbi::kDisabled),
    fDbiMapEntries(fPlaneMaps.GetNumRows()),
    fPlaneMapCache(0), 
    fFindWasInSteel(false), fPosInSteel(-999.,-999.,-9999.),
    fWasHigh(0), fWasLow(0), fSumHigh(0), fSumLow(0)
{
   InitFlags();
   fVldRange = fUgliGeomHandle.GetVldRange();
   MSG("Bfld",Msg::kDebug) << " BfldCache vldc " << vldc 
                           << "  geometry " << fVldRange << endl;
   const DbiValidityRec* dbivrec = fPlaneMaps.GetValidityRec();
   if (dbivrec) {
     fVldRange.TrimTo(dbivrec->GetVldRange());
     MSG("Bfld",Msg::kDebug) << "BfldCache range trimed to " 
                             << fVldRange << endl;
   }

   if ( Detector::kUnknown == fDetector ) {
     // this happens for the case of UseEverywhere
     VldTimeStamp start = VldTimeStamp(1970, 1, 1, 0, 0, 0);
     VldTimeStamp end   = VldTimeStamp(2038, 1,18,19,14, 6);
     fVldRange = VldRange(fDetector,
                          vldc.GetSimFlag(),
                          start,
                          end,
                          "fake");
     MSG("Bfld",Msg::kDebug) << "BfldCache ctor with no detector "  << endl;
   }
   
   FillSMZLimits(); 
}

//_____________________________________________________________________________
void BfldCache::InitFlags()
{
   // Initialize configuration flags from BfldLoanPool

   // get access to the BfldLoanPool
   BfldLoanPool *loanpool = BfldLoanPool::Instance();

   fDoLocalTransform = loanpool->GetDefaultDoLocalTransform();
   fInZTest          = loanpool->GetDefaultRequireInZTest();
   fZTolerance       = loanpool->GetDefaultZTolerance();
   MSG("Bfld",Msg::kSynopsis)
     << "BfldCache DoLocalTransform=" << fDoLocalTransform 
     << " InZTest=" << fInZTest
     << " ZTolerance=" << (fZTolerance/Munits::mm) << "mm"
     << endl;

}

//_____________________________________________________________________________
BfldCache::~BfldCache()
{
  if (fRef) {
    MSG("Bfld",Msg::kInfo) << "~BfldCache with " << fRef 
                           << " outstanding references." << endl
                           << "     using geometry" 
                           << fUgliGeomHandle.GetVldRange() << endl;

  }
  MSG("Bfld",Msg::kDebug) << "~BfldCache geometry" 
                          << fUgliGeomHandle.GetVldRange() << endl;

  Double_t avgHigh = 0;
  if (fWasHigh) avgHigh = fSumHigh/fWasHigh;
  Double_t avgLow = 0;
  if (fWasLow) avgLow = fSumLow/fWasLow;

  MSG("Bfld",Msg::kInfo) << "~BfldCache " 
                         << " z outside steel " 
                         << fWasHigh << " high avg=" << avgHigh << ", " 
                         << fWasLow << " low avg=" << avgLow << endl;


}

//_____________________________________________________________________________
const BfldDbiPlaneMap* BfldCache::FindPlaneMap(TVector3& position,
                                                     Bool_t isUVZ)
{
   SetSteelAndPlaneMapCache(position,isUVZ);
   return fPlaneMapCache;
}

#ifdef OBSOLETE_METHOD
//_____________________________________________________________________________
const BfldDbiPlaneMap* BfldCache::FindPlaneMapFromId(PlexPlaneId plnid)
{
   InvalidateCaches();
   SetPlaneMapCache(plnid);
   if (plnid.IsSteel() && ! plnid.IsNull() ) {
     fSteelCache = fUgliGeomHandle.GetSteelPlnHandle(plnid);
     SetSteelLimits();
   }
   return fPlaneMapCache;
}
#endif

//_____________________________________________________________________________
void BfldCache::SetSteelLimits()
{
   // Fill fSteelZmin/fSteelZmax for given fSteelCache
   if (fSteelCache.IsValid()) {
     Double_t z0 = fSteelCache.GetZ0();
     Double_t dz = fSteelCache.GetHalfThickness() + FLT_EPSILON + fZTolerance;
     fSteelZmin = z0 - dz;
     fSteelZmax = z0 + dz;
   }
   else {
     fSteelZmin = +9999.;
     fSteelZmax = -9999.;
   }
}

//_____________________________________________________________________________
void BfldCache::SetPlaneMapCache(PlexPlaneId steelid)
{
   // Fill fPlaneMapCache entry for this steel

   // avoid DBI lookup if possible
   if ( fPlaneMapCache && fPlaneMapCache->GetPlaneId() == steelid) return;

   fPlaneMapCache = 0;

   if (fDbiMapEntries > 0) {
     const BfldDbiPlaneMap* planemap = 
       fPlaneMaps.GetRowByIndex(steelid.GetPlane());
    
     if (planemap) fPlaneMapCache = planemap;
   }
}

//_____________________________________________________________________________
void BfldCache::InvalidateCaches()
{
   // Invalidate fSteelCache, fSteelZmin/max, fPlaneMapCache
   fSteelCache    = UgliSteelPlnHandle();
   fSteelZmin     = +9999.;
   fSteelZmax     = -9999.;
   fPlaneMapCache = 0;

}

//_____________________________________________________________________________
void BfldCache::SetSteelAndPlaneMapCache(TVector3& posGlobal, Bool_t isUVZ)
{
   // Fill fSteelCache and PlaneMap entry for this position.
   // Assumes that B is only important within steel
   // Assumes that the plane is perpendicular to z direction
   // i.e. no rotation out of the x-y plane

   fFindWasInSteel = true;
   Double_t z = posGlobal.z();

   // Try the cache first
   if ( fSteelCache.IsValid() && BFieldWithinZ(z,fSteelZmin,fSteelZmax) ) {
     //MSG("Bfld",Msg::kVerbose) 
     //  << "SetSteelAndPlaneMapCacheFromZ hit cache!" << endl;
   }
   else {
     // Ask the geometry
     fSteelCache = fUgliGeomHandle.GetNearestSteelPlnHandle(z);
     SetSteelLimits();
     SetPlaneMapCache(fSteelCache.GetPlexPlaneId());
     MSG("Bfld",Msg::kDebug) << "Asked geometry for NearestSteelPln " << endl;

     if ( ! fSteelCache.IsValid() ) 
       MSG("Bfld",Msg::kError) 
         << "Asked geometry for NearestSteelPln for z= " << z 
         << ", got invalid UgliGeomSteelPln" 
         << endl;
   }

   //MSG("Bfld",Msg::kInfo)
   //    << " posGlobal " << posGlobal << " fPosInSteel " 
   //    << fPosInSteel << endl;


   // InZTest:
   //   0 = no test that z is within range of a plane
   //   1 = return (nearest) valid UgliSteelPlnHandle, but complain
   //  >1 = return (nearest) valid UgliSteelPlnHandle (quietly)
   //       in this case set fPositionInSteel = false
   if ( fInZTest == 0) return;
   // DoLocalTransform:
   //   0 = no Global <--> Local transform
   //   1 = do only for shifts (ie.  Global->Local position but not field)
   //   2 = do for shifts and rotations (ie. transform returned field)
   if ( fDoLocalTransform <= 0 ) {
     fPosInSteel = posGlobal;
     if ( ! BFieldOutsideZ(z,fSteelZmin,fSteelZmax) ) return;
   } else {
     // use the in-the-steel local coord to determine if we're in "z"
     // note Ugli uses globalInXYZ, while BField uses isUVZ
     fPosInSteel = fSteelCache.GlobalToLocal(posGlobal,!isUVZ);
     Double_t hdzlimit = fSteelCache.GetHalfThickness()+FLT_EPSILON+fZTolerance;
     if ( TMath::Abs(fPosInSteel.z()) <= hdzlimit )  return;
   }

   // currently this returns "nearest steel" when outside the steel
   // but for coil we should always return either upstream or downstream
   // steel depending on Near/Far (which is which ... think hard)
   // be careful of falling off the ends
   fFindWasInSteel = false;
   
   if ( fInZTest < 2 ) {
     
     Double_t z0 = fSteelCache.GetZ0();
     Double_t ht = fSteelCache.GetHalfThickness();
     //Double_t dz = (z>z0) ? z-(z0+ht) : (z0-ht)-z;
     Double_t dz = 0;
     const char* side;
     if ( z > z0 ) {
       side = "high";
       dz = z-(z0+ht);
       fWasHigh++;
       fSumHigh += dz;
     } else {
       side = "low";
       dz = (z0-ht)-z;
       fWasLow++;
       fSumLow += dz;
     }
     
     MAXMSG("Bfld",Msg::kWarning,10) 
       << "Geom returned plane " << fSteelCache.GetPlexPlaneId() << endl
       << " but z " << z << " is not within " << ht << " of " << z0
       << " (" << dz << " too " << side << ")" << endl;
     
   }
   //if (fInZTest>1) InvalidateCaches();

}   

//_____________________________________________________________________________
void BfldCache::FillSMZLimits()
{
   // Determine z range of SM

   fSMZLimits[0] = fSMZLimits[2] = -1.0e6;
   fSMZLimits[1] = fSMZLimits[3] = +1.0e6;

   int ipln[4] = { 0, 0, 0, 0};

   switch (fDetector) {
   case Detector::kFar: 
     {
       ipln[0] = 0;
       ipln[1] = PlexPlaneId::LastPlaneFarSM0();
       ipln[2] = PlexPlaneId::LastPlaneFarSM0()+1;
       ipln[3] = PlexPlaneId::LastPlaneFarSM1();
       break;
     }
   case Detector::kNear:
     {
       ipln[0] = 0;
       ipln[1] = PlexPlaneId::LastPlaneNearSpect();
       ipln[2] = PlexPlaneId::LastPlaneNearSpect();
       ipln[3] = PlexPlaneId::LastPlaneNearSpect();
       break;
     }
   default:
     {
       return; // ! no SM at all, use defaults
     }
   }

   for (int indx = 0; indx < 4; ++indx) {
     PlexPlaneId steelid(fDetector,ipln[indx],true);
     UgliSteelPlnHandle uph = fUgliGeomHandle.GetSteelPlnHandle(steelid);
     if ( ! uph.IsValid() ) continue; // in case of no geometry
     Double_t dz = uph.GetHalfThickness();
     if ( 0 == indx || 2 == indx ) dz = -dz;
     fSMZLimits[indx] = uph.GetZ0() + dz;
   }

}

//_____________________________________________________________________________
Ugli::SMRegion_t BfldCache::InSMRegion(const TVector3& pos, Bool_t isUVZ )
{
   // Determine if the this position is within the SuperModules
   // If in XY region of "flange" then consider it part of SM if
   // within dz of SM.
   // Return:
   //      -1 = upstream
   //       1 = in SM1
   //       0 = in gap
   //       2 = in SM2
   //      -2 = downstream

   Double_t dz = 3.4 * Munits::cm;  // extra z in the area of the collar

   Double_t z = pos.Z();

   // definitely in a SM
   if ( z >= fSMZLimits[0] && z <= fSMZLimits[1] ) return Ugli::kInSM1; // 1
   if ( z >= fSMZLimits[2] && z <= fSMZLimits[3] ) return Ugli::kInSM2; // 2

   // definitely not in a SM
   if ( z < fSMZLimits[0]-dz ) return Ugli::kUpstream;   // -1
   if ( z > fSMZLimits[1]+dz && z < fSMZLimits[2]-dz ) return Ugli::kSMGap; // 0
   if ( z > fSMZLimits[3]+dz ) return Ugli::kDownstream; // -2

   // in region where flange might be part of SM
   bool inFlange = ( InXYRegion(pos,isUVZ) <= Ugli::kFlange );
   if ( z < fSMZLimits[0] ) {
     if ( inFlange ) return Ugli::kInSM1;      //  1
     else            return Ugli::kUpstream  ; // -1;
   }
   if ( Detector::kNear == fDetector ) {
     if ( inFlange ) return Ugli::kInSM1;      //  1
     else            return Ugli::kDownstream; // -2;
   } else { // far
     if ( z > fSMZLimits[3] ) {
       if ( inFlange ) return Ugli::kInSM2;      //  2;
       else            return Ugli::kDownstream; // -2;
     }
     // FarDet gap
     if ( ! inFlange ) return Ugli::kSMGap; // 0;
     if ( z < 0.5*(fSMZLimits[1]+fSMZLimits[2]) ) return Ugli::kInSM1; // 1
     else                                         return Ugli::kInSM2; // 2
   }
   return Ugli::kSMGap; // 0
}

//_____________________________________________________________________________
Ugli::XYRegion_t BfldCache::InXYRegion(const TVector3& pos, Bool_t isUVZ )
{
  // Determine if x-y is outside flange area ...
  // NearDet is square in UV, diamond in XY
  // FarDet is round

  // this should really be part of UgliGeom

  // inside "coil" return 0... not implemented yet
  const Double_t rCoil    = 0.0;  // inside this return 0
  const Double_t rThroat  = 13.0 * Munits::cm; // inside this return 1
  const Double_t rNeck    = 14.0 * Munits::cm; // inside this return 2
  const Double_t rHole    = 15.0 * Munits::cm; // inside this return 3
  const Double_t rFlange  = 17.0 * Munits::cm; // inside this return 4
  // outside return 5

  Double_t rlimit[5] = { rCoil, rThroat, rNeck, rHole, rFlange };
  int indx;

  switch (fDetector) {
  case Detector::kNear:
    {
      Double_t u, v;
      if (isUVZ) {
        u = pos.X();
        v = pos.Y();
      } else {
        Ugli::xy2uv(fDetector,pos.X(),pos.Y(),u,v);
      }
      for ( indx = Ugli::kCoil; indx <= Ugli::kFlange; ++indx ) {
        if ( TMath::Abs(u) < rlimit[indx] &&
             TMath::Abs(v) < rlimit[indx]     ) return (Ugli::XYRegion_t)indx;
      }
      break;
    }
  case Detector::kFar:
    {
      Double_t r = pos.Perp();
      for ( indx = Ugli::kCoil; indx <= Ugli::kFlange; ++indx ) {
        if ( r < rlimit[indx] ) return (Ugli::XYRegion_t)indx;
      }
    }
  default:
    return Ugli::kSteel;
  }
  return Ugli::kSteel;
}
//_____________________________________________________________________________
Int_t BfldCache::GetDefaultMapVariant(const VldContext& vldc)
{
   // Return default map (for context)

  // no real data ...
  Int_t indx = 0;
  Int_t variant;
  switch (vldc.GetDetector()) {
  case (Detector::kNear):
    {
      switch (vldc.GetSimFlag()) {
      case (SimFlag::kMC):
      case (SimFlag::kReroot):
        variant = (indx==0) ? 160 : 0; break;
      default:
        variant = (indx==0) ? 160 : 0; break;
      }
      break;
    }
  case (Detector::kFar):
    {
      switch (vldc.GetSimFlag()) {
      case (SimFlag::kMC):
      case (SimFlag::kReroot):
        variant = (indx==0) ? 208 : 0; break;
      default:
        variant = (indx==0) ? 208 : 0; break;
      }
      break;
    }
  default:                      variant =    0; break;
  } 

  return variant;
}

//_____________________________________________________________________________
Double_t BfldCache::GetDefaultScale(const VldContext& vldc)
{
   // Lookup scale factor for this plane

  // no real data ...
  Double_t scale;
  switch (vldc.GetDetector()) {
  case (Detector::kNear):   scale =  1.0; break;
  case (Detector::kFar):    scale =  1.0; break;
  default:                      scale =    0; break;
  } 

  return scale;
}

//_____________________________________________________________________________
void BfldCache::Print(Option_t * /* option */) const
{
  MSG("Bfld",Msg::kInfo) 
    << " RefCnt=" << setw(2) << fRef
    << " " << fVldRange << endl
    << "       DoLocalTransform=" << fDoLocalTransform
    << "    InZTest=" << fInZTest
    << "    ZTolerance=" << (fZTolerance/Munits::mm) << "mm"
    << endl;
}


//_____________________________________________________________________________

---