ParametrizedEngine/src/BFit3D.cc

0001 #include <cmath>
0002 #ifndef M_PI
0003 #define M_PI 3.1415926535897932385
0004 #endif
0005
0006 #include "BFit3D_data.h"
0007
0008 using namespace magfieldparam;
0009
0010 //_______________________________________________________________________________
0011 void BFit3D::SetCoeff_Linear(const double B) {
0012   unsigned jj, kk = 1;
0013   double w_0, w_1, B_mod = fabs(B);
0014   if (B_mod <= B_nom[0]) {
0015     w_0 = B / B_nom[0];
0016     for (jj = 0; jj < 360; ++jj) {
0017       C[jj] = w_0 * C0[jj][0];
0018     }
0019   } else if (B_mod >= B_nom[3]) {
0020     w_0 = B / B_nom[3];
0021     for (jj = 0; jj < 360; ++jj) {
0022       C[jj] = w_0 * C0[jj][3];
0023     }
0024   } else {
0025     while (B_nom[kk] < B_mod)
0026       ++kk;  //Simple linear search
0027     w_1 = (B_mod - B_nom[kk - 1]) / (B_nom[kk] - B_nom[kk - 1]);
0028     w_0 = 1.0 - w_1;
0029     if (B < 0.) {
0030       w_0 = -w_0;
0031       w_1 = -w_1;
0032     }
0033     for (jj = 0; jj < 360; ++jj) {
0034       C[jj] = w_0 * C0[jj][kk - 1] + w_1 * C0[jj][kk];
0035     }
0036   }
0037 }
0038
0039 //_______________________________________________________________________________
0040 void BFit3D::SetCoeff_Spline(const double B) {
0041   int jc, k0 = 0, k1 = 1;
0042   double dB2, dB = fabs(B);
0043   if (dB >= B_nom[3]) {  //Linear extrapolation for a large field
0044     dB -= B_nom[3];
0045     for (jc = 0; jc < 360; ++jc)
0046       C[jc] = C0[jc][3] + C1[jc][4] * dB;
0047   } else {
0048     if (dB < B_nom[0]) {
0049       dB2 = dB * dB / (3. * B_nom[0]);
0050       for (jc = 0; jc < 360; ++jc)
0051         C[jc] = (C2[jc][0] * dB2 + C1[jc][0]) * dB;
0052     } else {
0053       while (B_nom[k1] < dB)
0054         ++k1;  //Simple linear search
0055       k0 = k1 - 1;
0056       dB2 = (dB -= B_nom[k0]) / (3. * (B_nom[k1] - B_nom[k0]));
0057       if (k1 < 3) {  //"Regular" interval
0058         for (jc = 0; jc < 360; ++jc)
0059           C[jc] = (((C2[jc][k1] - C2[jc][k0]) * dB2 + C2[jc][k0]) * dB + C1[jc][k1]) * dB + C0[jc][k0];
0060       } else {  //The last interval
0061         dB2 = (1. - dB2) * dB;
0062         for (jc = 0; jc < 360; ++jc)
0063           C[jc] = (C2[jc][k0] * dB2 + C1[jc][k1]) * dB + C0[jc][k0];
0064       }
0065     }
0066   }
0067   if (B < 0)
0068     for (jc = 0; jc < 360; ++jc)
0069       C[jc] = -C[jc];
0070 }
0071
0072 //_______________________________________________________________________________
0073 void BFit3D::GetField(const double r, const double z, const double phi, double &Br, double &Bz, double &Bphi) {
0074   //Return field components in Br, Bz, Bphi. SetField must be called before use.
0075   //
0076   if (signed_rad || (r >= 0.))
0077     HB->SetPoint(r, z, phi);
0078   else
0079     HB->SetPoint(-r, z, phi + M_PI);
0080   HB->EvalBr();
0081   HB->EvalBz();
0082   HB->EvalBphi();
0083   Br = HB->GetBr(C);
0084   Bz = HB->GetBz(C);
0085   Bphi = HB->GetBphi(C);
0086 }
0087
0088 /////////////////////////////////////////////////////////////////////////////////
0089 //                                                                             //
0090 //                           T E S T  A R E A                                  //
0091 //                                                                             //
0092 /////////////////////////////////////////////////////////////////////////////////