source: trunk/examples/advanced/radioprotection/src/RemSimHadronicBertini.cc @ 1292

Last change on this file since 1292 was 1230, checked in by garnier, 14 years ago

update to geant4.9.3
File size: 12.0 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer                                           *
4// *                                                                  *
5// * The  Geant4 software  is  copyright of the Copyright Holders  of *
6// * the Geant4 Collaboration.  It is provided  under  the terms  and *
7// * conditions of the Geant4 Software License,  included in the file *
8// * LICENSE and available at  http://cern.ch/geant4/license .  These *
9// * include a list of copyright holders.                             *
10// *                                                                  *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work  make  any representation or  warranty, express or implied, *
14// * regarding  this  software system or assume any liability for its *
15// * use.  Please see the license in the file  LICENSE  and URL above *
16// * for the full disclaimer and the limitation of liability.         *
17// *                                                                  *
18// * This  code  implementation is the result of  the  scientific and *
19// * technical work of the GEANT4 collaboration.                      *
20// * By using,  copying,  modifying or  distributing the software (or *
21// * any work based  on the software)  you  agree  to acknowledge its *
22// * use  in  resulting  scientific  publications,  and indicate your *
23// * acceptance of all terms of the Geant4 Software license.          *
24// ********************************************************************
25//
26//    **************************************
27//    *                                    *
28//    *    RemSimHadronicBertini.cc        *
29//    *                                    *
30//    **************************************
31//
32// $Id: RemSimHadronicBertini.cc,v 1.6 2009/11/12 05:12:18 cirrone Exp $
33// GEANT4 tag $Name: geant4-09-03-cand-01 $
34//
35// Author : Susanna Guatelli, guatelli@ge.infn.it
36
37// Code review: MGP, 7 November 2006 (to be completed)
38//
39#include "RemSimHadronicBertini.hh"
40#include "G4BinaryLightIonReaction.hh"
41#include "G4TripathiCrossSection.hh"
42#include "G4IonsShenCrossSection.hh"
43#include "G4ParticleDefinition.hh"
44#include "G4ProcessManager.hh"
45#include "G4LElastic.hh"
46#include "G4CascadeInterface.hh"
47//#include "G4PreCompoundModel.hh"
48//#include "G4QGSModel.hh"
49#include "G4LEProtonInelastic.hh"
50#include "G4LENeutronInelastic.hh"
51#include "G4LEPionPlusInelastic.hh"
52#include "G4LEPionMinusInelastic.hh"
53#include "G4LEAlphaInelastic.hh"
54#include "G4LFission.hh"
55#include "G4LCapture.hh"
56//#include "G4TheoFSGenerator.hh"
57//#include "G4GeneratorPrecompoundInterface.hh"
58//#include "G4PreCompoundModel.hh"
59//#include "G4QGSMFragmentation.hh"
60//#include "G4ExcitedStringDecay.hh"
61#include "G4HadronElasticProcess.hh"
62#include "G4HadronFissionProcess.hh"
63#include "G4HadronCaptureProcess.hh"
64#include "G4ProtonInelasticProcess.hh"
65#include "G4NeutronInelasticProcess.hh"
66#include "G4PionPlusInelasticProcess.hh"
67#include "G4PionMinusInelasticProcess.hh"
68#include "G4AlphaInelasticProcess.hh"
69#include "G4ExcitationHandler.hh"
70
71//
72// BERTINI PHYSICS LIST
73//
74
75RemSimHadronicBertini::RemSimHadronicBertini(const G4String& name): 
76G4VPhysicsConstructor(name)
77{}
78
79RemSimHadronicBertini::~RemSimHadronicBertini()
80{}
81
82void RemSimHadronicBertini::ConstructProcess()
83{
84  // Physics for proton, neutron, pion+ and pion-
85 
86  // Elastic scattering: LElastic model
87  G4LElastic* elasticModel = new G4LElastic();
88  G4HadronElasticProcess* elasticScattering = new G4HadronElasticProcess();
89  elasticScattering->RegisterMe(elasticModel);
90
91  // Inelastic scattering: Bertini model up to 3.2 GeV
92  G4CascadeInterface* bertiniModel = new G4CascadeInterface();
93  // Energy limit of the Bertini model
94  G4double bertiniHighEnergyLimit = 3.2 * GeV;
95  bertiniModel->SetMaxEnergy(bertiniHighEnergyLimit);
96
97  // Inelastic scattering: LEP model between 2.8 * GeV and 25. * GeV
98  G4LEProtonInelastic* LEPProtonModel = new G4LEProtonInelastic();
99  G4LENeutronInelastic* LEPNeutronModel = new G4LENeutronInelastic();
100  G4LEPionPlusInelastic* LEPPionPlusModel = new G4LEPionPlusInelastic();
101  G4LEPionMinusInelastic* LEPPionMinusModel = new G4LEPionMinusInelastic();
102  // Set the energy limits
103  G4double LEPLowEnergyLimit = 2.8 * GeV;
104  G4double LEPHighEnergyLimit = 25. * GeV;
105  LEPProtonModel->SetMinEnergy(LEPLowEnergyLimit);
106  LEPProtonModel->SetMaxEnergy(LEPHighEnergyLimit);
107  LEPNeutronModel->SetMinEnergy(LEPLowEnergyLimit);
108  LEPNeutronModel->SetMaxEnergy(LEPHighEnergyLimit);
109  LEPPionPlusModel->SetMinEnergy(LEPLowEnergyLimit);
110  LEPPionPlusModel->SetMaxEnergy(LEPHighEnergyLimit);
111  LEPPionMinusModel->SetMinEnergy(LEPLowEnergyLimit);
112  LEPPionMinusModel->SetMaxEnergy(LEPHighEnergyLimit);
113
114  // Inelastic scattering: QGSP model between 20 GeV and 100 TeV
115  QGSPModel = new G4TheoFSGenerator();
116  // Set the energy limits
117  G4double QGSPLowEnergyLimit = 20.* GeV;
118  G4double QGSPHighEnergyLimit = 100.* GeV;
119  QGSPModel->SetMinEnergy(QGSPLowEnergyLimit);
120  QGSPModel->SetMaxEnergy(QGSPHighEnergyLimit);
121
122  theCascade = new G4GeneratorPrecompoundInterface();
123
124  thePreEquilib = new G4PreCompoundModel(&theHandler);
125  theCascade->SetDeExcitation(thePreEquilib);
126  QGSPModel->SetTransport(theCascade);
127
128  // Set the fragmentation
129  theFragmentation = new G4QGSMFragmentation();
130  theStringDecay = new G4ExcitedStringDecay(theFragmentation);
131  theStringModel.SetFragmentationModel(theStringDecay);
132  QGSPModel->SetHighEnergyGenerator(&theStringModel);
133
134  // ---------------------------------------------------------------------------------------------
135  // Proton processes
136  G4ParticleDefinition* proton = G4Proton::ProtonDefinition();
137  G4ProcessManager* protonProcessManager = proton->GetProcessManager();
138 
139  // Proton inelastic scattering
140  G4ProtonInelasticProcess* protonInelasticProcess = new G4ProtonInelasticProcess();
141  // Activate the cross-sections for proton nuclear scattering up to 20 GeV
142  //G4ProtonInelasticCrossSection protonCrossSection;
143  protonInelasticProcess->AddDataSet(&protonCrossSection);
144  // Set the models
145  protonInelasticProcess->RegisterMe(bertiniModel);
146  protonInelasticProcess->RegisterMe(LEPProtonModel);
147 
148  protonInelasticProcess->RegisterMe(QGSPModel);
149  // Activate the inelastic scattering
150  protonProcessManager->AddDiscreteProcess(protonInelasticProcess);
151  // Activate the elastic scattering
152  protonProcessManager->AddDiscreteProcess(elasticScattering);
153
154  //------------------------------------------------------
155  // Pion Plus processes
156  G4ParticleDefinition* piPlus = G4PionPlus::PionPlusDefinition();
157  G4ProcessManager* pionPlusProcessManager = piPlus->GetProcessManager();
158 
159  // Define the inelastic scattering for pion plus
160  G4PionPlusInelasticProcess* pionPlusInelasticProcess = new G4PionPlusInelasticProcess();
161  // Set the cross section
162 // G4PiNuclearCrossSection pionCrossSection;
163  pionPlusInelasticProcess->AddDataSet(&pionCrossSection);
164  // Register the models
165  pionPlusInelasticProcess->RegisterMe(bertiniModel);
166  pionPlusInelasticProcess->RegisterMe(LEPPionPlusModel);
167  pionPlusInelasticProcess->RegisterMe(QGSPModel);
168  // Activate the inelastic scattering
169  pionPlusProcessManager->AddDiscreteProcess(pionPlusInelasticProcess);
170  // Activate the elastic process
171  pionPlusProcessManager->AddDiscreteProcess(elasticScattering);
172
173  //------------------------------------------------------------
174  // Pion Minus processes
175  G4ParticleDefinition* piMinus = G4PionMinus::PionMinusDefinition();
176  G4ProcessManager* pionMinusProcessManager = piMinus->GetProcessManager();
177 
178  // Define the inelastic processes for pion minus
179  G4PionMinusInelasticProcess* pionMinusInelasticProcess = new G4PionMinusInelasticProcess();
180  // Set the cross section
181  pionMinusInelasticProcess->AddDataSet(&pionCrossSection);
182  // Register the models
183  pionMinusInelasticProcess->RegisterMe(bertiniModel);
184  pionMinusInelasticProcess->RegisterMe(LEPPionMinusModel);
185  pionMinusInelasticProcess->RegisterMe(QGSPModel);
186  // Activate the inelastic scattering
187  pionMinusProcessManager->AddDiscreteProcess(pionMinusInelasticProcess);
188  // Activate the elastic scattering
189  pionMinusProcessManager->AddDiscreteProcess(elasticScattering);
190 
191  //-----------------------------------------------------
192  // Neutron processes
193  G4ParticleDefinition* neutron = G4Neutron::NeutronDefinition();
194  G4ProcessManager* neutronProcessManager = neutron->GetProcessManager();
195
196  // Inelastic process
197  G4NeutronInelasticProcess* neutronInelasticProcess = new G4NeutronInelasticProcess();
198  // Set the cross section
199  //G4NeutronInelasticCrossSection neutronCrossSection;
200  neutronInelasticProcess->AddDataSet(&neutronCrossSection);
201  // Set the models
202  neutronInelasticProcess->RegisterMe(bertiniModel);
203  neutronInelasticProcess->RegisterMe(LEPNeutronModel);
204  neutronInelasticProcess->RegisterMe(QGSPModel);
205 // Activate the neutron inelastic scattering
206  neutronProcessManager->AddDiscreteProcess(neutronInelasticProcess);
207  // Activate the neutron elastic scattering
208  neutronProcessManager->AddDiscreteProcess(elasticScattering);
209
210  // Neutron capture process
211  G4HadronCaptureProcess* neutronCapture = new G4HadronCaptureProcess();
212  // Final state production model for capture of neutral hadrons in nuclei
213  G4LCapture* captureModel = new G4LCapture();
214  // Set the energy range for the capture model
215  G4double neutronLowEnergyLimit = 0. * MeV;
216  G4double neutronHighEnergyLimit = 100. * TeV;
217  captureModel->SetMinEnergy(neutronLowEnergyLimit);
218  captureModel->SetMaxEnergy(neutronHighEnergyLimit);
219  // Activate the neutron capture model
220  neutronCapture->RegisterMe(captureModel);
221  // Activate the neutron capture process
222  neutronProcessManager->AddDiscreteProcess(neutronCapture);
223
224  // Process for induced fission
225   G4HadronFissionProcess* fission = new G4HadronFissionProcess();
226  //Final state production model for induced fission
227  G4LFission* fissionModel = new G4LFission();
228  // Set the energy range for the fission model
229  fissionModel->SetMinEnergy(neutronLowEnergyLimit);
230  fissionModel->SetMaxEnergy(neutronHighEnergyLimit);
231  // Register the fission model
232  fission->RegisterMe(fissionModel); 
233  // Activate the fission process
234  neutronProcessManager->AddDiscreteProcess(fission); 
235
236  //--------------------------------------------------------
237  // Physics for alpha particles
238 
239  G4ParticleDefinition* alpha = G4Alpha::AlphaDefinition();
240  G4ProcessManager* alphaProcessManager = alpha->GetProcessManager();
241
242  // Cross section data sets
243 
244  // TRIPATHI CROSS SECTION
245  // Implementation of formulas taken from NASA technical paper 3621 by
246  // Tripathi, et al. Cross-sections for ion ion scattering
247  G4TripathiCrossSection* tripathi = new G4TripathiCrossSection();
248
249  // IONS SHEN CROSS SECTION
250  // Implementation of formulas
251  // Shen et al. Nuc. Phys. A 491 130 (1989)
252  // Total Reaction Cross Section for Heavy-Ion Collisions
253  G4IonsShenCrossSection* shen = new G4IonsShenCrossSection();
254
255  G4LEAlphaInelastic* LEPAlphaModel = new G4LEAlphaInelastic();
256  // Energy limit of the LEP model for alpha particles
257  G4double LEPAlphaHighLimit = 100 * MeV;
258  LEPAlphaModel->SetMaxEnergy(LEPAlphaHighLimit);
259
260  G4BinaryLightIonReaction* binaryIonModel = new G4BinaryLightIonReaction();
261  // Energy limit of the binary ion model
262  G4double binaryIonLowLimit = 80. * MeV;
263  G4double binaryIonHighLimit = 400. * GeV; 
264  binaryIonModel->SetMinEnergy(binaryIonLowLimit);
265  binaryIonModel->SetMaxEnergy(binaryIonHighLimit);
266
267  // Define the alpha inelastic scattering
268  G4AlphaInelasticProcess* alphaInelasticProcess = new G4AlphaInelasticProcess();
269  // Activate the Tripathi and Shen Cross Section
270  alphaInelasticProcess->AddDataSet(tripathi);
271  alphaInelasticProcess->AddDataSet(shen);
272  // Set the models
273  alphaInelasticProcess->RegisterMe(LEPAlphaModel);
274  alphaInelasticProcess->RegisterMe(binaryIonModel);
275  // Activate the inelastic scattering
276  alphaProcessManager->AddDiscreteProcess(alphaInelasticProcess);
277  // Activate the elastic scattering
278  alphaProcessManager->AddDiscreteProcess(elasticScattering);
279}
280
281
282
Note: See TracBrowser for help on using the repository browser.