| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 /// \file BiasingOperation.cc 26 /// \file BiasingOperation.cc
27 /// \brief Implementation of the BiasingOperat 27 /// \brief Implementation of the BiasingOperation class
28 // 28 //
29 // 29 //
30 30
31 //....oooOO0OOooo........oooOO0OOooo........oo 31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 //....oooOO0OOooo........oooOO0OOooo........oo 32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
33 33
34 #include "BiasingOperation.hh" 34 #include "BiasingOperation.hh"
35 <<
36 #include "G4BGGNucleonInelasticXS.hh" <<
37 #include "G4BGGPionInelasticXS.hh" <<
38 #include "G4BiasingProcessInterface.hh" 35 #include "G4BiasingProcessInterface.hh"
39 #include "G4CascadeInterface.hh" << 36 #include "G4VParticleChange.hh"
40 #include "G4CrossSectionDataStore.hh" << 37 #include "G4ProtonInelasticProcess.hh"
41 #include "G4ExcitedStringDecay.hh" << 38 #include "G4NeutronInelasticProcess.hh"
42 #include "G4FTFModel.hh" << 39 #include "G4PionPlusInelasticProcess.hh"
43 #include "G4GeneratorPrecompoundInterface.hh" << 40 #include "G4PionMinusInelasticProcess.hh"
44 #include "G4HadronInelasticProcess.hh" <<
45 #include "G4HadronicParameters.hh" 41 #include "G4HadronicParameters.hh"
46 #include "G4INCLXXInterface.hh" << 42 #include "G4FTFModel.hh"
47 #include "G4LundStringFragmentation.hh" 43 #include "G4LundStringFragmentation.hh"
48 #include "G4NeutronInelasticXS.hh" << 44 #include "G4ExcitedStringDecay.hh"
>> 45 #include "G4GeneratorPrecompoundInterface.hh"
49 #include "G4TheoFSGenerator.hh" 46 #include "G4TheoFSGenerator.hh"
50 #include "G4VParticleChange.hh" << 47 #include "G4HadronicInteractionRegistry.hh"
>> 48 #include "G4VPreCompoundModel.hh"
>> 49 #include "G4INCLXXInterface.hh"
>> 50 #include "G4HadronInelasticDataSet.hh"
51 51
52 //....oooOO0OOooo........oooOO0OOooo........oo 52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
53 53
54 BiasingOperation::BiasingOperation(G4String na << 54 BiasingOperation::BiasingOperation( G4String name ) : G4VBiasingOperation( name ) {
55 { << 55
56 // Create the inelastic processes for p , n << 56 // Create the inelastic processes for p , n , pi+ , pi-
57 fProtonInelasticProcess = new G4HadronInelas << 57 fProtonInelasticProcess = new G4ProtonInelasticProcess;
58 fNeutronInelasticProcess = << 58 fNeutronInelasticProcess = new G4NeutronInelasticProcess;
59 new G4HadronInelasticProcess("neutronInela << 59 fPionPlusInelasticProcess = new G4PionPlusInelasticProcess;
60 fPionPlusInelasticProcess = << 60 fPionMinusInelasticProcess = new G4PionMinusInelasticProcess;
61 new G4HadronInelasticProcess("pi+Inelastic <<
62 fPionMinusInelasticProcess = <<
63 new G4HadronInelasticProcess("pi-Inelastic <<
64 61
65 // Set the energy ranges 62 // Set the energy ranges
66 const G4double maxBERT = 41.0 * CLHEP::MeV; << 63 const G4double minPreco = 0.0;
67 const G4double maxINCLXX = 12.0 * CLHEP::GeV << 64 const G4double maxPreco = 2.0 * CLHEP::MeV;
68 const G4double minINCLXX = 40.0 * CLHEP::MeV << 65 const G4double maxBERT = 12.0 * CLHEP::GeV;
>> 66 const G4double minINCLXX = 1.0 * CLHEP::MeV;
69 const G4double minFTFP = 3.0 * CLHEP::GeV; 67 const G4double minFTFP = 3.0 * CLHEP::GeV;
70 const G4double maxFTFP = G4HadronicParameter 68 const G4double maxFTFP = G4HadronicParameters::Instance()->GetMaxEnergy();
71 69
72 // Create the hadronic models (to replace FT << 70 // Create the hadronic models (to replace FTFP_BERT with "FTFP_INCLXX",
73 // keeping the same energy ranges for the tr 71 // keeping the same energy ranges for the transition between models).
74 // Notice that it is better to create the mo 72 // Notice that it is better to create the models here from scratch,
75 // instead of reusing the existing ones, bec 73 // instead of reusing the existing ones, because we might pick up the
76 // existing ones associated to the wrong par 74 // existing ones associated to the wrong particles...
77 // --- FTFP model --- 75 // --- FTFP model ---
78 G4FTFModel* theStringModel = new G4FTFModel; 76 G4FTFModel* theStringModel = new G4FTFModel;
79 G4LundStringFragmentation* theLund = new G4L 77 G4LundStringFragmentation* theLund = new G4LundStringFragmentation;
80 G4ExcitedStringDecay* theStringDecay = new G << 78 G4ExcitedStringDecay* theStringDecay = new G4ExcitedStringDecay( theLund );
81 theStringModel->SetFragmentationModel(theStr << 79 theStringModel->SetFragmentationModel( theStringDecay );
82 G4GeneratorPrecompoundInterface* thePrecoInt 80 G4GeneratorPrecompoundInterface* thePrecoInterface = new G4GeneratorPrecompoundInterface;
83 G4TheoFSGenerator* theHighEnergyModel = new << 81 G4TheoFSGenerator* theHighEnergyModel = new G4TheoFSGenerator( "FTFP" );
84 theHighEnergyModel->SetHighEnergyGenerator(t << 82 theHighEnergyModel->SetHighEnergyGenerator( theStringModel );
85 theHighEnergyModel->SetTransport(thePrecoInt << 83 theHighEnergyModel->SetTransport( thePrecoInterface );
86 theHighEnergyModel->SetMinEnergy(minFTFP); << 84 theHighEnergyModel->SetMinEnergy( minFTFP );
87 theHighEnergyModel->SetMaxEnergy(maxFTFP); << 85 theHighEnergyModel->SetMaxEnergy( maxFTFP );
88 // Bertini : create a new model to be used b << 86 // Preco : create a new model to be used only for INCLXX for nucleons
89 G4CascadeInterface* theBertiniModel = new G4 << 87 G4VPreCompoundModel* thePreCompoundModel = new G4PreCompoundModel;
90 theBertiniModel->SetMinEnergy(0.0); << 88 thePreCompoundModel->SetMinEnergy( minPreco );
91 theBertiniModel->SetMaxEnergy(maxBERT); << 89 thePreCompoundModel->SetMaxEnergy( maxPreco );
92 // --- Preco --- 90 // --- Preco ---
93 // --- INCLXX model --- 91 // --- INCLXX model ---
94 // The instance for nucleons: 92 // The instance for nucleons:
95 G4INCLXXInterface* theInclxxModel = new G4IN << 93 G4INCLXXInterface* theInclxxModel = new G4INCLXXInterface( thePreCompoundModel );
96 theInclxxModel->SetMinEnergy(minINCLXX); << 94 theInclxxModel->SetMinEnergy( minINCLXX );
97 theInclxxModel->SetMaxEnergy(maxINCLXX); // << 95 theInclxxModel->SetMaxEnergy( maxBERT ); // Use the same as for Bertini
>> 96 // The instance for pions:
>> 97 G4INCLXXInterface* theInclxxModel_forPions = new G4INCLXXInterface;
>> 98 theInclxxModel_forPions->SetMinEnergy( minPreco );
>> 99 theInclxxModel_forPions->SetMaxEnergy( maxBERT ); // Use the same as for Bertini
98 100
99 // Register the models 101 // Register the models
100 fProtonInelasticProcess->RegisterMe(theHighE << 102 fProtonInelasticProcess->RegisterMe( theHighEnergyModel );
101 fProtonInelasticProcess->RegisterMe(theInclx << 103 fProtonInelasticProcess->RegisterMe( theInclxxModel );
102 fProtonInelasticProcess->RegisterMe(theBerti << 104 fProtonInelasticProcess->RegisterMe( thePreCompoundModel );
103 fNeutronInelasticProcess->RegisterMe(theHigh << 105 fNeutronInelasticProcess->RegisterMe( theHighEnergyModel );
104 fNeutronInelasticProcess->RegisterMe(theIncl << 106 fNeutronInelasticProcess->RegisterMe( theInclxxModel );
105 fNeutronInelasticProcess->RegisterMe(theBert << 107 fNeutronInelasticProcess->RegisterMe( thePreCompoundModel );
106 fPionPlusInelasticProcess->RegisterMe(theHig << 108 fPionPlusInelasticProcess->RegisterMe( theHighEnergyModel );
107 fPionPlusInelasticProcess->RegisterMe(theInc << 109 fPionPlusInelasticProcess->RegisterMe( theInclxxModel_forPions );
108 fPionPlusInelasticProcess->RegisterMe(theBer << 110 fPionMinusInelasticProcess->RegisterMe( theHighEnergyModel );
109 fPionMinusInelasticProcess->RegisterMe(theHi << 111 fPionMinusInelasticProcess->RegisterMe( theInclxxModel_forPions );
110 fPionMinusInelasticProcess->RegisterMe(theIn << 112
111 fPionMinusInelasticProcess->RegisterMe(theBe << 113 // Register the cross sections: this is mandatory starting from G4 10.6
112 << 114 // because the default Gheisha inelastic cross sections have been removed.
113 G4VCrossSectionDataSet* theProtonXSdata = ne << 115 // It is convenient to use the Gheisha inelastic cross sections here
114 theProtonXSdata->BuildPhysicsTable(*(G4Proto << 116 // because they do not require any special initialization.
115 fProtonInelasticProcess->AddDataSet(theProto << 117 fProtonInelasticProcess->AddDataSet( new G4HadronInelasticDataSet );
116 G4VCrossSectionDataSet* theNeutronXSdata = n << 118 fNeutronInelasticProcess->AddDataSet( new G4HadronInelasticDataSet );
117 theNeutronXSdata->BuildPhysicsTable(*(G4Neut << 119 fPionPlusInelasticProcess->AddDataSet( new G4HadronInelasticDataSet );
118 fNeutronInelasticProcess->AddDataSet(theNeut << 120 fPionMinusInelasticProcess->AddDataSet( new G4HadronInelasticDataSet );
119 G4VCrossSectionDataSet* thePionPlusXSdata = <<
120 thePionPlusXSdata->BuildPhysicsTable(*(G4Pio <<
121 fPionPlusInelasticProcess->AddDataSet(thePio <<
122 G4VCrossSectionDataSet* thePionMinusXSdata = <<
123 thePionMinusXSdata->BuildPhysicsTable(*(G4Pi <<
124 fPionMinusInelasticProcess->AddDataSet(thePi <<
125 } 121 }
126 122
127 //....oooOO0OOooo........oooOO0OOooo........oo 123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
128 124
129 BiasingOperation::~BiasingOperation() {} 125 BiasingOperation::~BiasingOperation() {}
130 126
131 //....oooOO0OOooo........oooOO0OOooo........oo 127 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
132 128
133 G4VParticleChange* BiasingOperation::ApplyFina << 129 G4VParticleChange* BiasingOperation::
134 << 130 ApplyFinalStateBiasing( const G4BiasingProcessInterface* ,
135 << 131 const G4Track* track, const G4Step* step, G4bool& ) {
136 { << 132 if ( track->GetParticleDefinition() == G4Proton::Definition() ) {
137 if (track->GetParticleDefinition() == G4Prot << 133 return fProtonInelasticProcess->PostStepDoIt( *track, *step );
138 auto particle = track->GetDynamicParticle( << 134 } else if ( track->GetParticleDefinition() == G4Neutron::Definition() ) {
139 auto material = track->GetMaterial(); << 135 return fNeutronInelasticProcess->PostStepDoIt( *track, *step );
140 fProtonInelasticProcess->GetCrossSectionDa << 136 } else if ( track->GetParticleDefinition() == G4PionPlus::Definition() ) {
141 return fProtonInelasticProcess->PostStepDo << 137 return fPionPlusInelasticProcess->PostStepDoIt( *track, *step );
142 } << 138 } else if ( track->GetParticleDefinition() == G4PionMinus::Definition() ) {
143 else if (track->GetParticleDefinition() == G << 139 return fPionMinusInelasticProcess->PostStepDoIt( *track, *step );
144 auto particle = track->GetDynamicParticle( << 140 } else {
145 auto material = track->GetMaterial(); <<
146 fNeutronInelasticProcess->GetCrossSectionD <<
147 return fNeutronInelasticProcess->PostStepD <<
148 } <<
149 else if (track->GetParticleDefinition() == G <<
150 auto particle = track->GetDynamicParticle( <<
151 auto material = track->GetMaterial(); <<
152 fPionPlusInelasticProcess->GetCrossSection <<
153 return fPionPlusInelasticProcess->PostStep <<
154 } <<
155 else if (track->GetParticleDefinition() == G <<
156 auto particle = track->GetDynamicParticle( <<
157 auto material = track->GetMaterial(); <<
158 fPionMinusInelasticProcess->GetCrossSectio <<
159 return fPionMinusInelasticProcess->PostSte <<
160 } <<
161 else { <<
162 G4cerr << "ERROR in BiasingOperation::Appl 141 G4cerr << "ERROR in BiasingOperation::ApplyFinalStateBiasing : unexpected particle = "
163 << track->GetParticleDefinition()-> 142 << track->GetParticleDefinition()->GetParticleName() << G4endl;
164 return 0; 143 return 0;
165 } 144 }
166 } 145 }
167 146
168 //....oooOO0OOooo........oooOO0OOooo........oo 147 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
169 148