From 801efaff4ad21ff5117e0c1840f3219311fb5ab7 Mon Sep 17 00:00:00 2001 From: Adrien Matta <matta@lpccaen.in2p3.fr> Date: Wed, 30 Mar 2022 11:39:17 +0200 Subject: [PATCH] * Initial commit of SuperX3 based on file provided by Hari Pai (ELI-NP) --- Examples/Example1/Analysis.cxx | 8 +- Examples/Example1/ShowResults.C | 2 +- Inputs/EventGenerator/np.source | 4 +- NPLib/Detectors/Minos/TMinosPhysics.cxx | 12 +- NPLib/Detectors/SuperX3/CMakeLists.txt | 6 + NPLib/Detectors/SuperX3/TSuperX3Data.cxx | 62 ++ NPLib/Detectors/SuperX3/TSuperX3Data.h | 138 ++++ NPLib/Detectors/SuperX3/TSuperX3Physics.cxx | 698 ++++++++++++++++++ NPLib/Detectors/SuperX3/TSuperX3Physics.h | 241 ++++++ NPLib/Detectors/SuperX3/TSuperX3Spectra.cxx | 66 ++ NPLib/Detectors/SuperX3/TSuperX3Spectra.h | 58 ++ NPSimulation/Detectors/SuperX3/CMakeLists.txt | 2 + NPSimulation/Detectors/SuperX3/SuperX3.cc | 366 +++++++++ NPSimulation/Detectors/SuperX3/SuperX3.hh | 174 +++++ Projects/S034/PhysicsListOption.txt | 2 +- Projects/Strasse/PhysicsListOption.txt | 2 +- .../Strasse/geometry/catana_only.detector | 2 +- .../Strasse/geometry/strasse_CAD.detector | 18 +- .../geometry/strasse_July2021.detector | 8 +- Projects/Strasse/reaction/C12_p2p.reaction | 2 +- Projects/Strasse/reaction/Ca54.source | 1 + Projects/SuperX3/Analysis.cxx | 83 +++ Projects/SuperX3/Analysis.h | 52 ++ Projects/SuperX3/CMakeLists.txt | 5 + Projects/SuperX3/SuperX3.detector | 20 + Projects/SuperX3/energy_loss/.gitignore | 0 Projects/SuperX3/project.config | 5 + Projects/SuperX3/root/analysis/.gitignore | 0 Projects/SuperX3/root/simulation/.gitignore | 0 29 files changed, 2005 insertions(+), 32 deletions(-) create mode 100644 NPLib/Detectors/SuperX3/CMakeLists.txt create mode 100644 NPLib/Detectors/SuperX3/TSuperX3Data.cxx create mode 100644 NPLib/Detectors/SuperX3/TSuperX3Data.h create mode 100644 NPLib/Detectors/SuperX3/TSuperX3Physics.cxx create mode 100644 NPLib/Detectors/SuperX3/TSuperX3Physics.h create mode 100644 NPLib/Detectors/SuperX3/TSuperX3Spectra.cxx create mode 100644 NPLib/Detectors/SuperX3/TSuperX3Spectra.h create mode 100644 NPSimulation/Detectors/SuperX3/CMakeLists.txt create mode 100644 NPSimulation/Detectors/SuperX3/SuperX3.cc create mode 100644 NPSimulation/Detectors/SuperX3/SuperX3.hh create mode 100644 Projects/SuperX3/Analysis.cxx create mode 100644 Projects/SuperX3/Analysis.h create mode 100644 Projects/SuperX3/CMakeLists.txt create mode 100644 Projects/SuperX3/SuperX3.detector create mode 100644 Projects/SuperX3/energy_loss/.gitignore create mode 100644 Projects/SuperX3/project.config create mode 100644 Projects/SuperX3/root/analysis/.gitignore create mode 100644 Projects/SuperX3/root/simulation/.gitignore diff --git a/Examples/Example1/Analysis.cxx b/Examples/Example1/Analysis.cxx index 0a60fe37d..02e68e4a8 100644 --- a/Examples/Example1/Analysis.cxx +++ b/Examples/Example1/Analysis.cxx @@ -61,10 +61,10 @@ void Analysis::Init(){ Si_Y_M2 = 0; TargetThickness = m_DetectorManager->GetTargetThickness(); // Energy loss table: the G4Table are generated by the simulation - He3CD2 = EnergyLoss("Example/He3_CD2.G4table","G4Table",100 ); - He3Al = EnergyLoss("Example/He3_Al.G4table","G4Table",10); - He3Si = EnergyLoss("Example/He3_Si.G4table","G4Table",10); - Li11CD2 = EnergyLoss("Example/Li11_CD2.G4table","G4Table",100); + He3CD2 = EnergyLoss("energy_loss/He3_CD2.G4table","G4Table",100 ); + He3Al = EnergyLoss("energy_loss/He3_Al.G4table","G4Table",10); + He3Si = EnergyLoss("energy_loss/He3_Si.G4table","G4Table",10); + Li11CD2 = EnergyLoss("energy_loss/Li11_CD2.G4table","G4Table",100); } //////////////////////////////////////////////////////////////////////////////// diff --git a/Examples/Example1/ShowResults.C b/Examples/Example1/ShowResults.C index e398a21e2..6d7c2ed41 100644 --- a/Examples/Example1/ShowResults.C +++ b/Examples/Example1/ShowResults.C @@ -16,7 +16,7 @@ void LoadCuts(){ //////////////////////////////////////////////////////////////////////////////// void LoadChain(){ chain = new TChain("PhysicsTree"); - chain->Add("../../Outputs/Analysis/Example1.root"); + chain->Add("root/analysis/Example1.root"); } //////////////////////////////////////////////////////////////////////////////// diff --git a/Inputs/EventGenerator/np.source b/Inputs/EventGenerator/np.source index 59867b6ea..79d252791 100644 --- a/Inputs/EventGenerator/np.source +++ b/Inputs/EventGenerator/np.source @@ -11,7 +11,7 @@ Isotropic x0= 0 mm y0= 0 mm z0= 0 mm - Multiplicity= 1 1 - Particle= neutron proton + Multiplicity= 1 1 1 + Particle= neutron proton 8He %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Supported particle type: proton, neutron, deuton, triton, He3 , alpha diff --git a/NPLib/Detectors/Minos/TMinosPhysics.cxx b/NPLib/Detectors/Minos/TMinosPhysics.cxx index ad776262f..cc2e729e7 100644 --- a/NPLib/Detectors/Minos/TMinosPhysics.cxx +++ b/NPLib/Detectors/Minos/TMinosPhysics.cxx @@ -101,7 +101,7 @@ void TMinosPhysics::PreTreat() { vector<unsigned short>* Time = m_EventData->GetTimePtr(i); m_utility.Calibrate(Time,Charge,i,T,Q); - if(T>0){ + if(T>0&&Q<77000){ PadNumber = m_EventData->GetPadNumber(i); double x_mm = m_X[PadNumber]; double y_mm = m_Y[PadNumber]; @@ -111,16 +111,12 @@ void TMinosPhysics::PreTreat() { double calV= cal->GetValue(cal_v,0); double calO=cal->GetValue(cal_o,0); - // regular case - //double z_mm = (T*m_TimeBin+cal->GetValue(cal_o,0))*cal->GetValue(cal_v,0); - - // testing dependency: - //calV+=-calV*30./100.;// -5% - //calO+=+calO*60./100; - double z_mm = (T*m_TimeBin+calO)*calV; + double z_mm = (T*m_TimeBin+calO/*-0.7-0.04*ring*/)*calV; + //cout << T*m_TimeBin+calO << endl; TVector3 Pos=TVector3(x_mm+m_Position.X(),y_mm+m_Position.Y(),z_mm+m_Position.Z()); Pos.RotateZ(m_ZRotation); // Calibrate the Pad: + X_Pad.push_back(Pos.X()); Y_Pad.push_back(Pos.Y()); Z_Pad.push_back(Pos.Z()); diff --git a/NPLib/Detectors/SuperX3/CMakeLists.txt b/NPLib/Detectors/SuperX3/CMakeLists.txt new file mode 100644 index 000000000..9615752a2 --- /dev/null +++ b/NPLib/Detectors/SuperX3/CMakeLists.txt @@ -0,0 +1,6 @@ +add_custom_command(OUTPUT TSuperX3PhysicsDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TSuperX3Physics.h TSuperX3PhysicsDict.cxx TSuperX3Physics.rootmap libNPSuperX3.dylib DEPENDS TSuperX3Physics.h) +add_custom_command(OUTPUT TSuperX3DataDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TSuperX3Data.h TSuperX3DataDict.cxx TSuperX3Data.rootmap libNPSuperX3.dylib DEPENDS TSuperX3Data.h) +add_library(NPSuperX3 SHARED TSuperX3Spectra.cxx TSuperX3Data.cxx TSuperX3Physics.cxx TSuperX3DataDict.cxx TSuperX3PhysicsDict.cxx ) +target_link_libraries(NPSuperX3 ${ROOT_LIBRARIES} NPCore) +install(FILES TSuperX3Data.h TSuperX3Physics.h TSuperX3Spectra.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY}) + diff --git a/NPLib/Detectors/SuperX3/TSuperX3Data.cxx b/NPLib/Detectors/SuperX3/TSuperX3Data.cxx new file mode 100644 index 000000000..15bed70e3 --- /dev/null +++ b/NPLib/Detectors/SuperX3/TSuperX3Data.cxx @@ -0,0 +1,62 @@ +/***************************************************************************** + * Copyright (C) 2009-2016 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: N. de Sereville contact address: deserevi@ipno.in2p3.fr * + * * + * Creation Date : january 2011 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: * + * This class holds the raw data storage for the SuperX3 detector from Micron * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ +#include "TSuperX3Data.h" + +#include <iostream> +using namespace std; + +ClassImp(TSuperX3Data) + + TSuperX3Data::TSuperX3Data() { + // Default constructor + Clear(); +} + +TSuperX3Data::~TSuperX3Data() {} + +void TSuperX3Data::Clear() { + // (Up, E) + fSuperX3_UpE_DetectorNbr.clear(); + fSuperX3_UpE_StripNbr.clear(); + fSuperX3_UpE_Energy.clear(); + // (Up, T) + fSuperX3_UpT_DetectorNbr.clear(); + fSuperX3_UpT_StripNbr.clear(); + fSuperX3_UpT_Time.clear(); + // (Down, E) + fSuperX3_DownE_DetectorNbr.clear(); + fSuperX3_DownE_StripNbr.clear(); + fSuperX3_DownE_Energy.clear(); + // (Down, T) + fSuperX3_DownT_DetectorNbr.clear(); + fSuperX3_DownT_StripNbr.clear(); + fSuperX3_DownT_Time.clear(); + + // (Back, E) + fSuperX3_BackE_DetectorNbr.clear(); + fSuperX3_BackE_Energy.clear(); + // (Back, T) + fSuperX3_BackT_DetectorNbr.clear(); + fSuperX3_BackT_Time.clear(); +} + +void TSuperX3Data::Dump() const {} diff --git a/NPLib/Detectors/SuperX3/TSuperX3Data.h b/NPLib/Detectors/SuperX3/TSuperX3Data.h new file mode 100644 index 000000000..b4cd61218 --- /dev/null +++ b/NPLib/Detectors/SuperX3/TSuperX3Data.h @@ -0,0 +1,138 @@ +#ifndef __SuperX3DATA__ +#define __SuperX3DATA__ +/***************************************************************************** + * Copyright (C) 2009-2016 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: N. de Sereville contact address: deserevi@ipno.in2p3.fr * + * * + * Creation Date : january 2011 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: * + * This class holds the raw data storage for the SuperX3 detector from Micron * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ +// C++ headers +#include <vector> +using namespace std; + +// ROOT headers +#include "TObject.h" + +class TSuperX3Data : public TObject { + private: + // Front + // Up + // Energy + vector<UShort_t> fSuperX3_UpE_DetectorNbr; + vector<UShort_t> fSuperX3_UpE_StripNbr; + vector<Double_t> fSuperX3_UpE_Energy; + // Time + vector<UShort_t> fSuperX3_UpT_DetectorNbr; + vector<UShort_t> fSuperX3_UpT_StripNbr; + vector<Double_t> fSuperX3_UpT_Time; + // Energy + vector<UShort_t> fSuperX3_DownE_DetectorNbr; + vector<UShort_t> fSuperX3_DownE_StripNbr; + vector<Double_t> fSuperX3_DownE_Energy; + // Time + vector<UShort_t> fSuperX3_DownT_DetectorNbr; + vector<UShort_t> fSuperX3_DownT_StripNbr; + vector<Double_t> fSuperX3_DownT_Time; + + // Back + // Energy + vector<Double_t> fSuperX3_BackE_Energy; + vector<UShort_t> fSuperX3_BackE_DetectorNbr; + // Time + vector<Double_t> fSuperX3_BackT_Time; + vector<UShort_t> fSuperX3_BackT_DetectorNbr; + + public: + TSuperX3Data(); + virtual ~TSuperX3Data(); + + void Clear(); + void Clear(const Option_t*){}; + void Dump() const; + + ///////////////////// SETTERS //////////////////////// + // up + void SetUpE(UShort_t DetNbr, UShort_t StripNbr, Double_t Energy) { + fSuperX3_UpE_DetectorNbr.push_back(DetNbr); + fSuperX3_UpE_StripNbr.push_back(StripNbr); + fSuperX3_UpE_Energy.push_back(Energy); + } + void SetUpT(UShort_t DetNbr, UShort_t StripNbr, Double_t Time) { + fSuperX3_UpT_DetectorNbr.push_back(DetNbr); + fSuperX3_UpT_StripNbr.push_back(StripNbr); + fSuperX3_UpT_Time.push_back(Time); + } + // down + void SetDownE(UShort_t DetNbr, UShort_t StripNbr, Double_t Energy) { + fSuperX3_DownE_DetectorNbr.push_back(DetNbr); + fSuperX3_DownE_StripNbr.push_back(StripNbr); + fSuperX3_DownE_Energy.push_back(Energy); + } + void SetDownT(UShort_t DetNbr, UShort_t StripNbr, Double_t Time) { + fSuperX3_DownT_DetectorNbr.push_back(DetNbr); + fSuperX3_DownT_StripNbr.push_back(StripNbr); + fSuperX3_DownT_Time.push_back(Time); + } + + // Back E + void SetBackE(UShort_t DetNbr, Double_t Energy) { + fSuperX3_BackE_DetectorNbr.push_back(DetNbr); + fSuperX3_BackE_Energy.push_back(Energy); + } + // Back T + void SetBackT(UShort_t DetNbr, Double_t Time) { + fSuperX3_BackT_DetectorNbr.push_back(DetNbr); + fSuperX3_BackT_Time.push_back(Time); + } + + ///////////////////// GETTERS //////////////////////// + // DSSD + // (Up, E) + UShort_t GetUpEMult() const { return fSuperX3_UpE_DetectorNbr.size(); } + UShort_t GetUpEDetectorNbr(Int_t i) const { return fSuperX3_UpE_DetectorNbr[i]; } + UShort_t GetUpEStripNbr(Int_t i) const { return fSuperX3_UpE_StripNbr[i]; } + Double_t GetUpEEnergy(Int_t i) const { return fSuperX3_UpE_Energy[i]; } + // (Up, T) + UShort_t GetUpTMult() const { return fSuperX3_UpT_DetectorNbr.size(); } + UShort_t GetUpTDetectorNbr(Int_t i) const { return fSuperX3_UpT_DetectorNbr[i]; } + UShort_t GetUpTStripNbr(Int_t i) const { return fSuperX3_UpT_StripNbr[i]; } + Double_t GetUpTTime(Int_t i) const { return fSuperX3_UpT_Time[i]; } + // (Down, E) + UShort_t GetDownEMult() const { return fSuperX3_DownE_DetectorNbr.size(); } + UShort_t GetDownEDetectorNbr(Int_t i) const { return fSuperX3_DownE_DetectorNbr[i]; } + UShort_t GetDownEStripNbr(Int_t i) const { return fSuperX3_DownE_StripNbr[i]; } + Double_t GetDownEEnergy(Int_t i) const { return fSuperX3_DownE_Energy[i]; } + // (Down, T) + UShort_t GetDownTMult() const { return fSuperX3_DownT_DetectorNbr.size(); } + UShort_t GetDownTDetectorNbr(Int_t i) const { return fSuperX3_DownT_DetectorNbr[i]; } + UShort_t GetDownTStripNbr(Int_t i) const { return fSuperX3_DownT_StripNbr[i]; } + Double_t GetDownTTime(Int_t i) const { return fSuperX3_DownT_Time[i]; } + + // (Back, E) + UShort_t GetBackEMult() const { return fSuperX3_BackE_DetectorNbr.size(); } + UShort_t GetBackEDetectorNbr(Int_t i) const { return fSuperX3_BackE_DetectorNbr[i]; } + Double_t GetBackEEnergy(Int_t i) const { return fSuperX3_BackE_Energy[i]; } + // (Back, T) + UShort_t GetBackTMult() const { return fSuperX3_BackT_DetectorNbr.size(); } + UShort_t GetBackTDetectorNbr(Int_t i) const { return fSuperX3_BackT_DetectorNbr[i]; } + Double_t GetBackTTime(Int_t i) const { return fSuperX3_BackT_Time[i]; } + + ClassDef(TSuperX3Data, 1) // TSuperX3Data raw data +}; + +#endif diff --git a/NPLib/Detectors/SuperX3/TSuperX3Physics.cxx b/NPLib/Detectors/SuperX3/TSuperX3Physics.cxx new file mode 100644 index 000000000..0d21f4073 --- /dev/null +++ b/NPLib/Detectors/SuperX3/TSuperX3Physics.cxx @@ -0,0 +1,698 @@ +/***************************************************************************** + * Copyright (C) 2009-2016 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: Adrien MATTA contact address: matta@lpccaen.in2p3.fr * + * * + * Creation Date : november 2009 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: * + * This class hold SuperX3 Physics * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ +// NPL +#include "TSuperX3Physics.h" +using namespace SuperX3_LOCAL; + +// STL +#include <cmath> +#include <fstream> +#include <iostream> +#include <limits> +#include <sstream> +#include <stdlib.h> +using namespace std; + +// NPTOOL +#include "NPDetectorFactory.h" +#include "NPOptionManager.h" +#include "RootInput.h" +#include "RootOutput.h" + +// ROOT +#include "TChain.h" + +ClassImp(TSuperX3Physics) + + /////////////////////////////////////////////////////////////////////////// + TSuperX3Physics::TSuperX3Physics() + : m_EventData(new TSuperX3Data), m_PreTreatedData(new TSuperX3Data), m_EventPhysics(this), m_Spectra(0), + m_nCounter(10), m_MaximumStripMultiplicityAllowed(1), // multiplidity 1 + m_StripEnergyMatchingSigma(0.060), // MeV + m_StripEnergyMatchingNumberOfSigma(5), // MeV + m_FrontE_Raw_Threshold(0), // adc channels + m_BackE_Raw_Threshold(0), // adc channels + m_FrontE_Calib_Threshold(0), // MeV + m_BackE_Calib_Threshold(0), // MeV + m_NumberOfDetectors(0), m_SiliconFace(49.6), // mm + m_NumberOfStrips(16) { + m_StripPitch = m_SiliconFace / (double)m_NumberOfStrips; + + for (Int_t i = 0; i < m_nCounter; ++i) { // loop on counters + m_Counter[i] = 0; + } // end loop on counters +} + +/////////////////////////////////////////////////////////////////////////// +// TSuperX3Physics::~TSuperX3Physics() +//{ +// delete m_EventData; +// delete m_PreTreatedData; +//} + +/////////////////////////////////////////////////////////////////////////// +void TSuperX3Physics::Clear() { + fEventType.clear(); + fDetectorNumber.clear(); + fFrontEnergy.clear(); + fBackEnergy.clear(); + fHalfEnergy.clear(); + fFrontTime.clear(); + fBackTime.clear(); + fFrontStrip.clear(); + fBackStrip.clear(); + + for (Int_t i = 0; i < m_nCounter; ++i) { // loop on counters + m_Counter[i] = 0; + } // end loop on counters +} + +/////////////////////////////////////////////////////////////////////////// +void TSuperX3Physics::ReadConfiguration(NPL::InputParser parser) { + vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("SuperX3"); + if (NPOptionManager::getInstance()->GetVerboseLevel()) + cout << "//// " << blocks.size() << " detectors found " << endl; + for (unsigned int i = 0; i < blocks.size(); i++) { + // Cartesian Case + vector<string> cart = {"X1_Y1", "X1_Y16", "X16_Y1", "X16_Y16"}; + // Spherical Case + vector<string> sphe = {"R", "THETA", "PHI", "BETA"}; + + if (blocks[i]->HasTokenList(cart)) { + cout << endl << "//// SuperX3 " << i + 1 << endl; + TVector3 A = blocks[i]->GetTVector3("X1_Y1", "mm"); + TVector3 B = blocks[i]->GetTVector3("X16_Y1", "mm"); + TVector3 C = blocks[i]->GetTVector3("X1_Y16", "mm"); + TVector3 D = blocks[i]->GetTVector3("X16_Y16", "mm"); + AddDetector(A, B, C, D); + } + + else if (blocks[i]->HasTokenList(sphe)) { + double Theta = blocks[i]->GetDouble("THETA", "deg"); + double Phi = blocks[i]->GetDouble("PHI", "deg"); + double R = blocks[i]->GetDouble("R", "mm"); + vector<double> beta = blocks[i]->GetVectorDouble("BETA", "deg"); + AddDetector(Theta, Phi, R, beta[0], beta[1], beta[2]); + } + + else { + cout << "ERROR: Missing token for SuperX3 blocks, check your input file" << endl; + exit(1); + } + } + + InitializeStandardParameters(); + ReadAnalysisConfig(); +} + +/////////////////////////////////////////////////////////////////////////// +void TSuperX3Physics::AddParameterToCalibrationManager() { + CalibrationManager* Cal = CalibrationManager::getInstance(); + + for (Int_t i = 0; i < m_NumberOfDetectors; i++) { + for (Int_t j = 0; j < m_NumberOfStrips; j++) { + // Energy + Cal->AddParameter("SuperX3", "D" + NPL::itoa(i + 1) + "_Front" + NPL::itoa(j) + "_E", + "SuperX3_D" + NPL::itoa(i + 1) + "_FRONT" + NPL::itoa(j) + "_E"); + Cal->AddParameter("SuperX3", "D" + NPL::itoa(i + 1) + "_Back" + NPL::itoa(j) + "_E", + "SuperX3_D" + NPL::itoa(i + 1) + "_BACK" + NPL::itoa(j) + "_E"); + // Time + Cal->AddParameter("SuperX3", "D" + NPL::itoa(i + 1) + "_Front" + NPL::itoa(j) + "_T", + "SuperX3_D" + NPL::itoa(i + 1) + "_FRONT" + NPL::itoa(j) + "_T"); + Cal->AddParameter("SuperX3", "D" + NPL::itoa(i + 1) + "_Back" + NPL::itoa(j) + "_T", + "SuperX3_D" + NPL::itoa(i + 1) + "_BACK" + NPL::itoa(j) + "_T"); + } + } +} + +/////////////////////////////////////////////////////////////////////////// +void TSuperX3Physics::InitializeRootInputRaw() { + TChain* inputChain = RootInput::getInstance()->GetChain(); + inputChain->SetBranchStatus("SuperX3", true); + inputChain->SetBranchStatus("fSuperX3_*", true); + inputChain->SetBranchAddress("SuperX3", &m_EventData); +} + +/////////////////////////////////////////////////////////////////////////// +void TSuperX3Physics::InitializeRootInputPhysics() { + TChain* inputChain = RootInput::getInstance()->GetChain(); + inputChain->SetBranchStatus("SuperX3", true); + inputChain->SetBranchStatus("fEventType", true); + inputChain->SetBranchStatus("fDetectorNumber", true); + inputChain->SetBranchStatus("fEnergy", true); + inputChain->SetBranchStatus("fTime", true); + inputChain->SetBranchStatus("fFrontStrip", true); + inputChain->SetBranchStatus("fBackStrip", true); + inputChain->SetBranchAddress("SuperX3", &m_EventPhysics); +} + +/////////////////////////////////////////////////////////////////////////// +void TSuperX3Physics::InitializeRootOutput() { + TTree* outputTree = RootOutput::getInstance()->GetTree(); + outputTree->Branch("SuperX3", "TSuperX3Physics", &m_EventPhysics); +} + +void TSuperX3Physics::AddDetector(TVector3 C_X1_Y1, TVector3 C_X16_Y1, TVector3 C_X1_Y16, TVector3 C_X16_Y16) { + m_NumberOfDetectors++; + + // remove warning using C_X16_Y16 + C_X16_Y16.Unit(); + + // Vector U on Module Face (paralelle to Y Strip) (NB: remember that Y strip are allong X axis) + TVector3 U = C_X16_Y1 - C_X1_Y1; + U = U.Unit(); + + // Vector V on Module Face (parallele to X Strip) + TVector3 V = C_X1_Y16 - C_X1_Y1; + V = V.Unit(); + + // Position Vector of Strip Center + TVector3 StripCenter = TVector3(0, 0, 0); + // Position Vector of X=1 Y=1 Strip + TVector3 Strip_1_1; + + // Buffer object to fill Position Array + vector<double> lineX; + vector<double> lineY; + vector<double> lineZ; + + vector<vector<double>> OneModuleStripPositionX; + vector<vector<double>> OneModuleStripPositionY; + vector<vector<double>> OneModuleStripPositionZ; + + // Moving StripCenter to 1.1 corner: + Strip_1_1 = C_X1_Y1 + (U + V) * (m_StripPitch / 2.); + + for (int i = 0; i < m_NumberOfStrips; i++) { + lineX.clear(); + lineY.clear(); + lineZ.clear(); + + for (int j = 0; j < m_NumberOfStrips; j++) { + StripCenter = Strip_1_1 + m_StripPitch * (i * U + j * V); + + lineX.push_back(StripCenter.X()); + lineY.push_back(StripCenter.Y()); + lineZ.push_back(StripCenter.Z()); + } + + OneModuleStripPositionX.push_back(lineX); + OneModuleStripPositionY.push_back(lineY); + OneModuleStripPositionZ.push_back(lineZ); + } + + m_StripPositionX.push_back(OneModuleStripPositionX); + m_StripPositionY.push_back(OneModuleStripPositionY); + m_StripPositionZ.push_back(OneModuleStripPositionZ); +} + +void TSuperX3Physics::AddDetector(double theta, double phi, double distance, double beta_u, double beta_v, + double beta_w) { + m_NumberOfDetectors++; + + // convert from degree to radian: + theta *= M_PI / 180; + phi *= M_PI / 180; + + // Vector U on Module Face (paralelle to Y Strip) (NB: remember that Y strip are allong X axis) + TVector3 U; + // Vector V on Module Face (parallele to X Strip) + TVector3 V; + // Vector W normal to Module Face (pointing CsI) + TVector3 W; + // Vector position of Module Face center + TVector3 C; + + C = TVector3(distance * sin(theta) * cos(phi), distance * sin(theta) * sin(phi), distance * cos(theta)); + + TVector3 YperpW = TVector3(cos(theta) * cos(phi), cos(theta) * sin(phi), -sin(theta)); + + W = C.Unit(); + U = W.Cross(YperpW); + V = W.Cross(U); + + U = U.Unit(); + V = V.Unit(); + + U.Rotate(beta_u * M_PI / 180, U); + V.Rotate(beta_u * M_PI / 180, U); + + U.Rotate(beta_v * M_PI / 180, V); + V.Rotate(beta_v * M_PI / 180, V); + + U.Rotate(beta_w * M_PI / 180, W); + V.Rotate(beta_w * M_PI / 180, W); + + // Position Vector of Strip Center + TVector3 StripCenter = TVector3(0, 0, 0); + // Position Vector of X=1 Y=1 Strip + TVector3 Strip_1_1; + + vector<double> lineX; + vector<double> lineY; + vector<double> lineZ; + + vector<vector<double>> OneModuleStripPositionX; + vector<vector<double>> OneModuleStripPositionY; + vector<vector<double>> OneModuleStripPositionZ; + + // Moving C to the 1.1 corner: + Strip_1_1 = C - 0.5 * (m_SiliconFace - m_StripPitch) * (U + V); + + for (int i = 0; i < m_NumberOfStrips; i++) { + lineX.clear(); + lineY.clear(); + lineZ.clear(); + + for (int j = 0; j < m_NumberOfStrips; j++) { + StripCenter = Strip_1_1 + m_StripPitch * (i * U + j * V); + + lineX.push_back(StripCenter.X()); + lineY.push_back(StripCenter.Y()); + lineZ.push_back(StripCenter.Z()); + } + + OneModuleStripPositionX.push_back(lineX); + OneModuleStripPositionY.push_back(lineY); + OneModuleStripPositionZ.push_back(lineZ); + } + + m_StripPositionX.push_back(OneModuleStripPositionX); + m_StripPositionY.push_back(OneModuleStripPositionY); + m_StripPositionZ.push_back(OneModuleStripPositionZ); +} + +TVector3 TSuperX3Physics::GetPositionOfInteraction(int i) { + TVector3 Position = TVector3(GetStripPositionX(fDetectorNumber[i], fFrontStrip[i], fBackStrip[i]), + GetStripPositionY(fDetectorNumber[i], fFrontStrip[i], fBackStrip[i]), + GetStripPositionZ(fDetectorNumber[i], fFrontStrip[i], fBackStrip[i])); + + return Position; +} + +TVector3 TSuperX3Physics::GetDetectorNormal(int i) { + TVector3 U = TVector3(GetStripPositionX(fDetectorNumber[i], m_NumberOfStrips, 1), + GetStripPositionY(fDetectorNumber[i], m_NumberOfStrips, 1), + GetStripPositionZ(fDetectorNumber[i], m_NumberOfStrips, 1)) + + - TVector3(GetStripPositionX(fDetectorNumber[i], 1, 1), GetStripPositionY(fDetectorNumber[i], 1, 1), + GetStripPositionZ(fDetectorNumber[i], 1, 1)); + + TVector3 V = TVector3(GetStripPositionX(fDetectorNumber[i], m_NumberOfStrips, m_NumberOfStrips), + GetStripPositionY(fDetectorNumber[i], m_NumberOfStrips, m_NumberOfStrips), + GetStripPositionZ(fDetectorNumber[i], m_NumberOfStrips, m_NumberOfStrips)) + + - TVector3(GetStripPositionX(fDetectorNumber[i], m_NumberOfStrips, 1), + GetStripPositionY(fDetectorNumber[i], m_NumberOfStrips, 1), + GetStripPositionZ(fDetectorNumber[i], m_NumberOfStrips, 1)); + + TVector3 Normal = U.Cross(V); + + return Normal.Unit(); +} + +void TSuperX3Physics::DumpStrippingScheme(Int_t detecNumber) { + cout << endl << "TSuperX3Physics::DumpStrippingScheme()" << endl; + cout << "Detector number " << detecNumber << endl; + + for (int i = 1; i < m_NumberOfStrips + 1; i++) { // front part + for (int j = 1; j < m_NumberOfStrips + 1; j++) { // back part + cout << "strips Front, Back: " << i << " " << j + << "\t--->\t (X,Y,Z) mm: " << GetStripPositionX(detecNumber, i, j) << "\t" + << GetStripPositionY(detecNumber, i, j) << "\t" << GetStripPositionZ(detecNumber, i, j) << endl; + } + } +} + +/////////////////////////////////////////////////////////////////////////// +void TSuperX3Physics::BuildPhysicalEvent() { BuildSimplePhysicalEvent(); } + +/////////////////////////////////////////////////////////////////////////// +void TSuperX3Physics::BuildSimplePhysicalEvent() { + // Select active channels and apply thresholds + PreTreat(); + /* m_Counter[0] = 1; + + // Begin treatement + Int_t evtType = EventType(); + + if (evtType == 1) { // case where multiplicity front = multiplicity back + m_Counter[1] = 1; + vector<TVector2> couple = Match_Front_Back(); + + for (UShort_t i = 0; i < couple.size(); i++) { // loop on selected events + Int_t DetecNbr = m_PreTreatedData->GetFrontEDetectorNbr(couple[i].X()); + Int_t StripFront = m_PreTreatedData->GetFrontEStripNbr(couple[i].X()); + Int_t StripBack = m_PreTreatedData->GetBackEStripNbr(couple[i].Y()); + Double_t EnergyFront = m_PreTreatedData->GetFrontEEnergy(couple[i].X()); + Double_t EnergyBack = m_PreTreatedData->GetBackEEnergy(couple[i].Y()); + + // Search for associate time + // Front + Double_t TimeFront = -1000; + for (UShort_t t = 0; t < m_PreTreatedData->GetFrontTMult(); t++) { + if (m_PreTreatedData->GetFrontTStripNbr(couple[i].X()) == m_PreTreatedData->GetFrontTStripNbr(t) && + m_PreTreatedData->GetFrontTDetectorNbr(couple[i].X()) == m_PreTreatedData->GetFrontTDetectorNbr(t)) { + TimeFront = m_PreTreatedData->GetFrontTTime(t); + m_Counter[4] = 1; + } + } + // Back + Double_t TimeBack = -1000; + for (UShort_t t = 0; t < m_PreTreatedData->GetBackTMult(); t++) { + if (m_PreTreatedData->GetBackTStripNbr(couple[i].Y()) == m_PreTreatedData->GetBackTStripNbr(t) && + m_PreTreatedData->GetBackTDetectorNbr(couple[i].Y()) == m_PreTreatedData->GetBackTDetectorNbr(t)) + TimeBack = m_PreTreatedData->GetBackTTime(t); + } + + // Fill TSuperX3Physics private members + fEventType.push_back(evtType); + fDetectorNumber.push_back(DetecNbr); + fFrontEnergy.push_back(EnergyFront); + fBackEnergy.push_back(EnergyBack); + fHalfEnergy.push_back((EnergyFront + EnergyBack) / 2); + fFrontTime.push_back(TimeFront); + fBackTime.push_back(TimeBack); + fFrontStrip.push_back(StripFront); + fBackStrip.push_back(StripBack); + } + }*/ +} + +/////////////////////////////////////////////////////////////////////////// +void TSuperX3Physics::PreTreat() { + // Clear pre treated object + ClearPreTreatedData(); + /* + // (Front, E) + for (UShort_t i = 0; i < m_EventData->GetFrontEMult(); i++) { + if (IsValidChannel("Front", m_EventData->GetFrontEDetectorNbr(i), m_EventData->GetFrontEStripNbr(i)) && + m_EventData->GetFrontEEnergy(i) > m_FrontE_Raw_Threshold) { + Double_t E = fSuperX3_Front_E(m_EventData, i); + if (E > m_FrontE_Calib_Threshold) { + m_PreTreatedData->SetFrontEDetectorNbr(m_EventData->GetFrontEDetectorNbr(i)); + m_PreTreatedData->SetFrontEStripNbr(m_EventData->GetFrontEStripNbr(i)); + m_PreTreatedData->SetFrontEEnergy(E); + } + } + } + // (Front, T) + for (UShort_t i = 0; i < m_EventData->GetFrontTMult(); i++) { + if (IsValidChannel("Front", m_EventData->GetFrontTDetectorNbr(i), m_EventData->GetFrontTStripNbr(i))) { + Double_t T = fSuperX3_Front_T(m_EventData, i); + m_PreTreatedData->SetFrontTDetectorNbr(m_EventData->GetFrontTDetectorNbr(i)); + m_PreTreatedData->SetFrontTStripNbr(m_EventData->GetFrontTStripNbr(i)); + m_PreTreatedData->SetFrontTTime(T); + } + } + + // (Back, E) + for (UShort_t i = 0; i < m_EventData->GetBackEMult(); i++) { + if (IsValidChannel("Back", m_EventData->GetBackEDetectorNbr(i), m_EventData->GetBackEStripNbr(i)) && + m_EventData->GetBackEEnergy(i) > m_BackE_Raw_Threshold) { + Double_t E = fSuperX3_Back_E(m_EventData, i); + if (E > m_BackE_Calib_Threshold) { + m_PreTreatedData->SetBackEDetectorNbr(m_EventData->GetBackEDetectorNbr(i)); + m_PreTreatedData->SetBackEStripNbr(m_EventData->GetBackEStripNbr(i)); + m_PreTreatedData->SetBackEEnergy(E); + } + } + } + // (Back, T) + for (UShort_t i = 0; i < m_EventData->GetBackTMult(); i++) { + if (IsValidChannel("Back", m_EventData->GetBackTDetectorNbr(i), m_EventData->GetBackTStripNbr(i))) { + Double_t T = fSuperX3_Back_T(m_EventData, i); + m_PreTreatedData->SetBackTDetectorNbr(m_EventData->GetBackTDetectorNbr(i)); + m_PreTreatedData->SetBackTStripNbr(m_EventData->GetBackTStripNbr(i)); + m_PreTreatedData->SetBackTTime(T); + } + }*/ +} + +/////////////////////////////////////////////////////////////////////////// +Int_t TSuperX3Physics::EventType() { + /* // Same multiplicity on front and back side + if (m_PreTreatedData->GetFrontEMult() == m_PreTreatedData->GetBackEMult()) { + return 1; + } + // Possibly interstrip + else if (m_PreTreatedData->GetFrontEMult() == m_PreTreatedData->GetBackEMult() + 1 || + m_PreTreatedData->GetFrontEMult() == m_PreTreatedData->GetBackEMult() - 1) { + return 2; + } + // Rejected event + else { + return -1; + }*/ + return 0; +} + +/////////////////////////////////////////////////////////////////////////// +vector<TVector2> TSuperX3Physics::Match_Front_Back() { + vector<TVector2> ArrayOfGoodCouple; + /* + // Select allowed multiplicity events. If multiplicity is too + // high, then return "empty" vector + if (m_PreTreatedData->GetFrontEMult() > m_MaximumStripMultiplicityAllowed || + m_PreTreatedData->GetBackEMult() > m_MaximumStripMultiplicityAllowed) + return ArrayOfGoodCouple; + + // Loop on front multiplicity + for (UShort_t i = 0; i < m_PreTreatedData->GetFrontEMult(); i++) { + // Loop on back multiplicity + for (UShort_t j = 0; j < m_PreTreatedData->GetBackEMult(); j++) { + // if same detector check energy + if (m_PreTreatedData->GetFrontEDetectorNbr(i) == m_PreTreatedData->GetBackEDetectorNbr(j)) { + m_Counter[2] = 1; + // Equal energy + if (abs((m_PreTreatedData->GetFrontEEnergy(i) - m_PreTreatedData->GetBackEEnergy(j)) / 2.) < + m_StripEnergyMatchingNumberOfSigma * m_StripEnergyMatchingSigma) { + ArrayOfGoodCouple.push_back(TVector2(i, j)); + m_Counter[3] = 1; + } // end test energy + } // end test same detector + } // end loop back multiplicity + } // end loop front multiplicity + + // Prevent treating event with ambiguous matchin beetween X and Y + if (ArrayOfGoodCouple.size() > m_PreTreatedData->GetFrontEMult()) + ArrayOfGoodCouple.clear(); + */ + return ArrayOfGoodCouple; +} + +/////////////////////////////////////////////////////////////////////////// +bool TSuperX3Physics::IsValidChannel(string Type, int detector, int channel) { + if (Type == "Front") { + return *(m_FrontChannelStatus[detector - 1].begin() + channel); + } + else if (Type == "Back") + return *(m_BackChannelStatus[detector - 1].begin() + channel); + + else + return false; +} + +/////////////////////////////////////////////////////////////////////////// +void TSuperX3Physics::InitializeStandardParameters() { + // Enable all channels + vector<bool> ChannelStatus; + m_FrontChannelStatus.clear(); + m_BackChannelStatus.clear(); + + ChannelStatus.resize(m_NumberOfStrips, true); + for (Int_t i = 0; i < m_NumberOfDetectors; i++) { + m_FrontChannelStatus[i] = ChannelStatus; + m_BackChannelStatus[i] = ChannelStatus; + } +} + +/////////////////////////////////////////////////////////////////////////// +void TSuperX3Physics::ReadAnalysisConfig() { + bool ReadingStatus = false; + + cout << "\t/////////// Reading ConfigSuperX3.dat file ///////////" << endl; + + // path to file + string FileName = "./configs/ConfigSuperX3.dat"; + + // open analysis config file + ifstream AnalysisConfigFile; + AnalysisConfigFile.open(FileName.c_str()); + + if (!AnalysisConfigFile.is_open()) { + cout << "\tNo ConfigSuperX3.dat found: default parameters loaded for Analysis " << FileName << endl; + return; + } + cout << "\tLoading user parameters from ConfigSuperX3.dat " << endl; + + // storing config file in the ROOT output file + TAsciiFile* asciiFile = RootOutput::getInstance()->GetAsciiFileAnalysisConfig(); + asciiFile->AppendLine("%% ConfigSuperX3.dat %%"); + asciiFile->Append(FileName.c_str()); + asciiFile->AppendLine(""); + + // read analysis config file + string LineBuffer, DataBuffer, whatToDo; + while (!AnalysisConfigFile.eof()) { + // Pick-up next line + getline(AnalysisConfigFile, LineBuffer); + + // search for "header" + if (LineBuffer.compare(0, 8, "ConfigSuperX3") == 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 == "MAX_STRIP_MULTIPLICITY") { + AnalysisConfigFile >> DataBuffer; + m_MaximumStripMultiplicityAllowed = atoi(DataBuffer.c_str()); + cout << "\t" << whatToDo << "\t" << m_MaximumStripMultiplicityAllowed << endl; + } + + else if (whatToDo == "FRONT_BACK_ENERGY_MATCHING_SIGMA") { + AnalysisConfigFile >> DataBuffer; + m_StripEnergyMatchingSigma = atof(DataBuffer.c_str()); + cout << "\t" << whatToDo << "\t" << m_StripEnergyMatchingSigma << endl; + } + + else if (whatToDo == "FRONT_BACK_ENERGY_MATCHING_NUMBER_OF_SIGMA") { + AnalysisConfigFile >> DataBuffer; + m_StripEnergyMatchingNumberOfSigma = atoi(DataBuffer.c_str()); + cout << "\t" << whatToDo << "\t" << m_StripEnergyMatchingNumberOfSigma << endl; + } + + else if (whatToDo == "FRONT_E_RAW_THRESHOLD") { + AnalysisConfigFile >> DataBuffer; + m_FrontE_Raw_Threshold = atoi(DataBuffer.c_str()); + cout << "\t" << whatToDo << "\t" << m_FrontE_Raw_Threshold << endl; + } + + else if (whatToDo == "BACK_E_RAW_THRESHOLD") { + AnalysisConfigFile >> DataBuffer; + m_BackE_Raw_Threshold = atoi(DataBuffer.c_str()); + cout << "\t" << whatToDo << "\t" << m_BackE_Raw_Threshold << endl; + } + + else if (whatToDo == "FRONT_E_CAL_THRESHOLD") { + AnalysisConfigFile >> DataBuffer; + m_FrontE_Calib_Threshold = atoi(DataBuffer.c_str()); + cout << "\t" << whatToDo << "\t" << m_FrontE_Calib_Threshold << endl; + } + + else if (whatToDo == "BACK_E_CAL_THRESHOLD") { + AnalysisConfigFile >> DataBuffer; + m_BackE_Calib_Threshold = atoi(DataBuffer.c_str()); + cout << "\t" << whatToDo << "\t" << m_BackE_Calib_Threshold << endl; + } + + else if (whatToDo == "DISABLE_ALL") { + AnalysisConfigFile >> DataBuffer; + cout << "\t" << whatToDo << "\t" << DataBuffer << endl; + Int_t Detector = atoi(DataBuffer.substr(4, 1).c_str()); + vector<bool> ChannelStatus; + ChannelStatus.resize(m_NumberOfStrips, false); + m_FrontChannelStatus[Detector - 1] = ChannelStatus; + m_BackChannelStatus[Detector - 1] = ChannelStatus; + } + + else if (whatToDo == "DISABLE_CHANNEL") { + AnalysisConfigFile >> DataBuffer; + cout << "\t" << whatToDo << "\t" << DataBuffer << endl; + int detector = atoi(DataBuffer.substr(4, 1).c_str()); + int channel = -1; + if (DataBuffer.compare(6, 5, "FRONT") == 0) { + channel = atoi(DataBuffer.substr(12).c_str()); + *(m_FrontChannelStatus[detector - 1].begin() + channel) = false; + } + else if (DataBuffer.compare(6, 4, "BACK") == 0) { + channel = atoi(DataBuffer.substr(11).c_str()); + *(m_BackChannelStatus[detector - 1].begin() + channel) = false; + } + } + else { + ReadingStatus = false; + } + } + } + cout << "\t/////////////////////////////////////////////////" << endl; +} + +/////////////////////////////////////////////////////////////////////////// +void TSuperX3Physics::InitSpectra() { m_Spectra = new TSuperX3Spectra(m_NumberOfDetectors); } + +/////////////////////////////////////////////////////////////////////////// +void TSuperX3Physics::FillSpectra() { + m_Spectra->FillRawSpectra(m_EventData); + m_Spectra->FillPreTreatedSpectra(m_PreTreatedData); + m_Spectra->FillPhysicsSpectra(m_EventPhysics); +} + +/////////////////////////////////////////////////////////////////////////// +void TSuperX3Physics::CheckSpectra() { m_Spectra->CheckSpectra(); } + +/////////////////////////////////////////////////////////////////////////// +void TSuperX3Physics::ClearSpectra() { + // To be done +} + +/////////////////////////////////////////////////////////////////////////// +void TSuperX3Physics::WriteSpectra() { m_Spectra->WriteSpectra(); } + +/////////////////////////////////////////////////////////////////////////// +map<string, TH1*> TSuperX3Physics::GetSpectra() { + if (m_Spectra) + return m_Spectra->GetMapHisto(); + else { + map<string, TH1*> empty; + return empty; + } +} + +//////////////////////////////////////////////////////////////////////////////// +// Construct Method to be pass to the DetectorFactory // +//////////////////////////////////////////////////////////////////////////////// +NPL::VDetector* TSuperX3Physics::Construct() { return (NPL::VDetector*)new TSuperX3Physics(); } + +//////////////////////////////////////////////////////////////////////////////// +// Registering the construct method to the factory // +//////////////////////////////////////////////////////////////////////////////// +extern "C" { +class proxy_w1 { + public: + proxy_w1() { + NPL::DetectorFactory::getInstance()->AddToken("SuperX3", "SuperX3"); + NPL::DetectorFactory::getInstance()->AddDetector("SuperX3", TSuperX3Physics::Construct); + } +}; + +proxy_w1 p; +} + diff --git a/NPLib/Detectors/SuperX3/TSuperX3Physics.h b/NPLib/Detectors/SuperX3/TSuperX3Physics.h new file mode 100644 index 000000000..f419d6f0c --- /dev/null +++ b/NPLib/Detectors/SuperX3/TSuperX3Physics.h @@ -0,0 +1,241 @@ +#ifndef __SuperX3Physics__ +#define __SuperX3Physics__ +/***************************************************************************** + * Copyright (C) 2009-2016 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: N. de Sereville contact address: deserevi@ipno.in2p3.fr * + * * + * Creation Date : january 2011 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: * + * This class holds the physics class for the SuperX3 detector from Micron * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ +// C++ headers +#include <map> +#include <vector> +using namespace std; + +// ROOT headers +#include "TCanvas.h" +#include "TH1.h" +#include "TObject.h" +#include "TVector2.h" +#include "TVector3.h" + +// NPTool headers +#include "NPCalibrationManager.h" +#include "NPInputParser.h" +#include "NPVDetector.h" +#include "TSuperX3Data.h" +#include "TSuperX3Spectra.h" + +// forward declaration +class TSuperX3Spectra; + +class TSuperX3Physics : public TObject, public NPL::VDetector { + public: // Constructor and Destructor + TSuperX3Physics(); + ~TSuperX3Physics(){}; + + public: + void Clear(); + void Clear(const Option_t*){}; + + private: // data obtained after BuildPhysicalEvent() and stored in ROOT output file + vector<Int_t> fEventType; + vector<Int_t> fDetectorNumber; + vector<Double_t> fFrontEnergy; + vector<Double_t> fBackEnergy; + vector<Double_t> fHalfEnergy; + vector<Double_t> fFrontTime; + vector<Double_t> fBackTime; + vector<Int_t> fFrontStrip; + vector<Int_t> fBackStrip; + + public: + // setters + void SetEventType(Int_t evtType) { fEventType.push_back(evtType); } + void SetDetectorNumber(Int_t moduleNbr) { fDetectorNumber.push_back(moduleNbr); } + void SetFrontEnergy(Double_t ener) { fFrontEnergy.push_back(ener); } + void SetBackEnergy(Double_t ener) { fBackEnergy.push_back(ener); } + void SethalfEnergy(Double_t ener) { fHalfEnergy.push_back(ener); } + void SetFrontTime(Double_t time) { fFrontTime.push_back(time); } + void SetBackTime(Double_t time) { fBackTime.push_back(time); } + void SetFrontStrip(Int_t x) { fFrontStrip.push_back(x); } + void SetBackStrip(Int_t y) { fBackStrip.push_back(y); } + + // getters + Int_t GetEventMultiplicity() { return fFrontEnergy.size(); } + Int_t GetEventType(Int_t i) { return fEventType[i]; } + Int_t GetDetectorNumber(Int_t i) { return fDetectorNumber[i]; } + Double_t GetFrontEnergy(Int_t i) { return fFrontEnergy[i]; } + Double_t GetBackEnergy(Int_t i) { return fBackEnergy[i]; } + Double_t GetHalfEnergy(Int_t i) { return fHalfEnergy[i]; } + Double_t GetFrontTime(Int_t i) { return fFrontTime[i]; } + Double_t GetBackTime(Int_t i) { return fBackTime[i]; } + Int_t GetFrontStrip(Int_t i) { return fFrontStrip[i]; } + Int_t GetBackStrip(Int_t i) { return fBackStrip[i]; } + + public: + Int_t m_nCounter; //! + Bool_t m_Counter[10]; //! + + public: // inherited from VDetector + // Read stream at ConfigFile to pick-up parameters of detector (Position,...) using Token + void ReadConfiguration(NPL::InputParser); + + // Add parameters to the CalibrationManger + void AddParameterToCalibrationManager(); + + // Activate associated branches and link them to the private member object m_EventData + void InitializeRootInputRaw(); + + // Activate associated branches and link them to the private member m_EventPhysics + void InitializeRootInputPhysics(); + + // Create associated branches and associated private member m_EventPhysics + void InitializeRootOutput(); + + // This method is called at each event read from the Input Tree. Aime is to build treat Raw dat in order to extract + // physical parameter. + void BuildPhysicalEvent(); + + // Same as above, but only the simplest event and/or simple method are used (low multiplicity, faster algorythm but + // less efficient ...). This method aimed to be used for analysis performed during experiment, when speed is + // requiered. NB: This method can eventually be the same as BuildPhysicalEvent. + void BuildSimplePhysicalEvent(); + + // Same as above but for online analysis + void BuildOnlinePhysicalEvent() { BuildPhysicalEvent(); }; + + // Clear raw and physics data + void ClearEventPhysics() { Clear(); } + void ClearEventData() { m_EventData->Clear(); } + + // Methods related to the TSuperX3Spectra classes + // Instantiate the TSuperX3Spectra class and the histograms + void InitSpectra(); + // Fill the spectra defined in TSuperX3Spectra + void FillSpectra(); + // Used for Online mainly, perform check on the histo and for example change their color if issues are found + void CheckSpectra(); + // Used for Online only, clear all the spectra hold by the Spectra class + void ClearSpectra(); + // Write Spectra to file + void WriteSpectra(); + + public: // Specific to SuperX3 + // Remove bad channel, calibrate the data and apply thresholds + void PreTreat(); + + // Clear the pre treated object + void ClearPreTreatedData() { m_PreTreatedData->Clear(); } + + // Return false if the channel is disabled by user + // Frist argument is either "Front" or "Back" + bool IsValidChannel(string Type, int detector, int channel); + + // Initialize the standard parameters for analysis, i.e.: all channel enable, + // maximum multiplicity for strip = number of telescope + void InitializeStandardParameters(); + + // Read the user configuration file; if no file found, load standard one + void ReadAnalysisConfig(); + + // Add detector using cartesian coordinates + void AddDetector(TVector3 C_X1_Y1, TVector3 C_X16_Y1, TVector3 C_X1_Y16, TVector3 C_X16_Y16); + + // Add detector using spherical coordinates + void AddDetector(double theta, double phi, double distance, double beta_u, double beta_v, double beta_w); + + // Give an external TSuperX3Data object to TSuperX3Physics. Needed for online analysis for example. + void SetRawDataPointer(TSuperX3Data* rawDataPointer) { m_EventData = rawDataPointer; } + + // Use for reading Calibration Run, very simple methods; only apply calibration, no condition + void ReadCalibrationRun(){}; + + public: // Methods used for event treatement + Int_t EventType(); + vector<TVector2> Match_Front_Back(); + + private: // Data not written in the tree + TSuperX3Data* m_EventData; //! + TSuperX3Data* m_PreTreatedData; //! + TSuperX3Physics* m_EventPhysics; //! + + public: + TSuperX3Data* GetRawData() const { return m_EventData; } + TSuperX3Data* GetPreTreatedData() const { return m_PreTreatedData; } + + private: // Map of activated Channel + map<int, vector<bool>> m_FrontChannelStatus; //! + map<int, vector<bool>> m_BackChannelStatus; //! + + private: // Parameters used in the analysis + // If multiplicity is greater than m_MaximumStripMultiplicityAllowed + // after PreTreat(), event is not treated + int m_MaximumStripMultiplicityAllowed; //! + + // Tolerance for front / back energy match + double m_StripEnergyMatchingSigma; //! + double m_StripEnergyMatchingNumberOfSigma; //! + + // Energy thresholds + // Raw Threshold + int m_FrontE_Raw_Threshold; //! + int m_BackE_Raw_Threshold; //! + // Calibrated Threshold + double m_FrontE_Calib_Threshold; //! + double m_BackE_Calib_Threshold; //! + + private: // Spatial Position of Strip Calculated on bases of detector position + int m_NumberOfDetectors; //! + vector<vector<vector<double>>> m_StripPositionX; //! + vector<vector<vector<double>>> m_StripPositionY; //! + vector<vector<vector<double>>> m_StripPositionZ; //! + + public: + double GetNumberOfDetectors() { return m_NumberOfDetectors; }; + double GetStripPositionX(int N, int Front, int Back) { return m_StripPositionX[N - 1][Front - 1][Back - 1]; }; + double GetStripPositionY(int N, int Front, int Back) { return m_StripPositionY[N - 1][Front - 1][Back - 1]; }; + double GetStripPositionZ(int N, int Front, int Back) { return m_StripPositionZ[N - 1][Front - 1][Back - 1]; }; + TVector3 GetPositionOfInteraction(int i); + TVector3 GetDetectorNormal(int i); + void DumpStrippingScheme(Int_t detecNumber); + + private: // Geometry and strip number + double m_SiliconFace; //! // mm + int m_NumberOfStrips; //! + double m_StripPitch; //! + + private: // Spectra Class + TSuperX3Spectra* m_Spectra; // ! + + public: // Spectra Getter + map<string, TH1*> GetSpectra(); + + public: // Static constructor to be passed to the Detector Factory + static NPL::VDetector* Construct(); + + ClassDef(TSuperX3Physics, 1) // TSuperX3Physics +}; + +namespace SuperX3_LOCAL { + Double_t fSuperX3_Front_E(TSuperX3Data* EventData, Int_t i); + Double_t fSuperX3_Front_T(TSuperX3Data* EventData, Int_t i); + Double_t fSuperX3_Back_E(TSuperX3Data* EventData, Int_t i); + Double_t fSuperX3_Back_T(TSuperX3Data* EventData, Int_t i); +} // namespace SuperX3_LOCAL + +#endif diff --git a/NPLib/Detectors/SuperX3/TSuperX3Spectra.cxx b/NPLib/Detectors/SuperX3/TSuperX3Spectra.cxx new file mode 100644 index 000000000..fe76f01a2 --- /dev/null +++ b/NPLib/Detectors/SuperX3/TSuperX3Spectra.cxx @@ -0,0 +1,66 @@ +/***************************************************************************** + * Copyright (C) 2009-2016 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: N. de Sereville address: deserevi@ipno.in2p3.fr * + * * + * Creation Date : November 2015 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: * + * This class holds all the online spectra needed for SuperX3 * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + *****************************************************************************/ + +// class header +#include "TSuperX3Spectra.h" + +// C++ headers +#include <iostream> +#include <string> +using namespace std; + +// NPTool headers +#include "NPOptionManager.h" + +//////////////////////////////////////////////////////////////////////////////// +TSuperX3Spectra::TSuperX3Spectra() + : fNumberOfDetectors(0), fNumberOfStripsFront(16), fNumberOfStripsBack(16), fNumberOfCounters(10) { + SetName("SuperX3"); +} + +//////////////////////////////////////////////////////////////////////////////// +TSuperX3Spectra::TSuperX3Spectra(unsigned int NumberOfDetectors) { + InitRawSpectra(); + InitPreTreatedSpectra(); + InitPhysicsSpectra(); +} + +//////////////////////////////////////////////////////////////////////////////// +TSuperX3Spectra::~TSuperX3Spectra() {} + +//////////////////////////////////////////////////////////////////////////////// +void TSuperX3Spectra::InitRawSpectra() {} + +//////////////////////////////////////////////////////////////////////////////// +void TSuperX3Spectra::InitPreTreatedSpectra() {} + +//////////////////////////////////////////////////////////////////////////////// +void TSuperX3Spectra::InitPhysicsSpectra() {} + +//////////////////////////////////////////////////////////////////////////////// +void TSuperX3Spectra::FillRawSpectra(TSuperX3Data* RawData) {} + +//////////////////////////////////////////////////////////////////////////////// +void TSuperX3Spectra::FillPreTreatedSpectra(TSuperX3Data* PreTreatedData) {} + +//////////////////////////////////////////////////////////////////////////////// +void TSuperX3Spectra::FillPhysicsSpectra(TSuperX3Physics* Physics) {} + diff --git a/NPLib/Detectors/SuperX3/TSuperX3Spectra.h b/NPLib/Detectors/SuperX3/TSuperX3Spectra.h new file mode 100644 index 000000000..0e9ce648e --- /dev/null +++ b/NPLib/Detectors/SuperX3/TSuperX3Spectra.h @@ -0,0 +1,58 @@ +#ifndef TSuperX3SPECTRA_H +#define TSuperX3SPECTRA_H +/***************************************************************************** + * Copyright (C) 2009-2016 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: n. de Sereville address: deserevi@ipno.in2p3.fr * + * * + * Creation Date : November 2015 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: * + * This class holds all the online spectra needed for SuperX3 * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + *****************************************************************************/ + +// NPLib headers +#include "NPVSpectra.h" +#include "TSuperX3Data.h" +#include "TSuperX3Physics.h" + +// ForwardDeclaration +class TSuperX3Physics; + +class TSuperX3Spectra : public VSpectra { + public: + // constructor and destructor + TSuperX3Spectra(); + TSuperX3Spectra(unsigned int NumberOfDetectors); + ~TSuperX3Spectra(); + + private: + // Initialization methods + void InitRawSpectra(); + void InitPreTreatedSpectra(); + void InitPhysicsSpectra(); + + public: + // Filling methods + void FillRawSpectra(TSuperX3Data*); + void FillPreTreatedSpectra(TSuperX3Data*); + void FillPhysicsSpectra(TSuperX3Physics*); + + private: // Information on SuperX3 + unsigned int fNumberOfDetectors; + unsigned int fNumberOfStripsFront; + unsigned int fNumberOfStripsBack; + Int_t fNumberOfCounters; +}; + +#endif diff --git a/NPSimulation/Detectors/SuperX3/CMakeLists.txt b/NPSimulation/Detectors/SuperX3/CMakeLists.txt new file mode 100644 index 000000000..d7cfe0cfc --- /dev/null +++ b/NPSimulation/Detectors/SuperX3/CMakeLists.txt @@ -0,0 +1,2 @@ +add_library(NPSSuperX3 SHARED SuperX3.cc ) +target_link_libraries(NPSSuperX3 NPSCore ${ROOT_LIBRARIES} ${Geant4_LIBRARIES} NPSuperX3) diff --git a/NPSimulation/Detectors/SuperX3/SuperX3.cc b/NPSimulation/Detectors/SuperX3/SuperX3.cc new file mode 100644 index 000000000..0307a5bec --- /dev/null +++ b/NPSimulation/Detectors/SuperX3/SuperX3.cc @@ -0,0 +1,366 @@ +/***************************************************************************** + * Copyright (C) 2009-2016 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: N. de Sereville contact address: deserevi@ipno.in2p3.fr * + * * + * Creation Date : 12/01/11 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: Define the SuperX3 detector from Micron * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ + +// C++ headers +#include <cmath> +#include <sstream> +#include <string> + +// G4 Geometry headers +#include "G4Box.hh" +#include "G4Tubs.hh" + +// G4 various headers +#include "G4Colour.hh" +#include "G4Material.hh" +#include "G4PVDivision.hh" +#include "G4PVPlacement.hh" +#include "G4RotationMatrix.hh" +#include "G4Transform3D.hh" +#include "G4VisAttributes.hh" + +// G4 sensitive +#include "G4MultiFunctionalDetector.hh" +#include "G4SDManager.hh" + +// NPTool headers +#include "DSSDScorers.hh" +#include "MaterialManager.hh" +#include "NPOptionManager.h" +#include "NPSDetectorFactory.hh" +#include "ObsoleteGeneralScorers.hh" +#include "RootOutput.h" +#include "SuperX3.hh" +#include "TSuperX3Data.h" + +// CLHEP +#include "CLHEP/Random/RandGauss.h" + +using namespace std; +using namespace CLHEP; +using namespace SuperX3SQUARE; + +SuperX3::SuperX3() : m_Event(new TSuperX3Data) { InitializeMaterials(); } + +SuperX3::~SuperX3() { delete m_Event; } + +void SuperX3::AddDetector(G4ThreeVector X1_Y1, G4ThreeVector X16_Y1, G4ThreeVector X1_Y16, G4ThreeVector X16_Y16) { + m_DefinitionType.push_back(true); + + m_X1_Y1.push_back(X1_Y1); + m_X16_Y1.push_back(X16_Y1); + m_X1_Y16.push_back(X1_Y16); + m_X16_Y16.push_back(X16_Y16); + + m_R.push_back(0); + m_Theta.push_back(0); + m_Phi.push_back(0); + m_beta_u.push_back(0); + m_beta_v.push_back(0); + m_beta_w.push_back(0); +} + +void SuperX3::AddDetector(G4double R, G4double Theta, G4double Phi, G4double beta_u, G4double beta_v, G4double beta_w) { + G4ThreeVector empty = G4ThreeVector(0, 0, 0); + + m_DefinitionType.push_back(false); + + m_R.push_back(R); + m_Theta.push_back(Theta); + m_Phi.push_back(Phi); + m_beta_u.push_back(beta_u); + m_beta_v.push_back(beta_v); + m_beta_w.push_back(beta_w); + + m_X1_Y1.push_back(empty); + m_X16_Y1.push_back(empty); + m_X1_Y16.push_back(empty); + m_X16_Y16.push_back(empty); +} + +void SuperX3::VolumeMaker(G4int DetecNumber, G4ThreeVector position, G4RotationMatrix* rotation, + G4LogicalVolume* world) { + G4double NbrTelescopes = DetecNumber; + G4String DetectorNumber; + ostringstream Number; + Number << NbrTelescopes; + DetectorNumber = Number.str(); + + //////////////////////////////////////////////////////////////// + ////////////// Starting Volume Definition ////////////////////// + //////////////////////////////////////////////////////////////// + G4String Name = "SuperX3Square" + DetectorNumber; + + // Definition of the volume containing the sensitive detector + G4Box* solidSuperX3 = new G4Box(Name, 0.5 * FaceFront, 0.5 * SiliconFaceLength * mm, 0.5 * Length); + G4LogicalVolume* logicSuperX3 = new G4LogicalVolume(solidSuperX3, m_MaterialVacuum, Name, 0, 0, 0); + + new G4PVPlacement(G4Transform3D(*rotation, position), logicSuperX3, Name, world, false, 0); + + logicSuperX3->SetVisAttributes(G4VisAttributes::Invisible); + if (m_non_sensitive_part_visiualisation) + logicSuperX3->SetVisAttributes(G4VisAttributes(G4Colour(0.90, 0.90, 0.90))); + + // Aluminium dead layers + G4ThreeVector positionAluStripFront = G4ThreeVector(0, 0, AluStripBack_PosZ); + G4ThreeVector positionAluStripBack = G4ThreeVector(0, 0, AluStripFront_PosZ); + + G4Box* solidAluStrip = new G4Box("AluBox", 0.5 * SiliconFaceWidth, 0.5 * SiliconFaceLength, 0.5 * AluStripThickness); + // G4LogicalVolume* logicAluStrip = new G4LogicalVolume(solidAluStrip, + // m_MaterialAluminium, "logicAluStrip", 0, 0, 0); + G4LogicalVolume* logicAluStrip = new G4LogicalVolume(solidAluStrip, m_MaterialVacuum, "logicAluStrip", 0, 0, 0); + + new G4PVPlacement(0, positionAluStripFront, logicAluStrip, Name + "_AluStripFront", logicSuperX3, false, 0); + new G4PVPlacement(0, positionAluStripBack, logicAluStrip, Name + "_AluStripBack", logicSuperX3, false, 0); + + logicAluStrip->SetVisAttributes(G4VisAttributes(G4Colour(1.90, 0.0, 0.0))); + + // Silicon detector itself + G4ThreeVector positionSilicon = G4ThreeVector(0, 0, Silicon_PosZ); + + G4Box* solidSilicon = + new G4Box("solidSilicon", 0.5 * SiliconFaceWidth, 0.5 * SiliconFaceLength * mm, 0.5 * SiliconThickness); + G4LogicalVolume* logicSilicon = new G4LogicalVolume(solidSilicon, m_MaterialSilicon, "logicSilicon", 0, 0, 0); + + new G4PVPlacement(0, positionSilicon, logicSilicon, Name + "_Silicon", logicSuperX3, false, 0); + + // Set Silicon strip sensible + logicSilicon->SetSensitiveDetector(m_Scorer); + + /// Visualisation of Silicon Strip + G4VisAttributes* SiliconVisAtt = new G4VisAttributes(G4Colour(0.0, 0.0, 0.9)); + logicSilicon->SetVisAttributes(SiliconVisAtt); +} + +// Read stream at Configfile to pick-up parameters of detector (Position,...) +// Called in DetecorConstruction::ReadDetextorConfiguration Method +void SuperX3::ReadConfiguration(NPL::InputParser parser) { + vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("SuperX3"); + if (NPOptionManager::getInstance()->GetVerboseLevel()) + cout << "//// " << blocks.size() << " detectors found " << endl; + for (unsigned int i = 0; i < blocks.size(); i++) { + // Cartesian Case + vector<string> cart = {"X1_Y1", "X1_Y16", "X16_Y1", "X16_Y16", "VIS"}; + // Spherical Case + vector<string> sphe = {"R", "THETA", "PHI", "BETA", "VIS"}; + + if (blocks[i]->HasTokenList(cart)) { + cout << endl << "//// SuperX3 " << i + 1 << endl; + G4ThreeVector A = NPS::ConvertVector(blocks[i]->GetTVector3("X1_Y1", "mm")); + G4ThreeVector B = NPS::ConvertVector(blocks[i]->GetTVector3("X16_Y1", "mm")); + G4ThreeVector C = NPS::ConvertVector(blocks[i]->GetTVector3("X1_Y16", "mm")); + G4ThreeVector D = NPS::ConvertVector(blocks[i]->GetTVector3("X16_Y16", "mm")); + if (blocks[i]->GetInt("VIS")) + m_non_sensitive_part_visiualisation = true; + AddDetector(A, B, C, D); + } + + else if (blocks[i]->HasTokenList(sphe)) { + cout << endl << "//// SuperX3 " << i + 1 << endl; + double Theta = blocks[i]->GetDouble("THETA", "deg"); + double Phi = blocks[i]->GetDouble("PHI", "deg"); + double R = blocks[i]->GetDouble("R", "mm"); + vector<double> beta = blocks[i]->GetVectorDouble("BETA", "deg"); + if (blocks[i]->GetInt("VIS")) + m_non_sensitive_part_visiualisation = true; + + AddDetector(Theta, Phi, R, beta[0], beta[1], beta[2]); + } + + else { + cout << "ERROR: Missing token for SuperX3 blocks, check your input file" << endl; + exit(1); + } + } +} + +// Construct detector and inialise sensitive part. +// Called After DetecorConstruction::AddDetector Method +void SuperX3::ConstructDetector(G4LogicalVolume* world) { + G4RotationMatrix* SuperX3rot = NULL; + G4ThreeVector SuperX3pos = G4ThreeVector(0, 0, 0); + G4ThreeVector SuperX3u = G4ThreeVector(0, 0, 0); + G4ThreeVector SuperX3v = G4ThreeVector(0, 0, 0); + G4ThreeVector SuperX3w = G4ThreeVector(0, 0, 0); + G4ThreeVector SuperX3Center = G4ThreeVector(0, 0, 0); + + G4int NumberOfDetector = m_DefinitionType.size(); + for (G4int i = 0; i < NumberOfDetector; i++) { + // By Point + if (m_DefinitionType[i]) { + // (u,v,w) unitary vector associated to telescope referencial + // (u,v) // to silicon plan + // w perpendicular to (u,v) plan and pointing ThirdStage + SuperX3u = m_X16_Y1[i] - m_X1_Y1[i]; + SuperX3u = SuperX3u.unit(); + + SuperX3v = m_X1_Y16[i] - m_X1_Y1[i]; + SuperX3v = SuperX3v.unit(); + + SuperX3w = SuperX3u.cross(SuperX3v); + SuperX3w = SuperX3w.unit(); + + SuperX3Center = (m_X1_Y1[i] + m_X1_Y16[i] + m_X16_Y1[i] + m_X16_Y16[i]) / 4; + + // Passage Matrix from Lab Referential to Telescope Referential + SuperX3rot = new G4RotationMatrix(SuperX3u, SuperX3v, SuperX3w); + // translation to place Telescope + SuperX3pos = SuperX3w * Length * 0.5 + SuperX3Center; + } + + // By Angle + else { + G4double Theta = m_Theta[i]; + G4double Phi = m_Phi[i]; + + // (u,v,w) unitary vector associated to telescope referencial + // (u,v) // to silicon plan + // w perpendicular to (u,v) plan and pointing ThirdStage + // Phi is angle between X axis and projection in (X,Y) plan + // Theta is angle between position vector and z axis + G4double wX = m_R[i] * sin(Theta / rad) * cos(Phi / rad); + G4double wY = m_R[i] * sin(Theta / rad) * sin(Phi / rad); + G4double wZ = m_R[i] * cos(Theta / rad); + SuperX3w = G4ThreeVector(wX, wY, wZ); + + // vector corresponding to the center of the module + SuperX3Center = SuperX3w; + + // vector parallel to one axis of silicon plane + G4double ii = cos(Theta / rad) * cos(Phi / rad); + G4double jj = cos(Theta / rad) * sin(Phi / rad); + G4double kk = -sin(Theta / rad); + G4ThreeVector Y = G4ThreeVector(ii, jj, kk); + + SuperX3w = SuperX3w.unit(); + SuperX3u = SuperX3w.cross(Y); + SuperX3v = SuperX3w.cross(SuperX3u); + SuperX3v = SuperX3v.unit(); + SuperX3u = SuperX3u.unit(); + + // Passage Matrix from Lab Referential to Telescope Referential + // MUST2 + SuperX3rot = new G4RotationMatrix(SuperX3u, SuperX3v, SuperX3w); + // Telescope is rotate of Beta angle around SuperX3v axis. + SuperX3rot->rotate(m_beta_u[i], SuperX3u); + SuperX3rot->rotate(m_beta_v[i], SuperX3v); + SuperX3rot->rotate(m_beta_w[i], SuperX3w); + // translation to place Telescope + SuperX3pos = SuperX3w * Length * 0.5 + SuperX3Center; + } + + VolumeMaker(i + 1, SuperX3pos, SuperX3rot, world); + } + + delete SuperX3rot; +} + +// Connect the GaspardTrackingData class to the output TTree +// of the simulation +void SuperX3::InitializeRootOutput() { + RootOutput* pAnalysis = RootOutput::getInstance(); + TTree* pTree = pAnalysis->GetTree(); + if (!pTree->FindBranch("SuperX3")) { + pTree->Branch("SuperX3", "TSuperX3Data", &m_Event); + } + pTree->SetBranchAddress("SuperX3", &m_Event); +} + +// Read sensitive part and fill the Root tree. +// Called at in the EventAction::EndOfEventAvtion +void SuperX3::ReadSensitive(const G4Event*) { + // Clear ROOT objects + m_Event->Clear(); + + auto resistive = (DSSDScorers::PS_Resistive*)m_Scorer->GetPrimitive(0); + auto sizeUp = resistive->GetUpMult(); + for (unsigned int i = 0; i < sizeUp; i++) { + double energy = resistive->GetEnergyUp(i); + double time = resistive->GetTimeUp(i); + int det = resistive->GetDetectorUp(i); + int strip = resistive->GetStripUp(i); + m_Event->SetUpE(det, strip, energy); + m_Event->SetUpT(det, strip, time); + } + auto sizeDown = resistive->GetDownMult(); + for (unsigned int i = 0; i < sizeDown; i++) { + double energy = resistive->GetEnergyDown(i); + double time = resistive->GetTimeDown(i); + int det = resistive->GetDetectorDown(i); + int strip = resistive->GetStripDown(i); + m_Event->SetDownE(det, strip, energy); + m_Event->SetDownT(det, strip, time); + } + auto sizeBack = resistive->GetBackMult(); + for (unsigned int i = 0; i < sizeBack; i++) { + double energy = resistive->GetEnergyBack(i); + double time = resistive->GetTimeBack(i); + int det = resistive->GetDetectorBack(i); + m_Event->SetBackE(det, energy); + m_Event->SetBackT(det, time); + } +} + +void SuperX3::InitializeMaterials() { + m_MaterialSilicon = MaterialManager::getInstance()->GetMaterialFromLibrary("Si"); + m_MaterialAluminium = MaterialManager::getInstance()->GetMaterialFromLibrary("Al"); + m_MaterialIron = MaterialManager::getInstance()->GetMaterialFromLibrary("Fe"); + m_MaterialVacuum = MaterialManager::getInstance()->GetMaterialFromLibrary("Vacuum"); +} + +void SuperX3::InitializeScorers() { + bool already_exist = false; + // Associate Scorer + m_Scorer = CheckScorer("ScorerSuperX3", already_exist); + if (already_exist) + return; + + //..... resistive starts.. + G4VPrimitiveScorer* resistivestrip = + new DSSDScorers::PS_Resistive("resistivestrip", 1, SiliconFaceLength, SiliconFaceWidth, NbStrips); + //... resistive ends...... + // and register it to the multifunctionnal detector + //.... resistive starts... + m_Scorer->RegisterPrimitive(resistivestrip); + //.....resistive ends... + + // Add All Scorer to the Global Scorer Manager + G4SDManager::GetSDMpointer()->AddNewDetector(m_Scorer); +} +//////////////////////////////////////////////////////////////////////////////// +// Construct Method to be pass to the DetectorFactory // +//////////////////////////////////////////////////////////////////////////////// +NPS::VDetector* SuperX3::Construct() { return (NPS::VDetector*)new SuperX3(); } + +//////////////////////////////////////////////////////////////////////////////// +// Registering the construct method to the factory // +//////////////////////////////////////////////////////////////////////////////// +extern "C" { +class proxy_nps_w1 { + public: + proxy_nps_w1() { + NPS::DetectorFactory::getInstance()->AddToken("SuperX3", "SuperX3"); + NPS::DetectorFactory::getInstance()->AddDetector("SuperX3", SuperX3::Construct); + } +}; + +proxy_nps_w1 p_nps_w1; +} diff --git a/NPSimulation/Detectors/SuperX3/SuperX3.hh b/NPSimulation/Detectors/SuperX3/SuperX3.hh new file mode 100644 index 000000000..bff28fdff --- /dev/null +++ b/NPSimulation/Detectors/SuperX3/SuperX3.hh @@ -0,0 +1,174 @@ +#ifndef SuperX3_h +#define SuperX3_h 1 +/***************************************************************************** + * Copyright (C) 2009-2016 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: N. de Sereville contact address: deserevi@ipno.in2p3.fr * + * * + * Creation Date : 12/01/11 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: Define the SuperX3 detector from Micron * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ + +// C++ headers +#include <vector> + +// NPTool header +#include "NPSVDetector.hh" + +// Geant4 headers +#include "G4MultiFunctionalDetector.hh" + +// NPTool - ROOT headers +#include "NPInputParser.h" +#include "TSuperX3Data.h" +using namespace std; +using namespace CLHEP; + +class SuperX3 : public NPS::VDetector { + //////////////////////////////////////////////////// + /////// Default Constructor and Destructor ///////// + //////////////////////////////////////////////////// + public: + SuperX3(); + virtual ~SuperX3(); + + //////////////////////////////////////////////////// + //////// Specific Function of this Class /////////// + //////////////////////////////////////////////////// + public: + // Detector positionning + // By Position Method + void AddDetector(G4ThreeVector TL, G4ThreeVector BL, G4ThreeVector BR, + G4ThreeVector CT); + // By Angle Method + void AddDetector(G4double R, G4double Theta, G4double Phi, G4double beta_u, + G4double beta_v, G4double beta_w); + + // Effectively construct Volume + // Avoid to have two time same code for Angle and Point definition + void VolumeMaker(G4int DetecNumber, G4ThreeVector pos, G4RotationMatrix* rot, + G4LogicalVolume* world); + + ///////////////////////////////////////// + //// Inherite from NPS::VDetector class ///// + ///////////////////////////////////////// + public: + // Read stream at Configfile to pick-up parameters of detector (Position,...) + // Called in DetecorConstruction::ReadDetextorConfiguration Method + void ReadConfiguration(NPL::InputParser); + + // Construct detector and inialise sensitive part. + // Called After DetecorConstruction::AddDetector Method + void ConstructDetector(G4LogicalVolume* world); + + // Add Detector branch to the EventTree. + // Called After DetecorConstruction::AddDetector Method + void InitializeRootOutput(); + + // Initialize all scorers necessary for the detector + void InitializeScorers(); + + // Read sensitive part and fill the Root tree. + // Called at in the EventAction::EndOfEventAvtion + void ReadSensitive(const G4Event* event); + + //////////////////////////////////////////////////// + ///////////Event class to store Data//////////////// + //////////////////////////////////////////////////// + private: + TSuperX3Data* m_Event; + + //////////////////////////////////////////////////// + //////////////////// Scorers /////////////////////// + //////////////////////////////////////////////////// + private: + G4MultiFunctionalDetector* m_Scorer; + + //////////////////////////////////////////////////// + //////////////////// Material ////////////////////// + //////////////////////////////////////////////////// + private: + // Declare all material used by theSuperX3 detector + void InitializeMaterials(); + + // Vacuum + G4Material* m_MaterialVacuum; + // Si + G4Material* m_MaterialSilicon; + // Al + G4Material* m_MaterialAluminium; + // Iron + G4Material* m_MaterialIron; + + //////////////////////////////////////////////////// + ///////////////Private intern Data////////////////// + //////////////////////////////////////////////////// + private: + // True if Define by Position, False is Define by angle + vector<bool> m_DefinitionType; + + // Used for "By Point Definition" + vector<G4ThreeVector> m_X1_Y1; // Top Left Corner Position Vector + vector<G4ThreeVector> m_X1_Y16; // Bottom Left Corner Position Vector + vector<G4ThreeVector> m_X16_Y1; // Bottom Right Corner Position Vector + vector<G4ThreeVector> m_X16_Y16; // Center Corner Position Vector + + // Used for "By Angle Definition" + vector<G4double> m_R; // | + vector<G4double> m_Theta; // > Spherical coordinate of Strips Silicium Plate + vector<G4double> m_Phi; // | + + vector<G4double> m_beta_u; // | + vector<G4double> m_beta_v; // > Tilt angle of the Telescope + vector<G4double> m_beta_w; // | + + // Set to true if you want to see Telescope Frame in your visualisation + bool m_non_sensitive_part_visiualisation; + + public: + static NPS::VDetector* Construct(); +}; + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +namespace SuperX3SQUARE { +// Energy/Time resolutions for the different layers +const G4double EnergyResolution = 12.8e-3; // 30 keV FWHM +const G4double TimeResolution = 0.638; // 1.5 ns (FWHM) + +// Geometry +const G4double FaceFront = 40 * mm; +const G4double Length = 10 * mm; + +// First stage +// const G4double AluStripThickness = 0.00000001*micrometer; +const G4double AluStripThickness = 0.4 * micrometer; +const G4double SiliconThickness = 1000 * micrometer; +const G4double SiliconFace = 40.0 * mm; +const G4double SiliconFaceWidth = 40.0 * mm; +const G4double SiliconFaceLength = 75.0 * mm; +// const G4double SiliconFaceWidtht = 75*mm; +// const G4double SiliconFaceLengtht = 39.6*mm; +// Characteristics +const G4int NbStrips = 4; + +// Starting at the front and going in direction of third stage +const G4double AluStripFront_PosZ = Length * -0.5 + 0.5 * AluStripThickness; +const G4double Silicon_PosZ = + AluStripFront_PosZ + 0.5 * AluStripThickness + 0.5 * SiliconThickness; +const G4double AluStripBack_PosZ = + Silicon_PosZ + 0.5 * SiliconThickness + 0.5 * AluStripThickness; +} // namespace SuperX3SQUARE + +#endif diff --git a/Projects/S034/PhysicsListOption.txt b/Projects/S034/PhysicsListOption.txt index 3b244f709..5a8ee25ad 100644 --- a/Projects/S034/PhysicsListOption.txt +++ b/Projects/S034/PhysicsListOption.txt @@ -1,5 +1,5 @@ EmPhysicsList Option4 -DefaultCutOff 1 +DefaultCutOff 10 IonBinaryCascadePhysics 0 NPIonInelasticPhysics 0 EmExtraPhysics 0 diff --git a/Projects/Strasse/PhysicsListOption.txt b/Projects/Strasse/PhysicsListOption.txt index ff38d11c7..495166f43 100644 --- a/Projects/Strasse/PhysicsListOption.txt +++ b/Projects/Strasse/PhysicsListOption.txt @@ -1,5 +1,5 @@ EmPhysicsList Option4 -DefaultCutOff 1000000 +DefaultCutOff 1000 IonBinaryCascadePhysics 0 NPIonInelasticPhysics 0 EmExtraPhysics 0 diff --git a/Projects/Strasse/geometry/catana_only.detector b/Projects/Strasse/geometry/catana_only.detector index dcbfc9e5d..3445cbaa5 100644 --- a/Projects/Strasse/geometry/catana_only.detector +++ b/Projects/Strasse/geometry/catana_only.detector @@ -11,7 +11,7 @@ Target %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1 Catana CSV - Path= geometry/Catana.csv + Path= geometry/CATANA_nptool_sim.csv Pos= 0 0 100 mm Rshift= 100 micrometer diff --git a/Projects/Strasse/geometry/strasse_CAD.detector b/Projects/Strasse/geometry/strasse_CAD.detector index 28aeff76d..8e917df8d 100644 --- a/Projects/Strasse/geometry/strasse_CAD.detector +++ b/Projects/Strasse/geometry/strasse_CAD.detector @@ -87,15 +87,15 @@ Strasse Outer Ref= 0 0 0 mm %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -%Strasse InactiveMaterial -% Chamber= ./geometry/STRASSE_Chamber.stl -% Stars= ./geometry/STRASSE_StarSupports.stl -% Base= ./geometry/STRASSE_Base.stl -% Blades= ./geometry/STRASSE_Blades.stl +Strasse InactiveMaterial + Chamber= ./geometry/STRASSE_Chamber.stl + Stars= ./geometry/STRASSE_StarSupports.stl + Base= ./geometry/STRASSE_Base.stl + Blades= ./geometry/STRASSE_Blades.stl %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1 -%Catana CSV -% Path= geometry/Catana.csv -% Pos= 0 0 100 mm -% Rshift= 100 micrometer +Catana CSV + Path= geometry/Catana.csv + Pos= 0 0 100 mm + Rshift= 100 micrometer diff --git a/Projects/Strasse/geometry/strasse_July2021.detector b/Projects/Strasse/geometry/strasse_July2021.detector index aa02386b2..fe5d254ca 100644 --- a/Projects/Strasse/geometry/strasse_July2021.detector +++ b/Projects/Strasse/geometry/strasse_July2021.detector @@ -93,9 +93,9 @@ Strasse Outer % Base= ./geometry/STRASSE_Base.stl % Blades= ./geometry/STRASSE_Blades.stl %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1 -%Catana CSV -% Path= geometry/Catana.csv -% Pos= 0 0 100 mm -% Rshift= 100 micrometer +Catana CSV + Path= geometry/Catana.csv + Pos= 0 0 100 mm + Rshift= 100 micrometer diff --git a/Projects/Strasse/reaction/C12_p2p.reaction b/Projects/Strasse/reaction/C12_p2p.reaction index 6e9e34a8c..2da7ec347 100755 --- a/Projects/Strasse/reaction/C12_p2p.reaction +++ b/Projects/Strasse/reaction/C12_p2p.reaction @@ -20,7 +20,7 @@ QFSReaction KnockedOut= 1H Heavy= 11B ExcitationEnergyBeam= 0.0 MeV - ExcitationEnergyHeavy= 0.0 MeV + ExcitationEnergyHeavy= 1.0 MeV MomentumSigma= 50.0 ShootHeavy= 1 ShootLight= 1 diff --git a/Projects/Strasse/reaction/Ca54.source b/Projects/Strasse/reaction/Ca54.source index e75757861..d5bfa4982 100755 --- a/Projects/Strasse/reaction/Ca54.source +++ b/Projects/Strasse/reaction/Ca54.source @@ -12,6 +12,7 @@ Beam MeanX= 0 mm MeanY= 0 mm ExcitationEnergy= 3680 keV + ZEmmission= 0 mm %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% LevelData 54Ca Path= reaction/54Ca.dat diff --git a/Projects/SuperX3/Analysis.cxx b/Projects/SuperX3/Analysis.cxx new file mode 100644 index 000000000..555fb5771 --- /dev/null +++ b/Projects/SuperX3/Analysis.cxx @@ -0,0 +1,83 @@ +/***************************************************************************** + * Copyright (C) 2009-2014 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: Adrien MATTA contact address: a.matta@surrey.ac.uk * + * * + * Creation Date : july 2020 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: * + * Class describing the property of an Analysis object * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ +#include <iostream> +using namespace std; +#include "Analysis.h" +#include "NPAnalysisFactory.h" +#include "NPDetectorManager.h" +#include "NPFunction.h" +#include "NPOptionManager.h" +#include "NPPhysicalConstants.h" +#include "TRandom3.h" + +//////////////////////////////////////////////////////////////////////////////// +Analysis::Analysis() {} +//////////////////////////////////////////////////////////////////////////////// +Analysis::~Analysis() {} + +//////////////////////////////////////////////////////////////////////////////// +void Analysis::Init() { + RC = new TReactionConditions; + + InitOutputBranch(); + InitInputBranch(); + + SuperX3 = (TSuperX3Physics*)m_DetectorManager->GetDetector("SuperX3"); +} + +//////////////////////////////////////////////////////////////////////////////// +void Analysis::TreatEvent() { + // Reinitiate calculated variable + ReInitValue(); +} + +//////////////////////////////////////////////////////////////////////////////// +void Analysis::End() {} +//////////////////////////////////////////////////////////////////////////////// +void Analysis::InitOutputBranch() { + // RootOutput::getInstance()->GetTree()->Branch("Ex",&Ex,"Ex/D"); +} + +//////////////////////////////////////////////////////////////////////////////// +void Analysis::InitInputBranch() { RootInput::getInstance()->GetChain()->SetBranchAddress("ReactionConditions", &RC); } +//////////////////////////////////////////////////////////////////////////////// +void Analysis::ReInitValue() { + // Ex=-1000; +} + +//////////////////////////////////////////////////////////////////////////////// +// Construct Method to be pass to the AnalysisFactory // +//////////////////////////////////////////////////////////////////////////////// +NPL::VAnalysis* Analysis::Construct() { return (NPL::VAnalysis*)new Analysis(); } + +//////////////////////////////////////////////////////////////////////////////// +// Registering the construct method to the factory // +//////////////////////////////////////////////////////////////////////////////// +extern "C" { +class proxy { + public: + proxy() { NPL::AnalysisFactory::getInstance()->SetConstructor(Analysis::Construct); } +}; + +proxy p; +} + diff --git a/Projects/SuperX3/Analysis.h b/Projects/SuperX3/Analysis.h new file mode 100644 index 000000000..3e244d183 --- /dev/null +++ b/Projects/SuperX3/Analysis.h @@ -0,0 +1,52 @@ +#ifndef Analysis_h +#define Analysis_h +/***************************************************************************** + * Copyright (C) 2009-2014 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: Adrien MATTA contact address: a.matta@surrey.ac.uk * + * * + * Creation Date : march 2025 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: * + * Class describing the property of an Analysis object * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ +#include "NPEnergyLoss.h" +#include "NPVAnalysis.h" +#include "RootInput.h" +#include "RootOutput.h" +#include "TReactionConditions.h" +#include "TSuperX3Physics.h" +#include <TLorentzVector.h> +#include <TMath.h> +#include <TRandom3.h> +#include <TVector3.h> + +class Analysis : public NPL::VAnalysis { + public: + Analysis(); + ~Analysis(); + + public: + void Init(); + void TreatEvent(); + void End(); + void InitOutputBranch(); + void InitInputBranch(); + void ReInitValue(); + static NPL::VAnalysis* Construct(); + + private: + TSuperX3Physics* SuperX3; +}; +#endif diff --git a/Projects/SuperX3/CMakeLists.txt b/Projects/SuperX3/CMakeLists.txt new file mode 100644 index 000000000..22c74affd --- /dev/null +++ b/Projects/SuperX3/CMakeLists.txt @@ -0,0 +1,5 @@ +cmake_minimum_required (VERSION 2.8) +# Setting the policy to match Cmake version +cmake_policy(VERSION ${CMAKE_MAJOR_VERSION}.${CMAKE_MINOR_VERSION}) +# include the default NPAnalysis cmake file +include("../../NPLib/ressources/CMake/NPAnalysis.cmake") diff --git a/Projects/SuperX3/SuperX3.detector b/Projects/SuperX3/SuperX3.detector new file mode 100644 index 000000000..3b716d6de --- /dev/null +++ b/Projects/SuperX3/SuperX3.detector @@ -0,0 +1,20 @@ +%%%%%%% +Target + THICKNESS= 1 micrometer + ANGLE= 0 deg + RADIUS= 15 mm + MATERIAL= CD2 + X= 0 mm + Y= 0 mm + Z= 0 mm + NbSlices= 10 + +%%%%%%% Detector 1 %%%%%%% +SuperX3 + THETA= 150 + PHI= 0 + R= 150 + BETA= 0 0 0 + VIS= all + + diff --git a/Projects/SuperX3/energy_loss/.gitignore b/Projects/SuperX3/energy_loss/.gitignore new file mode 100644 index 000000000..e69de29bb diff --git a/Projects/SuperX3/project.config b/Projects/SuperX3/project.config new file mode 100644 index 000000000..756657e98 --- /dev/null +++ b/Projects/SuperX3/project.config @@ -0,0 +1,5 @@ +Project SuperX3 + AnalysisOutput= ./root/analysis + SimulationOutput= ./root/simulation + EnergyLoss= ./energy_loss + diff --git a/Projects/SuperX3/root/analysis/.gitignore b/Projects/SuperX3/root/analysis/.gitignore new file mode 100644 index 000000000..e69de29bb diff --git a/Projects/SuperX3/root/simulation/.gitignore b/Projects/SuperX3/root/simulation/.gitignore new file mode 100644 index 000000000..e69de29bb -- GitLab