diff --git a/NPLib/Detectors/Sofia/CMakeLists.txt b/NPLib/Detectors/Sofia/CMakeLists.txt
index 10ce558b9212cb0f6a69dc28c96d7d485c039d58..e9efe35c7f17a91031d6e9b4dd9d87545d2c2b96 100644
--- a/NPLib/Detectors/Sofia/CMakeLists.txt
+++ b/NPLib/Detectors/Sofia/CMakeLists.txt
@@ -1,4 +1,5 @@
add_custom_command(OUTPUT TSofSciDataDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TSofSciData.h TSofSciDataDict.cxx TSofSciData.rootmap libNPSofia.dylib DEPENDS TSofSciData.h)
+add_custom_command(OUTPUT TSofSciPhysicsDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TSofSciPhysics.h TSofSciPhysicsDict.cxx TSofSciPhysics.rootmap libNPSofia.dylib DEPENDS TSofSciPhysics.h)
add_custom_command(OUTPUT TSofTofWDataDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TSofTofWData.h TSofTofWDataDict.cxx TSofTofWData.rootmap libNPSofia.dylib DEPENDS TSofTofWData.h)
add_custom_command(OUTPUT TSofTofWPhysicsDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TSofTofWPhysics.h TSofTofWPhysicsDict.cxx TSofTofWPhysics.rootmap libNPSofia.dylib DEPENDS TSofTofWPhysics.h)
@@ -12,9 +13,11 @@ add_custom_command(OUTPUT TSofTrimPhysicsDict.cxx COMMAND ${CMAKE_BINARY_DIR}/sc
add_custom_command(OUTPUT TSofTwimDataDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TSofTwimData.h TSofTwimDataDict.cxx TSofTwimData.rootmap libNPSofia.dylib DEPENDS TSofTwimData.h)
-add_library(NPSofia SHARED TSofSciData.cxx TSofSciDataDict.cxx TSofMwpcData.cxx TSofMwpcDataDict.cxx TSofAtData.cxx TSofAtDataDict.cxx TSofTrimData.cxx TSofTrimDataDict.cxx TSofTrimPhysics.cxx TSofTrimPhysicsDict.cxx TSofTwimData.cxx TSofTwimDataDict.cxx TSofTofWData.cxx TSofTofWDataDict.cxx TSofTofWPhysics.cxx TSofTofWPhysicsDict.cxx)
+add_custom_command(OUTPUT TSofBeamIDDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TSofBeamID.h TSofBeamIDDict.cxx TSofBeamID.rootmap libNPSofia.dylib DEPENDS TSofBeamID.h)
+
+add_library(NPSofia SHARED TSofSciData.cxx TSofSciDataDict.cxx TSofSciPhysics.cxx TSofSciPhysicsDict.cxx TSofMwpcData.cxx TSofMwpcDataDict.cxx TSofAtData.cxx TSofAtDataDict.cxx TSofTrimData.cxx TSofTrimDataDict.cxx TSofTrimPhysics.cxx TSofTrimPhysicsDict.cxx TSofTwimData.cxx TSofTwimDataDict.cxx TSofTofWData.cxx TSofTofWDataDict.cxx TSofTofWPhysics.cxx TSofTofWPhysicsDict.cxx TSofBeamID.cxx TSofBeamIDDict.cxx)
target_link_libraries(NPSofia ${ROOT_LIBRARIES} NPCore NPPhysics)
-install(FILES TSofSciData.h TSofMwpcData.h TSofAtData.h TSofTrimData.h TSofTrimPhysics.h TSofTwimData.h TSofTofWData.h TSofTofWPhysics.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
+install(FILES TSofSciData.h TSofSciPhysics.h TSofMwpcData.h TSofAtData.h TSofTrimData.h TSofTrimPhysics.h TSofTwimData.h TSofTofWData.h TSofTofWPhysics.h TSofBeamID.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
diff --git a/NPLib/Detectors/Sofia/TSofBeamID.cxx b/NPLib/Detectors/Sofia/TSofBeamID.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..fca40725bae41768c9b33f054a4daea835c51da0
--- /dev/null
+++ b/NPLib/Detectors/Sofia/TSofBeamID.cxx
@@ -0,0 +1,74 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2020 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: Pierre Morfouace contact address: pierre.morfouace2@cea.fr *
+ * Creation Date : May 2021 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class hold SofBeamID Raw data *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+#include "TSofBeamID.h"
+
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include <string>
+using namespace std;
+
+ClassImp(TSofBeamID)
+
+
+//////////////////////////////////////////////////////////////////////
+TSofBeamID::TSofBeamID() {
+}
+
+
+
+//////////////////////////////////////////////////////////////////////
+TSofBeamID::~TSofBeamID() {
+}
+
+
+
+//////////////////////////////////////////////////////////////////////
+void TSofBeamID::Clear() {
+ Zbeam = -1;
+ Qmax = -1;
+ AoQ = -1;
+ Abeam = -1;
+ Beta = -1;
+ Gamma = -1;
+ Brho = -1;
+ XS2 = -1000;
+ XCC = -1000;
+
+}
+
+
+
+//////////////////////////////////////////////////////////////////////
+void TSofBeamID::Dump() const {
+ // This method is very useful for debuging and worth the dev.
+ cout << "XXXXXXXXXXXXXXXXXXXXXXXX New Event [TSofBeamID::Dump()] XXXXXXXXXXXXXXXXX" << endl;
+
+ cout << "Zbeam: " << Zbeam << endl;
+ cout << "AoQ: " << AoQ << endl;
+ cout << "Abeam: " << Abeam << endl;
+ cout << "Beta: " << Beta << endl;
+ cout << "Gamma: " << Gamma << endl;
+ cout << "Brho: " << Brho << endl;
+ cout << "XS2: " << XS2 << endl;
+ cout << "XCC: " << XCC << endl;
+
+}
diff --git a/NPLib/Detectors/Sofia/TSofBeamID.h b/NPLib/Detectors/Sofia/TSofBeamID.h
new file mode 100644
index 0000000000000000000000000000000000000000..72b52ad7a2f4423e5ff37fcbd555d4a8b23fb86f
--- /dev/null
+++ b/NPLib/Detectors/Sofia/TSofBeamID.h
@@ -0,0 +1,95 @@
+#ifndef __SofBeamIDDATA__
+#define __SofBeamIDDATA__
+/*****************************************************************************
+ * Copyright (C) 2009-2020 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: Pierre Morfouace contact address: pierre.morfouace2@cea.fr *
+ * *
+ * Creation Date : May 2021 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class hold SofBeamID Raw data *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+
+// STL
+#include <vector>
+using namespace std;
+
+// ROOT
+#include "TObject.h"
+
+class TSofBeamID : public TObject {
+ //////////////////////////////////////////////////////////////
+ // data members are hold into vectors in order
+ // to allow multiplicity treatment
+ private:
+ double Zbeam;
+ double Qmax;
+ double AoQ;
+ double Abeam;
+ double Beta;
+ double Gamma;
+ double Brho;
+ double XS2;
+ double XCC;
+
+ //////////////////////////////////////////////////////////////
+ // Constructor and destructor
+ public:
+ TSofBeamID();
+ ~TSofBeamID();
+
+
+ //////////////////////////////////////////////////////////////
+ // Inherited from TObject and overriden to avoid warnings
+ public:
+ void Clear();
+ void Clear(const Option_t*) {};
+ void Dump() const;
+
+
+ //////////////////////////////////////////////////////////////
+ // Getters and Setters
+ // Prefer inline declaration to avoid unnecessary called of
+ // frequently used methods
+ // add //! to avoid ROOT creating dictionnary for the methods
+ public:
+ ////////////////////// SETTERS ////////////////////////
+ inline void SetZbeam(double val){Zbeam = val;};//!
+ inline void SetQmax(double val){Qmax = val;};//!
+ inline void SetAoQ(double val){AoQ = val;};//!
+ inline void SetAbeam(double val){Abeam = val;};//!
+ inline void SetBeta(double val){Beta = val;};//!
+ inline void SetGamma(double val){Gamma = val;};//!
+ inline void SetBrho(double val){Brho = val;};//!
+ inline void SetXS2(double val){XS2 = val;};//!
+ inline void SetXCC(double val){XCC = val;};//!
+
+ ////////////////////// GETTERS ////////////////////////
+ inline double GetZbeam() const {return Zbeam;}//!
+ inline double GetQmax() const {return Qmax;}//!
+ inline double GetAoQ() const {return AoQ;}//!
+ inline double GetAbeam() const {return Abeam;}//!
+ inline double GetBeta() const {return Beta;}//!
+ inline double GetGamma() const {return Gamma;}//!
+ inline double GetBrho() const {return Brho;}//!
+ inline double GetXS2() const {return XS2;}//!
+ inline double GetXCC() const {return XCC;}//!
+
+ //////////////////////////////////////////////////////////////
+ // Required for ROOT dictionnary
+ ClassDef(TSofBeamID,1) // SofBeamID structure
+};
+
+#endif
diff --git a/NPLib/Detectors/Sofia/TSofSciData.h b/NPLib/Detectors/Sofia/TSofSciData.h
index 460e11ab92b86f1aff073ae747a1bdec83db23fa..e17540e5ff9ff22edb4d2ed3cc9f62fc34a1a491 100644
--- a/NPLib/Detectors/Sofia/TSofSciData.h
+++ b/NPLib/Detectors/Sofia/TSofSciData.h
@@ -34,10 +34,10 @@ class TSofSciData : public TObject {
// data members are hold into vectors in order
// to allow multiplicity treatment
private:
- vector<int> fSofSci_DetNbr;
- vector<int> fSofSci_Pmt;
- vector<double> fSofSci_CT;
- vector<double> fSofSci_FT;
+ vector<int> fSofSci_DetNbr;
+ vector<int> fSofSci_Pmt;
+ vector<int> fSofSci_CT;
+ vector<int> fSofSci_FT;
//////////////////////////////////////////////////////////////
// Constructor and destructor
@@ -63,15 +63,15 @@ class TSofSciData : public TObject {
////////////////////// SETTERS ////////////////////////
inline void SetDetectorNbr(int det){fSofSci_DetNbr.push_back(det);};//!
inline void SetPmt(int pmt){fSofSci_Pmt.push_back(pmt);};//!
- inline void SetCoarseTime(double Time){fSofSci_CT.push_back(Time);};//!
- inline void SetFineTime(double Time){fSofSci_FT.push_back(Time);};//!
+ inline void SetCoarseTime(int Time){fSofSci_CT.push_back(Time);};//!
+ inline void SetFineTime(int Time){fSofSci_FT.push_back(Time);};//!
////////////////////// GETTERS ////////////////////////
inline int GetMultiplicity() const {return fSofSci_DetNbr.size();}//!
inline int GetDetectorNbr(const unsigned int &i) const {return fSofSci_DetNbr[i];}//!
inline int GetPmt(const unsigned int &i) const {return fSofSci_Pmt[i];}//!
- inline double GetCoarseTime(const unsigned int &i) const {return fSofSci_CT[i];}//!
- inline double GetFineTime(const unsigned int &i) const {return fSofSci_FT[i];}//!
+ inline int GetCoarseTime(const unsigned int &i) const {return fSofSci_CT[i];}//!
+ inline int GetFineTime(const unsigned int &i) const {return fSofSci_FT[i];}//!
//////////////////////////////////////////////////////////////
// Required for ROOT dictionnary
diff --git a/NPLib/Detectors/Sofia/TSofSciPhysics.cxx b/NPLib/Detectors/Sofia/TSofSciPhysics.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..0d869d13c500b146ccd604c457cb6fa2f5d5d26b
--- /dev/null
+++ b/NPLib/Detectors/Sofia/TSofSciPhysics.cxx
@@ -0,0 +1,485 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2020 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: Pierre Morfouace contact address: pierre.morfouace2@cea.fr *
+ * *
+ * Creation Date : November 2020 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class hold SofSci Treated data *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+
+#include "TSofSciPhysics.h"
+
+// STL
+#include <sstream>
+#include <iostream>
+#include <cmath>
+#include <stdlib.h>
+#include <limits>
+using namespace std;
+
+// NPL
+#include "RootInput.h"
+#include "RootOutput.h"
+#include "NPDetectorFactory.h"
+#include "NPOptionManager.h"
+#include "NPPhysicalConstants.h"
+#include "NPGlobalSystemOfUnits.h"
+
+// ROOT
+#include "TChain.h"
+
+ClassImp(TSofSciPhysics)
+
+
+ ///////////////////////////////////////////////////////////////////////////
+TSofSciPhysics::TSofSciPhysics()
+ : m_EventData(new TSofSciData),
+ m_PreTreatedData(new TSofSciData),
+ m_EventPhysics(this),
+ m_DET1_PosNs_Min(-20),
+ m_DET1_PosNs_Max(20),
+ m_DET2_PosNs_Min(-20),
+ m_DET2_PosNs_Max(20),
+ m_RawTof_Min(0),
+ m_RawTof_Max(2000),
+ m_NumberOfSignals(3),
+ m_NumberOfPmts(2),
+ m_NumberOfDetectors(0){
+ }
+
+///////////////////////////////////////////////////////////////////////////
+/// A usefull method to bundle all operation to add a detector
+void TSofSciPhysics::AddDetector(TVector3 ){
+ // In That simple case nothing is done
+ // Typically for more complex detector one would calculate the relevant
+ // positions (stripped silicon) or angles (gamma array)
+ m_NumberOfDetectors++;
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TSofSciPhysics::AddDetector(double R, double Theta, double Phi){
+ // Compute the TVector3 corresponding
+ TVector3 Pos(R*sin(Theta)*cos(Phi),R*sin(Theta)*sin(Phi),R*cos(Theta));
+ // Call the cartesian method
+ AddDetector(Pos);
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TSofSciPhysics::BuildSimplePhysicalEvent() {
+ BuildPhysicalEvent();
+}
+
+
+
+///////////////////////////////////////////////////////////////////////////
+void TSofSciPhysics::BuildPhysicalEvent() {
+ Clear();
+ // apply thresholds and calibration
+ //PreTreat();
+
+ const int N = m_NumberOfDetectors;
+
+ vector<double> S2_pmtL;
+ vector<double> S2_pmtR;
+ vector<double> S2_pmtTref;
+ vector<double> S2_pos;
+ vector<double> S2_time;
+
+ vector<double> CC_pmtL;
+ vector<double> CC_pmtR;
+ vector<double> CC_pmtTref;
+ vector<double> CC_pos;
+ vector<double> CC_time;
+
+ unsigned int mysizeE = m_EventData->GetMultiplicity();
+ for (UShort_t e = 0; e < mysizeE ; e++) {
+ int det = m_EventData->GetDetectorNbr(e);
+ int pmt = m_EventData->GetPmt(e);
+ int FT = m_EventData->GetFineTime(e);
+ int CT = m_EventData->GetCoarseTime(e);
+
+ double T = CalculateTimeNs(det, pmt, FT, CT);
+
+ if(m_NumberOfDetectors==2){
+ if(det==1){
+ if(pmt==1) S2_pmtR.push_back(T);
+ else if(pmt==2) S2_pmtL.push_back(T);
+ else if(pmt==3) S2_pmtTref.push_back(T);
+ }
+
+ else if(det==2){
+ if(pmt==1) CC_pmtR.push_back(T);
+ else if(pmt==2) CC_pmtL.push_back(T);
+ else if(pmt==3) CC_pmtTref.push_back(T);
+ }
+ }
+
+ else if(m_NumberOfDetectors==1){
+ if(pmt==1) CC_pmtR.push_back(T);
+ else if(pmt==2) CC_pmtL.push_back(T);
+ else if(pmt==3) CC_pmtTref.push_back(T);
+ }
+ }
+
+ if(m_NumberOfDetectors==2){
+ double CC_rawpos;
+ double S2_rawpos;
+ double CC_calpos;
+ double S2_calpos;
+ double CC_rawtime;
+ double S2_rawtime;
+ double rawtof;
+
+ static CalibrationManager* Cal = CalibrationManager::getInstance();
+ if(S2_pmtTref.size()==1 && CC_pmtTref.size()==1){
+ for(unsigned int i=0; i<CC_pmtR.size(); i++){
+ for(unsigned int j=0; j<CC_pmtL.size(); j++){
+ CC_rawpos = CC_pmtR[i] - CC_pmtL[j];
+ CC_calpos = Cal->ApplyCalibration("SofSci/DET2_POSPAR", CC_rawpos);
+ CC_rawtime = 0.5*(CC_pmtR[i] + CC_pmtL[j]);
+ if(CC_rawpos<m_DET2_PosNs_Min || CC_rawpos>m_DET2_PosNs_Max)
+ continue;
+ for(int p=0; p<S2_pmtR.size(); p++){
+ for(int k=0; k<S2_pmtL.size(); k++){
+ S2_rawpos = S2_pmtR[p] - S2_pmtL[k];
+ S2_calpos = Cal->ApplyCalibration("SofSci/DET1_POSPAR", S2_rawpos);
+ S2_rawtime = 0.5*(S2_pmtR[p] + S2_pmtL[k]);
+
+ if(S2_rawpos<m_DET1_PosNs_Min || S2_rawpos>m_DET1_PosNs_Max)
+ continue;
+
+ /*CC_pos.push_back(CC_rawpos);
+ CC_time.push_back(CC_rawtime);
+ S2_pos.push_back(S2_rawpos);
+ S2_time.push_back(S2_rawtime); */
+
+ rawtof = CC_rawtime - CC_pmtTref[0] - (S2_rawtime - S2_pmtTref[0]);
+
+ if(rawtof<m_RawTof_Min || rawtof>m_RawTof_Max)
+ continue;
+
+ double velocity_mns;
+ double caltof;
+ double betaS2;
+
+ /*cout << "*** Printing physics calibration parameter:" << endl;
+ cout << Cal->GetValue("SofSci/TOF2INV_V",0) << endl;
+ cout << Cal->GetValue("SofSci/TOF2INV_V",1) << endl;
+ cout << Cal->GetValue("SofSci/LENGTH_S2",0) << endl;*/
+
+ velocity_mns = 1./Cal->ApplyCalibration("SofSci/TOF2INV_V",rawtof);
+ caltof = Cal->GetValue("SofSci/LENGTH_S2",0) / velocity_mns;
+ betaS2 = velocity_mns * m/ns / NPUNITS::c_light;
+
+ // Filling ouput Tree;
+ DetectorNbr.push_back(1);
+ TimeNs.push_back(S2_rawtime);
+ PosNs.push_back(S2_rawpos);
+ PosMm.push_back(S2_calpos);
+
+ DetectorNbr.push_back(2);
+ TimeNs.push_back(CC_rawtime);
+ PosNs.push_back(CC_rawpos);
+ PosMm.push_back(CC_calpos);
+
+ RawTof.push_back(rawtof);
+ CalTof.push_back(caltof);
+ VelocityMNs.push_back(velocity_mns);
+ Beta.push_back(betaS2);
+ }
+ }
+ }
+ }
+ }
+ }
+ else if(m_NumberOfDetectors==1){
+ double CC_rawpos;
+ double CC_rawtime;
+ if(CC_pmtTref.size()==1){
+ for(unsigned int i=0; i<CC_pmtR.size(); i++){
+ for(unsigned int j=0; j<CC_pmtL.size(); j++){
+ CC_rawpos = CC_pmtR[i] - CC_pmtL[j];
+ CC_rawtime = 0.5*(CC_pmtR[i] + CC_pmtL[j]);
+
+ if(CC_rawpos<m_DET2_PosNs_Min || CC_rawpos>m_DET2_PosNs_Max)
+ continue;
+
+ //Filling output Tree;
+ DetectorNbr.push_back(1);
+ TimeNs.push_back(CC_rawtime);
+ PosNs.push_back(CC_rawpos);
+ }
+ }
+ }
+ }
+}
+
+///////////////////////////////////////////////////////////////////////////
+double TSofSciPhysics::CalculateTimeNs(int det, int pmt, int ft, int ct){
+
+ static CalibrationManager* Cal = CalibrationManager::getInstance();
+ //ft = ft+1;
+ 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 ift_ns;
+
+ if(r<0){
+ double par_prev = Cal->GetValue("SofSci/DET"+NPL::itoa(det)+"_SIGNAL"+NPL::itoa(pmt)+"_TIME",ft-1);
+ ift_ns = par + r*(par - par_prev);
+ }
+
+ else{
+ double par_next = Cal->GetValue("SofSci/DET"+NPL::itoa(det)+"_SIGNAL"+NPL::itoa(pmt)+"_TIME",ft+1);
+ ift_ns = par + r*(par_next - par);
+ }
+
+ double time_ns = (double)ict_ns - ift_ns;
+
+ return time_ns;
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TSofSciPhysics::PreTreat() {
+ // This method typically applies thresholds and calibrations
+ // Might test for disabled channels for more complex detector
+
+ // clear pre-treated object
+ ClearPreTreatedData();
+
+ // instantiate CalibrationManager
+ static CalibrationManager* Cal = CalibrationManager::getInstance();
+
+ unsigned int mysize = m_EventData->GetMultiplicity();
+ for (unsigned int i = 0; i < mysize ; ++i) {
+
+ m_PreTreatedData->SetDetectorNbr(m_EventData->GetDetectorNbr(i));
+ m_PreTreatedData->SetPmt(m_EventData->GetPmt(i));
+ m_PreTreatedData->SetCoarseTime(m_EventData->GetCoarseTime(i));
+ m_PreTreatedData->SetFineTime(m_EventData->GetFineTime(i));
+ }
+}
+
+
+
+///////////////////////////////////////////////////////////////////////////
+void TSofSciPhysics::ReadAnalysisConfig() {
+ bool ReadingStatus = false;
+
+ // path to file
+ string FileName = "./configs/ConfigSofSci.dat";
+
+ // open analysis config file
+ ifstream AnalysisConfigFile;
+ AnalysisConfigFile.open(FileName.c_str());
+
+ if (!AnalysisConfigFile.is_open()) {
+ cout << " No ConfigSofSci.dat found: Default parameter loaded for Analayis " << FileName << endl;
+ return;
+ }
+ cout << " Loading user parameter for Analysis from ConfigSofSci.dat " << endl;
+
+ // Save it in a TAsciiFile
+ TAsciiFile* asciiConfig = RootOutput::getInstance()->GetAsciiFileAnalysisConfig();
+ asciiConfig->AppendLine("%%% ConfigSofSci.dat %%%");
+ asciiConfig->Append(FileName.c_str());
+ asciiConfig->AppendLine("");
+ // read analysis config file
+ string LineBuffer,DataBuffer,whatToDo;
+ while (!AnalysisConfigFile.eof()) {
+ // Pick-up next line
+ getline(AnalysisConfigFile, LineBuffer);
+
+ // search for "header"
+ string name = "ConfigSofSci";
+ if (LineBuffer.compare(0, name.length(), name) == 0)
+ ReadingStatus = true;
+
+ // loop on tokens and data
+ while (ReadingStatus ) {
+ whatToDo="";
+ AnalysisConfigFile >> whatToDo;
+
+ // Search for comment symbol (%)
+ if (whatToDo.compare(0, 1, "%") == 0) {
+ AnalysisConfigFile.ignore(numeric_limits<streamsize>::max(), '\n' );
+ }
+
+ else if (whatToDo=="E_THRESHOLD") {
+ AnalysisConfigFile >> DataBuffer;
+ m_E_Threshold = atof(DataBuffer.c_str());
+ cout << whatToDo << " " << m_E_Threshold << endl;
+ }
+ else if (whatToDo=="DET1_POSNS_MIN") {
+ AnalysisConfigFile >> DataBuffer;
+ m_DET1_PosNs_Min = atof(DataBuffer.c_str());
+ cout << whatToDo << " " << m_DET1_PosNs_Min << endl;
+ }
+ else if (whatToDo=="DET1_POSNS_MAX") {
+ AnalysisConfigFile >> DataBuffer;
+ m_DET1_PosNs_Max = atof(DataBuffer.c_str());
+ cout << whatToDo << " " << m_DET1_PosNs_Max << endl;
+ }
+ else if (whatToDo=="DET2_POSNS_MIN") {
+ AnalysisConfigFile >> DataBuffer;
+ m_DET2_PosNs_Min = atof(DataBuffer.c_str());
+ cout << whatToDo << " " << m_DET2_PosNs_Min << endl;
+ }
+ else if (whatToDo=="DET2_POSNS_MAX") {
+ AnalysisConfigFile >> DataBuffer;
+ m_DET2_PosNs_Max = atof(DataBuffer.c_str());
+ cout << whatToDo << " " << m_DET2_PosNs_Max << endl;
+ }
+ else if (whatToDo=="RAWTOF_MIN") {
+ AnalysisConfigFile >> DataBuffer;
+ m_RawTof_Min = atof(DataBuffer.c_str());
+ cout << whatToDo << " " << m_RawTof_Min << endl;
+ }
+ else if (whatToDo=="RAWTOF_MAX") {
+ AnalysisConfigFile >> DataBuffer;
+ m_RawTof_Max = atof(DataBuffer.c_str());
+ cout << whatToDo << " " << m_RawTof_Max << endl;
+ }
+
+
+
+ else {
+ ReadingStatus = false;
+ }
+ }
+ }
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TSofSciPhysics::Clear() {
+ DetectorNbr.clear();
+ TimeNs.clear();
+ PosNs.clear();
+ PosMm.clear();
+ RawTof.clear();
+ CalTof.clear();
+ VelocityMNs.clear();
+ Beta.clear();
+}
+
+
+
+///////////////////////////////////////////////////////////////////////////
+void TSofSciPhysics::ReadConfiguration(NPL::InputParser parser) {
+ vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("SofSci");
+ if(NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << "//// " << blocks.size() << " detectors found " << endl;
+
+ vector<string> cart = {"POS"};
+ vector<string> sphe = {"R","Theta","Phi"};
+
+ for(unsigned int i = 0 ; i < blocks.size() ; i++){
+ if(blocks[i]->HasTokenList(cart)){
+ if(NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << endl << "//// SofSci " << i+1 << endl;
+
+ TVector3 Pos = blocks[i]->GetTVector3("POS","mm");
+ m_X = Pos.Z();
+ m_Y = Pos.Z();
+ m_Z = Pos.Z();
+ AddDetector(Pos);
+ }
+ else if(blocks[i]->HasTokenList(sphe)){
+ if(NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << endl << "//// SofSci " << i+1 << endl;
+ double R = blocks[i]->GetDouble("R","mm");
+ double Theta = blocks[i]->GetDouble("Theta","deg");
+ double Phi = blocks[i]->GetDouble("Phi","deg");
+ AddDetector(R,Theta,Phi);
+ }
+ else{
+ cout << "ERROR: check your input file formatting " << endl;
+ exit(1);
+ }
+ }
+
+ ReadAnalysisConfig();
+}
+
+
+///////////////////////////////////////////////////////////////////////////
+void TSofSciPhysics::AddParameterToCalibrationManager() {
+ CalibrationManager* Cal = CalibrationManager::getInstance();
+
+ Cal->AddParameter("SofSci","TOF2INV_V","SofSci_TOF2INV_V");
+ Cal->AddParameter("SofSci","LENGTH_S2","SofSci_LENGTH_S2");
+
+ for(int d = 0; d < m_NumberOfDetectors; d++){
+ Cal->AddParameter("SofSci","DET"+NPL::itoa(d+1)+"_POSPAR","SofSci_DET"+NPL::itoa(d+1)+"_POSPAR");
+ for(int s = 0; s < m_NumberOfSignals; s++){
+ Cal->AddParameter("SofSci","DET"+NPL::itoa(d+1)+"_SIGNAL"+NPL::itoa(s+1)+"_TIME","SofSci_DET"+NPL::itoa(d+1)+"_SIGNAL"+NPL::itoa(s+1)+"_TIME");
+ Cal->AddParameter("SofSci","DET"+NPL::itoa(d+1)+"_SIGNAL"+NPL::itoa(s+1)+"_CLOCKOFFSET","SofSci_DET"+NPL::itoa(d+1)+"_SIGNAL"+NPL::itoa(s+1)+"_CLOCKOFFSET");
+ }
+ }
+}
+
+
+
+///////////////////////////////////////////////////////////////////////////
+void TSofSciPhysics::InitializeRootInputRaw() {
+ TChain* inputChain = RootInput::getInstance()->GetChain();
+ inputChain->SetBranchStatus("SofSci", true );
+ inputChain->SetBranchAddress("SofSci", &m_EventData );
+}
+
+
+
+///////////////////////////////////////////////////////////////////////////
+void TSofSciPhysics::InitializeRootInputPhysics() {
+ TChain* inputChain = RootInput::getInstance()->GetChain();
+ inputChain->SetBranchAddress("SofSci", &m_EventPhysics);
+}
+
+
+
+///////////////////////////////////////////////////////////////////////////
+void TSofSciPhysics::InitializeRootOutput() {
+ TTree* outputTree = RootOutput::getInstance()->GetTree();
+ outputTree->Branch("SofSci", "TSofSciPhysics", &m_EventPhysics);
+}
+
+
+
+////////////////////////////////////////////////////////////////////////////////
+// Construct Method to be pass to the DetectorFactory //
+////////////////////////////////////////////////////////////////////////////////
+NPL::VDetector* TSofSciPhysics::Construct() {
+ return (NPL::VDetector*) new TSofSciPhysics();
+}
+
+
+
+////////////////////////////////////////////////////////////////////////////////
+// Registering the construct method to the factory //
+////////////////////////////////////////////////////////////////////////////////
+extern "C"{
+ class proxy_SofSci{
+ public:
+ proxy_SofSci(){
+ NPL::DetectorFactory::getInstance()->AddToken("SofSci","SofSci");
+ NPL::DetectorFactory::getInstance()->AddDetector("SofSci",TSofSciPhysics::Construct);
+ }
+ };
+
+ proxy_SofSci p_SofSci;
+}
+
diff --git a/NPLib/Detectors/Sofia/TSofSciPhysics.h b/NPLib/Detectors/Sofia/TSofSciPhysics.h
new file mode 100644
index 0000000000000000000000000000000000000000..fe74eca0f1b4df613db81ca54ee6a7908d568b99
--- /dev/null
+++ b/NPLib/Detectors/Sofia/TSofSciPhysics.h
@@ -0,0 +1,173 @@
+#ifndef TSofSciPHYSICS_H
+#define TSofSciPHYSICS_H
+/*****************************************************************************
+ * Copyright (C) 2009-2020 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: Pierre Morfouace contact address: pierre.morfouace2@cea.fr *
+ * *
+ * Creation Date : November 2020 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class hold SofSci Treated data *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+
+// C++ headers
+#include <vector>
+#include <map>
+#include <string>
+using namespace std;
+
+// ROOT headers
+#include "TObject.h"
+#include "TH1.h"
+#include "TVector3.h"
+#include "TSpline.h"
+#include "TRandom3.h"
+// NPTool headers
+#include "TSofSciData.h"
+#include "NPCalibrationManager.h"
+#include "NPVDetector.h"
+#include "NPInputParser.h"
+
+
+
+class TSofSciPhysics : public TObject, public NPL::VDetector {
+ //////////////////////////////////////////////////////////////
+ // constructor and destructor
+ public:
+ TSofSciPhysics();
+ ~TSofSciPhysics() {};
+
+
+ //////////////////////////////////////////////////////////////
+ // Inherited from TObject and overriden to avoid warnings
+ public:
+ void Clear();
+ void Clear(const Option_t*) {};
+
+
+ //////////////////////////////////////////////////////////////
+ // data obtained after BuildPhysicalEvent() and stored in
+ // output ROOT file
+ public:
+ vector<int> DetectorNbr;
+ vector<double> TimeNs;
+ vector<double> PosNs;
+ vector<double> PosMm;
+ vector<double> RawTof;
+ vector<double> CalTof;
+ vector<double> VelocityMNs;
+ vector<double> Beta;
+
+ /// A usefull method to bundle all operation to add a detector
+ void AddDetector(TVector3 POS);
+ void AddDetector(double R, double Theta, double Phi);
+
+ //////////////////////////////////////////////////////////////
+ // methods inherited from the VDetector ABC class
+ public:
+ // read stream from ConfigFile to pick-up detector parameters
+ void ReadConfiguration(NPL::InputParser);
+
+ // add parameters to the CalibrationManger
+ void AddParameterToCalibrationManager();
+
+ // method called event by event, aiming at extracting the
+ // physical information from detector
+ void BuildPhysicalEvent();
+
+ // same as BuildPhysicalEvent() method but with a simpler
+ // treatment
+ void BuildSimplePhysicalEvent();
+
+ // same as above but for online analysis
+ void BuildOnlinePhysicalEvent() {BuildPhysicalEvent();};
+
+ // activate raw data object and branches from input TChain
+ // in this method mother branches (Detector) AND daughter leaves
+ // (fDetector_parameter) have to be activated
+ void InitializeRootInputRaw();
+
+ // activate physics data object and branches from input TChain
+ // in this method mother branches (Detector) AND daughter leaves
+ // (fDetector_parameter) have to be activated
+ void InitializeRootInputPhysics();
+
+ // create branches of output ROOT file
+ void InitializeRootOutput();
+
+ // clear the raw and physical data objects event by event
+ void ClearEventPhysics() {Clear();}
+ void ClearEventData() {m_EventData->Clear();}
+
+ double GetDistance(){
+ return sqrt(m_X*m_X + m_Y*m_Y + m_Z*m_Z);
+ }
+ double CalculateTimeNs(int, int, int, int);
+ double GetNumberOfDetectors() {return m_NumberOfDetectors;}
+
+ //////////////////////////////////////////////////////////////
+ // specific methods to SofSci array
+ public:
+ // remove bad channels, calibrate the data and apply thresholds
+ void PreTreat();
+
+ // clear the pre-treated object
+ void ClearPreTreatedData() {m_PreTreatedData->Clear();}
+
+ // read the user configuration file. If no file is found, load standard one
+ void ReadAnalysisConfig();
+
+ // give and external TSofSciData object to TSofSciPhysics.
+ // needed for online analysis for example
+ void SetRawDataPointer(TSofSciData* rawDataPointer) {m_EventData = rawDataPointer;}
+
+ // objects are not written in the TTree
+ private:
+ TSofSciData* m_EventData; //!
+ TSofSciData* m_PreTreatedData; //!
+ TSofSciPhysics* m_EventPhysics; //!
+
+ // getters for raw and pre-treated data object
+ public:
+ TSofSciData* GetRawData() const {return m_EventData;}
+ TSofSciData* GetPreTreatedData() const {return m_PreTreatedData;}
+
+ // parameters used in the analysis
+ private:
+ double m_E_Threshold; //!
+ double m_DET1_PosNs_Min; //!
+ double m_DET1_PosNs_Max; //!
+ double m_DET2_PosNs_Min; //!
+ double m_DET2_PosNs_Max; //!
+ double m_RawTof_Min; //!
+ double m_RawTof_Max; //!
+ double m_X; //!
+ double m_Y; //!
+ double m_Z; //!
+ TRandom3 rand; //!
+
+ // number of detectors
+ private:
+ int m_NumberOfDetectors; //!
+ int m_NumberOfSignals; //!
+ int m_NumberOfPmts; //!
+
+ // Static constructor to be passed to the Detector Factory
+ public:
+ static NPL::VDetector* Construct();
+
+ ClassDef(TSofSciPhysics,1) // SofSciPhysics structure
+};
+#endif
diff --git a/NPLib/Detectors/Sofia/TSofTrimPhysics.cxx b/NPLib/Detectors/Sofia/TSofTrimPhysics.cxx
index 5d55976570e598be69d8c18f5c699af326b02439..3f3ecea24dc8d681fa056e05ffa5f280a8953175 100644
--- a/NPLib/Detectors/Sofia/TSofTrimPhysics.cxx
+++ b/NPLib/Detectors/Sofia/TSofTrimPhysics.cxx
@@ -48,7 +48,10 @@ TSofTrimPhysics::TSofTrimPhysics()
m_PreTreatedData(new TSofTrimData),
m_EventPhysics(this),
m_NumberOfDetectors(0),
+ m_Beta(0.838),
+ m_BetaNorm(0.838),
m_NumberOfSections(3),
+ m_NumberOfAnodesPaired(3),
m_NumberOfAnodesPerSection(6) {
}
@@ -80,12 +83,119 @@ void TSofTrimPhysics::BuildSimplePhysicalEvent() {
void TSofTrimPhysics::BuildPhysicalEvent() {
// apply thresholds and calibration
PreTreat();
+ if(m_PreTreatedData->GetMultiplicity() != 18)
+ return;
+
+ double Ep1[3], DTp1[3];
+ double Ep2[3], DTp2[3];
+ double Ep3[3], DTp3[3];
+ double Esec[3];
+ for(int i=0; i<m_NumberOfSections; i++){
+ Ep1[i] = 0;
+ Ep2[i] = 0;
+ Ep3[i] = 0;
+ DTp1[i] = 0;
+ DTp2[i] = 0;
+ DTp3[i] = 0;
+ Esec[i] = 0;
+ }
- // match energy and time together
unsigned int mysizeE = m_PreTreatedData->GetMultiplicity();
for (UShort_t e = 0; e < mysizeE ; e++) {
- //to do
+ int SectionNbr = m_PreTreatedData->GetSectionNbr(e);
+ int AnodeNbr = m_PreTreatedData->GetAnodeNbr(e);
+ double Energy = m_PreTreatedData->GetEnergy(e);
+ double DT = m_PreTreatedData->GetDriftTime(e);
+
+ if(AnodeNbr==1 || AnodeNbr==2){
+ Ep1[SectionNbr-1] += Energy;
+ DTp1[SectionNbr-1] += DT;
+ }
+ if(AnodeNbr==3 || AnodeNbr==4){
+ Ep2[SectionNbr-1] += Energy;
+ DTp2[SectionNbr-1] += DT;
+ }
+ if(AnodeNbr==5 || AnodeNbr==6){
+ Ep3[SectionNbr-1] += Energy;
+ DTp3[SectionNbr-1] += DT;
+ }
+ }
+
+
+ for(int i=0; i<m_NumberOfSections; i++){
+ DTp1[i] = 0.5*DTp1[i];
+ DTp2[i] = 0.5*DTp2[i];
+ DTp3[i] = 0.5*DTp3[i];
+
+ Ep1[i] = 0.5*Ep1[i];
+ Ep2[i] = 0.5*Ep2[i];
+ Ep3[i] = 0.5*Ep3[i];
+ }
+
+ static CalibrationManager* Cal = CalibrationManager::getInstance();
+ double Ddt = DTp2[2] - DTp2[0];
+ double p0_1[3], p0_2[3], p0_3[3];
+ double p1_1[3], p1_2[3], p1_3[3];
+
+ for(int i=0; i<m_NumberOfSections; i++){
+ // Energy Alignement of pairs per section
+ Ep1[i] = Cal->ApplyCalibration("SofTrim/SEC"+NPL::itoa(i+1)+"_ANODE1_ALIGN",Ep1[i]);
+ Ep2[i] = Cal->ApplyCalibration("SofTrim/SEC"+NPL::itoa(i+1)+"_ANODE2_ALIGN",Ep2[i]);
+ Ep3[i] = Cal->ApplyCalibration("SofTrim/SEC"+NPL::itoa(i+1)+"_ANODE3_ALIGN",Ep3[i]);
+
+ // Beta correction per pair: DE = [0] + [1]*pow(Beta, -5./3);
+ p0_1[i] = Cal->GetValue("SofTrim/SEC"+NPL::itoa(i+1)+"_ANODE1_BETA",0);
+ p0_2[i] = Cal->GetValue("SofTrim/SEC"+NPL::itoa(i+1)+"_ANODE2_BETA",0);
+ p0_3[i] = Cal->GetValue("SofTrim/SEC"+NPL::itoa(i+1)+"_ANODE3_BETA",0);
+ p1_1[i] = Cal->GetValue("SofTrim/SEC"+NPL::itoa(i+1)+"_ANODE1_BETA",1);
+ 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.));
+
+ // Angle correction per pair: spline
+ Ep1[i] = Ep1[i] / fcorr_EvsA[i][0]->Eval(Ddt) * fcorr_EvsA[i][0]->Eval(0);
+ Ep2[i] = Ep2[i] / fcorr_EvsA[i][1]->Eval(Ddt) * fcorr_EvsA[i][1]->Eval(0);
+ Ep3[i] = Ep3[i] / fcorr_EvsA[i][2]->Eval(Ddt) * fcorr_EvsA[i][2]->Eval(0);
+
+ // DT correction per pair: spline
+ Ep1[i] = Ep1[i] / fcorr_EvsDT[i][0]->Eval(DTp1[i]) * fcorr_EvsDT[i][0]->Eval(3000);
+ Ep2[i] = Ep2[i] / fcorr_EvsDT[i][1]->Eval(DTp2[i]) * fcorr_EvsDT[i][1]->Eval(3000);
+ Ep3[i] = Ep3[i] / fcorr_EvsDT[i][2]->Eval(DTp3[i]) * fcorr_EvsDT[i][2]->Eval(3000);
+ }
+
+ for(int i=0; i<m_NumberOfSections; i++){
+ Esec[i] = (Ep1[i] + Ep2[i] + Ep3[i])/3;
+
+ // 2nd DT correction per section: spline
+ Esec[i] = Esec[i] / fcorr_sec[i]->Eval(DTp2[i]) * fcorr_sec[i]->Eval(0);
+
+ // Section ALignement
+ Esec[i] = Cal->ApplyCalibration("SofTrim/SEC"+NPL::itoa(i+1)+"_ALIGN",Esec[i]);
}
+
+
+
+ // Filling output Tree //
+ for(int i=0; i<m_NumberOfSections; i++){
+ if(DTp2[i]!=0){
+ SectionNbr.push_back(i+1);
+ EnergyPair1.push_back(Ep1[i]);
+ EnergyPair2.push_back(Ep2[i]);
+ EnergyPair3.push_back(Ep3[i]);
+ DriftTimePair1.push_back(DTp1[i]);
+ DriftTimePair2.push_back(DTp2[i]);
+ DriftTimePair3.push_back(DTp3[i]);
+ EnergySection.push_back(Esec[i]);
+ Theta.push_back(DTp2[2]-DTp2[0]);
+ }
+ }
+
}
///////////////////////////////////////////////////////////////////////////
@@ -101,13 +211,13 @@ void TSofTrimPhysics::PreTreat() {
unsigned int mysize = m_EventData->GetMultiplicity();
for (unsigned int i = 0; i < mysize ; ++i) {
- Double_t Energy = Cal->ApplyCalibration("SofTrim/ENERGY_SEC"+NPL::itoa(m_EventData->GetSectionNbr(i))+"_ANODE"+NPL::itoa(m_EventData->GetAnodeNbr(i))+"_ENERGY",m_EventData->GetEnergy(i));
- Double_t Time = Cal->ApplyCalibration("SofTrim/TIME_SEC"+NPL::itoa(m_EventData->GetSectionNbr(i))+"_ANODE"+NPL::itoa(m_EventData->GetAnodeNbr(i))+"_TIME",m_EventData->GetDriftTime(i));
-
+ Double_t Energy = Cal->ApplyCalibration("SofTrim/SEC"+NPL::itoa(m_EventData->GetSectionNbr(i))+"_ANODE"+NPL::itoa(m_EventData->GetAnodeNbr(i))+"_ENERGY",m_EventData->GetEnergy(i));
+ Double_t DT = Cal->ApplyCalibration("SofTrim/SEC"+NPL::itoa(m_EventData->GetSectionNbr(i))+"_ANODE"+NPL::itoa(m_EventData->GetAnodeNbr(i))+"_TIME",m_EventData->GetDriftTime(i));
+
m_PreTreatedData->SetSectionNbr(m_EventData->GetSectionNbr(i));
m_PreTreatedData->SetAnodeNbr(m_EventData->GetAnodeNbr(i));
m_PreTreatedData->SetEnergy(Energy);
- m_PreTreatedData->SetDriftTime(Time);
+ m_PreTreatedData->SetDriftTime(DT);
m_PreTreatedData->SetPileUp(m_EventData->GetPileUp(i));
m_PreTreatedData->SetOverflow(m_EventData->GetOverflow(i));
}
@@ -158,17 +268,33 @@ void TSofTrimPhysics::ReadAnalysisConfig() {
AnalysisConfigFile.ignore(numeric_limits<streamsize>::max(), '\n' );
}
- /*else if (whatToDo=="E_RAW_THRESHOLD") {
+ else if (whatToDo=="SPLINE_PAIR_ANGLE_PATH") {
+ AnalysisConfigFile >> DataBuffer;
+ m_SPLINE_PAIR_ANGLE_PATH = DataBuffer;
+ cout << "*** Loading Spline for Angle correction per pair ***" << endl;
+ LoadSplinePairAngle();
+ }
+
+ else if (whatToDo=="SPLINE_PAIR_DT_PATH") {
+ AnalysisConfigFile >> DataBuffer;
+ m_SPLINE_PAIR_DT_PATH = DataBuffer;
+ cout << "*** Loading Spline for Dritf Time correction per pair ***" << endl;
+ LoadSplinePairDriftTime();
+ }
+
+ else if (whatToDo=="SPLINE_SECTION_DT_PATH") {
AnalysisConfigFile >> DataBuffer;
- m_E_RAW_Threshold = atof(DataBuffer.c_str());
- cout << whatToDo << " " << m_E_RAW_Threshold << endl;
- }
+ m_SPLINE_SECTION_DT_PATH = DataBuffer;
+ cout << "*** Loading Spline for Dritf Time correction per section ***" << endl;
+ LoadSplineSectionDriftTime();
+ }
- else if (whatToDo=="E_THRESHOLD") {
+
+ else if (whatToDo=="E_THRESHOLD") {
AnalysisConfigFile >> DataBuffer;
m_E_Threshold = atof(DataBuffer.c_str());
cout << whatToDo << " " << m_E_Threshold << endl;
- }*/
+ }
else {
ReadingStatus = false;
@@ -177,6 +303,77 @@ void TSofTrimPhysics::ReadAnalysisConfig() {
}
}
+///////////////////////////////////////////////////////////////////////////
+double TSofTrimPhysics::GetMaxEnergySection(){
+ double Emax=-1;
+
+ if(EnergySection.size() != 3){
+ cout << "WARNING! Size of EnergySection different than 3, size= " << EnergySection.size() << endl;
+ return Emax;
+ }
+
+ Emax = *max_element(EnergySection.begin(), EnergySection.end());
+
+ return Emax;
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TSofTrimPhysics::LoadSplinePairAngle(){
+ TString filename = m_SPLINE_PAIR_ANGLE_PATH;
+ TFile* ifile = new TFile(filename,"read");
+
+ if(ifile->IsOpen()){
+ cout << "Loading splines..." << endl;
+ for(int s=0; s<m_NumberOfSections; s++){
+ for(int a=0; a<3; a++){
+ TString splinename = Form("spline_EvsA_sec%i_anode%i",s+1,a+1);
+ fcorr_EvsA[s][a] = (TSpline3*) ifile->FindObjectAny(splinename);
+ cout << fcorr_EvsA[s][a]->GetName() << endl;
+ }
+ }
+ }
+ else
+ cout << "File " << filename << " not found!" << endl;
+ ifile->Close();
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TSofTrimPhysics::LoadSplinePairDriftTime(){
+ TString filename = m_SPLINE_PAIR_DT_PATH;
+ TFile* ifile = new TFile(filename,"read");
+
+ if(ifile->IsOpen()){
+ cout << "Loading splines..." << endl;
+ for(int s=0; s<m_NumberOfSections; s++){
+ for(int a=0; a<3; a++){
+ TString splinename = Form("spline_EvsDT_sec%i_anode%i",s+1,a+1);
+ fcorr_EvsDT[s][a] = (TSpline3*) ifile->FindObjectAny(splinename);
+ cout << fcorr_EvsDT[s][a]->GetName() << endl;
+ }
+ }
+ }
+ else
+ cout << "File " << filename << " not found!" << endl;
+ ifile->Close();
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TSofTrimPhysics::LoadSplineSectionDriftTime(){
+ TString filename = m_SPLINE_SECTION_DT_PATH;
+ TFile* ifile = new TFile(filename,"read");
+
+ if(ifile->IsOpen()){
+ cout << "Loading splines..." << endl;
+ for(int s=0; s<m_NumberOfSections; s++){
+ TString splinename = Form("spline_sec%i",s+1);
+ fcorr_sec[s] = (TSpline3*) ifile->FindObjectAny(splinename);
+ cout << fcorr_sec[s]->GetName() << endl;
+ }
+ }
+ else
+ cout << "File " << filename << " not found!" << endl;
+ ifile->Close();
+}
///////////////////////////////////////////////////////////////////////////
@@ -184,11 +381,12 @@ void TSofTrimPhysics::Clear() {
SectionNbr.clear();
EnergyPair1.clear();
EnergyPair2.clear();
- EnergyPair2.clear();
+ EnergyPair3.clear();
DriftTimePair1.clear();
DriftTimePair2.clear();
DriftTimePair3.clear();
- EnergySum.clear();
+ EnergySection.clear();
+ Theta.clear();
}
@@ -223,12 +421,15 @@ void TSofTrimPhysics::ReadConfiguration(NPL::InputParser parser) {
exit(1);
}
}
+
+ ReadAnalysisConfig();
}
///////////////////////////////////////////////////////////////////////////
void TSofTrimPhysics::AddParameterToCalibrationManager() {
CalibrationManager* Cal = CalibrationManager::getInstance();
+
for(int sec = 0; sec < m_NumberOfSections; sec++){
for(int anode = 0; anode < m_NumberOfAnodesPerSection; anode++){
Cal->AddParameter("SofTrim","SEC"+NPL::itoa(sec+1)+"_ANODE"+NPL::itoa(anode+1)+"_ENERGY","SofTrim_SEC"+NPL::itoa(sec+1)+"_ANODE"+NPL::itoa(anode+1)+"_ENERGY");
@@ -236,9 +437,16 @@ void TSofTrimPhysics::AddParameterToCalibrationManager() {
}
}
-}
-
+ for(int sec = 0; sec < m_NumberOfSections; sec++){
+ Cal->AddParameter("SofTrim","SEC"+NPL::itoa(sec+1)+"_ALIGN","SofTrim_SEC"+NPL::itoa(sec+1)+"_ALIGN");
+
+ for(int anode = 0; anode < m_NumberOfAnodesPaired; anode++){
+ Cal->AddParameter("SofTrim","SEC"+NPL::itoa(sec+1)+"_ANODE"+NPL::itoa(anode+1)+"_ALIGN","SofTrim_SEC"+NPL::itoa(sec+1)+"_ANODE"+NPL::itoa(anode+1)+"_ALIGN");
+ Cal->AddParameter("SofTrim","SEC"+NPL::itoa(sec+1)+"_ANODE"+NPL::itoa(anode+1)+"_BETA","SofTrim_SEC"+NPL::itoa(sec+1)+"_ANODE"+NPL::itoa(anode+1)+"_BETA");
+ }
+ }
+}
///////////////////////////////////////////////////////////////////////////
void TSofTrimPhysics::InitializeRootInputRaw() {
diff --git a/NPLib/Detectors/Sofia/TSofTrimPhysics.h b/NPLib/Detectors/Sofia/TSofTrimPhysics.h
index 759fff86d56b500de3f6766655701e592e2da0ce..923b46288ada001d9abb38d94f3c266a76eeaade 100644
--- a/NPLib/Detectors/Sofia/TSofTrimPhysics.h
+++ b/NPLib/Detectors/Sofia/TSofTrimPhysics.h
@@ -14,7 +14,7 @@
* Last update : *
*---------------------------------------------------------------------------*
* Decription: *
- * This class hold TofTofW Treated data *
+ * This class hold SofTrim Treated data *
* *
*---------------------------------------------------------------------------*
* Comment: *
@@ -32,6 +32,7 @@ using namespace std;
#include "TObject.h"
#include "TH1.h"
#include "TVector3.h"
+#include "TSpline.h"
// NPTool headers
#include "TSofTrimData.h"
#include "NPCalibrationManager.h"
@@ -66,7 +67,8 @@ class TSofTrimPhysics : public TObject, public NPL::VDetector {
vector<double> DriftTimePair1;
vector<double> DriftTimePair2;
vector<double> DriftTimePair3;
- vector<double> EnergySum;
+ vector<double> EnergySection;
+ vector<double> Theta;
/// A usefull method to bundle all operation to add a detector
void AddDetector(TVector3 POS);
@@ -125,7 +127,16 @@ class TSofTrimPhysics : public TObject, public NPL::VDetector {
// give and external TSofTrimData object to TSofTrimPhysics.
// needed for online analysis for example
void SetRawDataPointer(TSofTrimData* rawDataPointer) {m_EventData = rawDataPointer;}
-
+
+ void LoadSplinePairAngle();
+ void LoadSplinePairDriftTime();
+ void LoadSplineSectionDriftTime();
+
+ void SetBeta(double beta) {m_Beta = beta;}
+ double GetBeta() {return m_Beta;}
+
+ double GetMaxEnergySection();
+
// objects are not written in the TTree
private:
TSofTrimData* m_EventData; //!
@@ -139,13 +150,23 @@ class TSofTrimPhysics : public TObject, public NPL::VDetector {
// parameters used in the analysis
private:
- // thresholds
double m_E_Threshold; //!
+ double m_Beta; //!
+ double m_BetaNorm; //!
+ string m_SPLINE_PAIR_ANGLE_PATH; //!
+ string m_SPLINE_PAIR_DT_PATH; //!
+ string m_SPLINE_SECTION_DT_PATH; //!
+
+ private:
+ TSpline3* fcorr_EvsA[3][3]; //!
+ TSpline3* fcorr_EvsDT[3][3]; //!
+ TSpline3* fcorr_sec[3]; //!
// number of detectors
private:
int m_NumberOfDetectors; //!
int m_NumberOfSections; //!
+ int m_NumberOfAnodesPaired; //!
int m_NumberOfAnodesPerSection; //!
// Static constructor to be passed to the Detector Factory
diff --git a/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.cxx b/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.cxx
index 26d119ac77efcb47c6020b7254afe83576c1d22d..a551e0c5f3adb468aa3fa51613a25f6eb05f5d9b 100755
--- a/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.cxx
+++ b/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.cxx
@@ -14,6 +14,8 @@ double devE(const double* parameter){
TVector3 dir;
//Initilize histogram
+ //canv->Clear();
+ //canv->ResetDrawn();
h->Reset();
h1->Reset();
h2->Reset();
@@ -22,36 +24,15 @@ double devE(const double* parameter){
h5->Reset();
h7->Reset();
+ double FitResultMatrix[7][5];
+
//Loop over events
for(unsigned int i = 0 ; i < size ; i++){
//Particle path vector
dir=*(pos[i])-offset;
//Define normal vector for the MG# of detection
- //[[[[ UPDATE WITH NEW MG POSITIONS FROM SURVEY ]]]]
- switch(detnum[i]){
- case 1:
- MugastNormal.SetXYZ(-0.453915, +0.455463, -0.765842);
- break;
- case 2:
- MugastNormal.SetXYZ(-0.642828, +0.000000, -0.766010);
- break;
- case 3:
- MugastNormal.SetXYZ(-0.454594, -0.450670, -0.768271);
- break;
- case 4:
- MugastNormal.SetXYZ(-0.002437, -0.638751, -0.769409);
- break;
- case 5:
- MugastNormal.SetXYZ(+0.452429, -0.454575, -0.767248);
- break;
- case 7:
- MugastNormal.SetXYZ(+0.443072, +0.443265, -0.779232);
- break;
- default:
- cout << "ERROR:: Invalid DetNum " << detnum[i] << " at event " << i << endl;
- return 1; // Exit code
- }
+ DetectorSwitch(detnum[i], MugastNormal);
//Angle leaving target, angle entering MUGAST & energy deposited in MUGAST
double ThetaTarget = dir.Angle(TVector3(0,0,1));
@@ -74,51 +55,26 @@ double devE(const double* parameter){
//Fill histograms with Ex
h->Fill(Ex);
- switch(detnum[i]){
- case 1:
- h1->Fill(Ex);
- break;
- case 2:
- h2->Fill(Ex);
- break;
- case 3:
- h3->Fill(Ex);
- break;
- case 4:
- h4->Fill(Ex);
- break;
- case 5:
- h5->Fill(Ex);
- break;
- case 7:
- h7->Fill(Ex);
- break;
- default:
- cout << "ERROR:: Invalid DetNum " << detnum[i] << " at event " << i << endl;
- return 1; // Exit code
- }
+ DetectorSwitch(detnum[i], Ex);
}
- //End loop over events
-
+
+ //Initilise, Draw & Fit histograms
+ InitiliseCanvas(FitResultMatrix);
+
//Write vals to screen
- cout << "Mean: " << h->GetMean()
- << "\t StdDev: " << h->GetStdDev()
- << "\t Thickness: " << parameter[3] << " um"
- << endl;
-
- //Draw histogram(s)
- h->Draw();
- if(flagDraw){ InitiliseCanvas(); }
-
- /*
- cout << pow(h->GetMean()-refE,2) << " + "
- << pow(0.1*h->GetStdDev(),2) << " = "
- << pow(h->GetMean()-refE,2) + pow(0.1*h->GetStdDev(),2)
+ if(flagDraw){cout << "==================================================" << endl;}
+ cout << "Mean: " << FitResultMatrix[mgSelect][0]
+ << " +/- "
+ << FitResultMatrix[mgSelect][1]
+ << " StdDev: " << FitResultMatrix[mgSelect][2]
+ << " +/- "
+ << FitResultMatrix[mgSelect][3]
+ << " Thick: " << parameter[3] << " um"
+ << " Fit Chi2/NDF = " << FitResultMatrix[mgSelect][4]
<< endl;
- */
//Adapt the metric as needed
- return sqrt( pow(h->GetMean()-refE,2) + pow(0.1*h->GetStdDev(),2) );
+ return sqrt( pow(FitResultMatrix[mgSelect][0]-refE,2) + pow(0.1*FitResultMatrix[mgSelect][2],2) );
}
////////////////////////////////////////////////////////////////////////////////
void MinimizeBeamSpot(){
@@ -126,6 +82,19 @@ void MinimizeBeamSpot(){
//Read data in
LoadFile();
+ //Output formatting
+ cout << fixed << showpoint << setprecision(6) << showpos;
+
+ //Read in
+ cout << "==================================================" << endl;
+ cout << "=--------- SELECT TELESCOPE TO MINIMIZE ---------=" << endl;
+ cout << "= Type MG# of telescope metric to use, or type 0 =" << endl;
+ cout << "= to use the sum of all MG's =" << endl;
+ cout << "==================================================" << endl;
+ cin >> mgSelect;
+ if(mgSelect==7){mgSelect=6;} // Correct the input for MG7
+ cout << "==================================================" << endl;
+
//Start with beam (0,0,0) and 4.76um 0.5mg/c2 target
double parameter[4] = {0.0, 0.0, 0.0, 4.76};
devE(parameter);
@@ -133,9 +102,8 @@ void MinimizeBeamSpot(){
//Function with 4 parameter XYZ and Target thickness
auto func = ROOT::Math::Functor(&devE,4);
- //Minimizer
+ //Initilise minimizer
auto minim = ROOT::Math::Factory::CreateMinimizer("Minuit2","Migrad");
-
minim->SetPrintLevel(0);
minim->SetPrecision(1e-10);
@@ -150,20 +118,31 @@ void MinimizeBeamSpot(){
//Don't draw iterations of minimizer
flagDraw = 0;
-
+ //canv->SetBatch(kTRUE);
+ gROOT->SetBatch(kTRUE);
+
//Shrink it, babeyyy
minim->Minimize();
//Draw minimal value
flagDraw = 1;
+ gROOT->SetBatch(kFALSE);
//Pull values from minimizer
const double* x = minim->X();
- cout << "========================================" << endl;
- cout << "\t\tX =" << x[0] << endl;
- cout << "\t\tY =" << x[1] << endl;
- cout << "\t\tZ =" << x[2] << endl;
- cout << "\t\tT =" << x[3] << endl;
- devE(x);
- cout << "========================================" << endl;
+ cout << "==================================================" << endl;
+ cout << "=---------------- FINAL PEAK FITS ---------------=" << endl;
+ cout << "==================================================" << endl;
+ devE(x);
+// canv->DrawClone();
+ cout << "==================================================" << endl;
+ cout << "=------------ RESULTS OF MINIMIZATION -----------=" << endl;
+ if(mgSelect==6){mgSelect=7;} // Correct the input for MG7
+ cout << "=------------------- USING MG " << mgSelect << " -----------------=" << endl;
+ cout << "==================================================" << endl;
+ cout << "\t\tX = " << x[0] << " mm" << endl;
+ cout << "\t\tY = " << x[1] << " mm" << endl;
+ cout << "\t\tZ = " << x[2] << " mm" << endl;
+ cout << "\t\tT = " << x[3] << " um" << endl;
+ cout << "==================================================" << endl;
}
diff --git a/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.h b/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.h
index 5994406553c00aa43a6cce6ec12b5a62b08ea547..ef2448e5e3c27a6d5330e992fc5da810f0fba6af 100755
--- a/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.h
+++ b/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.h
@@ -9,12 +9,21 @@ double refE = 0.143; // the energy of the selected states
vector<TVector3*> pos;
vector<double> energy;
vector<int> detnum;
+unsigned int mgSelect = 10;
NPL::EnergyLoss CD2("proton_CD2.G4table","G4Table",100);
NPL::EnergyLoss Al("proton_Al.G4table","G4Table",100);
using namespace std;
bool flagDraw = 0;
-static auto h = new TH1D("h","All MG#'s", 60,-1.,1.);
+
+//double FitResultMatrix[7][5];
+// 7 => Sum in 0 and them MG's in 1-6
+// 5 => Mean, MeanErr, StdDev, StdDevErr, Chi2/NDF
+
+
+//TCanvas *canv = new TCanvas("canv","Ex Histograms",20,20,1600,800);
+
+static auto h = new TH1D("h","All MG#'s", 80,-1.,1.);
static auto h1 = new TH1D("h1","Individual MG#'s", 40,-1.,1.);
static auto h2 = new TH1D("h2","h2", 40,-1.,1.);
static auto h3 = new TH1D("h3","h3", 40,-1.,1.);
@@ -46,7 +55,104 @@ void LoadFile(){
file.close();
}
////////////////////////////////////////////////////////////////////////////////
-void InitiliseCanvas(){
+void FillMatrix(double* matrix, TFitResultPtr fit){
+ matrix[0] = fit->Parameter(1); //Mean
+ matrix[1] = fit->ParError(1);
+ matrix[2] = fit->Parameter(2); //StdDev
+ matrix[3] = fit->ParError(2);
+ matrix[4] = fit->Chi2()/fit->Ndf(); //Chi2/NDF
+
+ if(flagDraw){
+ cout << "\n Mean = " << matrix[0] << " +/- " << matrix[1] << endl;
+ cout << " StdDev = " << matrix[2] << " +/- " << matrix[3] << endl;
+ cout << " Chi2/NDF = " << matrix[4] << endl;
+ }
+}
+////////////////////////////////////////////////////////////////////////////////
+//[[[[ UPDATE WITH NEW MG POSITIONS FROM SURVEY ]]]]
+//Overloaded function definiton; this is for MUGAST Normal vectors
+void DetectorSwitch(int MG, TVector3& Normal ){
+ switch(MG){
+ case 1:
+ Normal.SetXYZ(-0.453915, +0.455463, -0.765842);
+ break;
+ case 2:
+ Normal.SetXYZ(-0.642828, +0.000000, -0.766010);
+ break;
+ case 3:
+ Normal.SetXYZ(-0.454594, -0.450670, -0.768271);
+ break;
+ case 4:
+ Normal.SetXYZ(-0.002437, -0.638751, -0.769409);
+ break;
+ case 5:
+ Normal.SetXYZ(+0.452429, -0.454575, -0.767248);
+ break;
+ case 7:
+ Normal.SetXYZ(+0.443072, +0.443265, -0.779232);
+ break;
+ default:
+ cout << "ERROR:: Invalid DetNum " << MG << endl;
+ return 1; // Exit code
+ }
+}
+////////////////////////////////////////////////////////////////////////////////
+//Overloaded function definiton; this is for filling individual Ex histograms
+void DetectorSwitch(int MG, double Ex){
+ switch(MG){
+ case 1:
+ h1->Fill(Ex);
+ break;
+ case 2:
+ h2->Fill(Ex);
+ break;
+ case 3:
+ h3->Fill(Ex);
+ break;
+ case 4:
+ h4->Fill(Ex);
+ break;
+ case 5:
+ h5->Fill(Ex);
+ break;
+ case 7:
+ h7->Fill(Ex);
+ break;
+ default:
+ cout << "ERROR:: Invalid DetNum " << MG << endl;
+ return 1; // Exit code
+ }
+}
+////////////////////////////////////////////////////////////////////////////////
+void DrawOneHistogram(TH1D* hist, int mg, int colour, int fill, double *FitResultMatrixMG){
+ //Hist settings
+ hist->SetStats(0);
+ hist->SetLineColor(colour);
+ hist->SetFillStyle(fill);
+ hist->SetFillColor(colour);
+ hist->Draw("same");
+
+ if (flagDraw){
+ //Header
+ cout << noshowpos;
+ cout << "\n==================================================" << endl;
+ if (mg==6){
+ cout << "=---------------------- MG7 ---------------------=" << endl;
+ } else if (mg==0) {
+ cout << "=---------------------- SUM ---------------------=" << endl;
+ } else {
+ cout << "=---------------------- MG" << mg << " ---------------------=" << endl;
+ }
+ cout << showpos;
+ }
+
+ TFitResultPtr fit = hist->Fit("gaus","WQS"); //N=stop drawing, Q=stop writing
+ FillMatrix(FitResultMatrixMG,fit);
+}
+////////////////////////////////////////////////////////////////////////////////
+void InitiliseCanvas(double FitResultMatrix[7][5]){
+
+ //Canvas setup
TCanvas *canv = new TCanvas("canv","Ex Histograms",20,20,1600,800);
gStyle->SetOptStat(0);
canv->Divide(2,1,0.005,0.005,0);
@@ -58,61 +164,22 @@ void InitiliseCanvas(){
canv->cd(2)->SetBottomMargin(0.15);
gPad->SetTickx();
gPad->SetTicky();
-
+
+ //Histogram setup - Individual
canv->cd(1);
h1->SetMaximum(75.);
h1->GetXaxis()->SetTitle("Ex [MeV]");
h1->GetYaxis()->SetTitle("Counts");
-
- // ----- MG1 -----
- h1->SetStats(0);
- h1->SetLineColor(kRed);
- h1->SetFillStyle(3244);
- h1->SetFillColor(kRed);
- h1->Draw();
- h1->Fit("gaus","WQ"); //add N to stop it drawing
- // ----- MG2 -----
- h2->SetStats(0);
- h2->SetLineColor(kOrange);
- h2->SetFillStyle(3244);
- h2->SetFillColor(kOrange);
- h2->Draw("same");
- h2->Fit("gaus","WQ"); //add N to stop it drawing
-
- // ----- MG3 -----
- h3->SetStats(0);
- h3->SetLineColor(kGreen);
- h3->SetFillStyle(3344);
- h3->SetFillColor(kGreen);
- h3->Draw("same");
- h3->Fit("gaus","WQ"); //add N to stop it drawing
-
- // ----- MG4 -----
- h4->SetStats(0);
- h4->SetLineColor(kTeal);
- h4->SetFillStyle(3444);
- h4->SetFillColor(kTeal);
- h4->Draw("same");
- h4->Fit("gaus","WQ"); //add N to stop it drawing
-
- // ----- MG5 -----
- h5->SetStats(0);
- h5->SetLineColor(kBlue);
- h5->SetFillStyle(3544);
- h5->SetFillColor(kBlue);
- h5->Draw("same");
- h5->Fit("gaus","WQ"); //add N to stop it drawing
-
- // ----- MG7 -----
- h7->SetStats(0);
- h7->SetLineColor(kViolet);
- h7->SetFillStyle(3644);
- h7->SetFillColor(kViolet);
- h7->Draw("same");
- h7->Fit("gaus","WQ"); //add N to stop it drawing
-
- // Format legend
+ //Histogram draw - Individual
+ DrawOneHistogram(h1, 1, 632, 3244, FitResultMatrix[1]);
+ DrawOneHistogram(h2, 2, 800, 3244, FitResultMatrix[2]);
+ DrawOneHistogram(h3, 3, 416, 3344, FitResultMatrix[3]);
+ DrawOneHistogram(h4, 4, 840, 3444, FitResultMatrix[4]);
+ DrawOneHistogram(h5, 5, 600, 3544, FitResultMatrix[5]);
+ DrawOneHistogram(h7, 6, 880, 3644, FitResultMatrix[6]);
+
+ //Format legend
auto legend = new TLegend(0.15,0.7,0.35,0.9);
legend->AddEntry(h1,"MUGAST 1","f");
legend->AddEntry(h2,"MUGAST 2","f");
@@ -122,15 +189,18 @@ void InitiliseCanvas(){
legend->AddEntry(h7,"MUGAST 7","f");
legend->Draw();
- // ----- ALL -----
+ //Histogram setup - Sum
canv->cd(2);
- h->SetStats(0);
- h->SetLineColor(kBlack);
h->GetXaxis()->SetTitle("Ex [MeV]");
h->GetYaxis()->SetTitle("Counts");
- h->Draw();
- h->Fit("gaus", "WQ");
+
+ //Histogram draw - Sum
+ DrawOneHistogram(h, 0, 1, 0, FitResultMatrix[0]);
+
+ //Refresh
gPad->Update();
+
+ if(!flagDraw){delete canv;}
}
////////////////////////////////////////////////////////////////////////////////