From 4f6105f307583aabd42511bad67b2def500be570 Mon Sep 17 00:00:00 2001
From: Pierre Morfouace <pierre.morfouace2@cea.fr>
Date: Fri, 7 Jan 2022 16:13:18 +0100
Subject: [PATCH] * Updating Sofia project to match TofWall and MWPC3
---
NPLib/Detectors/Sofia/TSofFissionFragment.cxx | 9 +-
NPLib/Detectors/Sofia/TSofFissionFragment.h | 23 +-
NPLib/Detectors/Sofia/TSofMwpcPhysics.cxx | 197 ++++++++++++++----
NPLib/Detectors/Sofia/TSofMwpcPhysics.h | 22 +-
NPLib/Detectors/Sofia/TSofTofWPhysics.cxx | 12 +-
NPLib/Detectors/Sofia/TSofTofWPhysics.h | 5 +-
Projects/s455/Analysis.cxx | 145 ++++++++++---
.../s455/calibration/SofSci/ClockOffset.cal | 4 +-
.../SofTrim/SofTrim_SectionAlign.cal | 4 +-
Projects/s455/macro/GetChargeDist.C | 28 ++-
10 files changed, 352 insertions(+), 97 deletions(-)
diff --git a/NPLib/Detectors/Sofia/TSofFissionFragment.cxx b/NPLib/Detectors/Sofia/TSofFissionFragment.cxx
index 5f107d4e3..84620fb29 100644
--- a/NPLib/Detectors/Sofia/TSofFissionFragment.cxx
+++ b/NPLib/Detectors/Sofia/TSofFissionFragment.cxx
@@ -44,7 +44,7 @@ TSofFissionFragment::~TSofFissionFragment() {
//////////////////////////////////////////////////////////////////////
void TSofFissionFragment::Clear() {
fFF_Z.clear();
- fFF_Qmax.clear();
+ fFF_iZ.clear();
fFF_AoQ.clear();
fFF_A.clear();
fFF_Beta.clear();
@@ -52,7 +52,12 @@ void TSofFissionFragment::Clear() {
fFF_Gamma.clear();
fFF_Brho.clear();
fFF_DT.clear();
-
+ fFF_ThetaIn.clear();
+ fFF_TofPosX.clear();
+ fFF_TofPosY.clear();
+ fFF_MwpcPosX.clear();
+ fFF_MwpcPosY.clear();
+
fFF_Zsum = -1;
fFF_IntZsum = -1;
}
diff --git a/NPLib/Detectors/Sofia/TSofFissionFragment.h b/NPLib/Detectors/Sofia/TSofFissionFragment.h
index 0fbed2f6e..c277c833d 100644
--- a/NPLib/Detectors/Sofia/TSofFissionFragment.h
+++ b/NPLib/Detectors/Sofia/TSofFissionFragment.h
@@ -35,7 +35,7 @@ class TSofFissionFragment : public TObject {
// to allow multiplicity treatment
private:
vector<double> fFF_Z;
- vector<double> fFF_Qmax;
+ vector<int> fFF_iZ;
vector<double> fFF_AoQ;
vector<double> fFF_A;
vector<double> fFF_Beta;
@@ -43,6 +43,11 @@ class TSofFissionFragment : public TObject {
vector<double> fFF_Gamma;
vector<double> fFF_Brho;
vector<double> fFF_DT;
+ vector<double> fFF_ThetaIn;
+ vector<double> fFF_TofPosX;
+ vector<double> fFF_TofPosY;
+ vector<double> fFF_MwpcPosX;
+ vector<double> fFF_MwpcPosY;
double fFF_Zsum;
int fFF_IntZsum;
@@ -71,7 +76,7 @@ class TSofFissionFragment : public TObject {
inline void SetZsum(double val){fFF_Zsum = val;};//!
inline void SetIntZsum(int val){fFF_IntZsum = val;};//!
inline void SetZ(double val){fFF_Z.push_back(val);};//!
- inline void SetQmax(double val){fFF_Qmax.push_back(val);};//!
+ inline void SetiZ(int val){fFF_iZ.push_back(val);};//!
inline void SetAoQ(double val){fFF_AoQ.push_back(val);};//!
inline void SetA(double val){fFF_A.push_back(val);};//!
inline void SetBeta(double val){fFF_Beta.push_back(val);};//!
@@ -79,13 +84,19 @@ class TSofFissionFragment : public TObject {
inline void SetGamma(double val){fFF_Gamma.push_back(val);};//!
inline void SetBrho(double val){fFF_Brho.push_back(val);};//!
inline void SetDT(double val){fFF_DT.push_back(val);};//!
+ inline void SetThetaIn(double val){fFF_ThetaIn.push_back(val);};//!
+ inline void SetTofPosX(double val){fFF_TofPosX.push_back(val);};//!
+ inline void SetTofPosY(double val){fFF_TofPosY.push_back(val);};//!
+ inline void SetMwpcPosX(double val){fFF_MwpcPosX.push_back(val);};//!
+ inline void SetMwpcPosY(double val){fFF_MwpcPosY.push_back(val);};//!
+
////////////////////// GETTERS ////////////////////////
int GetMult() {return fFF_Z.size();}//!
inline double GetZsum() const {return fFF_Zsum;}//!
inline int GetIntZsum() const {return fFF_IntZsum;}//!
inline double GetZ(int i) const {return fFF_Z[i];}//!
- inline double GetQmax(int i) const {return fFF_Qmax[i];}//!
+ inline int GetiZ(int i) const {return fFF_iZ[i];}//!
inline double GetAoQ(int i) const {return fFF_AoQ[i];}//!
inline double GetA(int i) const {return fFF_A[i];}//!
inline double GetBeta(int i) const {return fFF_Beta[i];}//!
@@ -93,6 +104,12 @@ class TSofFissionFragment : public TObject {
inline double GetGamma(int i) const {return fFF_Gamma[i];}//!
inline double GetBrho(int i) const {return fFF_Brho[i];}//!
inline double GetDT(int i) const {return fFF_DT[i];}//!
+ inline double GetThetaIn(int i) const {return fFF_ThetaIn[i];}//!
+ inline double GetTofPosX(int i) const {return fFF_TofPosX[i];}//!
+ inline double GetTofPosY(int i) const {return fFF_TofPosY[i];}//!
+ inline double GetMwpcPosX(int i) const {return fFF_MwpcPosX[i];}//!
+ inline double GetMwpcPosY(int i) const {return fFF_MwpcPosY[i];}//!
+
//////////////////////////////////////////////////////////////
// Required for ROOT dictionnary
diff --git a/NPLib/Detectors/Sofia/TSofMwpcPhysics.cxx b/NPLib/Detectors/Sofia/TSofMwpcPhysics.cxx
index 182332b55..a7351c1f4 100644
--- a/NPLib/Detectors/Sofia/TSofMwpcPhysics.cxx
+++ b/NPLib/Detectors/Sofia/TSofMwpcPhysics.cxx
@@ -41,6 +41,7 @@ using namespace std;
ClassImp(TSofMwpcPhysics)
+ bool compare_charge(pair<int,int> p1, pair<int,int> p2) {return p1.second < p2.second;}
///////////////////////////////////////////////////////////////////////////
TSofMwpcPhysics::TSofMwpcPhysics()
@@ -80,7 +81,6 @@ void TSofMwpcPhysics::BuildSimplePhysicalEvent() {
///////////////////////////////////////////////////////////////////////////
void TSofMwpcPhysics::BuildPhysicalEvent() {
-
PreTreat();
unsigned int mysizeE = m_PreTreatedData->GetMultiplicity();
@@ -95,6 +95,8 @@ void TSofMwpcPhysics::BuildPhysicalEvent() {
double ChargeMaxX2[4]={0,0,0,0};
double ChargeMaxY[4] ={0,0,0,0};
+ fPairX.reserve(50);
+ fPairY.reserve(50);
for(int i=0; i<m_NumberOfDetectors; i++){
Buffer_X1_Q.push_back(ChargeArray);
@@ -103,11 +105,10 @@ void TSofMwpcPhysics::BuildPhysicalEvent() {
}
for (UShort_t e = 0; e < mysizeE ; e++) {
- int det = m_PreTreatedData->GetDetectorNbr(e);
+ int det = m_PreTreatedData->GetDetectorNbr(e);//start at 1
int plane = m_PreTreatedData->GetPlane(e);
int strip = m_PreTreatedData->GetPad(e);// starts at 0
double charge = m_PreTreatedData->GetCharge(e);
-
// Xup for MWPC2 and 3 and X for MWPC1 and 4
if(plane==1){
Buffer_X1_Q[det-1][strip] = charge;
@@ -116,9 +117,14 @@ void TSofMwpcPhysics::BuildPhysicalEvent() {
ChargeMaxX1[det-1] = charge;
StripMaxX1[det-1] = strip;
}
+ if(det==4){
+ pair<int, int> hit_pair = make_pair(strip,charge);
+ fPairX.push_back(hit_pair);
+ }
}
+
// Xdown for MWPC2 and 3 and nothing for MWPC1 and 4
- else if(plane==2){
+ if(plane==2){
Buffer_X2_Q[det-1][strip] = charge;
if(charge>ChargeMaxX2[det-1]){
@@ -126,6 +132,7 @@ void TSofMwpcPhysics::BuildPhysicalEvent() {
StripMaxX2[det-1] = strip;
}
}
+
// Y for all MWPCx
if(plane==3){
Buffer_Y_Q[det-1][strip] = charge;
@@ -134,62 +141,167 @@ void TSofMwpcPhysics::BuildPhysicalEvent() {
ChargeMaxY[det-1] = charge;
StripMaxY[det-1] = strip;
}
+ if(det==4){
+ pair<int, int> hit_pair = make_pair(strip,charge);
+ fPairY.push_back(hit_pair);
+ }
}
}
- double X1 = -100;
- double X2 = -100;
- double Y = -100;
+ double X1 = -1000;
+ double X2 = -1000;
+ double Y = -1000;
for(int i=0; i<m_NumberOfDetectors; i++){
- double qleft = Buffer_X1_Q[i][StripMaxX1[i]-1];
- double qright = Buffer_X1_Q[i][StripMaxX1[i]+1];
- if(qleft>0 && qright>0){
- X1 = GetPositionX(ChargeMaxX1[i], StripMaxX1[i], qleft, qright);
- X1 = X1 - DetPosX[i];
- }
- qleft = Buffer_X2_Q[i][StripMaxX2[i]-1];
- qright = Buffer_X2_Q[i][StripMaxX2[i]+1];
- if(qleft>0 && qright>0){
- X2 = GetPositionX(ChargeMaxX2[i], StripMaxX2[i], qleft, qright);
- X2 = X2 - DetPosX[i];
- }
- double qdown = Buffer_Y_Q[i][StripMaxY[i]-1];
- double qup = Buffer_Y_Q[i][StripMaxY[i]+1];
- if(qdown>0 && qup>0){
- Y = GetPositionY(ChargeMaxY[i], StripMaxY[i], qdown, qup);
- Y = Y - DetPosY[i];
+ int det_num = i+1;
+ if(det_num<4){
+ double qleft = Buffer_X1_Q[i][StripMaxX1[i]-1];
+ double qright = Buffer_X1_Q[i][StripMaxX1[i]+1];
+ if(qleft>0 && qright>0){
+ X1 = GetPositionX(det_num, ChargeMaxX1[i], StripMaxX1[i], qleft, qright);
+ X1 = X1 - DetPosX[i];
+ }
+ qleft = Buffer_X2_Q[i][StripMaxX2[i]-1];
+ qright = Buffer_X2_Q[i][StripMaxX2[i]+1];
+ if(qleft>0 && qright>0){
+ X2 = GetPositionX(det_num, ChargeMaxX2[i], StripMaxX2[i], qleft, qright);
+ X2 = X2 - DetPosX[i];
+ }
+ double qdown = Buffer_Y_Q[i][StripMaxY[i]-1];
+ double qup = Buffer_Y_Q[i][StripMaxY[i]+1];
+ if(qdown>0 && qup>0){
+ Y = GetPositionY(det_num, ChargeMaxY[i], StripMaxY[i], qdown, qup);
+ Y = Y - DetPosY[i];
+ }
+
+ DetectorNbr.push_back(det_num);
+ PositionX1.push_back(X1);
+ PositionX2.push_back(X2);
+ PositionY.push_back(Y);
}
- DetectorNbr.push_back(i+1);
- PositionX1.push_back(X1);
- PositionX2.push_back(X2);
- PositionY.push_back(Y);
+ if(det_num==4 && fPairX.size()>0 && fPairY.size()>0){
+ sort(fPairX.begin(), fPairX.end());
+ sort(fPairY.begin(), fPairY.end());
+
+ fThresholdX = max(0.2 * GetChargeMax(fPairX), 500.);
+ fThresholdY = max(0.2 * GetChargeMax(fPairY), 500.);
+
+ // *** Finding X clusters *** //
+ while(GetChargeMax(fPairX) > fThresholdX && fPairX.size()>4){
+ const auto pelt = max_element(fPairX.begin(), fPairX.end(), compare_charge);
+ int indice = distance(fPairX.begin(), pelt);
+
+ if(indice>0 && indice<fPairX.size()-1){
+ fClusterX.push_back(FindCluster(fPairX));
+ }
+ else break;
+ }
+
+ // *** Finding Y clusters *** //
+ while(GetChargeMax(fPairY) > fThresholdY && fPairY.size() > 4){
+ const auto pelt = max_element(fPairY.begin(), fPairY.end(), compare_charge);
+ int indice = distance(fPairY.begin(), pelt);
+
+ if(indice>0 && indice<fPairY.size()-1){
+ fClusterY.push_back(FindCluster(fPairY));
+ }
+ else break;
+ }
+
+
+ if(fClusterX.size()==fClusterY.size()){
+ int size = fClusterX.size();
+ vector<double> Xpos;
+ vector<double> Ypos;
+ for(unsigned int i=0; i<size; i++){
+ // *** strip X *** //
+ vector<pair<int,int>> hitX;
+ hitX = fClusterX[i];
+ double x = -1000;
+ int qleft = hitX[0].second;
+ int qmax = hitX[1].second;
+ int qright = hitX[2].second;
+ int padmax = hitX[2].first;
+ if(padmax>0 && padmax+1<288 && qmax>0 && qleft>0 && qright>0){
+ x = GetPositionX(det_num, qmax, padmax, qleft, qright);
+ Xpos.push_back(x);
+ }
+
+ // *** strip Y *** //
+ vector<pair<int,int>> hitY;
+ hitY = fClusterY[i];
+ double y = -1000;
+ int qdown = hitY[0].second;
+ qmax = hitY[1].second;
+ int qup = hitY[2].second;
+ padmax = hitY[2].first;
+ if(padmax>0 && padmax+1<120 && qmax>0 && qdown>0 && qup>0){
+ y = GetPositionY(det_num, qmax, padmax, qdown, qup);
+ Ypos.push_back(y);
+ }
+ DetectorNbr.push_back(det_num);
+ PositionX1.push_back(x);
+ PositionX2.push_back(-1000);
+ PositionY.push_back(y);
+ }
+ }
+ }
}
}
///////////////////////////////////////////////////////////////////////////
-double TSofMwpcPhysics::GetPositionX(double qmax, int padmax, double qleft, double qright){
- double fwx = 3.125;
- double fSize = 200.0;
+vector<pair<int,int>> TSofMwpcPhysics::FindCluster(vector<pair<int,int>>& p1){
+ vector<pair<int,int>> hit;
+ hit.clear();
+ pair<int,int> couple;
+
+ const auto pelt = max_element(p1.begin(), p1.end(), compare_charge);
+ int padmax = pelt->first;
+ int qmax = pelt->second;
+ int imax = distance(p1.begin(), pelt);
+
+ couple = {p1[imax - 1].first, p1[imax - 1].second};
+ //cout << "pad= " << p1[imax-1].first << " / Q= " << p1[imax-1].second << endl;
+ hit.push_back(couple);
+ couple = {p1[imax].first, p1[imax].second};
+ //cout << "pad= " << p1[imax].first << " / Q= " << p1[imax].second << endl;
+ hit.push_back(couple);
+ couple = {p1[imax + 1].first, p1[imax + 1].second};
+ //cout << "pad= " << p1[imax+1].first << " / Q= " << p1[imax+1].second << endl;
+ hit.push_back(couple);
+ // Removing those three points from fPair
+ p1.erase(p1.begin() + imax -1, p1.begin() + imax);
+ p1.erase(p1.begin() + imax -1, p1.begin() + imax);
+ p1.erase(p1.begin() + imax -1, p1.begin() + imax);
+
+
+ return hit;
+}
+///////////////////////////////////////////////////////////////////////////
+double TSofMwpcPhysics::GetChargeMax(vector<pair<int,int>> p1){
+ const auto comp = max_element(p1.begin(), p1.end(), compare_charge);
+ double Qmax = comp->second;
+ return Qmax;
+}
- double a3 = TMath::Pi() * fwx / (TMath::ACosH(0.5 * (TMath::Sqrt(qmax / qleft) + TMath::Sqrt(qmax / qright))));
+///////////////////////////////////////////////////////////////////////////
+double TSofMwpcPhysics::GetPositionX(int det, double qmax, int padmax, double qleft, double qright){
+ double a3 = TMath::Pi() * fwx[det-1] / (TMath::ACosH(0.5 * (TMath::Sqrt(qmax / qleft) + TMath::Sqrt(qmax / qright))));
Double_t a2 = (a3 / TMath::Pi()) * TMath::ATanH((TMath::Sqrt(qmax / qleft) - TMath::Sqrt(qmax / qright)) /
- (2 * TMath::SinH(TMath::Pi() * fwx / a3)));
+ (2 * TMath::SinH(TMath::Pi() * fwx[det-1] / a3)));
- return (-1. * padmax * fwx + (fSize / 2) - (fwx / 2) - a2); // Left is positive and right negative
+ return (-1. * padmax * fwx[det-1] + (fSizeX[det-1] / 2) - (fwx[det-1] / 2) - a2); // Left is positive and right negative
}
///////////////////////////////////////////////////////////////////////////
-double TSofMwpcPhysics::GetPositionY(double qmax, int padmax, double qdown, double qup){
- double fwy = 3.125;
- double fSize = 200.0;
- double a3 = TMath::Pi() * fwy / (TMath::ACosH(0.5 * (TMath::Sqrt(qmax / qdown) + TMath::Sqrt(qmax / qup))));
+double TSofMwpcPhysics::GetPositionY(int det, double qmax, int padmax, double qdown, double qup){
+ double a3 = TMath::Pi() * fwy[det-1] / (TMath::ACosH(0.5 * (TMath::Sqrt(qmax / qdown) + TMath::Sqrt(qmax / qup))));
Double_t a2 = (a3 / TMath::Pi()) * TMath::ATanH((TMath::Sqrt(qmax / qdown) - TMath::Sqrt(qmax / qup)) /
- (2 * TMath::SinH(TMath::Pi() * fwy / a3)));
+ (2 * TMath::SinH(TMath::Pi() * fwy[det-1] / a3)));
- return (padmax * fwy - (fSize / 2) + (fwy / 2) + a2);
+ return (padmax * fwy[det-1] - (fSizeY[det-1] / 2) + (fwy[det-1] / 2) + a2);
}
///////////////////////////////////////////////////////////////////////////
@@ -281,6 +393,11 @@ void TSofMwpcPhysics::Clear() {
Buffer_X1_Q.clear();
Buffer_X2_Q.clear();
Buffer_Y_Q.clear();
+
+ fPairX.clear();
+ fPairY.clear();
+ fClusterX.clear();
+ fClusterY.clear();
}
diff --git a/NPLib/Detectors/Sofia/TSofMwpcPhysics.h b/NPLib/Detectors/Sofia/TSofMwpcPhysics.h
index 14935bd46..ee4b6a18c 100644
--- a/NPLib/Detectors/Sofia/TSofMwpcPhysics.h
+++ b/NPLib/Detectors/Sofia/TSofMwpcPhysics.h
@@ -131,9 +131,25 @@ class TSofMwpcPhysics : public TObject, public NPL::VDetector {
// needed for online analysis for example
void SetRawDataPointer(TSofMwpcData* rawDataPointer) {m_EventData = rawDataPointer;}
- double GetPositionX(double qmax, int padmax, double qleft, double qright);
- double GetPositionY(double qmax, int padmax, double qdown, double qup);
-
+ double fwx[4]={3.125, 3.125, 3.125, 3.125};
+ double fwy[4]={3.125, 5.000, 5.000, 5.000};
+ double fSizeX[4]={200, 200, 200, 900};
+ double fSizeY[4]={200, 200, 200, 600};
+
+ vector<pair<int, int>> fPairX;
+ vector<pair<int, int>> fPairY;
+ vector<vector<pair<int,int>>> fClusterX;
+ vector<vector<pair<int,int>>> fClusterY;
+
+ double fThresholdX;
+ double fThresholdY;
+
+ double GetChargeMax(vector<pair<int,int>> p1);
+ vector<pair<int,int>> FindCluster(vector<pair<int,int>>& p1);
+
+ double GetPositionX(int det, double qmax, int padmax, double qleft, double qright);
+ double GetPositionY(int det, double qmax, int padmax, double qdown, double qup);
+
// objects are not written in the TTree
private:
TSofMwpcData* m_EventData; //!
diff --git a/NPLib/Detectors/Sofia/TSofTofWPhysics.cxx b/NPLib/Detectors/Sofia/TSofTofWPhysics.cxx
index 2eb7514bc..99e1bb88f 100644
--- a/NPLib/Detectors/Sofia/TSofTofWPhysics.cxx
+++ b/NPLib/Detectors/Sofia/TSofTofWPhysics.cxx
@@ -120,16 +120,18 @@ void TSofTofWPhysics::BuildPhysicalEvent() {
static CalibrationManager* Cal = CalibrationManager::getInstance();
for(int p=0; p<m_NumberOfPlastics; p++){
if(mult1[p]==1 && mult2[p]==1){
+ double calposX = -(p*PlasticWidth + rand.Uniform(0,PlasticWidth) - (double)m_NumberOfPlastics/2*PlasticWidth);
+
double time_ns = 0.5*(T1[p][0] + T2[p][0]);
- double rawpos = T1[p][0] - T2[p][0];
- double calpos = Cal->ApplyCalibration("SofTofW/TOFW"+NPL::itoa(p+1)+"_POSPAR",rawpos);
+ double rawposY = T1[p][0] - T2[p][0];
+ double calposY = Cal->ApplyCalibration("SofTofW/TOFW"+NPL::itoa(p+1)+"_POSPAR",rawposY);
double rawtof = time_ns - m_StartTime;
double caltof = Cal->ApplyCalibration("SofTofW/TOFW"+NPL::itoa(p+1)+"_TOFPAR",rawtof) + m_TofAlignedValue;
PlasticNbr.push_back(p+1);
TimeNs.push_back(time_ns);
- RawPosY.push_back(rawpos);
- CalPosY.push_back(calpos);
+ CalPosX.push_back(calposX);
+ CalPosY.push_back(calposY);
RawTof.push_back(rawtof);
CalTof.push_back(caltof);
}
@@ -256,7 +258,7 @@ void TSofTofWPhysics::ReadAnalysisConfig() {
void TSofTofWPhysics::Clear() {
PlasticNbr.clear();
TimeNs.clear();
- RawPosY.clear();
+ CalPosX.clear();
CalPosY.clear();
RawTof.clear();
CalTof.clear();
diff --git a/NPLib/Detectors/Sofia/TSofTofWPhysics.h b/NPLib/Detectors/Sofia/TSofTofWPhysics.h
index 94a402ec7..75ce72a12 100644
--- a/NPLib/Detectors/Sofia/TSofTofWPhysics.h
+++ b/NPLib/Detectors/Sofia/TSofTofWPhysics.h
@@ -39,7 +39,8 @@ using namespace std;
#include "NPVDetector.h"
#include "NPInputParser.h"
-
+#define PlasticWidth 32
+#define PlasticHeight 660
class TSofTofWPhysics : public TObject, public NPL::VDetector {
//////////////////////////////////////////////////////////////
@@ -62,7 +63,7 @@ class TSofTofWPhysics : public TObject, public NPL::VDetector {
public:
vector<int> PlasticNbr;
vector<double> TimeNs;
- vector<double> RawPosY;
+ vector<double> CalPosX;
vector<double> CalPosY;
vector<double> RawTof;
vector<double> CalTof;
diff --git a/Projects/s455/Analysis.cxx b/Projects/s455/Analysis.cxx
index 76eb904db..6fe048409 100644
--- a/Projects/s455/Analysis.cxx
+++ b/Projects/s455/Analysis.cxx
@@ -77,7 +77,6 @@ void Analysis::FissionFragmentAnalysis(){
double TOF_CC[2];
double Plastic[2];
- double PosY[2];
double Plastic_left = -1;
double Plastic_right = -1;
double TOF_left = -1;
@@ -102,18 +101,76 @@ void Analysis::FissionFragmentAnalysis(){
double Beta3 = -1;
double Beta4 = -1;
double Beta_norm = 0.745;
+ double Gamma1 = -1;
+ double Gamma2 = -1;
+ double Beta_Z1 = -1;
+ double Beta_Z2 = -1;
+ double Brho1 = -1;
+ double Brho2 = -1;
+ double AoQ1 = -1;
+ double AoQ2 = -1;
+ double A1 = -1;
+ double A2 = -1;
+ double Z1 = -1;
+ double Z2 = -1;
+ double Zsum = -1;
+ int iZ1 = -1;
+ int iZ2 = -1;
+ int iZsum = -1;
+ double ThetaIn1 = -1000;
+ double ThetaIn2 = -1000;
for(int i = 0; i<2; i++){
TOF_CC[i] = -1;
Plastic[i] = -1;
- PosY[i] = -1;
}
+ vector<double> PosY;
+ vector<double> PosX;
if(softofw_size==2){
for(unsigned int i=0; i<softofw_size; i++){
TOF_CC[i] = SofTofW->CalTof[i];
Plastic[i] = SofTofW->PlasticNbr[i];
- PosY[i] = SofTofW->CalPosY[i];
+ PosX.push_back(SofTofW->CalPosX[i]);
+ PosY.push_back(SofTofW->CalPosY[i]);
+
+ //SofFF->SetTofPosX(SofTofW->CalPosX[i]);
+ //SofFF->SetTofPosY(SofTofW->CalPosY[i]);
+ }
+ }
+
+ vector<double> Xmwpc4;
+ vector<double> Ymwpc4;
+ for(unsigned int i=0; i<SofMwpc->DetectorNbr.size(); i++){
+ if(SofMwpc->DetectorNbr[i]==4){
+ Xmwpc4.push_back(SofMwpc->PositionX1[i]);
+ Ymwpc4.push_back(SofMwpc->PositionY[i]);
+
+ //SofFF->SetMwpcPosX(SofMwpc->PositionX1[i]);
+ //SofFF->SetMwpcPosY(SofMwpc->PositionY[i]);
+ }
+ }
+
+ vector<double> good_posx;
+ vector<double> good_posy;
+ for(unsigned int i=0; i<PosX.size(); i++){
+ double tofx = PosX[i];
+ double tofy = PosY[i];
+ for(unsigned int k=0; k<Xmwpc4.size(); k++){
+ double posx = Xmwpc4[k];
+ if(abs(posx-tofx) < 100){
+ good_posx.push_back(posx);
+ SofFF->SetMwpcPosX(posx);
+ SofFF->SetTofPosX(tofx);
+ }
+ }
+ for(unsigned int p=0; p<Ymwpc4.size(); p++){
+ double posy = Ymwpc4[p];
+ if(abs(posy-tofy) < 20){
+ good_posy.push_back(posy);
+ SofFF->SetMwpcPosY(posy);
+ SofFF->SetTofPosY(tofy);
+ }
}
}
@@ -181,13 +238,15 @@ void Analysis::FissionFragmentAnalysis(){
TOF_right = TOF_CC[0];
}
- if(PosY[0]>PosY[1]){
- TOF_up = TOF_CC[0];
- TOF_down = TOF_CC[1];
- }
- if(PosY[0]<PosY[1]){
- TOF_up = TOF_CC[1];
- TOF_down = TOF_CC[0];
+ if(PosY.size()==2){
+ if(PosY[0]>PosY[1]){
+ TOF_up = TOF_CC[0];
+ TOF_down = TOF_CC[1];
+ }
+ if(PosY[0]<PosY[1]){
+ TOF_up = TOF_CC[1];
+ TOF_down = TOF_CC[0];
+ }
}
double velocity_left = L_CC/TOF_left;
@@ -245,34 +304,41 @@ void Analysis::FissionFragmentAnalysis(){
// Z calibration //
- double Z1=-1;
- double Z2=-1;
- double Zsum=-1;
- int iZsum=-1;
-
if(E1>0 && E2>0 && E3==-1 && E4==-1){
Z1 = E1;
Z2 = E2;
+ Beta_Z1 = Beta1;
+ Beta_Z2 = Beta2;
}
if(E1>0 && E2==-1 && E3>0 && E4==-1){
Z1 = E1;
Z2 = E3;
+ Beta_Z1 = Beta1;
+ Beta_Z2 = Beta3;
}
if(E1>0 && E2==-1 && E3==-1 && E4>0){
Z1 = E1;
Z2 = E4;
+ Beta_Z1 = Beta1;
+ Beta_Z2 = Beta4;
}
if(E1==-1 && E2>0 && E3>0 && E4==-1){
Z1 = E2;
Z2 = E3;
+ Beta_Z1 = Beta2;
+ Beta_Z2 = Beta3;
}
if(E1==-1 && E2>0 && E3==-1 && E4>0){
Z1 = E2;
Z2 = E4;
+ Beta_Z1 = Beta2;
+ Beta_Z2 = Beta4;
}
if(E1==-1 && E2==-1 && E3>0 && E4>0){
Z1 = E3;
Z2 = E4;
+ Beta_Z1 = Beta3;
+ Beta_Z2 = Beta4;
}
if(Z1>0 && Z2>0){
@@ -284,20 +350,37 @@ void Analysis::FissionFragmentAnalysis(){
Zsum = Z1+Z2;
- int iZ1 = (int) round(Z1);
- int iZ2 = (int) round(Z2);
+ iZ1 = (int) round(Z1);
+ iZ2 = (int) round(Z2);
iZsum = iZ1 + iZ2;
}
+ Gamma1 = 1. / sqrt(1 - Beta_Z1 * Beta_Z1);
+ Gamma2 = 1. / sqrt(1 - Beta_Z2 * Beta_Z2);
+
+ AoQ1 = Brho1 / (3.10761 * Beta_Z1 * Gamma1);
+ AoQ2 = Brho2 / (3.10761 * Beta_Z2 * Gamma2);
+
+ A1 = AoQ1 * iZ1;
+ A2 = AoQ2 * iZ2;
+
+ //*** Filling the Fission Fragment Tree ***//
SofFF->SetTOF(TOF_left);
SofFF->SetTOF(TOF_right);
- SofFF->SetBeta(Beta1);
- SofFF->SetBeta(Beta2);
- SofFF->SetBeta(Beta3);
- SofFF->SetBeta(Beta4);
-
+ SofFF->SetBeta(Beta_Z1);
+ SofFF->SetBeta(Beta_Z2);
+ SofFF->SetGamma(Gamma1);
+ SofFF->SetGamma(Gamma2);
+ SofFF->SetiZ(iZ1);
+ SofFF->SetiZ(iZ2);
SofFF->SetZ(Z1);
SofFF->SetZ(Z2);
+ SofFF->SetAoQ(AoQ1);
+ SofFF->SetAoQ(AoQ1);
+ SofFF->SetA(A1);
+ SofFF->SetA(A2);
+ SofFF->SetBrho(Brho1);
+ SofFF->SetBrho(Brho2);
SofFF->SetDT(DT1);
SofFF->SetDT(DT2);
@@ -325,12 +408,10 @@ void Analysis::BeamAnalysis(){
double Anode3 = SofTrim->EnergySection[2];
if(fRunID==13)
Zbeam = max(Anode1, Anode2);
- else if(fRunID==14)
- Zbeam = max(Anode2, Anode3);
else
Zbeam = SofTrim->GetMaxEnergySection();
- Qmax = DetermineQmax();
+ //Qmax = DetermineQmax();
Theta = SofTrim->Theta[0];
double TofFromS2 = SofSci->CalTof[0];
@@ -371,7 +452,7 @@ void Analysis::BeamAnalysis(){
// Filling Beam tree
SofBeamID->SetZbeam(Zbeam);
- SofBeamID->SetQmax(rand.Gaus(Qmax,0.15));
+ //SofBeamID->SetQmax(rand.Gaus(Qmax,0.15));
SofBeamID->SetAoQ(AoQ);
SofBeamID->SetAbeam(A);
SofBeamID->SetBeta(Beta);
@@ -493,7 +574,7 @@ void Analysis::InitParameter(){
fDCC = -10000;
fK_LS2 = -30e-8;
- fRunID = 13;
+ fRunID = 6;
// Beam parameter //
fZBeta_p0 = 1;
@@ -503,14 +584,14 @@ void Analysis::InitParameter(){
fZff_p0 = 2.80063;
fZff_p1 = 6.91985e-2;
fZff_p2 = 1.01598e-7;
-
+
if(fRunID==1 || fRunID==2){
//fBrho0 = 10.6813; // 180Hg
fBrho0 = 10.6955; // 180Hg
fZbeam_p0 = -5303.06;
fZbeam_p1 = 0.674945;
fZbeam_p2 = -8.32085e-6;
-
+
fZBeta_p0 = 72.946;
fZBeta_p1 = 6.0644;
}
@@ -543,7 +624,7 @@ void Analysis::InitParameter(){
fZbeam_p0 = 1590.66;
fZbeam_p1 = 0.0956455;
fZbeam_p2 = 3.84585e-6;
-
+
fZBeta_p0 = 74.6063;
fZBeta_p1 = 6.4635;
}
@@ -561,7 +642,7 @@ void Analysis::InitParameter(){
fZbeam_p0 = 4122.94;
fZbeam_p1 = -0.119867;
fZbeam_p2 = 8.29115e-6;
-
+
fZBeta_p0 = 63.9575;
fZBeta_p1 = 25.1988;
}
@@ -588,7 +669,7 @@ void Analysis::InitParameter(){
fZbeam_p0 = -116.425;
fZbeam_p1 = 0.218256;
fZbeam_p2 = 1.62399e-6;
-
+
fZBeta_p0 = 61.3889;
fZBeta_p1 = 25.8908;
}
diff --git a/Projects/s455/calibration/SofSci/ClockOffset.cal b/Projects/s455/calibration/SofSci/ClockOffset.cal
index ec8ed402e..c29a83575 100644
--- a/Projects/s455/calibration/SofSci/ClockOffset.cal
+++ b/Projects/s455/calibration/SofSci/ClockOffset.cal
@@ -1,6 +1,6 @@
SofSci_DET1_SIGNAL1_CLOCKOFFSET 0
SofSci_DET1_SIGNAL2_CLOCKOFFSET 0
SofSci_DET1_SIGNAL3_CLOCKOFFSET 0
-SofSci_DET2_SIGNAL1_CLOCKOFFSET 0
-SofSci_DET2_SIGNAL2_CLOCKOFFSET 0
+SofSci_DET2_SIGNAL1_CLOCKOFFSET 1
+SofSci_DET2_SIGNAL2_CLOCKOFFSET 1
SofSci_DET2_SIGNAL3_CLOCKOFFSET 0
diff --git a/Projects/s455/calibration/SofTrim/SofTrim_SectionAlign.cal b/Projects/s455/calibration/SofTrim/SofTrim_SectionAlign.cal
index df3a40198..8a3e3df45 100644
--- a/Projects/s455/calibration/SofTrim/SofTrim_SectionAlign.cal
+++ b/Projects/s455/calibration/SofTrim/SofTrim_SectionAlign.cal
@@ -1,3 +1,3 @@
-SofTrim_SEC1_ALIGN 912.97 0.987851
+SofTrim_SEC1_ALIGN 310.216 1.00792
SofTrim_SEC2_ALIGN 0 1
-SofTrim_SEC3_ALIGN -2002.22 1.15275
+SofTrim_SEC3_ALIGN -1382.53 1.13744
diff --git a/Projects/s455/macro/GetChargeDist.C b/Projects/s455/macro/GetChargeDist.C
index 37588e0df..0806d1275 100644
--- a/Projects/s455/macro/GetChargeDist.C
+++ b/Projects/s455/macro/GetChargeDist.C
@@ -90,7 +90,7 @@ void GetChargeDist(string nucleus_name){
cut_beam->SetVarY("fBeam_Z");
double Rx=0.005;
double Ry=0.5;
- AoQ = AoQ;// + 0.002;
+ AoQ = AoQ + 0.003;
double Zcenter = Z;
for(int i=0; i<360; i++){
double x = AoQ+Rx*cos(i*TMath::Pi()/180);
@@ -202,8 +202,13 @@ void GetChargeDist(string nucleus_name){
void LoadRootFiles(string iso){
chain = new TChain("PhysicsTree");
+ //*** Pa isotopes ***//
+ if(iso=="221Pa" || iso=="222Pa"){
+ chain->Add("../../../Outputs/Analysis/221Pa.root");
+ }
+
//*** Th isotopes ***//
- if(iso=="217Th" || iso=="216Th" || iso=="209Ra" || iso=="210Ra" || iso=="211Ra" || iso=="212Ra" || iso=="213Ra" || iso=="214Ra"){
+ if(iso=="217Th" || iso=="216Th"){
chain->Add("../../../Outputs/Analysis/216Th_1.root");
}
if(iso=="219Th" || iso=="220Th"){
@@ -211,12 +216,15 @@ void LoadRootFiles(string iso){
}
//*** Ac isotopes ***//
- if(iso=="212Ac" || iso=="213Ac" || iso=="214Ac" || iso=="215Ac" || iso=="216Ac"){
+ if(iso=="212Ac" || iso=="213Ac" || iso=="214Ac" || iso=="215Ac" || iso=="216Ac"){
chain->Add("../../../Outputs/Analysis/216Th_1.root");
}
-
+ if(iso=="214Ac" || iso=="215Ac" || iso=="216Ac" || iso=="217Ac" || iso=="218Ac"){
+ chain->Add("../../../Outputs/Analysis/221Pa.root");
+ }
+
//*** Ra isotopes ***//
- if(iso=="209Ra" || iso=="210Ra" || iso=="211Ra" || iso=="212Ra" || iso=="213Ra" || iso=="214Ra"){
+ if(iso=="209Ra" || iso=="210Ra" || iso=="211Ra" || iso=="212Ra" || iso=="213Ra" || iso=="214Ra"){
chain->Add("../../../Outputs/Analysis/216Th_1.root");
}
@@ -354,12 +362,20 @@ void LoadRootFiles(string iso){
}
//*** Ir isotopes ***//
- if(iso=="172Ir" || iso=="173Ir" || iso=="174Ir" || iso=="175Ir" || iso=="176It"){
+ if(iso=="171Ir" ||iso=="172Ir" || iso=="173Ir" || iso=="174Ir" || iso=="175Ir" || iso=="176Ir"){
chain->Add("../../../Outputs/Analysis/175Pt.root");
chain->Add("../../../Outputs/Analysis/180Hg_1.root");
chain->Add("../../../Outputs/Analysis/180Hg_2.root");
}
+ //*** Os isotopes ***//
+ if(iso=="169Os" ||iso=="170Os" || iso=="171Os" || iso=="172Os" || iso=="173Os" || iso=="174Os"){
+ chain->Add("../../../Outputs/Analysis/175Pt.root");
+ chain->Add("../../../Outputs/Analysis/180Hg_1.root");
+ chain->Add("../../../Outputs/Analysis/180Hg_2.root");
+ }
+
+
cout << "Number of entries: " << chain->GetEntries() << endl;
}
--
GitLab