diff --git a/NPLib/Detectors/Sofia/GladFieldMap.cxx b/NPLib/Detectors/Sofia/GladFieldMap.cxx
index 3277b93e94d5dbdcd12aae9567884eda24d496cb..c011841489d85fe21e8ecb6763c3c03e9d42d51f 100644
--- a/NPLib/Detectors/Sofia/GladFieldMap.cxx
+++ b/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;
+
}
diff --git a/NPLib/Detectors/Sofia/GladFieldMap.h b/NPLib/Detectors/Sofia/GladFieldMap.h
index d303692eee92b428e6a3a638715e1f1782a07d60..d94501fdcb976ca2cce1b912f24d7e2b7e3e1025 100644
--- a/NPLib/Detectors/Sofia/GladFieldMap.h
+++ b/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);
diff --git a/Projects/s455/Analysis.cxx b/Projects/s455/Analysis.cxx
index efb950f10f929f369143b945ef53c642a7901929..7c33c0b3a883706cf0ef3cce1c207aea7c5ca32d 100644
--- a/Projects/s455/Analysis.cxx
+++ b/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);