diff --git a/Inputs/DetectorConfiguration/Vamos.detector b/Inputs/DetectorConfiguration/Vamos.detector new file mode 100644 index 0000000000000000000000000000000000000000..188e5613668a2c99fb99518b73ec6f4ee70e2f5b --- /dev/null +++ b/Inputs/DetectorConfiguration/Vamos.detector @@ -0,0 +1,98 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +Target + THICKNESS= 5 micrometer + RADIUS= 20 mm + MATERIAL= CD2 + ANGLE= 0 deg + X= 0 mm + Y= 0 mm + Z= 0 mm + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Vamos Configuration +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +Vamos + R= 1 m + Theta= 0 deg + +Vamos DC + Z= 60 mm + Gas= iC4H10 + Pressure= 0.01 bar + Temperature= 295 kelvin + +Vamos DC + Z= 200 mm + Gas= iC4H10 + Pressure= 0.01 bar + Temperature= 295 kelvin + +Vamos BeamCatcher + Material= BC400 + Width= 30 mm + Length= 80 mm + Thickness= 30 mm + Pos= 0 0 600 mm + +Vamos MWPPAC + Z= 750 mm + Gas= iC4H10 + Pressure= 0.01 bar + Temperature= 295 kelvin + +Vamos DC + Z= 940 mm + Gas= iC4H10 + Pressure= 0.01 bar + Temperature= 295 kelvin + +Vamos DC + Z= 1060 mm + Gas= iC4H10 + Pressure= 0.01 bar + Temperature= 295 kelvin + + +Vamos IC + Z= 1200 mm + Thickness= 50 mm + Gas= CF4 + Pressure= 0.2 bar + Temperature= 295 kelvin + +Vamos IC + Z= 1250 mm + Thickness= 50 mm + Gas= CF4 + Pressure= 0.2 bar + Temperature= 295 kelvin + +Vamos IC + Z= 1300 mm + Thickness= 50 mm + Gas= CF4 + Pressure= 0.2 bar + Temperature= 295 kelvin + +Vamos IC + Z= 1375 mm + Thickness= 100 mm + Gas= CF4 + Pressure= 0.2 bar + Temperature= 295 kelvin + +Vamos IC + Z= 1525 mm + Thickness= 200 mm + Gas= CF4 + Pressure= 0.2 bar + Temperature= 295 kelvin + +Vamos IC + Z= 1725 mm + Thickness= 200 mm + Gas= CF4 + Pressure= 0.2 bar + Temperature= 295 kelvin +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% diff --git a/NPLib/Detectors/MUST2/TMust2Physics.cxx b/NPLib/Detectors/MUST2/TMust2Physics.cxx index 40e114fc83407229e6c4b92ca8318320f38f751d..14d18e037da057d2ce4c55d77bed100503d20b65 100644 --- a/NPLib/Detectors/MUST2/TMust2Physics.cxx +++ b/NPLib/Detectors/MUST2/TMust2Physics.cxx @@ -74,108 +74,104 @@ ClassImp(TMust2Physics) m_SiLi_MatchingX.resize(16, 0); m_SiLi_MatchingY.resize(16, 0); - for (int i = 0; i < 16; ++i) { - m_SiLi_MatchingX[0] = 112; - m_SiLi_MatchingY[1] = 112; + m_SiLi_MatchingX[0] = 112; + m_SiLi_MatchingY[1] = 112; - m_SiLi_MatchingX[1] = 112; - m_SiLi_MatchingY[1] = 80; + m_SiLi_MatchingX[1] = 112; + m_SiLi_MatchingY[1] = 80; - m_SiLi_MatchingX[2] = 80; - m_SiLi_MatchingY[2] = 112; + m_SiLi_MatchingX[2] = 80; + m_SiLi_MatchingY[2] = 112; - m_SiLi_MatchingX[3] = 80; - m_SiLi_MatchingY[3] = 80; - // - m_SiLi_MatchingX[4] = 48; - m_SiLi_MatchingY[4] = 80; + m_SiLi_MatchingX[3] = 80; + m_SiLi_MatchingY[3] = 80; + // + m_SiLi_MatchingX[4] = 48; + m_SiLi_MatchingY[4] = 80; - m_SiLi_MatchingX[5] = 48; - m_SiLi_MatchingY[5] = 112; + m_SiLi_MatchingX[5] = 48; + m_SiLi_MatchingY[5] = 112; - m_SiLi_MatchingX[6] = 16; - m_SiLi_MatchingY[6] = 80; + m_SiLi_MatchingX[6] = 16; + m_SiLi_MatchingY[6] = 80; - m_SiLi_MatchingX[7] = 16; - m_SiLi_MatchingY[7] = 112; - // - m_SiLi_MatchingX[8] = 112; - m_SiLi_MatchingY[8] = 16; + m_SiLi_MatchingX[7] = 16; + m_SiLi_MatchingY[7] = 112; + // + m_SiLi_MatchingX[8] = 112; + m_SiLi_MatchingY[8] = 16; - m_SiLi_MatchingX[9] = 112; - m_SiLi_MatchingY[9] = 48; + m_SiLi_MatchingX[9] = 112; + m_SiLi_MatchingY[9] = 48; - m_SiLi_MatchingX[10] = 80; - m_SiLi_MatchingY[10] = 16; + m_SiLi_MatchingX[10] = 80; + m_SiLi_MatchingY[10] = 16; - m_SiLi_MatchingX[11] = 80; - m_SiLi_MatchingY[11] = 48; - // - m_SiLi_MatchingX[12] = 48; - m_SiLi_MatchingY[12] = 48; + m_SiLi_MatchingX[11] = 80; + m_SiLi_MatchingY[11] = 48; + // + m_SiLi_MatchingX[12] = 48; + m_SiLi_MatchingY[12] = 48; - m_SiLi_MatchingX[13] = 48; - m_SiLi_MatchingY[13] = 16; + m_SiLi_MatchingX[13] = 48; + m_SiLi_MatchingY[13] = 16; - m_SiLi_MatchingX[14] = 16; - m_SiLi_MatchingY[14] = 48; + m_SiLi_MatchingX[14] = 16; + m_SiLi_MatchingY[14] = 48; - m_SiLi_MatchingX[15] = 16; - m_SiLi_MatchingY[15] = 16; - } + m_SiLi_MatchingX[15] = 16; + m_SiLi_MatchingY[15] = 16; m_CsI_Size = 32; m_CsI_MatchingX.resize(16, 0); m_CsI_MatchingY.resize(16, 0); - for (int i = 0; i < 16; ++i) { - m_CsI_MatchingX[0] = 112; - m_CsI_MatchingY[0] = 112; + m_CsI_MatchingX[0] = 112; + m_CsI_MatchingY[0] = 112; - m_CsI_MatchingX[1] = 112; - m_CsI_MatchingY[1] = 80; + m_CsI_MatchingX[1] = 112; + m_CsI_MatchingY[1] = 80; - m_CsI_MatchingX[2] = 112; - m_CsI_MatchingY[2] = 48; + m_CsI_MatchingX[2] = 112; + m_CsI_MatchingY[2] = 48; - m_CsI_MatchingX[3] = 112; - m_CsI_MatchingY[3] = 16; - // - m_CsI_MatchingX[4] = 80; - m_CsI_MatchingY[4] = 16; + m_CsI_MatchingX[3] = 112; + m_CsI_MatchingY[3] = 16; + // + m_CsI_MatchingX[4] = 80; + m_CsI_MatchingY[4] = 16; - m_CsI_MatchingX[5] = 80; - m_CsI_MatchingY[5] = 48; + m_CsI_MatchingX[5] = 80; + m_CsI_MatchingY[5] = 48; - m_CsI_MatchingX[6] = 80; - m_CsI_MatchingY[6] = 80; + m_CsI_MatchingX[6] = 80; + m_CsI_MatchingY[6] = 80; - m_CsI_MatchingX[7] = 80; - m_CsI_MatchingY[7] = 112; - // - m_CsI_MatchingX[8] = 48; - m_CsI_MatchingY[8] = 16; + m_CsI_MatchingX[7] = 80; + m_CsI_MatchingY[7] = 112; + // + m_CsI_MatchingX[8] = 48; + m_CsI_MatchingY[8] = 16; - m_CsI_MatchingX[9] = 48; - m_CsI_MatchingY[9] = 48; + m_CsI_MatchingX[9] = 48; + m_CsI_MatchingY[9] = 48; - m_CsI_MatchingX[10] = 48; - m_CsI_MatchingY[10] = 80; + m_CsI_MatchingX[10] = 48; + m_CsI_MatchingY[10] = 80; - m_CsI_MatchingX[11] = 48; - m_CsI_MatchingY[11] = 112; - // - m_CsI_MatchingX[12] = 16; - m_CsI_MatchingY[12] = 16; + m_CsI_MatchingX[11] = 48; + m_CsI_MatchingY[11] = 112; + // + m_CsI_MatchingX[12] = 16; + m_CsI_MatchingY[12] = 16; - m_CsI_MatchingX[13] = 16; - m_CsI_MatchingY[13] = 48; + m_CsI_MatchingX[13] = 16; + m_CsI_MatchingY[13] = 48; - m_CsI_MatchingX[14] = 16; - m_CsI_MatchingY[14] = 80; + m_CsI_MatchingX[14] = 16; + m_CsI_MatchingY[14] = 80; - m_CsI_MatchingX[15] = 16; - m_CsI_MatchingY[15] = 112; - } + m_CsI_MatchingX[15] = 16; + m_CsI_MatchingY[15] = 112; } /////////////////////////////////////////////////////////////////////////// diff --git a/NPLib/Detectors/Vamos/CMakeLists.txt b/NPLib/Detectors/Vamos/CMakeLists.txt index 4b065cf72f16a658df070e6364441b9f18222f73..04ae96ae27bc9b59cd0710bfa2bab2aea72338c6 100644 --- a/NPLib/Detectors/Vamos/CMakeLists.txt +++ b/NPLib/Detectors/Vamos/CMakeLists.txt @@ -1,8 +1,6 @@ -add_custom_command(OUTPUT TVamosPlasticDataDict.cxx COMMAND ../../scripts/build_dict.sh TVamosPlasticData.h TVamosPlasticDataDict.cxx TVamosPlasticData.rootmap libNPVamos.dylib DEPENDS TVamosPlasticData.h) -add_custom_command(OUTPUT TVamosFingerDataDict.cxx COMMAND ../../scripts/build_dict.sh TVamosFingerData.h TVamosFingerDataDict.cxx TVamosFingerData.rootmap libNPVamos.dylib DEPENDS TVamosFingerData.h) -add_custom_command(OUTPUT TVamosDCDataDict.cxx COMMAND ../../scripts/build_dict.sh TVamosDCData.h TVamosDCDataDict.cxx TVamosDCData.rootmap libNPVamos.dylib DEPENDS TVamosDCData.h) -add_custom_command(OUTPUT TVamosCHIODataDict.cxx COMMAND ../../scripts/build_dict.sh TVamosCHIOData.h TVamosCHIODataDict.cxx TVamosCHIOData.rootmap libNPVamos.dylib DEPENDS TVamosCHIOData.h) -add_library(NPVamos SHARED TVamosCHIOData.cxx TVamosDCData.cxx TVamosPlasticData.cxx TVamosFingerData.cxx TVamosCHIODataDict.cxx TVamosDCDataDict.cxx TVamosPlasticDataDict.cxx TVamosFingerDataDict.cxx ) +add_custom_command(OUTPUT TVamosPhysicsDict.cxx COMMAND ../../scripts/build_dict.sh TVamosPhysics.h TVamosPhysicsDict.cxx TVamosPhysics.rootmap libNPVamos.dylib DEPENDS TVamosPhysics.h) +add_custom_command(OUTPUT TVamosDataDict.cxx COMMAND ../../scripts/build_dict.sh TVamosData.h TVamosDataDict.cxx TVamosData.rootmap libNPVamos.dylib DEPENDS TVamosData.h) +add_library(NPVamos SHARED TVamosData.cxx TVamosPhysics.cxx TVamosDataDict.cxx TVamosPhysicsDict.cxx ) target_link_libraries(NPVamos ${ROOT_LIBRARIES} NPCore) -install(FILES TVamosCHIOData.h TVamosDCData.h TVamosPlasticData.h TVamosFingerData.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY}) +install(FILES TVamosData.h TVamosPhysics.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY}) diff --git a/NPLib/Detectors/Vamos/TVamosData.cxx b/NPLib/Detectors/Vamos/TVamosData.cxx new file mode 100644 index 0000000000000000000000000000000000000000..1f4662067d6069dea85b6c188bc72a342da61a0f --- /dev/null +++ b/NPLib/Detectors/Vamos/TVamosData.cxx @@ -0,0 +1,80 @@ +/***************************************************************************** + * Copyright (C) 2009-2018 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: Cyril Lenain contact address: lenain@lpccaen.in2p3.fr * + * * + * Creation Date : octobre 2018 * + * Last update : 09/01/2019 * + *---------------------------------------------------------------------------* + * Decription: * + * This class hold Vamos Raw data * + *---------------------------------------------------------------------------* + * Comment: * + * * + *****************************************************************************/ +#include "TVamosData.h" + +#include <iostream> +#include <fstream> +#include <sstream> +#include <string> +using namespace std; + +ClassImp(TVamosData) + + + ////////////////////////////////////////////////////////////////////// + TVamosData::TVamosData() { + } + + + +////////////////////////////////////////////////////////////////////// +TVamosData::~TVamosData() { +} + + + +////////////////////////////////////////////////////////////////////// +void TVamosData::Clear() { + // Energy + fVamos_E_DetectorNbr.clear(); + fVamos_Energy.clear(); + // Time + fVamos_T_DetectorNbr.clear(); + fVamos_Time.clear(); + // Pos + fVamos_Drift_DetectorNbr.clear(); + fVamos_DriftTime.clear(); + fVamos_X.clear(); +} + + +////////////////////////////////////////////////////////////////////// +void TVamosData::Dump() const { + // This method is very useful for debuging and worth the dev. + cout << "XXXXXXXXXXXXXXXXXXXXXXXX New Event [TVamosData::Dump()] XXXXXXXXXXXXXXXXX" << endl; + + // Energy + size_t mysize = fVamos_E_DetectorNbr.size(); + cout << "Vamos_E_Mult: " << mysize << endl; + + for (size_t i = 0 ; i < mysize ; i++){ + cout << "DetNbr: " << fVamos_E_DetectorNbr[i] + << " Energy: " << fVamos_Energy[i]; + } + + // Time + mysize = fVamos_T_DetectorNbr.size(); + cout << "Vamos_T_Mult: " << mysize << endl; + + for (size_t i = 0 ; i < mysize ; i++){ + cout << "DetNbr: " << fVamos_T_DetectorNbr[i] + << " Time: " << fVamos_Time[i]; + } +} diff --git a/NPLib/Detectors/Vamos/TVamosData.h b/NPLib/Detectors/Vamos/TVamosData.h new file mode 100644 index 0000000000000000000000000000000000000000..b5957639b9de2e0ad1e6ad86c44ccb0b5b5eac6d --- /dev/null +++ b/NPLib/Detectors/Vamos/TVamosData.h @@ -0,0 +1,123 @@ +#ifndef __VamosDATA__ +#define __VamosDATA__ +/***************************************************************************** + * Copyright (C) 2009-2018 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: Cyril Lenain contact address: lenain@lpccaen.in2p3.fr * + * * + * Creation Date : octobre 2018 * + * Last update : 09/01/2019 * + *---------------------------------------------------------------------------* + * Decription: * + * This class hold Vamos Raw data * + *---------------------------------------------------------------------------* + * Comment: * + * * + *****************************************************************************/ + +// STL +#include <vector> +using namespace std; + +// ROOT +#include "TObject.h" + +class TVamosData : public TObject { + ////////////////////////////////////////////////////////////// + // data members are hold into vectors in order + // to allow multiplicity treatment + private: + // Energy + vector<UShort_t> fVamos_E_DetectorNbr; + vector<Double_t> fVamos_Energy; + + // Time + vector<UShort_t> fVamos_T_DetectorNbr; + vector<Double_t> fVamos_Time; + + // DriftTime in DC + vector<UShort_t> fVamos_Drift_DetectorNbr; + vector<Double_t> fVamos_DriftTime; + vector<Double_t> fVamos_X; + + ////////////////////////////////////////////////////////////// + // Constructor and destructor + public: + TVamosData(); + ~TVamosData(); + + ////////////////////////////////////////////////////////////// + // Inherited from TObject and overriden to avoid warnings + public: + void Clear(); + void Clear(const Option_t*){}; + void Dump() const; + + ////////////////////////////////////////////////////////////// + // Getters and Setters + // Prefer inline declaration to avoid unnecessary called of + // frequently used methods + // add //! to avoid ROOT creating dictionnary for the methods + public: + ////////////////////// SETTERS //////////////////////// + // Energy + inline void SetEnergy(const UShort_t& DetNbr, const Double_t& Energy) { + fVamos_E_DetectorNbr.push_back(DetNbr); + fVamos_Energy.push_back(Energy); + }; //! + + // Time + inline void SetTime(const UShort_t& DetNbr, const Double_t& Time) { + fVamos_T_DetectorNbr.push_back(DetNbr); + fVamos_Time.push_back(Time); + }; //! + + // Position DriftTime and X in DC + inline void SetDrift(const UShort_t& DetNbr, const Double_t& DriftTime, + const Double_t& X) { + fVamos_Drift_DetectorNbr.push_back(DetNbr); + fVamos_DriftTime.push_back(DriftTime); + fVamos_X.push_back(X); + }; //! + + ////////////////////// GETTERS //////////////////////// + // Energy + inline UShort_t GetMultEnergy() const { return fVamos_E_DetectorNbr.size(); } + inline UShort_t GetE_DetectorNbr(const unsigned int& i) const { + return fVamos_E_DetectorNbr[i]; + } //! + inline Double_t Get_Energy(const unsigned int& i) const { + return fVamos_Energy[i]; + } //! + + // Time + inline UShort_t GetMultTime() const { return fVamos_T_DetectorNbr.size(); } + inline UShort_t GetT_DetectorNbr(const unsigned int& i) const { + return fVamos_T_DetectorNbr[i]; + } //! + inline Double_t Get_Time(const unsigned int& i) const { + return fVamos_Time[i]; + } //! + + // Position + inline UShort_t GetMultDrift() const { return fVamos_Drift_DetectorNbr.size(); } + inline UShort_t GetDrift_DetectorNbr(const unsigned int& i) const { + return fVamos_Drift_DetectorNbr[i]; + } //! + inline Double_t Get_DriftTime(const unsigned int& i) const { + return fVamos_DriftTime[i]; + } //! + + inline Double_t Get_X(const unsigned int& i) const { return fVamos_X[i]; } //! + + ////////////////////////////////////////////////////////////// + // Required for ROOT dictionnary + ClassDef(TVamosData, 1) // VamosData structure +}; + +#endif diff --git a/NPLib/Detectors/Vamos/TVamosPhysics.cxx b/NPLib/Detectors/Vamos/TVamosPhysics.cxx new file mode 100644 index 0000000000000000000000000000000000000000..0ea9b62bc94a575fbdbe0992cc84d0ef1269054d --- /dev/null +++ b/NPLib/Detectors/Vamos/TVamosPhysics.cxx @@ -0,0 +1,382 @@ +/***************************************************************************** + * Copyright (C) 2009-2018 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: Cyril Lenain contact address: lenain@lpccaen.in2p3.fr * + * * + * Creation Date : octobre 2018 * + * Last update : 09/01/2019 * + *---------------------------------------------------------------------------* + * Decription: * + * This class hold Vamos Treated data * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ + +#include "TVamosPhysics.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 "NPDetectorManager.h" +#include "NPOptionManager.h" +#include "NPSystemOfUnits.h" + +// ROOT +#include "TChain.h" + +ClassImp(TVamosPhysics) + + using namespace NPUNITS; + + /////////////////////////////////////////////////////////////////////////// +TVamosPhysics::TVamosPhysics() + : m_EventData(new TVamosData), + m_EventPhysics(this), + m_NumberOfDetectors(0) { + } + +/////////////////////////////////////////////////////////////////////////// +/// A usefull method to bundle all operation to add a detector +/////////////////////////////////////////////////////////////////////////// +void TVamosPhysics::BuildSimplePhysicalEvent() { + BuildPhysicalEvent(); +} + +/////////////////////////////////////////////////////////////////////////// +void TVamosPhysics::BuildPhysicalEvent() { + + double DriftSpeed = 1*centimeter/microsecond; + double Z_t = m_TargetThickness*0.5 ; + double E = 0 ; + double DriftTime =0; + double Raw_Y = 0; + double t_para =0; + int DetectorNumber = 0; + + // match energy and time together + unsigned int mysizeE = m_EventData->GetMultEnergy(); + unsigned int mysizeT = m_EventData->GetMultTime(); + for (UShort_t e = 0; e < mysizeE ; e++) { + /* for (UShort_t t = 0; t < mysizeT ; t++) { */ + /* if (m_EventData->GetE_DetectorNbr(e) == m_EventData->GetT_DetectorNbr(t)) { */ + if (e!= (mysizeE-1)){ + E+=m_EventData->Get_Energy(e); + + } + else { + /* DetectorNumber.push_back(m_EventData->GetE_DetectorNbr(e)); */ + E+=m_EventData->Get_Energy(e); + Energy.push_back(E); + + /* Time.push_back(m_EventData->Get_Time(t)); */ + + /* } */ + /* } */ + } + } + + // extract X and Y from DriftChamber + unsigned int sizeDrift = m_EventData->GetMultDrift(); + vector<double> X_Event, Y_Event; + + for (UShort_t d=0; d<sizeDrift ; d++){ + DetectorNumber = m_EventData->GetDrift_DetectorNbr(d); + DriftDetectorNumber.push_back(DetectorNumber); + Drift_X.push_back(m_EventData->Get_X(d)); + X_Event.push_back(m_EventData->Get_X(d)); + + DriftTime = m_EventData->Get_DriftTime(d)*microsecond; + + // 7.5 centimeter is half times the height of Chamber + // to separate negative Y from positive Y + if (DriftTime > 7.5*centimeter/DriftSpeed) { + Raw_Y = (DriftTime-7.5*centimeter/DriftSpeed)*DriftSpeed; + Drift_Y.push_back(Raw_Y); + Y_Event.push_back(Raw_Y); + } + + else { + Raw_Y = (7.5*centimeter/DriftSpeed-DriftTime)*DriftSpeed; + Drift_Y.push_back(-Raw_Y); + Y_Event.push_back(-Raw_Y); + } + + } + + // calculate + if (sizeDrift >=2){ + + double X1, X2, Y1, Y2, t, r; + double Y_t_Event, X_t_Event, Phi_t_Event, Theta_t_Event; + + X1=X_Event[0]; + Y1=Y_Event[0]; + X2=X_Event[1]; + Y2=Y_Event[1]; + + // Calculate X,Y on the target + t =(Z_t-Z1)/(Z2-Z1); + X_t_Event = t*(X2-X1)+X1; + Y_t_Event = t*(Y2-Y1)+Y1; + + // Theta and Phi out of the target + r = sqrt(pow(X2-X1,2) + pow(Y2-Y1,2) + pow(Z2-Z1,2)); + Phi_t_Event = acos((Y2-Y1)/r); + Theta_t_Event = acos((X2-X1)/(r*sin(Phi_t_Event))); + + X_target.push_back(X_t_Event); + Y_target.push_back(Y_t_Event); + + Theta_target.push_back((Theta_t_Event*180)/3.14159); + Phi_target.push_back((Phi_t_Event*180)/3.14159); + + // in anticipation to use macro + TargetThickness.push_back(m_TargetThickness); + + } + + Y_Event.clear(); + X_Event.clear(); + +} + +/////////////////////////////////////////////////////////////////////////// +void TVamosPhysics::ReadAnalysisConfig() { + bool ReadingStatus = false; + + // path to file + string FileName = "./configs/ConfigVamos.dat"; + + // open analysis config file + ifstream AnalysisConfigFile; + AnalysisConfigFile.open(FileName.c_str()); + + if (!AnalysisConfigFile.is_open()) { + cout << " No ConfigVamos.dat found: Default parameter loaded for Analayis " << FileName << endl; + return; + } + cout << " Loading user parameter for Analysis from ConfigVamos.dat " << endl; + + // Save it in a TAsciiFile + TAsciiFile* asciiConfig = RootOutput::getInstance()->GetAsciiFileAnalysisConfig(); + asciiConfig->AppendLine("%%% ConfigVamos.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 = "ConfigVamos"; + 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_RAW_THRESHOLD") { + AnalysisConfigFile >> DataBuffer; + m_E_RAW_Threshold = atof(DataBuffer.c_str()); + cout << whatToDo << " " << m_E_RAW_Threshold << endl; + } + + else if (whatToDo=="E_THRESHOLD") { + AnalysisConfigFile >> DataBuffer; + m_E_Threshold = atof(DataBuffer.c_str()); + cout << whatToDo << " " << m_E_Threshold << endl; + } + + else { + ReadingStatus = false; + } + } + } +} + +/////////////////////////////////////////////////////////////////////////// +void TVamosPhysics::Clear() { + DetectorNumber.clear(); + Energy.clear(); + Time.clear(); + DriftDetectorNumber.clear(); + Drift_X.clear(); + Drift_Y.clear(); + X_target.clear(); + Y_target.clear(); + Theta_target.clear(); + Phi_target.clear(); + TargetThickness.clear(); +} + +/////////////////////////////////////////////////////////////////////////// +void TVamosPhysics::ReadConfiguration(NPL::InputParser parser) { + + vector<NPL::InputBlock*> blocks1= parser.GetAllBlocksWithToken("Target"); + if (NPOptionManager::getInstance()->GetVerboseLevel()) + + for (unsigned int i=0; i< blocks1.size(); i++) { + vector<string> TokenTarget = {"THICKNESS", "RADIUS","MATERIAL","ANGLE","X","Y","Z"}; + if (blocks1[i]->HasTokenList(TokenTarget)) { + if (NPOptionManager::getInstance()->GetVerboseLevel()) + cout << endl << "//// Target " << i + 1 << endl; + m_TargetThickness = blocks1[i]->GetDouble("THICKNESS", "micrometer"); + } + } + + vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("Vamos"); + if (NPOptionManager::getInstance()->GetVerboseLevel()) + cout << "//// " << blocks.size() << " detectors found " << endl; + + vector<string> TokenBeamCatcher + = {"Material", "Width", "Length", "Thickness", "Pos"}; + vector<string> sphe = {"R", "Theta"}; + vector<string> TokenMWPPAC = {"Z", "Gas", "Pressure", "Temperature"}; + vector<string> TokenDC = {"Z", "Gas", "Pressure", "Temperature"}; + vector<string> TokenIC = {"Z", "Thickness", "Gas", "Pressure", "Temperature"}; + + for (unsigned int i = 0; i < blocks.size(); i++) { + + if (blocks[i]->HasTokenList(sphe)) { + if (NPOptionManager::getInstance()->GetVerboseLevel()) + cout << endl << "//// Vamos " << i + 1 << endl; + double R = blocks[i]->GetDouble("R", "mm"); + double Theta = blocks[i]->GetDouble("Theta", "deg"); + m_NumberOfDetectors++; + DistanceTargetVamos = R; + } + + else if (blocks[i]->GetMainValue() == "BeamCatcher" + && blocks[i]->HasTokenList(TokenBeamCatcher)) { + if (NPOptionManager::getInstance()->GetVerboseLevel()) + cout << endl << "//// BeamCatcher" << i + 1 << endl; + string Material = blocks[i]->GetString("Material"); + double Width = blocks[i]->GetDouble("Width", "mm"); + double Length = blocks[i]->GetDouble("Length", "mm"); + double Thickness = blocks[i]->GetDouble("Thickness", "mm"); + m_NumberOfDetectors++; + } + + else if (blocks[i]->GetMainValue() == "MWPPAC" + && blocks[i]->HasTokenList(TokenMWPPAC)) { + if (NPOptionManager::getInstance()->GetVerboseLevel()) + cout << endl << "//// MWPPAC" << i + 1 << endl; + double Z = blocks[i]->GetDouble("Z", "mm"); + string Gas = blocks[i]->GetString("Gas"); + double Pressure = blocks[i]->GetDouble("Pressure", "bar"); + double Temperature = blocks[i]->GetDouble("Temperature", "kelvin"); + m_NumberOfDetectors++; + } + + else if (blocks[i]->GetMainValue() == "DC" + && blocks[i]->HasTokenList(TokenDC)) { + if (NPOptionManager::getInstance()->GetVerboseLevel()) + cout << endl << "//// DC" << i + 1 << endl; + double Z = blocks[i]->GetDouble("Z", "mm"); + string Gas = blocks[i]->GetString("Gas"); + double Pressure = blocks[i]->GetDouble("Pressure", "bar"); + double Temperature = blocks[i]->GetDouble("Temperature", "kelvin"); + Z_DriftChambers.push_back(Z); + m_NumberOfDetectors++; + } + + else if (blocks[i]->GetMainValue() == "IC" + && blocks[i]->HasTokenList(TokenIC)) { + if (NPOptionManager::getInstance()->GetVerboseLevel()) + cout << endl << "//// IC" << i + 1 << endl; + double Z = blocks[i]->GetDouble("Z", "mm"); + double Thickness = blocks[i]->GetDouble("Thickness", "mm"); + string Gas = blocks[i]->GetString("Gas"); + double Pressure = blocks[i]->GetDouble("Pressure", "bar"); + double Temperature = blocks[i]->GetDouble("Temperature", "kelvin"); + m_NumberOfDetectors++; + } + + else { + cout << "ERROR: check your input file formatting " << endl; + exit(1); + } + } + + Z1 = Z_DriftChambers[0]+ DistanceTargetVamos; + Z2 = Z_DriftChambers[1]+ DistanceTargetVamos; + /* OriginalBeamEnergy = myReaction.GetBeamEnergy(); */ + + /* cout<<OriginalBeamEnergy <<endl; */ + +} + +/////////////////////////////////////////////////////////////////////////// +void TVamosPhysics::AddParameterToCalibrationManager() { + CalibrationManager* Cal = CalibrationManager::getInstance(); + for (int i = 0; i < m_NumberOfDetectors; ++i) { + Cal->AddParameter("Vamos", "D"+ NPL::itoa(i+1)+"_ENERGY","Vamos_D"+ NPL::itoa(i+1)+"_ENERGY"); + Cal->AddParameter("Vamos", "D"+ NPL::itoa(i+1)+"_TIME","Vamos_D"+ NPL::itoa(i+1)+"_TIME"); + } +} + +/////////////////////////////////////////////////////////////////////////// +void TVamosPhysics::InitializeRootInputRaw() { + TChain* inputChain = RootInput::getInstance()->GetChain(); + inputChain->SetBranchStatus("Vamos", true ); + inputChain->SetBranchAddress("Vamos", &m_EventData ); +} + +/////////////////////////////////////////////////////////////////////////// +void TVamosPhysics::InitializeRootInputPhysics() { + TChain* inputChain = RootInput::getInstance()->GetChain(); + inputChain->SetBranchAddress("Vamos", &m_EventPhysics); +} + +/////////////////////////////////////////////////////////////////////////// +void TVamosPhysics::InitializeRootOutput() { + TTree* outputTree = RootOutput::getInstance()->GetTree(); + outputTree->Branch("Vamos", "TVamosPhysics", &m_EventPhysics); +} + +//////////////////////////////////////////////////////////////////////////////// +// Construct Method to be pass to the DetectorFactory // +//////////////////////////////////////////////////////////////////////////////// +NPL::VDetector* TVamosPhysics::Construct() { + return (NPL::VDetector*) new TVamosPhysics(); +} + +//////////////////////////////////////////////////////////////////////////////// +// Registering the construct method to the factory // +//////////////////////////////////////////////////////////////////////////////// +extern "C"{ + class proxy_Vamos{ + public: + proxy_Vamos(){ + NPL::DetectorFactory::getInstance()->AddToken("Vamos","Vamos"); + NPL::DetectorFactory::getInstance()->AddDetector("Vamos",TVamosPhysics::Construct); + } + }; + + proxy_Vamos p_Vamos; +} + diff --git a/NPLib/Detectors/Vamos/TVamosPhysics.h b/NPLib/Detectors/Vamos/TVamosPhysics.h new file mode 100644 index 0000000000000000000000000000000000000000..ce91ad978e0f933dfc35e4b01f95bdd0bf664436 --- /dev/null +++ b/NPLib/Detectors/Vamos/TVamosPhysics.h @@ -0,0 +1,150 @@ +#ifndef TVamosPHYSICS_H +#define TVamosPHYSICS_H +/***************************************************************************** + * Copyright (C) 2009-2018 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: Cyril Lenain contact address: lenain@lpccaen.in2p3.fr * + * * + * Creation Date : octobre 2018 * + * Last update : 09/01/2019 * + *---------------------------------------------------------------------------* + * Decription: * + * This class hold Vamos Treated data * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ + +// C++ headers +#include <vector> +#include <map> +#include <string> +using namespace std; + +// ROOT headers +#include "TObject.h" +#include "TH1.h" +#include "TVector3.h" +// NPTool headers +#include "TVamosData.h" +#include "NPCalibrationManager.h" +#include "NPVDetector.h" +#include "NPInputParser.h" +#include"NPEnergyLoss.h" +#include"NPReaction.h" + + +class TVamosPhysics : public TObject, public NPL::VDetector { + ////////////////////////////////////////////////////////////// + // constructor and destructor + public: + TVamosPhysics(); + ~TVamosPhysics() {}; + + + ////////////////////////////////////////////////////////////// + // 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> DetectorNumber; + vector<int> DriftDetectorNumber; + vector<double> Energy; + vector<double> Time; + vector<double> Drift_X; + vector<double> Drift_Y; + vector<double> X_target; + vector<double> Y_target; + vector<double> Theta_target; + vector<double> Phi_target; + vector<double> TargetThickness; + + ////////////////////////////////////////////////////////////// + // 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 Vamos array + public: + + // read the user configuration file. If no file is found, load standard one + void ReadAnalysisConfig(); + + // objects are not written in the TTree + private: + TVamosData* m_EventData; //! + /* TVamosData* m_PreTreatedData; //! */ + TVamosPhysics* m_EventPhysics; //! + /* NPL::Reaction myReaction; //! */ + + double m_TargetThickness; + double DistanceTargetVamos; + vector<double> Z_DriftChambers;// + double Z1, Z2 ; + + // getters for raw and pre-treated data object + public: + TVamosData* GetRawData() const {return m_EventData;} + + // parameters used in the analysis + private: + // thresholds + double m_E_RAW_Threshold; //! + double m_E_Threshold; //! + + private: + int m_NumberOfDetectors; //! + + // Static constructor to be passed to the Detector Factory + public: + static NPL::VDetector* Construct(); + + ClassDef(TVamosPhysics,1) // VamosPhysics structure +}; +#endif diff --git a/NPLib/Physics/TInitialConditions.cxx b/NPLib/Physics/TInitialConditions.cxx index 4618cca002f5f61161bb7f8feb730c7960ca0a97..0987d268813fafbf8a612953615ab89735d1e462 100644 --- a/NPLib/Physics/TInitialConditions.cxx +++ b/NPLib/Physics/TInitialConditions.cxx @@ -106,3 +106,12 @@ TVector3 TInitialConditions::GetParticleDirection (const int &i) const { fIC_Momentum_Direction_Z[i]); } +double TInitialConditions::GetParticlePositionX (const int &i, const double &Z) const { + double t = (Z-fIC_Incident_Position_Z) / fIC_Momentum_Direction_Z[i]; + return t*fIC_Momentum_Direction_X[i] + fIC_Incident_Position_X ; +} + +double TInitialConditions::GetParticlePositionY (const int &i, const double &Z) const { + double t = (Z-fIC_Incident_Position_Z) / fIC_Momentum_Direction_Z[i]; + return t*fIC_Momentum_Direction_Y[i] + fIC_Incident_Position_Y; +} diff --git a/NPLib/Physics/TInitialConditions.h b/NPLib/Physics/TInitialConditions.h index c2a235bb4ef2061264178ddf51e19102f9c8f78b..045dd431980f9d3abec0276a442d39c2faa1eeaa 100644 --- a/NPLib/Physics/TInitialConditions.h +++ b/NPLib/Physics/TInitialConditions.h @@ -128,6 +128,10 @@ public: TVector3 GetBeamDirection () const ; TVector3 GetParticleDirection (const int &i) const ; + double GetParticlePositionX (const int &i, const double &Z) const; + double GetParticlePositionY (const int &i, const double &Z) const; + + double GetThetaLab_WorldFrame (const int &i) const { return (GetParticleDirection(i).Angle(TVector3(0,0,1)))/deg; } diff --git a/NPSimulation/Core/MaterialManager.cc b/NPSimulation/Core/MaterialManager.cc index 1ed6000c297625d9c562914b91aef5a5dd384990..914809ff261903b5d9c99a6bcc12df5055fd3785 100644 --- a/NPSimulation/Core/MaterialManager.cc +++ b/NPSimulation/Core/MaterialManager.cc @@ -861,96 +861,90 @@ G4Element* MaterialManager::GetElementFromLibrary(string Name) { } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... -// -G4Material* MaterialManager::GetGasFromLibrary(string Name, double Pressure, - double Temperature) { - ostringstream oss; - oss << Name << "_" << Pressure << "_" << Temperature; - string newName = oss.str(); - map<string, G4Material*>::iterator it; - it = m_Material.find(Name); - double density = 0; - - G4double Vm = 0.08206 * Temperature * atmosphere / (Pressure * kelvin); - - // The element is not found - if (it == m_Material.end()) { - if (Name == "CF4") { // 52 torr - density = 3.72 * kg / m3; - double refTemp = (273.15 + 15) * kelvin; - double refPres = 1.01325 * bar; - density = density * (refTemp / Temperature) / (refPres / Pressure); - G4Material* material = new G4Material("NPS_" + newName, density, 2, - kStateGas, Temperature, Pressure); - material->AddElement(GetElementFromLibrary("C"), 1); - material->AddElement(GetElementFromLibrary("F"), 4); - m_Material[newName] = material; - return material; - } - - if (Name == "He") { - density = (4.0026 / Vm) * mg / cm3; - G4Material* material = new G4Material("NPS_" + newName, density, 1, - kStateGas, Temperature, Pressure); - material->AddElement(GetElementFromLibrary("He"), 1); - m_Material[newName] = material; - return material; - } - - if (Name == "iC4H10" || Name == "Isobutane" || Name == "isobutane") { - density = ((4 * 12.0107 + 10 * 1.00794) / Vm) * mg / cm3; - G4Material* material = new G4Material("NPS_" + newName, density, 2, - kStateGas, Temperature, Pressure); - material->AddElement(GetElementFromLibrary("C"), 4); - material->AddElement(GetElementFromLibrary("H"), 10); - m_Material[newName] = material; - return material; - } - - if (Name == "CH4") { - density = ((12.0107 + 4 * 1.00794) / Vm) * mg / cm3; - G4Material* material = new G4Material("NPS_" + newName, density, 2, - kStateGas, Temperature, Pressure); - material->AddElement(GetElementFromLibrary("C"), 1); - material->AddElement(GetElementFromLibrary("H"), 4); - m_Material[newName] = material; - return material; - } - - if (Name == "CO2") { - density = ((12.0107 + 2 * 16) / Vm) * mg / cm3; - G4Material* material = new G4Material("NPS_" + newName, density, 2, - kStateGas, Temperature, Pressure); - material->AddElement(GetElementFromLibrary("C"), 1); - material->AddElement(GetElementFromLibrary("O"), 2); - m_Material[newName] = material; - return material; - } - - if (Name == "H2") { - density = (2 * 1.00794 / Vm) * mg / cm3; - G4Material* material = new G4Material("NPS_" + newName, density, 1, - kStateGas, Temperature, Pressure); - material->AddElement(GetElementFromLibrary("H"), 2); - // material->AddElement(GetElementFromLibrary("H"), 1); - m_Material[newName] = material; - return material; - } - - if (Name == "D2") { - density = (2 * 2.0140 / Vm) * mg / cm3; - G4Material* material = new G4Material("NPS_" + newName, density, 1, - kStateGas, Temperature, Pressure); - material->AddElement(GetElementFromLibrary("D"), 2); - // material->AddElement(GetElementFromLibrary("D"), 1); - m_Material[newName] = material; - return material; - } - - else { - exit(1); - } - } +// +G4Material* MaterialManager::GetGasFromLibrary(string Name, double Pressure, double Temperature){ + ostringstream oss; + oss << Name<< "_"<<Pressure<<"_"<<Temperature; + string newName= oss.str(); + map<string,G4Material*>::iterator it; + it = m_Material.find(Name); + + double density = 0 ; + + G4double Vm=0.08206*Temperature*atmosphere/(Pressure*kelvin); + + // The element is not found + if(it==m_Material.end()){ + if(Name == "CF4"){ // 52 torr + density = 3.72*kg/m3; + double refTemp= (273.15+15)*kelvin; + double refPres= 1.01325*bar; + density = density*(refTemp/Temperature)/(refPres/Pressure); + G4Material* material = new G4Material("NPS_"+newName,density,2,kStateGas,Temperature,Pressure); + material->AddElement(GetElementFromLibrary("C"), 1); + material->AddElement(GetElementFromLibrary("F"), 4); + m_Material[newName]=material; + return material; + } + + if(Name == "He"){ + density = (4.0026/Vm)*mg/cm3; + G4Material* material = new G4Material("NPS_"+newName,density,1,kStateGas,Temperature,Pressure); + material->AddElement(GetElementFromLibrary("He"), 1); + m_Material[newName]=material; + return material; + } + + if(Name == "iC4H10" || Name == "Isobutane" || Name == "isobutane"){ + density = ((4*12.0107+10*1.00794)/Vm)*mg/cm3; + G4Material* material = new G4Material("NPS_"+newName,density,2,kStateGas,Temperature,Pressure); + material->AddElement(GetElementFromLibrary("C"), 4); + material->AddElement(GetElementFromLibrary("H"), 10); + m_Material[newName]=material; + return material; + } + + if(Name == "CH4"){ + density = ((12.0107+4*1.00794)/Vm)*mg/cm3; + G4Material* material = new G4Material("NPS_"+newName,density,2,kStateGas,Temperature,Pressure); + material->AddElement(GetElementFromLibrary("C"), 1); + material->AddElement(GetElementFromLibrary("H"), 4); + m_Material[newName]=material; + return material; + } + + if(Name == "CO2"){ + density = ((12.0107+2*16)/Vm)*mg/cm3; + G4Material* material = new G4Material("NPS_"+newName,density,2,kStateGas,Temperature,Pressure); + material->AddElement(GetElementFromLibrary("C"), 1); + material->AddElement(GetElementFromLibrary("O"), 2); + m_Material[newName]=material; + return material; + } + + if(Name == "H2"){ + density = (2*1.00794/Vm)*mg/cm3; + G4Material* material = new G4Material("NPS_"+newName,density,1,kStateGas,Temperature,Pressure); + material->AddElement(GetElementFromLibrary("H"), 2); + //material->AddElement(GetElementFromLibrary("H"), 1); + m_Material[newName]=material; + return material; + } + + if(Name == "D2"){ + density = (2*2.0140/Vm)*mg/cm3; + G4Material* material = new G4Material("NPS_"+newName,density,1,kStateGas,Temperature,Pressure); + material->AddElement(GetElementFromLibrary("D"), 2); + //material->AddElement(GetElementFromLibrary("D"), 1); + m_Material[newName]=material; + return material; + } + + + else{ + exit(1); + } + } return NULL; } diff --git a/NPSimulation/Core/PrimaryGeneratorAction.cc b/NPSimulation/Core/PrimaryGeneratorAction.cc index 725041c272860dbe95fa66bcb672c2c382f75a0f..8a5ab2947fd7d720b9629536b6c7d340f3d05508 100644 --- a/NPSimulation/Core/PrimaryGeneratorAction.cc +++ b/NPSimulation/Core/PrimaryGeneratorAction.cc @@ -41,7 +41,7 @@ // Event Generator Class #include "EventGeneratorIsotropic.hh" -#include "EventGeneratorCosmic.hh" +/* #include "EventGeneratorCosmic.hh" */ #include "EventGeneratorMultipleParticle.hh" #include "EventGeneratorBeam.hh" @@ -111,13 +111,13 @@ void PrimaryGeneratorAction::ReadEventGeneratorFile(string Path){ m_EventGenerator.push_back(myEventGenerator); } blocks.clear(); - blocks = parser.GetAllBlocksWithToken("Cosmic"); - if (blocks.size()>0) { - NPS::VEventGenerator* myEventGenerator = new EventGeneratorCosmic(); - myEventGenerator->ReadConfiguration(parser); - myEventGenerator->InitializeRootOutput(); - m_EventGenerator.push_back(myEventGenerator); - } + /* blocks = parser.GetAllBlocksWithToken("Cosmic"); */ + /* if (blocks.size()>0) { */ + /* NPS::VEventGenerator* myEventGenerator = new EventGeneratorCosmic(); */ + /* myEventGenerator->ReadConfiguration(parser); */ + /* myEventGenerator->InitializeRootOutput(); */ + /* m_EventGenerator.push_back(myEventGenerator); */ + /* } */ m_Target=m_detector->GetTarget(); if(m_Target!=NULL){ diff --git a/NPSimulation/Detectors/Vamos/CMakeLists.txt b/NPSimulation/Detectors/Vamos/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..e1307310c822d5aff8449dfcbc51c1d892029fd2 --- /dev/null +++ b/NPSimulation/Detectors/Vamos/CMakeLists.txt @@ -0,0 +1,2 @@ +add_library(NPSVamos SHARED Vamos.cc) +target_link_libraries(NPSVamos NPSCore ${ROOT_LIBRARIES} ${Geant4_LIBRARIES} ${NPLib_LIBRARIES} -lNPVamos) diff --git a/NPSimulation/Detectors/Vamos/Vamos.cc b/NPSimulation/Detectors/Vamos/Vamos.cc new file mode 100644 index 0000000000000000000000000000000000000000..ad808109e2f52f09d70d810bd019869c3b958354 --- /dev/null +++ b/NPSimulation/Detectors/Vamos/Vamos.cc @@ -0,0 +1,640 @@ +/***************************************************************************** + * Copyright (C) 2009-2018 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: Cyril Lenain contact address: lenain@lpccaen.in2p3.fr * + * * + * Creation Date : octobre 2018 * + * Last update : 09/01/2019 * + *---------------------------------------------------------------------------* + * Decription: * + * This class describe the Vamos spectrometer * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + *****************************************************************************/ + +// C++ headers +#include <cmath> +#include <limits> +#include <sstream> +// G4 Geometry object +#include "G4Box.hh" +#include "G4Tubs.hh" + +// G4 sensitive +#include "G4MultiFunctionalDetector.hh" +#include "G4SDManager.hh" + +// G4 various object +#include "G4Colour.hh" +#include "G4Material.hh" +#include "G4PVPlacement.hh" +#include "G4Transform3D.hh" +#include "G4UnionSolid.hh" +#include "G4VisAttributes.hh" + +// NPTool header +#include "CalorimeterScorers.hh" +#include "DriftChamberScorers.hh" +#include "InteractionScorers.hh" +#include "MaterialManager.hh" +#include "NPOptionManager.h" +#include "NPSDetectorFactory.hh" +#include "NPSHitsMap.hh" +#include "RootOutput.h" +#include "Vamos.hh" + +// CLHEP header +#include "CLHEP/Random/RandGauss.h" + +using namespace std; +using namespace CLHEP; + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +namespace Vamos_NS { + + // Energy,time and position Resolution + const G4double EnergyThreshold = -10 * MeV; + const G4double ResoTime = 4.5 * ns; + const G4double ResoEnergy = 0.0001 * MeV; + const G4double ResoDriftTime = 0.0001 * ns; + const G4double ResoPosX = 0.0001 * mm; + + // Drift features + const G4double DriftSpeed = 1 * cm / microsecond; + const G4ThreeVector DriftDir = G4ThreeVector(0, 1, 0); + + // Geometry + const G4double E_DCWidth = 600 * mm; // Entrance_DriftChamber + const G4double E_DCLength = 150 * mm; + const G4double E_DCThickness = 60 * mm; + + const G4double MagnetWidth = 1000 * mm; + const G4double MagnetLenght = 150 * mm; + const G4double MagnetThickness = 200 * mm; + const G4double DipolThickness = 600 * mm; + + const G4double ChamberWidth = 1000 * mm; + const G4double ChamberLength = 150 * mm; + const G4double ChamberThickness = 120 * mm; + + // Mother Volume of Vamos + const G4double Phi = 0 * deg; + const G4double VamosVolumeWidth = 4500 * mm; + const G4double VamosVolumeLength = 1000 * mm; + const G4double VamosVolumeThickness = 5000 * mm; + + // SubVolume of detection + const G4double DetectionVolumeThickness = 1300 * mm; + const G4double DetectionVolumeWidth = 1000 * mm; + const G4double DetectionVolumeLength = 1000 * mm; + const G4double Det_Theta = 0 * deg; + +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +// Vamos Specific Method +Vamos::Vamos() { + + m_Event = new TVamosData(); + m_CalorimeterScorer = 0; + m_DCScorer = 0; + m_Quad1 = 0; + m_Quad2 = 0; + m_Dipol = 0; + m_BeamCatcher = 0; + m_MWPPAC = 0; + m_DC3 = 0; + m_DC4 = 0; + m_DC1 = 0; + m_DC2 = 0; + + ICcounter = 0; + + // RGB Color + Transparency + m_VisQuad = new G4VisAttributes(G4Colour(1.0, 0.0, 0.0)); + m_VisVolumeVamos = new G4VisAttributes(G4Colour(1.0, 0.0, 0.0, 0.1)); + m_VisDC = new G4VisAttributes(G4Colour(0.0, 1.0, 0.0, 0.2)); + m_VisCatcher = new G4VisAttributes(G4Colour(0.0, 0.0, 1.0)); + m_VisGasC4H10 = new G4VisAttributes(G4Colour(0.0, 1.0, 0.0, 0.2)); + m_VisGasCF4 = new G4VisAttributes(G4Colour(0.0, 0.0, 1.0, 0.2)); +} + +Vamos::~Vamos() {} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +using namespace Vamos_NS; + +void Vamos::AddVamos(G4double R, double Theta) { + m_R = R; + m_Theta = Theta; +} + +void Vamos::AddBeamCatcher(string Material, G4double Width, double Length, + double Thickness, G4ThreeVector Pos) { + CatcherMaterial = Material; + CatcherWidth = Width; + CatcherLength = Length; + CatcherThickness = Thickness; + R_Catcher = Pos[2]; + Pos[2] = -DetectionVolumeThickness * 0.5 + CatcherThickness * 0.5; + m_PosCatcher = Pos; +} + +// To add DriftChambers and the MWPPAC +void Vamos::AddDetector(G4double Z, string Gas, double Pressure, + double Temperature) { + m_Z.push_back(Z); + m_Gas.push_back(Gas); + m_Pressure.push_back(Pressure); + m_Temperature.push_back(Temperature); +} + +void Vamos::AddIC(G4double Z, double Thickness, string Gas, double Pressure, + double Temperature) { + m_ZIC.push_back(Z); + m_ThicknessIC.push_back(Thickness); + m_GasIC.push_back(Gas); + m_PressureIC.push_back(Pressure); + m_TemperatureIC.push_back(Temperature); +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +// The two entry DriftChambers +G4LogicalVolume* Vamos::BuildDC1() { + if (!m_DC1) { + G4Box* box = new G4Box("Vamos_DC1", E_DCWidth * 0.5, E_DCLength * 0.5, + E_DCThickness * 0.5); + + G4Material* DetectorMaterial + = MaterialManager::getInstance()->GetGasFromLibrary( + m_Gas[0], m_Pressure[0], m_Temperature[0]); + + m_DC1 = new G4LogicalVolume(box, DetectorMaterial, "logic_Vamos_DC1", 0, 0, + 0); + m_DC1->SetVisAttributes(m_VisGasC4H10); + m_DC1->SetSensitiveDetector(m_DCScorer); + } + return m_DC1; +} + +G4LogicalVolume* Vamos::BuildDC2() { + if (!m_DC2) { + G4Box* box = new G4Box("Vamos_DC2", E_DCWidth * 0.5, E_DCLength * 0.5, + E_DCThickness * 0.5); + + G4Material* DetectorMaterial + = MaterialManager::getInstance()->GetGasFromLibrary( + m_Gas[1], m_Pressure[1], m_Temperature[1]); + + m_DC2 = new G4LogicalVolume(box, DetectorMaterial, "logic_Vamos_DC2", 0, 0, + 0); + m_DC2->SetVisAttributes(m_VisGasC4H10); + m_DC2->SetSensitiveDetector(m_DCScorer); + } + return m_DC2; +} + +// Quadruples and Dipole just to make the visualisation nice +G4LogicalVolume* Vamos::BuildQuad1() { + if (!m_Quad1) { + G4Box* box = new G4Box("Vamos_Box", Vamos_NS::MagnetWidth * 0.5, + Vamos_NS::MagnetLenght * 0.5, + Vamos_NS::MagnetThickness * 0.5); + + G4Material* VamosMaterial + = MaterialManager::getInstance()->GetMaterialFromLibrary("Vacuum"); + m_Quad1 = new G4LogicalVolume(box, VamosMaterial, "logic_Quad1", 0, 0, 0); + m_Quad1->SetVisAttributes(m_VisQuad); + } + return m_Quad1; +} + +G4LogicalVolume* Vamos::BuildQuad2() { + if (!m_Quad2) { + G4Box* box = new G4Box("Vamos_Quad2", Vamos_NS::MagnetWidth * 0.5, + Vamos_NS::MagnetLenght * 0.5, + Vamos_NS::MagnetThickness * 0.5); + + G4Material* VamosMaterial + = MaterialManager::getInstance()->GetMaterialFromLibrary("Vacuum"); + m_Quad1 = new G4LogicalVolume(box, VamosMaterial, "logic_Quad1", 0, 0, 0); + m_Quad1->SetVisAttributes(m_VisQuad); + } + return m_Quad1; +} + +G4LogicalVolume* Vamos::BuildDipol() { + if (!m_Dipol) { + G4Box* box = new G4Box("Vamos_Dipol", Vamos_NS::MagnetWidth * 0.5, + Vamos_NS::MagnetLenght * 0.5, + Vamos_NS::DipolThickness * 0.5); + + G4Material* VamosMaterial + = MaterialManager::getInstance()->GetMaterialFromLibrary("Vacuum"); + m_Dipol = new G4LogicalVolume(box, VamosMaterial, "logic_Dipol", 0, 0, 0); + m_Dipol->SetVisAttributes(m_VisQuad); + } + return m_Dipol; +} + +// Detection at the end of Vamos +G4LogicalVolume* Vamos::BuildBeamCatcher() { + if (!m_BeamCatcher) { + G4Box* box = new G4Box("Vamos_Catcher", CatcherWidth * 0.5, + CatcherLength * 0.5, CatcherThickness * 0.5); + + G4Material* Material + = MaterialManager::getInstance()->GetMaterialFromLibrary( + CatcherMaterial); + m_BeamCatcher + = new G4LogicalVolume(box, Material, "logic_Vamos_Catcher", 0, 0, 0); + m_BeamCatcher->SetVisAttributes(m_VisCatcher); + m_BeamCatcher->SetSensitiveDetector(m_CalorimeterScorer); + } + return m_BeamCatcher; +} + +G4LogicalVolume* Vamos::BuildMWPPAC() { + if (!m_MWPPAC) { + G4Box* box = new G4Box("Vamos_MWPPAC", ChamberWidth * 0.5, + ChamberLength * 0.5, ChamberThickness * 0.5); + + G4Material* DetectorMaterial + = MaterialManager::getInstance()->GetGasFromLibrary( + m_Gas[2], m_Pressure[2], m_Temperature[2]); + m_MWPPAC = new G4LogicalVolume(box, DetectorMaterial, "logic_Vamos_MWPPAC", + 0, 0, 0); + m_MWPPAC->SetVisAttributes(m_VisGasC4H10); + } + return m_MWPPAC; +} + +G4LogicalVolume* Vamos::BuildDC3() { + if (!m_DC3) { + G4Box* box = new G4Box("Vamos_DC3", ChamberWidth * 0.5, ChamberLength * 0.5, + ChamberThickness * 0.5); + + G4Material* DetectorMaterial + = MaterialManager::getInstance()->GetGasFromLibrary( + m_Gas[3], m_Pressure[3], m_Temperature[3]); + + m_DC3 = new G4LogicalVolume(box, DetectorMaterial, "logic_Vamos_DC3", 0, 0, + 0); + m_DC3->SetVisAttributes(m_VisGasC4H10); + m_DC3->SetSensitiveDetector(m_DCScorer); + } + return m_DC3; +} + +G4LogicalVolume* Vamos::BuildDC4() { + if (!m_DC4) { + G4Box* box = new G4Box("Vamos_DC4", ChamberWidth * 0.5, ChamberLength * 0.5, + ChamberThickness * 0.5); + + G4Material* DetectorMaterial + = MaterialManager::getInstance()->GetGasFromLibrary( + m_Gas[4], m_Pressure[4], m_Temperature[4]); + + m_DC4 = new G4LogicalVolume(box, DetectorMaterial, "logic_Vamos_DC4", 0, 0, + 0); + m_DC4->SetVisAttributes(m_VisGasC4H10); + m_DC4->SetSensitiveDetector(m_DCScorer); + } + return m_DC4; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +// In anticipation of use a macro +void Vamos::ClearGeometry(){ + + m_Z.clear(); + m_Gas.clear(); + m_Pressure.clear(); + m_Temperature.clear(); + + m_ZIC.clear(); + m_ThicknessIC.clear(); + m_PressureIC.clear(); + m_TemperatureIC.clear(); + m_GasIC.clear(); + +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +// Virtual Method of NPS:VDetector class +// Read stream at Configfile to pick-up parameters of detector (Position,...) +// Called in DetecorConstruction::ReadDetectorConfiguration Method + +void Vamos::ReadConfiguration(NPL::InputParser parser) { + + vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("Vamos"); + if (NPOptionManager::getInstance()->GetVerboseLevel()) + cout << "//// " << blocks.size() << " detectors found " << endl; + + vector<string> TokenBeamCatcher + = {"Material", "Width", "Length", "Thickness", "Pos"}; + vector<string> sphe = {"R", "Theta"}; + vector<string> TokenMWPPAC = {"Z", "Gas", "Pressure", "Temperature"}; + vector<string> TokenDC = {"Z", "Gas", "Pressure", "Temperature"}; + vector<string> TokenIC = {"Z", "Thickness", "Gas", "Pressure", "Temperature"}; + + for (unsigned int i = 0; i < blocks.size(); i++) { + if (blocks[i]->HasTokenList(sphe)) { + if (NPOptionManager::getInstance()->GetVerboseLevel()) + cout << endl << "//// Vamos " << i + 1 << endl; + G4double R = blocks[i]->GetDouble("R", "mm"); + G4double Theta = blocks[i]->GetDouble("Theta", "deg"); + AddVamos(R, Theta); + } + + else if (blocks[i]->GetMainValue() == "BeamCatcher" + && blocks[i]->HasTokenList(TokenBeamCatcher)) { + if (NPOptionManager::getInstance()->GetVerboseLevel()) + cout << endl << "//// BeamCatcher" << i + 1 << endl; + string Material = blocks[i]->GetString("Material"); + G4double Width = blocks[i]->GetDouble("Width", "mm"); + G4double Length = blocks[i]->GetDouble("Length", "mm"); + G4double Thickness = blocks[i]->GetDouble("Thickness", "mm"); + G4ThreeVector Pos + = NPS::ConvertVector(blocks[i]->GetTVector3("Pos", "mm")); + AddBeamCatcher(Material, Width, Length, Thickness, Pos); + } + + else if (blocks[i]->GetMainValue() == "MWPPAC" + && blocks[i]->HasTokenList(TokenMWPPAC)) { + if (NPOptionManager::getInstance()->GetVerboseLevel()) + cout << endl << "//// MWPPAC" << i + 1 << endl; + G4double Z = blocks[i]->GetDouble("Z", "mm"); + string Gas = blocks[i]->GetString("Gas"); + G4double Pressure = blocks[i]->GetDouble("Pressure", "bar"); + G4double Temperature = blocks[i]->GetDouble("Temperature", "kelvin"); + AddDetector(Z, Gas, Pressure, Temperature); + } + + else if (blocks[i]->GetMainValue() == "DC" + && blocks[i]->HasTokenList(TokenDC)) { + if (NPOptionManager::getInstance()->GetVerboseLevel()) + cout << endl << "//// DC" << i + 1 << endl; + G4double Z = blocks[i]->GetDouble("Z", "mm"); + string Gas = blocks[i]->GetString("Gas"); + G4double Pressure = blocks[i]->GetDouble("Pressure", "bar"); + G4double Temperature = blocks[i]->GetDouble("Temperature", "kelvin"); + AddDetector(Z, Gas, Pressure, Temperature); + } + + else if (blocks[i]->GetMainValue() == "IC" + && blocks[i]->HasTokenList(TokenIC)) { + if (NPOptionManager::getInstance()->GetVerboseLevel()) + cout << endl << "//// IC" << ICcounter+1 << endl; + G4double Z = blocks[i]->GetDouble("Z", "mm"); + G4double Thickness = blocks[i]->GetDouble("Thickness", "mm"); + string Gas = blocks[i]->GetString("Gas"); + G4double Pressure = blocks[i]->GetDouble("Pressure", "bar"); + G4double Temperature = blocks[i]->GetDouble("Temperature", "kelvin"); + AddIC(Z, Thickness, Gas, Pressure, Temperature); + ICcounter++; + } + + else { + cout << "ERROR: check your input file formatting " << endl; + exit(1); + } + + } +} + +// Construct detector and inialise sensitive part. +// Called After DetecorConstruction::AddDetector Method +void Vamos::ConstructDetector(G4LogicalVolume* world) { + + // Mother Volume of Vamos + G4double R = m_R + VamosVolumeThickness * 0.5; + + G4double X = R * sin(m_Theta) * cos(Phi); + G4double Y = R * sin(m_Theta) * sin(Phi); + G4double Z = R * cos(m_Theta); + G4ThreeVector Det_pos = G4ThreeVector(X, Y, Z); + + G4RotationMatrix* Rot1 = new G4RotationMatrix(); + Rot1->rotateY(m_Theta); + + G4Box* MotherSolid + = new G4Box("MotherVolume", VamosVolumeWidth * 0.5, + VamosVolumeLength * 0.5, VamosVolumeThickness * 0.5); + + G4Material* VolumeMaterial + = MaterialManager::getInstance()->GetMaterialFromLibrary("Vacuum"); + G4LogicalVolume* MotherVolume = new G4LogicalVolume( + MotherSolid, VolumeMaterial, "MotherVolume", 0, 0, 0); + + new G4PVPlacement(G4Transform3D(*Rot1, Det_pos), MotherVolume, "MotherVolume", + world, false, 0); + MotherVolume->SetVisAttributes(m_VisVolumeVamos); + + // SubVolume of Detection at the end of Vamos + // The position of the subvolume is defined by the position of the BeamCatcher + G4double R2 + = R_Catcher - CatcherThickness * 0.5 + DetectionVolumeThickness * 0.5; + + G4double X2 = R2 * sin(Det_Theta) * cos(0); + G4double Z2 = R2 * cos(Det_Theta); + G4ThreeVector Det_pos2 = G4ThreeVector(X2, 0, Z2); + + G4RotationMatrix* Rot2 = new G4RotationMatrix(); + Rot2->rotateY(Det_Theta); + + G4Box* MotherDetectorSolid + = new G4Box("MotherDetector", DetectionVolumeWidth * 0.5, + DetectionVolumeLength * 0.5, DetectionVolumeThickness * 0.5); + + G4LogicalVolume* MotherDetector = new G4LogicalVolume( + MotherDetectorSolid, VolumeMaterial, "MotherDetector", 0, 0, 0); + + new G4PVPlacement(G4Transform3D(*Rot2, Det_pos2), MotherDetector, + "MotherDetector", MotherVolume, false, 0); + MotherDetector->SetVisAttributes(m_VisVolumeVamos); + + + // Position the entry DCs and the magnets in the MotherVolume + new G4PVPlacement(0, G4ThreeVector(0, 0, -VamosVolumeThickness * 0.5 + m_Z[0]), + BuildDC1(), "Entrance_DC1", MotherVolume, false, 1); + + new G4PVPlacement(0, G4ThreeVector(0, 0, -VamosVolumeThickness * 0.5 + m_Z[1]), + BuildDC2(), "Entrance_DC2", MotherVolume, false, 2); + + new G4PVPlacement( + 0, + G4ThreeVector(0, 0, (-VamosVolumeThickness + MagnetThickness) * 0.5 + 400), + BuildQuad1(), "Vamos", MotherVolume, false, 0); + + new G4PVPlacement( + 0, + G4ThreeVector(0, 0, + (-VamosVolumeThickness + MagnetThickness) * 0.5 + 700 * mm), + BuildQuad2(), "Vamos", MotherVolume, false, 0); + + new G4PVPlacement( + 0, + G4ThreeVector(0, 0, + (-VamosVolumeThickness + MagnetThickness) * 0.5 + 1500 * mm), + BuildDipol(), "Vamos", MotherVolume, false, 0); + + // Position the system of detection at the end of Vamos in the sub Volume + new G4PVPlacement(0, m_PosCatcher, BuildBeamCatcher(), "BeamCatcher", + MotherDetector, false, 3); + + new G4PVPlacement( + 0, + G4ThreeVector(0, 0, + -DetectionVolumeThickness * 0.5 + + (m_Z[2] - R_Catcher + CatcherThickness * 0.5)), + BuildMWPPAC(), "MWPPAC", MotherDetector, false, 4); + + new G4PVPlacement( + 0, + G4ThreeVector(0, 0, + -DetectionVolumeThickness * 0.5 + + (m_Z[3] - R_Catcher + CatcherThickness * 0.5)), + BuildDC3(), "DC", MotherDetector, false, 5); + + new G4PVPlacement( + 0, + G4ThreeVector(0, 0, + -DetectionVolumeThickness * 0.5 + + (m_Z[4] - R_Catcher + CatcherThickness * 0.5)), + BuildDC4(), "DC", MotherDetector, false, 6); + + // Construct and position the Ionisations Chambers + for (int i = 0; i < ICcounter; i++) { + G4Box* box = new G4Box("Vamos_IC", ChamberWidth * 0.5, ChamberLength * 0.5, + m_ThicknessIC[i] * 0.5); + G4Material* GasIC = MaterialManager::getInstance()->GetGasFromLibrary( + m_GasIC[i], m_PressureIC[i], m_TemperatureIC[i]); + G4LogicalVolume* IC + = new G4LogicalVolume(box, GasIC, "logic_Vamos_IC", 0, 0, 0); + IC->SetVisAttributes(m_VisGasCF4); + IC->SetSensitiveDetector(m_CalorimeterScorer); + new G4PVPlacement( + 0, + G4ThreeVector(0, 0, + -DetectionVolumeThickness * 0.5 + + (m_ZIC[i] - R_Catcher + CatcherThickness * 0.5)), + IC, "IC", MotherDetector, false, i + 7); + } +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +// Add Detector branch to the EventTree. +// Called After DetecorConstruction::AddDetector Method + +void Vamos::InitializeRootOutput() { + RootOutput* pAnalysis = RootOutput::getInstance(); + TTree* pTree = pAnalysis->GetTree(); + if (!pTree->FindBranch("Vamos")) { + pTree->Branch("Vamos", "TVamosData", &m_Event); + } + pTree->SetBranchAddress("Vamos", &m_Event); +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +// Read sensitive part and fill the Root Three +// Called at in the EventAction::EndOfEventAvtion +void Vamos::ReadSensitive(const G4Event*) { + + m_Event->Clear(); + + /////////// + // Calorimeter scorer + CalorimeterScorers::PS_Calorimeter* Scorer + = (CalorimeterScorers::PS_Calorimeter*)m_CalorimeterScorer->GetPrimitive( + 0); + unsigned int size = Scorer->GetMult(); + for (unsigned int i = 0; i < size; i++) { + vector<unsigned int> level = Scorer->GetLevel(i); + G4double Energy + = RandGauss::shoot(Scorer->GetEnergy(i), Vamos_NS::ResoEnergy); + if (Energy > Vamos_NS::EnergyThreshold) { + G4double Time = RandGauss::shoot(Scorer->GetTime(i), Vamos_NS::ResoTime); + int DetectorNbr = level[0]; + m_Event->SetEnergy(DetectorNbr, Energy); + m_Event->SetTime(DetectorNbr, Time); + } + } + /////////// + // DriftChamber scorer + DriftChamberScorers::PS_DriftChamber* Scorer2 + = (DriftChamberScorers::PS_DriftChamber*)m_DCScorer->GetPrimitive(0); + unsigned int size2 = Scorer2->GetMult(); + for (unsigned int i = 0; i < size2; i++) { + vector<unsigned int> level = Scorer2->GetLevel(i); + G4double DriftTime + = RandGauss::shoot(Scorer2->GetDriftTime(i)/Scorer2->GetCounter(i), Vamos_NS::ResoDriftTime); + G4double X + = RandGauss::shoot(Scorer2->GetX(i)/Scorer2->GetCounter(i), ResoPosX); + int DetectorNbr = level[0]; + m_Event->SetDrift(DetectorNbr, DriftTime, X); + } + +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//////////////////////////////////////////////////////////////// +void Vamos::InitializeScorers() { + // This check is necessary in case the geometry is reloaded + bool already_exist = false; + m_DCScorer = CheckScorer("DCScorer", already_exist); + m_CalorimeterScorer = CheckScorer("VamosScorer", already_exist); + if (already_exist) + return; + + // Otherwise the scorer is initialised + vector<int> level; + level.push_back(0); + + G4VPrimitiveScorer* Calorimeter + = new CalorimeterScorers::PS_Calorimeter("Calorimeter", level, 0); + m_CalorimeterScorer->RegisterPrimitive(Calorimeter); + + G4VPrimitiveScorer* Drift = new DriftChamberScorers::PS_DriftChamber( + "Drift", level, DriftDir, DriftSpeed, 0); + m_DCScorer->RegisterPrimitive(Drift); + + // and register it to the multifunctionnal detector + G4SDManager::GetSDMpointer()->AddNewDetector(m_DCScorer); + G4SDManager::GetSDMpointer()->AddNewDetector(m_CalorimeterScorer); +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//////////////////////////////////////////////////////////////////////////////// +// Construct Method to be pass to the DetectorFactory // +//////////////////////////////////////////////////////////////////////////////// +NPS::VDetector* Vamos::Construct() { return (NPS::VDetector*)new Vamos(); } + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//////////////////////////////////////////////////////////////////////////////// +// Registering the construct method to the factory // +//////////////////////////////////////////////////////////////////////////////// +extern "C" { + class proxy_nps_Vamos { + public: + proxy_nps_Vamos() { + NPS::DetectorFactory::getInstance()->AddToken("Vamos", "Vamos"); + NPS::DetectorFactory::getInstance()->AddDetector("Vamos", Vamos::Construct); + } + }; + + proxy_nps_Vamos p_nps_Vamos; +} diff --git a/NPSimulation/Detectors/Vamos/Vamos.hh b/NPSimulation/Detectors/Vamos/Vamos.hh new file mode 100644 index 0000000000000000000000000000000000000000..d4e6b3dc0372b2e92ccef1261623cb02a75cc9fc --- /dev/null +++ b/NPSimulation/Detectors/Vamos/Vamos.hh @@ -0,0 +1,172 @@ +#ifndef Vamos_h +#define Vamos_h +/***************************************************************************** + * Copyright (C) 2009-2018 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: Cyril Lenain contact address: lenain@lpccaen.in2p3.fr * + * * + * Creation Date : octobre 2018 * + * Last update : 09/01/2019 * + *---------------------------------------------------------------------------* + * Decription: * + * This class describe The Vamos spectrometer * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + *****************************************************************************/ + +// C++ header +#include <string> +#include <vector> +using namespace std; + +// G4 headers +#include "G4ThreeVector.hh" +#include "G4RotationMatrix.hh" +#include "G4LogicalVolume.hh" +#include "G4MultiFunctionalDetector.hh" +#include "G4UnionSolid.hh" + +// NPTool header +#include "NPSVDetector.hh" +#include "TVamosData.h" +#include "NPInputParser.h" + +class Vamos : public NPS::VDetector{ + //////////////////////////////////////////////////// + /////// Default Constructor and Destructor ///////// + //////////////////////////////////////////////////// + public: + Vamos() ; + virtual ~Vamos() ; + + //////////////////////////////////////////////////// + /////// Specific Function of this Class /////////// + //////////////////////////////////////////////////// + public: + + void AddVamos(G4double R,double Theta); + void AddBeamCatcher(string Material, G4double Width, double Length, double Thickness, G4ThreeVector Pos); + void AddDetector(G4double Z, string Gas, double Pressure, double Temperature); + void AddIC(G4double Z, double Thickness, string Gas, double Pressure, double Temperature); + + G4LogicalVolume* BuildDC1(); + G4LogicalVolume* BuildDC2(); + + G4LogicalVolume* BuildQuad1(); + G4LogicalVolume* BuildQuad2(); + G4LogicalVolume* BuildDipol(); + + G4LogicalVolume* BuildBeamCatcher(); + G4LogicalVolume* BuildMWPPAC(); + G4LogicalVolume* BuildDC3(); + G4LogicalVolume* BuildDC4(); + G4LogicalVolume* BuildIC(); + + private: + + G4LogicalVolume* m_DC1; + G4LogicalVolume* m_DC2; + + G4LogicalVolume* m_Quad1; + G4LogicalVolume* m_Quad2; + G4LogicalVolume* m_Dipol; + + G4LogicalVolume* m_BeamCatcher; + G4LogicalVolume* m_MWPPAC; + G4LogicalVolume* m_DC3; + G4LogicalVolume* m_DC4; + G4LogicalVolume* m_IC; + + G4double ICcounter; + + //////////////////////////////////////////////////// + ////// Inherite from NPS::VDetector class ///////// + //////////////////////////////////////////////////// + public: + // Read stream at Configfile to pick-up parameters of detector (Position,...) + + void ClearGeometry(); + + // 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() ; + + // Read sensitive part and fill the Root tree. + // Called at in the EventAction::EndOfEventAvtion + void ReadSensitive(const G4Event* event) ; + + public: + // Scorer + // Initialize all Scorer used by Vamos + void InitializeScorers() ; + + // Associated Scorer + G4MultiFunctionalDetector* m_CalorimeterScorer ; + G4MultiFunctionalDetector* m_DCScorer ; + + //////////////////////////////////////////////////// + ///////////Event class to store Data//////////////// + //////////////////////////////////////////////////// + + private: + TVamosData* m_Event ; + + //////////////////////////////////////////////////// + ///////////////Private intern Data////////////////// + //////////////////////////////////////////////////// + + private: + + // Geometry + + G4double R_Catcher = 0; + G4double CatcherWidth = 0; + G4double CatcherLength = 0; + G4double CatcherThickness = 0; + G4double m_R = 0; // distance Target- Entrance of the Mother Volume + G4double m_Theta = 0; + // Detector Coordinate + + string CatcherMaterial; + + G4ThreeVector m_PosCatcher; + vector<G4double> m_Z ; + vector<string> m_Gas; + vector<G4double> m_Pressure; + vector<G4double> m_Temperature; + + vector<G4double> m_ZIC; + vector<G4double> m_ThicknessIC; + vector<G4double> m_PressureIC; + vector<G4double> m_TemperatureIC; + vector<string> m_GasIC; + + // Shape type + + // Visualisation Attribute + G4VisAttributes* m_VisQuad; + G4VisAttributes* m_VisDC; + G4VisAttributes* m_VisVolumeVamos; + G4VisAttributes* m_VisCatcher; + G4VisAttributes* m_VisGasC4H10; + G4VisAttributes* m_VisGasCF4; + + public: + static NPS::VDetector* Construct(); +}; + +#endif diff --git a/NPSimulation/Detectors/beam_dump/cmake_install.cmake b/NPSimulation/Detectors/beam_dump/cmake_install.cmake index be0d02e91221e2f41fea379a06699fee73960373..0974316586c4759d2b0697085160a601e85d6888 100644 --- a/NPSimulation/Detectors/beam_dump/cmake_install.cmake +++ b/NPSimulation/Detectors/beam_dump/cmake_install.cmake @@ -1,4 +1,4 @@ -# Install script for directory: /scratch/nptool/NPSimulation/Detectors/beam_dump +# Install script for directory: /home/lenain/nptool/NPSimulation/Detectors/beam_dump # Set the install prefix if(NOT DEFINED CMAKE_INSTALL_PREFIX) diff --git a/NPSimulation/EventGenerator/CMakeLists.txt b/NPSimulation/EventGenerator/CMakeLists.txt index 3bcc1361a6f2285bca1cbc0c19b6d88e650c4752..5a478e5119b7dfa61806f03673015a48091b008d 100644 --- a/NPSimulation/EventGenerator/CMakeLists.txt +++ b/NPSimulation/EventGenerator/CMakeLists.txt @@ -1,2 +1,3 @@ -add_library(NPSEventGenerator OBJECT EventGeneratorBeam.cc EventGeneratorMultipleParticle.cc EventGeneratorCosmic.cc EventGeneratorIsotropic.cc VEventGenerator.cc) +#add_library(NPSEventGenerator OBJECT EventGeneratorBeam.cc EventGeneratorMultipleParticle.cc EventGeneratorCosmic.cc EventGeneratorIsotropic.cc VEventGenerator.cc) +add_library(NPSEventGenerator OBJECT EventGeneratorBeam.cc EventGeneratorMultipleParticle.cc EventGeneratorIsotropic.cc VEventGenerator.cc) #target_link_libraries(NPSEventGenerator ${ROOT_LIBRARIES} ${Geant4_LIBRARIES} ${NPLib_LIBRARIES} -lNPInitialConditions -lNPInteractionCoordinates ) diff --git a/NPSimulation/Scorers/CMakeLists.txt b/NPSimulation/Scorers/CMakeLists.txt index 83972fed277ca34a8e2963b26912283ad7e971ac..0e285e53a64a9043ddc806641b5918128ce85c81 100644 --- a/NPSimulation/Scorers/CMakeLists.txt +++ b/NPSimulation/Scorers/CMakeLists.txt @@ -1,2 +1,2 @@ -add_library(NPSScorers SHARED NPSHitsMap.hh CalorimeterScorers.cc InteractionScorers.cc DSSDScorers.cc SiliconScorers.cc PhotoDiodeScorers.cc ObsoleteGeneralScorers.cc DriftElectronScorers.cc TPCScorers.cc MDMScorer.cc NeutronDetectorScorers.cc ) +add_library(NPSScorers SHARED DriftChamberScorers.cc NPSHitsMap.hh CalorimeterScorers.cc InteractionScorers.cc DSSDScorers.cc SiliconScorers.cc PhotoDiodeScorers.cc ObsoleteGeneralScorers.cc DriftElectronScorers.cc TPCScorers.cc MDMScorer.cc NeutronDetectorScorers.cc ) target_link_libraries(NPSScorers ${ROOT_LIBRARIES} ${Geant4_LIBRARIES} ${NPLib_LIBRARIES} -lNPInitialConditions -lNPInteractionCoordinates) diff --git a/NPSimulation/Scorers/DriftChamberScorers.cc b/NPSimulation/Scorers/DriftChamberScorers.cc new file mode 100644 index 0000000000000000000000000000000000000000..5073878a1dd5ed6c7362839374454cbd5beacd9a --- /dev/null +++ b/NPSimulation/Scorers/DriftChamberScorers.cc @@ -0,0 +1,128 @@ +/***************************************************************************** + * 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: Cyril Lenain contact adresse: lenain@lpccaen.in2p3.fr * + * * + * Creation Date : June 2018 * + * Last update : 09/01/19 * + *---------------------------------------------------------------------------* + * Decription: * + * Scorer for Drift Chamber based on cathode strips for X position and * + * drift time for Y * + * * + *---------------------------------------------------------------------------* + * Comment: * + *****************************************************************************/ + +#include "DriftChamberScorers.hh" +#include "G4UnitsTable.hh" +#include "G4VPhysicalVolume.hh" +using namespace DriftChamberScorers; + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +unsigned int DCData::CalculateIndex(const vector<unsigned int>& level) { + + unsigned int size = level.size(); + unsigned int result = 0; + unsigned int multiplier = 1; + for (unsigned i = 0; i < size; i++) { + result += level[i] * multiplier; + multiplier *= 1000; + } + return result; +} + +vector<DCData>::iterator DCDataVector::find(const unsigned int& index) { + for (vector<DCData>::iterator it = m_Data.begin(); it != m_Data.end(); it++) { + if ((*it).GetIndex() == index) + return it; + } + return m_Data.end(); +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +PS_DriftChamber::PS_DriftChamber(G4String name, vector<G4int> NestingLevel, + const G4ThreeVector DriftDir, + const double DriftSpeed, G4int depth) + : G4VPrimitiveScorer(name, depth) { + m_NestingLevel = NestingLevel; + m_Dir = DriftDir; + m_DriftSpeed = DriftSpeed; + } + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... +PS_DriftChamber::~PS_DriftChamber() {} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +G4bool PS_DriftChamber::ProcessHits(G4Step* aStep, G4TouchableHistory*) { + + static unsigned int mysize = m_NestingLevel.size(); + + // Extract X and DriftTime + // We use PostStep and PreStep to have at least 2 points + t_PrePosition = aStep->GetPreStepPoint()->GetPosition(); + t_PostPosition = aStep->GetPostStepPoint()->GetPosition(); + G4VSolid* t_solid = aStep->GetPreStepPoint() + ->GetTouchableHandle() + ->GetVolume() + ->GetLogicalVolume() + ->GetSolid(); + + t_postout = t_solid->DistanceToOut(t_PostPosition, -m_Dir); + t_preout = t_solid->DistanceToOut(t_PrePosition, -m_Dir); + t_DriftTime = ((t_preout + t_postout) / m_DriftSpeed) * 0.5; + t_X = (t_PrePosition.x() + t_PostPosition.x()) * 0.5; + + t_Level.clear(); + + for (unsigned int i = 0; i < mysize; i++) + { + t_Level.push_back( + aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber( + m_NestingLevel[i])); } + + vector<DCData>::iterator it; + it = m_Data.find(DCData::CalculateIndex(t_Level)); + if (it != m_Data.end()) { + it->Add(t_DriftTime / microsecond, t_X); + } + + else { + m_Data.Set(t_DriftTime / microsecond, t_X , t_Level); + } + + return TRUE; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +void PS_DriftChamber::Initialize(G4HCofThisEvent*) { clear(); } + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +void PS_DriftChamber::EndOfEvent(G4HCofThisEvent*) { + + // An average is taken + for (vector<DCData>::iterator it = m_Data.begin(); it != m_Data.end(); it++) { + + it->Moy((*it).GetCounter()); + } + +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +void PS_DriftChamber::clear() { + m_Data.clear(); + t_Level.clear(); +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +void PS_DriftChamber::DrawAll() {} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +void PS_DriftChamber::PrintAll() {} +//....oooOO0OOooo.......oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... diff --git a/NPSimulation/Scorers/DriftChamberScorers.hh b/NPSimulation/Scorers/DriftChamberScorers.hh new file mode 100644 index 0000000000000000000000000000000000000000..2c0c7c6755cef1b49481743310855f1e4ccb8e13 --- /dev/null +++ b/NPSimulation/Scorers/DriftChamberScorers.hh @@ -0,0 +1,154 @@ +#ifndef DriftChamberScorers_h +#define DriftChamberScorers_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: Cyril Lenain contact adress: lenain@lpccaen.in2p3.fr * + * * + * Creation Date : July 2018 * + * Last update : 09/01/19 * + *---------------------------------------------------------------------------* + * Decription: * + * Scorer for Drift Chamber based on cathode strips for X position and * + * drift time for Y * + *---------------------------------------------------------------------------* + * Comment: * + *****************************************************************************/ +#include "G4VFastSimulationModel.hh" +#include "G4VPrimitiveScorer.hh" +#include "NPSHitsMap.hh" +#include <map> +using namespace std; +using namespace CLHEP; + +namespace DriftChamberScorers { + + //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + class DCData { + public: + DCData(const double& DriftTime, const double& X, + const vector<unsigned int>& Nesting) { + m_Index = CalculateIndex(Nesting); + m_Level = Nesting; + m_DriftTime = DriftTime; + m_X = X; + m_Counter = 1; + }; + ~DCData(){}; + + private: + unsigned int m_Index; + vector<unsigned int> m_Level; + double m_DriftTime; + double m_X; + double m_Counter; + public: + static unsigned int CalculateIndex(const vector<unsigned int>& Nesting); + + public: + inline double GetIndex() const { return m_Index;} + inline vector<unsigned int> GetLevel() const {return m_Level;}; + inline double GetDriftTime() const { return m_DriftTime; }; + inline double GetX() const { return m_X; }; + inline double GetCounter() {return m_Counter;}; + + public: + void Add(const double& DriftTime, const double& X) { + m_DriftTime += DriftTime; + m_X += X; + m_Counter ++; + }; + + void Moy(const double& Counter){ + m_DriftTime /= Counter; + m_X /= Counter; + }; + }; + + //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + class DCDataVector { + public: + DCDataVector(){}; + ~DCDataVector(){}; + + private: + vector<DCData> m_Data; + + public: + vector<DCData>::iterator find(const unsigned int& index); + inline void clear() { m_Data.clear(); }; + inline vector<DCData>::iterator end() { return m_Data.end(); }; + inline vector<DCData>::iterator begin() { return m_Data.begin(); }; + inline unsigned int size() { return m_Data.size(); }; + inline void Add(const unsigned int& index, const double& DriftTime, const double& X) { + find(index)->Add(DriftTime, X); + }; + + inline void Moy( const unsigned int& index, const double& Counter){ + find(index)->Moy(Counter); + }; + + void Set(const double& DriftTime, const double& X, + const vector<unsigned int>& Nesting) { + m_Data.push_back(DCData(DriftTime, X, Nesting)); + }; + DCData* operator[](const unsigned int& i) { return &m_Data[i]; }; + }; + + //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + class PS_DriftChamber : public G4VPrimitiveScorer { + + public: // with description + PS_DriftChamber(G4String name,vector<G4int> NestingLevel, const G4ThreeVector DriftDir, + const double DriftSpeed, G4int depth = 0); + ~PS_DriftChamber(); + + protected: // with description + G4bool ProcessHits(G4Step*, G4TouchableHistory*); + + public: + void Initialize(G4HCofThisEvent*); + void EndOfEvent(G4HCofThisEvent*); + void clear(); + void DrawAll(); + void PrintAll(); + + private: + vector<G4int> m_NestingLevel; + G4int m_Index; + + + private: // Needed for intermediate calculation (avoid multiple instantiation + // in Processing Hit) + DCDataVector m_Data; + vector<unsigned int> t_Level; + unsigned int t_Detector; + G4ThreeVector m_Dir; + double m_DriftSpeed; + double t_preout; + double t_postout; + double t_DriftTime; + double t_X; + G4ThreeVector t_PrePosition; + G4ThreeVector t_PostPosition; + + public: + inline unsigned int GetMult() { return m_Data.size(); }; + inline double GetDriftTime(const unsigned int& i) { + return m_Data[i]->GetDriftTime(); + }; + inline double GetX(const unsigned int& i) { return m_Data[i]->GetX(); }; + inline double GetCounter(const unsigned int& i) { return m_Data[i]->GetCounter(); }; + inline vector<unsigned int> GetLevel(const unsigned int& i) { + return m_Data[i]->GetLevel(); + }; + }; + +} // namespace DriftChamberScorers + +#endif