Merge branch 'NPTool.2.dev' of gitlab.in2p3.fr:np/nptool into NPTool.2.dev

Examples/Example1/RunToTreat.txt

 TTreeName 
 	SimulatedTree
 RootFileName 
-	../../Outputs/Simulation/Example1.root
+	./root/simulation/Example1.root
 

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");
 }

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

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

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})
 

NPLib/Detectors/PISTA/TFPMWPhysics.cxx

+/*****************************************************************************
+ * 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;
+}
+

NPLib/Detectors/PISTA/TFPMWPhysics.h

+#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

NPLib/Detectors/PISTA/TICPhysics.cxx

+/*****************************************************************************
+ * 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;
+}
+

NPLib/Detectors/PISTA/TICPhysics.h

+#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

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;
 

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);