Updating GladFieldMap

NPLib/Detectors/Sofia/GladFieldMap.cxx

@@ -32,31 +32,31 @@ using namespace std;
 
 ClassImp(GladFieldMap)
 
-
-//////////////////////////////////////////////////////////////////////
-GladFieldMap::GladFieldMap() {
-
-  m_BrhoScan=NULL;
-  m_min = ROOT::Math::Factory::CreateMinimizer("Minuit2","Migrad");
-  m_func = ROOT::Math::Functor(this,&GladFieldMap::Delta,1);
-  m_min->SetFunction(m_func);
-  m_min->SetPrintLevel(-1);
-  m_Bfield = TVector3(0,1.5/1000,0);
-  m_Z_Glad = 4.434*m;
-  m_Leff = 2.*m;
-  m_Tilt = 14.*deg;
-
-  m_Zmax = 9.*m;
-  m_Limit = 1000;
-  m_dt = 0.1*nanosecond;
-
-  m_CentralTheta = 20.*deg;
-  m_X_MWPC3 = -1436.;
-  m_Z_MWPC3 = 8380.;
-
-  m_InitPos = TVector3(0,0,0);
-  m_InitDir = TVector3(0,0,1);
-}
+  //////////////////////////////////////////////////////////////////////
+  GladFieldMap::GladFieldMap() {
+
+    m_BrhoScan=NULL;
+    m_min = ROOT::Math::Factory::CreateMinimizer("Minuit2","Migrad");
+    m_func = ROOT::Math::Functor(this,&GladFieldMap::Delta,1);
+    m_min->SetFunction(m_func);
+    m_min->SetPrintLevel(-1);
+    m_bin = 50;
+    m_Scale = -2135./3583.81;
+    m_Z_Glad = 4.434*m;
+    m_Leff = 2.*m;
+    m_Tilt = 14.*deg;
+
+    m_Zmax = 9.*m;
+    m_Limit = 1000;
+    m_dt = 0.1*nanosecond;
+
+    m_CentralTheta = 20.*deg;
+    m_X_MWPC3 = -1436.;
+    m_Z_MWPC3 = 8380.;
+
+    m_InitPos = TVector3(0,0,0);
+    m_InitDir = TVector3(0,0,1);
+  }
 
 
 
@@ -187,11 +187,12 @@ void GladFieldMap::func(NPL::Particle& N, TVector3 Pos, TVector3 Imp, TVector3&
   xk.SetY(vy);
   xk.SetZ(vz);
 
-  TVector3 B = LoadMap(Pos);
+  vector<double> position = {Pos.X(),Pos.Y(),Pos.Z()};
+  vector<double> B = InterpolateB(position);
   double Bx, By, Bz;
-  Bx = B.X();
-  By = B.Y();
-  Bz = B.Z();
+  Bx = B[0];
+  By = B[1];
+  Bz = B[2];
   double q = N.GetZ()*eplus;
   pk.SetX(q*(vy*Bz - vz*By));
   pk.SetY(q*(vz*Bx - vx*Bz));
@@ -202,7 +203,7 @@ void GladFieldMap::func(NPL::Particle& N, TVector3 Pos, TVector3 Imp, TVector3&
 
 //////////////////////////////////////////////////////////////////////
 TVector3 GladFieldMap::CalculateIntersectionPoint(vector<TVector3> vPos){
- 
+
   unsigned int size = vPos.size();
 
   // Track equation Z = a0*X + b0
@@ -225,22 +226,188 @@ TVector3 GladFieldMap::CalculateIntersectionPoint(vector<TVector3> vPos){
 }
 
 //////////////////////////////////////////////////////////////////////
-TVector3 GladFieldMap::LoadMap(TVector3 Pos) {
-
-  double Bx = 0;
-  double By = 0;
-  double Bz = 0;
-
-  double Zinf = m_Z_Glad + Pos.X()*tan(m_Tilt) - m_Leff/2;
-  double Zsup = m_Z_Glad + Pos.X()*tan(m_Tilt) + m_Leff/2;
-  
-  if(Pos.Z()>Zinf && Pos.Z()<Zsup){
-    Bx = m_Bfield.X();
-    By = m_Bfield.Y();
-    Bz = m_Bfield.Z();
+void GladFieldMap::LoadMap(string filename) {
+  ifstream ifile(filename.c_str());
+  if(!ifile.is_open()){
+    cout << "Error: failed to load Glad field map: " << filename << endl;
+    exit(1);
+  }
+
+  cout << "///////// Loading ASCII Glad field map: " << filename << endl;
+  double x,y,z,Bx,By,Bz;
+
+  unsigned int count=0;
+  while(!ifile.eof()){
+    ifile >> x >> y >> z >> Bx >> By >> Bz;
+    if(++count%10000==0)
+      cout << "\r - Loading " << count << " values " << flush;
+
+    x = x*10;
+    y = y*10;
+    z = z*10+m_Z_Glad;
+
+    Bx *= m_Scale;
+    By *= m_Scale;
+    Bz *= m_Scale;
+
+    vector<double> p = {x,y,z};
+    Bx*=tesla;
+    By*=tesla;
+    Bz*=tesla;
+    vector<double> B = {Bx,By,Bz};
+    m_field[p] = B;
+
+    if(x<m_x_min)
+      m_x_min=x;
+    if(y<m_y_min)
+      m_y_min=y;
+    if(z<m_z_min)
+      m_z_min=z;
+    if(x>m_x_max)
+      m_x_max=x;
+    if(y>m_y_max)
+      m_y_max=y;
+    if(z>m_z_max)
+      m_z_max=z;
+
+  }
+
+  cout << "///////// ASCII file loaded"<< endl;
+
+}
+
+//////////////////////////////////
+vector<double> GladFieldMap::InterpolateB(const vector<double>& pos)
+{
+  static vector<double> nullv = {0,0,0};
+
+  double x,y,z;
+  x=pos[0];
+  y=pos[1];
+  z=pos[2];
+
+  // out of bound
+  if(x<m_x_min || x>m_x_max)
+    return nullv;
+  if(y<m_y_min || y>m_y_max)
+    return nullv;
+  if(z<m_z_min || z>m_z_max)
+    return nullv;
+
+  double x0 = (double)((int)(x)/m_bin)*m_bin;
+  if(x<=x0)
+    x0=(double)((int)(x-m_bin)/m_bin)*m_bin; 
+
+  double x1 = (double)((int)(x)/m_bin)*m_bin;
+  if(x>=x1)
+    x1=(double)((int)(x+m_bin)/m_bin)*m_bin; 
+
+  double y0 = (double)((int)(y)/m_bin)*m_bin;
+  if(y<=y0)
+    y0=(double)((int)(y-m_bin)/m_bin)*m_bin; 
+
+  double y1 = (double)((int)(y)/m_bin)*m_bin;
+  if(y>=y1)
+    y1=(double)((int)(y+m_bin)/m_bin)*m_bin; 
+
+  double z0 = (double)((int)(z)/m_bin)*m_bin;
+  if(z<=z0)
+    z0=(double)((int)(z-m_bin)/m_bin)*m_bin; 
+
+  double z1 = (double)((int)(z)/m_bin)*m_bin;
+  if(z>=z1)
+    z1=(double)((int)(z+m_bin)/m_bin)*m_bin; 
+
+  //vector<double> X={xm,ym,zm};
+  vector<double> X000={x0,y0,z0};
+  vector<double> X111={x1,y1,z1};
+  vector<double> X100={x1,y0,z0};
+  vector<double> X010={x0,y1,z0};
+  vector<double> X001={x0,y0,z1};
+  vector<double> X101={x1,y0,z1};
+  vector<double> X011={x0,y1,z1};
+  vector<double> X110={x1,y1,z0};
+
+  vector<double> C000;
+  vector<double> C111;
+  vector<double> C100;
+  vector<double> C010;
+  vector<double> C001;
+  vector<double> C101;
+  vector<double> C011;
+  vector<double> C110;
+
+  if(m_field.lower_bound(X000)!=m_field.end())
+    C000=m_field.lower_bound(X000)->second;
+  else 
+    return nullv;
+
+  if(m_field.lower_bound(X111)!=m_field.end())
+    C111=m_field.lower_bound(X111)->second;
+  else 
+    return nullv;
+
+  if(m_field.lower_bound(X100)!=m_field.end())
+    C100=m_field.lower_bound(X100)->second;
+  else 
+    return nullv;
+
+  if(m_field.lower_bound(X010)!=m_field.end())
+    C010=m_field.lower_bound(X010)->second;
+  else 
+    return nullv;
+
+  if(m_field.lower_bound(X001)!=m_field.end())
+    C001=m_field.lower_bound(X001)->second;
+  else 
+    return nullv;
+
+  if(m_field.lower_bound(X101)!=m_field.end())
+    C101=m_field.lower_bound(X101)->second;
+  else 
+    return nullv;
+
+  if(m_field.lower_bound(X011)!=m_field.end())
+    C011=m_field.lower_bound(X011)->second;
+  else 
+    return nullv;
+
+  if(m_field.lower_bound(X110)!=m_field.end())
+    C110=m_field.lower_bound(X110)->second;
+  else 
+    return nullv;
+
+  double xd = (x-x0)/(x1-x0);    
+  double yd = (y-y0)/(y1-y0);    
+  double zd = (z-z0)/(z1-z0);    
+  double alphaX = 1-xd;
+  double alphaY = 1-yd;
+  double alphaZ = 1-zd;
+  // X 
+  vector<double> C00 = {C000[0]*alphaX+C100[0]*xd,C000[1]*alphaX+C100[1]*xd,C000[2]*alphaX+C100[2]*xd};
+  vector<double> C01 = {C001[0]*alphaX+C101[0]*xd,C001[1]*alphaX+C101[1]*xd,C001[2]*alphaX+C101[2]*xd};
+  vector<double> C10 = {C010[0]*alphaX+C110[0]*xd,C010[1]*alphaX+C110[1]*xd,C010[2]*alphaX+C110[2]*xd};
+  vector<double> C11 = {C011[0]*alphaX+C111[0]*xd,C011[1]*alphaX+C111[1]*xd,C011[2]*alphaX+C111[2]*xd};
+  // Y
+  vector<double> C0  = {C00[0] *alphaY+C10[0] *yd,C00[1] *alphaY+C10[1] *yd,C00[2] *alphaY+C10[2] *yd};
+  vector<double> C1  = {C01[0] *alphaY+C11[0] *yd,C01[1] *alphaY+C11[1] *yd,C01[2] *alphaY+C11[2] *yd};
+  // Z
+  vector<double> res = {C0[0]  *alphaZ+C1[0]  *zd,C0[1]  *alphaZ+C1[1]  *zd,C0[2]  *alphaZ+C1[2]  *zd};
+
+  return res;
+}
+
+//////////////////////////////////
+vector<double> GladFieldMap::GetB(const vector<double>& pos)
+{
+  static vector<double> nullv = {0,0,0};
+
+  auto it = m_field.find(pos);
+  if(it!=m_field.end()){
+    return it->second;
   }
-  
-  TVector3 B = TVector3(Bx,By,Bz);
-  return B;
+  else
+    return nullv;
+
 }
 

NPLib/Detectors/Sofia/GladFieldMap.h

@@ -24,6 +24,7 @@
 
 // STL
 #include <vector>
+#include <map>
 using namespace std;
 
 // ROOT
@@ -35,6 +36,7 @@ using namespace std;
 
 #include "NPParticle.h"
 
+
 class GladFieldMap{
   public: 
     GladFieldMap();
@@ -42,10 +44,18 @@ class GladFieldMap{
     
   private:
     // GLAD parameters
-    TVector3 m_Bfield;
+    map<vector<double>, vector<double>> m_field;
     double m_Z_Glad;
     double m_Leff;
     double m_Tilt;
+    double m_x_min;
+    double m_y_min;
+    double m_z_min;
+    double m_x_max;
+    double m_y_max;
+    double m_z_max;
+    double m_Scale;
+    int m_bin;
     // MWPC3 paramters
     double m_CentralTheta;
     double m_X_MWPC3;
@@ -61,7 +71,6 @@ class GladFieldMap{
 
 
   public:
-    void SetBfield(TVector3 vec) {m_Bfield = vec;}
     void SetZGlad(double val) {m_Z_Glad = val;}
     void SetLeff(double val) {m_Leff = val;}
     void SetGladTiltAngle(double val) {m_Tilt = val;}
@@ -78,7 +87,9 @@ class GladFieldMap{
     void SetInitDir(TVector3 Dir) {m_InitPos = Dir;}
 
   public:
-    TVector3 LoadMap(TVector3 pos);
+    void LoadMap(string filename);
+    vector<double> InterpolateB(const vector<double>& pos);
+    vector<double> GetB(const vector<double>& pos);
     TGraph* BrhoScan(double Brho_min, double Brho_max, double Brho_step);
     TVector3 CalculateIntersectionPoint(vector<TVector3> vPos);
     vector<TVector3> Propagate(double Brho, TVector3 Pos, TVector3 Dir);

Projects/s455/Analysis.cxx

@@ -77,7 +77,15 @@ void Analysis::Init(){
   SofTofW= (TSofTofWPhysics*) m_DetectorManager->GetDetector("SofTofW");
   SofAt= (TSofAtPhysics*) m_DetectorManager->GetDetector("SofAt");
   SofMwpc= (TSofMwpcPhysics*) m_DetectorManager->GetDetector("SofMwpc");
+  
   m_GladField = new GladFieldMap();
+  m_GladField->SetZGlad(4.434*m);
+  m_GladField->SetLeff(2.067*m);
+  m_GladField->SetGladTiltAngle(14.*deg);
+  m_GladField->SetCentralTheta(20.*deg);
+  m_GladField->SetX_MWPC3(-1.436*m);
+  m_GladField->SetZ_MWPC3(8.380*m);
+  m_GladField->LoadMap("GladFieldMap.dat");
 
   InitParameter();
   InitOutputBranch();
@@ -241,17 +249,17 @@ void Analysis::FissionFragmentAnalysis(){
     if(SofMwpc->DetectorNbr[i]==4){
       if(SofMwpc->PositionX1[i]!=-1000)
         X3.push_back(SofMwpc->PositionX1[i]);
-      
+
       if(SofMwpc->PositionY[i]!=-1000)
         Y3.push_back(SofMwpc->PositionY[i]);
     }
   }
-  
+
 
   for(unsigned int i=0; i<2; i++){
     double tofx = TofHit[i].x;
 
-   for(unsigned int k=0; k<X3.size(); k++){
+    for(unsigned int k=0; k<X3.size(); k++){
       double posx = X3[k];
       double posy = -1000;
       if(Y3.size()==X3.size())
@@ -351,7 +359,7 @@ void Analysis::FissionFragmentAnalysis(){
 
           TofHit[0].theta_in = Theta2;
           TofHit[1].theta_in = Theta1;
- 
+
           TofHit[0].DT = DT2;
           TofHit[1].DT = DT1;
 
@@ -364,7 +372,7 @@ void Analysis::FissionFragmentAnalysis(){
 
           TofHit[0].theta_in = Theta1;
           TofHit[1].theta_in = Theta2;
- 
+
           TofHit[0].DT = DT1;
           TofHit[1].DT = DT2;
 
@@ -381,7 +389,7 @@ void Analysis::FissionFragmentAnalysis(){
 
           TofHit[0].theta_in = Theta3;
           TofHit[1].theta_in = Theta1;
- 
+
           TofHit[0].DT = DT3;
           TofHit[1].DT = DT1;
 
@@ -397,7 +405,7 @@ void Analysis::FissionFragmentAnalysis(){
 
           TofHit[0].theta_in = Theta1;
           TofHit[1].theta_in = Theta3;
- 
+
           TofHit[0].section = 1;
           TofHit[1].section = 3; 
         }
@@ -411,7 +419,7 @@ void Analysis::FissionFragmentAnalysis(){
 
           TofHit[0].theta_in = Theta1;
           TofHit[1].theta_in = Theta4;
-          
+
           TofHit[0].DT = DT1;
           TofHit[1].DT = DT4;
 
@@ -424,7 +432,7 @@ void Analysis::FissionFragmentAnalysis(){
 
           TofHit[0].theta_in = Theta4;
           TofHit[1].theta_in = Theta1;
-           
+
           TofHit[0].DT = DT4;
           TofHit[1].DT = DT1;
 
@@ -441,7 +449,7 @@ void Analysis::FissionFragmentAnalysis(){
 
           TofHit[0].theta_in = Theta2;
           TofHit[1].theta_in = Theta3;
-           
+
           TofHit[0].DT = DT2;
           TofHit[1].DT = DT3;
 
@@ -454,7 +462,7 @@ void Analysis::FissionFragmentAnalysis(){
 
           TofHit[0].theta_in = Theta3;
           TofHit[1].theta_in = Theta2;
- 
+
           TofHit[0].DT = DT3;
           TofHit[1].DT = DT2;
 
@@ -471,7 +479,7 @@ void Analysis::FissionFragmentAnalysis(){
 
           TofHit[0].theta_in = Theta2;
           TofHit[1].theta_in = Theta4;
- 
+
           TofHit[0].DT = DT2;
           TofHit[1].DT = DT4;
 
@@ -484,7 +492,7 @@ void Analysis::FissionFragmentAnalysis(){
 
           TofHit[0].theta_in = Theta4;
           TofHit[1].theta_in = Theta2;
- 
+
           TofHit[0].DT = DT4;
           TofHit[1].DT = DT2;
 
@@ -501,7 +509,7 @@ void Analysis::FissionFragmentAnalysis(){
 
           TofHit[0].theta_in = Theta3;
           TofHit[1].theta_in = Theta4;
- 
+
           TofHit[0].DT = DT3;
           TofHit[1].DT = DT4;
 
@@ -514,7 +522,7 @@ void Analysis::FissionFragmentAnalysis(){
 
           TofHit[0].theta_in = Theta4;
           TofHit[1].theta_in = Theta3;
- 
+
           TofHit[0].DT = DT4;
           TofHit[1].DT = DT3;
 
@@ -529,10 +537,10 @@ void Analysis::FissionFragmentAnalysis(){
         int section = TofHit[i].section;
 
         /*double drift_time;
-        if(section==1) drift_time = DT1;
-        if(section==2) drift_time = DT2;
-        if(section==3) drift_time = DT3;
-        if(section==4) drift_time = DT4;*/
+          if(section==1) drift_time = DT1;
+          if(section==2) drift_time = DT2;
+          if(section==3) drift_time = DT3;
+          if(section==4) drift_time = DT4;*/
 
         double DT_eval;
         if(section<3) DT_eval=55;
@@ -577,7 +585,7 @@ void Analysis::FissionFragmentAnalysis(){
         InitPos[i]  = TVector3(XA,0,ZA);
         InitDir[i]  = TVector3(sin(TofHit[i].theta_in),0,cos(TofHit[i].theta_in));
         FinalPos[i] = TVector3(X3lab,0,Z3lab);
-        
+
         vC    = TVector3(XC,0,ZC);
         vOut  = TVector3(X3lab-XC,0,Z3lab-ZC);
 
@@ -616,23 +624,15 @@ void Analysis::FissionFragmentAnalysis(){
       TVector3 B = TVector3(Bx,By,Bz);
       double Leff = 2.067;
       /*double rho1 = Leff/abs(2*sin(0.5*(TofHit[0].theta_out-TofHit[0].theta_in))*cos(Tilt-0.5*(TofHit[0].theta_out-TofHit[0].theta_in)));
-      double rho2 = Leff/abs(2*sin(0.5*(TofHit[1].theta_out-TofHit[1].theta_in))*cos(Tilt-0.5*(TofHit[1].theta_out-TofHit[1].theta_in)));
-      double Brho1 = MagB*rho1;
-      double Brho2 = MagB*rho2;*/
-
-      m_GladField->SetBfield(B);
-      m_GladField->SetZGlad(4.434*m);
-      m_GladField->SetLeff(2.067*m);
-      m_GladField->SetGladTiltAngle(14.*deg);
-      m_GladField->SetCentralTheta(20.*deg);
-      m_GladField->SetX_MWPC3(-1.436*m);
-      m_GladField->SetZ_MWPC3(8.380*m);
-     
+        double rho2 = Leff/abs(2*sin(0.5*(TofHit[1].theta_out-TofHit[1].theta_in))*cos(Tilt-0.5*(TofHit[1].theta_out-TofHit[1].theta_in)));
+        double Brho1 = MagB*rho1;
+        double Brho2 = MagB*rho2;*/
+
       double Brho1 = 0;
       double Brho2 = 0;
       //Brho1 = m_GladField->FindBrho(InitPos[0], InitDir[0], FinalPos[0]);
       //Brho2 = m_GladField->FindBrho(InitPos[1], InitDir[1], FinalPos[1]);
-      
+
       Beta_Z1 = TofHit[0].beta;
       Beta_Z2 = TofHit[1].beta;
       Gamma1 = 1. / sqrt(1 - Beta_Z1 * Beta_Z1);
@@ -662,7 +662,7 @@ void Analysis::FissionFragmentAnalysis(){
       SofFF->SetPosX3(TofHit[1].x3);
       SofFF->SetPosY3(TofHit[0].y3);
       SofFF->SetPosY3(TofHit[1].y3);
- 
+
       //SofFF->SetPosX3(TofHit[0].x3lab);
       //SofFF->SetPosX3(TofHit[1].x3lab);