From cf0dcf9c762d3b6ecbebbed4871ccda460801786 Mon Sep 17 00:00:00 2001
From: Pierre Morfouace <pierre.morfouace2@cea.fr>
Date: Fri, 11 Jun 2021 16:20:09 +0200
Subject: [PATCH] Updating SofTofW Physics analysis

---
 NPLib/Detectors/Sofia/TSofSciPhysics.cxx  |   2 +-
 NPLib/Detectors/Sofia/TSofTofWPhysics.cxx | 118 ++++++++++++++++++----
 NPLib/Detectors/Sofia/TSofTofWPhysics.h   |  21 +++-
 Projects/s455/Analysis.cxx                |  16 ++-
 Projects/s455/calibration.txt             |   2 +
 Projects/s455/s455.detector               |   9 +-
 6 files changed, 135 insertions(+), 33 deletions(-)

diff --git a/NPLib/Detectors/Sofia/TSofSciPhysics.cxx b/NPLib/Detectors/Sofia/TSofSciPhysics.cxx
index 0d869d13c..da7b7ad1c 100644
--- a/NPLib/Detectors/Sofia/TSofSciPhysics.cxx
+++ b/NPLib/Detectors/Sofia/TSofSciPhysics.cxx
@@ -235,7 +235,7 @@ double TSofSciPhysics::CalculateTimeNs(int det, int pmt, int ft, int ct){
   double par = Cal->GetValue("SofSci/DET"+NPL::itoa(det)+"_SIGNAL"+NPL::itoa(pmt)+"_TIME",ft);
   double r = (double)rand.Rndm()-0.5;
   double ClockOffset = Cal->GetValue("SofSci/DET"+NPL::itoa(det)+"_SIGNAL"+NPL::itoa(pmt)+"_CLOCKOFFSET",0);
-  double ict_ns = ((double)ct - ClockOffset) * 5.; // to do... take care of the clock offset
+  double ict_ns = ((double)ct - ClockOffset) * 5.; 
   double ift_ns;
 
   if(r<0){
diff --git a/NPLib/Detectors/Sofia/TSofTofWPhysics.cxx b/NPLib/Detectors/Sofia/TSofTofWPhysics.cxx
index f4d85066f..8e98fb7c5 100644
--- a/NPLib/Detectors/Sofia/TSofTofWPhysics.cxx
+++ b/NPLib/Detectors/Sofia/TSofTofWPhysics.cxx
@@ -49,6 +49,9 @@ TSofTofWPhysics::TSofTofWPhysics()
      m_EventPhysics(this),
      m_E_RAW_Threshold(0), // adc channels
      m_E_Threshold(0),     // MeV
+     m_NumberOfPlastics(28),
+     m_StartTime(-1),
+     m_TofAlignedValue(0), // ns
      m_NumberOfDetectors(0) {
 }
 
@@ -79,14 +82,88 @@ void TSofTofWPhysics::BuildSimplePhysicalEvent() {
 ///////////////////////////////////////////////////////////////////////////
 void TSofTofWPhysics::BuildPhysicalEvent() {
   // apply thresholds and calibration
+  if(m_StartTime == -1)
+    return;
+
   PreTreat();
 
-  // match energy and time together
+  double T1[28][10], T2[28][10];
+  int mult1[28], mult2[28];
+  for(int i=0; i<28; i++){
+    mult1[i] = 0;
+    mult2[i] = 0;
+    for(int j=0; j<10; j++){
+      T1[i][j] = 0;
+      T2[i][j] = 0;
+    }
+  }
+
   unsigned int mysizeE = m_PreTreatedData->GetMultiplicity();
-  for (UShort_t e = 0; e < mysizeE ; e++) {
-    PlasticNumber.push_back(m_PreTreatedData->GetPlasticNbr(e));
-    Energy.push_back(m_PreTreatedData->GetEnergy(e));
+  for (UShort_t i = 0; i < mysizeE ; i++) {
+    int plastic = m_PreTreatedData->GetPlasticNbr(i);
+    int pmt     = m_PreTreatedData->GetPmt(i);
+    int FT      = m_PreTreatedData->GetFineTime(i);
+    int CT      = m_PreTreatedData->GetCoarseTime(i);
+    
+    double T = CalculateTimeNs(plastic, pmt, FT, CT);
+
+    if(pmt==1){
+      T1[plastic-1][mult1[plastic-1]] = T;
+      mult1[plastic-1]++;
+    }
+    else if(pmt==2){
+      T2[plastic-1][mult2[plastic-1]] = T;    
+      mult2[plastic-1]++;
+    }
+  }
+  
+  static CalibrationManager* Cal = CalibrationManager::getInstance();
+  for(int p=0; p<m_NumberOfPlastics; p++){
+    if(mult1[p]==1 && mult2[p]==1){
+      for(int i=0; i<mult1[p]; i++){
+        for(int j=0; j<mult2[p]; j++){
+          double time_ns = 0.5*(T1[p][i] + T2[p][j]);
+          double rawpos = T1[p][i] - T2[p][j];
+          double calpos = Cal->ApplyCalibration("SofTofW/TOFW"+NPL::itoa(p+1)+"_POSPAR",rawpos);
+          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);
+          RawTof.push_back(rawtof);
+          CalTof.push_back(caltof);
+        }
+      }
+    }
   }
+  m_StartTime = -1;
+}
+
+///////////////////////////////////////////////////////////////////////////
+double TSofTofWPhysics::CalculateTimeNs(int det, int pmt, int ft, int ct){
+
+  static CalibrationManager* Cal = CalibrationManager::getInstance();
+  double par = Cal->GetValue("SofTofW/TOFW"+NPL::itoa(det)+"_PMT"+NPL::itoa(pmt)+"_TIME",ft);
+  double r = (double)rand.Rndm()-0.5;
+  double ClockOffset = Cal->GetValue("SofTofW/TOFW"+NPL::itoa(det)+"_PMT"+NPL::itoa(pmt)+"_CLOCKOFFSET",0);
+  double ict_ns = ((double)ct - ClockOffset) * 5.; 
+  double ift_ns;
+
+  if(r<0){
+    double par_prev = Cal->GetValue("SofTofW/TOFW"+NPL::itoa(det)+"_PMT"+NPL::itoa(pmt)+"_TIME",ft-1);
+    ift_ns = par + r*(par - par_prev);
+  }
+
+  else{
+    double par_next = Cal->GetValue("SofSci/TOFW"+NPL::itoa(det)+"_PMT"+NPL::itoa(pmt)+"_TIME",ft+1);
+    ift_ns = par + r*(par_next - par);
+  }
+
+  double time_ns = (double)ict_ns - ift_ns;
+
+  return time_ns;
 }
 
 ///////////////////////////////////////////////////////////////////////////
@@ -103,13 +180,11 @@ void TSofTofWPhysics::PreTreat() {
   // Energy
   unsigned int mysize = m_EventData->GetMultiplicity();
   for (UShort_t i = 0; i < mysize ; ++i) {
-    if (m_EventData->GetEnergy(i) > m_E_RAW_Threshold) {
-      Double_t Energy = Cal->ApplyCalibration("SofTofW/ENERGY"+NPL::itoa(m_EventData->GetPlasticNbr(i)),m_EventData->GetEnergy(i));
-      if (Energy > m_E_Threshold) {
-        m_PreTreatedData->SetPlasticNbr(m_EventData->GetPlasticNbr(i));
-        m_PreTreatedData->SetEnergy(Energy);
-      }
-    }
+    m_PreTreatedData->SetPlasticNbr(m_EventData->GetPlasticNbr(i));
+    m_PreTreatedData->SetPmt(m_EventData->GetPmt(i));
+    m_PreTreatedData->SetCoarseTime(m_EventData->GetCoarseTime(i));
+    m_PreTreatedData->SetFineTime(m_EventData->GetFineTime(i));
+    m_PreTreatedData->SetWhichFlag(m_EventData->GetWhichFlag(i));
   }
 }
 
@@ -181,10 +256,12 @@ void TSofTofWPhysics::ReadAnalysisConfig() {
 
 ///////////////////////////////////////////////////////////////////////////
 void TSofTofWPhysics::Clear() {
-  PlasticNumber.clear();
-  Energy.clear();
-  Time.clear();
-  PosY.clear();
+  PlasticNbr.clear();
+  TimeNs.clear();
+  RawPosY.clear();
+  CalPosY.clear();
+  RawTof.clear();
+  CalTof.clear();
 }
 
 
@@ -225,9 +302,14 @@ void TSofTofWPhysics::ReadConfiguration(NPL::InputParser parser) {
 ///////////////////////////////////////////////////////////////////////////
 void TSofTofWPhysics::AddParameterToCalibrationManager() {
   CalibrationManager* Cal = CalibrationManager::getInstance();
-  for (int i = 0; i < m_NumberOfDetectors; ++i) {
-    Cal->AddParameter("SofTofW", "D"+ NPL::itoa(i+1)+"_ENERGY","SofTofW_D"+ NPL::itoa(i+1)+"_ENERGY");
-    Cal->AddParameter("SofTofW", "D"+ NPL::itoa(i+1)+"_TIME","SofTofW_D"+ NPL::itoa(i+1)+"_TIME");
+  for (int i = 0; i < m_NumberOfPlastics; ++i) {
+    Cal->AddParameter("SofTofW", "TOFW"+ NPL::itoa(i+1)+"_POSPAR","SofTofW_TOFW"+ NPL::itoa(i+1)+"_POSPAR");
+    Cal->AddParameter("SofTofW", "TOFW"+ NPL::itoa(i+1)+"_TOFPAR","SofTofW_TOFW"+ NPL::itoa(i+1)+"_TOFPAR");
+    
+    for(int j = 0; j < 2; j++){
+      Cal->AddParameter("SofTofW", "TOFW"+ NPL::itoa(i+1)+"_PMT"+NPL::itoa(j+1)+"_TIME","SofTofW_TOFW"+ NPL::itoa(i+1)+"_PMT"+NPL::itoa(j+1)+"_TIME");
+      Cal->AddParameter("SofTofW", "TOFW"+ NPL::itoa(i+1)+"_PMT"+NPL::itoa(j+1)+"_CLOCKOFFSET","SofTofW_TOFW"+ NPL::itoa(i+1)+"_PMT"+NPL::itoa(j+1)+"_CLOCKOFFSET");
+    }
   }
 }
 
diff --git a/NPLib/Detectors/Sofia/TSofTofWPhysics.h b/NPLib/Detectors/Sofia/TSofTofWPhysics.h
index 3b8ac0370..94a402ec7 100644
--- a/NPLib/Detectors/Sofia/TSofTofWPhysics.h
+++ b/NPLib/Detectors/Sofia/TSofTofWPhysics.h
@@ -32,6 +32,7 @@ using namespace std;
 #include "TObject.h"
 #include "TH1.h"
 #include "TVector3.h"
+#include "TRandom3.h"
 // NPTool headers
 #include "TSofTofWData.h"
 #include "NPCalibrationManager.h"
@@ -59,10 +60,12 @@ class TSofTofWPhysics : public TObject, public NPL::VDetector {
   // data obtained after BuildPhysicalEvent() and stored in
   // output ROOT file
   public:
-    vector<int>      PlasticNumber;
-    vector<double>   Energy;
-    vector<double>   Time;
-    vector<double>   PosY;
+    vector<int>      PlasticNbr;
+    vector<double>   TimeNs;
+    vector<double>   RawPosY;
+    vector<double>   CalPosY;
+    vector<double>   RawTof;
+    vector<double>   CalTof;
 
   /// A usefull method to bundle all operation to add a detector
   void AddDetector(TVector3 POS); 
@@ -105,6 +108,11 @@ class TSofTofWPhysics : public TObject, public NPL::VDetector {
     void ClearEventPhysics() {Clear();}      
     void ClearEventData()    {m_EventData->Clear();}   
 
+    double CalculateTimeNs(int, int, int, int);
+    double GetStartTime() {return m_StartTime;}
+    double GetTofAlignedValue() {return m_TofAlignedValue;}
+    void SetStartTime(double val) {m_StartTime = val;}
+    void SetTofAlignedValue(double val) {m_TofAlignedValue = val;}
 
   //////////////////////////////////////////////////////////////
   // specific methods to SofTofW array
@@ -135,10 +143,13 @@ class TSofTofWPhysics : public TObject, public NPL::VDetector {
 
   // parameters used in the analysis
   private:
-    // thresholds
+    double m_StartTime; //!
+    double m_TofAlignedValue; //!
+    int m_NumberOfPlastics; //!
     double m_E_RAW_Threshold; //!
     double m_E_Threshold;     //!
 
+    TRandom3 rand; //!
   // number of detectors
   private:
     int m_NumberOfDetectors;  //!
diff --git a/Projects/s455/Analysis.cxx b/Projects/s455/Analysis.cxx
index 13471dc5b..1f8b2a72e 100644
--- a/Projects/s455/Analysis.cxx
+++ b/Projects/s455/Analysis.cxx
@@ -40,7 +40,7 @@ void Analysis::Init(){
   SofSci= (TSofSciPhysics*) m_DetectorManager->GetDetector("SofSci");
   SofTrim= (TSofTrimPhysics*) m_DetectorManager->GetDetector("SofTrim");
   SofTwim= (TSofTwimPhysics*) m_DetectorManager->GetDetector("SofTwim");
-  //SofTofW= (TSofTofWPhysics*) m_DetectorManager->GetDetector("SofTofW");
+  SofTofW= (TSofTofWPhysics*) m_DetectorManager->GetDetector("SofTofW");
 
   InitParameter();
   InitOutputBranch();
@@ -54,7 +54,13 @@ void Analysis::TreatEvent(){
   ReInitValue();
   //cout << "************" << endl;
   BeamAnalysis();
-
+  unsigned int sofsci_size = SofSci->DetectorNbr.size();
+  if(sofsci_size==2){
+    double start_time = SofSci->TimeNs[1];
+    SofTofW->SetTofAlignedValue(36);
+    SofTofW->SetStartTime(start_time);
+    SofTofW->BuildPhysicalEvent();
+  }
 }
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -165,8 +171,10 @@ void Analysis::End(){
 ////////////////////////////////////////////////////////////////////////////////
 void Analysis::InitParameter(){
   fLS2_0 = 136.3706933;
-  fDS2   = 9500;
-  fDCC   = -30000;
+  //fDS2   = 9500;
+  fDS2   = 10000;
+  //fDCC   = -30000;
+  fDCC   = -40000;
   fK_LS2 = -2.5e-8;
   fBrho0 = 10.8183; // run401 -> 182Hg
 
diff --git a/Projects/s455/calibration.txt b/Projects/s455/calibration.txt
index 383018a10..d77958d2c 100644
--- a/Projects/s455/calibration.txt
+++ b/Projects/s455/calibration.txt
@@ -14,6 +14,8 @@ CalibrationFilePath
 ./calibration/SofSci/SofSci_physics.cal
 ./calibration/SofTwim/SofTwim_Energy.cal
 ./calibration/SofTwim/SofTwim_Time.cal
+./calibration/SofTofW/ClockOffset.cal
+./calibration/SofTofW/SofTofW_physics.cal
 ./calibration/SofTofW/VFTX_TOFW1_PMT1.cal
 ./calibration/SofTofW/VFTX_TOFW2_PMT1.cal
 ./calibration/SofTofW/VFTX_TOFW3_PMT1.cal
diff --git a/Projects/s455/s455.detector b/Projects/s455/s455.detector
index b30b0fa9f..3cb3ec2ce 100644
--- a/Projects/s455/s455.detector
+++ b/Projects/s455/s455.detector
@@ -20,10 +20,9 @@ SofTrim
 SofTwim
  POS= 0 0 2 m
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-%SofTofW
-% R= 8 m
-% THETA= -9.5 deg
-% PHI= 0 deg
+SofTofW
+ R= 8 m
+ THETA= -9.5 deg
+ PHI= 0 deg
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-- 
GitLab