LBNFMagneticFieldPolyconeHorn.hh

Go to the documentation of this file.

// $Id: LBNEMagneticField.hh,v 1.2.2.8 2013/12/24 20:30:46 lebrun Exp $
// --------------------------------------------------------------
// LBNEMagneticField.hh modified by Yuki 2004/7/16
// modified by Yuki 8/2/04
#ifndef LBNFMagneticFieldPolyconeHorn_H
#define LBNFMagneticFieldPolyconeHorn_H 1
#include <fstream>
#include <vector>
#include <map>

#include "globals.hh"
#include "G4MagneticField.hh"
#include "G4RotationMatrix.hh"
#include <fstream>

class LBNFMagneticFieldPolyconeHorn : public G4MagneticField
{
  public:
    explicit LBNFMagneticFieldPolyconeHorn(size_t iH);
    ~LBNFMagneticFieldPolyconeHorn();
    
  virtual void GetFieldValue( const double Point[3], double *Bfield ) const;

  private:
   size_t hornNumber; // from 0 to 4, as Horn1 magentic field is still serviced by LBNEMagneticField
   G4double fHornCurrent;
   // data specific for the Polyhorn, and the skin depth effect.  
   std::vector<double> fZDCBegin;
   std::vector<double> fRadsInnerC;
   std::vector<double> fRadsOuterC;
   bool fHornIsMisaligned; // Then, we will rotate the field 
   G4RotationMatrix fRotMatrixHornContainer;
   double fEffectiveLength;
   double fOuterRadius, fOuterRadiusEff; // Not sure we need both of them..
  // Data for the World coordinate transform to the Horn coordinate system. 
   // We linearize the tranformation, leave Z unchanged, except for a translation.
   G4double fSkinDepth;
   G4double fSkinDepthInnerRad;
   G4double  fSkinDepthIIRInv;
   double fDeltaEccentricityIO; // March 2016: the delta parameter that set the eccentricity between the inner and Outer conductors
                                // See Note NuMI-B-471, Feb. 15 1999, page 31
   double fDeltaEllipticityI; // March 2016: the delta parameter that set the ellipticity of the inner conductor
                                // See Note NuMI-B-471, Feb. 15 1999, page 33
   double fCurrentEqualizerLongAbsLength; // the characteristic length at which the azimutal anisotropy decreases. 
                                          // Units: inverse mm. Default value: 100 mm 
  
   double fCurrentEqualizerQuadAmpl;  // Quadrupolar moment amplitude of the azimuthal anisotropy, at Z = zMax. 
                                      // unit: pure number Default value = 0.                                     
   double fCurrentEqualizerOctAmpl;  // Octupole moment amplitude of the azimuthal anisotropy, at Z = zMax. 
                                      // unit: pure number Default value = 0.                                     
  
   mutable std::ofstream fOutTraj;
   //
   // More variable for debugging and checking field values from the stepping Action 
   //
   mutable bool fDebugIsOn;
   mutable double fEffectiveInnerRad;
   mutable double fEffectiveOuterRad; // of the inner conductor. 
   mutable double fEffectiveSkinDepthFact;
   //
   // For possible (albeit futuer, not implemented yet) misalignment studies. 
   //
   std::vector<double> fShift;
   std::vector<double> fShiftSlope;
   mutable G4double fZShiftDrawingCoordinate; // same, might not be needed...
   
   std::string fFileNameFieldMapForCE;
   mutable unsigned int fNumZPtsForCE; 
   mutable unsigned int fNumRPtsForCE; 
   mutable unsigned int fNumPhiPtsForCE;
   // a simple Cartesian grid in Z, R square, and phi. 
   mutable double fField3DMapCurrentEqualizerZMin; 
   mutable double fField3DMapCurrentEqualizerZStep;
   mutable double fField3DMapCurrentEqualizerPhiStep;
   
   mutable std::map<unsigned int, std::pair<float, float> > fField3DMapCurrentEqualizer;
   mutable std::vector<double> fField3DMapCurrentEqualizerRadSqAtZ;   
   mutable std::vector<double> fField3DMapCurrentEqualizerRadSqRStepAtZ;   
   // In global coodinate system. ( but rotation and misalignment taken out.  ) 
   void getFieldFrom3DMapCurrentEqualizer(const double Point[3], double *Bfield) const;
   // No skin depth effect using this map.. 
   void fill3DMapCurrentEqualizer() const ; // pseudo const. 
   
   inline void fillHornPolygonRadii(double z, double *rIn, double *rOut) const {
     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 fillTrajectories(const double Point[3], double bx, double by) const;
   // 
   // Implement the Skin depth for the Inner conductors. 
   // References:  Zarko's Thesis (MINOS Document 5694-v1), section 5.2.4, p150 - 156)
   // DeltaR : difference from the radial coordinate of the track and the inner radius of the inner conductor.
   //
   // Note: Approximate, because this reasoning is based on a cylindical (not tubular) conductor.
   //
   inline double getNormCurrDensityInnerC(double deltaR, double rOut) const {
     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);
   }
   
  public:
    
   inline void SetHornCurrent(G4double ihorn) {fHornCurrent = ihorn;}
   inline G4double GetHornCurrent() const { return fHornCurrent;}
   inline void SetSkinDepthInnerRad(G4double f) {fSkinDepthInnerRad = f; fSkinDepthIIRInv= std::sqrt(2.0)/f;}
   inline void SetZShiftDrawingCoordinate(double z) { fZShiftDrawingCoordinate = z;}
   inline void SetEffectiveLength(double z) { fEffectiveLength = z;}
   inline void SetDeltaEccentricityIO(G4double f) {fDeltaEccentricityIO = f;}
   inline void SetDeltaEllipticityI(G4double f) {fDeltaEllipticityI = f;}
   inline double GetEffectiveInnerRad() const { return fEffectiveInnerRad;} 
   inline double GetEffectiveOuterRad() const { return fEffectiveOuterRad;}
   inline double GetEffectiveSkinDepthFact() const { return  fEffectiveSkinDepthFact;}
   
   inline void SetCurrentEqualizerLongAbsLength(double v) {fCurrentEqualizerLongAbsLength = v;}
   inline void SetCurrentEqualizerQuadAmpl(double v) {fCurrentEqualizerQuadAmpl = v;}
   inline void SetCurrentEqualizerOctAmpl(double v) {fCurrentEqualizerOctAmpl = v;}
   inline void SetFileNameFieldMapForCE(std::string s) { fFileNameFieldMapForCE = s;}
   inline double GetCurrentEqualizerQuadAmpl() const {return fCurrentEqualizerQuadAmpl;}
   inline std::string GetFileNameFieldMapForCE() const { return fFileNameFieldMapForCE;}
  
   void dumpField() const;
   void TestDivergence(int opt) const;
   void TestEccentricityBerkeley() const;
   void writeFieldMapCurrentEqualizer() const;
   void readFieldMapCurrentEqualizer() const; // not really, evidently.. 
};

#endif