| 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 // neutron_hp -- source file 26 // neutron_hp -- source file
27 // J.P. Wellisch, Nov-1996 27 // J.P. Wellisch, Nov-1996
28 // A prototype of the low energy neutron trans 28 // A prototype of the low energy neutron transport model.
29 // 29 //
30 // P. Arce, June-2014 Conversion neutron_hp to 30 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
31 // 31 //
32 #include "G4ParticleHPFSFissionFS.hh" 32 #include "G4ParticleHPFSFissionFS.hh"
33 <<
34 #include "G4Alpha.hh" <<
35 #include "G4Deuteron.hh" <<
36 #include "G4LorentzVector.hh" <<
37 #include "G4Nucleus.hh" <<
38 #include "G4ParticleHPDataUsed.hh" <<
39 #include "G4ParticleHPManager.hh" 33 #include "G4ParticleHPManager.hh"
40 #include "G4Poisson.hh" <<
41 #include "G4Proton.hh" <<
42 #include "G4ReactionProduct.hh" 34 #include "G4ReactionProduct.hh"
43 #include "G4ThreeVector.hh" << 35 #include "G4Nucleus.hh"
>> 36 #include "G4Proton.hh"
>> 37 #include "G4Deuteron.hh"
44 #include "G4Triton.hh" 38 #include "G4Triton.hh"
>> 39 #include "G4Alpha.hh"
>> 40 #include "G4ThreeVector.hh"
>> 41 #include "G4Poisson.hh"
>> 42 #include "G4LorentzVector.hh"
>> 43 #include "G4ParticleHPDataUsed.hh"
45 44
46 void G4ParticleHPFSFissionFS::Init(G4double A, << 45 void G4ParticleHPFSFissionFS::Init (G4double A, G4double Z, G4int M, G4String & dirName, G4String &, G4ParticleDefinition* )
47 const G4Str <<
48 { <<
49 G4String tString = "/FS/"; <<
50 G4bool dbool; <<
51 const G4ParticleHPDataUsed& aFile = <<
52 theNames.GetName(static_cast<G4int>(A), st <<
53 const G4String& filename = aFile.GetName(); <<
54 SetAZMs(A, Z, M, aFile); <<
55 if (!dbool) { <<
56 hasAnyData = false; <<
57 hasFSData = false; <<
58 hasXsec = false; <<
59 return; <<
60 } <<
61 <<
62 std::istringstream theData(std::ios::in); <<
63 G4ParticleHPManager::GetInstance()->GetDataS <<
64 <<
65 G4int infoType, dataType; <<
66 hasFSData = false; <<
67 while (theData >> infoType) // Loop checkin <<
68 { 46 {
69 hasFSData = true; << 47 G4String tString = "/FS/";
70 theData >> dataType; << 48 G4bool dbool;
71 switch (infoType) { << 49 G4ParticleHPDataUsed aFile = theNames.GetName(static_cast<G4int>(A), static_cast<G4int>(Z), M, dirName, tString, dbool);
72 case 1: << 50 G4String filename = aFile.GetName();
73 if (dataType == 4) theNeutronAngularDi << 51 SetAZMs( A, Z, M, aFile );
74 if (dataType == 5) thePromptNeutronEnD << 52 if(!dbool)
75 if (dataType == 12) theFinalStatePhoto << 53 {
76 if (dataType == 14) theFinalStatePhoto << 54 hasAnyData = false;
77 if (dataType == 15) theFinalStatePhoto << 55 hasFSData = false;
78 break; << 56 hasXsec = false;
79 case 2: << 57 return;
80 if (dataType == 1) theFinalStateNeutro <<
81 break; <<
82 case 3: <<
83 if (dataType == 1) theFinalStateNeutro <<
84 if (dataType == 5) theDelayedNeutronEn <<
85 break; <<
86 case 4: <<
87 if (dataType == 1) theFinalStateNeutro <<
88 break; <<
89 case 5: <<
90 if (dataType == 1) theEnergyRelease.In <<
91 break; <<
92 default: <<
93 G4cout << "G4ParticleHPFSFissionFS::In <<
94 throw G4HadronicException(__FILE__, __ <<
95 "G4ParticleH <<
96 break; <<
97 } 58 }
>> 59 //std::ifstream theData(filename, std::ios::in);
>> 60 std::istringstream theData(std::ios::in);
>> 61 G4ParticleHPManager::GetInstance()->GetDataStream(filename,theData);
>> 62
>> 63 // here it comes
>> 64 G4int infoType, dataType;
>> 65 hasFSData = false;
>> 66 while (theData >> infoType) // Loop checking, 11.05.2015, T. Koi
>> 67 {
>> 68 hasFSData = true;
>> 69 theData >> dataType;
>> 70 switch(infoType)
>> 71 {
>> 72 case 1:
>> 73 if(dataType==4) theNeutronAngularDis.Init(theData);
>> 74 if(dataType==5) thePromptNeutronEnDis.Init(theData);
>> 75 if(dataType==12) theFinalStatePhotons.InitMean(theData);
>> 76 if(dataType==14) theFinalStatePhotons.InitAngular(theData);
>> 77 if(dataType==15) theFinalStatePhotons.InitEnergies(theData);
>> 78 break;
>> 79 case 2:
>> 80 if(dataType==1) theFinalStateNeutrons.InitMean(theData);
>> 81 break;
>> 82 case 3:
>> 83 if(dataType==1) theFinalStateNeutrons.InitDelayed(theData);
>> 84 if(dataType==5) theDelayedNeutronEnDis.Init(theData);
>> 85 break;
>> 86 case 4:
>> 87 if(dataType==1) theFinalStateNeutrons.InitPrompt(theData);
>> 88 break;
>> 89 case 5:
>> 90 if(dataType==1) theEnergyRelease.Init(theData);
>> 91 break;
>> 92 default:
>> 93 G4cout << "G4ParticleHPFSFissionFS::Init: unknown data type"<<dataType<<G4endl;
>> 94 throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPFSFissionFS::Init: unknown data type");
>> 95 break;
>> 96 }
>> 97 }
>> 98 //targetMass = theFinalStateNeutrons.GetTargetMass();
>> 99 //theData.close();
98 } 100 }
99 } << 101
100 << 102
101 G4DynamicParticleVector* G4ParticleHPFSFission << 103 G4DynamicParticleVector * G4ParticleHPFSFissionFS::ApplyYourself(G4int nPrompt,
102 << 104 G4int nDelayed, G4double * theDecayConst)
103 { << 105 {
104 G4int i; << 106 G4int i;
105 auto aResult = new G4DynamicParticleVector; << 107 G4DynamicParticleVector * aResult = new G4DynamicParticleVector;
106 G4ReactionProduct boosted; << 108 G4ReactionProduct boosted;
107 boosted.Lorentz(*(fCache.Get().theNeutronRP) << 109 boosted.Lorentz( *(fCache.Get().theNeutronRP) , *(fCache.Get().theTarget) );
108 G4double eKinetic = boosted.GetKineticEnergy << 110 G4double eKinetic = boosted.GetKineticEnergy();
109 << 111
110 // Build neutrons << 112 // Build neutrons
111 std::vector<G4ReactionProduct> theNeutrons; << 113 G4ReactionProduct * theNeutrons = new G4ReactionProduct[nPrompt+nDelayed];
112 for (i = 0; i < nPrompt + nDelayed; ++i) { << 114 for(i=0; i<nPrompt+nDelayed; i++)
113 theNeutrons.emplace_back(); << 115 {
114 theNeutrons[i].SetDefinition(G4Neutron::Ne << 116 theNeutrons[i].SetDefinition(G4Neutron::Neutron());
115 } << 117 }
116 << 118
117 // sample energies << 119 // sample energies
118 G4int it, dummy; << 120 G4int it, dummy;
119 G4double tempE; << 121 G4double tempE;
120 for (i = 0; i < nPrompt; ++i) { << 122 for(i=0; i<nPrompt; i++)
121 tempE = << 123 {
122 thePromptNeutronEnDis.Sample(eKinetic, d << 124 tempE = thePromptNeutronEnDis.Sample(eKinetic, dummy); // energy distribution (file5) always in lab
123 theNeutrons[i].SetKineticEnergy(tempE); << 125 theNeutrons[i].SetKineticEnergy(tempE);
124 } << 126 }
125 for (i = nPrompt; i < nPrompt + nDelayed; ++ << 127 for(i=nPrompt; i<nPrompt+nDelayed; i++)
126 theNeutrons[i].SetKineticEnergy(theDelayed << 128 {
127 if (it == 0) theNeutrons[i].SetKineticEner << 129 theNeutrons[i].SetKineticEnergy(theDelayedNeutronEnDis.Sample(eKinetic, it)); // dito
128 theDecayConst[i - nPrompt] = theFinalState << 130 if(it==0) theNeutrons[i].SetKineticEnergy(thePromptNeutronEnDis.Sample(eKinetic, dummy));
129 } << 131 theDecayConst[i-nPrompt] = theFinalStateNeutrons.GetDecayConstant(it); // this is returned
130 << 132 }
131 // sample neutron angular distribution << 133
132 for (i = 0; i < nPrompt + nDelayed; ++i) { << 134 // sample neutron angular distribution
133 theNeutronAngularDis.SampleAndUpdate( << 135 for(i=0; i<nPrompt+nDelayed; i++)
134 theNeutrons[i]); // angular comes back << 136 {
135 } << 137 theNeutronAngularDis.SampleAndUpdate(theNeutrons[i]); // angular comes back in lab automatically
136 << 138 }
137 // already in lab. Add neutrons to dynamic p << 139
138 for (i = 0; i < nPrompt + nDelayed; ++i) { << 140 // already in lab. Add neutrons to dynamic particle vector
139 auto dp = new G4DynamicParticle; << 141 for(i=0; i<nPrompt+nDelayed; i++)
140 dp->SetDefinition(theNeutrons[i].GetDefini << 142 {
141 dp->SetMomentum(theNeutrons[i].GetMomentum << 143 G4DynamicParticle * dp = new G4DynamicParticle;
142 aResult->push_back(dp); << 144 dp->SetDefinition(theNeutrons[i].GetDefinition());
>> 145 dp->SetMomentum(theNeutrons[i].GetMomentum());
>> 146 aResult->push_back(dp);
>> 147 }
>> 148 delete [] theNeutrons;
>> 149 // return the result
>> 150 return aResult;
143 } 151 }
144 return aResult; <<
145 } <<
146 152
147 void G4ParticleHPFSFissionFS::SampleNeutronMul << 153 void G4ParticleHPFSFissionFS::SampleNeutronMult(G4int&all, G4int&Prompt, G4int&delayed, G4double eKinetic, G4int off)
148 <<
149 { 154 {
150 G4double promptNeutronMulti = 0; << 155 G4double promptNeutronMulti = 0;
151 promptNeutronMulti = theFinalStateNeutrons.G << 156 promptNeutronMulti = theFinalStateNeutrons.GetPrompt(eKinetic);
152 G4double delayedNeutronMulti = 0; << 157 G4double delayedNeutronMulti = 0;
153 delayedNeutronMulti = theFinalStateNeutrons. << 158 delayedNeutronMulti = theFinalStateNeutrons.GetDelayed(eKinetic);
154 << 159
155 if (delayedNeutronMulti == 0 && promptNeutro << 160 if(delayedNeutronMulti==0&&promptNeutronMulti==0)
156 Prompt = 0; << 161 {
157 delayed = 0; << 162 Prompt = 0;
158 G4double totalNeutronMulti = theFinalState << 163 delayed = 0;
159 all = (G4int)G4Poisson(totalNeutronMulti - << 164 G4double totalNeutronMulti = theFinalStateNeutrons.GetMean(eKinetic);
160 all += off; << 165 all = G4Poisson(totalNeutronMulti-off);
161 } << 166 all += off;
162 else { << 167 }
163 Prompt = (G4int)G4Poisson(promptNeutronMul << 168 else
164 Prompt += off; << 169 {
165 delayed = (G4int)G4Poisson(delayedNeutronM << 170 Prompt = G4Poisson(promptNeutronMulti-off);
166 all = Prompt + delayed; << 171 Prompt += off;
167 } << 172 delayed = G4Poisson(delayedNeutronMulti);
>> 173 all = Prompt+delayed;
>> 174 }
168 } 175 }
169 176
170 G4DynamicParticleVector* G4ParticleHPFSFission << 177 G4DynamicParticleVector * G4ParticleHPFSFissionFS::GetPhotons()
171 { 178 {
172 // sample photons << 179 // sample photons
173 G4ReactionProductVector* temp; << 180 G4ReactionProductVector * temp;
174 G4ReactionProduct boosted; << 181 G4ReactionProduct boosted;
175 << 182 // the photon distributions are in the Nucleus rest frame.
176 // the photon distributions are in the Nucle << 183 boosted.Lorentz( *(fCache.Get().theNeutronRP) , *(fCache.Get().theTarget) );
177 boosted.Lorentz(*(fCache.Get().theNeutronRP) << 184 G4double anEnergy = boosted.GetKineticEnergy();
178 G4double anEnergy = boosted.GetKineticEnergy << 185 temp = theFinalStatePhotons.GetPhotons(anEnergy);
179 temp = theFinalStatePhotons.GetPhotons(anEne << 186 if(temp == 0) { return 0; }
180 if (temp == nullptr) { << 187
181 return nullptr; << 188 // lorentz transform, and add photons to final state
182 } << 189 unsigned int i;
183 << 190 G4DynamicParticleVector * result = new G4DynamicParticleVector;
184 // lorentz transform, and add photons to fin << 191 for(i=0; i<temp->size(); i++)
185 unsigned int i; << 192 {
186 auto result = new G4DynamicParticleVector; << 193 // back to lab
187 for (i = 0; i < temp->size(); ++i) { << 194 temp->operator[](i)->Lorentz(*(temp->operator[](i)), -1.* (*(fCache.Get().theTarget)) );
188 // back to lab << 195 G4DynamicParticle * theOne = new G4DynamicParticle;
189 temp->operator[](i)->Lorentz(*(temp->opera << 196 theOne->SetDefinition(temp->operator[](i)->GetDefinition());
190 auto theOne = new G4DynamicParticle; << 197 theOne->SetMomentum(temp->operator[](i)->GetMomentum());
191 theOne->SetDefinition(temp->operator[](i)- << 198 result->push_back(theOne);
192 theOne->SetMomentum(temp->operator[](i)->G << 199 delete temp->operator[](i);
193 result->push_back(theOne); << 200 }
194 delete temp->operator[](i); << 201 delete temp;
195 } << 202 return result;
196 delete temp; <<
197 return result; <<
198 } 203 }
199 204