LBNFMagneticFieldPolyconeHorn.cc

Go to the documentation of this file.

//----------------------------------------------------------------------
//LBNFMagneticFieldPolyconeHorn 
// $Id
//----------------------------------------------------------------------

#include "LBNFMagneticFieldPolyconeHorn.hh"

#include "G4VPhysicalVolume.hh"
#include "G4ios.hh"
#include "G4TransportationManager.hh"
#include "LBNEVolumePlacements.hh"
#include "LBNESurveyor.hh"
#include "Randomize.hh"
#include "LBNERunManager.hh"
#include <fstream>

//magnetic field between conductors ====================================================

LBNFMagneticFieldPolyconeHorn::LBNFMagneticFieldPolyconeHorn(size_t iH) : 
hornNumber(iH), // indexed 0 to 4..
fHornCurrent(0.),
fEffectiveLength(0.),
fOuterRadius(1.0e9),
fOuterRadiusEff(1.0e9),
fSkinDepth(4.1*CLHEP::mm),
fSkinDepthInnerRad(1.0e10*CLHEP::mm), // Unrealistic high, because that was the default sometime ago. 
fSkinDepthIIRInv(std::sqrt(2)/fSkinDepthInnerRad),
fDeltaEccentricityIO(0.),
fDeltaEllipticityI(0.),
fShift(2, 0.),
fShiftSlope(2,0.)
{
// 
// Magnetic field for the simple Polycone horn.  This field is attached to container volume 
// LBNFSimpleHornxContainer, which is a box. The field is defined the Polycone, we keep our own geometry 
// such we can implement field error in the future, should it be necessary..
//
  const LBNEVolumePlacements *aPlacementHandler = LBNEVolumePlacements::Instance();
  fZShiftDrawingCoordinate = 0.; // Check with Geantino here... 
  fZDCBegin.clear();
  fRadsInnerC.clear();
  fRadsOuterC.clear();
  double zzMinTmp = 1.0e20;
  double zzMaxTmp = -1.0e20;
  fOuterRadius = aPlacementHandler->GetHornsPolyOuterRadius(hornNumber + 1);
  fZShiftDrawingCoordinate = aPlacementHandler->GetHornsPolyZStartPos(hornNumber + 1);
  if (aPlacementHandler->IsHornPnParabolic(hornNumber + 1)) {
    size_t nPts = aPlacementHandler->GetHornParabolicNumInnerPts(iH);
    const double epsil = 0.050*CLHEP::mm; // from LBNEVolumePlacements::UpdateParamsForHornMotherPolyNum
    for (size_t k=0; k !=  nPts; k++) {
      const double zzz = aPlacementHandler->GetHornParabolicRZCoord(iH, k) + fZShiftDrawingCoordinate + epsil;
      zzMinTmp = std::min(zzMinTmp, zzz);      
      zzMaxTmp = std::max(zzMaxTmp, zzz);
      fZDCBegin.push_back(zzz);
      const double rIn = aPlacementHandler->GetHornParabolicRIn(iH, k) - epsil;
      fRadsInnerC.push_back(rIn);
      const double rOut = rIn + aPlacementHandler->GetHornParabolicThick(iH, k);
      fRadsOuterC.push_back(rOut);
    } 
  } else {
    if ((hornNumber == 0) && aPlacementHandler->GetUseLBNFOptimConceptDesignHornA()) {
     fZDCBegin = aPlacementHandler->GetHornsPolyInnerConductorZsHornARev2();
     std::cerr << " First fZDCBegin for horn A " << fZDCBegin[0]+ fZShiftDrawingCoordinate << std::endl;
     // We shift in global coordinate.. 
     std::cerr << "LBNFMagneticFieldPolyconeHorn::LBNFMagneticFieldPolyconeHorn, CD HornA, fZShiftDrawingCoordinate " 
               << fZShiftDrawingCoordinate << std::endl;
//             std::cerr << " And quit after Horn A " << std::endl; exit(2);
     for (size_t k=0; k !=fZDCBegin.size(); k++) { fZDCBegin[k] += fZShiftDrawingCoordinate; }
     fRadsInnerC = aPlacementHandler->GetHornsPolyInnerConductorRadiiHornARev2();
     fRadsOuterC = fRadsInnerC;
     for (size_t k=0; k != fRadsInnerC.size(); k++) { fRadsOuterC[k] += aPlacementHandler->GetThickICDRevisedHornA(); }
     std::cerr << " ....    /Rin/Rout map, global coordinate system " << std::endl;
     for (size_t k=0; k !=fZDCBegin.size(); k++) 
       std::cerr << " " << k << "  " 
                 << fZDCBegin[k] << " " << fRadsInnerC[k] << " " << fRadsOuterC[k] << std::endl;
     zzMinTmp = fZDCBegin[0];
     zzMaxTmp = fZDCBegin[fZDCBegin.size()-1];
    } else if ((hornNumber == 1) && aPlacementHandler->GetUseLBNFOptimConceptDesignHornB()) {
     fZDCBegin = aPlacementHandler->GetHornsPolyInnerConductorZsHornBRev2();
     // We shift in global coordinate.. 
     std::cerr << "LBNFMagneticFieldPolyconeHorn::LBNFMagneticFieldPolyconeHorn, CD HornB, fZShiftDrawingCoordinate " 
               << fZShiftDrawingCoordinate << std::endl;
     for (size_t k=0; k !=fZDCBegin.size(); k++) { fZDCBegin[k] += fZShiftDrawingCoordinate; }
     fRadsInnerC = aPlacementHandler->GetHornsPolyInnerConductorRadiiHornBRev2();
     fRadsOuterC = fRadsInnerC;
     for (size_t k=0; k != fRadsInnerC.size(); k++) { fRadsOuterC[k] += aPlacementHandler->GetThickICDRevisedHornB(); }
     std::cerr << " ....    /Rin/Rout map, global coordinate system " << std::endl;
     for (size_t k=0; k !=fZDCBegin.size(); k++) 
       std::cerr << " " << k << "  " 
                 << fZDCBegin[k] << " " << fRadsInnerC[k] << " " << fRadsOuterC[k] << std::endl;
     
     zzMinTmp = fZDCBegin[0];
     zzMaxTmp = fZDCBegin[fZDCBegin.size()-1];
    } else if ((hornNumber == 2) && aPlacementHandler->GetUseLBNFOptimConceptDesignHornC()) {
     fZDCBegin = aPlacementHandler->GetHornsPolyInnerConductorZsHornCRev2();
     // We shift in global coordinate.. 
     std::cerr << "LBNFMagneticFieldPolyconeHorn::LBNFMagneticFieldPolyconeHorn, CD HornC, fZShiftDrawingCoordinate " 
               << fZShiftDrawingCoordinate << std::endl;
     for (size_t k=0; k !=fZDCBegin.size(); k++) { fZDCBegin[k] += fZShiftDrawingCoordinate; }
     fRadsInnerC = aPlacementHandler->GetHornsPolyInnerConductorRadiiHornCRev2();
     fRadsOuterC = fRadsInnerC;
     for (size_t k=0; k != fRadsInnerC.size(); k++) { fRadsOuterC[k] += aPlacementHandler->GetThickICDRevisedHornC(); }
     std::cerr << " ....    /Rin/Rout map, global coordinate system " << std::endl;
     for (size_t k=0; k !=fZDCBegin.size(); k++) 
       std::cerr << " " << k << "  " 
                 << fZDCBegin[k] << " " << fRadsInnerC[k] << " " << fRadsOuterC[k] << std::endl;
     
     zzMinTmp = fZDCBegin[0];
     zzMaxTmp = fZDCBegin[fZDCBegin.size()-1];
    } else { // simplified geometry for the horns.
      for (size_t k=0; k!= static_cast<size_t>(aPlacementHandler->GetUseHornsPolyNumInnerPts(hornNumber + 1)); k++) {
        const double zzz = aPlacementHandler->GetHornsPolyInnerConductorZCoord(hornNumber + 1, k) + fZShiftDrawingCoordinate;
        zzMinTmp = std::min(zzMinTmp, zzz);      
        zzMaxTmp = std::max(zzMaxTmp, zzz);
        fZDCBegin.push_back(zzz);
        double rIn =  aPlacementHandler->GetHornsPolyInnerConductorRadius(hornNumber + 1, k);     
        double rOut = rIn + aPlacementHandler->GetHornsPolyInnerConductorThickness(hornNumber + 1, k);
        fRadsInnerC.push_back(rIn);
        fRadsOuterC.push_back(rOut);
     }
//       std::cerr << " k " << k << " z " << fZDCBegin[fZDCBegin.size()-1] << " InnerC " << rIn << " out " << rOut << std::endl;
   }
  }
  fEffectiveLength = zzMaxTmp - zzMinTmp;
  
//       std::cerr << " And quit for now.. I say so.. " << std::endl; exit(2);
//
  
// Use now the survey data (simulated for now.. ) to establish the coordinate transform..  
//
  LBNESurveyor *theSurvey= LBNESurveyor::Instance();
  std::ostringstream hornIndexStrStr; hornIndexStrStr << hornNumber+1;
  std::string hornIndexStr(hornIndexStrStr.str());
  std::string upstrStr("Upstream");std::string downstrStr("Downstream");
  std::string leftStr("Left"); std::string rightStr("Right"); 
  std::string ballStr("BallHorn");
  std::vector<LBNESurveyedPt>::const_iterator itUpLeft=
    theSurvey->GetPoint(G4String(upstrStr + leftStr + ballStr + hornIndexStr));
  std::vector<LBNESurveyedPt>::const_iterator itUpRight=
    theSurvey->GetPoint(G4String(upstrStr + rightStr + ballStr + hornIndexStr));
  std::vector<LBNESurveyedPt>::const_iterator itDwLeft=
    theSurvey->GetPoint(G4String(downstrStr + leftStr + ballStr + hornIndexStr));
  std::vector<LBNESurveyedPt>::const_iterator itDwRight=
    theSurvey->GetPoint(G4String(downstrStr + rightStr + ballStr + hornIndexStr));
  const G4ThreeVector pUpLeft = itUpLeft->GetPosition();
  const G4ThreeVector pUpRight = itUpRight->GetPosition();
  const G4ThreeVector pDwLeft = itDwLeft->GetPosition();
  const G4ThreeVector pDwRight = itDwRight->GetPosition();
  std::cerr << " For Horn " << hornNumber+1 << " pUpLeft " 
            << pUpLeft << " right " <<  pUpRight << " pDwLeft " 
            << pDwLeft << " effLength " << fEffectiveLength << std::endl;
  fHornIsMisaligned  = false;
  for (size_t k=0; k!=2; k++) {
    fShift[k] = -0.50*(pUpLeft[k] + pUpRight[k]); // minus sign comes from the global to local. 
                                                  // Transverse shoft at the upstream  entrance of Horn1 
    fShiftSlope[k] = 0.5*(pUpLeft[k] + pUpRight[k] -  pDwLeft[k] - pDwRight[k])/fEffectiveLength; 
    if (std::abs(fShiftSlope[k]) > 1.0e-8) fHornIsMisaligned = true;
  }
  if (fHornIsMisaligned) {
    std::ostringstream aNameStrStr; aNameStrStr << "LBNFSimpleHorn" << (hornNumber+1) << "Container";
    std::string namePlacedVol(aNameStrStr.str()); 
    if (hornNumber == 0) namePlacedVol = std::string("Horn1PolyM1");
    const LBNEVolumePlacementData* plDat= 
       aPlacementHandler->Find("BFieldInitialization", namePlacedVol.c_str(), "BFieldInitialization");
    fRotMatrixHornContainer = plDat->fRotation;
    std::cerr << " Horn " << hornNumber+1 << " is misaligned " 
              <<  " xz term " << fRotMatrixHornContainer[0][2] << " yz " << fRotMatrixHornContainer[1][2] << std::endl;
//    fHornIsMisaligned = false;
//    std::cerr << " !!!! Back door setting of the HornisMialigned, we will not rotate the field " << std::endl; 
  }
  
//
// Open a file to study the value of the field for displaced trajectory... 
//  
// if (!fOutTraj.is_open()) {
//   G4String fNameStr("./FieldTrajectories_Horn");
//   if (amHorn1) fNameStr += G4String("1_1.txt");
//   else  fNameStr += G4String("2_1.txt");
//   fOutTraj.open(fNameStr.c_str()); // we will place a GUI interface at a later stage... 
//   fOutTraj << " id x y z bx by " << std::endl;
// }
//   if (hornNumber == 1) this->TestDivergence(0); Not to be done here, the correction term are not yet declared. 
// See LBNEDetectorConstruction,  at debug line "Defining Polycone Mag field, usedNumberOfHornsPoly"
// 
// Build the field if the Quadrupole amplitude greater than 1. 
//
    fField3DMapCurrentEqualizer.clear(); // Optional fill below.. Once we know the parameters for the simulation.. 
    fDebugIsOn = false;
    fCurrentEqualizerQuadAmpl = 0.;
    fCurrentEqualizerOctAmpl = 0.;
    fFileNameFieldMapForCE = std::string("?");
}

void LBNFMagneticFieldPolyconeHorn::GetFieldValue(const double Point[3],double *Bfield) const
{
//   std::cerr << " LBNFMagneticFieldPolyconeHorn::GetFieldValue, z = " << Point[2] << std::endl;
   for (size_t k=0; k!=3; ++k) Bfield[k] = 0.;
   const double zGlobal = Point[2];
   // Intialization of the coordinate transform. 
//   const LBNEVolumePlacements *aPlacementHandler = LBNEVolumePlacements::Instance(); // no longer needed here. See above.. 
   std::vector<double> ptTrans(3,0.);
   for (size_t k=0; k != 2; k++) 
     ptTrans[k] = Point[k] + fShift[k] + (zGlobal-fZDCBegin[0])*fShiftSlope[k]; 
   if ((fCurrentEqualizerQuadAmpl > 1.) && (fField3DMapCurrentEqualizer.size() == 0)) {
     fDebugIsOn = false; // O.K. for Now, until debugged. 
     this->fill3DMapCurrentEqualizer();
     fDebugIsOn = false; 
   }
   if (fCurrentEqualizerQuadAmpl > 1.) {
   // From field map.. 
      double bFieldFromMap[3];
      ptTrans[2] = zGlobal;
      this->getFieldFrom3DMapCurrentEqualizer(&ptTrans[0], bFieldFromMap);
      Bfield[0] = bFieldFromMap[0];
      Bfield[1] = bFieldFromMap[1];
      if (fHornIsMisaligned) { 
        std::vector<double> bFieldLocal(3); 
        for(size_t k=0; k != 3; k++) { bFieldLocal[k] = Bfield[k]; Bfield[k] = 0.; } 
        for (size_t k1=0; k1 !=3; k1++) {
          for (size_t k2=0; k2 !=3; k2++) Bfield[k1] += fRotMatrixHornContainer[k1][k2]*bFieldLocal[k2];
        }
      }
      return; // No Correction for skin depth effect.  Will be included in the map, if need be...  
   }                                    
   G4double current = fHornCurrent/CLHEP::ampere/1000.;
//   std::cerr << " LBNFMagneticFieldPolyconeHorn::GetFieldValue Z gobal  " << Point[2] 
//            << " fZDCBegin, first " << fZDCBegin[0] << " last " <<  fZDCBegin[fZDCBegin.size() -1]
//             << " for horn " << (hornNumber+1) << " current " << current << std::endl;
   const double r = std::sqrt(ptTrans[0]*ptTrans[0] + ptTrans[1]*ptTrans[1]); 
   const double phi = std::atan2(Point[1], Point[0]);
   if (r > fOuterRadius) return;
//   if ( r < 1.0e-3) return; // The field should go smoothly to zero at the center of the horm No longer valid if Eccentricity 

   double magBField = current / (5.*r/CLHEP::cm)/10*CLHEP::tesla; //B(kG)=i(kA)/[5*r(cm)], 1T=10kG
   // Implement the azimuthal anisotropy due to imperfection of the current equalizer section.
   const bool doCurrentEqualizerCorr = ((std::abs(fCurrentEqualizerQuadAmpl) > 1.0e-9) || 
                                       (std::abs(fCurrentEqualizerOctAmpl) > 1.0e-9) );
   if (doCurrentEqualizerCorr) { 
      double zFactAzi = 1.0e-6;
      switch(hornNumber) {
        case 0:
        case 2:
          zFactAzi = ((zGlobal - fZDCBegin[0]) > 0.) ? (zGlobal - fZDCBegin[0])/fCurrentEqualizerLongAbsLength : 1.0e-6;
          break;
        case 1:
          zFactAzi = ((fZDCBegin[fZDCBegin.size() -1] - zGlobal) > 0.) ? 
                         (fZDCBegin[fZDCBegin.size() -1] - zGlobal)/fCurrentEqualizerLongAbsLength : 1.0e-6; 
          break;
        default:
         zFactAzi = 1.0e-6; // undefined, only 3 horns   
      }
      const double factCE = std::exp(-1.0*std::abs(zFactAzi)) * (fCurrentEqualizerQuadAmpl*std::sin(2.0*phi) 
                                                   + fCurrentEqualizerOctAmpl*std::sin(4.0*phi));
      magBField *= (1.0 + factCE);                                         
//      std::cerr << " ... Current Equalizer factor " << factCE << " Z " << zFactAzi << std::endl;                                          
   }
   

   double radIC=0.;
   double radOC=0.;   
   this->fillHornPolygonRadii(zGlobal, &radIC, &radOC);
   if (std::abs(radIC) > 1.0e6) {
      fEffectiveInnerRad = -1.0;
      fEffectiveOuterRad = -1.0;
//      std::cerr << " No effective radius at z " << zGlobal << std::endl;
      return;
   } else {
//      std::cerr << " At zGlobal " << zGlobal << " zLocal " << ptTrans[2] << " radIC " 
//                << radIC << " radOC " << radOC <<  std::endl; 
   } 
   fEffectiveInnerRad = radIC;
   fEffectiveOuterRad = radOC;
   fEffectiveSkinDepthFact = 1.0;
     // Study of systematic effects. Define an effective radIC wich is large than the physics radIC.
     // We assume here that the finite skin depth matter, and the current density is not 
     // uniform. Simply increase the radIC
   if (r < radIC && (fDeltaEccentricityIO < 1.0e-20)) {
       return; // zero field region (Amps law). 
    }
    // We assume that the eccentricity is along the X axis.  The resulting dipole will then be vertical.  It does not depend on r 
    // ( the variable magField goes like 1/r ) 
    // Ref: NuMI note # 710, December 14 2000. 
    //
    double dipoleAtInnerInner = 0.;
    if (std::abs(fDeltaEccentricityIO) > 1.0e-20) { 
       // July 22 2017: flip this sign.  
      dipoleAtInnerInner  = (-1.0)*magBField*r*fDeltaEccentricityIO/(radOC*radOC - radIC*radIC);
      if (r < radIC && (std::abs(fDeltaEccentricityIO) > 1.0e-20)) {
       Bfield[0] = 0.; 
       Bfield[1] = dipoleAtInnerInner; // Is this realy indepedent of phi? Or an approximation
       // We skip the eventual additional complexity of the mis-aligned horn.. 
       return; 
     }
    } else {
      if (r < radIC) return;
    }
//     const double deltaRIC = radOC - radIC;
// Version from Dec 18-19 2013.      
//     radIC += deltaRIC*std::min(0.99, fSkinDepthCorrInnerRad); 
     //
//     if ((r > radIC) && (r < radOC)) {
//      const double surfCyl = radOC*radOC - radIC*radIC;
//      const double deltaRSq = (r*r - radIC*radIC);
//      magBField *= deltaRSq/surfCyl; // Assume uniform current density and apply Amps law. 
//     }
// Skin Depth effect, current disctribution in the conductor itself, 
//  and Version from Zarko's thesis. ( (MINOS Document 5694-v1) ) 
// 
     double magFieldPhi = magBField;
     double magFieldR = 0.;
     
     double factRInAlum = 1.0;  
     if ((r > radIC) && (r < radOC)) {
        const double surfCyl = radOC*radOC - radIC*radIC;
        const double deltaRSq = (r*r - radIC*radIC);
        factRInAlum = deltaRSq/surfCyl; // Inside the alumimum  
        fEffectiveSkinDepthFact = factRInAlum;
       if (fSkinDepthInnerRad > 10.*CLHEP::mm) { // default assume very large skin depth, assume constant current 
                                      // current density .
         magFieldPhi *= factRInAlum; // Assume uniform current density and apply Amps law. 
       } else { // Realistic skin depth. 
         magFieldPhi *= factRInAlum*getNormCurrDensityInnerC((r-radIC), (radOC-radIC))
                            /getNormCurrDensityInnerC((radOC-radIC), (radOC-radIC));
       }
     } 
// Sin and cos factor bphi -> bx, by 
//     std::cerr << " ....check r " << r << " magFieldPhi " 
//               << magFieldPhi << " magBField " << magBField << std::endl;
// 
// See informal document given to Alberto M and Laura F. regarding the solution of 
// divergence equation, for the multipole expansion. Sept 29, 2016. 
//
     double r12Fact = 1.0;
     if (std::abs(fDeltaEllipticityI) > 1.0e-20) { 
       // NuMi note 0710 and/or NuMI 0471 
       const double rDelta = sqrt((ptTrans[0] - fDeltaEccentricityIO)*(ptTrans[0] - fDeltaEccentricityIO) +
                              ptTrans[1]*ptTrans[1]);
       r12Fact = fEffectiveInnerRad*fEffectiveInnerRad/
           (fEffectiveOuterRad*fEffectiveOuterRad - fEffectiveInnerRad*fEffectiveInnerRad);
       double deltaMagFieldPhiEc = magBField * r12Fact * (fDeltaEccentricityIO / rDelta) * std::cos(phi);
       double magFieldREc = magBField * r12Fact * (fDeltaEccentricityIO / rDelta) * std::sin(phi);
       double deltaMagFieldPhiEl = magBField * r12Fact * fEffectiveInnerRad * (fDeltaEllipticityI /(r*r)) * std::cos(2.0*phi);
       double magFieldREl = magBField * r12Fact * fEffectiveInnerRad * (fDeltaEllipticityI /(r * r)) * std::sin(2.0*phi);
       // Reintroduce the bug, by setting 
//       factRInAlum = 1.;
       magFieldPhi = factRInAlum*(magBField - deltaMagFieldPhiEl - deltaMagFieldPhiEc);
       magFieldR = factRInAlum*(magFieldREl + magFieldREc);       
//       magFieldPhi = (magBField - deltaMagFieldPhiEl - deltaMagFieldPhiEc);
//       magFieldR = (magFieldREl + magFieldREc); 
 // To be checked for ellipticity.       
     }
     Bfield[0] = -magFieldPhi*ptTrans[1]/r +  magFieldR*ptTrans[0]/r; 
     Bfield[1] = magFieldPhi*ptTrans[0]/r +  magFieldR*ptTrans[1]/r;
     //
     // Analytical calculation, in cartesian coordinate, for the eccentricity. 
     //
     if (std::abs(fDeltaEccentricityIO) > 1.0e-20)  {
        // Two Cartesian coordinate system, the first one, centered on Outer surface of the IC 
        // The second one, translated by delta, say to the left, and with critical radius radIn 
        const double ITot = magBField * r; // Up to a constant.. 
        const double JTot = ITot/(radOC*radOC - radIC*radIC);
        const double BField1Phi = (r < radOC) ? ( JTot * r) : (JTot * radOC*radOC/ r);
        const double BField1X = -1.0*BField1Phi*std::sin(phi); 
        const double BField1Y = BField1Phi*std::cos(phi);
        const double h1 = r * std::sin(phi);
        const double l2 = fDeltaEccentricityIO + r*std::cos(phi);
        const double r2 = std::sqrt(l2*l2 + h1*h1); // with respect to the displaced inner radius. 
        const double phi2 = std::atan2(h1, l2);
        const double BField2Phi = (r2 < radIC) ? ( JTot * r2) : (JTot * radIC*radIC/ r2);
        const double BField2X = -1.0*BField2Phi*std::sin(phi2); 
        const double BField2Y = BField2Phi*std::cos(phi2);
        Bfield[0] = BField1X - BField2X;
        Bfield[1] = BField1Y - BField2Y;
     }
//     std::cerr << " Bfield .. " << Bfield[0]/CLHEP::tesla << " / " << Bfield[1]/CLHEP::tesla  
//               << " r " << r << " radIC " << radIC << " radIC-eff " << fEffectiveInnerRad << std::endl;
     if (doCurrentEqualizerCorr) { 
        const double kOfPhiEQ = 
            (+1.0) * ( 2.0*fCurrentEqualizerQuadAmpl*std::cos(2.0*phi) + 
                                 4.0*fCurrentEqualizerOctAmpl*std::cos(4.0*phi)) * std::log(r/radIC); 
                                 // See 2-page solution
                                 // The sign has been flipped, because the sign of magBField has been flipped 
                                 //for the main component.. 
                                 // Based on the arrow plot, and the divergence test. 
         Bfield[0] += (magFieldPhi*ptTrans[0]/r) * kOfPhiEQ;
         Bfield[1] += (magFieldPhi*ptTrans[1]/r) * kOfPhiEQ; 
         // We don't bother iwth the dadial componenent, too small.. 
     }
//
// Important correction.. May 27 2017... ( a bug from a long time ago !?)
// One needs to rotate the field as well.
    if (fHornIsMisaligned) { 
       std::vector<double> bFieldLocal(3); 
       for(size_t k=0; k != 3; k++) { bFieldLocal[k] = Bfield[k]; Bfield[k] = 0.; } 
       for (size_t k1=0; k1 !=3; k1++) {
         for (size_t k2=0; k2 !=3; k2++) Bfield[k1] += fRotMatrixHornContainer[k1][k2]*bFieldLocal[k2];
       }
//       std::cerr << " Misaligned Horn" << hornNumber+1 << ", rotation of the field ... " << std::endl;
//       for(size_t k=0; k != 3; k++) std::cerr << " Component " << k 
//                                              << " local " << bFieldLocal[k] << " global " << Bfield[k] << std::endl;
    }
}
void LBNFMagneticFieldPolyconeHorn::dumpField() const {
//  this->TestEccentricityBerkeley();
  fDebugIsOn = true;
  std::ostringstream fNameStr;  fNameStr << "./FieldMapHornTmpV3" << (hornNumber +1);
  if (std::abs(fDeltaEccentricityIO) > 1.0e-20) fNameStr << "Ecc.txt";
  else fNameStr << ".txt";
  std::string fName(fNameStr.str());
  std::ofstream fOut(fName.c_str());
//  const double zStart = fZShiftDrawingCoordinate;
//  const double zEnd = zStart + fEffectiveLength + 1000. ;
//  Study Horn A only.
  const double zStart = 2219.;  
//  const double zEnd = 2200.0000001;  
  const double zEnd = 2219.0000001;  
  std::cerr << "LBNFMagneticFieldPolyconeHorn::dumpField .. fZShiftDrawingCoordinate " 
            << fZShiftDrawingCoordinate << " length " << fEffectiveLength 
            << " Z start/end " << zStart << "/" << zEnd << std::endl;
//  if ((std::abs(fDeltaEllipticityI) < 1.0e-20) && (std::abs(fDeltaEccentricityIO) < 1.0e-20) ) {
// Splitting this method in two parts was a bad idea, as it prevented us to make maps with the exactly the same grid,
// but different option.. I'll leave the code there for now, in case we want to generate simple maps, but disable it..  
   if (std::abs(zEnd) < -1.e10) { // to ensure the same map. 
    fOut << " z zd r bphi br  " << std::endl;
    double rMax = fOuterRadius + 5.0*CLHEP::mm;
//  if (amHorn1) std::cerr << " Horn1, Dumpfield, r Max = " << rMax << std::endl;
//  else std::cerr << " Horn2, Dumpfield, r Max = " << rMax << std::endl;
    const int numZStep = 200;
    const int numRStep= 200;
    const double zStep = (zEnd-zStart)/numZStep;
    const double rStep = rMax/numRStep;
    double point[3];
    for (int iZ=0; iZ !=numZStep; iZ++) { 
      double z = zStart +  iZ*zStep;
      point[2] = z;
      const double zLocD = point[2];
      double radIC = 0.; double radOC = 0.;
      this->fillHornPolygonRadii(zLocD, &radIC, &radOC);
      double r = radIC - 0.25*CLHEP::mm;
    // inside the conductor... 
      while (r < rMax/2.)  {
        const double phi = 2.0*M_PI*G4RandFlat::shoot();
        point[0] = r*std::cos(phi);
        point[1] = r*std::sin(phi);
        double bf[3];
        this->GetFieldValue(point, &bf[0]);
        const double bphi = bf[0]*std::sin(phi) - bf[1]*std::cos(phi);
        const double br = bf[0]*std::cos(phi) + bf[1]*std::sin(phi);
        fOut << " " << z << "  " << zLocD << " " << r << " " 
           << bphi/CLHEP::tesla << " " << br/CLHEP::tesla << std::endl;
        if (r < (radIC + 5.0*CLHEP::mm)) r += 0.1*CLHEP::mm;
        else  r += rStep;
     }
    }
  } else { 
    fOut << " z zd x y r bphi br bx by " << std::endl;
//  if (amHorn1) std::cerr << " Horn1, Dumpfield, r Max = " << rMax << std::endl;
//  else std::cerr << " Horn2, Dumpfield, r Max = " << rMax << std::endl;
//    const int numZStep = 200;
    const int numZStep = 1;
    const int numRStep= 200;
    const int numPhiStep= 100;
    const double zStep = (zEnd-zStart)/numZStep;
    double point[3];
    const double rInit = 20.;
//    const double rMax = radIC + 10.0*CLHEP::mm;
    const double rMax = 50.0*CLHEP::mm;
//    const double rMax = radIC - 0.1*CLHEP::mm;
      const double rStep = (rMax - rInit)/numRStep;
    for (int iZ=0; iZ !=numZStep; iZ++) { 
      double z = zStart +  iZ*zStep;
      point[2] = z;
      const double zLocD = point[2];
      double radIC = 0.; double radOC = 0.;
      this->fillHornPolygonRadii(zLocD, &radIC, &radOC);
    // inside the conductor... 
      double r = rInit;
      while (r < rMax)  {
        for (size_t kPhi=0; kPhi!= (numPhiStep); kPhi++) { 
          const double phi = (2.0*M_PI*kPhi)/numPhiStep;
          point[0] = r*std::cos(phi);
          point[1] = r*std::sin(phi);
          double bf[3];
          this->GetFieldValue(point, &bf[0]);
          const double bphi = bf[0]*std::sin(phi) - bf[1]*std::cos(phi);
          const double br = bf[0]*std::cos(phi) +  bf[1]*std::sin(phi);
          fOut << " " << z << "  " << zLocD << " " << point[0] << " " << point[1] << " " << r << " " 
           << bphi/CLHEP::tesla << " " << br/CLHEP::tesla 
           << " " << bf[0]/CLHEP::tesla << " " << bf[1]/CLHEP::tesla << std::endl;
         }
         r += rStep;
      }
    }
  }  
   std::cerr << " And quit for now... " << std::endl; exit(2);
  fOut.close();
  //
  // Provide the arrow plot to visualize the correction to the Current equalizer section.. 
  //
  std::ostringstream fNameAStr;  
  fNameAStr << "./FieldMapHorn" << (hornNumber +1) 
            << "_ArrowPlot_" << fNumZPtsForCE << "_" << fNumRPtsForCE << "_" << fNumPhiPtsForCE << "_V2.txt";
  std::string fNameA(fNameAStr.str());
  fOut.open(fNameA.c_str());
  fOut << " x y r bx by " << std::endl;
  const int numR = 20;
  const int numPhi = 10.; 
  const double stepPhi = (M_PI/2.) / numPhi;
  G4ThreeVector dir;
  double pos[3]; 
  size_t iPtZSel = fRadsInnerC.size()/2;
  pos[2] = fZDCBegin[iPtZSel]; 
  const double zLoc = pos[2];
  double radICAr=0.; double radOCAr=0.;
  this->fillHornPolygonRadii(zLoc, &radICAr, &radOCAr);
  const double rStep = (fOuterRadius - radOCAr)/numR;

  double bField[3];  
  for (int iR = 0; iR != numR; iR++) { 
    double r = radOCAr + iR*rStep;
    if (iR == 0) r = radOCAr + 5.0*CLHEP::mm;
    else r = radOCAr + 2.0*CLHEP::mm + iR*rStep; // for graphical clarity ..
    for (int iPhi = 0; iPhi != numPhi; iPhi++) { 
      const double phi = 0.3333*stepPhi + iPhi*stepPhi;
      pos[0] = r * std::cos(phi); pos[1] = r * std::sin(phi);
      this->GetFieldValue(pos, bField);
      fOut << " " << pos[0] << " " << pos[1] << " " << r << " " 
           << bField[0]/CLHEP::tesla << " " << bField[1]/CLHEP::tesla << std::endl;
    }
  }
  fOut.close();
  //
  // Special map to debug the CurrentEqualizer induced azimuthal variation in the current. 
  
  if (fCurrentEqualizerQuadAmpl > 1.) {
    std::ostringstream fNameBStr;  fNameBStr << "./FieldMapHornCEQ_" << (hornNumber +1) << "_Debug.txt";
    std::string fNameB(fNameBStr.str());
    fOut.open(fNameB.c_str());
    fOut << " x y z bx by bNorm br bphi" << std::endl;
    const int numRCE = 50;
    const int numPhiCE = 100.; 
    const double stepPhiCE = (2.0*M_PI) / numPhiCE;
    double posCE[3]; 
    const double zGrid = fField3DMapCurrentEqualizerZMin + (2*fNumZPtsForCE/3)*fField3DMapCurrentEqualizerZStep;
    double radIC; double radOC; 
    this->fillHornPolygonRadii(zGrid, &radIC, &radOC);
    double bFieldCE[3];  
    for (int iZCE=0; iZCE != 10; iZCE++) { 
       posCE[2] = zGrid + (iZCE-1)*fField3DMapCurrentEqualizerZStep/5.;
       const double rStepCE = (fOuterRadius - radOC- 2.5*CLHEP::mm)/numRCE;
       for (int iR = 0; iR != numRCE; iR++) { 
          const double r = radOC -1.*CLHEP::mm + iR*rStepCE;
          for (int iPhi = 0; iPhi != numPhiCE; iPhi++) { 
            const double phi = iPhi*stepPhiCE;
             posCE[0] = r * std::cos(phi); posCE[1] = r * std::sin(phi);
//           if ((iR == 15) && (iPhi > 20) && (iPhi < 25)) fDebugIsOn = true;
//           if ((r < 225.) && (iPhi > 20) && (iPhi < 25)) fDebugIsOn = true;
             if ((std::abs(r - 205.997) < 1.) && (std::abs(phi) < 1.0)) fDebugIsOn = true;
             else fDebugIsOn = false;
             this->GetFieldValue(posCE, bFieldCE);
             const double bNorm  = sqrt(bFieldCE[0]*bFieldCE[0] + bFieldCE[1]*bFieldCE[1]);
             const double bR = bFieldCE[0]*std::cos(phi) + bFieldCE[1]*std::sin(phi);
             const double bPhi = bFieldCE[1]*std::cos(phi) - bFieldCE[0]*std::sin(phi);
             
             fOut << " " << posCE[0] << " " << posCE[1] << " " <<  posCE[2] << " " 
                << bFieldCE[0]/CLHEP::tesla << " " << bFieldCE[1]/CLHEP::tesla << " " << bNorm/CLHEP::tesla 
                << " " << bR/CLHEP::tesla << " " << bPhi/CLHEP::tesla << std::endl;
      }
    }
  }
 }
 fOut.close(); 
}  

void LBNFMagneticFieldPolyconeHorn::fillTrajectories(const double Point[3], double bx, double by) const {
 if (!fOutTraj.is_open()) return; 
 LBNERunManager* pRunManager = dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());
 fOutTraj << " " << pRunManager->GetCurrentEvent()->GetEventID();
 for (size_t k=0; k!= 3; k++) fOutTraj << " " << Point[k];
 fOutTraj << " " << bx/CLHEP::tesla << " " << " " << by/CLHEP::tesla << std::endl;
 // we flush, the destructor never called.. 
 fOutTraj.flush();
}

void LBNFMagneticFieldPolyconeHorn::TestDivergence(int opt) const {
//
// Test that it satisfy Maxwell equation, div B  = 0
//
// Since we compute in polar coordinate system, do it Cartesian.. 
//
   G4ThreeVector dir;
   double pos[3]; 
   size_t iPtZSel = fRadsInnerC.size()/2;
   pos[2] = fZDCBegin[iPtZSel]; 
   double radIC=0.;
   double radOC=0.;   
   radOC += 1.5*CLHEP::mm;
   this->fillHornPolygonRadii(pos[2], &radIC, &radOC);
   std::string fNameOut;
   switch (opt) {
     case 0:
        dir = G4ThreeVector(std::sqrt(1.0-0.01), 0.1, 0.);
        pos[0] = 1.1*radOC; pos[1] = 0.1; 
        fNameOut=std::string("./LBNFMagneticFieldPolyconeHorn_TestDivergence_X_v1.txt");
        break;
     case 1: 
        dir = G4ThreeVector(0.1, std::sqrt(1.0-0.01), 0.);
        pos[1] = 1.1*radOC; pos[0] = 0.1; 
        fNameOut=std::string("./LBNFMagneticFieldPolyconeHorn_TestDivergence_Y_v1.txt");
        break;
      case 2:
        dir = G4ThreeVector(0.480729, sqrt(1.0 - (0.480729*0.480729)), 0.);
        pos[0] = 1.1*radIC/std::sqrt(2.); pos[1] = pos[0]; 
        fNameOut=std::string("./LBNFMagneticFieldPolyconeHorn_TestDivergence_U_v1.txt");
        break;
      default: 
         std::cerr << "LBNFMagneticFieldPolyconeHorn::TestDivergence Valid switches so far are 0, 1 or 2" << std::endl;
         exit(1);
   }
   std::ofstream fOut(fNameOut.c_str());
   fOut << " x y r bx by b dbdx dbdy obnablab " << std::endl;
   const double step = 0.002*CLHEP::mm;
   double deltaMaxCart = 0.;
   const double range = fOuterRadius - fRadsOuterC[iPtZSel]- 2.5*CLHEP::mm; 
   const int numPts = static_cast<int> (range/step);
   std::cerr << " LBNFMagneticFieldPolyconeHorn::TestDivergence, number of pts " 
             << numPts << " Z location for scan " << pos[2] << " radii " 
             << radIC << " / " << fOuterRadius << std::endl;
   double bFieldPrev[3]; double bField[3]; double posPrevious[3];
   for (size_t kk=0; kk != 3; kk++) { posPrevious[kk] = pos[kk];  bFieldPrev[kk] = 0.; }
   double dbVal[3]; dbVal[2] = 0.; // assume Bz = 0.
   double bdlIntegral = 0.;
   for (int iStep = 0; iStep != numPts; iStep++) {
     this->GetFieldValue(pos, bField);
     const double bNorm = std::sqrt( bField[0]*bField[0] + bField[1]*bField[1]);
     
     if (bNorm < 1.0e-10) {
       for (size_t k=0; k!=3; k++) {
         bFieldPrev[k] = 0.;
         posPrevious[k] = pos[k];
         pos[k] += step*dir[k];
       }
        continue;
     }
     bdlIntegral += bField[0]*step*dir[1] - bField[1]*step*dir[0];
     const double phi = std::atan2(pos[1], pos[0]);
     if (iStep != 0) {
       dbVal[0] = std::cos(phi) * (bField[0] - bFieldPrev[0])/(pos[0] - posPrevious[0]);
       dbVal[1] = std::sin(phi) * (bField[1] - bFieldPrev[1])/(pos[1] - posPrevious[1]);
       const double dd = (1.0/bNorm)*(dbVal[0] + dbVal[1]);
       fOut << " " << pos[0] << " " << pos[1] << " " << std::sqrt(pos[0]*pos[0] + pos[1]*pos[1]) 
            << " " << bField[0]/CLHEP::tesla << " " << bField[1]/CLHEP::tesla 
            << " " << bNorm/CLHEP::tesla << " " 
            << dbVal[0]/bNorm << " " << dbVal[1]/bNorm << " " << dd << std::endl;
       deltaMaxCart = std::max(std::abs(dd), deltaMaxCart);
     }
     for (size_t k=0; k!=3; k++) {
       bFieldPrev[k] = bField[k];
       posPrevious[k] = pos[k];
       pos[k] += step*dir[k];
     }
  }
  std::cerr << " LBNFMagneticFieldPolyconeHorn::TestDivergence, final max deviation (rel) " 
            << deltaMaxCart  << " Integral bdl " << bdlIntegral/(CLHEP::tesla*CLHEP::meter) << std::endl;
  fOut.close();     
//  std::cerr << " And quit for now " << std::endl; exit(2);

}
void LBNFMagneticFieldPolyconeHorn::getFieldFrom3DMapCurrentEqualizer(
                              const double Point[3], double *Bfield) const {
  
  if (fDebugIsOn) std::cerr << " LBNFMagneticFieldPolyconeHorn::getFieldFrom3DMapCurrentEqualizer xyz  " 
                            << Point[0] << " / " << Point[1] << " / " <<  Point[2] << std::endl;
  for (int k=0; k !=3; k++) Bfield[k] = 0.;
  if (fField3DMapCurrentEqualizer.size() == 0) return; // should not happen
  if (Point[2] < fField3DMapCurrentEqualizerZMin) return;
  const unsigned int iZ1 = 
    static_cast<unsigned int>(std::floor((Point[2] - fField3DMapCurrentEqualizerZMin)/fField3DMapCurrentEqualizerZStep));
  if( iZ1 >= fNumZPtsForCE) return;
  const unsigned int iZ2 = (iZ1 == (fNumZPtsForCE-1)) ? iZ1 : (iZ1 + 1); 
  if(fDebugIsOn) std::cerr << " ..... Z " << Point[2] << " indices " << iZ1 << " / " <<  iZ2 << std::endl;
  double phi = std::atan2(Point[1], Point[0]); 
  if (phi < 0.) phi += 2.0*M_PI;
  const unsigned int iQuadrant = (unsigned int) std::floor(2.0*phi/M_PI);
  const double phiFirst = phi - iQuadrant*M_PI/2.;
  const unsigned int iPhi1 = 
    static_cast<unsigned int>(phiFirst/fField3DMapCurrentEqualizerPhiStep);
  const unsigned int iPhi2 = (iPhi1 == (fNumPhiPtsForCE-1)) ? 0 : (iPhi1 + 1); 
    
  if (fDebugIsOn) std::cerr << " ..... phi " << phi << " phiFirst " << phiFirst 
                            << " indices " << iPhi1 << " / " <<  iPhi2 << std::endl;
  const double rSqRequested = Point[0]*Point[0] + Point[1]*Point[1];
  const double rSqZ1 =  rSqRequested - fField3DMapCurrentEqualizerRadSqAtZ[iZ1]; 
  const double rSqZ2 =  rSqRequested - fField3DMapCurrentEqualizerRadSqAtZ[iZ2]; 
  if ((rSqZ1 < 0.) && (rSqZ2 < 0.)) return;
  double radOCAtZ = 0.; 
  double radICAtZ = 0.; // irrelevant if no OC crossing betwen Z1 and Z2
  if ((rSqZ1 < 0.) || (rSqZ2 < 0.))  {
       this->fillHornPolygonRadii(Point[2], &radICAtZ,  &radOCAtZ);
       if (rSqRequested < (radOCAtZ*radOCAtZ)) return;
  }  
  const unsigned int iR1Z1 = (rSqZ1 < 0.) ? 0 : std::floor(rSqZ1 /fField3DMapCurrentEqualizerRadSqRStepAtZ[iZ1]);
  const unsigned int iR2Z1 = (rSqZ1 < 0.) ? 0 : ((iR1Z1 == (fNumRPtsForCE-1)) ? iR1Z1 : (iR1Z1 + 1));   
  const unsigned int iR1Z2 = (rSqZ2 < 0.) ? 0 : std::floor(rSqZ2 /fField3DMapCurrentEqualizerRadSqRStepAtZ[iZ2]);
  const unsigned int iR2Z2 = (rSqZ2 < 0.) ? 0 : ((iR1Z2 == (fNumRPtsForCE-1)) ? iR1Z2 : (iR1Z2 + 1)); 
  if ((iR1Z1 >=  fNumRPtsForCE) || (iR1Z2 >=  fNumRPtsForCE)) return;
  // bilinear interposlation on the phiXRsq plane, the linear interpolation on Z 
  if (fDebugIsOn) std::cerr << " ..... RsQ, Z1 " << rSqZ1 << " at Z2 " << rSqZ2 << " indices " 
                            << iR1Z1 << " / " <<  iR2Z1 << " / " << iR1Z2 << " / " <<  iR2Z2 << std::endl;
  
  const unsigned int ii1 = fNumZPtsForCE*(fNumRPtsForCE*iPhi1 + iR1Z1) + iZ1;
  const unsigned int ii2 = fNumZPtsForCE*(fNumRPtsForCE*iPhi1 + iR2Z1) + iZ1;
  const unsigned int ii3 = fNumZPtsForCE*(fNumRPtsForCE*iPhi2 + iR2Z1) + iZ1;
  const unsigned int ii4 = fNumZPtsForCE*(fNumRPtsForCE*iPhi2 + iR1Z1) + iZ1;
  const double ui = (iPhi2 == 0) ?  (phiFirst - (fNumPhiPtsForCE - 1)* 
                                   fField3DMapCurrentEqualizerPhiStep)/ fField3DMapCurrentEqualizerPhiStep : 
                     (phiFirst - iPhi1*fField3DMapCurrentEqualizerPhiStep)/fField3DMapCurrentEqualizerPhiStep;
                    
  if (fDebugIsOn) std::cerr << " ..... ii indices " << ii1 << " / " << ii2 << " / " <<  ii3 << " / " << ii4 
            << " ui " << ui << std::endl;
  double bxIz1 = 0.; double byIz1 = 0.;
  // Same at the point iZ2.. 
  //
  const unsigned int jj1 = fNumZPtsForCE*(fNumRPtsForCE*iPhi1 + iR1Z2) + iZ2;
  const unsigned int jj2 = fNumZPtsForCE*(fNumRPtsForCE*iPhi1 + iR2Z2) + iZ2;
  const unsigned int jj3 = fNumZPtsForCE*(fNumRPtsForCE*iPhi2 + iR2Z2) + iZ2;
  const unsigned int jj4 = fNumZPtsForCE*(fNumRPtsForCE*iPhi2 + iR1Z2) + iZ2;
  double bxIz2 = 0.; double byIz2 = 0.;
  if (((iR2Z1 == iR1Z1) || (iR2Z2 == iR1Z2)) && (iR1Z1 != 0) &&
       (iR1Z2 != 0) && (iR2Z1 != 0) && (iR2Z2 != 0) ) { // only do the interpolation in phi . Far away from the IC, 
                                             // does not matter which interpolation table section we use.. Field too small.
    std::map<unsigned int, std::pair<float, float> >::const_iterator itm1 = 
         fField3DMapCurrentEqualizer.find(ii1);
     std::map<unsigned int, std::pair<float, float> >::const_iterator itm4 = 
         fField3DMapCurrentEqualizer.find(ii4);
     const double bix1 = static_cast<double>(itm1->second.first);
     const double biy1 = static_cast<double>(itm1->second.second);     
     const double bix4 = static_cast<double>(itm4->second.first);
     const double biy4 = static_cast<double>(itm4->second.second);
     bxIz1 = (1. - ui)*bix1 + ui*bix4;
     byIz1 = (1. - ui)*biy1 + ui*biy4;
     if (iZ1 == iZ2) {
       Bfield[0] = bxIz1;
       Bfield[1] = byIz1;
     } else { 
       std::map<unsigned int, std::pair<float, float> >::const_iterator jtm1 = 
       fField3DMapCurrentEqualizer.find(jj1);
       std::map<unsigned int, std::pair<float, float> >::const_iterator jtm4 = 
       fField3DMapCurrentEqualizer.find(jj4);
       const double bjx1 = static_cast<double>(jtm1->second.first);
       const double bjy1 = static_cast<double>(jtm1->second.second);     
       const double bjx4 = static_cast<double>(jtm4->second.first);
       const double bjy4 = static_cast<double>(jtm4->second.second);
       bxIz2 = (1. - ui)*bjx1 + ui*bjx4;
       byIz2 = (1. - ui)*bjy1 + ui*bjy4;
       const double v = 
          (Point[2] - fField3DMapCurrentEqualizerZMin 
                    - iZ1*fField3DMapCurrentEqualizerZStep)/fField3DMapCurrentEqualizerZStep;
       Bfield[0] = (1.0 - v)*bxIz1 + v*bxIz2;
       Bfield[1] = (1.0 - v)*byIz1 + v*byIz2;
     }
     if (fDebugIsOn) std::cerr << " ....Phi/Z interpolation only.  " 
                              << " Bfield " << Bfield[0] << " / " << Bfield[1] << std::endl;
   } else {
     if (rSqZ1 < 0.) {
         std::map<unsigned int, std::pair<float, float> >::const_iterator jtm1 = 
          fField3DMapCurrentEqualizer.find(jj1);
         std::map<unsigned int, std::pair<float, float> >::const_iterator jtm2 = 
          fField3DMapCurrentEqualizer.find(jj2);
         std::map<unsigned int, std::pair<float, float> >::const_iterator jtm3 = 
           fField3DMapCurrentEqualizer.find(jj3);
         std::map<unsigned int, std::pair<float, float> >::const_iterator jtm4 = 
           fField3DMapCurrentEqualizer.find(jj4);
         const double bjx1 = static_cast<double>(jtm1->second.first);
         const double bjy1 = static_cast<double>(jtm1->second.second);     
         const double bjx2 = static_cast<double>(jtm2->second.first);
         const double bjy2 = static_cast<double>(jtm2->second.second);     
         const double bjx3 = static_cast<double>(jtm3->second.first);
         const double bjy3 = static_cast<double>(jtm3->second.second);     
         const double bjx4 = static_cast<double>(jtm4->second.first);
         const double bjy4 = static_cast<double>(jtm4->second.second); 
         const double tj = 
          std::max(std::abs((rSqZ2 - radOCAtZ*radOCAtZ)/fField3DMapCurrentEqualizerRadSqRStepAtZ[iZ2]), 1.); // Approximate!. 
         bxIz2 = (1. - ui)*(1.0 - tj)*bjx1 + (1.0 - ui)*tj*bjx2 + ui*tj*bjx3 + (1.0 - tj)*ui*bjx4 ;
         byIz2 = (1. - ui)*(1.0 - tj)*bjy1 + (1.0 - ui)*tj*bjy2 + ui*tj*bjy3 + (1.0 - tj)*ui*bjy4 ;
         if (fDebugIsOn) {
           std::cerr << " tj " << tj << " bjxs " << bjx1 << " / " << bjx2 << " / " << bjx3 << " / " << bjx4 << std::endl;
           std::cerr << " bjys " << bjy1 << " / " << bjy2 << " / " << bjy3 << " / " << bjy4 << std::endl;
           std::cerr << " bxIz2 " << bxIz2 << " byIz2 " << byIz2 << std::endl;
         }
         const double v = 
          (Point[2] - fField3DMapCurrentEqualizerZMin 
                    - iZ1*fField3DMapCurrentEqualizerZStep)/fField3DMapCurrentEqualizerZStep;
         Bfield[0] =  v*bxIz2;
         Bfield[1] =  v*byIz2; 
         if (fDebugIsOn) std::cerr << " ....    v " << v << " Bfield " << Bfield[0] << " / " << Bfield[1] << std::endl;
     } else if (rSqZ2 < 0.) { 
       std::map<unsigned int, std::pair<float, float> >::const_iterator itm1 = 
         fField3DMapCurrentEqualizer.find(ii1);
       std::map<unsigned int, std::pair<float, float> >::const_iterator itm2 = 
         fField3DMapCurrentEqualizer.find(ii2);
       std::map<unsigned int, std::pair<float, float> >::const_iterator itm3 = 
         fField3DMapCurrentEqualizer.find(ii3);
       std::map<unsigned int, std::pair<float, float> >::const_iterator itm4 = 
         fField3DMapCurrentEqualizer.find(ii4);
       const double bix1 = static_cast<double>(itm1->second.first);
       const double biy1 = static_cast<double>(itm1->second.second);     
       const double bix2 = static_cast<double>(itm2->second.first);
       const double biy2 = static_cast<double>(itm2->second.second);     
       const double bix3 = static_cast<double>(itm3->second.first);
       const double biy3 = static_cast<double>(itm3->second.second);     
       const double bix4 = static_cast<double>(itm4->second.first);
       const double biy4 = static_cast<double>(itm4->second.second);
       const double ti =  
        std::max(std::abs((rSqZ1 - radOCAtZ*radOCAtZ)/fField3DMapCurrentEqualizerRadSqRStepAtZ[iZ1]) , 1.);
       bxIz1 = (1. - ui)*(1.0 - ti)*bix1 + (1.0 - ui)*ti*bix2 + ui*ti*bix3 + (1.0 - ti)*ui*bix4 ;
       byIz1 = (1. - ui)*(1.0 - ti)*biy1 + (1.0 - ui)*ti*biy2 + ui*ti*biy3 + (1.0 - ti)*ui*biy4 ;
       if (fDebugIsOn) { 
         std::cerr << " ui " << ui << " ti " << ti << std::endl;
         std::cerr << " bixs " << bix1 << " / " << bix2 << " / " << bix3 << " / " << bix4 << std::endl;
         std::cerr << " biys " << biy1 << " / " << biy2 << " / " << biy3 << " / " << biy4 << std::endl;
         std::cerr << " bxIz1 " << bxIz1 << " byIz1 " << byIz1 << std::endl;
       }
         const double v = 
          (Point[2] - fField3DMapCurrentEqualizerZMin 
                    - iZ1*fField3DMapCurrentEqualizerZStep)/fField3DMapCurrentEqualizerZStep;
       Bfield[0] = (1.0 - v)*bxIz1;
       Bfield[1] = (1.0 - v)*byIz1; 
       if (fDebugIsOn) std::cerr << " ....    v " << v << " Bfield " << Bfield[0] << " / " << Bfield[1] << std::endl;
     
     } else {
       std::map<unsigned int, std::pair<float, float> >::const_iterator itm1 = 
         fField3DMapCurrentEqualizer.find(ii1);
       std::map<unsigned int, std::pair<float, float> >::const_iterator itm2 = 
         fField3DMapCurrentEqualizer.find(ii2);
       std::map<unsigned int, std::pair<float, float> >::const_iterator itm3 = 
         fField3DMapCurrentEqualizer.find(ii3);
       std::map<unsigned int, std::pair<float, float> >::const_iterator itm4 = 
         fField3DMapCurrentEqualizer.find(ii4);
       const double bix1 = static_cast<double>(itm1->second.first);
       const double biy1 = static_cast<double>(itm1->second.second);     
       const double bix2 = static_cast<double>(itm2->second.first);
       const double biy2 = static_cast<double>(itm2->second.second);     
       const double bix3 = static_cast<double>(itm3->second.first);
       const double biy3 = static_cast<double>(itm3->second.second);     
       const double bix4 = static_cast<double>(itm4->second.first);
       const double biy4 = static_cast<double>(itm4->second.second);
       const double ti = 
        (rSqZ1 - iR1Z1*fField3DMapCurrentEqualizerRadSqRStepAtZ[iZ1])/fField3DMapCurrentEqualizerRadSqRStepAtZ[iZ1];
       bxIz1 = (1. - ui)*(1.0 - ti)*bix1 + (1.0 - ui)*ti*bix2 + ui*ti*bix3 + (1.0 - ti)*ui*bix4 ;
       byIz1 = (1. - ui)*(1.0 - ti)*biy1 + (1.0 - ui)*ti*biy2 + ui*ti*biy3 + (1.0 - ti)*ui*biy4 ;
       if (fDebugIsOn) { 
         std::cerr << " ui " << ui << " ti " << ti << std::endl;
         std::cerr << " bixs " << bix1 << " / " << bix2 << " / " << bix3 << " / " << bix4 << std::endl;
         std::cerr << " biys " << biy1 << " / " << biy2 << " / " << biy3 << " / " << biy4 << std::endl;
         std::cerr << " bxIz1 " << bxIz1 << " byIz1 " << byIz1 << std::endl;
       }
       if (iZ1 == iZ2) {
         Bfield[0] = bxIz1;
         Bfield[1] = byIz1;
       } else { 
         std::map<unsigned int, std::pair<float, float> >::const_iterator jtm1 = 
          fField3DMapCurrentEqualizer.find(jj1);
         std::map<unsigned int, std::pair<float, float> >::const_iterator jtm2 = 
          fField3DMapCurrentEqualizer.find(jj2);
         std::map<unsigned int, std::pair<float, float> >::const_iterator jtm3 = 
           fField3DMapCurrentEqualizer.find(jj3);
         std::map<unsigned int, std::pair<float, float> >::const_iterator jtm4 = 
           fField3DMapCurrentEqualizer.find(jj4);
         const double bjx1 = static_cast<double>(jtm1->second.first);
         const double bjy1 = static_cast<double>(jtm1->second.second);     
         const double bjx2 = static_cast<double>(jtm2->second.first);
         const double bjy2 = static_cast<double>(jtm2->second.second);     
         const double bjx3 = static_cast<double>(jtm3->second.first);
         const double bjy3 = static_cast<double>(jtm3->second.second);     
         const double bjx4 = static_cast<double>(jtm4->second.first);
         const double bjy4 = static_cast<double>(jtm4->second.second); 
         const double tj =
          (rSqZ2 - iR1Z2*fField3DMapCurrentEqualizerRadSqRStepAtZ[iZ1])/fField3DMapCurrentEqualizerRadSqRStepAtZ[iZ2];
         bxIz2 = (1. - ui)*(1.0 - tj)*bjx1 + (1.0 - ui)*tj*bjx2 + ui*tj*bjx3 + (1.0 - tj)*ui*bjx4 ;
         byIz2 = (1. - ui)*(1.0 - tj)*bjy1 + (1.0 - ui)*tj*bjy2 + ui*tj*bjy3 + (1.0 - tj)*ui*bjy4 ;
         if (fDebugIsOn) {
           std::cerr << " tj " << tj << " bjxs " << bjx1 << " / " << bjx2 << " / " << bjx3 << " / " << bjx4 << std::endl;
           std::cerr << " bjys " << bjy1 << " / " << bjy2 << " / " << bjy3 << " / " << bjy4 << std::endl;
           std::cerr << " bxIz2 " << bxIz2 << " byIz2 " << byIz2 << std::endl;
         }
         const double v = 
          (Point[2] - fField3DMapCurrentEqualizerZMin 
                    - iZ1*fField3DMapCurrentEqualizerZStep)/fField3DMapCurrentEqualizerZStep;
         Bfield[0] = (1.0 - v)*bxIz1 + v*bxIz2;
         Bfield[1] = (1.0 - v)*byIz1 + v*byIz2;
         if (fDebugIsOn) std::cerr << " ....    v " << v << " Bfield " << Bfield[0] << " / " << Bfield[1] << std::endl;
       }
     } // both rSqZ1 and rSqZ2 positiv  
  }
  if (iQuadrant == 0) return;
  double bfQ[2]; bfQ[0] = Bfield[0]; bfQ[1] = Bfield[1];
  switch (iQuadrant ) {
    case 1:
       Bfield[0] = -bfQ[1]; Bfield[1] = bfQ[0]; break;
    case 2:
       Bfield[0] = -bfQ[0]; Bfield[1] = -bfQ[1]; break;
    case 3:
       Bfield[0] = bfQ[1]; Bfield[1] = -bfQ[0]; break;
  }     
  return;
}
void LBNFMagneticFieldPolyconeHorn::fill3DMapCurrentEqualizer() const {
   //
   if (fFileNameFieldMapForCE != std::string("?")) {
     this->readFieldMapCurrentEqualizer();
     fDebugIsOn = false;
//     this->dumpField();
     return;
   }  
   
   if (fDebugIsOn) std::cerr << " LBNFMagneticFieldPolyconeHorn::fill3DMapCurrentEqualizer, starting .. " << std::endl;
   const G4double current = fHornCurrent/CLHEP::ampere/1000.;
   double magBField = current / (5./CLHEP::cm)/10*CLHEP::tesla; //B(kG)=i(kA)/[5*r(cm)], 1T=10kG 
  // The 1/r dependence is included in the map. 
   fField3DMapCurrentEqualizerRadSqAtZ.clear();
   const double AQuad = fCurrentEqualizerQuadAmpl/100.; // temporary convention.. 
   fNumZPtsForCE = 400; // a bit arbitrary... 
   fNumRPtsForCE = 200;
   fNumPhiPtsForCE = 128;
   fField3DMapCurrentEqualizerZMin = fZDCBegin[0];
   fField3DMapCurrentEqualizerZStep = fEffectiveLength/(fNumZPtsForCE-1);
   fField3DMapCurrentEqualizerPhiStep = (M_PI/2.)/(fNumPhiPtsForCE-1); // Only one quadrant.. 
   const int numWires = 500; // Really arbitrary.. 
   const double phiWStep = 2.0*M_PI/numWires; // Over two pi, full sources. 
   for (unsigned int iZ = 0; iZ != fNumZPtsForCE; iZ++) {
//     std::cerr << " At Z index " << iZ << std::endl;
     const double z = fField3DMapCurrentEqualizerZMin + iZ * fField3DMapCurrentEqualizerZStep;
     double radIC=0.;
     double radOC=0.;   
     this->fillHornPolygonRadii(z, &radIC, &radOC);
     fField3DMapCurrentEqualizerRadSqAtZ.push_back(radOC*radOC);
     double stepRadial = (fOuterRadius*fOuterRadius - radOC*radOC)/(fNumRPtsForCE-1);
     fField3DMapCurrentEqualizerRadSqRStepAtZ.push_back(stepRadial);
     double zFactAzi = 1.0e-6;
     switch(hornNumber) {
        case 0:
        case 2:
          zFactAzi = ((z - fZDCBegin[0]) > 0.) ? (z - fZDCBegin[0])/fCurrentEqualizerLongAbsLength : 1.0e-6;
          break;
        case 1:
          zFactAzi = ((fZDCBegin[fZDCBegin.size() -1] - z) > 0.) ? 
                         (fZDCBegin[fZDCBegin.size() -1] - z)/fCurrentEqualizerLongAbsLength : 1.0e-6; 
          break;
        default:
         zFactAzi = 1.0e-6; // undefined, only 3 horns   
      }
      const double factCE = std::exp(-1.0*std::abs(zFactAzi));
     if (fDebugIsOn) {
         std::cerr << " At z = " << z << "radIC " << radIC << " radOC " 
                               << radOC << " zFactAzi " <<  factCE << " " << magBField << std::endl;
         std::cerr << " Size of the map " << fField3DMapCurrentEqualizer.size() << std::endl;
      }
     for (unsigned int iR=0; iR != fNumRPtsForCE; iR++) {
     // Choose a quadratic spacing... 
       const double r = 1.0e-12 + sqrt(radOC*radOC + iR*stepRadial); // 1 microns to avoid numerics 
       for (unsigned int iPhi=0; iPhi != fNumPhiPtsForCE; iPhi++) {
         const double phi = iPhi*fField3DMapCurrentEqualizerPhiStep;
         const double xP = r * std::cos(phi);
         const double yP = r * std::sin(phi);
         double bX = 0.; double bY = 0.;
         // Discretize the current density on the conductor. Normalize to the unform current density, 
         //  
         double sumCurrent = 0.;
         for (int iW=0;  iW != numWires; iW++) {
           const double phiW =  1.0e-4 + phiWStep*(iW-1);
           const double xW = radOC*std::cos(phiW);
           const double yW = radOC*std::sin(phiW);
           const double xWf = xP - xW;
           const double yWf = yP - yW;
           const double phiff = atan2(yWf, xWf);
           const double radff = std::sqrt(xWf*xWf + yWf*yWf);
           const double iDens = 1.0 + AQuad*factCE*std::abs(sin(2.0*phiW));
           sumCurrent += iDens;
           const double iDensByR = iDens/radff;
           bX -= iDensByR*std::sin(phiff);
           bY += iDensByR*std::cos(phiff);
           if (std::isnan(bX) || std::isnan(bY) || std::isinf(bX) || std::isinf(bY)) {
             std::cerr << " Numerics problem, xP " << xP << " yP " << yP << " iW " 
                       << iW << " xWf " << xWf << " yWf " << yWf << std::endl;
                std::cerr << " And quit here ! " << std::endl;       
             exit(2);       
           }
         }
         if (sumCurrent < 1.0e-10) { bX = 0.; bY = 0.; } else { 
            bX *=magBField/sumCurrent; bY *= magBField/sumCurrent; 
         }
         const unsigned int ii = fNumZPtsForCE*(fNumRPtsForCE*iPhi + iR) + iZ;
         fField3DMapCurrentEqualizer.emplace(ii, 
           std::pair<float, float>(static_cast<float>(bX), static_cast<float>(bY)));
       } // On phi;
    } // on R 
  } // on Z
  /*
   std::cerr << " ... size of final Field map .. " << fField3DMapCurrentEqualizer.size() << std::endl;//
   std::map<unsigned int, std::pair<float, float> >::const_iterator itCheck0 = 
        fField3DMapCurrentEqualizer.find(84029);
   std::cerr << " Check map at location 84029 " << itCheck0->second.first << " / " << itCheck0->second.second << std::endl;
   std::map<unsigned int, std::pair<float, float> >::const_iterator itCheck1 = 
        fField3DMapCurrentEqualizer.find(84009);
   std::cerr << " Check map at location 84009 " << itCheck1->second.first << " / " << itCheck1->second.second << std::endl;
   std::map<unsigned int, std::pair<float, float> >::const_iterator itCheck2 = 
        fField3DMapCurrentEqualizer.find(86009);
   std::cerr << " Check map at location 86009 " << itCheck2->second.first << " / " << itCheck2->second.second << std::endl;
        
  */    
   fDebugIsOn = false;
//   this->dumpField();
   //
   // Write the map out .. 
   //
   this->writeFieldMapCurrentEqualizer();
//   std::cerr << " And quit !. " << std::endl; exit(2);
   
}

void LBNFMagneticFieldPolyconeHorn::writeFieldMapCurrentEqualizer() const {

   std::ostringstream fNameOutStrStr; 
   fNameOutStrStr << "./FieldMapCEQ_Horn_" << (hornNumber +1) 
                  << "_LAbs_" << static_cast<int>(fCurrentEqualizerLongAbsLength) 
                  << "_Quad_" << static_cast<int>(fCurrentEqualizerQuadAmpl) 
                  << "_NZ_" << fNumZPtsForCE << "_NR_" << fNumRPtsForCE 
                  << "_NP_" << fNumPhiPtsForCE << ".dat";
   std::string fNameOutStr(fNameOutStrStr.str()); 
   std::ofstream fOutDat (fNameOutStr.c_str(), std::ios::out | std::ios::binary);
   const char *pTmp1 = (const char*) &fCurrentEqualizerLongAbsLength;
   fOutDat.write(pTmp1, sizeof(double));
   const char *pTmp2 = (const char*) &fCurrentEqualizerQuadAmpl;
   fOutDat.write(pTmp2, sizeof(double));
   char *pTmp = (char*) &fNumZPtsForCE;
   fOutDat.write(pTmp, sizeof(unsigned int));
   pTmp = (char*) &fNumRPtsForCE;
   fOutDat.write(pTmp, sizeof(unsigned int));
   pTmp = (char*) &fNumPhiPtsForCE;
   fOutDat.write(pTmp, sizeof(unsigned int));
   size_t lTotMap =   fField3DMapCurrentEqualizer.size();                 
   pTmp = (char*) &lTotMap;
   fOutDat.write(pTmp, sizeof(size_t));                   
   for (std::map<unsigned int, std::pair<float, float> >::const_iterator it = fField3DMapCurrentEqualizer.begin(); 
                     it != fField3DMapCurrentEqualizer.end(); it++) {
        unsigned int ii = it->first;
        pTmp = (char*) &ii;
        fOutDat.write(pTmp, sizeof(unsigned int));                
        float bb[2];
        bb[0] = it->second.first;
        bb[1] = it->second.second;
        pTmp = (char*) &bb[0];
        fOutDat.write(pTmp, 2*sizeof(float));                     
   }
   fOutDat.close();                  

}

void LBNFMagneticFieldPolyconeHorn::readFieldMapCurrentEqualizer() const {
   
   std::ifstream fInDat(fFileNameFieldMapForCE.c_str(), std::ios::in | std::ios::binary);
   if (!fInDat.is_open()) {
     std::string msg(" LBNFMagneticFieldPolyconeHorn::readFieldMapCurrentEqualizer, could not open file ");
     msg +=  fFileNameFieldMapForCE;
     G4Exception("LBNFMagneticFieldPolyconeHorn::readFieldMapCurrentEqualizer", " ",  RunMustBeAborted, msg.c_str());           
   }
   char *pTmp1 = (char*) &fCurrentEqualizerLongAbsLength;
   fInDat.read(pTmp1, sizeof(double));
   char *pTmp2 = (char*) &fCurrentEqualizerQuadAmpl;
   fInDat.read(pTmp2, sizeof(double));
   char *pTmp = (char*) &fNumZPtsForCE;
   fInDat.read(pTmp, sizeof(unsigned int));
   pTmp = (char*) &fNumRPtsForCE;
   fInDat.read(pTmp, sizeof(unsigned int));
   pTmp = (char*) &fNumPhiPtsForCE;
   fInDat.read(pTmp, sizeof(unsigned int));
   size_t lTotMap = 0;
   pTmp = (char*) &lTotMap;
   fInDat.read(pTmp, sizeof(size_t));
   std::cerr << " LBNFMagneticFieldPolyconeHorn::readFieldMapCurrentEqualizer " 
             << std::endl << " ...... Size of map in file " 
             <<  lTotMap << " Grid " << fNumZPtsForCE << " /  " << fNumRPtsForCE << " / " << fNumPhiPtsForCE << std::endl;                
   for (size_t i=0; i!= lTotMap; i++) {
        unsigned int ii = 0;
        pTmp = (char*) &ii;
        fInDat.read(pTmp, sizeof(unsigned int));                  
        float bb[2];
        pTmp = (char*) &bb[0];
        fInDat.read(pTmp, 2*sizeof(float)); 
        fField3DMapCurrentEqualizer.emplace(ii, std::pair<float, float>(bb[0], bb[1]));                   
   }
   fField3DMapCurrentEqualizerZMin = fZDCBegin[0];
   fField3DMapCurrentEqualizerZStep = fEffectiveLength/(fNumZPtsForCE - 1);
   fField3DMapCurrentEqualizerPhiStep = (M_PI/2.)/(fNumPhiPtsForCE-1); // Only one quadrant.. 
   fField3DMapCurrentEqualizerRadSqRStepAtZ.clear();
   fField3DMapCurrentEqualizerRadSqAtZ.clear();
   for (unsigned int iZ = 0; iZ != fNumZPtsForCE; iZ++) {
//     std::cerr << " At Z index " << iZ << std::endl;
     const double z = fField3DMapCurrentEqualizerZMin + iZ * fField3DMapCurrentEqualizerZStep;
     double radIC=0.;
     double radOC=0.;   
     this->fillHornPolygonRadii(z, &radIC, &radOC);
     fField3DMapCurrentEqualizerRadSqAtZ.push_back(radOC*radOC);
     double stepRadial = (fOuterRadius*fOuterRadius - radOC*radOC)/(fNumRPtsForCE-1);
     fField3DMapCurrentEqualizerRadSqRStepAtZ.push_back(stepRadial);
   }
   fInDat.close();
   this->dumpField();
   std::cerr << " And quit !. " << std::endl; exit(2);
   
} 
void LBNFMagneticFieldPolyconeHorn::TestEccentricityBerkeley() const {
//
// July 25 2017.. 
//
// Alberto M.  suggested the following exercise: Let us take the algorithm in the GetFieldValue method, 
// fill in the radOC and radIC value such that it corresponds to the problem ??? (needs reference )
//
  std::vector<double> Bfield(2, 0.);
  std::vector<double> BfieldB(2, 0.);
  const double radOC = 6.0*CLHEP::mm;
  const double radIC = 2.0*CLHEP::mm;
  const double aBer = radIC;
  const double deltaE = radIC;
  std::ofstream fOut("./TestEccentricityBerkelyV1.txt");
  fOut << " x y Bx By BxB ByB " << std::endl;
  const int numStepY = 50;
  const double stepY = 2*(radIC - 0.01)/numStepY; 
  const int numStepX = 200;
  const double stepX = 3*radOC/numStepX;
  G4double current = fHornCurrent/CLHEP::ampere/1000.;
  const double magBField = current / (5./CLHEP::cm)/10*CLHEP::tesla; //B(kG)=i(kA)/[5*r(cm)], 1T=10kG 
//   double magBField = current / (5.*r/CLHEP::cm)/10*CLHEP::tesla; //B(kG)=i(kA)/[5*r(cm)], 1T=10kG ..
// from GetFieldValue
  const double magBFieldB = 2.0*M_PI*magBField; // adjusted.. Pretty sure the constant above is correct!. 
  for (int iy =0; iy != numStepY+1; iy++) {
    const double y = (iy == numStepY) ? 0. : -radIC + 0.5*stepY + iy*stepY;
    // last loop around is for the test along the X axis 
    for (int ix =0; ix != numStepX; ix++) {
      const double x = -1.5*radOC + 0.5*stepX + ix*stepX;
      const double r = std::sqrt(x*x + y*y);
      const double phi = std::atan2(y, x);
      // cloned from GetFieldValue.. 
      
          // Two Cartesian coordinate system, the first one, centered on Outer surface of the IC 
        // The second one, translated by delta, say to the left, and with critical radius radIn 
//      const double ITot = magBField * r; // Up to a constant.. We don't multiply by r, see above the 
        const double ITot = magBField; // Up to a constant.. We don't multiply by r, see above the 
        // definition of constant magBField.. 
        const double JTot = ITot/(radOC*radOC - radIC*radIC);
        const double BField1Phi = (r < radOC) ? ( JTot * r) : (JTot * radOC*radOC/ r);
        const double BField1X = -1.0*BField1Phi*std::sin(phi); 
        const double BField1Y = BField1Phi*std::cos(phi);
        const double h1 = r * std::sin(phi);
//      const double l2 = fDeltaEccentricityIO + r*std::cos(phi);
        const double l2 = (-1.0*deltaE) + r*std::cos(phi); // sign convention on delta change..
        const double r2 = std::sqrt(l2*l2 + h1*h1); // with respect to the displaced inner radius. 
        const double phi2 = std::atan2(h1, l2);
        const double BField2Phi = (r2 < radIC) ? ( JTot * r2) : (JTot * radIC*radIC/ r2);
        const double BField2X = -1.0*BField2Phi*std::sin(phi2); 
        const double BField2Y = BField2Phi*std::cos(phi2);
        Bfield[0] = BField1X - BField2X;
        Bfield[1] = BField1Y - BField2Y;
      
      // Now the BerKeley soution.. 
      BfieldB[0] = 0.;
      BfieldB[1] = 0.;
      const double r1 = r2; // notation in the Berkeley book 
      if (r1 < radIC) {
        BfieldB[1] = magBFieldB/ (16.0*M_PI*aBer); 
      } else {
        // valid only along the X - axis 
        if (iy == numStepY) {
           if (std::abs(x) > radOC) {
             BfieldB[1] = (8.0*x - 9.0*aBer)*magBFieldB/(16.0*M_PI*x*(x - aBer));
           } else if ((((2.0*aBer <= x) && (x <= radOC))) || (((-radOC <= x) && (x <= 0)))) {
             BfieldB[1] = (x*x -aBer*x -aBer*aBer)*magBFieldB/(16.0*M_PI*aBer*aBer*(x - aBer));
           } else {
             std::cerr << " Unexpected case in LBNFMagneticFieldPolyconeHorn::TestEccentricityBerkeley " << 
                          " x " << x << " y " << y << " r1 " << r1 << " iy " << iy<< std::endl;
                          exit(2);
           }
        }
      }
     fOut << " " << x << " " << y << " " 
          << Bfield[0] << " " << Bfield[1] << " " << BfieldB[0] << " " << BfieldB[1] << std::endl;
    }  
  }
  fOut.close();
  std::cerr << " And quit for now.. " << std::endl; exit(2);

}

LBNFMagneticFieldPolyconeHorn::~LBNFMagneticFieldPolyconeHorn() {
   std::cerr << " Closing Output Trajectory file in LBNFMagneticFieldPolyconeHorn " << std::endl;
   if (fOutTraj.is_open()) fOutTraj.close();
}