LBNFMagneticFieldPolyconeHorn Class Reference

#include <LBNFMagneticFieldPolyconeHorn.hh>

Inheritance diagram for LBNFMagneticFieldPolyconeHorn:

Public Member Functions
	LBNFMagneticFieldPolyconeHorn (size_t iH)
	~LBNFMagneticFieldPolyconeHorn ()
virtual void	GetFieldValue (const double Point[3], double *Bfield) const
void	SetHornCurrent (G4double ihorn)
G4double	GetHornCurrent () const
void	SetSkinDepthInnerRad (G4double f)
void	SetZShiftDrawingCoordinate (double z)
void	SetEffectiveLength (double z)
void	SetDeltaEccentricityIO (G4double f)
void	SetDeltaEllipticityI (G4double f)
double	GetEffectiveInnerRad () const
double	GetEffectiveOuterRad () const
double	GetEffectiveSkinDepthFact () const
void	SetCurrentEqualizerLongAbsLength (double v)
void	SetCurrentEqualizerQuadAmpl (double v)
void	SetCurrentEqualizerOctAmpl (double v)
void	SetFileNameFieldMapForCE (std::string s)
double	GetCurrentEqualizerQuadAmpl () const
std::string	GetFileNameFieldMapForCE () const
void	dumpField () const
void	TestDivergence (int opt) const
void	TestEccentricityBerkeley () const
void	writeFieldMapCurrentEqualizer () const
void	readFieldMapCurrentEqualizer () const

Private Member Functions
void	getFieldFrom3DMapCurrentEqualizer (const double Point[3], double *Bfield) const
void	fill3DMapCurrentEqualizer () const
void	fillHornPolygonRadii (double z, double *rIn, double *rOut) const
void	fillTrajectories (const double Point[3], double bx, double by) const
double	getNormCurrDensityInnerC (double deltaR, double rOut) const

Private Attributes
size_t	hornNumber
G4double	fHornCurrent
std::vector< double >	fZDCBegin
std::vector< double >	fRadsInnerC
std::vector< double >	fRadsOuterC
bool	fHornIsMisaligned
G4RotationMatrix	fRotMatrixHornContainer
double	fEffectiveLength
double	fOuterRadius
double	fOuterRadiusEff
G4double	fSkinDepth
G4double	fSkinDepthInnerRad
G4double	fSkinDepthIIRInv
double	fDeltaEccentricityIO
double	fDeltaEllipticityI
double	fCurrentEqualizerLongAbsLength
double	fCurrentEqualizerQuadAmpl
double	fCurrentEqualizerOctAmpl
std::ofstream	fOutTraj
bool	fDebugIsOn
double	fEffectiveInnerRad
double	fEffectiveOuterRad
double	fEffectiveSkinDepthFact
std::vector< double >	fShift
std::vector< double >	fShiftSlope
G4double	fZShiftDrawingCoordinate
std::string	fFileNameFieldMapForCE
unsigned int	fNumZPtsForCE
unsigned int	fNumRPtsForCE
unsigned int	fNumPhiPtsForCE
double	fField3DMapCurrentEqualizerZMin
double	fField3DMapCurrentEqualizerZStep
double	fField3DMapCurrentEqualizerPhiStep
std::map< unsigned int, std::pair< float, float > >	fField3DMapCurrentEqualizer
std::vector< double >	fField3DMapCurrentEqualizerRadSqAtZ
std::vector< double >	fField3DMapCurrentEqualizerRadSqRStepAtZ

Detailed Description

Definition at line 17 of file LBNFMagneticFieldPolyconeHorn.hh.

Constructor & Destructor Documentation

LBNFMagneticFieldPolyconeHorn::LBNFMagneticFieldPolyconeHorn ( size_t  iH )

explicit

Definition at line 19 of file LBNFMagneticFieldPolyconeHorn.cc.

                                                                      : 
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("?");
}

LBNFMagneticFieldPolyconeHorn::~LBNFMagneticFieldPolyconeHorn

(

)

Definition at line 1195 of file LBNFMagneticFieldPolyconeHorn.cc.

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

Member Function Documentation

void LBNFMagneticFieldPolyconeHorn::dumpField

(

)

const

Definition at line 403 of file LBNFMagneticFieldPolyconeHorn.cc.

                                                    {
//  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::fill3DMapCurrentEqualizer ( ) const

private

Definition at line 908 of file LBNFMagneticFieldPolyconeHorn.cc.

                                                                    {
   //
   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::fillHornPolygonRadii ( double  z, double *  rIn, double *  rOut  ) const

inlineprivate

Definition at line 85 of file LBNFMagneticFieldPolyconeHorn.hh.

                                                                               {
     size_t kLow = 0;
     size_t kHigh = fZDCBegin.size() - 1;
     for (size_t k=0; k != fZDCBegin.size() - 1; k++) {
       if ((z >=  fZDCBegin[k]) && (z < fZDCBegin[k+1])) {
         kLow = k;
         kHigh = k+1;
         const double dzFrac = (z - fZDCBegin[kLow])/(fZDCBegin[kHigh] - fZDCBegin[kLow]);
         *rIn = fRadsInnerC[kLow] + dzFrac*(fRadsInnerC[kHigh] - fRadsInnerC[kLow]);
         *rOut = fRadsOuterC[kLow] + dzFrac*(fRadsOuterC[kHigh] - fRadsOuterC[kLow]);
//       if ((*rIn) < 160.) std::cerr << " .......... z " << z << " k " << k << " rIn " << (*rIn) << " rOut  " << (*rOut) << std::endl;
         return;
       }
     }
     *rIn = 1.0e10;
     *rOut = 2.0e10;
   }

void LBNFMagneticFieldPolyconeHorn::fillTrajectories ( const double  Point[3], double  bx, double  by  ) const

private

Definition at line 575 of file LBNFMagneticFieldPolyconeHorn.cc.

                                                                                                      {
 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();
}

double LBNFMagneticFieldPolyconeHorn::GetCurrentEqualizerQuadAmpl ( ) const

inline

Definition at line 142 of file LBNFMagneticFieldPolyconeHorn.hh.

{return fCurrentEqualizerQuadAmpl;}

double LBNFMagneticFieldPolyconeHorn::GetEffectiveInnerRad ( ) const

inline

Definition at line 134 of file LBNFMagneticFieldPolyconeHorn.hh.

{ return fEffectiveInnerRad;} 

double LBNFMagneticFieldPolyconeHorn::GetEffectiveOuterRad ( ) const

inline

Definition at line 135 of file LBNFMagneticFieldPolyconeHorn.hh.

{ return fEffectiveOuterRad;}

double LBNFMagneticFieldPolyconeHorn::GetEffectiveSkinDepthFact ( ) const

inline

Definition at line 136 of file LBNFMagneticFieldPolyconeHorn.hh.

{ return  fEffectiveSkinDepthFact;}

void LBNFMagneticFieldPolyconeHorn::getFieldFrom3DMapCurrentEqualizer ( const double  Point[3], double *  Bfield  ) const

private

Definition at line 669 of file LBNFMagneticFieldPolyconeHorn.cc.

                                                                           {
  
  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::GetFieldValue ( const double  Point[3], double *  Bfield  ) const

virtual

Definition at line 194 of file LBNFMagneticFieldPolyconeHorn.cc.

{
//   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;
    }
}

std::string LBNFMagneticFieldPolyconeHorn::GetFileNameFieldMapForCE ( ) const

inline

Definition at line 143 of file LBNFMagneticFieldPolyconeHorn.hh.

{ return fFileNameFieldMapForCE;}

G4double LBNFMagneticFieldPolyconeHorn::GetHornCurrent ( ) const

inline

Definition at line 128 of file LBNFMagneticFieldPolyconeHorn.hh.

{ return fHornCurrent;}

double LBNFMagneticFieldPolyconeHorn::getNormCurrDensityInnerC ( double  deltaR, double  rOut  ) const

inlineprivate

Definition at line 111 of file LBNFMagneticFieldPolyconeHorn.hh.

                                                                            {
     if (deltaR < 0.) return 0.;
     const double xby2 = fSkinDepthIIRInv*deltaR/2.; // x must be less < 1, otherwise Taylor expansion fails.. 
     const double xby2Out = fSkinDepthIIRInv*rOut/2.; // x must be less < 1, otherwise Taylor expansion fails.. 
     // keep only 3 terms.. 
     const double ber = 1.0 - std::pow(xby2, 4.)/4. + std::pow(xby2, 8.)/576.;
     const double bei = std::pow(xby2, 2.) - std::pow(xby2, 6.)/36 + std::pow(xby2, 10.)/14400.;
     const double beNSq = ber*ber + bei*bei;
     const double berOut = 1.0 - std::pow(xby2Out, 4.)/4. + std::pow(xby2Out, 8.)/576.;
     const double beiOut = std::pow(xby2Out, 2.) - std::pow(xby2Out, 6.)/36 + std::pow(xby2Out, 10.)/14400.;
     const double beNSqOut = berOut*berOut + beiOut*beiOut;
     return std::sqrt(beNSq/beNSqOut);
   }

void LBNFMagneticFieldPolyconeHorn::readFieldMapCurrentEqualizer

(

)

const

Definition at line 1060 of file LBNFMagneticFieldPolyconeHorn.cc.

                                                                       {
   
   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::SetCurrentEqualizerLongAbsLength ( double  v )

inline

Definition at line 138 of file LBNFMagneticFieldPolyconeHorn.hh.

{fCurrentEqualizerLongAbsLength = v;}

void LBNFMagneticFieldPolyconeHorn::SetCurrentEqualizerOctAmpl ( double  v )

inline

Definition at line 140 of file LBNFMagneticFieldPolyconeHorn.hh.

{fCurrentEqualizerOctAmpl = v;}

void LBNFMagneticFieldPolyconeHorn::SetCurrentEqualizerQuadAmpl ( double  v )

inline

Definition at line 139 of file LBNFMagneticFieldPolyconeHorn.hh.

{fCurrentEqualizerQuadAmpl = v;}

void LBNFMagneticFieldPolyconeHorn::SetDeltaEccentricityIO ( G4double  f )

inline

Definition at line 132 of file LBNFMagneticFieldPolyconeHorn.hh.

{fDeltaEccentricityIO = f;}

void LBNFMagneticFieldPolyconeHorn::SetDeltaEllipticityI ( G4double  f )

inline

Definition at line 133 of file LBNFMagneticFieldPolyconeHorn.hh.

{fDeltaEllipticityI = f;}

void LBNFMagneticFieldPolyconeHorn::SetEffectiveLength ( double  z )

inline

Definition at line 131 of file LBNFMagneticFieldPolyconeHorn.hh.

{ fEffectiveLength = z;}

void LBNFMagneticFieldPolyconeHorn::SetFileNameFieldMapForCE ( std::string  s )

inline

Definition at line 141 of file LBNFMagneticFieldPolyconeHorn.hh.

{ fFileNameFieldMapForCE = s;}

void LBNFMagneticFieldPolyconeHorn::SetHornCurrent ( G4double  ihorn )

inline

Definition at line 127 of file LBNFMagneticFieldPolyconeHorn.hh.

{fHornCurrent = ihorn;}

void LBNFMagneticFieldPolyconeHorn::SetSkinDepthInnerRad ( G4double  f )

inline

Definition at line 129 of file LBNFMagneticFieldPolyconeHorn.hh.

{fSkinDepthInnerRad = f; fSkinDepthIIRInv= std::sqrt(2.0)/f;}

void LBNFMagneticFieldPolyconeHorn::SetZShiftDrawingCoordinate ( double  z )

inline

Definition at line 130 of file LBNFMagneticFieldPolyconeHorn.hh.

{ fZShiftDrawingCoordinate = z;}

void LBNFMagneticFieldPolyconeHorn::TestDivergence

(

int

opt

)

const

Definition at line 585 of file LBNFMagneticFieldPolyconeHorn.cc.

                                                                {
//
// 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::TestEccentricityBerkeley

(

)

const

Definition at line 1113 of file LBNFMagneticFieldPolyconeHorn.cc.

                                                                   {
//
// 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);

}

void LBNFMagneticFieldPolyconeHorn::writeFieldMapCurrentEqualizer

(

)

const

Definition at line 1022 of file LBNFMagneticFieldPolyconeHorn.cc.

                                                                        {

   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();                  

}

Member Data Documentation

double LBNFMagneticFieldPolyconeHorn::fCurrentEqualizerLongAbsLength

private

Definition at line 45 of file LBNFMagneticFieldPolyconeHorn.hh.

double LBNFMagneticFieldPolyconeHorn::fCurrentEqualizerOctAmpl

private

Definition at line 50 of file LBNFMagneticFieldPolyconeHorn.hh.

double LBNFMagneticFieldPolyconeHorn::fCurrentEqualizerQuadAmpl

private

Definition at line 48 of file LBNFMagneticFieldPolyconeHorn.hh.

bool LBNFMagneticFieldPolyconeHorn::fDebugIsOn

mutableprivate

Definition at line 57 of file LBNFMagneticFieldPolyconeHorn.hh.

double LBNFMagneticFieldPolyconeHorn::fDeltaEccentricityIO

private

Definition at line 41 of file LBNFMagneticFieldPolyconeHorn.hh.

double LBNFMagneticFieldPolyconeHorn::fDeltaEllipticityI

private

Definition at line 43 of file LBNFMagneticFieldPolyconeHorn.hh.

double LBNFMagneticFieldPolyconeHorn::fEffectiveInnerRad

mutableprivate

Definition at line 58 of file LBNFMagneticFieldPolyconeHorn.hh.

double LBNFMagneticFieldPolyconeHorn::fEffectiveLength

private

Definition at line 34 of file LBNFMagneticFieldPolyconeHorn.hh.

double LBNFMagneticFieldPolyconeHorn::fEffectiveOuterRad

mutableprivate

Definition at line 59 of file LBNFMagneticFieldPolyconeHorn.hh.

double LBNFMagneticFieldPolyconeHorn::fEffectiveSkinDepthFact

mutableprivate

Definition at line 60 of file LBNFMagneticFieldPolyconeHorn.hh.

std::map<unsigned int, std::pair<float, float> > LBNFMagneticFieldPolyconeHorn::fField3DMapCurrentEqualizer

mutableprivate

Definition at line 77 of file LBNFMagneticFieldPolyconeHorn.hh.

double LBNFMagneticFieldPolyconeHorn::fField3DMapCurrentEqualizerPhiStep

mutableprivate

Definition at line 75 of file LBNFMagneticFieldPolyconeHorn.hh.

std::vector<double> LBNFMagneticFieldPolyconeHorn::fField3DMapCurrentEqualizerRadSqAtZ

mutableprivate

Definition at line 78 of file LBNFMagneticFieldPolyconeHorn.hh.

std::vector<double> LBNFMagneticFieldPolyconeHorn::fField3DMapCurrentEqualizerRadSqRStepAtZ

mutableprivate

Definition at line 79 of file LBNFMagneticFieldPolyconeHorn.hh.

double LBNFMagneticFieldPolyconeHorn::fField3DMapCurrentEqualizerZMin

mutableprivate

Definition at line 73 of file LBNFMagneticFieldPolyconeHorn.hh.

double LBNFMagneticFieldPolyconeHorn::fField3DMapCurrentEqualizerZStep

mutableprivate

Definition at line 74 of file LBNFMagneticFieldPolyconeHorn.hh.

std::string LBNFMagneticFieldPolyconeHorn::fFileNameFieldMapForCE

private

Definition at line 68 of file LBNFMagneticFieldPolyconeHorn.hh.

G4double LBNFMagneticFieldPolyconeHorn::fHornCurrent

private

Definition at line 27 of file LBNFMagneticFieldPolyconeHorn.hh.

bool LBNFMagneticFieldPolyconeHorn::fHornIsMisaligned

private

Definition at line 32 of file LBNFMagneticFieldPolyconeHorn.hh.

unsigned int LBNFMagneticFieldPolyconeHorn::fNumPhiPtsForCE

mutableprivate

Definition at line 71 of file LBNFMagneticFieldPolyconeHorn.hh.

unsigned int LBNFMagneticFieldPolyconeHorn::fNumRPtsForCE

mutableprivate

Definition at line 70 of file LBNFMagneticFieldPolyconeHorn.hh.

unsigned int LBNFMagneticFieldPolyconeHorn::fNumZPtsForCE

mutableprivate

Definition at line 69 of file LBNFMagneticFieldPolyconeHorn.hh.

double LBNFMagneticFieldPolyconeHorn::fOuterRadius

private

Definition at line 35 of file LBNFMagneticFieldPolyconeHorn.hh.

double LBNFMagneticFieldPolyconeHorn::fOuterRadiusEff

private

Definition at line 35 of file LBNFMagneticFieldPolyconeHorn.hh.

std::ofstream LBNFMagneticFieldPolyconeHorn::fOutTraj

mutableprivate

Definition at line 53 of file LBNFMagneticFieldPolyconeHorn.hh.

std::vector<double> LBNFMagneticFieldPolyconeHorn::fRadsInnerC

private

Definition at line 30 of file LBNFMagneticFieldPolyconeHorn.hh.

std::vector<double> LBNFMagneticFieldPolyconeHorn::fRadsOuterC

private

Definition at line 31 of file LBNFMagneticFieldPolyconeHorn.hh.

G4RotationMatrix LBNFMagneticFieldPolyconeHorn::fRotMatrixHornContainer

private

Definition at line 33 of file LBNFMagneticFieldPolyconeHorn.hh.

std::vector<double> LBNFMagneticFieldPolyconeHorn::fShift

private

Definition at line 64 of file LBNFMagneticFieldPolyconeHorn.hh.

std::vector<double> LBNFMagneticFieldPolyconeHorn::fShiftSlope

private

Definition at line 65 of file LBNFMagneticFieldPolyconeHorn.hh.

G4double LBNFMagneticFieldPolyconeHorn::fSkinDepth

private

Definition at line 38 of file LBNFMagneticFieldPolyconeHorn.hh.

G4double LBNFMagneticFieldPolyconeHorn::fSkinDepthIIRInv

private

Definition at line 40 of file LBNFMagneticFieldPolyconeHorn.hh.

G4double LBNFMagneticFieldPolyconeHorn::fSkinDepthInnerRad

private

Definition at line 39 of file LBNFMagneticFieldPolyconeHorn.hh.

std::vector<double> LBNFMagneticFieldPolyconeHorn::fZDCBegin

private

Definition at line 29 of file LBNFMagneticFieldPolyconeHorn.hh.

G4double LBNFMagneticFieldPolyconeHorn::fZShiftDrawingCoordinate

mutableprivate

Definition at line 66 of file LBNFMagneticFieldPolyconeHorn.hh.

size_t LBNFMagneticFieldPolyconeHorn::hornNumber

private

Definition at line 26 of file LBNFMagneticFieldPolyconeHorn.hh.

---

The documentation for this class was generated from the following files:

• g4lbne/include/LBNFMagneticFieldPolyconeHorn.hh
• g4lbne/src/LBNFMagneticFieldPolyconeHorn.cc