| Geant4 Cross Reference |
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 /// \file BiasingOperation.cc
27 /// \brief Implementation of the BiasingOperation class
28 //
29 //
30
31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
33
34 #include "BiasingOperation.hh"
35
36 #include "G4BGGNucleonInelasticXS.hh"
37 #include "G4BGGPionInelasticXS.hh"
38 #include "G4BiasingProcessInterface.hh"
39 #include "G4CascadeInterface.hh"
40 #include "G4CrossSectionDataStore.hh"
41 #include "G4ExcitedStringDecay.hh"
42 #include "G4FTFModel.hh"
43 #include "G4GeneratorPrecompoundInterface.hh"
44 #include "G4HadronInelasticProcess.hh"
45 #include "G4HadronicParameters.hh"
46 #include "G4INCLXXInterface.hh"
47 #include "G4LundStringFragmentation.hh"
48 #include "G4NeutronInelasticXS.hh"
49 #include "G4TheoFSGenerator.hh"
50 #include "G4VParticleChange.hh"
51
52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
53
54 BiasingOperation::BiasingOperation(G4String name) : G4VBiasingOperation(name)
55 {
56 // Create the inelastic processes for p , n , pi+ , pi-
57 fProtonInelasticProcess = new G4HadronInelasticProcess("protonInelastic", G4Proton::Definition());
58 fNeutronInelasticProcess =
59 new G4HadronInelasticProcess("neutronInelastic", G4Neutron::Definition());
60 fPionPlusInelasticProcess =
61 new G4HadronInelasticProcess("pi+Inelastic", G4PionPlus::Definition());
62 fPionMinusInelasticProcess =
63 new G4HadronInelasticProcess("pi-Inelastic", G4PionMinus::Definition());
64
65 // Set the energy ranges
66 const G4double maxBERT = 41.0 * CLHEP::MeV;
67 const G4double maxINCLXX = 12.0 * CLHEP::GeV;
68 const G4double minINCLXX = 40.0 * CLHEP::MeV;
69 const G4double minFTFP = 3.0 * CLHEP::GeV;
70 const G4double maxFTFP = G4HadronicParameters::Instance()->GetMaxEnergy();
71
72 // Create the hadronic models (to replace FTFP_BERT with "FTFP_INCLXX",
73 // keeping the same energy ranges for the transition between models).
74 // Notice that it is better to create the models here from scratch,
75 // instead of reusing the existing ones, because we might pick up the
76 // existing ones associated to the wrong particles...
77 // --- FTFP model ---
78 G4FTFModel* theStringModel = new G4FTFModel;
79 G4LundStringFragmentation* theLund = new G4LundStringFragmentation;
80 G4ExcitedStringDecay* theStringDecay = new G4ExcitedStringDecay(theLund);
81 theStringModel->SetFragmentationModel(theStringDecay);
82 G4GeneratorPrecompoundInterface* thePrecoInterface = new G4GeneratorPrecompoundInterface;
83 G4TheoFSGenerator* theHighEnergyModel = new G4TheoFSGenerator("FTFP");
84 theHighEnergyModel->SetHighEnergyGenerator(theStringModel);
85 theHighEnergyModel->SetTransport(thePrecoInterface);
86 theHighEnergyModel->SetMinEnergy(minFTFP);
87 theHighEnergyModel->SetMaxEnergy(maxFTFP);
88 // Bertini : create a new model to be used below INCLXX limit
89 G4CascadeInterface* theBertiniModel = new G4CascadeInterface();
90 theBertiniModel->SetMinEnergy(0.0);
91 theBertiniModel->SetMaxEnergy(maxBERT);
92 // --- Preco ---
93 // --- INCLXX model ---
94 // The instance for nucleons:
95 G4INCLXXInterface* theInclxxModel = new G4INCLXXInterface();
96 theInclxxModel->SetMinEnergy(minINCLXX);
97 theInclxxModel->SetMaxEnergy(maxINCLXX); // Use the same as for FTFP_BERT
98
99 // Register the models
100 fProtonInelasticProcess->RegisterMe(theHighEnergyModel);
101 fProtonInelasticProcess->RegisterMe(theInclxxModel);
102 fProtonInelasticProcess->RegisterMe(theBertiniModel);
103 fNeutronInelasticProcess->RegisterMe(theHighEnergyModel);
104 fNeutronInelasticProcess->RegisterMe(theInclxxModel);
105 fNeutronInelasticProcess->RegisterMe(theBertiniModel);
106 fPionPlusInelasticProcess->RegisterMe(theHighEnergyModel);
107 fPionPlusInelasticProcess->RegisterMe(theInclxxModel);
108 fPionPlusInelasticProcess->RegisterMe(theBertiniModel);
109 fPionMinusInelasticProcess->RegisterMe(theHighEnergyModel);
110 fPionMinusInelasticProcess->RegisterMe(theInclxxModel);
111 fPionMinusInelasticProcess->RegisterMe(theBertiniModel);
112
113 G4VCrossSectionDataSet* theProtonXSdata = new G4BGGNucleonInelasticXS(G4Proton::Definition());
114 theProtonXSdata->BuildPhysicsTable(*(G4Proton::Definition()));
115 fProtonInelasticProcess->AddDataSet(theProtonXSdata);
116 G4VCrossSectionDataSet* theNeutronXSdata = new G4NeutronInelasticXS;
117 theNeutronXSdata->BuildPhysicsTable(*(G4Neutron::Definition()));
118 fNeutronInelasticProcess->AddDataSet(theNeutronXSdata);
119 G4VCrossSectionDataSet* thePionPlusXSdata = new G4BGGPionInelasticXS(G4PionPlus::Definition());
120 thePionPlusXSdata->BuildPhysicsTable(*(G4PionPlus::Definition()));
121 fPionPlusInelasticProcess->AddDataSet(thePionPlusXSdata);
122 G4VCrossSectionDataSet* thePionMinusXSdata = new G4BGGPionInelasticXS(G4PionMinus::Definition());
123 thePionMinusXSdata->BuildPhysicsTable(*(G4PionMinus::Definition()));
124 fPionMinusInelasticProcess->AddDataSet(thePionMinusXSdata);
125 }
126
127 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
128
129 BiasingOperation::~BiasingOperation() {}
130
131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
132
133 G4VParticleChange* BiasingOperation::ApplyFinalStateBiasing(const G4BiasingProcessInterface*,
134 const G4Track* track,
135 const G4Step* step, G4bool&)
136 {
137 if (track->GetParticleDefinition() == G4Proton::Definition()) {
138 auto particle = track->GetDynamicParticle();
139 auto material = track->GetMaterial();
140 fProtonInelasticProcess->GetCrossSectionDataStore()->ComputeCrossSection(particle, material);
141 return fProtonInelasticProcess->PostStepDoIt(*track, *step);
142 }
143 else if (track->GetParticleDefinition() == G4Neutron::Definition()) {
144 auto particle = track->GetDynamicParticle();
145 auto material = track->GetMaterial();
146 fNeutronInelasticProcess->GetCrossSectionDataStore()->ComputeCrossSection(particle, material);
147 return fNeutronInelasticProcess->PostStepDoIt(*track, *step);
148 }
149 else if (track->GetParticleDefinition() == G4PionPlus::Definition()) {
150 auto particle = track->GetDynamicParticle();
151 auto material = track->GetMaterial();
152 fPionPlusInelasticProcess->GetCrossSectionDataStore()->ComputeCrossSection(particle, material);
153 return fPionPlusInelasticProcess->PostStepDoIt(*track, *step);
154 }
155 else if (track->GetParticleDefinition() == G4PionMinus::Definition()) {
156 auto particle = track->GetDynamicParticle();
157 auto material = track->GetMaterial();
158 fPionMinusInelasticProcess->GetCrossSectionDataStore()->ComputeCrossSection(particle, material);
159 return fPionMinusInelasticProcess->PostStepDoIt(*track, *step);
160 }
161 else {
162 G4cerr << "ERROR in BiasingOperation::ApplyFinalStateBiasing : unexpected particle = "
163 << track->GetParticleDefinition()->GetParticleName() << G4endl;
164 return 0;
165 }
166 }
167
168 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
169