From ae4a521be158a871010841a61169cc4e92dfbc4f Mon Sep 17 00:00:00 2001
From: matta <matta@npt>
Date: Tue, 17 Jul 2012 14:22:18 +0000
Subject: [PATCH] Gamma and Particle Decay working
---
Inputs/EventGenerator/10He.reaction | 31 +-
.../include/EventGeneratorGammaDecay.hh | 77 ++++
.../include/EventGeneratorParticleDecay.hh | 89 ++++
.../include/PrimaryGeneratorAction.hh | 5 +-
NPSimulation/include/VEventGenerator.hh | 3 +-
NPSimulation/src/EventGeneratorGammaDecay.cc | 398 ++++++++++++++++
.../src/EventGeneratorParticleDecay.cc | 434 ++++++++++++++++++
NPSimulation/src/EventGeneratorTransfert.cc | 13 +-
NPSimulation/src/PrimaryGeneratorAction.cc | 247 +++++-----
9 files changed, 1171 insertions(+), 126 deletions(-)
create mode 100644 NPSimulation/include/EventGeneratorGammaDecay.hh
create mode 100644 NPSimulation/include/EventGeneratorParticleDecay.hh
create mode 100644 NPSimulation/src/EventGeneratorGammaDecay.cc
create mode 100644 NPSimulation/src/EventGeneratorParticleDecay.cc
diff --git a/Inputs/EventGenerator/10He.reaction b/Inputs/EventGenerator/10He.reaction
index fd541ea4a..95cc7b9f6 100644
--- a/Inputs/EventGenerator/10He.reaction
+++ b/Inputs/EventGenerator/10He.reaction
@@ -2,25 +2,42 @@
%%%%%%%%% Reaction file for 11Li(d,3He)10He reaction %%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Beam energy given in MeV ; Excitation in MeV
-TransfertToResonance
+Transfert
Beam= 11Li
Target= 2H
Light= 3He
Heavy= 10He
ExcitationEnergyLight= 0.0
- ExcitationEnergyHeavy= 1.0
+ ExcitationEnergyHeavy= 10.0
BeamEnergy= 550
BeamEnergySpread= 0
SigmaThetaX= 0.01
SigmaPhiY= 0.01
SigmaX= 0.01
SigmaY= 0.01
- ResonanceWidth= 0
- ResonanceDecayZ= 2
- ResonanceDecayA= 8
CrossSectionPath= 11Li(d,3He)10He.txt
ShootLight= 1
- ShootHeavy= 0
- ShootDecayProduct= 0
+ ShootHeavy= 1
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ParticleDecay 10He
+ Daughter= 9He n
+ ExcitationEnergy= 0.5 0
+ DifferentialCrossSection= flat.txt
+ shoot= 1 1
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+GammaDecay 9He
+ Cascade
+ BranchingRatio= 30
+ Energies= 0.1
+ DifferentialCrossSection= flat.txt
+ Cascade
+ BranchingRatio= 70
+ Energies= 0.2 0.3
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ParticleDecay 9He
+ Daughter= 8He n
+ DifferentialCrossSection= flat.txt
+ shoot= 1 1
\ No newline at end of file
diff --git a/NPSimulation/include/EventGeneratorGammaDecay.hh b/NPSimulation/include/EventGeneratorGammaDecay.hh
new file mode 100644
index 000000000..ae030178e
--- /dev/null
+++ b/NPSimulation/include/EventGeneratorGammaDecay.hh
@@ -0,0 +1,77 @@
+#ifndef EventGeneratorGammaDecay_H
+#define EventGeneratorGammaDecay_H
+/*****************************************************************************
+ * Copyright (C) 2009-2010 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: matta@ipno.in2p3.fr *
+ * *
+ * Creation Date : May 2012 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This event Generator is used to simulated a gamma decay of nuclei genera-*
+ * ted by previous event generator. Multiple cases are supported: *
+ * - Only one gamma is emmited, in this case a Cross section can be given *
+ * - A cascade decay, in this case the CS is ignore *
+ * - If more than one cascade are given, Branching Ratio could be given *
+ * *
+ * Using the "Source Token" the class can be used to simulate gamma ray sou-*
+ * rce and Doppler shifted gamma rays. *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ * *
+ *****************************************************************************/
+
+// STL
+#include <string>
+#include <iostream>
+using namespace std;
+
+// NPSimulation
+#include "VEventGenerator.hh"
+#include "Target.hh"
+#include "ParticleStack.hh"
+
+class EventGeneratorGammaDecay : public VEventGenerator
+{
+ public: // Constructor and destructor
+ EventGeneratorGammaDecay();
+ ~EventGeneratorGammaDecay();
+
+ public: // Inherit from VEventGenerator class
+ void ReadConfiguration(string,int);
+ void GenerateEvent(G4Event*);
+ void SetTarget(Target* Target) ;
+
+ private: // Target Parameter
+ Target* m_Target;
+
+
+ private: // The decaying nuclei
+ string m_NucleiName;
+ private: // the gamma rays property
+ vector<double> m_BranchingRatio;
+ vector<string> m_CrossSectionPath;
+ vector<double*> m_CrossSectionArray;
+ vector<double> m_CrossSectionThetaMin;
+ vector<double> m_CrossSectionThetaMax;
+ vector<double> m_CrossSectionSize;
+ vector<double> m_CascadeTotalEnergy;
+ vector< vector<double> > m_Energies;
+
+ private: // Pointer to the Particle stack for faster acces
+ ParticleStack* m_ParticleStack;
+ public: // Managing the different cascade
+ // Add a new cascade
+ void AddCascade(vector<double> Energies, double BranchingRatio=-1, string CrossSectionPath="_void_");
+ // Read all the added cscade en instentiate every thing that is needed
+ void PrepareCascade();
+};
+#endif
\ No newline at end of file
diff --git a/NPSimulation/include/EventGeneratorParticleDecay.hh b/NPSimulation/include/EventGeneratorParticleDecay.hh
new file mode 100644
index 000000000..2c344524a
--- /dev/null
+++ b/NPSimulation/include/EventGeneratorParticleDecay.hh
@@ -0,0 +1,89 @@
+#ifndef EventGeneratorParticleDecay_H
+#define EventGeneratorParticleDecay_H
+/*****************************************************************************
+ * Copyright (C) 2009-2010 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: matta@ipno.in2p3.fr *
+ * *
+ * Creation Date : May 2012 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This event Generator is used to simulated the particle decay of nuclei *
+ * generated by previous event generator. Multiple cases are supported: *
+ * - Only one particle is emmited, in this case a Cross section can be given*
+ * - If Multiple particle are emitted, a Phase Space generator is used *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ * *
+ *****************************************************************************/
+
+// STL
+#include <string>
+#include <iostream>
+using namespace std;
+
+// ROOT
+#include "TGenPhaseSpace.h"
+
+// G4 headers including CLHEP headers
+// for generating random numbers
+#include "Randomize.hh"
+
+// NPSimulation
+#include "VEventGenerator.hh"
+#include "Target.hh"
+#include "ParticleStack.hh"
+
+// Geant4
+#include "G4ParticleDefinition.hh"
+
+class EventGeneratorParticleDecay : public VEventGenerator
+{
+ public: // Constructor and destructor
+ EventGeneratorParticleDecay();
+ ~EventGeneratorParticleDecay();
+
+ public: // Inherit from VEventGenerator class
+ void ReadConfiguration(string,int);
+ void GenerateEvent(G4Event*);
+ void SetTarget(Target* Target) ;
+
+ private: // Target Parameter
+ Target* m_Target;
+
+ private: // The decaying nuclei
+ // General used
+ string m_MotherNucleiName;
+ vector<G4ParticleDefinition*> m_DaughterNuclei;
+ vector<string> m_DaughterName;
+ vector<double> m_ExcitationEnergy;
+ vector<bool> m_shoot;
+
+ // Used for cross section
+ string m_DifferentialCrossSection;
+ double* m_CrossSectionArray;
+ double m_CrossSectionThetaMin;
+ double m_CrossSectionThetaMax;
+ double m_CrossSectionSize;
+ RandGeneral* m_CrossSectionShoot;
+
+ // Used for the phase space
+ TGenPhaseSpace m_TPhaseSpace ;
+ double* m_Masses;
+
+ private: // Pointer to the Particle stack for faster acces
+ ParticleStack* m_ParticleStack;
+ public: // Managing the decay
+ // Set everything for the decay
+ void SetDecay(vector<string> DaughterName, vector<bool> shoot, vector<double> ExcitationEnergy, string CSPath);
+};
+#endif
\ No newline at end of file
diff --git a/NPSimulation/include/PrimaryGeneratorAction.hh b/NPSimulation/include/PrimaryGeneratorAction.hh
index dfc2cce7e..0313fa172 100644
--- a/NPSimulation/include/PrimaryGeneratorAction.hh
+++ b/NPSimulation/include/PrimaryGeneratorAction.hh
@@ -27,7 +27,6 @@
// G4 headers
#include "G4VUserPrimaryGeneratorAction.hh"
-#include "G4ParticleGun.hh"
#include "G4Event.hh"
// NPTool headers
@@ -46,15 +45,13 @@ public:
public:
void GeneratePrimaries(G4Event*);
- G4ParticleGun* GetParticleGun() {return m_particleGun;}
public:
void ReadEventGeneratorFile(string Path);
private:
- G4ParticleGun* m_particleGun;
DetectorConstruction* m_detector;
- VEventGenerator* m_EventGenerator;
+ vector<VEventGenerator*> m_EventGenerator;
};
#endif
diff --git a/NPSimulation/include/VEventGenerator.hh b/NPSimulation/include/VEventGenerator.hh
index b158b6ae8..f65a0b019 100644
--- a/NPSimulation/include/VEventGenerator.hh
+++ b/NPSimulation/include/VEventGenerator.hh
@@ -54,7 +54,8 @@ public:
public:
virtual void ReadConfiguration(string) {};
- virtual void GenerateEvent(G4Event*, G4ParticleGun*) {};
+ virtual void ReadConfiguration(string,int) {};
+ virtual void GenerateEvent(G4Event*) {};
virtual void InitializeRootOutput() {};
// Used in some case to generate event inside the target
diff --git a/NPSimulation/src/EventGeneratorGammaDecay.cc b/NPSimulation/src/EventGeneratorGammaDecay.cc
new file mode 100644
index 000000000..788c90fd5
--- /dev/null
+++ b/NPSimulation/src/EventGeneratorGammaDecay.cc
@@ -0,0 +1,398 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2010 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: matta@ipno.in2p3.fr *
+ * *
+ * Creation Date : May 2012 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This event Generator is used to simulated a gamma decay of nuclei genera-*
+ * ted by previous event generator. Multiple cases are supported: *
+ * - Only one gamma is emmited, in this case a Cross section can be given *
+ * - A cascade decay, in this case the CS is ignore *
+ * - If more than one cascade are given, Branching Ratio could be given *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ * *
+ *****************************************************************************/
+
+#include "EventGeneratorGammaDecay.hh"
+
+// NPS
+#include "Particle.hh"
+
+// G4 headers including CLHEP headers
+// for generating random numbers
+#include "Randomize.hh"
+
+// G4
+#include "G4ParticleTable.hh"
+
+// ROOT
+#include "TLorentzVector.h"
+#include "TVector3.h"
+EventGeneratorGammaDecay::EventGeneratorGammaDecay(){
+ m_ParticleStack = ParticleStack::getInstance();
+}
+
+EventGeneratorGammaDecay::~EventGeneratorGammaDecay(){
+
+}
+
+void EventGeneratorGammaDecay::ReadConfiguration(string Path, int Occurence){
+ ////////General Reading needs////////
+ string LineBuffer;
+ string DataBuffer;
+ istringstream LineStream;
+ int TokenOccurence = 0 ;
+ //////// Setting needs///////
+ unsigned int NumberOfCascade = 0;
+ bool ReadingStatusGammaDecay = false ;
+ bool CascadeStatus = false ;
+
+ bool check_E = false;
+ bool check_BranchingRatio = false;
+ bool check_CSPath = false ;
+ bool check_created = false;
+
+ //////////////////////////////////////////////////////////////////////////////////////////
+ ifstream InputFile;
+ InputFile.open(Path.c_str());
+
+ if (InputFile.is_open()) {} else {
+ return;
+ }
+
+ while (!InputFile.eof() && !check_created) {
+ //Pick-up next line
+ getline(InputFile, LineBuffer);
+
+ if (LineBuffer.compare(0, 10, "GammaDecay") == 0) {
+ TokenOccurence++;
+ if (TokenOccurence == Occurence) {
+ ReadingStatusGammaDecay = true ;
+ G4cout << "///////////////////////////////////////// " << G4endl;
+ // Get the nuclei name
+ LineStream.clear();
+ LineStream.str(LineBuffer);
+ LineStream >> DataBuffer;
+ DataBuffer.erase();
+ LineStream >> DataBuffer;
+ m_NucleiName = DataBuffer ;
+ G4cout << "Gamma Decay for " << m_NucleiName << G4endl;
+ }
+ }
+
+ ///////////////////////////////
+ /// Gamma Decay case
+ while(ReadingStatusGammaDecay){
+ InputFile >> DataBuffer;
+ //Search for comment Symbol %
+ if (DataBuffer.compare(0, 1, "%") == 0) { InputFile.ignore ( std::numeric_limits<std::streamsize>::max(), '\n' );}
+
+ else if (DataBuffer == "Cascade") {
+ CascadeStatus = true ;
+ NumberOfCascade++;
+ G4cout << " Cascade " << NumberOfCascade << G4endl;
+
+ LineStream.clear();
+ LineStream.str(LineBuffer);
+
+ // Instantiate new variable for the up coming cascade
+ check_E = false;
+ check_BranchingRatio = false;
+ check_CSPath = false ;
+
+ double BranchingRatio = -1;
+ vector<double> E ;
+ string CSPath;
+
+ while (CascadeStatus) {
+ getline(InputFile, LineBuffer);
+ LineStream.clear();
+ LineStream.str(LineBuffer);
+ LineStream >> DataBuffer ;
+
+ // G4cout << DataBuffer << G4endl;
+
+ //Search for comment Symbol %
+ if (DataBuffer.compare(0, 1, "%") == 0) { InputFile.ignore ( std::numeric_limits<std::streamsize>::max(), '\n' );}
+
+ else if(DataBuffer == "BranchingRatio=") {
+ check_BranchingRatio = true;
+ LineStream.clear();
+ LineStream.str(LineBuffer);
+ LineStream >> DataBuffer ;
+ LineStream >> DataBuffer ;
+ BranchingRatio = atof(DataBuffer.c_str());
+ G4cout << " Branching Ratio: " << atof(DataBuffer.c_str()) << G4endl;
+
+ }
+
+ else if(DataBuffer == "Energies=") {
+ check_E = true;
+ G4cout << " Energies: " ;
+ LineStream.clear();
+ LineStream.str(LineBuffer);
+ LineStream >> DataBuffer;
+ while(LineStream >> DataBuffer){
+ E.push_back(atof(DataBuffer.c_str()));
+ G4cout << atof(DataBuffer.c_str()) << " ";
+
+ }
+ G4cout << G4endl;
+ }
+
+ else if(DataBuffer == "DifferentialCrossSection="){
+ LineStream.clear();
+ LineStream.str(LineBuffer);
+ LineStream >> DataBuffer >> DataBuffer ;
+ CSPath = DataBuffer;
+ G4cout << " Cross Section Path: " << DataBuffer << G4endl;
+ check_CSPath = true;
+ }
+
+ // Cascade ended
+ if(check_E && check_BranchingRatio && E.size()>1){
+ AddCascade(E, BranchingRatio, "_void_");
+ CascadeStatus = false;
+ }
+
+ if(check_E && check_BranchingRatio && (E.size()<2 && check_CSPath)){ AddCascade(E, BranchingRatio, CSPath);
+ CascadeStatus = false;
+ }
+
+ }
+
+ }
+
+ //////////////////////////////////////////////////////
+ // If no Token and no comment, toggle out //
+ // else
+ // {ReadingStatusGammaDecay = false; G4cout << "WARNING : Wrong Token Sequence: Getting out " << G4endl ;}
+
+
+ if(InputFile.eof()) {ReadingStatusGammaDecay=false;check_created=true;}
+ }
+ }
+ G4cout << "///////////////////////////////////////// " << G4endl;
+ InputFile.close();
+ PrepareCascade();
+}
+
+
+void EventGeneratorGammaDecay::GenerateEvent(G4Event*){
+ // Choose a Cascade to follow
+ int ChoosenCascade = -1;
+ double RandomNumber = RandFlat::shoot();
+
+ for (unsigned int i = 1; i<m_BranchingRatio.size(); i++) {
+ if(RandomNumber > m_BranchingRatio[i-1] && RandomNumber< m_BranchingRatio[i])
+ ChoosenCascade=i;
+ }
+
+ if (ChoosenCascade==-1) ChoosenCascade=0;
+
+ // Look for the decaying nucleus
+ Particle decayingParticle = m_ParticleStack->SearchAndRemoveParticle(m_NucleiName);
+
+ if(decayingParticle.GetParticleDefinition()==NULL){
+ G4cout << "Gamma Decay Warning: The decaying particle " << m_NucleiName
+ << " was not found in the particle stack " << G4endl;
+ return ;
+ }
+ // Check for energies conservation (i.e: Cascade Energies lower than Excitation energie)
+ string ExcitationString = decayingParticle.GetParticleDefinition()->GetParticleName();
+ ExcitationString.erase(0,m_NucleiName.length()+1);
+ ExcitationString.erase(ExcitationString.length()-1,ExcitationString.length());
+ double ExcitationEnergy = atof(ExcitationString.c_str())*keV;
+
+ // Compute the final excitation energy of the decaying nuclei
+ double FinalExcitationEnergy = ExcitationEnergy-m_CascadeTotalEnergy[ChoosenCascade];
+ if(FinalExcitationEnergy<0){
+ G4cout << "Gamma Decay Warning: The cascade energy exceed the excitation energy of the decaying nuclei: "
+ << G4endl << " Excitation Energy : " << ExcitationEnergy
+ << G4endl << " Cascade Energy : " << m_CascadeTotalEnergy[ChoosenCascade] << G4endl;
+ FinalExcitationEnergy=0;
+ }
+
+ // Put back the decaying nucleus with its new excitation energy
+ G4ParticleDefinition* FinalParticleDefition
+ = G4ParticleTable::GetParticleTable()->GetIon(decayingParticle.GetParticleDefinition()->GetAtomicNumber(), decayingParticle.GetParticleDefinition()->GetAtomicMass(), FinalExcitationEnergy*MeV);
+
+ Particle FinalParticle = Particle( FinalParticleDefition,
+ decayingParticle.GetParticleKineticEnergy(),
+ decayingParticle.GetParticleMomentumDirection(),
+ decayingParticle.GetParticlePosition(),
+ decayingParticle.GetShootStatus());
+
+ m_ParticleStack->AddParticleToStack(FinalParticle);
+
+ // Instantiate and add the gamma to the particle stack
+ for (unsigned int i = 0; i < m_Energies[ChoosenCascade].size(); i++) {
+ G4ParticleDefinition* gammaDefinition = G4ParticleTable::GetParticleTable()->FindParticle("gamma");
+ G4ThreeVector gammaDirection;
+ double theta=0;
+ double phi=0;
+ // If more than one gamma shoot no cross section to follow
+ if(m_Energies[ChoosenCascade].size()>1){
+
+ // Shoot flat in cos(theta) and Phi to have isotropic emission
+ double cos_theta = RandFlat::shoot();
+ theta = acos(cos_theta);
+ phi = RandFlat::shoot()*2.*pi;
+
+ gammaDirection= G4ThreeVector( cos(phi)*sin(theta),
+ sin(phi)*sin(theta),
+ cos(theta));
+ }
+
+ // Only one gamma to shoot, use the given cross section
+ else{
+ RandGeneral CrossSectionShoot(m_CrossSectionArray[ChoosenCascade], m_CrossSectionSize[ChoosenCascade]);
+ theta = ( m_CrossSectionThetaMin[ChoosenCascade]
+ + CrossSectionShoot.shoot()
+ * ( m_CrossSectionThetaMax[ChoosenCascade] - m_CrossSectionThetaMin[ChoosenCascade]) ) * deg;
+ phi = RandFlat::shoot() * 2. *pi;
+ gammaDirection= G4ThreeVector( cos(phi)*sin(theta),
+ sin(phi)*sin(theta),
+ cos(theta));
+ }
+
+ // Doppler shifted gamma emission
+ decayingParticle.GetParticleMomentumDirection();
+ double gammaEnergy = m_Energies[ChoosenCascade][i];
+ TLorentzVector GammaLV( gammaEnergy*cos(phi)*sin(theta),
+ gammaEnergy*sin(phi)*sin(theta),
+ gammaEnergy*cos(theta),
+ gammaEnergy);
+
+ double NucleiEnergy= decayingParticle.GetParticleKineticEnergy()+FinalParticleDefition->GetPDGMass();
+ double NucleiMomentum= sqrt(NucleiEnergy*NucleiEnergy-FinalParticleDefition->GetPDGMass()*FinalParticleDefition->GetPDGMass());
+ TLorentzVector NuvleiLV( NucleiMomentum*decayingParticle.GetParticleMomentumDirection().x(),
+ NucleiMomentum*decayingParticle.GetParticleMomentumDirection().y(),
+ NucleiMomentum*decayingParticle.GetParticleMomentumDirection().z(),
+ NucleiEnergy);
+
+ GammaLV.Boost(NuvleiLV.BoostVector());
+ gammaDirection= G4ThreeVector( GammaLV.Px(),
+ GammaLV.Py(),
+ GammaLV.Pz());
+ // Add the gamma to the stack
+ Particle gammaParticle(gammaDefinition,GammaLV.E(),gammaDirection, decayingParticle.GetParticlePosition());
+ m_ParticleStack->AddParticleToStack(gammaParticle);
+ }
+}
+
+
+void EventGeneratorGammaDecay::SetTarget(Target* Target){
+ m_Target = Target;
+}
+
+void EventGeneratorGammaDecay::AddCascade(vector<double> Energies, double BranchingRatio, string CrossSectionPath){
+ m_BranchingRatio.push_back(BranchingRatio);
+ m_CrossSectionPath.push_back(CrossSectionPath);
+ m_Energies.push_back(Energies);
+}
+
+void EventGeneratorGammaDecay::PrepareCascade(){
+ // Change the given branching ratio so total is one (always have a decay during the event)
+ double TotalRatio=0;
+ for (unsigned int i = 0; i < m_BranchingRatio.size(); i++) {
+ TotalRatio+=m_BranchingRatio[i]/100.;
+ }
+
+ // Check that the total ratio is not over 100% (below is allowed)
+ if(TotalRatio>1) {
+ G4cout << "Gamma Decay Error: Sum of branching ratio is over 100%" << endl;
+ exit(1);
+ }
+
+ for (unsigned int i = 0; i < m_BranchingRatio.size(); i++) {
+ m_BranchingRatio[i]=(m_BranchingRatio[i]/TotalRatio)/100.;
+ }
+
+ // Shift the Branching ratio for faster shooting during event generation
+ for (unsigned int i = 1; i < m_BranchingRatio.size(); i++) {
+ m_BranchingRatio[i]=m_BranchingRatio[i-1]+m_BranchingRatio[i];
+ }
+
+ // Compute the total energy of the cascade
+ double TotalEnergy=0;
+ for (unsigned int i = 0 ; i < m_Energies.size(); i++) {
+ TotalEnergy=0;
+ for (unsigned int j = 0; j < m_Energies[i].size(); j++) {
+ TotalEnergy+=m_Energies[i][j];
+ }
+ m_CascadeTotalEnergy.push_back(TotalEnergy);
+ }
+
+ // Transform the particle name to G4 standard: i.e: 10He -> He10
+ m_NucleiName = m_ParticleStack->ChangeNameToG4Standard(m_NucleiName);
+ if (m_NucleiName=="proton" ||
+ m_NucleiName=="deuteron" ||
+ m_NucleiName=="triton" ||
+ m_NucleiName=="3He" ||
+ m_NucleiName=="alpha" ){
+ G4cout << "Gamma Decay Error: Gamma Decay not allowed for light particles" << endl;
+ exit(1);
+ }
+
+ /// Load the differential cross section
+ for (unsigned int i = 0; i < m_BranchingRatio.size(); i++) {
+ unsigned int CrossSectionSize = 0;
+ double* CrossSection = new double[CrossSectionSize] ;
+ double thetamin = 0;
+ double thetamax = 0;
+
+ if(m_CrossSectionPath[i]!="_void_"){
+ string GlobalPath = getenv("NPTOOL");
+ string StandardPath = GlobalPath + "/Inputs/CrossSection/" + m_CrossSectionPath[i];
+ ifstream CSFile;
+ CSFile.open( StandardPath.c_str() );
+
+ if(CSFile.is_open()) cout << "Reading Cross Section File " << m_CrossSectionPath[i] << endl;
+
+ // In case the file is not found in the standard path, the programm try to interpret the file name as an absolute or relative file path.
+ else{
+ CSFile.open( m_CrossSectionPath[i].c_str() );
+ if(CSFile.is_open()) { cout << "Reading Cross Section File " << m_CrossSectionPath[i] << endl;}
+
+ else {cout << "Cross Section File " << m_CrossSectionPath[i] << " not found" << endl;return;}
+ }
+
+ double CSBuffer,AngleBuffer;
+ vector<double> CrossSectionBuffer;
+ thetamin = 200;
+ thetamax = -10;
+ while(!CSFile.eof()) {
+ CSFile >> AngleBuffer;
+ CSFile >> CSBuffer;
+ double CSFinal = CSBuffer*sin(AngleBuffer*deg);
+ CrossSectionBuffer.push_back(CSFinal);
+ if (AngleBuffer < thetamin) thetamin = AngleBuffer;
+ if (AngleBuffer > thetamax) thetamax = AngleBuffer;
+ }
+
+ CSFile.close();
+ CrossSectionSize = CrossSectionBuffer.size();
+ delete CrossSection;
+ CrossSection = new double[CrossSectionSize] ;
+ for(unsigned int i = 0 ; i <CrossSectionSize ; i++ ) CrossSection[i] = CrossSectionBuffer[i];
+ }
+ m_CrossSectionArray.push_back(CrossSection);
+ m_CrossSectionSize.push_back(CrossSectionSize);
+ m_CrossSectionThetaMax.push_back(thetamax);
+ m_CrossSectionThetaMin.push_back(thetamin);
+ }
+}
+
diff --git a/NPSimulation/src/EventGeneratorParticleDecay.cc b/NPSimulation/src/EventGeneratorParticleDecay.cc
new file mode 100644
index 000000000..813d00211
--- /dev/null
+++ b/NPSimulation/src/EventGeneratorParticleDecay.cc
@@ -0,0 +1,434 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2010 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: matta@ipno.in2p3.fr *
+ * *
+ * Creation Date : May 2012 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This event Generator is used to simulated the particle decay of nuclei *
+ * generated by previous event generator. Multiple cases are supported: *
+ * - Only one particle is emmited, in this case a Cross section can be given*
+ * - If Multiple particle are emitted, a Phase Space generator is used *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ * *
+ *****************************************************************************/
+
+#include "EventGeneratorParticleDecay.hh"
+
+// NPS
+#include "Particle.hh"
+
+// G4
+#include "G4ParticleTable.hh"
+
+// NPL
+#include "NPNucleus.h"
+using namespace NPL;
+
+// ROOT
+#include "TLorentzVector.h"
+#include "TVector3.h"
+EventGeneratorParticleDecay::EventGeneratorParticleDecay(){
+ m_ParticleStack = ParticleStack::getInstance();
+}
+
+EventGeneratorParticleDecay::~EventGeneratorParticleDecay(){
+
+}
+
+void EventGeneratorParticleDecay::ReadConfiguration(string Path,int Occurence){
+ ////////General Reading needs////////
+ string LineBuffer;
+ string DataBuffer;
+ istringstream LineStream;
+ int TokkenOccurence = 0;
+ //////// Setting needs///////
+ bool ReadingStatusParticleDecay = false ;
+
+ bool check_Daughter = false ;
+ bool check_CrossSection = false ;
+ bool check_ExcitationEnergy = false ;
+ bool check_shoot = false ;
+ bool check_created = false ;
+
+ // Instantiate new variable for the up coming Particle
+ vector<string> DaughterName;
+ vector<double> ExcitationEnergy;
+ vector<bool> shoot;
+ string CSPath = "TGenPhaseSpace";
+
+ //////////////////////////////////////////////////////////////////////////////////////////
+ ifstream InputFile;
+ InputFile.open(Path.c_str());
+
+ if (InputFile.is_open()) {}
+
+ else {
+ return;
+ }
+
+ while (!InputFile.eof()&& !check_created) {
+ //Pick-up next line
+ getline(InputFile, LineBuffer);
+
+ if (LineBuffer.compare(0, 13, "ParticleDecay") == 0) {
+ TokkenOccurence++;
+ if(TokkenOccurence==Occurence){
+ ReadingStatusParticleDecay = true ;
+ G4cout << "///////////////////////////////////////// " << G4endl;
+ // Get the nuclei name
+ LineStream.clear();
+ LineStream.str(LineBuffer);
+ LineStream >> DataBuffer;
+ DataBuffer.erase();
+ LineStream >> DataBuffer;
+ m_MotherNucleiName = DataBuffer ;
+ G4cout << "Particle Decay for " << m_MotherNucleiName << G4endl;
+ }
+ }
+
+ ///////////////////////////////
+ /// Gamma Decay case
+ while(ReadingStatusParticleDecay){
+
+ InputFile >> DataBuffer;
+ //Search for comment Symbol %
+ if (DataBuffer.compare(0, 1, "%") == 0) {
+ InputFile.ignore ( std::numeric_limits<std::streamsize>::max(), '\n' );
+ }
+
+ else if (DataBuffer == "Daughter=") {
+ check_Daughter = true ;
+ LineStream.clear();
+ LineStream.str(LineBuffer);
+
+ getline(InputFile, LineBuffer);
+ LineStream.clear();
+ LineStream.str(LineBuffer);
+ G4cout << " Daughter: " ;
+ while(LineStream >> DataBuffer){
+ DaughterName.push_back(DataBuffer);
+ G4cout << DataBuffer << " " ;
+ }
+
+ G4cout << G4endl;
+
+ }
+
+ else if(DataBuffer == "ExcitationEnergy=") {
+ check_ExcitationEnergy = true;
+ LineStream.clear();
+ LineStream.str(LineBuffer);
+
+ getline(InputFile, LineBuffer);
+ LineStream.clear();
+ LineStream.str(LineBuffer);
+ G4cout << " Excitation Energy: " ;
+ while(LineStream >> DataBuffer){
+ ExcitationEnergy.push_back( atof(DataBuffer.c_str()) );
+ G4cout << DataBuffer << " " ;
+ }
+
+ G4cout << G4endl;
+ }
+
+ else if(DataBuffer == "DifferentialCrossSection=") {
+ check_CrossSection = true;
+ InputFile >> CSPath ;
+ G4cout << " Cross Section: " << CSPath << G4endl;
+ }
+
+ else if(DataBuffer == "shoot=") {
+ check_shoot = true;
+ LineStream.clear();
+ LineStream.str(LineBuffer);
+
+ getline(InputFile, LineBuffer);
+ LineStream.clear();
+ LineStream.str(LineBuffer);
+ G4cout << " Shoot Particle: " ;
+ while(LineStream >> DataBuffer){
+ shoot.push_back( atof(DataBuffer.c_str()) );
+ G4cout << DataBuffer << " " ;
+ }
+
+ G4cout << G4endl;
+ }
+
+
+
+ //////////////////////////////////////////////////////
+ // If no Token and no comment, toggle out //
+ else
+ {ReadingStatusParticleDecay = false; G4cout << "ERROR : Wrong Token Sequence: Getting out " << G4endl ;
+ exit(1);
+ }
+
+ // Decay ended
+ if(check_Daughter && check_shoot){
+ SetDecay(DaughterName,shoot,ExcitationEnergy,CSPath);
+ ReadingStatusParticleDecay = false;
+ check_created=true;
+ }
+ }
+ }
+
+ G4cout << "///////////////////////////////////////// " << G4endl;
+ InputFile.close();
+}
+
+
+void EventGeneratorParticleDecay::GenerateEvent(G4Event* anEvent){
+ // Look for the decaying nucleus
+ Particle decayingParticle = m_ParticleStack->SearchAndRemoveParticle(m_MotherNucleiName);
+ if(decayingParticle.GetParticleDefinition()==NULL){
+ G4cout << "Gamma Decay Warning: The decaying particle " << m_MotherNucleiName
+ << " was not found in the particle stack " << G4endl;
+ return ;
+ }
+
+ G4ParticleDefinition* decayingParticleDefinition = decayingParticle.GetParticleDefinition();
+ // Build the decaying particle four momenta vector:
+
+ double NucleiEnergy= decayingParticle.GetParticleKineticEnergy()+decayingParticleDefinition->GetPDGMass();
+
+ double NucleiMomentum=sqrt(NucleiEnergy*NucleiEnergy -
+ decayingParticleDefinition->GetPDGMass()*decayingParticleDefinition->GetPDGMass());
+
+ G4ThreeVector Momentum = decayingParticle.GetParticleMomentumDirection().unit();
+
+ // Case of one particle decaying with a user given cross section
+ if(m_DifferentialCrossSection!="TGenPhaseSpace"){
+ TLorentzVector NucleiLV( NucleiMomentum*Momentum.x(),
+ NucleiMomentum*Momentum.y(),
+ NucleiMomentum*Momentum.z(),
+ NucleiEnergy);
+ // Shoot the angle in Center of Mass (CM) frame
+ G4double ThetaCM = (m_CrossSectionThetaMin + m_CrossSectionShoot->shoot() * (m_CrossSectionThetaMax - m_CrossSectionThetaMin)) * deg;
+ G4double phi = RandFlat::shoot()*2.*pi;
+
+ // Build daughter particule CM LV
+ // Pre compute variable for the decay
+ double M = decayingParticleDefinition->GetPDGMass();
+ double m1 = m_DaughterNuclei[0]->GetPDGMass();
+ double m2 = m_DaughterNuclei[1]->GetPDGMass();
+
+ if(M<(m1+m2))
+ cout << "Warning: Particle Decay forbiden by kinematic, no particle emitted "<<endl;
+
+ else {
+ double Energy = ( 1./(2.*M) )*( M*M + m1*m1 - m2*m2);
+ double Momentum = sqrt(Energy*Energy - m1*m1);
+
+ TVector3 FirstDaughterMomentum = Momentum * TVector3( sin(ThetaCM) * cos(phi),
+ sin(ThetaCM) * sin(phi),
+ cos(ThetaCM));
+
+ TLorentzVector FirstDaughterLV(FirstDaughterMomentum,Energy);
+
+ FirstDaughterLV.Boost( NucleiLV.BoostVector() );
+ TLorentzVector SecondDaughterLV = NucleiLV - FirstDaughterLV;
+
+ G4ThreeVector DaughterDirection = G4ThreeVector( FirstDaughterLV.X() ,
+ FirstDaughterLV.Y() ,
+ FirstDaughterLV.Z() );
+
+ Particle FirstDaughterParticle( m_DaughterNuclei[0],
+ FirstDaughterLV.E()-m_DaughterNuclei[0]->GetPDGMass(),
+ DaughterDirection.unit(),
+ decayingParticle.GetParticlePosition(),
+ m_shoot[0]);
+
+ DaughterDirection = G4ThreeVector( SecondDaughterLV.X() ,
+ SecondDaughterLV.Y() ,
+ SecondDaughterLV.Z() );
+
+ Particle SecondDaughterParticle( m_DaughterNuclei[1],
+ SecondDaughterLV.E()-m_DaughterNuclei[1]->GetPDGMass(),
+ DaughterDirection.unit(),
+ decayingParticle.GetParticlePosition(),
+ m_shoot[1]);
+
+ ParticleStack::getInstance()->AddParticleToStack(FirstDaughterParticle);
+ ParticleStack::getInstance()->AddParticleToStack(SecondDaughterParticle);
+ }
+ }
+
+ // Case of a TGenPhaseSpace
+ else{
+ TLorentzVector NucleiLV( NucleiMomentum*Momentum.x()/GeV,
+ NucleiMomentum*Momentum.y()/GeV,
+ NucleiMomentum*Momentum.z()/GeV,
+ NucleiEnergy/GeV);
+
+ if( !m_TPhaseSpace.SetDecay(NucleiLV, m_DaughterNuclei.size(), m_Masses) )
+ cout << "Warning: Phase Space Decay forbiden by kinematic, or more than 18 particles, no particle emitted "<<endl;
+
+ else{
+ m_TPhaseSpace.Generate();
+
+ TLorentzVector* daughterLV ;
+ double KineticEnergy;
+
+ for (unsigned int i = 0 ; i < m_DaughterNuclei.size(); i++) {
+
+ daughterLV = m_TPhaseSpace.GetDecay(i);
+ G4ThreeVector daughterDirection = G4ThreeVector( daughterLV->X() ,
+ daughterLV->Y() ,
+ daughterLV->Z() );
+
+ KineticEnergy = daughterLV->E()-m_Masses[i] ;
+
+ Particle daughterParticle( m_DaughterNuclei[i],
+ KineticEnergy*GeV,
+ daughterDirection.unit(),
+ decayingParticle.GetParticlePosition(),
+ m_shoot[i]);
+ ParticleStack::getInstance()->AddParticleToStack(daughterParticle);
+ }
+ }
+ }
+}
+
+
+void EventGeneratorParticleDecay::SetTarget(Target* Target){
+ m_Target = Target;
+}
+
+void EventGeneratorParticleDecay::SetDecay(vector<string> DaughterName, vector<bool> shoot, vector<double> ExcitationEnergy, string CSPath){
+
+ // Check the validity of the given data:
+ if (DaughterName.size() != shoot.size() || (DaughterName.size() != ExcitationEnergy.size() && ExcitationEnergy.size()!=0) ) {
+ G4cout << "ERROR : Missmatching information: Getting out " << G4endl ;
+ exit(1);
+ }
+
+ if ( DaughterName.size() != 2 && CSPath!="TGenPhaseSpace" ) {
+ G4cout << "ERROR 2: Missmatching information: Getting out " << G4endl ;
+ exit(1);
+ }
+
+ m_shoot = shoot ;
+ m_ExcitationEnergy= ExcitationEnergy ;
+
+ // If the Excitation Energy Token was omitted, then it is set to zero
+ if(m_ExcitationEnergy.size()==0)
+ for (unsigned int i = 0 ; i < DaughterName.size(); i++) {
+ m_ExcitationEnergy.push_back(0);
+ }
+
+ // used check for mass and charge conservation
+ Nucleus* myNucleus = new Nucleus(m_MotherNucleiName);
+ int InitialCharge = myNucleus->GetZ() ; int FinalCharge = 0 ;
+ int InitialMass = myNucleus->GetA() ; int FinalMass = 0 ;
+ delete myNucleus;
+ for (unsigned int i = 0 ; i< DaughterName.size(); i++) {
+ if(DaughterName[i] == "p"){
+ m_DaughterNuclei.push_back(G4ParticleTable::GetParticleTable()->FindParticle("proton"));
+ FinalMass++;FinalCharge++;
+ }
+
+ else if (DaughterName[i] == "n"){
+ m_DaughterNuclei.push_back(G4ParticleTable::GetParticleTable()->FindParticle("neutron"));
+ FinalMass++;
+ }
+
+ else{
+ Nucleus* myNucleus = new Nucleus(DaughterName[i]);
+ m_DaughterNuclei.push_back(G4ParticleTable::GetParticleTable()->GetIon(myNucleus->GetZ(),
+ myNucleus->GetA(),
+ m_ExcitationEnergy[i]*MeV));
+ FinalMass+=myNucleus->GetA();
+ FinalCharge+=myNucleus->GetZ();
+ delete myNucleus;
+ }
+ }
+
+ // Check mass and charge conservation
+ if (InitialMass!=FinalMass || InitialCharge!=FinalCharge) {
+ G4cout << "ERROR: Mass and charge are not conserved." << G4endl;
+ exit(1);
+ }
+
+ m_DaughterName = DaughterName;
+
+ m_DifferentialCrossSection = CSPath;
+ if(CSPath!="TGenPhaseSpace") {
+ unsigned int CrossSectionSize = 0;
+ double* m_CrossSection = new double[CrossSectionSize] ;
+ double m_CrossSectionThetaMin = 0;
+ double m_CrossSectionThetaMax = 0;
+
+ string GlobalPath = getenv("NPTOOL");
+ string StandardPath = GlobalPath + "/Inputs/CrossSection/" + m_DifferentialCrossSection;
+ ifstream CSFile;
+ CSFile.open( StandardPath.c_str() );
+
+ if(CSFile.is_open()) cout << "Reading Cross Section File " << m_DifferentialCrossSection << endl;
+
+ // In case the file is not found in the standard path, the programm try to interpret the file name as an absolute or relative file path.
+ else{
+ CSFile.open( m_DifferentialCrossSection.c_str() );
+ if(CSFile.is_open()) {
+ cout << "Reading Cross Section File " << m_DifferentialCrossSection << endl;
+ }
+
+ else {
+ cout << "ERROR : Cross Section File " << m_DifferentialCrossSection << " not found" << endl;
+ exit(1);
+ }
+ }
+
+ double CSBuffer,AngleBuffer;
+ vector<double> CrossSectionBuffer;
+ m_CrossSectionThetaMin = 200;
+ m_CrossSectionThetaMax = -10;
+ while(!CSFile.eof()) {
+ CSFile >> AngleBuffer;
+ CSFile >> CSBuffer;
+ double CSFinal = CSBuffer*sin(AngleBuffer*deg);
+ CrossSectionBuffer.push_back(CSFinal);
+ if (AngleBuffer < m_CrossSectionThetaMin) m_CrossSectionThetaMin = AngleBuffer;
+ if (AngleBuffer > m_CrossSectionThetaMax) m_CrossSectionThetaMax = AngleBuffer;
+ }
+
+ CSFile.close();
+
+ m_CrossSectionSize = CrossSectionBuffer.size();
+ m_CrossSectionArray = new double[CrossSectionBuffer.size()] ;
+
+ for(unsigned int i = 0 ; i < m_CrossSectionSize ; i++ ) {
+ m_CrossSectionArray[i] = CrossSectionBuffer[i];
+ }
+
+ m_CrossSectionShoot = new RandGeneral(m_CrossSectionArray, m_CrossSectionSize);
+
+ }
+
+
+ else{
+ // Set up the array of masses
+ m_Masses = new double[m_DaughterNuclei.size()];
+
+ // Mass of the daugther nuclei are set once
+ for (unsigned int i = 0 ; i < m_DaughterNuclei.size(); i++) {
+ m_Masses[i] = m_DaughterNuclei[i]->GetPDGMass()/GeV;
+ }
+
+ }
+
+ // Change the name of the decaying nucleus to G4 standard
+ m_MotherNucleiName = m_ParticleStack->ChangeNameToG4Standard(m_MotherNucleiName);
+
+}
\ No newline at end of file
diff --git a/NPSimulation/src/EventGeneratorTransfert.cc b/NPSimulation/src/EventGeneratorTransfert.cc
index feb67d4bf..435b5f687 100644
--- a/NPSimulation/src/EventGeneratorTransfert.cc
+++ b/NPSimulation/src/EventGeneratorTransfert.cc
@@ -185,7 +185,11 @@ void EventGeneratorTransfert::ReadConfiguration(string Path)
- if (LineBuffer.compare(0, 9, "Transfert") == 0) { ReadingStatus = true ;}
+ if (LineBuffer.compare(0, 9, "Transfert") == 0) {
+ ReadingStatus = true ;
+
+ G4cout << "////////// Transfert ///////////" << G4endl;
+ }
while(ReadingStatus){
@@ -323,10 +327,11 @@ void EventGeneratorTransfert::ReadConfiguration(string Path)
// If all Token found toggle out
if(check_Beam && check_Target && check_Light && check_Heavy && check_ExcitationEnergyLight && check_ExcitationEnergyHeavy
&& check_BeamEnergy && check_BeamEnergySpread && check_FWHMX && check_FWHMY && check_EmmitanceTheta
- && check_EmmitancePhi && check_CrossSectionPath && check_ShootLight && check_ShootHeavy)
- ReadingStatus = false ;
-
+ && check_EmmitancePhi && check_CrossSectionPath && check_ShootLight && check_ShootHeavy){
+ ReadingStatus = false ;
+ G4cout << "////////////////////////////////" << G4endl;
}
+ }
}
diff --git a/NPSimulation/src/PrimaryGeneratorAction.cc b/NPSimulation/src/PrimaryGeneratorAction.cc
index 50f0ae0df..20d7b1902 100644
--- a/NPSimulation/src/PrimaryGeneratorAction.cc
+++ b/NPSimulation/src/PrimaryGeneratorAction.cc
@@ -20,151 +20,178 @@
* Comment: *
* *
*****************************************************************************/
+// STL
+#include "math.h"
+#include <cstdlib>
-#include "PrimaryGeneratorAction.hh"
-#include "PhysicsList.hh"
+// G4
#include "G4ParticleTable.hh"
#include "G4ParticleDefinition.hh"
+
+// NPL
+#include "PrimaryGeneratorAction.hh"
+#include "PhysicsList.hh"
#include "Randomize.hh"
-#include "math.h"
+#include "ParticleStack.hh"
+
+// CLHEP
#include "CLHEP/Random/RandGauss.h"
#include "CLHEP/Random/RandGeneral.h"
// Event Generator Class
#include "EventGeneratorTransfert.hh"
-#include "EventGeneratorTransfertToResonance.hh"
#include "EventGeneratorIsotropic.hh"
#include "EventGeneratorBeam.hh"
#include "EventGeneratorPhaseSpace.hh"
+#include "EventGeneratorGammaDecay.hh"
+#include "EventGeneratorParticleDecay.hh"
-#include <cstdlib>
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
PrimaryGeneratorAction::~PrimaryGeneratorAction()
{
- delete m_particleGun;
- delete m_EventGenerator;
+ for (unsigned int i = 0 ; i < m_EventGenerator.size(); i++) {
+ delete m_EventGenerator[i];
+ }
+ m_EventGenerator.clear();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
PrimaryGeneratorAction::PrimaryGeneratorAction(DetectorConstruction* det)
: m_detector(det)
{
- m_EventGenerator = NULL;
-
- m_particleGun = new G4ParticleGun(1);
- G4ParticleDefinition* particle
- = G4ParticleTable::GetParticleTable()->FindParticle("proton");
- m_particleGun->SetParticleDefinition(particle);
- m_particleGun->SetParticleEnergy(10*MeV);
- m_particleGun->SetParticlePosition(G4ThreeVector(0., 0., 0.));
- m_particleGun->SetParticleMomentumDirection(G4ThreeVector(0., 0., 1.));
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
{
- if (m_EventGenerator != NULL)m_EventGenerator->GenerateEvent(anEvent, m_particleGun) ;
- else {
- G4cout << "No Event Generator!" << G4endl ;
- return ;
- }
+ for (unsigned int i = 0 ; i < m_EventGenerator.size(); i++) {
+ m_EventGenerator[i]->GenerateEvent(anEvent);
+ }
+ ParticleStack::getInstance()->ShootAllParticle(anEvent);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void PrimaryGeneratorAction::ReadEventGeneratorFile(string Path)
{
- bool check_TransfertToResonance = false;
- bool check_PhaseSpace = false;
- bool check_Isotropic = false;
- bool check_Transfert = false;
- bool check_Beam = false;
-
- string LineBuffer;
- ifstream EventGeneratorFile;
- EventGeneratorFile.open(Path.c_str());
-
- if (EventGeneratorFile.is_open()) { // should always be true
- cout << "Event Generator file " << Path << " loading " << endl ;
- }
- else {
- cout << "Error, Event Generator file " << Path << " found" << endl;
- }
-
-
- while (!EventGeneratorFile.eof()) {
- //Pick-up next line
- getline(EventGeneratorFile, LineBuffer);
-
- //Search for comment Symbol %
- if (LineBuffer.compare(0, 1, "%") == 0) { /*do nothing*/
- ;
- }
-
-
- //Search for Isotropic source
- else if (LineBuffer.compare(0, 9, "Isotropic") == 0 && !check_Isotropic) {
- check_Isotropic = true;
- VEventGenerator* myEventGenerator = new EventGeneratorIsotropic();
- EventGeneratorFile.close();
- myEventGenerator->ReadConfiguration(Path);
- EventGeneratorFile.open(Path.c_str());
- myEventGenerator->InitializeRootOutput();
- m_EventGenerator = myEventGenerator;
- }
-
- //Search for Beam
- else if (LineBuffer.compare(0, 4, "Beam") == 0 && !check_Beam) {
- check_Beam = true;
- VEventGenerator* myEventGenerator = new EventGeneratorBeam();
- EventGeneratorFile.close();
- myEventGenerator->ReadConfiguration(Path);
- EventGeneratorFile.open(Path.c_str());
- myEventGenerator->InitializeRootOutput();
- myEventGenerator->SetTarget(m_detector->GetTarget());
- m_EventGenerator = myEventGenerator;
- }
-
-
- //Search for Transfert
- else if (LineBuffer.compare(0, 9, "Transfert") == 0 && !check_Transfert && LineBuffer.compare(0, 11, "TransfertTo") != 0) {
- check_Transfert = true;
- VEventGenerator* myEventGenerator = new EventGeneratorTransfert();
- EventGeneratorFile.close();
- myEventGenerator->ReadConfiguration(Path);
- EventGeneratorFile.open(Path.c_str());
- myEventGenerator->InitializeRootOutput();
- myEventGenerator->SetTarget(m_detector->GetTarget());
- m_EventGenerator = myEventGenerator;
- }
-
- //Search for Transfert To Resonance
- else if (LineBuffer.compare(0, 21, "TransfertToResonance") == 0 && !check_TransfertToResonance) {
- check_TransfertToResonance = true;
- VEventGenerator* myEventGenerator = new EventGeneratorTransfertToResonance();
- EventGeneratorFile.close();
- myEventGenerator->ReadConfiguration(Path);
- EventGeneratorFile.open(Path.c_str());
- myEventGenerator->InitializeRootOutput();
- myEventGenerator->SetTarget(m_detector->GetTarget());
- m_EventGenerator = myEventGenerator;
- }
-
- //Search for Transfert To Resonance
- else if (LineBuffer.compare(0, 10, "PhaseSpace") == 0 && !check_PhaseSpace) {
- check_PhaseSpace = true;
- VEventGenerator* myEventGenerator = new EventGeneratorPhaseSpace();
- EventGeneratorFile.close();
- myEventGenerator->ReadConfiguration(Path);
- EventGeneratorFile.open(Path.c_str());
- myEventGenerator->InitializeRootOutput();
- myEventGenerator->SetTarget(m_detector->GetTarget());
- m_EventGenerator = myEventGenerator;
- }
- }
- EventGeneratorFile.close();
+ bool check_TransfertToResonance = false;
+ bool check_PhaseSpace = false;
+ bool check_Isotropic = false;
+ bool check_Transfert = false;
+ bool check_Beam = false;
+
+ // You can have more than one of those
+ int alreadyiInstantiate_GammaDecay = 0;
+ int seenToken_GammaDecay = 0;
+ int alreadyiInstantiate_ParticleDecay = 0;
+ int seenToken_ParticleDecay = 0;
+
+ string LineBuffer;
+ ifstream EventGeneratorFile;
+ EventGeneratorFile.open(Path.c_str());
+
+ if (EventGeneratorFile.is_open()) { // should always be true
+ cout << "Event Generator file " << Path << " loading " << endl ;
+ }
+ else {
+ cout << "Error, Event Generator file " << Path << " found" << endl;
+ }
+
+
+ while (!EventGeneratorFile.eof()) {
+ //Pick-up next line
+ getline(EventGeneratorFile, LineBuffer);
+
+ //Search for comment Symbol %
+ if (LineBuffer.compare(0, 1, "%") == 0) { /*do nothing*/
+ ;
+ }
+
+
+ //Search for Isotropic source
+ else if (LineBuffer.compare(0, 9, "Isotropic") == 0 && !check_Isotropic) {
+ check_Isotropic = true;
+ VEventGenerator* myEventGenerator = new EventGeneratorIsotropic();
+ EventGeneratorFile.close();
+ myEventGenerator->ReadConfiguration(Path);
+ EventGeneratorFile.open(Path.c_str());
+ myEventGenerator->InitializeRootOutput();
+ m_EventGenerator.push_back(myEventGenerator);
+ }
+
+ //Search for Beam
+ else if (LineBuffer.compare(0, 4, "Beam") == 0 && !check_Beam) {
+ check_Beam = true;
+ VEventGenerator* myEventGenerator = new EventGeneratorBeam();
+ EventGeneratorFile.close();
+ myEventGenerator->ReadConfiguration(Path);
+ EventGeneratorFile.open(Path.c_str());
+ myEventGenerator->InitializeRootOutput();
+ myEventGenerator->SetTarget(m_detector->GetTarget());
+ m_EventGenerator.push_back(myEventGenerator);
+ }
+
+ //Search for Transfert
+ else if (LineBuffer.compare(0, 9, "Transfert") == 0 && !check_Transfert && LineBuffer.compare(0, 11, "TransfertTo") != 0) {
+ check_Transfert = true;
+ VEventGenerator* myEventGenerator = new EventGeneratorTransfert();
+ EventGeneratorFile.close();
+ myEventGenerator->ReadConfiguration(Path);
+ EventGeneratorFile.open(Path.c_str());
+ myEventGenerator->InitializeRootOutput();
+ myEventGenerator->SetTarget(m_detector->GetTarget());
+ m_EventGenerator.push_back(myEventGenerator);
+ }
+
+ //Search for GammaDecay
+ else if ( LineBuffer.compare(0, 10, "GammaDecay") == 0 ) {
+ seenToken_GammaDecay++;
+ if (seenToken_GammaDecay>alreadyiInstantiate_GammaDecay) {
+ alreadyiInstantiate_GammaDecay++;
+ VEventGenerator* myEventGenerator = new EventGeneratorGammaDecay();
+ EventGeneratorFile.close();
+ myEventGenerator->ReadConfiguration(Path,alreadyiInstantiate_GammaDecay);
+ EventGeneratorFile.open(Path.c_str());
+ myEventGenerator->InitializeRootOutput();
+ myEventGenerator->SetTarget(m_detector->GetTarget());
+ m_EventGenerator.push_back(myEventGenerator);
+ seenToken_GammaDecay=0;
+ }
+
+ }
+
+ //Search for ParticleDecay
+ else if ( LineBuffer.compare(0, 13, "ParticleDecay") == 0 ) {
+ seenToken_ParticleDecay++;
+ if(seenToken_ParticleDecay>alreadyiInstantiate_ParticleDecay){
+ alreadyiInstantiate_ParticleDecay++;
+ VEventGenerator* myEventGenerator = new EventGeneratorParticleDecay();
+ EventGeneratorFile.close();
+ myEventGenerator->ReadConfiguration(Path,alreadyiInstantiate_ParticleDecay);
+ EventGeneratorFile.open(Path.c_str());
+ myEventGenerator->InitializeRootOutput();
+ myEventGenerator->SetTarget(m_detector->GetTarget());
+ m_EventGenerator.push_back(myEventGenerator);
+ seenToken_ParticleDecay=0;
+ }
+ }
+
+
+ //Search for Transfert To Resonance
+ else if (LineBuffer.compare(0, 10, "PhaseSpace") == 0 && !check_PhaseSpace) {
+ check_PhaseSpace = true;
+ VEventGenerator* myEventGenerator = new EventGeneratorPhaseSpace();
+ EventGeneratorFile.close();
+ myEventGenerator->ReadConfiguration(Path);
+ EventGeneratorFile.open(Path.c_str());
+ myEventGenerator->InitializeRootOutput();
+ myEventGenerator->SetTarget(m_detector->GetTarget());
+ m_EventGenerator.push_back(myEventGenerator);
+ }
+ }
+ EventGeneratorFile.close();
}
--
GitLab