[PISTA] Adding IC calibration in PISTA 2/2

NPLib/Detectors/PISTA/CMakeLists.txt

@@ -21,5 +21,6 @@ add_library(NPPISTA SHARED TPISTASpectra.cxx TPISTAData.cxx TPISTAPhysics.cxx
 
 target_link_libraries(NPPISTA ${ROOT_LIBRARIES} NPCore) 
 
-install(FILES TPISTAData.h TPISTAPhysics.h TPISTASpectra.h TFPMWData.h TFPMWPhysics.h TICData.h TICPhysics.h TVamosReconstruction.h TTimeData.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
+install(FILES TPISTAData.h TPISTAPhysics.h TPISTASpectra.h TFPMWData.h TFPMWPhysics.h TICData.h TICPhysics.h
+    TVamosReconstruction.h TTimeData.h TProfileEvaluator.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
 

NPLib/Detectors/PISTA/TICPhysics.cxx

@@ -9,10 +9,10 @@
  * Original Author: P. Morfouace  contact address: pierre.morfouace@cea.fr   *
  *                                                                           *
  * Creation Date  : October 2023                                             *
- * Last update    :                                                          *
+ * Last update    : 22/01/24                                                 *
  *---------------------------------------------------------------------------*
  * Decription:                                                               *
- *  This class hold IC Treated  data                                    *
+ *  This class hold IC Treated  data                                         *
  *                                                                           *
  *---------------------------------------------------------------------------*
  * Comment:                                                                  *
@@ -38,6 +38,7 @@ using namespace std;
 
 //   ROOT
 #include "TChain.h"
+#include "TKey.h"
 
 ClassImp(TICPhysics);
 
@@ -47,7 +48,6 @@ TICPhysics::TICPhysics()
   m_PreTreatedData(new TICData),
   m_EventPhysics(this),
   m_FPMW_Section(-1),
-  m_number_of_splines(34),
   m_Eres_Threshold(3000),
   m_Z_SPLINE_CORRECTION(false),
   m_NumberOfDetectors(0){
@@ -94,11 +94,26 @@ void TICPhysics::BuildPhysicalEvent() {
     double GainInit = Cal->GetValue("IC/INIT_SEG"+NPL::itoa(segment)+"_ALIGN",0);
 
     fIC_raw[i] = m_PreTreatedData->GetIC_Charge(i);
+
     fIC[i] = gain*m_PreTreatedData->GetIC_Charge(i);
+    
+    if(i < 4){
+        fIC_PID[i] = 0.5* m_PreTreatedData->GetIC_Charge(i);
+    }
+    else if(i >= 6){
+        fIC_PID[i] = 2* m_PreTreatedData->GetIC_Charge(i);
+    }
+    else {
+        fIC_PID[i] = m_PreTreatedData->GetIC_Charge(i);
+    }
+
     fIC_Init[i] = GainInit * m_PreTreatedData->GetIC_Charge(i);
     fIC_TS.push_back(m_PreTreatedData->GetIC_TS(i));
   }
 
+  if(m_Y_SPLINE_CORRECTION && m_XY0_SPLINE_CORRECTION){
+      ApplyXYCorrections(); 
+       }
 
   if (fIC_Init[1]>0 && fIC_Init[5]>0) {
     EtotInit = 0 ;
@@ -121,10 +136,26 @@ void TICPhysics::BuildPhysicalEvent() {
     EtotInit = -100 ;
   }
 
-  if(fIC[1]>0 && fIC[5]>0){
-    DE = 0.5*(fIC_raw[0] + fIC_raw[1] + fIC_raw[2] + fIC_raw[3]) + fIC_raw[4];
-    Eres = fIC_raw[5] + fIC_raw[6] + fIC_raw[7] + fIC_raw[8] + fIC_raw[9];
+    DE = (fIC_PID[0] + fIC_PID[1] + fIC_PID[2] + fIC_PID[3]) + fIC_PID[4];
+    Eres = fIC_PID[5] + (fIC_PID[6] + fIC_PID[7] + fIC_PID[8] + fIC_PID[9]);
+
+    if (m_DE_SPLINE_CORRECTION && m_Y_SPLINE_CORRECTION){
+        if (m_TimeData->GetMWPC13Mult() ==1 && fIC_TS.size()>=8 ){ //only mult 1 event
+            UShort_t FPMW_Section = m_FPMW_Section;
+            double TempY;
+
+            //Data year sensitive loading
+            if (m_Data_Year == 2024){
+               TempY  =  10* (fIC_TS.at(0) - m_TimeData->GetTS_MWPC13(0)) - ((m_TimeData->GetTime_MWPC13(0)+m_TimeData->GetToff_DT13(FPMW_Section))) ;
+            }
+            else if (m_Data_Year == 2023){
+                TempY =  m_Y ;
+            }
+            DE = DE * m_DEspline.at(0)->Eval(0) / m_DEspline.at(0)->Eval(TempY) ;
+        } // end if mult 1
+    } // end DE correction
 
+  if(fIC[1]>0 && fIC[5]>0){
     double scalor = Cal->GetValue("IC/ETOT_SCALING_SEC"+NPL::itoa(m_FPMW_Section),0);
 
     for(int i=0; i<10; i++){
@@ -232,17 +263,39 @@ void TICPhysics::ReadAnalysisConfig() {
         m_Data_Year = atof(DataBuffer.c_str());
         cout << whatToDo << " " << m_Data_Year << endl;
       }
+      else if (whatToDo=="LOAD_Y_SPLINE") {
+        AnalysisConfigFile >> DataBuffer;
+        m_Y_SPLINE_PATH = DataBuffer;
+        cout << "*** Loading Y spline ***" << endl;
+        m_Y_SPLINE_CORRECTION = LoadSpline(m_Yspline,m_number_Y_spline,m_Y_SPLINE_PATH);
+      }
+      else if (whatToDo=="LOAD_XY0_PROFILE") {
+        AnalysisConfigFile >> DataBuffer;
+        m_XY0_PROFILE_PATH = DataBuffer;
+        cout << "*** Loading XY0 profile ***" << endl;
+        TString PathTemp = m_XY0_PROFILE_PATH;
+        TFile* ifile = new TFile(PathTemp,"read");
+        if(ifile->IsOpen() && !ifile->IsZombie()){
+            m_XY0_SPLINE_CORRECTION = true;
+            m_IC0_Profile.LoadProfile(PathTemp,"ICOneZeroProfile");
+        }
+        else {
+            m_XY0_SPLINE_CORRECTION = false ;
+        }
+        ifile->Close();
+      }
 
-      else if(whatToDo=="NUMBER_OF_SPLINES"){
+      else if (whatToDo=="LOAD_DE_SPLINE") {
         AnalysisConfigFile >> DataBuffer;
-        m_number_of_splines = atoi(DataBuffer.c_str());
-        cout << whatToDo << " " << m_number_of_splines << endl;
+        m_DE_SPLINE_PATH = DataBuffer;
+        cout << "*** Loading DE spline ***" << endl;
+        m_DE_SPLINE_CORRECTION = LoadSpline(m_DEspline,m_number_DE_spline,m_DE_SPLINE_PATH);
       }
       else if (whatToDo=="LOAD_Z_SPLINE"){
         AnalysisConfigFile >> DataBuffer;
         m_Z_SPLINE_PATH = DataBuffer;
         cout << "*** Loading Z spline ***" << endl;
-        LoadZSpline();
+        m_Z_SPLINE_CORRECTION = LoadSpline(m_Zspline,m_number_zspline,m_Z_SPLINE_PATH);
       }
       else {
         ReadingStatus = false;
@@ -253,35 +306,51 @@ void TICPhysics::ReadAnalysisConfig() {
 
 
 ///////////////////////////////////////////////////////////////////////////
-void TICPhysics::LoadZSpline(){
-  TString filename = m_Z_SPLINE_PATH;
+bool TICPhysics::LoadSpline(vector<TSpline3*> &iSpline, int &NSpline , string Path){
+  TString filename = Path;
   TFile* ifile = new TFile(filename,"read");
+  if(ifile->IsOpen() && !ifile->IsZombie()){
+
+    // Get number of spline
+    TIter next(ifile->GetListOfKeys());
+    TKey* key;
+    NSpline = 0 ;
 
-  if(ifile->IsOpen()){
-    m_Z_SPLINE_CORRECTION = true;
-    for(int i=0; i<m_number_of_splines; i++){
-      m_Zspline[i] = (TSpline3*) ifile->FindObjectAny(Form("fspline_%d",i+1));
-      cout << "*** " << m_Zspline[i]->GetName() << " is loaded!" << endl;
+    while ((key=(TKey*)next())){
+        if (std::string(key->GetClassName()) == "TSpline3"){
+            NSpline ++;
+        }
+    }
+    cout << "This file contains  " << NSpline << " splines "  << endl;
+    // Load Spline
+    for(int i=0; i<NSpline; i++){
+      iSpline.at(i) = (TSpline3*) ifile->FindObjectAny(Form("fspline_%d",i+1));
+      iSpline.at(i)->SetName(Form("fspline_%s_%d",Path.c_str(),i+1));
+      cout  << iSpline.at(i)->GetName() << " is loaded!" << endl;
     }
 
-  }
-  else
-    cout << "File " << filename << " not found!" << endl;
   ifile->Close();
+  return true;
+  }
+  else{
+      cout << "File " << filename << " not found!" << endl;
+      ifile->Close();
+      return false;
+  }
 }
 
 ///////////////////////////////////////////////////////////////////////////
 double TICPhysics::ApplyZSpline(){
   double DEcorr;
-  double FF_DEcorr0[m_number_of_splines];
-  for(int i=0; i<m_number_of_splines; i++){
+  double FF_DEcorr0[m_number_zspline];
+  for(int i=0; i<m_number_zspline; i++){
     FF_DEcorr0[i] = m_Zspline[i]->Eval(8500);
   }
 
   double DEspline0;
   double Eval_DEspline;
   int index=0;
-  for(int i=0; i<m_number_of_splines; i++){
+  for(int i=0; i<m_number_zspline; i++){
     Eval_DEspline = m_Zspline[i]->Eval(Eres);
     if(DE<Eval_DEspline) break;
     index = i;
@@ -290,7 +359,7 @@ double TICPhysics::ApplyZSpline(){
   Eval_DEspline = m_Zspline[index]->Eval(Eres);
   DEspline0 = FF_DEcorr0[index];
   double dmin, dsup;
-  if(index<(m_number_of_splines-1) && DE>m_Zspline[0]->Eval(Eres)){
+  if(index<(m_number_zspline-1) && DE>m_Zspline[0]->Eval(Eres)){
     dmin = DE - m_Zspline[index]->Eval(Eres);
     dsup = m_Zspline[index+1]->Eval(Eres) - DE;
 
@@ -299,7 +368,7 @@ double TICPhysics::ApplyZSpline(){
 
     DEcorr        = DEspline0 * DE / Eval_DEspline;
   }
-  else if(index==m_number_of_splines-1){
+  else if(index==m_number_zspline-1){
     Eval_DEspline = m_Zspline[index]->Eval(Eres);
     DEspline0     = FF_DEcorr0[index];
 
@@ -308,6 +377,66 @@ double TICPhysics::ApplyZSpline(){
 
   return DEcorr;
 }
+///////////////////////////////////////////////////////////////////////////
+void TICPhysics::ApplyXYCorrections(){
+    vector<double> ICcorr_Y(11), ICcorr_X(11); 
+    double FF_DriftTime ;
+    if (m_TimeData->GetMWPC13Mult() ==1 && fIC_TS.size()>=8){ //only mult 1 event
+        UShort_t FPMW_Section = m_FPMW_Section;
+
+        // ***************************Different def of correction depending on year**********************************
+        if (m_Data_Year == 2024){
+            FF_DriftTime =  10* (fIC_TS.at(0) - m_TimeData->GetTS_MWPC13(0)) - ((m_TimeData->GetTime_MWPC13(0)+m_TimeData->GetToff_DT13(FPMW_Section))) ;
+        }
+        else if (m_Data_Year == 2023){
+            FF_DriftTime = m_Y;
+        }
+        else {
+            return ;
+        }
+
+        //**************************  Correction of section 1 to 4 ***************************************************
+        for (int seg = 2; seg < fIC_TS.size() ; seg++) { // loop on multiplicity of event
+
+            if (m_Yspline.at(seg-2)==0){
+                ICcorr_Y.at(seg) = fIC_PID[seg]; 
+            }
+
+            else {
+                ICcorr_Y.at(seg) = fIC_PID[seg] * m_Yspline.at(seg-2)->Eval(0)/ m_Yspline.at(seg-2)->Eval(FF_DriftTime);
+                if (!(ICcorr_Y.at(seg)==ICcorr_Y.at(seg))) ICcorr_Y.at(seg) = 0;
+            } //end if non empty
+            if (seg == 0) break;
+        }//endloop seg
+
+
+        //**************************  Correction of section 0 ***************************************************
+
+        Double_t PolX =  1.37622;
+
+        Double_t ICRatio =  fIC_PID[1]/fIC_PID[0];
+        Double_t ICRatioCorr =  ICRatio * PolX / m_IC0_Profile.Evaluate(m_X,FF_DriftTime,false);
+        Double_t ICcorr_Y0 = fIC_PID[0] / PolX * m_IC0_Profile.Evaluate(m_X,FF_DriftTime,false);
+
+        if ( ICRatioCorr<1.4 || ICRatioCorr >1.5){
+            ICRatioCorr =  ICRatio * PolX / m_IC0_Profile.Evaluate(m_X,FF_DriftTime,false);
+            ICcorr_Y0 =     fIC_PID[0] / PolX * m_IC0_Profile.Evaluate(m_X,FF_DriftTime,false);
+            if (ICRatioCorr >100) {
+                ICcorr_Y0 = fIC_PID[0];
+            }
+        }
+        
+
+        //**************************  Overwrite ICRAW ***************************************************
+        //Overwrite ic_raw
+        for (int i = 1 ; i<fIC_TS.size() ;i++){
+            if (ICcorr_Y.at(i) != 0 ){
+                fIC_PID[i] = ICcorr_Y.at(i);
+            }
+        }
+        fIC_PID[0] = ICcorr_Y0;
+    }
+}
 
 ///////////////////////////////////////////////////////////////////////////
 void TICPhysics::Clear() {
@@ -320,6 +449,7 @@ void TICPhysics::Clear() {
     fIC[i] = 0;
     fIC_raw[i] = 0;
     fIC_Init[i] = 0;
+    fIC_PID[i] = 0;
   }
 
   fIC_TS.clear();

NPLib/Detectors/PISTA/TICPhysics.h

@@ -34,8 +34,11 @@ using namespace std;
 #include "TH1.h"
 #include "TVector3.h"
 #include "TSpline.h"
+#include "TFile.h"
 // NPTool headers
 #include "TICData.h"
+#include "TTimeData.h"
+#include "TProfileEvaluator.h"
 #include "NPCalibrationManager.h"
 #include "NPVDetector.h"
 #include "NPInputParser.h"
@@ -71,6 +74,7 @@ class TICPhysics : public TObject, public NPL::VDetector {
 
     double fIC[11];
     double fIC_raw[11]; 
+    double fIC_PID[11]; 
     double fIC_Init[11];//! 
   private:
 
@@ -132,10 +136,22 @@ class TICPhysics : public TObject, public NPL::VDetector {
     // needed for online analysis for example
     void SetRawDataPointer(TICData* rawDataPointer) {m_EventData = rawDataPointer;}
 
+    // Setter and getter for section of the FPMW.
     void SetFPMWSection(int section) {m_FPMW_Section = section;}
     int GetFPMWSection() {return m_FPMW_Section;}
-    void LoadZSpline();
+    
+    //Setter and Getter for TTimeData used for drift time calculation and XY
+    void SetTTimeData(TTimeData *newTime) {m_TimeData = newTime;}
+    TTimeData* GetTTimeData() {return m_TimeData;}
+    
+    void SetX(double PosX) {m_X = PosX;}
+    double GetX() {return m_X;}
+    void SetY(double PosY) {m_Y = PosY;}
+    double GetY() {return m_Y;}
+
+    bool LoadSpline(vector<TSpline3*> &iSpline, int &NSpline, string Path);
     double ApplyZSpline();
+    void ApplyXYCorrections();
 
   // objects are not written in the TTree
   private:
@@ -152,14 +168,36 @@ class TICPhysics : public TObject, public NPL::VDetector {
   private:
     int m_NumberOfDetectors;  //!
     int m_FPMW_Section; //!
-    int m_number_of_splines; //!
+    int m_number_zspline; //!
     string m_Z_SPLINE_PATH; //!
     bool m_Z_SPLINE_CORRECTION; //!
-    TSpline3* m_Zspline[50]; //!
+    vector<TSpline3*>  m_Zspline = vector<TSpline3*>(50); //!
     double m_Eres_Threshold; //!
 
-   //Year of acquisition
+    //Time temporary variable
+    TTimeData*      m_TimeData;  //!
+    //Correction in XY of IC
+    string m_Y_SPLINE_PATH; //!
+    string m_XY0_PROFILE_PATH; //!
+    string m_DE_SPLINE_PATH; //!
+
+    bool m_Y_SPLINE_CORRECTION = false; //!
+    bool m_DE_SPLINE_CORRECTION = false; //!
+    bool m_XY0_SPLINE_CORRECTION = false; //!
+
+    int m_number_Y_spline; //!
+    int m_number_XY0_spline; //!
+    int m_number_DE_spline; //!
+
+    double m_X; //!
+    double m_Y; //!
+    // Be careful this array contains only the correction from IC1 to 4
+    vector<TSpline3*> m_Yspline = vector<TSpline3*>(11); //! 
+    vector<TSpline3*> m_DEspline = vector<TSpline3*>(11); //! 
+    // The following contains correction for IC0
+    ProfileEvaluator m_IC0_Profile;//!
   private:
+    //Year of acquisition
     double m_Data_Year; //!
   // Static constructor to be passed to the Detector Factory
   public:

NPLib/Detectors/PISTA/TTimeData.cxx

@@ -54,6 +54,9 @@ void TTimeData::Clear() {
   fTime_MWPC23.clear();
   fTime_MWPC24.clear();
 
+  fToff_DT13.clear();
+  fToff_DT14.clear();
+
   fSection_MWPC3.clear();
   fSection_MWPC4.clear();
 

NPLib/Detectors/PISTA/TTimeData.h

@@ -44,6 +44,9 @@ class TTimeData : public TObject {
         vector<float> fTime_MWPC23;
         vector<float> fTime_MWPC24;
 
+        vector<double> fToff_DT13;
+        vector<double> fToff_DT14;
+
         vector<short> fSection_MWPC3;
         vector<short> fSection_MWPC4;
 
@@ -82,6 +85,9 @@ class TTimeData : public TObject {
         inline void SetTime_MWPC23(float time ){fTime_MWPC23.push_back(time);};//!
         inline void SetTime_MWPC24(float time ){fTime_MWPC24.push_back(time);};//!
 
+        inline void SetToff_DT13(double offset ){fToff_DT13.push_back(offset);};//!
+        inline void SetToff_DT14(double offset ){fToff_DT14.push_back(offset);};//!
+
         inline void SetSection_MWPC3(short section ){fSection_MWPC3.push_back(section);};//!
         inline void SetSection_MWPC4(short section ){fSection_MWPC4.push_back(section);};//!
 
@@ -107,6 +113,11 @@ class TTimeData : public TObject {
         inline float GetTime_MWPC24(const unsigned int &i) const
         {return fTime_MWPC24.at(i) ;}//!
         
+        inline double GetToff_DT13(const unsigned int &i) const
+        {return fToff_DT13.at(i) ;}//!
+        inline double GetToff_DT14(const unsigned int &i) const
+        {return fToff_DT14.at(i) ;}//!
+          
         inline short GetSection_MWPC3(const unsigned int &i) const
         {return fSection_MWPC3.at(i) ;}//!
         inline short GetSection_MWPC4(const unsigned int &i) const
@@ -123,7 +134,7 @@ class TTimeData : public TObject {
 
         //////////////////////////////////////////////////////////////
         // Required for ROOT dictionnary
-        ClassDef(TTimeData,2)  // TimeData structure
+        ClassDef(TTimeData,3)  // TimeData structure
 };
 
 #endif