From 9353ae97de8a24a546c34cefc14cd4ef4b8c177c Mon Sep 17 00:00:00 2001
From: matta adrien <matta@lpccaen.in2p3.fr>
Date: Wed, 9 Nov 2016 11:24:25 +0100
Subject: [PATCH] * Removing NPLib dependency on libMathMore from Root
- the interpolation is done by a TSpline3 instead of specialized
GSL interpolator - No impact on running time or results
---
NPLib/Detectors/Helios2/THelios2Physics.cxx | 9 ++++-
NPLib/Physics/CMakeLists.txt | 6 ++--
NPLib/Physics/NPEnergyLoss.cxx | 40 ++++++++-------------
NPLib/Physics/NPEnergyLoss.h | 20 +++++------
NPLib/scripts/Style_nptool.C | 4 +--
NPSimulation/Core/ParticleStack.cc | 3 +-
NPSimulation/Detectors/Helios2/Helios2.cc | 23 +++++++-----
NPSimulation/Detectors/Sharc/Sharc.cc | 1 +
NPSimulation/Detectors/Sharc/Sharc.hh | 4 +--
9 files changed, 56 insertions(+), 54 deletions(-)
diff --git a/NPLib/Detectors/Helios2/THelios2Physics.cxx b/NPLib/Detectors/Helios2/THelios2Physics.cxx
index 5306e3462..35352db0f 100644
--- a/NPLib/Detectors/Helios2/THelios2Physics.cxx
+++ b/NPLib/Detectors/Helios2/THelios2Physics.cxx
@@ -201,7 +201,14 @@ double THelios2Physics::ThetaLab(unsigned int i, double m, double B, double q){
double v_par = Z[i]/T;
double v = sqrt(2*(Energy[i]*NPUNITS::MeV)
/(m*NPUNITS::MeV/(NPUNITS::c_squared)));
- return acos(v_par/v);
+ double result = acos(v_par/v);
+
+ // in case ThetaLab is not a number (nan)
+ // return -1000;
+ if(result!=result){
+ result= -1000;
+ }
+ return result;
}
///////////////////////////////////////////////////////////////////////////
diff --git a/NPLib/Physics/CMakeLists.txt b/NPLib/Physics/CMakeLists.txt
index 3d19e9ec7..71610f68c 100644
--- a/NPLib/Physics/CMakeLists.txt
+++ b/NPLib/Physics/CMakeLists.txt
@@ -4,12 +4,12 @@ add_custom_command(OUTPUT TInitialConditionsDict.cxx COMMAND ../scripts/build_di
add_custom_command(OUTPUT TInteractionCoordinatesDict.cxx COMMAND ../scripts/build_dict.sh TInteractionCoordinates.h TInteractionCoordinatesDict.cxx TInteractionCoordinates.rootmap libNPInteractionCoordinates.so DEPENDS TInteractionCoordinates.h)
add_library(NPPhysics SHARED NPBeam.cxx NPEnergyLoss.cxx NPFunction.cxx NPNucleus.cxx NPReaction.cxx NPReactionDict.cxx NPEnergyLossDict.cxx )
-target_link_libraries(NPPhysics ${ROOT_LIBRARIES} MathMore NPCore)
+target_link_libraries(NPPhysics ${ROOT_LIBRARIES} NPCore)
add_library(NPInitialConditions SHARED TInitialConditions.cxx TInitialConditionsDict.cxx )
-target_link_libraries(NPInitialConditions ${ROOT_LIBRARIES} MathMore NPCore)
+target_link_libraries(NPInitialConditions ${ROOT_LIBRARIES} NPCore)
add_library(NPInteractionCoordinates SHARED TInteractionCoordinates.cxx TInteractionCoordinatesDict.cxx)
-target_link_libraries(NPInteractionCoordinates ${ROOT_LIBRARIES} MathMore NPCore)
+target_link_libraries(NPInteractionCoordinates ${ROOT_LIBRARIES} )
install(FILES NPBeam.h NPEnergyLoss.h NPFunction.h NPNucleus.h NPReaction.h TInitialConditions.h TInteractionCoordinates.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
diff --git a/NPLib/Physics/NPEnergyLoss.cxx b/NPLib/Physics/NPEnergyLoss.cxx
index ee6abe535..6f618565a 100644
--- a/NPLib/Physics/NPEnergyLoss.cxx
+++ b/NPLib/Physics/NPEnergyLoss.cxx
@@ -39,8 +39,6 @@
using namespace std;
#include "NPEnergyLoss.h"
-#include "TGraph.h"
-#include "TSpline.h"
#include "TAxis.h"
// NPL
@@ -51,6 +49,7 @@ using namespace NPL;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
EnergyLoss::EnergyLoss() {
fInter = NULL ;
+ fSpline = NULL;
fMax = -10000;
fMin = 1e12;
}
@@ -168,21 +167,14 @@ EnergyLoss::EnergyLoss(string Path , string Source, int NumberOfSlice=100 , int
exit(1);
}
- fInter = new Interpolator( fEnergy , fdEdX_Total ) ;
+ fInter = new TGraph(fEnergy.size(), &fEnergy[0], &fdEdX_Total[0]);
+ fInter->Sort();
+ fSpline = new TSpline3 ("Energy Loss Spline", fInter->GetX(), fInter->GetY(), fInter->GetN());
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void EnergyLoss::Draw() const{
- /*TGraph *gr = new TGraph(fDim, fEnergy, fDedx_Tot);
- gr->Draw("A*");
- gr->GetXaxis()->SetTitle("E (MeV)");
- gr->GetYaxis()->SetTitle("dE/dx (MeV / (mg/cm^{2})");
- gr->Draw("A");
-
- // use a cubic spline to smooth the graph
- TSpline3 *s = new TSpline3("grs",gr) ;
- s->SetLineColor(kRed) ;
- s->Draw("same") ;*/
+void EnergyLoss::Draw(string option) const{
+ fInter->Draw(option.c_str());
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -371,19 +363,17 @@ double EnergyLoss::EvaluateEnergyFromDeltaE( double DeltaE , // Ene
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
double EnergyLoss::Eval(double energy) const {
- if(energy > fMax){ // out of the table, make a first order extrapolation
- return (fInter->Eval(fMax)+(energy-fMax)*fInter->Deriv(fMax));
+ if(energy < 0){
+ cout << "WARNING: negative energy given to NPL::EnergyLoss" << endl;
+ return 0;
}
- else if(energy < 0){
- cout << "WARNING: negative energy given to energy loss" << endl;
- return energy;
+ else if (energy!=energy){
+ cout << "WARNING: nan energy given to NPL::EnergyLoss" << endl;
+ return 0;
}
- else if(energy < fMin){ // out of the table, make a first order extrapolation
- return (fInter->Eval(fMin)+(fMin-energy)*fInter->Deriv(fMin));
- }
- else
- return fInter->Eval(energy);
-
+
+
+ return fInter->Eval(energy,fSpline);
}
diff --git a/NPLib/Physics/NPEnergyLoss.h b/NPLib/Physics/NPEnergyLoss.h
index 84d35a0f1..07612f018 100644
--- a/NPLib/Physics/NPEnergyLoss.h
+++ b/NPLib/Physics/NPEnergyLoss.h
@@ -37,7 +37,8 @@
using namespace std ;
// ROOT
-//#include "TObject.h"
+#include "TGraph.h"
+#include "TSpline.h"
// Use CLHEP System of unit and Physical Constant
#include "NPGlobalSystemOfUnits.h"
@@ -46,17 +47,12 @@ using namespace std ;
using namespace NPUNITS;
#endif
-// ROOT
-#include "Math/InterpolationTypes.h"
-#include "Math/Interpolator.h"
-using namespace ROOT::Math;
-
namespace NPL{
class EnergyLoss{
public : // Constructor
EnergyLoss();
- EnergyLoss( string Path , // Path of dE/dX table file
+ EnergyLoss( string Path , // Path of dE/dX table file
string Source , // Type of file : Geant4,Lise,SRIM
int NumberOfSlice , // Low number = Faster, High Number = more accurate / typical: 100 to 1000
int LiseColumns=0 , // Indicate which model to read in a lise File, set to 0 (Default value) for a SRIM / Geant4 file
@@ -72,11 +68,13 @@ namespace NPL{
vector<double> fdEdX_Nuclear ; // Nuclear Stopping Power
vector<double> fdEdX_Electronic ; // Electronic Stopping Power
vector<double> fdEdX_Total ; // Total Stopping Power
- Interpolator* fInter ; // Interpolator Used to evaluate Energy loss at given energy
-
- double Eval(double ener) const; // return the evaluated energy
+ //Interpolator* fInter ; // Interpolator Used to evaluate Energy loss at given energy
+ TGraph* fInter ; // Graph use to perform interpolation
+ TSpline3* fSpline ; // Spline 3rd order used to perfom the interpolation
+
public : // General Function on dE/dX table
+ double Eval(double ener) const; // return the evaluated energy
double EvaluateNuclearLoss (double ener) const;
double EvaluateElectronicLoss (double ener) const;
double EvaluateTotalLoss (double ener) const;
@@ -124,7 +122,7 @@ namespace NPL{
// Display parameter
void Print() const;
// Draw (CERN ROOT)
- void Draw() const;
+ void Draw(string option="") const;
};
}
diff --git a/NPLib/scripts/Style_nptool.C b/NPLib/scripts/Style_nptool.C
index ccf35f6be..6026c490c 100644
--- a/NPLib/scripts/Style_nptool.C
+++ b/NPLib/scripts/Style_nptool.C
@@ -29,8 +29,8 @@ void Style_nptool(){
// Pad
style->SetPadBottomMargin(0.10);
style->SetPadLeftMargin(0.10);
- style->SetPadTopMargin(0.15);
- style->SetPadRightMargin(0.15);
+ style->SetPadTopMargin(0.05);
+ style->SetPadRightMargin(0.05);
style->SetPadBorderMode(0);
style->SetPadBorderSize(1);
// style->SetPadColor(kWhite);
diff --git a/NPSimulation/Core/ParticleStack.cc b/NPSimulation/Core/ParticleStack.cc
index 704e20333..1d7ac6c17 100644
--- a/NPSimulation/Core/ParticleStack.cc
+++ b/NPSimulation/Core/ParticleStack.cc
@@ -70,7 +70,8 @@ ParticleStack::~ParticleStack(){
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void ParticleStack::AttachInitialConditions(){
// Reasssigned the branch address
- RootOutput::getInstance()->GetTree()->SetBranchAddress("InitialConditions",&m_InitialConditions);
+ if(RootOutput::getInstance()->GetTree()->FindBranch("InitialConditions"))
+ RootOutput::getInstance()->GetTree()->SetBranchAddress("InitialConditions",&m_InitialConditions);
}
diff --git a/NPSimulation/Detectors/Helios2/Helios2.cc b/NPSimulation/Detectors/Helios2/Helios2.cc
index 06a4c0490..d880bb859 100644
--- a/NPSimulation/Detectors/Helios2/Helios2.cc
+++ b/NPSimulation/Detectors/Helios2/Helios2.cc
@@ -63,8 +63,8 @@ namespace Helios2_NS{
// Energy and time Resolution
const double EnergyThreshold = 100*keV;
const double ResoTime = 1*ns ;
- const double ResoEnergy = 21.27*keV ;
-
+ const double ResoEnergyFront = 50*keV ;
+ const double ResoEnergyBack = 24*keV;
const double MagnetInnerRadius = 46*cm;
const double MagnetOutterRadius = 1*m;
const double MagnetLength = 2.35*m;
@@ -76,6 +76,7 @@ namespace Helios2_NS{
const double WaferLength = 56 *mm ;
const double ActiveWaferWidth = 9*mm ;
const double ActiveWaferLength = 50.5 *mm ;
+ const double AluThicness = 0.3*micrometer;
const double WaferThickness = 700* micrometer;
}
@@ -132,14 +133,16 @@ G4LogicalVolume* Helios2::BuildSquareTube(){
G4LogicalVolume* Helios2::BuildSiliconWafer(){
if(!m_SiliconWafer){
G4Box* box1 = new G4Box("Helios2_Box1",Helios2_NS::WaferWidth*0.5,
- Helios2_NS::WaferThickness*0.5,Helios2_NS::WaferLength*0.5);
+ Helios2_NS::WaferThickness*0.5+Helios2_NS::AluThicness,Helios2_NS::WaferLength*0.5);
G4Box* box2 = new G4Box("Helios2_Box2",Helios2_NS::ActiveWaferWidth*0.5,
Helios2_NS::WaferThickness*0.5,Helios2_NS::ActiveWaferLength*0.5);
G4Material* Si= MaterialManager::getInstance()->GetMaterialFromLibrary("Si");
- m_SiliconWafer= new G4LogicalVolume(box1,Si,"logic_Helios2_Wafer",0,0,0);
+ G4Material* Al= MaterialManager::getInstance()->GetMaterialFromLibrary("Al");
+
+ m_SiliconWafer= new G4LogicalVolume(box1,Al,"logic_Helios2_Wafer",0,0,0);
m_ActiveWafer= new G4LogicalVolume(box2,Si,"logic_Helios2_ActiveWafer",0,0,0);
G4ThreeVector AWPos(0,0,0);
@@ -307,6 +310,10 @@ void Helios2::ConstructDetector(G4LogicalVolume* world){
fieldMgr->CreateChordFinder(magField);
BuildMagnet()->SetFieldManager(fieldMgr,true);
+ fieldMgr->SetMinimumEpsilonStep( 1*mm);
+ fieldMgr->SetMaximumEpsilonStep( 10*m );
+ fieldMgr->SetDeltaOneStep( 1 * mm );
+
// Place detectors and support inside it
for (unsigned short i = 0 ; i < m_Z.size() ; i++) {
G4ThreeVector DetPos;
@@ -384,11 +391,11 @@ void Helios2::ReadSensitive(const G4Event* event){
// Loop on the Resistive map
for (Resistive_itr = ResistiveHitMap->GetMap()->begin() ; Resistive_itr != ResistiveHitMap->GetMap()->end() ; Resistive_itr++){
G4double* Info = *(Resistive_itr->second);
- double EBack = RandGauss::shoot(Info[0]+Info[1],Helios2_NS::ResoEnergy);
- double TBack = RandGauss::shoot(Info[2],Helios2_NS::ResoEnergy);
- double EUp = RandGauss::shoot(Info[1],Helios2_NS::ResoEnergy);
+ double EBack = RandGauss::shoot(Info[0]+Info[1],Helios2_NS::ResoEnergyBack);
+ double TBack = RandGauss::shoot(Info[2],Helios2_NS::ResoTime);
+ double EUp = RandGauss::shoot(Info[1],Helios2_NS::ResoEnergyFront);
double TUp = RandGauss::shoot(Info[2],Helios2_NS::ResoTime);
- double EDw = RandGauss::shoot(Info[0],Helios2_NS::ResoEnergy);
+ double EDw = RandGauss::shoot(Info[0],Helios2_NS::ResoEnergyFront);
double TDw = RandGauss::shoot(Info[2],Helios2_NS::ResoTime);
if(EBack>Helios2_NS::EnergyThreshold){
diff --git a/NPSimulation/Detectors/Sharc/Sharc.cc b/NPSimulation/Detectors/Sharc/Sharc.cc
index 96270641d..31c1fcc8e 100644
--- a/NPSimulation/Detectors/Sharc/Sharc.cc
+++ b/NPSimulation/Detectors/Sharc/Sharc.cc
@@ -776,6 +776,7 @@ void Sharc::ReadSensitive(const G4Event* event){
m_Event->SetBack_DetectorNbr(DetNbr);
m_Event->SetBack_StripNbr(BOX_Wafer_Back_NumberOfStrip-StripBack+1);
+
m_Event->SetBack_Energy(RandGauss::shoot(Energy, ResoEnergy));
m_Event->SetBack_TimeCFD(RandGauss::shoot(Time, ResoTime));
m_Event->SetBack_TimeLED(RandGauss::shoot(Time, ResoTime));
diff --git a/NPSimulation/Detectors/Sharc/Sharc.hh b/NPSimulation/Detectors/Sharc/Sharc.hh
index ec2f33d0e..9e2f0d3ea 100644
--- a/NPSimulation/Detectors/Sharc/Sharc.hh
+++ b/NPSimulation/Detectors/Sharc/Sharc.hh
@@ -46,9 +46,7 @@ using namespace CLHEP;
namespace SHARC{
// Energy and time Resolution
const G4double ResoTime = 0 ;
- //const G4double ResoEnergy = 0.035*MeV ;// = zzkeV of Resolution // Unit is MeV/2.35
- const G4double ResoEnergy = 0.042*MeV ;// = zzkeV of Resolution // Unit is MeV/2.35
-
+ const G4double ResoEnergy = 0.035*MeV ;// = zzkeV of Resolution // Unit is MeV/2.35
const G4double EnergyThreshold = 0.1*MeV;
// Geometry
--
GitLab