LBNEMagneticField.cc

Go to the documentation of this file.

//----------------------------------------------------------------------
//LBNEMagneticField 
// $Id
//----------------------------------------------------------------------

#include "LBNEMagneticField.hh"

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

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

LBNEMagneticFieldHorn::LBNEMagneticFieldHorn(bool isOne) : 
amHorn1(isOne),
fHornCurrent(0.),
fCoordinateSet(false),  
fShift(2, 0.),
fShiftSlope(2,0.),
fZShiftDrawingCoordinate(-1.0e13),
fZShiftUpstrWorldToLocal(-1.0e13),
fEffectiveLength(0.),
fHornNeckOuterRadius(-1.0e13),
fHornNeckInnerRadius(-1.0e13),
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. 

{
 // 
 // We rely on the Volume Placement utility and the surveyor data to compute the coordinate transfers. 
 // Two distinct ones:
 //
 //  (i) Going from G4 World to the local coordinate system.  This is done in two steps: first a Z shift, fixed.
 //         Key variable: fZShiftUpstrWorldToLocal
 // second, a linearized rotation due misalignement, in both X and Y. 
 //
 //  (ii) A Z translation to go from G4 World the 8875.112 MD 363xxx drawing coordinate system, such that 
 //      we can used the exact parabolic equations (Outer radius one)
 //    fZShiftDrawingCoordinate 
 //
 // This strategy has the drawback of the danger of making mistakes in the code below, but has 
 // has the advantage of less coordinate transforms. (inverse one done in the G4 Propagation. )
 // It is assumed here that the misalignment are small, i.e., cos^2(tilt) ~ 0.  
// 
// Having the field manager assigned to the top level of the Horn volume instead of just 
// the conductor polycones volumes, and the volume in between the two conductors, 
//  allows us to add fringe fields if need be, and avoid creating the volumes in between 
// the conductors..
//
  const LBNEVolumePlacements *aPlacementHandler = LBNEVolumePlacements::Instance();
// Radial equation interface.  Take the outer equations..
  
  fEqnIndicesInner.clear();
  fEqnIndicesOuter.clear();
  fZDCBegin.clear();
  fZDCEnd.clear();
  fUseHorn1Polycone = false;
  if(amHorn1) { // Use the Placement data to get the relevant Z coordinate to compute the distance to the inner or outer conductor.  
     const LBNEVolumePlacementData* plDatH1 = aPlacementHandler->Find("Bfield", 
                            "Horn1PolyM1", "LBNEMagneticField::LBNEMagneticField");
     fRotMatrixHornContainer = plDatH1->fRotation;
     const double zCDNeckStart = aPlacementHandler->GetHorn1NeckZPosition() - aPlacementHandler->GetHorn1NeckLength()/2.;
     const double zCDNeckEnd = aPlacementHandler->GetHorn1NeckZPosition() + aPlacementHandler->GetHorn1NeckLength()/2.;
     const double zCDFirstRing = 0.544*zCDNeckStart; // transition from drawing equation 1 to 2 .
                                                   // Drawing 8875.112-MD 363104
//     std::cerr << " Horn1 Magnetized region, check, zCDNeckStart " << zCDNeckStart 
//               <<  " zCDNeckEnd " << zCDNeckEnd << std::endl;
      //
     // No Split geometry, August September 2014. Needs Checking..!!!..... 
     //
    fZDCBegin.push_back(0.); fZDCEnd.push_back(zCDFirstRing); 
     fEqnIndicesInner.push_back(0); fEqnIndicesOuter.push_back(5);
     fZDCBegin.push_back(zCDFirstRing); fZDCEnd.push_back(zCDNeckStart); 
     fEqnIndicesInner.push_back(1); fEqnIndicesOuter.push_back(5);
     fZDCBegin.push_back(zCDNeckStart); fZDCEnd.push_back(zCDNeckEnd);
     fEqnIndicesInner.push_back(99); fEqnIndicesOuter.push_back(99); // No equations, use the radius
     fHornNeckOuterRadius = aPlacementHandler->GetHorn1NeckOuterRadius(); 
     fHornNeckInnerRadius = aPlacementHandler->GetHorn1NeckInnerRadius(); 
     const double zEndNR = aPlacementHandler->GetHorn1ZDEndNeckRegion();
     fZDCBegin.push_back(zCDNeckEnd); fZDCEnd.push_back(zEndNR); 
     fEqnIndicesInner.push_back(3); fEqnIndicesOuter.push_back(7);
     const double zEnd = aPlacementHandler->GetHorn1ZEndIC();
     fZDCBegin.push_back(zEndNR); fZDCEnd.push_back(zEnd);
     fEqnIndicesInner.push_back(4); fEqnIndicesOuter.push_back(7);
     fEffectiveLength = aPlacementHandler->GetHorn1EffectiveLength();
     fOuterRadius = aPlacementHandler->GetHorn1OuterTubeOuterRad();
     fZShiftUpstrWorldToLocal = 0.; // by new definition (beginning of Oct. 2013 ) of the geometry. 
     fZShiftDrawingCoordinate = 30.0*CLHEP::mm*aPlacementHandler->GetHorn1LongRescale(); // Drawing 8875.112-MD-363097
     if (aPlacementHandler->GetUseHorn1Polycone() || (aPlacementHandler->GetNumberOfHornsPolycone() > 0)) {
       fUseHorn1Polycone = true;
       fZShiftDrawingCoordinate = 0.; // Check with Geantino here... 
       fZDCBegin.clear();
       fRadsInnerC.clear();
       fRadsOuterC.clear();
//       std::cerr << " Check coordinates of inner radius of Horn1 Polycone inner conductor... fOuterRadius " << fOuterRadius << std::endl;
       if (aPlacementHandler->IsHornPnParabolic(1)) {
         size_t nPts = aPlacementHandler->GetHornParabolicNumInnerPts(0);
         const double epsil = 0.050*CLHEP::mm; // from LBNEVolumePlacements::UpdateParamsForHornMotherPolyNum
         for (size_t k=0; k !=  nPts; k++) {
           fZDCBegin.push_back(aPlacementHandler->GetHornParabolicRZCoord(0, k) + fZShiftDrawingCoordinate + epsil);
           const double rIn = aPlacementHandler->GetHornParabolicRIn(0, k) - epsil;
           fRadsInnerC.push_back(rIn);
           const double rOut = rIn + aPlacementHandler->GetHornParabolicThick(0, k);
           fRadsOuterC.push_back(rOut);
        }
      } else {
        for (size_t k=0; k!= static_cast<size_t>(aPlacementHandler->GetUseHorn1PolyNumInnerPts()); k++) {
          fZDCBegin.push_back(aPlacementHandler->GetHorn1PolyInnerConductorZCoord(k) + fZShiftDrawingCoordinate);
          const double rIn =  aPlacementHandler->GetHorn1PolyInnerConductorRadius(k);    
          fRadsInnerC.push_back(rIn);
          const double rOut = rIn + aPlacementHandler->GetHorn1PolyInnerConductorThickness(k);
          fRadsOuterC.push_back(rOut);
//       std::cerr << " k " << k << " z " << fZDCBegin[fZDCBegin.size()-1] << " InnerC " << rIn << " out " << rOut << std::endl;
        }
//        std::cerr << " Horn 1 magentic field decalred.. And quit for now.. I say so.. " << std::endl; exit(2);
      }
    } // Polycone, Parabolic style...or straight...   
  } else {
     const  LBNEVolumePlacementData* plDatH2 = aPlacementHandler->Find("Bfield", 
                            "Horn2TopLevel", "LBNEMagneticField::LBNEMagneticField");
//     std::cerr << " Number of equation changes " << aPlacementHandler->GetHorn2ZEqnChangeNumEqn() << std::endl;
//     for (size_t k=0; k != aPlacementHandler->GetHorn2ZEqnChangeNumEqn(); k++) {
//       std::cerr << " Horn2 zone " << k << " at z " 
//                 << aPlacementHandler->GetHorn2ZEqnChange(k)/25.4 << std::endl;
//     }
     fEffectiveLength = plDatH2->fParams[2];
     const double z1 =  aPlacementHandler->GetHorn2ZEqnChange(0);           
     fZDCBegin.push_back(0.); fZDCEnd.push_back(z1); 
     fEqnIndicesInner.push_back(6); fEqnIndicesOuter.push_back(10); //IO trasnsition, Part 1
      
     const double z2 =  aPlacementHandler->GetHorn2ZEqnChange(1);           
     fZDCBegin.push_back(z1); fZDCEnd.push_back(z2); 
     fEqnIndicesInner.push_back(6);  fEqnIndicesOuter.push_back(0); // Getting to the neck region

     const double z3 =  aPlacementHandler->GetHorn2ZEqnChange(2);           
     fZDCBegin.push_back(z2); fZDCEnd.push_back(z3); 
     fEqnIndicesInner.push_back(6);  fEqnIndicesOuter.push_back(0); // Getting to the neck, neck region 

     const double z4 =  aPlacementHandler->GetHorn2ZEqnChange(3);           
     fZDCBegin.push_back(z3); fZDCEnd.push_back(z4); 
     fEqnIndicesInner.push_back(6);  fEqnIndicesOuter.push_back(1); // Just before the neck

     const double z5 =  aPlacementHandler->GetHorn2ZEqnChange(4);           
     fZDCBegin.push_back(z4); fZDCEnd.push_back(z5);     
     fEqnIndicesInner.push_back(99); fEqnIndicesOuter.push_back(99); //The neck, fixe radius
     fHornNeckOuterRadius = aPlacementHandler->GetHorn2NeckOuterRadius();
     fHornNeckInnerRadius = aPlacementHandler->GetHorn2NeckInnerRadius();

     const double z6 =  aPlacementHandler->GetHorn2ZEqnChange(5);           
     fZDCBegin.push_back(z5); fZDCEnd.push_back(z6); 
     fEqnIndicesInner.push_back(7);  fEqnIndicesOuter.push_back(2); // Just after the neck

     const double z7 =  aPlacementHandler->GetHorn2ZEqnChange(6);           
     fZDCBegin.push_back(z6); fZDCEnd.push_back(z7); 
     fEqnIndicesInner.push_back(7);  fEqnIndicesOuter.push_back(3); // moving along positive Z, neck region 

      const double z8 =  aPlacementHandler->GetHorn2ZEqnChange(7);  //  Next physical section     
     fZDCBegin.push_back(z7); fZDCEnd.push_back(z8); 
     fEqnIndicesInner.push_back(7);  fEqnIndicesOuter.push_back(3);

      const double z9 =  aPlacementHandler->GetHorn2ZEqnChange(8);  //      
     fZDCBegin.push_back(z8); fZDCEnd.push_back(z9); 
     fEqnIndicesInner.push_back(8);  fEqnIndicesOuter.push_back(4);

      const double z10 =  aPlacementHandler->GetHorn2ZEqnChange(9);         
     fZDCBegin.push_back(z9); fZDCEnd.push_back(z10); 
     fEqnIndicesInner.push_back(9);  fEqnIndicesOuter.push_back(5);
     
     fOuterRadius = aPlacementHandler->GetHorn2OuterTubeOuterRad();
     fZShiftUpstrWorldToLocal = aPlacementHandler->GetHorn2LongPosition() + 
        aPlacementHandler->GetHorn2LongRescale()*
        aPlacementHandler->GetHorn2DeltaZEntranceToZOrigin(); // by new definition (beginning of Oct. 2013 ) of the geometry. 
     fZShiftDrawingCoordinate = -1.0*aPlacementHandler->GetHorn2LongRescale()*
         aPlacementHandler->GetHorn2DeltaZEntranceToZOrigin(); 
     const  LBNEVolumePlacementData* plDatH2b = aPlacementHandler->Find("Bfield", 
                            "Horn2Hall", "LBNEMagneticField::LBNEMagneticField");
     fRotMatrixHornContainer = plDatH2b->fRotation;
  }
  std::cerr << " Magnetic field coordinate definition done.  fZShiftUpstrWorldToLocal " << 
                    fZShiftUpstrWorldToLocal << " fZShiftDrawingCoordinate " << fZShiftDrawingCoordinate << std::endl;
//  fCoordinateSet = true; obsolete. Except for use of dumping the field... 
  fOuterRadiusEff = fOuterRadius - 2.0*fSkinDepth; // skin depth at 0.43 kHz 
  // Z coordinate change not yet initialize, done at the first track.
/*
  std::cerr << " Table of Z position and equations for ";
  if (amHorn1) std::cerr << " Horn1 " ;
  else  std::cerr << " Horn2 " ;
  std::cerr << std::endl <<  " Z-Start      Z-End      Innner Eqn    Outer Eqn rIC thick" << std::endl;
  for(size_t k=0; k!= fZDCBegin.size(); k++) {
    std::cerr << " " << fZDCBegin[k] << " " << fZDCEnd[k] 
              << " " << fEqnIndicesInner[k] << " " << fEqnIndicesOuter[k];
    const double zMid = 0.5 *( fZDCBegin[k] + fZDCEnd[k]);
    double rOut = fHornNeckOuterRadius; 
    double rIn =  fHornNeckInnerRadius; 
    if (fEqnIndicesInner[k]  != 99) {
      rOut = (amHorn1)  ? aPlacementHandler->GetConductorRadiusHorn1(zMid, fEqnIndicesOuter[k]) :
              aPlacementHandler->GetConductorRadiusHorn2(zMid, fEqnIndicesOuter[k]);
      rIn = (amHorn1)  ? aPlacementHandler->GetConductorRadiusHorn1(zMid, fEqnIndicesInner[k]) :
              aPlacementHandler->GetConductorRadiusHorn2(zMid, fEqnIndicesInner[k]);
    }
    std::cerr << " " << rIn << " " << rOut - rIn  << std::endl;
  }
*/
//
// Use now the survey data (simulated for now.. ) to establish the coordinate transform..  
//
  LBNESurveyor *theSurvey= LBNESurveyor::Instance();
  std::string hornIndexStr = (amHorn1) ? std::string("1") : std::string("2"); 
  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();
  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) {
    if (amHorn1) std::cerr << " Horn1 ";
    else std::cerr << " Horn2 ";
    std::cerr << " is misaligned " 
              <<  " xz term " << fRotMatrixHornContainer[0][2] << " yz " << fRotMatrixHornContainer[1][2] << std::endl;
  }
/*
  if (!amHorn1) {
    std::cerr << " fShift X " << fShift[0] << " Y " << fShift[1] << std::endl;
    std::cerr << " fShiftSlope X " << fShiftSlope[0] << " Y " << fShiftSlope[1] << std::endl;
    std::cerr << " UpLeft X " << pUpLeft[0] << " Y " << pUpLeft[1]  << std::endl;
    std::cerr << " UpRight X " << pUpRight[0] << " Y " << pUpRight[1]  << std::endl;
    std::cerr << " DwLeft X " << pDwLeft[0] << " Y " << pDwLeft[1]  << std::endl;
    std::cerr << " DwRight X " << pDwRight[0] << " Y " << pDwRight[1] << " eff L " << fEffectiveLength << std::endl;
    exit(2); 
  } 
*/
//
// 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;
// }
}

void LBNEMagneticFieldHorn::GetFieldValue(const double Point[3],double *Bfield) const
{

  
   G4double current = fHornCurrent/CLHEP::ampere/1000.;
//   std::cerr << " LBNEMagneticFieldHorn::GetFieldValue Z " << Point[2] << " current " << current << std::endl;
   for (size_t k=0; k!=3; ++k) Bfield[k] = 0.;
   // Intialization of the coordinate transform. 
   const LBNEVolumePlacements *aPlacementHandler = LBNEVolumePlacements::Instance();
   /* Obsolete... 
   if (!fCoordinateSet) {

     G4Navigator* numinavigator=new G4Navigator(); //geometry navigator
     G4Navigator* theNavigator=G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking();
     numinavigator->SetWorldVolume(theNavigator->GetWorldVolume());
     G4ThreeVector Position=G4ThreeVector(Point[0],Point[1],Point[2]); 
     G4VPhysicalVolume* myVolume = numinavigator->LocateGlobalPointAndSetup(Position);
     G4TouchableHistoryHandle aTouchable = numinavigator->CreateTouchableHistoryHandle();
     G4ThreeVector localPosition = aTouchable->GetHistory()->GetTopTransform().TransformPoint(Position);
     delete numinavigator;
     G4String vName = myVolume->GetLogicalVolume()->GetName();
     if (amHorn1 ) {
       if (vName != G4String("Horn1Hall") || (localPosition[2] > 0.)) return; // not there yet.., or getting by the back side.
        fZShiftUpstrWorldToLocal = Point[2]; // to correct for the misalignment, via the first coordinate transform
         // Note: because the logival Horn1 container volume is split into to of them, 
         // the offset for the field might a bit different than for the actual volume. This is a small correction. 
         // What matter most the field is the tilt, which is correct. 
        fZShiftDrawingCoordinate = Point[2] - aPlacementHandler->GetHorn1DeltaZEntranceToZOrigin(); // Check the sign ! 
     } else {
        if (vName != G4String("Horn2TopLevel") || (localPosition[2] > 0.)) return; 
        fZShiftUpstrWorldToLocal = Point[2]; // to correct for the misalignment, via the first coordinate transform
        fZShiftDrawingCoordinate = Point[2] - aPlacementHandler->GetHorn2DeltaZEntranceToZOrigin();
     }
     fCoordinateSet = true;
      this->dumpField();
     std::cerr << " Coordinate transform, Z shifts at Z World  " << Point[2] << " in  " << vName << std::endl;
//      if (!amHorn1) { std::cerr << " And quit  !!!! " << std::endl; exit(2); }
   } // Initialization of Z coordinate transform 
   if (!fCoordinateSet) return;
   */
   std::vector<double> ptTrans(2,0.);
   const double zLocal = Point[2] - fZShiftUpstrWorldToLocal; // from the start of the relevant volume 
   for (size_t k=0; k != 2; k++) 
     ptTrans[k] = Point[k] + fShift[k] + zLocal*fShiftSlope[k]; 
   const double r = std::sqrt(ptTrans[0]*ptTrans[0] + ptTrans[1]*ptTrans[1]); 
   if (r > fOuterRadius) return;
   if ( r < 1.0e-3) return; // The field should go smoothly to zero at the center of the horm 

   double magBField = current / (5.*r/CLHEP::cm)/10*CLHEP::tesla; //B(kG)=i(kA)/[5*r(cm)], 1T=10kG
   double radIC = fHornNeckInnerRadius;
   double radOC = fHornNeckOuterRadius;
   if ( r > fOuterRadiusEff) {
     const double dr = r - fOuterRadiusEff;
     const double attLength  = std::exp(-dr/fSkinDepth);
     magBField *= attLength;
   } else { 
     if (fUseHorn1Polycone) { // Custom Horn1, based on polygons 
       this->fillHorn1PolygonRadii(zLocal, &radIC, &radOC);
       if (std::abs(radIC) > 1.0e6) return;
     } else { // Numi style Horn1 (or Horn2) 
       size_t kSelZ = fEqnIndicesInner.size();
       const double zLocD = zLocal - fZShiftDrawingCoordinate;
       for (size_t k=0; k != fEqnIndicesInner.size(); ++k) {
         if (zLocD < fZDCBegin[k]) break;
         if (zLocD > fZDCEnd[k]) continue; // They are Z ordered.. 
         kSelZ = k; 
         break;
       }
       if (kSelZ ==  fEqnIndicesInner.size()) return;
       if (fEqnIndicesInner[kSelZ] != 99) {
         radIC = (amHorn1) ? aPlacementHandler->GetConductorRadiusHorn1(zLocD,fEqnIndicesInner[kSelZ] ) : 
                        aPlacementHandler->GetConductorRadiusHorn2(zLocD, fEqnIndicesInner[kSelZ]);
         radOC = (amHorn1) ? aPlacementHandler->GetConductorRadiusHorn1(zLocD,fEqnIndicesOuter[kSelZ] ) : 
                        aPlacementHandler->GetConductorRadiusHorn2(zLocD, fEqnIndicesOuter[kSelZ]);
       }
       if (amHorn1 && aPlacementHandler->IsHorn1RadiusBigEnough()) { // Obsolete.. (Aug, Sept. 2014) 
         radIC = std::max(radIC, aPlacementHandler->GetHorn1RMinAllBG()); 
         radOC = std::max(radOC, aPlacementHandler->GetHorn1RMaxAllBG()); 
       }
       if ((radIC < 1.0e-3) || (radOC < 1.0e-3)) {
         std::ostringstream mStrStr;
         mStrStr << " Wrong equation index at Z " << Point[2] << std::endl;
         G4String mStr(mStrStr.str());
         G4Exception("LBNEMagneticFieldHorn::GetFieldValue", " ",  FatalErrorInArgument, mStr.c_str()); 
         return; 
       }
     }
     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) return; // zero field region (Amps law). 
//     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. 
//     }
// Version from Zarko's thesis. ( (MINOS Document 5694-v1) ) 
// 
     if ((r > radIC) && (r < radOC)) {
        const double surfCyl = radOC*radOC - radIC*radIC;
        const double deltaRSq = (r*r - radIC*radIC);
        const double factR = deltaRSq/surfCyl;     
        fEffectiveSkinDepthFact = factR;
       if (fSkinDepthInnerRad > 10.*CLHEP::mm) { // default assume very large skin depth, assume constant current 
                                      // current density .
         magBField *= factR; // Assume uniform current density and apply Amps law. 
       } else { // Realistic skin depth. 
         magBField *= factR*getNormCurrDensityInnerC((r-radIC), (radOC-radIC))
                            /getNormCurrDensityInnerC((radOC-radIC), (radOC-radIC));
       }
     }
// Sin and cos factor bphi -> bx, by 
     Bfield[0] = -magBField*ptTrans[1]/r;
     Bfield[1] = magBField*ptTrans[0]/r;
//
// Important correction.. ( a bug from a long time ago !)
// One needs to rotate the field as well. It is expected to be given in global coordinate  
// 
    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;
    }
     if (!fCoordinateSet) { 
      fCoordinateSet = true; //done once and only once. 
//      this->dumpField();
      }
//     this->fillTrajectories(Point, Bfield[0],  Bfield[1]);
//     std::cerr << " Field region at Z " << Point[2] << " r = " << r 
//         << " radIC " << radIC << " radOC " 
//         << radOC << "zLocD " << zLocD << " magBField " 
//         <<  magBField/tesla << " Bx " << Bfield[0]/tesla << std::endl;
   }
}
void LBNEMagneticFieldHorn::dumpField() const {

  std::string fName = (amHorn1) ? std::string("./FieldMapNuMIHorn1.txt") : 
                                  std::string("./FieldMapNuMIHorn2.txt");
  std::ofstream fOut(fName.c_str());
  fOut << " z zd r bphi " << std::endl;
  double zStart = -500.;
  double zEnd = 4000.;
  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;
  if (!amHorn1) { zStart = 6000.; zEnd = 11000.; }
  const int numZStep = 100;
  const int numRStep= 2000;
  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] - fZShiftDrawingCoordinate;
    double radIC = fHornNeckInnerRadius;
//    size_t kSelZ = fEqnIndicesInner.size();
//    for (size_t k=0; k != fEqnIndicesInner.size(); ++k) {
//       if (zLocD < fZDCBegin[k]) break;
//       if (zLocD > fZDCEnd[k]) continue; // They are Z ordered.. 
//       kSelZ = k; 
//       break;
//    }
//    if (kSelZ ==  fEqnIndicesInner.size()) continue;
//    if (fEqnIndicesInner[kSelZ] != 99) {
//       radIC = (amHorn1) ? aPlacementHandler->GetConductorRadiusHorn1(zLocD,fEqnIndicesInner[kSelZ] ) : 
//                        aPlacementHandler->GetConductorRadiusHorn2(zLocD, fEqnIndicesInner[kSelZ]);
//    }
    double r = std::max(0., (radIC - 5.*CLHEP::mm)); // make sure we cover the Horn1Poly volume
    // inside the conductor... 
    while (r < rMax)  {
      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]);
      fOut << " " << z << "  " << zLocD << " " << r << " " 
           << std::sqrt(bf[0]*bf[0] + bf[1]*bf[1])/CLHEP::tesla << std::endl;
      if (r < (radIC + 5.0*CLHEP::mm)) r += 0.1*CLHEP::mm;
      else  r += rStep;
    }
  } 
  fOut.close();
}  

void LBNEMagneticFieldHorn::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();
}
LBNEMagneticFieldHorn::~LBNEMagneticFieldHorn() {
   std::cerr << " Closing Output Trajectory file in LBNEMagneticFieldHorn " << std::endl;
   if (fOutTraj.is_open()) fOutTraj.close();
}
//
// Decay pipe field.. 
//
LBNEMagneticFieldDecayPipe::LBNEMagneticFieldDecayPipe(bool isToroidal) : 
isToroid(isToroidal),
fWireCurrent(0.),
fWireRadius(0.)
{
}
void LBNEMagneticFieldDecayPipe::GetFieldValue(const double Point[3],double *Bfield) const
{

  
   G4double current = fWireCurrent/CLHEP::ampere/1000.;
   for (size_t k=0; k!=3; ++k) Bfield[k] = 0.;
//
// No misalignement of this wire..
//
   const double r = std::sqrt(Point[0]*Point[0] + Point[1]*Point[1]); 
   if (r < fWireRadius) return;
   
   double magBField = current / (5.*r/CLHEP::cm)/10*CLHEP::tesla; //B(kG)=i(kA)/[5*r(cm)], 1T=10kG
   Bfield[0] = -magBField*Point[1]/r;
   Bfield[1] = magBField*Point[0]/r;
   
//   std::cerr << " LBNEMagneticFieldWire::GetFieldValue Z " << Point[2] 
//             << " current " << current << " r " << r << " magfield " << magBField << std::endl;
   
}