diff --git a/NPAnalysis/Hira2/Analysis.cxx b/NPAnalysis/Hira2/Analysis.cxx new file mode 100644 index 0000000000000000000000000000000000000000..3b566c8fb55fde8eef3987620d11515e2f6fbc22 --- /dev/null +++ b/NPAnalysis/Hira2/Analysis.cxx @@ -0,0 +1,197 @@ +/***************************************************************************** + * 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<iostream> +using namespace std; +#include"Analysis.h" +#include"NPAnalysisFactory.h" +#include"NPDetectorManager.h" +#include"NPOptionManager.h" +#include"RootOutput.h" +#include"RootInput.h" +//////////////////////////////////////////////////////////////////////////////// +Analysis::Analysis(){ +} +//////////////////////////////////////////////////////////////////////////////// +Analysis::~Analysis(){ +} + +//////////////////////////////////////////////////////////////////////////////// +void Analysis::Init(){ + Hira = (THiraPhysics*) m_DetectorManager->GetDetector("HIRAArray"); + InitialConditions=new TInitialConditions(); + InitOutputBranch(); + InitInputBranch(); + Rand = TRandom3(); + TransferReaction= new NPL::Reaction(); + TransferReaction->ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile()); + + E_ThinSi = 0; + E_ThickSi = 0; + E_CsI = 0; + ELab = 0; + ThetaLab = 0; + PhiLab = 0; + ThetaLab_simu = 0; + ThetaCM = 0; + ExcitationEnergy = 0; + X = 0; + Y = 0; + Z = 0; + TelescopeNumber = 0; + EnergyThreshold = 1; + + TargetThickness = m_DetectorManager->GetTargetThickness()*micrometer; + // Energy loss table: the G4Table are generated by the simulation + He3CD2 = EnergyLoss("He3_CD2.G4table","G4Table",100 ); + He3Al = EnergyLoss("He3_Al.G4table","G4Table",10); + He3Si = EnergyLoss("He3_Si.G4table","G4Table",10); + //Li11CD2 = EnergyLoss("Li11[0.0]_CD2.G4table","G4Table",100); +} + +//////////////////////////////////////////////////////////////////////////////// +void Analysis::TreatEvent(){ + // Reinitiate calculated variable + ReInitValue(); + + + //double BeamEnergy = Rand.Gaus(Initial->GetIncidentInitialKineticEnergy(),4.5); + double BeamEnergy = InitialConditions->GetIncidentInitialKineticEnergy(); + //BeamEnergy = Li11CD2.Slow(BeamEnergy,TargetThickness/2.,0); + + TransferReaction->SetBeamEnergy(BeamEnergy); + //////////////////////////// LOOP on Hira Hit ////////////////// + if(Hira -> ThickSi_E.size() == 1){ + for(unsigned int countHira = 0 ; countHira < Hira->ThickSi_E.size() ; countHira++){ + TelescopeNumber = Hira->TelescopeNumber[countHira]; + + TargetThickness = m_DetectorManager->GetTargetThickness(); + + X = Hira->GetPositionOfInteraction(0).X(); + Y = Hira->GetPositionOfInteraction(0).Y(); + Z = Hira->GetPositionOfInteraction(0).Z(); + + TVector3 PositionOnHira = TVector3(X,Y,Z); + TVector3 ZUnit = TVector3(0,0,1); + + double X_target = InitialConditions->GetIncidentPositionX(); + double Y_target = InitialConditions->GetIncidentPositionY(); + double Z_target = InitialConditions->GetIncidentPositionZ(); + + TVector3 PositionOnTarget = TVector3(X_target,Y_target,Z_target); + TVector3 HitDirection = PositionOnHira - PositionOnTarget; + TVector3 HitDirectionUnit = HitDirection.Unit(); + + TVector3 BeamDirection = InitialConditions->GetBeamDirection(); + double XBeam = BeamDirection.X(); + double YBeam = BeamDirection.Y(); + double ZBeam = BeamDirection.Z(); + + ThetaLab = BeamDirection.Angle(HitDirection); + ThetaLab = ThetaLab/deg; + + PhiLab = HitDirection.Phi()/deg; + + E_ThickSi = Hira->ThickSi_E[countHira]; + E_ThinSi = Hira->ThinSi_E[countHira]; + //ELab = E_ThinSi + E_ThickSi; + if(Hira->CsI_E.size() == 1){ + for(int countCsI =0; countCsI<Hira->CsI_E.size(); countCsI++){ + E_CsI = Hira->CsI_E[countCsI]; + //ELab += E_CsI; + if(E_CsI>EnergyThreshold) ELab = E_ThinSi + E_ThickSi + E_CsI; + } + } + else ELab = -1000; + + // ********************** Angle in the CM frame ***************************** + TransferReaction -> SetNuclei3(ELab, ThetaLab*deg); + ThetaCM = TransferReaction->GetThetaCM()/deg; + ExcitationEnergy = TransferReaction->GetExcitation4(); + }//end loop Hira + }//end if Hira +} + +//////////////////////////////////////////////////////////////////////////////// +void Analysis::End(){ +} + +//////////////////////////////////////////////////////////////////////////////// +void Analysis::InitOutputBranch() { + RootOutput::getInstance()->GetTree()->Branch( "ThetaLab" , &ThetaLab , "ThetaLab/D" ); + RootOutput::getInstance()->GetTree()->Branch( "PhiLab" , &PhiLab , "PhiLab/D" ) ; + RootOutput::getInstance()->GetTree()->Branch("ThetaCM", &ThetaCM,"ThetaCM/D") ; + RootOutput::getInstance()->GetTree()->Branch( "E_ThinSi" , &E_ThinSi , "E_ThinSi/D" ) ; + RootOutput::getInstance()->GetTree()->Branch( "E_ThickSi" , &E_ThickSi , "E_ThickSi/D" ) ; + RootOutput::getInstance()->GetTree()->Branch( "E_CsI" , &E_CsI , "E_CsI/D" ) ; + RootOutput::getInstance()->GetTree()->Branch( "ELab" , &ELab , "ELab/D" ) ; + RootOutput::getInstance()->GetTree()->Branch("ExcitationEnergy", &ExcitationEnergy,"ExcitationEnergy/D") ; + RootOutput::getInstance()->GetTree()->Branch( "X" , &X , "X/D" ) ; + RootOutput::getInstance()->GetTree()->Branch( "Y" , &Y , "Y/D" ) ; + RootOutput::getInstance()->GetTree()->Branch( "Z" , &Z , "Z/D" ) ; + RootOutput::getInstance()->GetTree()->Branch( "TelescopeNumber" , &TelescopeNumber , "TelescopeNumber/D" ) ; + RootOutput::getInstance()->GetTree()->Branch("InteractionCoordinates","TInteractionCoordinates",&InteractionCoordinates); + RootOutput::getInstance()->GetTree()->Branch("InitialConditions","TInitialConditions",&InitialConditions); +} + +//////////////////////////////////////////////////////////////////////////////// +void Analysis::InitInputBranch(){ + RootInput:: getInstance()->GetChain()->SetBranchStatus("InitialConditions",true ); + RootInput:: getInstance()->GetChain()->SetBranchStatus("fIC_*",true ); + RootInput:: getInstance()->GetChain()->SetBranchAddress("InitialConditions",&InitialConditions); +} + +//////////////////////////////////////////////////////////////////////////////// +void Analysis::ReInitValue(){ + E_ThinSi = -1000; + E_ThickSi = -1000; + E_CsI = -1000; + ELab = -1000; + ThetaLab = -1000; + PhiLab = -1000; + ThetaCM = -1000; + ExcitationEnergy = -1000; + X = -1000; + Y = -1000; + Z = -1000; + TelescopeNumber = -1; +} + +//////////////////////////////////////////////////////////////////////////////// +// Construct Method to be pass to the DetectorFactory // +//////////////////////////////////////////////////////////////////////////////// +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; +} \ No newline at end of file diff --git a/NPAnalysis/Hira2/Analysis.h b/NPAnalysis/Hira2/Analysis.h new file mode 100644 index 0000000000000000000000000000000000000000..6cb0822d419cfa63fdd1cb8c66e5995ed7ed464e --- /dev/null +++ b/NPAnalysis/Hira2/Analysis.h @@ -0,0 +1,77 @@ +#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"NPVAnalysis.h" +#include"THiraPhysics.h" +#include "TInitialConditions.h" +#include "TInteractionCoordinates.h" +#include "NPEnergyLoss.h" +#include "NPReaction.h" +#include "TRandom3.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: + double TargetThickness; + double ExcitationEnergy; + double ELab; + double E_ThinSi; + double E_ThickSi; + double E_CsI; + double PhiLab; + double ThetaLab; + double ThetaLab_simu; + double ThetaCM; + double X,Y,Z; + double TelescopeNumber; + double EnergyThreshold; + + + + NPL::Reaction* TransferReaction; + + // intermediate variable + TRandom3 Rand; + + + NPL::EnergyLoss He3CD2 ; + NPL::EnergyLoss He3Al ; + NPL::EnergyLoss He3Si ; + NPL::EnergyLoss Li11CD2 ; + + THiraPhysics* Hira; + TInitialConditions* InitialConditions; + TInteractionCoordinates* InteractionCoordinates; +}; +#endif \ No newline at end of file diff --git a/NPAnalysis/Hira2/CMakeLists.txt b/NPAnalysis/Hira2/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..ceb5ebea2e999f49b4aae22dc7bcef23eeb6bb2d --- /dev/null +++ b/NPAnalysis/Hira2/CMakeLists.txt @@ -0,0 +1,33 @@ +cmake_minimum_required (VERSION 2.8) +#Finding NPTool +set(NPTOOL "$ENV{NPTOOL}") +set(NPLIB "${NPTOOL}/NPLib") +set(NPTOOL_INCLUDE_DIR "${NPLIB}/include") +set(NPTOOL_LIB_DIR "${NPLIB}/lib") + +include("${NPLIB}/FindROOT.cmake") + +project (NPAnalysis) +set(CMAKE_BUILD_TYPE Release) + + +# Add root to the link and include directories +include_directories( ${ROOT_INCLUDE_DIR}) +link_directories( ${ROOT_LIBRARY_DIR}) +include_directories( ${NPTOOL_INCLUDE_DIR}) +link_directories( ${NPTOOL_LIB_DIR}) + +# Get the compilator flag from root to assure consistancy +EXEC_PROGRAM(${ROOT_CONFIG_EXECUTABLE} + ARGS "--cflags" + OUTPUT_VARIABLE root_cflags ) + +set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${root_cflags}") + +# If the compiler is Clang, silence the unrecognised flags +if(${CMAKE_CXX_COMPILER_ID} MATCHES ".*Clang.*") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Qunused-arguments -undefined dynamic_lookup") +endif() + +add_library(NPAnalysis SHARED Analysis.cxx) +target_link_libraries(NPAnalysis ${ROOT_LIBRARIES} -L${NPLIB}/lib -lNPCore -lNPPhysics -lNPInitialConditions -lNPInteractionCoordinates) diff --git a/NPAnalysis/Hira2/RunToTreat.txt b/NPAnalysis/Hira2/RunToTreat.txt new file mode 100755 index 0000000000000000000000000000000000000000..82908222e590c8760e20a5e47347834912b6c19c --- /dev/null +++ b/NPAnalysis/Hira2/RunToTreat.txt @@ -0,0 +1,5 @@ +TTreeName + SimulatedTree +RootFileName + ../../Outputs/Simulation/Example1.root + diff --git a/NPLib/Sharc/TSharcPhysics.cxx b/NPLib/Sharc/TSharcPhysics.cxx index c62cdbb9e1476751e42b42cb2365e57bdc955b39..03679bb3d31b0a3a9452bce8546c8a427959f4a3 100644 --- a/NPLib/Sharc/TSharcPhysics.cxx +++ b/NPLib/Sharc/TSharcPhysics.cxx @@ -960,14 +960,14 @@ NPL::VDetector* TSharcPhysics::Construct(){ // Registering the construct method to the factory // //////////////////////////////////////////////////////////////////////////////// extern "C"{ -class proxy{ +class proxy_sharc{ public: - proxy(){ + proxy_sharc(){ NPL::DetectorFactory::getInstance()->AddToken("Sharc","Sharc"); NPL::DetectorFactory::getInstance()->AddDetector("Sharc",TSharcPhysics::Construct); } }; -proxy p; +proxy_sharc p; } diff --git a/NPLib/Tigress/TTigressPhysics.cxx b/NPLib/Tigress/TTigressPhysics.cxx index 6ff43e8a50e0c65bae5f2dada533df2ef417ef71..bd10da52313ea9aefc39fe96d5757ca8bfc5221e 100644 --- a/NPLib/Tigress/TTigressPhysics.cxx +++ b/NPLib/Tigress/TTigressPhysics.cxx @@ -267,14 +267,14 @@ NPL::VDetector* TTigressPhysics::Construct(){ // Registering the construct method to the factory // //////////////////////////////////////////////////////////////////////////////// extern "C"{ -class proxy{ +class proxy_tigress{ public: - proxy(){ + proxy_tigress(){ NPL::DetectorFactory::getInstance()->AddToken("Tigress","Tigress"); NPL::DetectorFactory::getInstance()->AddDetector("Tigress",TTigressPhysics::Construct); } }; -proxy p; +proxy_tigress p; } diff --git a/NPLib/Trifoil/TTrifoilPhysics.cxx b/NPLib/Trifoil/TTrifoilPhysics.cxx index a403f851a484e08d43129d39beb8b62868951cd7..78ed7f32a4142baabbb05a323a191c71784a6465 100644 --- a/NPLib/Trifoil/TTrifoilPhysics.cxx +++ b/NPLib/Trifoil/TTrifoilPhysics.cxx @@ -123,14 +123,14 @@ NPL::VDetector* TTrifoilPhysics::Construct(){ // Registering the construct method to the factory // //////////////////////////////////////////////////////////////////////////////// extern "C"{ -class proxy{ +class proxy_trifoil{ public: - proxy(){ - NPL::DetectorFactory::getInstance()->AddToken("Trifoil","Trifoil"); + proxy_trifoil(){ + NPL::DetectorFactory::getInstance()->AddToken("Trifoil","Trifoil"); NPL::DetectorFactory::getInstance()->AddDetector("Trifoil",TTrifoilPhysics::Construct); - } + } }; -proxy p; +proxy_trifoil p; } diff --git a/NPSimulation/Core/PhysicsList.cc b/NPSimulation/Core/PhysicsList.cc index 2f9098136c0a97f10231175247e3de29ec067cea..dfc439efc167bc879682745900103c271876a6cf 100644 --- a/NPSimulation/Core/PhysicsList.cc +++ b/NPSimulation/Core/PhysicsList.cc @@ -239,6 +239,20 @@ void PhysicsList::ConstructEM(){ opt.SetMaxEnergy(50.*GeV) ; opt.SetDEDXBinning(5000) ; opt.SetLambdaBinning(5000) ; + + //energy loss + opt.SetLinearLossLimit(1.e-9); + opt.SetStepFunction(0.001, 10.*um); + + //Multiple scattering + opt.SetMscLateralDisplacement(true); + opt.SetLossFluctuations(true); + //opt.SetMscStepLimitation(fMinimal); + //opt.SetMscStepLimitation(fUseDistanceToBoundary); + opt.SetMscStepLimitation(fUseSafety); + + //ionization + opt.SetSubCutoff(false); } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... @@ -251,9 +265,11 @@ void PhysicsList::SetCuts(){ SetCutsWithDefault(); // for gamma-rays - // SetCutValue(0.01*mm, "gamma"); - // SetCutValue(0.01*mm, "e-"); - // SetCutValue(0.01*mm, "e+"); + /* G4double em_cuts = 0.01*mm; // Use 10cm to avoid generating delta-rays + //G4double em_cuts = 10.*cm; // Use 10cm to avoid generating delta-rays + SetCutValue(em_cuts, "gamma"); + SetCutValue(em_cuts, "e-"); + SetCutValue(em_cuts, "e+");*/ // Retrieve verbose level SetVerboseLevel(temp);