From f1eb934f9fdd89cf7e18077f7a29b6af5dd8c32b Mon Sep 17 00:00:00 2001
From: adrien-matta <a.matta@surrey.ac.uk>
Date: Tue, 20 Oct 2015 12:00:05 +0100
Subject: [PATCH] * Adding support for neutron beam * adding Uranium target for
fission exp * adding timestamp to Trifoil
---
NPAnalysis/S1554/Analysis.cxx | 150 +++++++++++++++---------
NPAnalysis/S1554/Analysis.h | 3 +-
NPLib/Physics/nubtab03.asc | 2 +-
NPLib/Trifoil/TTrifoilData.cxx | 2 +
NPLib/Trifoil/TTrifoilData.h | 7 +-
NPLib/Trifoil/TTrifoilPhysics.cxx | 11 +-
NPLib/Trifoil/TTrifoilPhysics.h | 5 +-
NPSimulation/Core/EventGeneratorBeam.cc | 8 +-
NPSimulation/Core/MaterialManager.cc | 7 ++
NPSimulation/Core/PhysicsList.cc | 114 ++++++++++--------
NPSimulation/Core/Target.cc | 4 +
11 files changed, 196 insertions(+), 117 deletions(-)
diff --git a/NPAnalysis/S1554/Analysis.cxx b/NPAnalysis/S1554/Analysis.cxx
index 0942b870b..a0718bfc5 100644
--- a/NPAnalysis/S1554/Analysis.cxx
+++ b/NPAnalysis/S1554/Analysis.cxx
@@ -43,13 +43,13 @@ void Analysis::Init(){
LightCD2 = EnergyLoss("deuteron_CD2.G4table","G4Table",10);
LightAl = EnergyLoss("deuteron_Al.G4table","G4Table",10);
-// LightAl = EnergyLoss("alpha_Al.G4table","G4Table",10);
+ // LightAl = EnergyLoss("alpha_Al.G4table","G4Table",10);
BeamCD2 = EnergyLoss("Si28_CD2.G4table","G4Table",10);
myReaction = new NPL::Reaction();
myReaction->ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
TargetThickness = m_DetectorManager->GetTargetThickness()*micrometer;
-// TargetThickness = 0;
+ // TargetThickness = 0;
OriginalBeamEnergy = myReaction->GetBeamEnergy();
Rand = TRandom3();
DetectorNumber = 0 ;
@@ -72,31 +72,43 @@ void Analysis::Init(){
myReaction->SetBeamEnergy(finalEnergy);
cout << "Set Beam energy to: " << finalEnergy << " MeV" << endl;
- // Load cut for PAD cal
- TFile* file ;
- file = new TFile("cuts/deuton_d10.root","READ");
- cut_deuton_d10 = (TCutG*) file->FindObjectAny("deuton_d10");
- file = new TFile("cuts/proton_d10.root","READ");
- cut_proton_d10 = (TCutG*) file->FindObjectAny("proton_d10");
- file = new TFile("cuts/deuton_d12.root","READ");
- cut_deuton_d12 = (TCutG*) file->FindObjectAny("deuton_d12");
- file = new TFile("cuts/proton_d12.root","READ");
- cut_proton_d12 = (TCutG*) file->FindObjectAny("proton_d12");
+ // Load cut for second Pad cal
+/* TFile* file ;
+ file = new TFile("cuts/cut_elgs_det10.root","READ");
+ cut_deuton_d10 = (TCutG*) file->FindObjectAny("cut_elgs_det10");
+ file = new TFile("cuts/cut_elgs_det12.root","READ");
+ cut_deuton_d12 = (TCutG*) file->FindObjectAny("cut_elgs_det12");
+ file = new TFile("cuts/Cut_dd.root","READ");
+ cut_dd = (TCutG*) file->FindObjectAny("Cut_dd");
PADFile.open("PADEvent.txt");
+*/
+
-/* // Load the cut for alignement
- for(unsigned int i = 0 ; i < 8 ; i++){
- TFile* file = new TFile(Form("cuts/Ex5_D%i.root",i+1),"READ");
- cut_ex5.push_back((TCutG*) file->FindObjectAny(Form("Ex5_D%i",i+1)));
- if(! file->FindObjectAny(Form("Ex5_D%i",i+1))){
+ // Load cut for PAD cal
+/* TFile* file ;
+ file = new TFile("cuts/deuton_d10.root","READ");
+ cut_deuton_d10 = (TCutG*) file->FindObjectAny("deuton_d10");
+ file = new TFile("cuts/proton_d10.root","READ");
+ cut_proton_d10 = (TCutG*) file->FindObjectAny("proton_d10");
+ file = new TFile("cuts/deuton_d12.root","READ");
+ cut_deuton_d12 = (TCutG*) file->FindObjectAny("deuton_d12");
+ file = new TFile("cuts/proton_d12.root","READ");
+ cut_proton_d12 = (TCutG*) file->FindObjectAny("proton_d12");
+ PADFile.open("PADEvent.txt");
+ */
+ /* // Load the cut for alignement
+ for(unsigned int i = 0 ; i < 8 ; i++){
+ TFile* file = new TFile(Form("cuts/Ex5_D%i.root",i+1),"READ");
+ cut_ex5.push_back((TCutG*) file->FindObjectAny(Form("Ex5_D%i",i+1)));
+ if(! file->FindObjectAny(Form("Ex5_D%i",i+1))){
cout << "Fail to Load " << Form("Ex5_D%i",i+1) << endl;
exit(1);
- }
- }
+ }
+ }
- box_pos.open("BoxPos.txt");
- qqq_pos.open("QQQPos.txt");
-*/
+ box_pos.open("BoxPos.txt");
+ qqq_pos.open("QQQPos.txt");
+ */
}
////////////////////////////////////////////////////////////////////////////////
@@ -113,7 +125,8 @@ void Analysis::TreatEvent(){
TVector3 HitDirection = Sharc -> GetPositionOfInteraction(0)-TargetPosition;
ThetaLab = HitDirection.Angle( BeamDirection );
ThetaNormalTarget = HitDirection.Angle( TVector3(0,0,1) ) ;
- ThetaDetector = HitDirection.Angle(Sharc->GetDetectorNormal(0));
+ ThetaDetector = HitDirection.Angle(-Sharc->GetDetectorNormal(0));
+ Distance = HitDirection.Mag();
/************************************************/
/************************************************/
@@ -129,46 +142,72 @@ void Analysis::TreatEvent(){
// Target Correction
ELab = LightCD2.EvaluateInitialEnergy( Energy ,TargetThickness*0.5, ThetaNormalTarget);
-
+
/************************************************/
- int DetectorNumber = Sharc->DetectorNumber[0];
- if(DetectorNumber==10 && Sharc->PAD_E[0]>0){
- if(cut_proton_d10->IsInside(Sharc->PAD_E[0],-Sharc->Strip_E[0]*cos(ThetaDetector)))
- PADFile << "10 proton " << Sharc->Strip_E[0] << " " << Sharc->PAD_E[0] << " " << ThetaDetector << endl;
- else if(cut_deuton_d10->IsInside(Sharc->PAD_E[0],-Sharc->Strip_E[0]*cos(ThetaDetector)))
- PADFile << "10 deuton " << Sharc->Strip_E[0] << " " << Sharc->PAD_E[0] << " " << ThetaDetector << endl;
- }
+ /* int DetectorNumber = Sharc->DetectorNumber[0];
+ if(DetectorNumber==10 && Sharc->PAD_E[0]>0){
+ if(cut_proton_d10->IsInside(Sharc->PAD_E[0],Sharc->Strip_E[0]*cos(ThetaDetector))){
+ PADFile << "10 proton " << Sharc->Strip_E[0] << " " << Sharc->PAD_E[0] << " " << ThetaDetector << endl;
+ }
+ else if(cut_deuton_d10->IsInside(Sharc->PAD_E[0],Sharc->Strip_E[0]*cos(ThetaDetector)))
+ PADFile << "10 deuton " << Sharc->Strip_E[0] << " " << Sharc->PAD_E[0] << " " << ThetaDetector << endl;
+ }
+
+ if(DetectorNumber==12 && Sharc->PAD_E[0]>0){
+ if(cut_proton_d12->IsInside(Sharc->PAD_E[0],Sharc->Strip_E[0]*cos(ThetaDetector)))
+ PADFile << "12 proton " << Sharc->Strip_E[0] << " " << Sharc->PAD_E[0] << " " << ThetaDetector << endl;
+ else if(cut_deuton_d12->IsInside(Sharc->PAD_E[0], Sharc->Strip_E[0]*cos(ThetaDetector)))
+ PADFile << "12 deuton " << Sharc->Strip_E[0] << " " << Sharc->PAD_E[0] << " " << ThetaDetector << endl;
+ }
+ */
- if(DetectorNumber==12 && Sharc->PAD_E[0]>0){
- if(cut_proton_d12->IsInside(Sharc->PAD_E[0],-Sharc->Strip_E[0]*cos(ThetaDetector)))
- PADFile << "12 proton " << Sharc->Strip_E[0] << " " << Sharc->PAD_E[0] << " " << ThetaDetector << endl;
- else if(cut_deuton_d12->IsInside(Sharc->PAD_E[0], -Sharc->Strip_E[0]*cos(ThetaDetector)))
- PADFile << "12 deuton " << Sharc->Strip_E[0] << " " << Sharc->PAD_E[0] << " " << ThetaDetector << endl;
- }
+ /************************************************/
+ /* int DetectorNumber = Sharc->DetectorNumber[0];
+ if(cut_dd->IsInside(Sharc->PAD_E[0],Sharc->Strip_E[0]*cos(ThetaDetector))){
+ if(DetectorNumber == 10){
+ if(cut_deuton_d10->IsInside(ThetaLab/deg,ELab)){
+ double a = 9.23351e-13 ;
+ double b = 8.35809e-06 ;
+ double c = 5.77842e-01-Sharc->PAD_E[0];
+ double delta = b*b-4*a*c;
+ double RawE = (-b + sqrt(delta))/(2*a);
+ PADFile << "10 deuton " << Sharc->Strip_E[0] << " " << RawE << " " << ThetaLab << " " << ThetaNormalTarget << endl;
+ }
+ }
+ else if(DetectorNumber == 12){
+ if(cut_deuton_d12->IsInside(ThetaLab/deg,ELab)){
+ double a = 1.63514e-12;
+ double b = 7.52050e-06;
+ double c = 5.93687e-01-Sharc->PAD_E[0];
+ double delta = b*b-4*a*c;
+ double RawE = (-b + sqrt(delta))/(2*a);
+ PADFile << "12 deuton " << Sharc->Strip_E[0] << " " << RawE<< " " << ThetaLab << " " << ThetaNormalTarget << endl;
+ }
+ }
- /************************************************/
-
- /************************************************/
-/* int DetectorNumber = Sharc->DetectorNumber[0];
- bool checkT = false;
- if(Tigress->AddBack_DC.size()>0){
- if(Tigress->AddBack_DC[0]/1000.> 4.6 && Tigress->AddBack_DC[0]/1000.< 5.4)
- checkT=true;
- }
+ }*/
- if(checkT){
- if(DetectorNumber<5){
- if(cut_ex5[DetectorNumber-1]->IsInside(ThetaLab/deg,ELab))
+ /************************************************/
+ /* int DetectorNumber = Sharc->DetectorNumber[0];
+ bool checkT = false;
+ if(Tigress->AddBack_DC.size()>0){
+ if(Tigress->AddBack_DC[0]/1000.> 4.6 && Tigress->AddBack_DC[0]/1000.< 5.4)
+ checkT=true;
+ }
+
+ if(checkT){
+ if(DetectorNumber<5){
+ if(cut_ex5[DetectorNumber-1]->IsInside(ThetaLab/deg,ELab))
qqq_pos << DetectorNumber << " " << Energy << " " << HitDirection.X() << " " << HitDirection.Y() << " " << HitDirection.Z() << endl;
- }
+ }
- else if(DetectorNumber<9){
- if(cut_ex5[DetectorNumber-1]->IsInside(ThetaLab/deg,ELab))
+ else if(DetectorNumber<9){
+ if(cut_ex5[DetectorNumber-1]->IsInside(ThetaLab/deg,ELab))
box_pos << DetectorNumber << " " << Energy << " " << HitDirection.X() << " " << HitDirection.Y() << " " << HitDirection.Z() << endl;
- }
- }
-*/
+ }
+ }
+ */
/************************************************/
/************************************************/
@@ -196,6 +235,7 @@ void Analysis::InitOutputBranch(){
RootOutput::getInstance()->GetTree()->Branch("ThetaLab",&ThetaLab,"ThetaLab/D");
RootOutput::getInstance()->GetTree()->Branch("ThetaCM",&ThetaCM,"ThetaCM/D");
RootOutput::getInstance()->GetTree()->Branch("ThetaDetector",&ThetaDetector,"ThetaDetector/D");
+ RootOutput::getInstance()->GetTree()->Branch("Distance",&Distance,"Distance/D");
RootOutput::getInstance()->GetTree()->Branch("RunNumber",&RunNumber,"RunNumber/I");
RootOutput::getInstance()->GetTree()->Branch("RunNumberMinor",&RunNumberMinor,"RunNumberMinor/I");
diff --git a/NPAnalysis/S1554/Analysis.h b/NPAnalysis/S1554/Analysis.h
index 5b14b830c..22d78fdf0 100644
--- a/NPAnalysis/S1554/Analysis.h
+++ b/NPAnalysis/S1554/Analysis.h
@@ -85,13 +85,14 @@ TInitialConditions* myInit ;
double ThetaDetector;
double Si_E_Sharc ;
double E_Sharc ;
+ double Distance;
TSharcPhysics* Sharc;
TTigressPhysics* Tigress;
std::vector<TCutG*> cut_ex5;
std::ofstream box_pos;
std::ofstream qqq_pos;
-
+ TCutG* cut_dd;
TCutG* cut_deuton_d10;
TCutG* cut_proton_d10;
TCutG* cut_deuton_d12;
diff --git a/NPLib/Physics/nubtab03.asc b/NPLib/Physics/nubtab03.asc
index 596f1bc09..ad13518a4 100644
--- a/NPLib/Physics/nubtab03.asc
+++ b/NPLib/Physics/nubtab03.asc
@@ -1,4 +1,4 @@
-001 0000 1 n 8071.3171 0.0005 613.9 s 0.6 1/2+ 00 02PaDGt B-=100
+001 0000 1n 8071.3171 0.0005 613.9 s 0.6 1/2+ 00 02PaDGt B-=100
001 0010 1H 7288.9705 0.0001 stbl 1/2+ 00 98Ro45d IS=99.9885 70
002 0010 2H 13135.7216 0.0003 stbl 1+ 99 IS=0.0115 70
003 0010 3H 14949.8060 0.0023 12.32 y 0.02 1/2+ 00 B-=100
diff --git a/NPLib/Trifoil/TTrifoilData.cxx b/NPLib/Trifoil/TTrifoilData.cxx
index 2b24e37a1..2830f8a85 100644
--- a/NPLib/Trifoil/TTrifoilData.cxx
+++ b/NPLib/Trifoil/TTrifoilData.cxx
@@ -48,6 +48,8 @@ void TTrifoilData::Clear(){
fTrifoil_Waveform.clear();
fTrifoil_TimeCFD.clear();
fTrifoil_TimeLED.clear();
+ fTrifoil_TimeStamp.clear();
+
}
/////////////////////////
void TTrifoilData::Dump() const{
diff --git a/NPLib/Trifoil/TTrifoilData.h b/NPLib/Trifoil/TTrifoilData.h
index 8f5ebbc78..d114bbdb9 100644
--- a/NPLib/Trifoil/TTrifoilData.h
+++ b/NPLib/Trifoil/TTrifoilData.h
@@ -39,7 +39,8 @@ private:
vector<TH1F> fTrifoil_Waveform;
vector<Double_t> fTrifoil_TimeCFD;
vector<Double_t> fTrifoil_TimeLED;
-
+ vector<Double_t> fTrifoil_TimeStamp;
+
public:
TTrifoilData();
virtual ~TTrifoilData();
@@ -52,11 +53,13 @@ public:
inline void SetWaveform(const TH1F &Waveform) {fTrifoil_Waveform.push_back(Waveform);}
inline void SetTimeCFD(const Double_t &TimeCFD) {fTrifoil_TimeCFD.push_back(TimeCFD);}
inline void SetTimeLED(const Double_t &TimeLED) {fTrifoil_TimeLED.push_back(TimeLED);}
-
+ inline void SetTimeStamp(const Double_t &TimeStamp) {fTrifoil_TimeStamp.push_back(TimeStamp);}
+
///////////////////// GETTERS ////////////////////////
inline TH1F* GetWaveform(const unsigned int &i) {return &(fTrifoil_Waveform[i]);}
inline Double_t GetTimeCFD(const unsigned int &i) {return fTrifoil_TimeCFD[i];}
inline Double_t GetTimeLED(const unsigned int &i) {return fTrifoil_TimeLED[i];}
+ inline Double_t GetTimeStamp(const unsigned int &i) {return fTrifoil_TimeStamp[i];}
inline unsigned int GetMultiplicity() {return fTrifoil_Waveform.size();}
ClassDef(TTrifoilData,1) // TrifoilData structure
diff --git a/NPLib/Trifoil/TTrifoilPhysics.cxx b/NPLib/Trifoil/TTrifoilPhysics.cxx
index e8da764dc..06bb3c5ae 100644
--- a/NPLib/Trifoil/TTrifoilPhysics.cxx
+++ b/NPLib/Trifoil/TTrifoilPhysics.cxx
@@ -52,13 +52,15 @@ void TTrifoilPhysics::BuildSimplePhysicalEvent(){
///////////////////////////////////////////////////////////////////////////
void TTrifoilPhysics::BuildPhysicalEvent(){
unsigned int mysize = m_EventData->GetMultiplicity();
-
+ double base,maxi,mytime;
for (unsigned int i = 0 ; i < mysize ; i++){
TH1F* h = m_EventData->GetWaveform(i);
- double base = h->GetBinContent(h->GetMinimumBin());
- double maxi = h->GetBinContent(h->GetMaximumBin());
+ base = h->GetBinContent(h->GetMinimumBin());
+ maxi = h->GetBinContent(h->GetMaximumBin());
if(maxi>2000 && base>-300){
- Time.push_back(h->GetMaximumBin());
+ mytime = (h->GetMaximumBin());
+ TimeStamp.push_back(m_EventData->GetTimeStamp(i)*10);
+ Time.push_back(mytime);
Energy.push_back(maxi);
}
}
@@ -68,6 +70,7 @@ void TTrifoilPhysics::BuildPhysicalEvent(){
void TTrifoilPhysics::Clear(){
Energy.clear() ;
Time.clear() ;
+ TimeStamp.clear();
}
///////////////////////////////////////////////////////////////////////////
diff --git a/NPLib/Trifoil/TTrifoilPhysics.h b/NPLib/Trifoil/TTrifoilPhysics.h
index ba92c5f75..58ff3b226 100644
--- a/NPLib/Trifoil/TTrifoilPhysics.h
+++ b/NPLib/Trifoil/TTrifoilPhysics.h
@@ -48,9 +48,10 @@ class TTrifoilPhysics : public TObject, public NPL::VDetector
public:
// EventType is True if the wave form analysis return a valid Trifoil event
- vector<double> Energy;
+ vector<double> Energy;
vector<double> Time ;
-
+ vector<double> TimeStamp ;
+
public: // Innherited from VDetector Class
// Read stream at ConfigFile to pick-up parameters of detector (Position,...) using Token
diff --git a/NPSimulation/Core/EventGeneratorBeam.cc b/NPSimulation/Core/EventGeneratorBeam.cc
index d9121520e..b4e0fe3e7 100644
--- a/NPSimulation/Core/EventGeneratorBeam.cc
+++ b/NPSimulation/Core/EventGeneratorBeam.cc
@@ -75,7 +75,13 @@ void EventGeneratorBeam::GenerateEvent(G4Event* anEvent){
if( anEvent->GetEventID()==0){
// Define the particle to be shoot
- m_particle = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(m_Beam->GetZ(), m_Beam->GetA() ,m_Beam->GetExcitationEnergy());
+ if(m_Beam->GetZ()==0 && m_Beam->GetA()==1){
+ m_particle = G4ParticleTable::GetParticleTable()->FindParticle("neutron");
+ }
+
+ else
+ m_particle = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(m_Beam->GetZ(), m_Beam->GetA() ,m_Beam->GetExcitationEnergy());
+
}
///////////////////////////////////////////////////////////////////////
diff --git a/NPSimulation/Core/MaterialManager.cc b/NPSimulation/Core/MaterialManager.cc
index 0df814159..666f6df87 100644
--- a/NPSimulation/Core/MaterialManager.cc
+++ b/NPSimulation/Core/MaterialManager.cc
@@ -243,6 +243,13 @@ G4Material* MaterialManager::GetMaterialFromLibrary(string Name){
m_Material[Name]=material;
return material;
}
+
+ else if(Name == "NaturalUranium"){
+ G4Material* material = new G4Material(Name, 19.1*g/cm3,1);
+ material->AddElement(GetElementFromLibrary("U"),1);
+ m_Material[Name]=material;
+ return material;
+ }
else if(Name == "CsI_Scintillator"){
G4Material* material = new G4Material(Name, 4.51*g/cm3,2);
diff --git a/NPSimulation/Core/PhysicsList.cc b/NPSimulation/Core/PhysicsList.cc
index dfc439efc..bf5576c3c 100644
--- a/NPSimulation/Core/PhysicsList.cc
+++ b/NPSimulation/Core/PhysicsList.cc
@@ -44,12 +44,15 @@
#include "G4eIonisation.hh"
#include "G4eBremsstrahlung.hh"
#include "G4eplusAnnihilation.hh"
-
+#include "G4NeutronHPFission.hh"
#include "G4MuIonisation.hh"
#include "G4MuMultipleScattering.hh"
#include "G4MuBremsstrahlung.hh"
#include "G4MuPairProduction.hh"
+//neutron
+#include "G4HadronFissionProcess.hh"
+
#include "G4hIonisation.hh"
#include "G4ionIonisation.hh"
#include "G4hMultipleScattering.hh"
@@ -79,7 +82,9 @@
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
PhysicsList::PhysicsList(){
// ie: no secondaries
- defaultCutValue = 1000 * pc;
+ //defaultCutValue = 1000 * pc;
+ defaultCutValue = 1*mm;
+
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
PhysicsList::~PhysicsList(){
@@ -91,22 +96,22 @@ void PhysicsList::ConstructParticle(){
// for all particles which you want to use.
// This ensures that objects of these particle types will be
// created in the program.
-
+
//Usefull to test geometry
G4Geantino::GeantinoDefinition();
-
+
//Usefull for Gamma
ConstructBosons();
-
+
//Usefull for betaDecay
ConstructLeptons();
-
+
//Needed by G4 (4.9.2) to run on mac os X ;-)
ConstructMesons();
-
+
//usefull for p and n
ConstructBaryons();
-
+
//Usefull of course :p
ConstructIons();
}
@@ -160,28 +165,28 @@ void PhysicsList::ConstructIons(){
G4Deuteron::DeuteronDefinition() ;
G4Triton::TritonDefinition() ;
G4Alpha::AlphaDefinition() ;
-
+
G4IonConstructor iConstructor ;
iConstructor.ConstructParticle() ;
-
+
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void PhysicsList::ConstructProcess(){
AddTransportation() ;
ConstructEM() ;
-
+
SetCuts() ;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void PhysicsList::ConstructEM(){
-
+
theParticleIterator->reset();
-
+
while ((*theParticleIterator)()) {
G4ParticleDefinition* particle = theParticleIterator->value() ;
G4ProcessManager* pmanager = particle->GetProcessManager() ;
G4String particleName = particle->GetParticleName() ;
-
+
if (particleName == "gamma") {
// gamma
//standard Geant4
@@ -196,7 +201,7 @@ void PhysicsList::ConstructEM(){
pmanager->AddProcess(new G4eIonisation , -1, 2, 2) ;
pmanager->AddProcess(new G4eBremsstrahlung , -1, -1, 3) ;
}
-
+
else if (particleName == "e+") {
//positron
pmanager->AddProcess(new G4eMultipleScattering , -1, 1, 1 );
@@ -205,54 +210,61 @@ void PhysicsList::ConstructEM(){
pmanager->AddProcess(new G4eBremsstrahlung , -1, -1, 3 );
pmanager->AddProcess(new G4eplusAnnihilation , 0, -1, 4 );
}
-
+
else if (particleName == "mu+" ||
- particleName == "mu-") {
+ particleName == "mu-") {
}
-
+
else if ( particleName == "GenericIon"
- || particleName == "He3"
- || particleName == "alpha"
- || particleName == "deuteron"
- || particleName == "triton") {
+ || particleName == "He3"
+ || particleName == "alpha"
+ || particleName == "deuteron"
+ || particleName == "triton") {
pmanager->AddProcess(new G4hMultipleScattering(), -1, 1, 1) ;
G4ionIonisation* iI = new G4ionIonisation ;
iI->ActivateStoppingData(true) ;
pmanager->AddProcess(iI , -1, 2, 2) ;
}
-
+
+ else if (particleName == "neutron"){
+ G4NeutronHPFission* fissionModel = new G4NeutronHPFission();
+ G4HadronFissionProcess* fissionProcess = new G4HadronFissionProcess();
+ fissionProcess->RegisterMe(fissionModel);
+ pmanager->AddDiscreteProcess(fissionProcess);
+ }
+
else if ((!particle->IsShortLived()) &&
- (particle->GetPDGCharge() != 0.0) &&
- (particleName != "chargedgeantino")) {
-
+ (particle->GetPDGCharge() != 0.0) &&
+ (particleName != "chargedgeantino")) {
+
G4hIonisation* hI = new G4hIonisation ;
pmanager->AddProcess(new G4hMultipleScattering , -1, 1, 1) ;
pmanager->AddProcess(hI , -1, 2, 2) ;
-
-
+
+
}
}
-
+
G4EmProcessOptions opt ;
opt.SetSubCutoff(true) ;
opt.SetMinEnergy(0.001*eV) ;
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);
+
+ //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......
@@ -263,20 +275,20 @@ void PhysicsList::SetCuts(){
// " G4VUserPhysicsList::SetCutsWithDefault" method sets
// the default cut value for all particle types
SetCutsWithDefault();
-
+
// for gamma-rays
/* 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+");*/
-
+ //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);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void PhysicsList::ConstructDecay(){
-
+
// Add Decay Process
G4Decay* theDecayProcess = new G4Decay() ;
theParticleIterator->reset() ;
@@ -291,7 +303,7 @@ void PhysicsList::ConstructDecay(){
}
}
//end Add Decay Process
-
+
}
void PhysicsList::MyOwnConstruction(){
diff --git a/NPSimulation/Core/Target.cc b/NPSimulation/Core/Target.cc
index df9eaccbc..803f58d24 100644
--- a/NPSimulation/Core/Target.cc
+++ b/NPSimulation/Core/Target.cc
@@ -501,6 +501,10 @@ G4double Target::SlowDownBeam(G4ParticleDefinition* Beam,
G4double ZInteraction,
G4double IncidentTheta){
+ if(Beam->GetParticleName()=="neutron"){
+ return IncidentEnergy;
+ }
+
G4double ThicknessBeforeInteraction =
abs(ZInteraction - 0.5*m_EffectiveThickness) / cos(m_TargetAngle);
--
GitLab