diff --git a/NPLib/Detectors/Sofia/GladFieldMap.cxx b/NPLib/Detectors/Sofia/GladFieldMap.cxx
index 8c8719def0ef8f882e33a641c574f885ed257b83..794d309a9c2144148ad7ebec9c972ed0b0595472 100644
--- a/NPLib/Detectors/Sofia/GladFieldMap.cxx
+++ b/NPLib/Detectors/Sofia/GladFieldMap.cxx
@@ -48,9 +48,9 @@ GladFieldMap::GladFieldMap() {
m_bin = 50;
m_Current = 2135.;
m_Scale = m_Current/3583.81;
- m_Glad_Entrance = TVector3(0,0,2774.);
- m_Glad_TurningPoint = TVector3(0,0,2774.+1654);
- m_Tilt = 14.*deg;
+ m_Glad_Entrance = TVector3(0,0,-1.135*m);
+ m_Glad_TurningPoint = TVector3(0,0,0);
+ m_Tilt = -14.*deg;
m_B = m_Scale*m_Bmax;
for(int i=0; i<81; i++){
for(int j=0; j<41; j++){
@@ -71,9 +71,9 @@ GladFieldMap::GladFieldMap() {
m_Nx= 0;
m_Ny= 0;
m_Nz= 0;
- m_CentralTheta = 20.*deg;
- m_MWPC3_POS = TVector3(-1436.,0,7852);
- m_Angle_MWPC3 = 20.*deg;
+ m_CentralTheta = -20.*deg;
+ m_MWPC3_POS = TVector3(-1436.,0,3402);
+ m_Angle_MWPC3 = -20.*deg;
m_R_MWPC3 = 4199.;
m_InitPos = TVector3(0,0,0);
@@ -92,7 +92,7 @@ double GladFieldMap::Delta(const double* parameter){
static vector<TVector3> pos;
pos = Propagate(parameter[0], m_InitPos, m_InitDir, false);
- //pos.back().RotateY(-m_Angle_MWPC3+m_Tilt);
+ pos.back().RotateY(-m_CentralTheta+m_Tilt);
diff = pos.back() - m_FinalPos;
//cout << pos.back().X() << " " << pos.back().Z() << endl;
@@ -103,7 +103,7 @@ double GladFieldMap::Delta(const double* parameter){
double GladFieldMap::FindBrho(TVector3 Pos_init, TVector3 Dir_init, TVector3 Pos_final){
if(!m_BrhoScan)
- BrhoScan(1.,20,0.2,TVector3(0,0,0), TVector3(0,0,1));
+ BrhoScan(6,11,0.1,TVector3(0,0,-2500), TVector3(0,0,1));
m_InitPos = Pos_init;
m_InitDir = Dir_init;
@@ -119,7 +119,7 @@ double GladFieldMap::FindBrho(TVector3 Pos_init, TVector3 Dir_init, TVector3 Pos
m_min->Clear();
m_min->SetPrecision(1e-6);
m_min->SetMaxFunctionCalls(1000);
- m_min->SetLimitedVariable(0,"B",param[0],0.1,1,20);
+ m_min->SetLimitedVariable(0,"B",param[0],0.1,6,11);
m_min->Minimize();
return m_min->X()[0];
@@ -137,11 +137,11 @@ TGraph* GladFieldMap::BrhoScan(double Brho_min, double Brho_max, double Brho_ste
unsigned int i=0;
//TVector3 pos = TVector3(0,0,0);
//TVector3 dir = TVector3(0,0,1);
- //pos.RotateY(m_Tilt);
- //dir.RotateY(m_Tilt);
+ pos.RotateY(m_Tilt);
+ dir.RotateY(m_Tilt);
for(double Brho=Brho_min; Brho<Brho_max; Brho+=Brho_step){
vector<TVector3> vPos = Propagate(Brho, pos, dir, false);
- //vPos.back().RotateY(-m_Angle_MWPC3);
+ vPos.back().RotateY(-m_CentralTheta);
m_BrhoScan->SetPoint(i++,vPos.back().X(),Brho);
//cout << vPos.back().X() << " " << Brho << endl;
diff --git a/Projects/s455/Analysis.cxx b/Projects/s455/Analysis.cxx
index 182eeecd558ecb1c350f22e1e295522033c7f31c..b275e42e599d4715ec43f207a40cb989a8b9b12b 100644
--- a/Projects/s455/Analysis.cxx
+++ b/Projects/s455/Analysis.cxx
@@ -103,14 +103,18 @@ void Analysis::Init(){
m_GladField = new GladFieldMap();
m_GladField->SetCurrent(2135.);
- m_GladField->SetGladEntrance(0, 0.02*m, 2.774*m + 0.5405*m);
- m_GladField->SetGladTurningPoint(0, 0.02*m, 2.774*m + 0.5405*m + 1.135*m);
- m_GladField->SetGladTiltAngle(14.*deg);
- m_GladField->LoadMap("GladFieldMap.dat");
- m_GladField->SetCentralTheta(20.*deg);
+ //m_GladField->SetGladEntrance(0, 0.02*m, 2.774*m + 0.5405*m);
+ m_GladField->SetGladEntrance(0, 0, -1.135*m);
+ //m_GladField->SetGladTurningPoint(0, 0.02*m, 2.774*m + 0.5405*m + 1.135*m);
+ m_GladField->SetGladTurningPoint(0, 0, 0);
+ m_GladField->SetGladTiltAngle(-14.*deg);
+ m_GladField->LoadMap("GladFieldMap_50mm.dat");
+ m_GladField->SetBin(50);
+ m_GladField->SetCentralTheta(-20.*deg);
- double Z_MWPC3 = 7.852*m;
- double X_MWPC3 = -(Z_MWPC3 - m_GladField->GetGladTurningPoint().Z())*tan(m_GladField->GetCentralTheta());
+ //double Z_MWPC3 = 7.852*m;
+ double Z_MWPC3 = 3.402*m;
+ double X_MWPC3 = (Z_MWPC3 - m_GladField->GetGladTurningPoint().Z())*tan(m_GladField->GetCentralTheta());
m_GladField->Set_MWPC3_Position(X_MWPC3,0,Z_MWPC3);
InitParameter();
@@ -586,16 +590,19 @@ void Analysis::FissionFragmentAnalysis(){
// *** Calculation Theta_out *** //
- double Theta0 = m_GladField->GetCentralTheta();
+ double DistanceStartToG = 2.774*m + 0.5405*m +1.135*m;
+ double DistanceStartToA = 2.3155*m;
+ double DistanceStartToMW2 = 2.651*m;
+ double Theta0 = 20.*deg;//m_GladField->GetCentralTheta();
double XA = 0;
- double ZA = 2315.5;
+ double ZA = DistanceStartToA - DistanceStartToG;
double ZG = m_GladField->GetGladTurningPoint().Z();
double ZMW3 = m_GladField->Get_MWPC3_Position().Z();
double XMW3 = m_GladField->Get_MWPC3_Position().X();
- double ZMW2 = 2651;
+ double ZMW2 = DistanceStartToMW2 - DistanceStartToG;
double X3lab = 0;
double Z3lab = 0;;
- double Tilt = m_GladField->GetGladTiltAngle();;//14.*deg;
+ double Tilt = 14.*deg;//;
TVector3 vZ = TVector3(0,0,1);
double XC = 0;
double ZC = 0;
@@ -613,18 +620,19 @@ void Analysis::FissionFragmentAnalysis(){
XC = (XA+(ZG-ZA)*tan(TofHit[i].theta_in)) / (1-tan(Tilt)*tan(TofHit[i].theta_in));
ZC = ZG + XC*tan(Tilt);
TVector3 vC = TVector3(XC,0,ZC);
-
TofHit[i].xc = XC;
TofHit[i].yc = (ZC/8592.)*TofHit[i].y;
TofHit[i].zc = ZC;
int ix, iy;
- ix = (int)(TofHit[0].xc - m_GladField->GetXmin())/50;
- iy = (int)(TofHit[0].yc -20. - m_GladField->GetYmin())/50;
+ ix = (int)(TofHit[0].xc - m_GladField->GetXmin())/m_GladField->GetBin();
+ iy = (int)(TofHit[0].yc - m_GladField->GetYmin())/m_GladField->GetBin();
TofHit[i].Leff = m_GladField->GetLeff(ix,iy);
X3lab = TofHit[i].x3*cos(Theta0) + XMW3;
Z3lab = TofHit[i].x3*sin(Theta0) + ZMW3;
+ TVector3 vE = TVector3(X3lab, 0, Z3lab);
+ TVector3 dir = TVector3(sin(TofHit[i].theta_in), 0, cos(TofHit[i].theta_in));
TofHit[i].x3lab = X3lab;
TofHit[i].z3lab = Z3lab;
TVector3 vOut = TVector3(X3lab-XC,0,Z3lab-ZC);
@@ -632,7 +640,8 @@ void Analysis::FissionFragmentAnalysis(){
TofHit[i].theta_out = angle;
TofHit[i].psi = TofHit[i].theta_in - TofHit[i].theta_out;
TofHit[i].rho = TofHit[i].Leff/(2*sin(0.5*TofHit[i].psi)*cos(Tilt-0.5*TofHit[i].psi));
- TofHit[i].Brho = MagB*TofHit[i].rho;
+ //TofHit[i].Brho = MagB*TofHit[i].rho;
+ TofHit[i].Brho = m_GladField->FindBrho(vA,dir,vE);
// *** Extrapolate to B position *** //
double ZI = ZG - TofHit[i].Leff/(2*cos(Tilt));