From 8233e66c4b6869e5a82c403dc4bed009bd9c8bbd Mon Sep 17 00:00:00 2001
From: Morfouace <pierre.morfouace@cea.fr>
Date: Thu, 8 Dec 2022 11:14:33 +0100
Subject: [PATCH] Updating Sofia project
---
NPLib/Detectors/Sofia/GladFieldMap.cxx | 7 ++-
NPLib/Detectors/Sofia/GladFieldMap.h | 1 +
NPLib/Detectors/Sofia/TSofFissionFragment.cxx | 2 +
NPLib/Detectors/Sofia/TSofFissionFragment.h | 6 +++
NPLib/Detectors/Sofia/TSofTrimPhysics.cxx | 14 +++---
Projects/s455/Analysis.cxx | 49 ++++++++++++-------
6 files changed, 53 insertions(+), 26 deletions(-)
diff --git a/NPLib/Detectors/Sofia/GladFieldMap.cxx b/NPLib/Detectors/Sofia/GladFieldMap.cxx
index faae6633b..0ffcc2ca3 100644
--- a/NPLib/Detectors/Sofia/GladFieldMap.cxx
+++ b/NPLib/Detectors/Sofia/GladFieldMap.cxx
@@ -44,7 +44,7 @@ GladFieldMap::GladFieldMap() {
m_func = ROOT::Math::Functor(this,&GladFieldMap::Delta,1);
m_min->SetFunction(m_func);
m_min->SetPrintLevel(-1);
- m_Bmax = 2.2;
+ m_Bmax = 2.25;
m_bin = 50;
m_Current = 2135.;
m_Scale = m_Current/3583.81;
@@ -358,8 +358,11 @@ void GladFieldMap::LoadMap(string filename) {
//m_Leff[ix][iy] += abs(By)*m_bin;
// Need to fill this TGraph before scaling the field to get the proper Leff //
- gBy->SetPoint(iz,z,abs(By));
+ gBy->SetPoint(iz,z,By);
+ x += m_Glad_Entrance.X();
+ y += m_Glad_Entrance.Y();
+
z = z + x*sin(m_Tilt);
z += m_Glad_Entrance.Z();
diff --git a/NPLib/Detectors/Sofia/GladFieldMap.h b/NPLib/Detectors/Sofia/GladFieldMap.h
index 462638c7e..0071bee2c 100644
--- a/NPLib/Detectors/Sofia/GladFieldMap.h
+++ b/NPLib/Detectors/Sofia/GladFieldMap.h
@@ -111,6 +111,7 @@ class GladFieldMap{
public:
double GetLeff(int ix, int iy) {return m_Leff[ix][iy];}
+ double GetGladTiltAngle() {return m_Tilt;}
TVector3 GetGladEntrance() {return m_Glad_Entrance;}
TVector3 GetGladTurningPoint() {return m_Glad_TurningPoint;}
double GetB() {return m_B;}
diff --git a/NPLib/Detectors/Sofia/TSofFissionFragment.cxx b/NPLib/Detectors/Sofia/TSofFissionFragment.cxx
index e7bbbb108..f79cc7528 100644
--- a/NPLib/Detectors/Sofia/TSofFissionFragment.cxx
+++ b/NPLib/Detectors/Sofia/TSofFissionFragment.cxx
@@ -52,6 +52,8 @@ void TSofFissionFragment::Clear() {
fFF_TOF.clear();
fFF_Gamma.clear();
fFF_Brho.clear();
+ fFF_Brho_X.clear();
+ fFF_Brho_Z.clear();
fFF_Rho.clear();
fFF_Omega.clear();
fFF_DT.clear();
diff --git a/NPLib/Detectors/Sofia/TSofFissionFragment.h b/NPLib/Detectors/Sofia/TSofFissionFragment.h
index 08943934a..785e7eef5 100644
--- a/NPLib/Detectors/Sofia/TSofFissionFragment.h
+++ b/NPLib/Detectors/Sofia/TSofFissionFragment.h
@@ -43,6 +43,8 @@ class TSofFissionFragment : public TObject {
vector<double> fFF_TOF;
vector<double> fFF_Gamma;
vector<double> fFF_Brho;
+ vector<double> fFF_Brho_X;
+ vector<double> fFF_Brho_Z;
vector<double> fFF_Rho;
vector<double> fFF_Omega;
vector<double> fFF_DT;
@@ -103,6 +105,8 @@ class TSofFissionFragment : public TObject {
inline void SetTOF(double val){fFF_TOF.push_back(val);};//!
inline void SetGamma(double val){fFF_Gamma.push_back(val);};//!
inline void SetBrho(double val){fFF_Brho.push_back(val);};//!
+ inline void SetBrhoX(double val){fFF_Brho_X.push_back(val);};//!
+ inline void SetBrhoZ(double val){fFF_Brho_Z.push_back(val);};//!
inline void SetRho(double val){fFF_Rho.push_back(val);};//!
inline void SetOmega(double val){fFF_Omega.push_back(val);};//!
inline void SetDT(double val){fFF_DT.push_back(val);};//!
@@ -146,6 +150,8 @@ class TSofFissionFragment : public TObject {
inline double GetTOF(int i) const {return fFF_TOF[i];}//!
inline double GetGamma(int i) const {return fFF_Gamma[i];}//!
inline double GetBrho(int i) const {return fFF_Brho[i];}//!
+ inline double GetBrhoX(int i) const {return fFF_Brho_X[i];}//!
+ inline double GetBrhoZ(int i) const {return fFF_Brho_Z[i];}//!
inline double GetRho(int i) const {return fFF_Rho[i];}//!
inline double GetOmega(int i) const {return fFF_Omega[i];}//!
inline double GetDT(int i) const {return fFF_DT[i];}//!
diff --git a/NPLib/Detectors/Sofia/TSofTrimPhysics.cxx b/NPLib/Detectors/Sofia/TSofTrimPhysics.cxx
index b2f606d34..e47aa40fd 100644
--- a/NPLib/Detectors/Sofia/TSofTrimPhysics.cxx
+++ b/NPLib/Detectors/Sofia/TSofTrimPhysics.cxx
@@ -188,12 +188,14 @@ void TSofTrimPhysics::BuildPhysicalEvent() {
p1_2[i] = Cal->GetValue("SofTrim/SEC"+NPL::itoa(i+1)+"_ANODE2_BETA",1);
p1_3[i] = Cal->GetValue("SofTrim/SEC"+NPL::itoa(i+1)+"_ANODE3_BETA",1);
- double norm1 = p0_1[i] + p1_1[i]*TMath::Power(m_BetaNorm, -5./3.);
- double norm2 = p0_2[i] + p1_2[i]*TMath::Power(m_BetaNorm, -5./3.);
- double norm3 = p0_3[i] + p1_3[i]*TMath::Power(m_BetaNorm, -5./3.);
- Ep1[i] = norm1 * Ep1[i] / (p0_1[i] + p1_1[i]*TMath::Power(m_Beta, -5./3.));
- Ep2[i] = norm2 * Ep2[i] / (p0_2[i] + p1_2[i]*TMath::Power(m_Beta, -5./3.));
- Ep3[i] = norm3 * Ep3[i] / (p0_3[i] + p1_3[i]*TMath::Power(m_Beta, -5./3.));
+ if(p0_1[i]!=0 && p0_2[i]!=0 && p0_3[i]!=0 && p1_1[i]!=0 && p1_2[i]!=0 && p1_3[i]!=0){
+ double norm1 = p0_1[i] + p1_1[i]*TMath::Power(m_BetaNorm, -5./3.);
+ double norm2 = p0_2[i] + p1_2[i]*TMath::Power(m_BetaNorm, -5./3.);
+ double norm3 = p0_3[i] + p1_3[i]*TMath::Power(m_BetaNorm, -5./3.);
+ Ep1[i] = norm1 * Ep1[i] / (p0_1[i] + p1_1[i]*TMath::Power(m_Beta, -5./3.));
+ Ep2[i] = norm2 * Ep2[i] / (p0_2[i] + p1_2[i]*TMath::Power(m_Beta, -5./3.));
+ Ep3[i] = norm3 * Ep3[i] / (p0_3[i] + p1_3[i]*TMath::Power(m_Beta, -5./3.));
+ }
// Angle correction per pair: spline
/*Ep1[i] = Ep1[i] / fcorr_EvsA[i][0]->Eval(Ddt) * fcorr_EvsA[i][0]->Eval(0);
diff --git a/Projects/s455/Analysis.cxx b/Projects/s455/Analysis.cxx
index 547849d13..182eeecd5 100644
--- a/Projects/s455/Analysis.cxx
+++ b/Projects/s455/Analysis.cxx
@@ -75,6 +75,8 @@ struct TofPair
double Leff = 0;
double rho = 0;
double Brho = 0;
+ double BrhoX = 0;
+ double BrhoZ = 0;
double omega = 0;
double deff1 = 0;
double deff2 = 0;
@@ -101,8 +103,8 @@ void Analysis::Init(){
m_GladField = new GladFieldMap();
m_GladField->SetCurrent(2135.);
- m_GladField->SetGladEntrance(0,0,2.774*m);
- m_GladField->SetGladTurningPoint(0,0,m_GladField->GetGladEntrance().Z() + 1.654*m);
+ 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);
@@ -587,16 +589,13 @@ void Analysis::FissionFragmentAnalysis(){
double Theta0 = m_GladField->GetCentralTheta();
double XA = 0;
double ZA = 2315.5;
- int ix = (int) (-m_GladField->GetXmin()/50);
- int iy = (int) (-m_GladField->GetYmin()/50);
- double Leff_init = m_GladField->GetLeff(ix,iy);
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 X3lab = 0;
double Z3lab = 0;;
- double Tilt = 14.*deg;
+ double Tilt = m_GladField->GetGladTiltAngle();;//14.*deg;
TVector3 vZ = TVector3(0,0,1);
double XC = 0;
double ZC = 0;
@@ -613,22 +612,22 @@ void Analysis::FissionFragmentAnalysis(){
// *** Extroplate to C position *** //
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 = 0;//TofHit[i].y*0.5;
+ 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 - m_GladField->GetYmin())/50;
+ iy = (int)(TofHit[0].yc -20. - m_GladField->GetYmin())/50;
TofHit[i].Leff = m_GladField->GetLeff(ix,iy);
X3lab = TofHit[i].x3*cos(Theta0) + XMW3;
Z3lab = TofHit[i].x3*sin(Theta0) + ZMW3;
TofHit[i].x3lab = X3lab;
TofHit[i].z3lab = Z3lab;
- TVector3 vC = TVector3(TofHit[i].xc,0,TofHit[i].zc);
- TVector3 vOut = TVector3(X3lab-TofHit[i].xc,0,Z3lab-TofHit[i].zc);
+ TVector3 vOut = TVector3(X3lab-XC,0,Z3lab-ZC);
double angle = -vZ.Angle(vOut);
TofHit[i].theta_out = angle;
TofHit[i].psi = TofHit[i].theta_in - TofHit[i].theta_out;
@@ -646,15 +645,21 @@ void Analysis::FissionFragmentAnalysis(){
// *** Extrapolate to D position *** //
//XD = XC + TofHit[i].Leff/2*tan(angle+Tilt)*cos(Tilt);
//ZD = ZC + TofHit[i].Leff/2*cos(Tilt);
- double phi = -angle -Tilt;
- double psi = TMath::Pi()/2+angle;
+ double phi = abs(angle) - Tilt;
+ double psi = TMath::Pi()/2-abs(angle);
double l = TofHit[i].Leff/(2*cos(phi));
XD = XC - l*cos(psi);
ZD = ZC + l*sin(psi);
TofHit[i].xd = XD;
TofHit[i].zd = ZD;
TVector3 vD = TVector3(XD,0,ZD);
-
+
+ // *** Extrapolate position of the rho point *** //
+ TofHit[i].BrhoX = XB - TofHit[i].rho*cos(TofHit[i].theta_in);
+ TofHit[i].BrhoZ = ZB + TofHit[i].rho*sin(TofHit[i].theta_in);
+ TVector3 vBrho = TVector3(TofHit[i].BrhoX, 0, TofHit[i].BrhoZ);
+
+
TVector3 v3lab = TVector3(X3lab,0,Z3lab);
TVector3 v1 = TVector3(XB,0,ZB);
TVector3 v3 = TVector3(X3lab-XD,0,Z3lab-ZD);
@@ -679,11 +684,17 @@ void Analysis::FissionFragmentAnalysis(){
TVector3 vCG = vG - vC;
TVector3 vCA = vA - vC;
- TofHit[i].deff1 = vCG.Angle(vCA)*180./TMath::Pi();
- TofHit[i].deff2 = vCG.Angle(vOut)*180./TMath::Pi();
+ //TofHit[i].deff1 = vCG.Angle(vCA)*180./TMath::Pi();
+ //TofHit[i].deff2 = vCG.Angle(vOut)*180./TMath::Pi();
+
+ TVector3 vBD = vD-vB;
+ double vBD_angle = vZ.Angle(vBD);
- //TofHit[i].deff1 = (vC-vB).Mag();
- //TofHit[i].deff2 = (vC-vD).Mag();
+ double XO = XB - TofHit[i].Leff/(2*cos(vBD_angle - Tilt))*cos(TMath::Pi()/2-vBD_angle);
+ double ZO = ZB + TofHit[i].Leff/(2*cos(vBD_angle - Tilt))*sin(TMath::Pi()/2-vBD_angle);
+ TVector3 vO = TVector3(XO,0,ZO);
+ TofHit[i].deff1 = (vO-vB).Mag();
+ TofHit[i].deff2 = (vD-vO).Mag();
TofHit[i].gamma = 1. / sqrt(1 - pow(TofHit[i].beta,2));
TofHit[i].AoQ = TofHit[i].Brho / (3.10716 * TofHit[i].beta * TofHit[i].gamma);
@@ -717,6 +728,8 @@ void Analysis::FissionFragmentAnalysis(){
SofFF->SetAoQ(TofHit[i].AoQ);
SofFF->SetA(TofHit[i].A);
SofFF->SetBrho(TofHit[i].Brho);
+ SofFF->SetBrhoX(TofHit[i].BrhoX);
+ SofFF->SetBrhoZ(TofHit[i].BrhoZ);
SofFF->SetRho(TofHit[i].rho);
SofFF->SetOmega(TofHit[i].omega);
SofFF->SetDT(TofHit[i].DT);
--
GitLab