From ba668fe760e6a857712377798888f2d403e33557 Mon Sep 17 00:00:00 2001
From: Morfouace <pierre.morfouace@cea.fr>
Date: Fri, 27 Jan 2023 17:57:05 +0100
Subject: [PATCH] Updating sofia project
---
Projects/Sofia/Analysis.cxx | 8 +++++---
Projects/s455/Analysis.cxx | 31 +++++++++++++++++++------------
2 files changed, 24 insertions(+), 15 deletions(-)
diff --git a/Projects/Sofia/Analysis.cxx b/Projects/Sofia/Analysis.cxx
index 0bc490065..452ca17e3 100644
--- a/Projects/Sofia/Analysis.cxx
+++ b/Projects/Sofia/Analysis.cxx
@@ -90,7 +90,7 @@ void Analysis::FFAnalysis(){
for(int i=0; i<mult; i++){
double Time = ran.Gaus(InteractionCoordinates->GetTime(i) - BeamTimeOffset,0.02);
m_XE.push_back(ran.Gaus(InteractionCoordinates->GetDetectedPositionX(i),1));
- m_YE.push_back(ran.Gaus(InteractionCoordinates->GetDetectedPositionY(i),0));
+ m_YE.push_back(ran.Gaus(InteractionCoordinates->GetDetectedPositionY(i),1));
m_ZE.push_back(ran.Gaus(InteractionCoordinates->GetDetectedPositionZ(i),1));
TVector3 vE = TVector3(m_XE[i],m_YE[i],m_ZE[i]);
//TVector3 vE = TVector3(m_XE[i],0,m_ZE[i]);
@@ -99,10 +99,12 @@ void Analysis::FFAnalysis(){
int Adet = InteractionCoordinates->GetA(i);
//int Z = FissionConditions->GetFragmentZ(i);
int Z = InteractionCoordinates->GetZ(i);
- double ThetaIn = FissionConditions->GetFragmentTheta(i);
+ double dtheta = 0.2e-3*180./3.1415; // 0.2 mrad
+ double dphi = 10e-3*180./3.1415; // 10 mrad
+ double ThetaIn = ran.Gaus(FissionConditions->GetFragmentTheta(i),dtheta);
//if(FissionConditions->GetFragmentMomentumX(i)<0)
//ThetaIn = -ThetaIn;
- double Phi = FissionConditions->GetFragmentPhi(i);
+ double Phi = ran.Gaus(FissionConditions->GetFragmentPhi(i),dphi);
double Brho = FissionConditions->GetFragmentBrho(i);
TVector3 dir = TVector3(sin(ThetaIn*deg)*cos(Phi*deg), sin(ThetaIn*deg)*sin(Phi*deg), cos(ThetaIn*deg));
diff --git a/Projects/s455/Analysis.cxx b/Projects/s455/Analysis.cxx
index 3dd2bda17..c74327938 100644
--- a/Projects/s455/Analysis.cxx
+++ b/Projects/s455/Analysis.cxx
@@ -39,6 +39,7 @@ struct TofPair
double beta = -1;
double gamma = -1;
double theta_in = -10;
+ double phi_in = -10;
double theta_out = -10;
double psi = -10;
int plastic = -1;
@@ -118,17 +119,17 @@ void Analysis::Init(){
m_GladField = new GladFieldMap();
m_GladField->SetCurrent(fGladCurrent);
//m_GladField->SetGladEntrance(0, 0.02*m, 2.774*m + 0.5405*m);
- m_GladField->SetGladEntrance(0, 0, -1.1135*m);
- //m_GladField->SetGladTurningPoint(0, 0.02*m, 2.774*m + 0.5405*m + 1.1135*m);
+ m_GladField->SetGladEntrance(0, 0, -1113.5*mm);
m_GladField->SetGladTurningPoint(0, 0, 0);
m_GladField->SetGladTiltAngle(-14.*deg);
m_GladField->LoadMap("GladFieldMap_50mm.dat");
m_GladField->SetBin(50);
m_GladField->SetTimeStep(0.8);
m_GladField->SetCentralTheta(-20.*deg);
-
- double Z_MWPC3 = fDistanceGToMW3;// * cos(m_GladField->GetCentralTheta());
- double X_MWPC3 = (Z_MWPC3 - m_GladField->GetGladTurningPoint().Z())*tan(m_GladField->GetCentralTheta());
+
+ double R_MWPC3 = fDistanceGToMW3;
+ double Z_MWPC3 = R_MWPC3 * cos(m_GladField->GetCentralTheta());
+ double X_MWPC3 = R_MWPC3 * sin(m_GladField->GetCentralTheta());
m_GladField->Set_MWPC3_Position(X_MWPC3,0,Z_MWPC3);
}
@@ -620,8 +621,12 @@ void Analysis::FissionFragmentAnalysis(int which_cathode){
double ATToG = fDistanceStartToG - fDistancePlasticToCathode[which_cathode-1];
+ double ATToA = fDistanceStartToA - fDistancePlasticToCathode[which_cathode-1];
+ double GToTof = (fDistanceGToMW3 + fDistanceMW3ToToF) * cos(20*deg);
+ double Zdistance = ATToG + GToTof;
double Theta0 = 20.*deg;//m_GladField->GetCentralTheta();
double XA = 0;
+ double YA = 0;
double ZA = fDistanceStartToA - fDistanceStartToG;
double ZG = m_GladField->GetGladTurningPoint().Z();
double ZMW3 = m_GladField->Get_MWPC3_Position().Z();
@@ -640,10 +645,13 @@ void Analysis::FissionFragmentAnalysis(int which_cathode){
double MagB = m_GladField->GetB()/1000;
for(int i=0; i<2; i++){
XA = TofHit[i].DT;
- if(XA != -1e6 && SofBeamID->GetBeta()>0){
+ YA = (ATToA/Zdistance)*TofHit[i].y;
+ TofHit[i].phi_in = atan(TofHit[i].y/Zdistance);
+ if(XA != -1e6 && SofBeamID->GetBeta()>0){
+
TVector3 vG = TVector3(0,0,ZG);
- TVector3 vA = TVector3(XA,0,ZA);
+ TVector3 vA = TVector3(XA,YA,ZA);
// *** 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);
@@ -652,15 +660,15 @@ void Analysis::FissionFragmentAnalysis(int which_cathode){
TofHit[i].yc = (ZC/8592.)*TofHit[i].y;
TofHit[i].zc = ZC;
- int ix, iy;
+ /*int ix, iy;
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);
+ 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, TofHit[i].y, Z3lab);
- TVector3 dir = TVector3(sin(TofHit[i].theta_in), 0, cos(TofHit[i].theta_in));
+ TVector3 dir = TVector3(sin(TofHit[i].theta_in)*cos(TofHit[i].phi_in), sin(TofHit[i].theta_in)*sin(TofHit[i].phi_in), cos(TofHit[i].theta_in));
TofHit[i].x3lab = X3lab;
TofHit[i].z3lab = Z3lab;
TVector3 vOut = TVector3(X3lab-XC,0,Z3lab-ZC);
@@ -671,6 +679,7 @@ void Analysis::FissionFragmentAnalysis(int which_cathode){
//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));
XB = (XA+(ZI-ZA)*tan(TofHit[i].theta_in)) / (1-tan(Tilt)*tan(TofHit[i].theta_in));
@@ -757,8 +766,6 @@ void Analysis::FissionFragmentAnalysis(int which_cathode){
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();
TVector3 vBD = vD-vB;
double vBD_angle = vZ.Angle(vBD);
--
GitLab