diff --git a/Inputs/EventGenerator/neutron.source b/Inputs/EventGenerator/neutron.source index bcfd9b541581c4f146a69bdd70853ca3bc26e7f9..4523cdc9a07541c2eb57f0a92c08d71cb391ed5f 100644 --- a/Inputs/EventGenerator/neutron.source +++ b/Inputs/EventGenerator/neutron.source @@ -7,11 +7,11 @@ Isotropic EnergyLow= 0.1 EnergyHigh= 12 %EnergyDistribution= flat - EnergyDistribution= Watt + %EnergyDistribution= Watt %EnergyDistribution= 0.38*sqrt(x)*exp(-x/0.847212) - %EnergyDistribution= -0.00372440431*pow(x,6)+0.387617479*pow(x,5)-14.3752948*pow(x,4)+225.888082*pow(x,3)-1555.60583*pow(x,2)+7983.24902*pow(x,1)+9069.96435 - %EnergyDistribution= 0.619676*TMath::SinH(sqrt(1.07777*x))*exp(-0.847212*x) %EnergyDistribution= 1.5*TMath::SinH(sqrt(1.3*x))*exp(-0.89*x) + EnergyDistribution= FromHisto + EnergyDistributionHist= hEx.root hEx HalfOpenAngleMin= 0 HalfOpenAngleMax= 180 x0= 0 diff --git a/NPSimulation/Detectors/PISTA/PISTA.cc b/NPSimulation/Detectors/PISTA/PISTA.cc index 4e985ebe1cc97fd03dccb86c6a7123e03c7e7cd9..b37da3d47bb530bb5ac28cde0a30737378c292bb 100644 --- a/NPSimulation/Detectors/PISTA/PISTA.cc +++ b/NPSimulation/Detectors/PISTA/PISTA.cc @@ -68,7 +68,7 @@ namespace PISTA_NS{ const double TrapezoidBaseSmall = 41.0*mm; const double TrapezoidHeight = 57.7*mm; const double TrapezoidLength = 1*cm; - const double FirstStageThickness = 100*um; + const double FirstStageThickness = 200*um; const double SecondStageThickness = 1.5*mm; const double DistanceBetweenSi = 4*mm; const double FirstStageNbrOfStrips = 91; diff --git a/NPSimulation/EventGenerator/EventGeneratorIsotropic.cc b/NPSimulation/EventGenerator/EventGeneratorIsotropic.cc index c6d56afc903394e191dec0e79b28074dda811b82..863aa6687e2255b746bba1ed8aa4111ea493a4b2 100644 --- a/NPSimulation/EventGenerator/EventGeneratorIsotropic.cc +++ b/NPSimulation/EventGenerator/EventGeneratorIsotropic.cc @@ -36,6 +36,7 @@ #include "RootOutput.h" #include "NPNucleus.h" #include "NPOptionManager.h" +#include "NPFunction.h" using namespace CLHEP; @@ -43,6 +44,8 @@ using namespace CLHEP; EventGeneratorIsotropic::EventGeneratorIsotropic(){ m_ParticleStack = ParticleStack::getInstance(); event_ID=0; + + m_EnergyDistributionHist = NULL; } @@ -82,6 +85,10 @@ void EventGeneratorIsotropic::ReadConfiguration(NPL::InputParser parser){ it->m_EnergyHigh =blocks[i]->GetDouble("EnergyHigh","MeV"); if(blocks[i]->HasToken("EnergyDistribution")) it->m_EnergyDistribution=blocks[i]->GetString("EnergyDistribution"); + if(blocks[i]->HasToken("EnergyDistributionHist")){ + vector<string> file = blocks[i]->GetVectorString("EnergyDistributionHist"); + m_EnergyDistributionHist = NPL::Read1DProfile(file[0],file[1]); + } it->m_HalfOpenAngleMin =blocks[i]->GetDouble("HalfOpenAngleMin","deg"); it->m_HalfOpenAngleMax =blocks[i]->GetDouble("HalfOpenAngleMax","deg"); it->m_x0 =blocks[i]->GetDouble("x0","mm"); @@ -165,14 +172,16 @@ void EventGeneratorIsotropic::GenerateEvent(G4Event*){ particle_energy = par.m_EnergyLow + RandFlat::shoot() * (par.m_EnergyHigh - par.m_EnergyLow) ; event_ID++; } - else if(par.m_EnergyDistribution=="Watt"){ - particle_energy = fEnergyDist->GetRandom(); + else if(par.m_EnergyDistribution=="FromHisto"){ + if(m_EnergyDistributionHist) + particle_energy = m_EnergyDistributionHist->GetRandom(); } else{ particle_energy = fEnergyDist->GetRandom(); } + // Direction of particle, energy and laboratory angle G4double momentum_x = sin(theta) * cos(phi) ; G4double momentum_y = sin(theta) * sin(phi) ; diff --git a/NPSimulation/EventGenerator/EventGeneratorIsotropic.hh b/NPSimulation/EventGenerator/EventGeneratorIsotropic.hh index b6c2db9a4a19ee7b83606fc8d9a342b59cae52d9..8c63b026d25d0cb9db77f00f510916c04917892b 100644 --- a/NPSimulation/EventGenerator/EventGeneratorIsotropic.hh +++ b/NPSimulation/EventGenerator/EventGeneratorIsotropic.hh @@ -39,6 +39,7 @@ using namespace CLHEP; // ROOT headers #include "TString.h" #include "TF1.h" +#include "TH1F.h" class EventGeneratorIsotropic : public NPS::VEventGenerator{ public: // Constructor and destructor @@ -56,7 +57,7 @@ private: // Source parameter from input file SourceParameters() ; G4double m_EnergyLow ; // Lower limit of energy range G4double m_EnergyHigh ; // Upper limit of energy range - TString m_EnergyDistribution; + TString m_EnergyDistribution; G4double m_HalfOpenAngleMin ; // Min Half open angle of the source G4double m_HalfOpenAngleMax ; // Max Half open angle of the source G4double m_x0 ; // Vertex Position X @@ -72,5 +73,6 @@ private: // Source parameter from input file vector<SourceParameters> m_Parameters ; ParticleStack* m_ParticleStack ; TF1* fEnergyDist; + TH1F* m_EnergyDistributionHist; }; #endif