Update Isotropic source to take into account energy distriubtion from

histogram

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 

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;

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)  ;

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