LBNEMagneticFieldHorn Class Reference

#include <LBNEMagneticField.hh>

Inheritance diagram for LBNEMagneticFieldHorn:

Public Member Functions
	LBNEMagneticFieldHorn (bool isHorn1)
	~LBNEMagneticFieldHorn ()
virtual void	GetFieldValue (const double Point[3], double *Bfield) const
void	SetHornCurrent (G4double ihorn)
G4double	GetHornCurrent () const
void	SetSkinDepthInnerRad (G4double f)
double	GetEffectiveInnerRad () const
double	GetEffectiveOuterRad () const
double	GetEffectiveSkinDepthFact () const
void	dumpField () const

Private Member Functions
void	fillHorn1PolygonRadii (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
G4bool	amHorn1
G4double	fHornCurrent
G4bool	fCoordinateSet
bool	fHornIsMisaligned
G4RotationMatrix	fRotMatrixHornContainer
std::vector< double >	fShift
std::vector< double >	fShiftSlope
G4double	fZShiftDrawingCoordinate
G4double	fZShiftUpstrWorldToLocal
G4double	fEffectiveLength
G4double	fHornNeckOuterRadius
G4double	fHornNeckInnerRadius
G4double	fOuterRadius
G4double	fOuterRadiusEff
G4double	fSkinDepth
std::vector< size_t >	fEqnIndicesOuter
std::vector< size_t >	fEqnIndicesInner
std::vector< double >	fZDCBegin
std::vector< double >	fZDCEnd
G4double	fNeckRadius
G4double	fSkinDepthInnerRad
G4double	fSkinDepthIIRInv
std::ofstream	fOutTraj
double	fEffectiveInnerRad
double	fEffectiveOuterRad
double	fEffectiveSkinDepthFact
bool	fUseHorn1Polycone
std::vector< double >	fRadsInnerC
std::vector< double >	fRadsOuterC

Detailed Description

Definition at line 17 of file LBNEMagneticField.hh.

Constructor & Destructor Documentation

LBNEMagneticFieldHorn::LBNEMagneticFieldHorn ( bool  isHorn1 )

explicit

Definition at line 18 of file LBNEMagneticField.cc.

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

LBNEMagneticFieldHorn::~LBNEMagneticFieldHorn

(

)

Definition at line 472 of file LBNEMagneticField.cc.

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

Member Function Documentation

void LBNEMagneticFieldHorn::dumpField

(

)

const

Definition at line 412 of file LBNEMagneticField.cc.

                                            {

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

inlineprivate

Definition at line 77 of file LBNEMagneticField.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]);
//       std::cerr << " .......... z " << z << " k " << k << " rIn " << (*rIn) << " rOut  " << (*rOut) << std::endl;
         return;
       }
     }
     *rIn = 1.0e10;
     *rOut = 2.0e10;
   }

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

private

Definition at line 463 of file LBNEMagneticField.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 LBNEMagneticFieldHorn::GetEffectiveInnerRad ( ) const

inline

Definition at line 123 of file LBNEMagneticField.hh.

{ return fEffectiveInnerRad;} 

double LBNEMagneticFieldHorn::GetEffectiveOuterRad ( ) const

inline

Definition at line 124 of file LBNEMagneticField.hh.

{ return fEffectiveOuterRad;}

double LBNEMagneticFieldHorn::GetEffectiveSkinDepthFact ( ) const

inline

Definition at line 125 of file LBNEMagneticField.hh.

{ return  fEffectiveSkinDepthFact;}

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

virtual

Definition at line 266 of file LBNEMagneticField.cc.

{

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

G4double LBNEMagneticFieldHorn::GetHornCurrent ( ) const

inline

Definition at line 121 of file LBNEMagneticField.hh.

{ return fHornCurrent;}

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

inlineprivate

Definition at line 104 of file LBNEMagneticField.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 LBNEMagneticFieldHorn::SetHornCurrent ( G4double  ihorn )

inline

Definition at line 120 of file LBNEMagneticField.hh.

{fHornCurrent = ihorn;}

void LBNEMagneticFieldHorn::SetSkinDepthInnerRad ( G4double  f )

inline

Definition at line 122 of file LBNEMagneticField.hh.

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

Member Data Documentation

G4bool LBNEMagneticFieldHorn::amHorn1

private

Definition at line 26 of file LBNEMagneticField.hh.

G4bool LBNEMagneticFieldHorn::fCoordinateSet

mutableprivate

Definition at line 30 of file LBNEMagneticField.hh.

double LBNEMagneticFieldHorn::fEffectiveInnerRad

mutableprivate

Definition at line 65 of file LBNEMagneticField.hh.

G4double LBNEMagneticFieldHorn::fEffectiveLength

private

Definition at line 43 of file LBNEMagneticField.hh.

double LBNEMagneticFieldHorn::fEffectiveOuterRad

mutableprivate

Definition at line 66 of file LBNEMagneticField.hh.

double LBNEMagneticFieldHorn::fEffectiveSkinDepthFact

mutableprivate

Definition at line 67 of file LBNEMagneticField.hh.

std::vector<size_t> LBNEMagneticFieldHorn::fEqnIndicesInner

private

Definition at line 52 of file LBNEMagneticField.hh.

std::vector<size_t> LBNEMagneticFieldHorn::fEqnIndicesOuter

private

Definition at line 51 of file LBNEMagneticField.hh.

G4double LBNEMagneticFieldHorn::fHornCurrent

private

Definition at line 27 of file LBNEMagneticField.hh.

bool LBNEMagneticFieldHorn::fHornIsMisaligned

private

Definition at line 31 of file LBNEMagneticField.hh.

G4double LBNEMagneticFieldHorn::fHornNeckInnerRadius

private

Definition at line 45 of file LBNEMagneticField.hh.

G4double LBNEMagneticFieldHorn::fHornNeckOuterRadius

private

Definition at line 44 of file LBNEMagneticField.hh.

G4double LBNEMagneticFieldHorn::fNeckRadius

private

Definition at line 55 of file LBNEMagneticField.hh.

G4double LBNEMagneticFieldHorn::fOuterRadius

private

Definition at line 46 of file LBNEMagneticField.hh.

G4double LBNEMagneticFieldHorn::fOuterRadiusEff

private

Definition at line 47 of file LBNEMagneticField.hh.

std::ofstream LBNEMagneticFieldHorn::fOutTraj

mutableprivate

Definition at line 61 of file LBNEMagneticField.hh.

std::vector<double> LBNEMagneticFieldHorn::fRadsInnerC

private

Definition at line 74 of file LBNEMagneticField.hh.

std::vector<double> LBNEMagneticFieldHorn::fRadsOuterC

private

Definition at line 75 of file LBNEMagneticField.hh.

G4RotationMatrix LBNEMagneticFieldHorn::fRotMatrixHornContainer

private

Definition at line 32 of file LBNEMagneticField.hh.

std::vector<double> LBNEMagneticFieldHorn::fShift

private

Definition at line 33 of file LBNEMagneticField.hh.

std::vector<double> LBNEMagneticFieldHorn::fShiftSlope

private

Definition at line 34 of file LBNEMagneticField.hh.

G4double LBNEMagneticFieldHorn::fSkinDepth

private

Definition at line 48 of file LBNEMagneticField.hh.

G4double LBNEMagneticFieldHorn::fSkinDepthIIRInv

private

Definition at line 60 of file LBNEMagneticField.hh.

G4double LBNEMagneticFieldHorn::fSkinDepthInnerRad

private

Definition at line 59 of file LBNEMagneticField.hh.

bool LBNEMagneticFieldHorn::fUseHorn1Polycone

private

Definition at line 73 of file LBNEMagneticField.hh.

std::vector<double> LBNEMagneticFieldHorn::fZDCBegin

private

Definition at line 53 of file LBNEMagneticField.hh.

std::vector<double> LBNEMagneticFieldHorn::fZDCEnd

private

Definition at line 54 of file LBNEMagneticField.hh.

G4double LBNEMagneticFieldHorn::fZShiftDrawingCoordinate

mutableprivate

Definition at line 35 of file LBNEMagneticField.hh.

G4double LBNEMagneticFieldHorn::fZShiftUpstrWorldToLocal

mutableprivate

Definition at line 38 of file LBNEMagneticField.hh.

---

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

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