diff --git a/Examples/Example1/RunToTreat.txt b/Examples/Example1/RunToTreat.txt
index 82908222e590c8760e20a5e47347834912b6c19c..4baf9cb3ab57671534853ca67caf68476db3a2d5 100644
--- a/Examples/Example1/RunToTreat.txt
+++ b/Examples/Example1/RunToTreat.txt
@@ -1,5 +1,5 @@
TTreeName
SimulatedTree
RootFileName
- ../../Outputs/Simulation/Example1.root
+ ./root/simulation/Example1.root
diff --git a/Examples/Example1/ShowResults.C b/Examples/Example1/ShowResults.C
index 6d7c2ed410b92cf005913800d9d94f8a5d95abd2..78853c616d53ac12f89a77969464b0fae51a44fb 100644
--- a/Examples/Example1/ShowResults.C
+++ b/Examples/Example1/ShowResults.C
@@ -6,6 +6,7 @@ TChain* chain=NULL ;
////////////////////////////////////////////////////////////////////////////////
void LoadCuts(){
+ std::cout << "Loading Cuts" << std::endl;
TFile* File_ETOF = new TFile("cuts/ETOF.root","READ");
ETOF = (TCutG*) File_ETOF->FindObjectAny("ETOF");
@@ -15,6 +16,8 @@ void LoadCuts(){
////////////////////////////////////////////////////////////////////////////////
void LoadChain(){
+
+ std::cout << "Loading Chain" << std::endl;
chain = new TChain("PhysicsTree");
chain->Add("root/analysis/Example1.root");
}
diff --git a/Examples/Example1/configs/ConfigMust2.dat b/Examples/Example1/configs/ConfigMust2.dat
index bebb0bf9ea475fdc1d492417228af5e96dbe37c5..4a76b40c82a19111e2d98809ddb225f33be4f602 100755
--- a/Examples/Example1/configs/ConfigMust2.dat
+++ b/Examples/Example1/configs/ConfigMust2.dat
@@ -14,4 +14,4 @@ ConfigMust2
DISABLE_ALL MM6
DISABLE_ALL MM7
DISABLE_ALL MM8
- CSI_SIZE 256
+ CSI_SIZE 32
diff --git a/Inputs/DetectorConfiguration/Example1.detector b/Inputs/DetectorConfiguration/Example1.detector
index e67ae572d0e1b77285c08345d01c6aaa5e07f11f..0b9645962390ac415160145e323619de434cef93 100644
--- a/Inputs/DetectorConfiguration/Example1.detector
+++ b/Inputs/DetectorConfiguration/Example1.detector
@@ -94,8 +94,6 @@ M2Telescope
CSI= 0
VIS= all
-%%%%%%%%%%%%%%%%%%%%%
-SSSDArray
%%%%%%%%%% Det 1 %%%%%%%%
SSSD
A= 17.61 9.85 104.11 mm
diff --git a/NPLib/Detectors/PISTA/CMakeLists.txt b/NPLib/Detectors/PISTA/CMakeLists.txt
index d2d58fc5dc5acac95e36eb1ed0279536c17d0d43..921869f3e6e7b2351ff87502dac0161a1a09e118 100644
--- a/NPLib/Detectors/PISTA/CMakeLists.txt
+++ b/NPLib/Detectors/PISTA/CMakeLists.txt
@@ -4,9 +4,14 @@ add_custom_command(OUTPUT TPISTADataDict.cxx COMMAND ../../scripts/build_dict.sh
add_custom_command(OUTPUT TFPMWDataDict.cxx COMMAND ../../scripts/build_dict.sh TFPMWData.h TFPMWDataDict.cxx TFPMWData.rootmap libNPPISTA.dylib DEPENDS TFPMWData.h)
+add_custom_command(OUTPUT TFPMWPhysicsDict.cxx COMMAND ../../scripts/build_dict.sh TFPMWPhysics.h TFPMWPhysicsDict.cxx TFPMWPhysics.rootmap libNPPISTA.dylib DEPENDS TFPMWPhysics.h)
+
add_custom_command(OUTPUT TICDataDict.cxx COMMAND ../../scripts/build_dict.sh TICData.h TICDataDict.cxx TICData.rootmap libNPPISTA.dylib DEPENDS TICData.h)
-add_library(NPPISTA SHARED TPISTASpectra.cxx TPISTAData.cxx TPISTAPhysics.cxx TPISTADataDict.cxx TPISTAPhysicsDict.cxx TFPMWData.cxx TFPMWDataDict.cxx TICData.cxx TICDataDict.cxx)
+add_custom_command(OUTPUT TICPhysicsDict.cxx COMMAND ../../scripts/build_dict.sh TICPhysics.h TICPhysicsDict.cxx TICPhysics.rootmap libNPPISTA.dylib DEPENDS TICPhysics.h)
+
+add_library(NPPISTA SHARED TPISTASpectra.cxx TPISTAData.cxx TPISTAPhysics.cxx TPISTADataDict.cxx TPISTAPhysicsDict.cxx TFPMWData.cxx TFPMWDataDict.cxx TFPMWPhysics.cxx TFPMWPhysicsDict.cxx TICData.cxx TICDataDict.cxx TICPhysics.cxx TICPhysicsDict.cxx)
+
target_link_libraries(NPPISTA ${ROOT_LIBRARIES} NPCore)
-install(FILES TPISTAData.h TPISTAPhysics.h TPISTASpectra.h TFPMWData.h TICData.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
+install(FILES TPISTAData.h TPISTAPhysics.h TPISTASpectra.h TFPMWData.h TFPMWPhysics.h TICData.h TICPhysics.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
diff --git a/NPLib/Detectors/PISTA/TFPMWPhysics.cxx b/NPLib/Detectors/PISTA/TFPMWPhysics.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..520a445673628113880ba9a873b0e9a3b36ed85a
--- /dev/null
+++ b/NPLib/Detectors/PISTA/TFPMWPhysics.cxx
@@ -0,0 +1,393 @@
+/*****************************************************************************
+ * 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: P. Morfouace contact address: pierre.morfouace@cea.fr *
+ * *
+ * Creation Date : October 2023 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class hold FPMW Treated data *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+
+#include "TFPMWPhysics.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"
+
+// ROOT
+#include "TChain.h"
+
+ClassImp(TFPMWPhysics);
+
+///////////////////////////////////////////////////////////////////////////
+TFPMWPhysics::TFPMWPhysics()
+ : m_EventData(new TFPMWData),
+ m_PreTreatedData(new TFPMWData),
+ m_EventPhysics(this),
+ m_NumberOfDetectors(0){
+ }
+
+///////////////////////////////////////////////////////////////////////////
+/// A usefull method to bundle all operation to add a detector
+void TFPMWPhysics::AddDetector(TVector3 Pos){
+ // In That simple case nothing is done
+ // Typically for more complex detector one would calculate the relevant
+ // positions (stripped silicon) or angles (gamma array)
+ DetPosX.push_back(Pos.X());
+ DetPosY.push_back(Pos.Y());
+ DetPosZ.push_back(Pos.Z());
+
+ m_NumberOfDetectors++;
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TFPMWPhysics::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 TFPMWPhysics::BuildSimplePhysicalEvent() {
+ BuildPhysicalEvent();
+}
+
+
+
+///////////////////////////////////////////////////////////////////////////
+void TFPMWPhysics::BuildPhysicalEvent() {
+ PreTreat();
+
+ unsigned int sizeX = m_PreTreatedData->GetFPMWMultX();
+ for(unsigned int i=0; i<sizeX; i++){
+ DetectorHitX.insert(m_PreTreatedData->GetFPMW_DetX(i));
+ }
+
+ unsigned int sizeY = m_PreTreatedData->GetFPMWMultY();
+ for(unsigned int i=0; i<sizeY; i++){
+ if(DetectorHitX.find(m_PreTreatedData->GetFPMW_DetY(i)) != DetectorHitX.end()){
+ DetectorHit.insert(m_PreTreatedData->GetFPMW_DetY(i));
+ }
+ }
+ unsigned int sizeDet = DetectorHit.size();
+
+
+ for(unsigned int i=0; i<sizeX; i++){
+ int StrX = m_PreTreatedData->GetFPMW_StripX(i);
+ int NX = m_PreTreatedData->GetFPMW_DetX(i);
+ double QX = m_PreTreatedData->GetFPMW_ChargeX(i);
+ if(DetectorHit.find(NX) != DetectorHit.end()){
+ MapX[NX].push_back(std::make_pair(StrX,QX));
+ QSumX[NX] += QX;
+ if(MaxQX.find(NX)==MaxQX.end() || MaxQX[NX].second<QX){
+ MaxQX[NX] = make_pair(StrX,QX);
+ }
+ }
+ }
+ for(unsigned int i=0; i<sizeY; i++){
+ int StrY = m_PreTreatedData->GetFPMW_StripY(i);
+ int NY = m_PreTreatedData->GetFPMW_DetY(i);
+ double QY = m_PreTreatedData->GetFPMW_ChargeY(i);
+ if(DetectorHit.find(NY) != DetectorHit.end()){
+ MapY[NY].push_back(std::make_pair(StrY,QY));
+ QSumY[NY] += QY;
+ if(MaxQY.find(NY)==MaxQY.end() || MaxQY[NY].second<QY){
+ MaxQY[NY] = make_pair(StrY,QY);
+ }
+ }
+ }
+
+ for(auto &DetN : DetectorHit){
+ double PosX = AnalyticHyperbolicSecant(MaxQX[DetN],MapX[DetN]);
+ double PosY = AnalyticHyperbolicSecant(MaxQY[DetN],MapY[DetN]);
+
+ int sx0 = (int) PosX;
+ int sx1 = sx0+1;
+ int sy0 = (int) PosY;
+ int sy1 = sy0+1;
+
+ DetectorNbr.push_back(DetN);
+ ChargeX.push_back(QSumX[DetN]);
+ ChargeY.push_back(QSumY[DetN]);
+ PositionX.push_back(PosX);
+ PositionY.push_back(PosY);
+ }
+
+}
+
+/////////////////////////////////////////////////////////////////////
+double TFPMWPhysics::AnalyticHyperbolicSecant(std::pair<int,double>& MaxQ,std::vector<std::pair<int,double>>& Map){
+ double sech = -1000 ;
+
+ // std::cout << "test AnH 1" << std::endl;
+ if(MaxQ.second > 0 && MaxQ.first > 2 && MaxQ.first<996){
+ // if(Buffer_Q[MaxQ.first-1+1]==0||Buffer_Q[MaxQ.first-1-1]==0)
+ // return sech;
+ // std::cout << "test AnH 2" << std::endl;
+ double q2 = MaxQ.second;
+ double q1 = 0,q3 = 0;
+ for(auto &strip : Map){
+ if(strip.first == MaxQ.first - 1){
+ q1 = strip.second;
+ }
+ else if(strip.first == MaxQ.first + 1){
+ q3 = strip.second;
+ }
+ }
+ //std::cout << "test q " << q1 << " " << q2 << " " << q3 << std::endl;
+ double vs[6];
+ // std::cout << "test AnH 3" << std::endl;
+ if(q1 > 0 && q3 > 0)
+ {
+ // QsumSample[DetNum].push_back(QSum);
+ vs[0] = sqrt(q2/q3);
+ vs[1] = sqrt(q2/q1);
+ vs[2] = 0.5*(vs[0] + vs[1]);
+ vs[3] = log( vs[2] + sqrt(vs[2]*vs[2]-1.0) );
+ vs[4] = abs((vs[0] - vs[1])/(2.0*sinh(vs[3])));
+ vs[5] = 0.5*log( (1.0+vs[4])/(1.0-vs[4]) ) ;
+ // std::cout << "test AnH 4" << std::endl;
+
+ if ( q3>q1 )
+ sech = MaxQ.first + vs[5]/vs[3] ;
+ else
+ sech = MaxQ.first - vs[5]/vs[3] ;
+ //std::cout << "test sech " << sech << std::endl;
+ }
+ }
+
+ return sech ;
+}
+
+
+///////////////////////////////////////////////////////////////////////////
+void TFPMWPhysics::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 mysizeX = m_EventData->GetFPMWMultX();
+ for (unsigned int i = 0; i < mysizeX ; ++i) {
+ int det = m_EventData->GetFPMW_DetX(i);
+ int strip = m_EventData->GetFPMW_StripX(i);
+ double QX = m_EventData->GetFPMW_ChargeX(i);
+ double Qcal = Cal->ApplyCalibration("FPMW/DET"+NPL::itoa(det)+"_STRIP"+NPL::itoa(strip),QX);
+ m_PreTreatedData->SetFPMW_DetX(det);
+ m_PreTreatedData->SetFPMW_StripX(strip);
+ m_PreTreatedData->SetFPMW_ChargeX(Qcal);
+ }
+
+ unsigned int mysizeY = m_EventData->GetFPMWMultY();
+ for (unsigned int i = 0; i < mysizeY ; ++i) {
+ int det = m_EventData->GetFPMW_DetY(i);
+ int strip = m_EventData->GetFPMW_StripY(i);
+ double QY = m_EventData->GetFPMW_ChargeY(i);
+ double Qcal = Cal->ApplyCalibration("FPMW/DET"+NPL::itoa(det)+"_STRIP"+NPL::itoa(strip),QY);
+ m_PreTreatedData->SetFPMW_DetY(det);
+ m_PreTreatedData->SetFPMW_StripY(strip);
+ m_PreTreatedData->SetFPMW_ChargeY(Qcal);
+ }
+
+}
+
+
+
+///////////////////////////////////////////////////////////////////////////
+void TFPMWPhysics::ReadAnalysisConfig() {
+ bool ReadingStatus = false;
+
+ // path to file
+ string FileName = "./configs/ConfigFPMW.dat";
+
+ // open analysis config file
+ ifstream AnalysisConfigFile;
+ AnalysisConfigFile.open(FileName.c_str());
+
+ if (!AnalysisConfigFile.is_open()) {
+ cout << " No ConfigFPMW.dat found: Default parameter loaded for Analayis " << FileName << endl;
+ return;
+ }
+ cout << " Loading user parameter for Analysis from ConfigFPMW.dat " << endl;
+
+ // Save it in a TAsciiFile
+ TAsciiFile* asciiConfig = RootOutput::getInstance()->GetAsciiFileAnalysisConfig();
+ asciiConfig->AppendLine("%%% ConfigFPMW.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 = "ConfigFPMW";
+ 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 {
+ ReadingStatus = false;
+ }
+ }
+ }
+}
+
+
+///////////////////////////////////////////////////////////////////////////
+void TFPMWPhysics::Clear() {
+ DetectorNbr.clear();
+ PositionX.clear();
+ PositionY.clear();
+ ChargeX.clear();
+ ChargeY.clear();
+
+ DetectorHitX.clear();
+ DetectorHit.clear();
+
+ MapX.clear();
+ MapY.clear();
+ MaxQX.clear();
+ MaxQY.clear();
+}
+
+
+
+///////////////////////////////////////////////////////////////////////////
+void TFPMWPhysics::ReadConfiguration(NPL::InputParser parser) {
+ vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("FPMW");
+ 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 << "//// FPMW " << i+1 << endl;
+
+ TVector3 Pos = blocks[i]->GetTVector3("POS","mm");
+ AddDetector(Pos);
+ }
+ else if(blocks[i]->HasTokenList(sphe)){
+ if(NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << endl << "//// FPMW " << 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 TFPMWPhysics::AddParameterToCalibrationManager() {
+ CalibrationManager* Cal = CalibrationManager::getInstance();
+
+ for(int i = 0; i < m_NumberOfDetectors; i++){
+ for(int s = 0; s < 996; s++){
+ Cal->AddParameter("FPMW","DET"+NPL::itoa(i)+"_STRIP"+NPL::itoa(s),"FPMW_DET"+NPL::itoa(i)+"_STRIP"+NPL::itoa(s));
+ }
+ }
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TFPMWPhysics::InitializeRootInputRaw() {
+ TChain* inputChain = RootInput::getInstance()->GetChain();
+ inputChain->SetBranchStatus("FPMW", true );
+ inputChain->SetBranchAddress("FPMW", &m_EventData );
+}
+
+
+
+///////////////////////////////////////////////////////////////////////////
+void TFPMWPhysics::InitializeRootInputPhysics() {
+ TChain* inputChain = RootInput::getInstance()->GetChain();
+ inputChain->SetBranchAddress("FPMW", &m_EventPhysics);
+}
+
+
+
+///////////////////////////////////////////////////////////////////////////
+void TFPMWPhysics::InitializeRootOutput() {
+ TTree* outputTree = RootOutput::getInstance()->GetTree();
+ outputTree->Branch("FPMW", "TFPMWPhysics", &m_EventPhysics);
+}
+
+
+
+////////////////////////////////////////////////////////////////////////////////
+// Construct Method to be pass to the DetectorFactory //
+////////////////////////////////////////////////////////////////////////////////
+NPL::VDetector* TFPMWPhysics::Construct() {
+ return (NPL::VDetector*) new TFPMWPhysics();
+}
+
+
+
+////////////////////////////////////////////////////////////////////////////////
+// Registering the construct method to the factory //
+////////////////////////////////////////////////////////////////////////////////
+extern "C"{
+ class proxy_FPMW{
+ public:
+ proxy_FPMW(){
+ NPL::DetectorFactory::getInstance()->AddToken("FPMW","FPMW");
+ NPL::DetectorFactory::getInstance()->AddDetector("FPMW",TFPMWPhysics::Construct);
+ }
+ };
+
+ proxy_FPMW p_FPMW;
+}
+
diff --git a/NPLib/Detectors/PISTA/TFPMWPhysics.h b/NPLib/Detectors/PISTA/TFPMWPhysics.h
new file mode 100644
index 0000000000000000000000000000000000000000..d78ac457498bd5b0af2cc94063c14e0b4f432789
--- /dev/null
+++ b/NPLib/Detectors/PISTA/TFPMWPhysics.h
@@ -0,0 +1,176 @@
+#ifndef TFPMWPHYSICS_H
+#define TFPMWPHYSICS_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: P. Morfouace contact address: pierre.morfouace@cea.fr *
+ * *
+ * Creation Date : October 2023 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class hold FPMW Treated data *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+
+// C++ headers
+#include <vector>
+#include <map>
+#include <set>
+#include <string>
+using namespace std;
+
+// ROOT headers
+#include "TObject.h"
+#include "TH1.h"
+#include "TVector3.h"
+#include "TSpline.h"
+// NPTool headers
+#include "TFPMWData.h"
+#include "NPCalibrationManager.h"
+#include "NPVDetector.h"
+#include "NPInputParser.h"
+
+
+
+class TFPMWPhysics : public TObject, public NPL::VDetector {
+ //////////////////////////////////////////////////////////////
+ // constructor and destructor
+ public:
+ TFPMWPhysics();
+ ~TFPMWPhysics() {};
+
+
+ //////////////////////////////////////////////////////////////
+ // 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> PositionX;
+ vector<double> PositionY;
+ vector<double> ChargeX;
+ vector<double> ChargeY;
+
+ private:
+ std::set<int> DetectorHitX;//!
+ std::set<int> DetectorHit;//!
+ vector<double> DetPosX;//!
+ vector<double> DetPosY;//!
+ vector<double> DetPosZ;//!
+ vector<vector<double>> Buffer_X_Q;//!
+ vector<vector<double>> Buffer_Y_Q;//!
+
+ map<int, vector<pair<int,double>>> MapX;//!
+ map<int, vector<pair<int,double>>> MapY;//!
+ map<int, pair<int, double>> MaxQX;//!
+ map<int, pair<int, double>> MaxQY;//!
+ map<int, double> QSumX;//!
+ map<int, double> QSumY;//!
+
+
+ /// 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();}
+
+
+ //////////////////////////////////////////////////////////////
+ // specific methods to FPMW 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 TFPMWData object to TFPMWPhysics.
+ // needed for online analysis for example
+ void SetRawDataPointer(TFPMWData* rawDataPointer) {m_EventData = rawDataPointer;}
+
+ double fThresholdX;
+ double fThresholdY;
+
+ double AnalyticHyperbolicSecant(std::pair<int,double>& MaxQ, std::vector<std::pair<int, double>>& Map);
+
+
+ // objects are not written in the TTree
+ private:
+ TFPMWData* m_EventData; //!
+ TFPMWData* m_PreTreatedData; //!
+ TFPMWPhysics* m_EventPhysics; //!
+
+ // getters for raw and pre-treated data object
+ public:
+ TFPMWData* GetRawData() const {return m_EventData;}
+ TFPMWData* GetPreTreatedData() const {return m_PreTreatedData;}
+
+ // parameters used in the analysis
+ private:
+ double m_E_Threshold; //!
+
+ // number of detectors
+ private:
+ int m_NumberOfDetectors; //!
+
+ // Static constructor to be passed to the Detector Factory
+ public:
+ static NPL::VDetector* Construct();
+
+ ClassDef(TFPMWPhysics,1) // FPMWPhysics structure
+};
+#endif
diff --git a/NPLib/Detectors/PISTA/TICPhysics.cxx b/NPLib/Detectors/PISTA/TICPhysics.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..34219350aa7e020e94fe566801df939efcd95f33
--- /dev/null
+++ b/NPLib/Detectors/PISTA/TICPhysics.cxx
@@ -0,0 +1,274 @@
+/*****************************************************************************
+ * 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: P. Morfouace contact address: pierre.morfouace@cea.fr *
+ * *
+ * Creation Date : October 2023 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class hold IC Treated data *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+
+#include "TICPhysics.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"
+
+// ROOT
+#include "TChain.h"
+
+ClassImp(TICPhysics);
+
+///////////////////////////////////////////////////////////////////////////
+TICPhysics::TICPhysics()
+ : m_EventData(new TICData),
+ m_PreTreatedData(new TICData),
+ m_EventPhysics(this),
+ m_NumberOfDetectors(0){
+ }
+
+///////////////////////////////////////////////////////////////////////////
+/// A usefull method to bundle all operation to add a detector
+void TICPhysics::AddDetector(TVector3 Pos){
+ // In That simple case nothing is done
+ // Typically for more complex detector one would calculate the relevant
+ // positions (stripped silicon) or angles (gamma array)
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TICPhysics::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 TICPhysics::BuildSimplePhysicalEvent() {
+ BuildPhysicalEvent();
+}
+
+
+
+///////////////////////////////////////////////////////////////////////////
+void TICPhysics::BuildPhysicalEvent() {
+ PreTreat();
+
+ double fIC[11];
+
+ int size = m_PreTreatedData->GetICMult();
+ for(int i=0; i<size; i++){
+ fIC[i] = m_PreTreatedData->GetIC_Charge(i);
+ }
+
+ if(fIC[1]>0 && fIC[5]>0){
+ DE = 0.5*(fIC[0] + fIC[1] + fIC[2] + fIC[3]) + fIC[4];
+ Eres = fIC[5] + fIC[6] + fIC[7] + fIC[8] + fIC[9];
+ }
+ else{
+ DE = -100;
+ Eres = -100;
+ }
+
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TICPhysics::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->GetICMult();
+ for (unsigned int i = 0; i < mysize ; ++i) {
+ m_PreTreatedData->SetIC_Charge(m_EventData->GetIC_Charge(i));
+ m_PreTreatedData->SetIC_Section(m_EventData->GetIC_Section(i));
+ }
+}
+
+
+
+///////////////////////////////////////////////////////////////////////////
+void TICPhysics::ReadAnalysisConfig() {
+ bool ReadingStatus = false;
+
+ // path to file
+ string FileName = "./configs/ConfigIC.dat";
+
+ // open analysis config file
+ ifstream AnalysisConfigFile;
+ AnalysisConfigFile.open(FileName.c_str());
+
+ if (!AnalysisConfigFile.is_open()) {
+ cout << " No ConfigIC.dat found: Default parameter loaded for Analayis " << FileName << endl;
+ return;
+ }
+ cout << " Loading user parameter for Analysis from ConfigIC.dat " << endl;
+
+ // Save it in a TAsciiFile
+ TAsciiFile* asciiConfig = RootOutput::getInstance()->GetAsciiFileAnalysisConfig();
+ asciiConfig->AppendLine("%%% ConfigIC.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 = "ConfigIC";
+ 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 {
+ ReadingStatus = false;
+ }
+ }
+ }
+}
+
+
+///////////////////////////////////////////////////////////////////////////
+void TICPhysics::Clear() {
+ DE = -100;
+ Eres = -100;
+}
+
+
+
+///////////////////////////////////////////////////////////////////////////
+void TICPhysics::ReadConfiguration(NPL::InputParser parser) {
+ vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("IC");
+ 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 << "//// IC " << i+1 << endl;
+
+ TVector3 Pos = blocks[i]->GetTVector3("POS","mm");
+ AddDetector(Pos);
+ }
+ else if(blocks[i]->HasTokenList(sphe)){
+ if(NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << endl << "//// IC " << 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 TICPhysics::AddParameterToCalibrationManager() {
+ CalibrationManager* Cal = CalibrationManager::getInstance();
+
+ for(int sec = 0; sec < m_NumberOfDetectors; sec++){
+ Cal->AddParameter("IC","SEC"+NPL::itoa(sec+1)+"_ALIGN","IC_SEC"+NPL::itoa(sec+1)+"_ALIGN");
+ }
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TICPhysics::InitializeRootInputRaw() {
+ TChain* inputChain = RootInput::getInstance()->GetChain();
+ inputChain->SetBranchStatus("IC", true );
+ inputChain->SetBranchAddress("IC", &m_EventData );
+}
+
+
+
+///////////////////////////////////////////////////////////////////////////
+void TICPhysics::InitializeRootInputPhysics() {
+ TChain* inputChain = RootInput::getInstance()->GetChain();
+ inputChain->SetBranchAddress("IC", &m_EventPhysics);
+}
+
+
+
+///////////////////////////////////////////////////////////////////////////
+void TICPhysics::InitializeRootOutput() {
+ TTree* outputTree = RootOutput::getInstance()->GetTree();
+ outputTree->Branch("IC", "TICPhysics", &m_EventPhysics);
+}
+
+
+
+////////////////////////////////////////////////////////////////////////////////
+// Construct Method to be pass to the DetectorFactory //
+////////////////////////////////////////////////////////////////////////////////
+NPL::VDetector* TICPhysics::Construct() {
+ return (NPL::VDetector*) new TICPhysics();
+}
+
+
+
+////////////////////////////////////////////////////////////////////////////////
+// Registering the construct method to the factory //
+////////////////////////////////////////////////////////////////////////////////
+extern "C"{
+ class proxy_IC{
+ public:
+ proxy_IC(){
+ NPL::DetectorFactory::getInstance()->AddToken("IC","IC");
+ NPL::DetectorFactory::getInstance()->AddDetector("IC",TICPhysics::Construct);
+ }
+ };
+
+ proxy_IC p_IC;
+}
+
diff --git a/NPLib/Detectors/PISTA/TICPhysics.h b/NPLib/Detectors/PISTA/TICPhysics.h
new file mode 100644
index 0000000000000000000000000000000000000000..beac07d38f693281f85f4015b023b6ae33521c2a
--- /dev/null
+++ b/NPLib/Detectors/PISTA/TICPhysics.h
@@ -0,0 +1,152 @@
+#ifndef TICPHYSICS_H
+#define TICPHYSICS_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: P. Morfouace contact address: pierre.morfouace@cea.fr *
+ * *
+ * Creation Date : October 2023 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class hold IC Treated data *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+
+// C++ headers
+#include <vector>
+#include <map>
+#include <set>
+#include <string>
+using namespace std;
+
+// ROOT headers
+#include "TObject.h"
+#include "TH1.h"
+#include "TVector3.h"
+#include "TSpline.h"
+// NPTool headers
+#include "TICData.h"
+#include "NPCalibrationManager.h"
+#include "NPVDetector.h"
+#include "NPInputParser.h"
+
+
+
+class TICPhysics : public TObject, public NPL::VDetector {
+ //////////////////////////////////////////////////////////////
+ // constructor and destructor
+ public:
+ TICPhysics();
+ ~TICPhysics() {};
+
+
+ //////////////////////////////////////////////////////////////
+ // 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:
+ double DE;
+ double Eres;
+
+ private:
+
+ /// 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();}
+
+
+ //////////////////////////////////////////////////////////////
+ // specific methods to IC 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 TICData object to TICPhysics.
+ // needed for online analysis for example
+ void SetRawDataPointer(TICData* rawDataPointer) {m_EventData = rawDataPointer;}
+
+ // objects are not written in the TTree
+ private:
+ TICData* m_EventData; //!
+ TICData* m_PreTreatedData; //!
+ TICPhysics* m_EventPhysics; //!
+
+ // getters for raw and pre-treated data object
+ public:
+ TICData* GetRawData() const {return m_EventData;}
+ TICData* GetPreTreatedData() const {return m_PreTreatedData;}
+
+ // parameters used in the analysis
+ private:
+ double m_E_Threshold; //!
+
+ // number of detectors
+ private:
+ int m_NumberOfDetectors; //!
+
+ // Static constructor to be passed to the Detector Factory
+ public:
+ static NPL::VDetector* Construct();
+
+ ClassDef(TICPhysics,1) // ICPhysics structure
+};
+#endif
diff --git a/NPLib/Physics/NPReaction.cxx b/NPLib/Physics/NPReaction.cxx
index b548dc168d24c69a70cafa4742fadda43c165b97..46aefd5692a7b4f768dab5145a18007890b6ed40 100644
--- a/NPLib/Physics/NPReaction.cxx
+++ b/NPLib/Physics/NPReaction.cxx
@@ -261,10 +261,10 @@ void Reaction::KineRelativistic(double& ThetaLab3, double& KineticEnergyLab3, do
// case of inverse kinematics
double theta = fThetaCM;
- /*if (m1 > m2){
+ if (m1 > m2){
theta = M_PI - fThetaCM;
//fThetaCM = M_PI - fThetaCM;
- }*/
+ }
fEnergyImpulsionCM_3 = TLorentzVector(pCM_3 * sin(theta), 0, pCM_3 * cos(theta), ECM_3);
fEnergyImpulsionCM_4 = fTotalEnergyImpulsionCM - fEnergyImpulsionCM_3;
diff --git a/NPSimulation/Core/SteppingAction.cc b/NPSimulation/Core/SteppingAction.cc
index 2069571a8c818c8b9c22889e9a09327f81f38ca8..d82804352761b4c0e2b141f5ff5c0e265e436063 100644
--- a/NPSimulation/Core/SteppingAction.cc
+++ b/NPSimulation/Core/SteppingAction.cc
@@ -28,18 +28,18 @@
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
SteppingAction::SteppingAction() {
- m_record_track = false;
- m_tree = RootOutput::getInstance()->GetTree();
+ m_record_track = NPOptionManager::getInstance()->GetRecordTrack();
+ m_tree = RootOutput::getInstance()->GetTree();
- m_First = true;
- ParticleName = "";
+ m_First = true;
+ ParticleName = "";
KineticEnergy = -1000;
- Theta = -1000;
- Phi = -1000;
- Mass = -1000;
- Charge = -1000;
- Z = -1000;
- A = -1000;
+ Theta = -1000;
+ Phi = -1000;
+ Mass = -1000;
+ Charge = -1000;
+ Z = -1000;
+ A = -1000;
Momentum.setX(-1000);
Momentum.setY(-1000);
Momentum.setZ(-1000);
@@ -47,15 +47,14 @@ SteppingAction::SteppingAction() {
Position.setY(-1000);
Position.setZ(-1000);
LastStepNumber = -1000;
- VolumeName = "";
- nClear = 0;
- TrackID = -1000;
+ VolumeName = "";
+ nClear = 0;
+ TrackID = -1000;
m_TrackInfo = new TTrackInfo();
AttachTrackInfo();
if (!RootOutput::getInstance()->GetTree()->FindBranch("TrackInfo"))
- RootOutput::getInstance()->GetTree()->Branch("TrackInfo", "TTrackInfo",
- &m_TrackInfo);
+ RootOutput::getInstance()->GetTree()->Branch("TrackInfo", "TTrackInfo", &m_TrackInfo);
// Attach track info class to m_tree
}
@@ -64,8 +63,7 @@ SteppingAction::SteppingAction() {
void SteppingAction::AttachTrackInfo() {
// Reasssigned the branch address
if (RootOutput::getInstance()->GetTree()->FindBranch("TrackInfo"))
- RootOutput::getInstance()->GetTree()->SetBranchAddress("TrackInfo",
- &m_TrackInfo);
+ RootOutput::getInstance()->GetTree()->SetBranchAddress("TrackInfo", &m_TrackInfo);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -82,33 +80,31 @@ void SteppingAction::TrackRecording(const G4Step* step) {
m_TrackInfo->Clear();
m_First = false;
- G4Track* track = step->GetTrack();
- int StepNumber = track->GetCurrentStepNumber();
+ G4Track* track = step->GetTrack();
+ int StepNumber = track->GetCurrentStepNumber();
- if (eventID
- < G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID()) {
+ if (eventID < G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID()) {
m_TrackInfo->Clear();
nClear++;
}
if (StepNumber == 1) {
- ParticleName = track->GetParticleDefinition()->GetParticleName();
+ ParticleName = track->GetParticleDefinition()->GetParticleName();
KineticEnergy = track->GetDynamicParticle()->GetKineticEnergy();
- Mass = track->GetDynamicParticle()->GetMass();
- Charge = track->GetDynamicParticle()->GetCharge();
- Z = track->GetParticleDefinition()->GetAtomicNumber();
- A = track->GetParticleDefinition()->GetAtomicMass();
- Momentum = track->GetDynamicParticle()->GetMomentum();
- Theta = Momentum.theta() * 180. / M_PI;
- Phi = Momentum.phi() * 180. / M_PI;
- Position = track->GetPosition();
+ Mass = track->GetDynamicParticle()->GetMass();
+ Charge = track->GetDynamicParticle()->GetCharge();
+ Z = track->GetParticleDefinition()->GetAtomicNumber();
+ A = track->GetParticleDefinition()->GetAtomicMass();
+ Momentum = track->GetDynamicParticle()->GetMomentum();
+ Theta = Momentum.theta() * 180. / M_PI;
+ Phi = Momentum.phi() * 180. / M_PI;
+ Position = track->GetPosition();
G4VPhysicalVolume* volume = track->GetVolume();
- VolumeName = volume->GetName();
- TrackID = track->GetTrackID();
+ VolumeName = volume->GetName();
+ TrackID = track->GetTrackID();
double c_light = 299.792458; // To go from T.m to MeV/e
- double Brho = sqrt(KineticEnergy * KineticEnergy + 2 * KineticEnergy * Mass)
- / (c_light * Charge);
+ double Brho = sqrt(KineticEnergy * KineticEnergy + 2 * KineticEnergy * Mass) / (c_light * Charge);
m_TrackInfo->SetKineticEnergy(KineticEnergy);
m_TrackInfo->SetTheta(Theta);