| 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 // particle_hp -- source file 26 // particle_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 // 070523 bug fix for G4FPE_DEBUG on by A. How 30 // 070523 bug fix for G4FPE_DEBUG on by A. Howard ( and T. Koi)
31 // 070606 bug fix and migrate to enable to Par << 31 // 070606 bug fix and migrate to enable to Partial cases by T. Koi
32 // 080603 bug fix for Hadron Hyper News #932 b << 32 // 080603 bug fix for Hadron Hyper News #932 by T. Koi
33 // 080612 bug fix contribution from Benoit Pir 33 // 080612 bug fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #4,6
34 // 080717 bug fix of calculation of residual m 34 // 080717 bug fix of calculation of residual momentum by T. Koi
35 // 080801 protect negative available energy by << 35 // 080801 protect negative avalable energy by T. Koi
36 // introduce theNDLDataA,Z which has A 36 // introduce theNDLDataA,Z which has A and Z of NDL data by T. Koi
37 // 081024 G4NucleiPropertiesTable:: to G4Nucle 37 // 081024 G4NucleiPropertiesTable:: to G4NucleiProperties::
38 // 090514 Fix bug in IC electron emission case << 38 // 090514 Fix bug in IC electron emission case
39 // Contribution from Chao Zhang (Chao.Z 39 // Contribution from Chao Zhang (Chao.Zhang@usd.edu) and Dongming Mei(Dongming.Mei@usd.edu)
40 // 100406 "nothingWasKnownOnHadron=1" then sam << 40 // 100406 "nothingWasKnownOnHadron=1" then sample mu isotropic in CM
41 // add two_body_reaction 41 // add two_body_reaction
42 // 100909 add safty << 42 // 100909 add safty
43 // 101111 add safty for _nat_ data case in Bin << 43 // 101111 add safty for _nat_ data case in Binary reaction, but break conservation
44 // 110430 add Reaction Q value and break up fl 44 // 110430 add Reaction Q value and break up flag (MF3::QI and LR)
45 // 45 //
46 // P. Arce, June-2014 Conversion neutron_hp to 46 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
47 // June-2019 - E. Mendoza - re-build "two_body <<
48 // (now isotropic emission in the CMS). Also r <<
49 // developments). Add photon emission when no <<
50 // 47 //
51 // nresp71_m03.hh and nresp71_m02.hh are alike <<
52 // is in the total carbon cross section that i <<
53 // These data are not used in nresp71_m0*.hh. <<
54 // <<
55 <<
56 #include "G4ParticleHPInelasticCompFS.hh" 48 #include "G4ParticleHPInelasticCompFS.hh"
57 << 49 #include "G4ParticleHPManager.hh"
>> 50 #include "G4Nucleus.hh"
>> 51 #include "G4NucleiProperties.hh"
>> 52 #include "G4He3.hh"
58 #include "G4Alpha.hh" 53 #include "G4Alpha.hh"
59 #include "G4Electron.hh" 54 #include "G4Electron.hh"
60 #include "G4He3.hh" <<
61 #include "G4IonTable.hh" <<
62 #include "G4NRESP71M03.hh" <<
63 #include "G4NucleiProperties.hh" <<
64 #include "G4Nucleus.hh" <<
65 #include "G4ParticleHPDataUsed.hh" 55 #include "G4ParticleHPDataUsed.hh"
66 #include "G4ParticleHPManager.hh" << 56 #include "G4IonTable.hh"
67 #include "G4RandomDirection.hh" << 57 #include "G4Pow.hh"
68 #include "G4SystemOfUnits.hh" <<
69 <<
70 #include <fstream> <<
71 <<
72 G4ParticleHPInelasticCompFS::G4ParticleHPInela <<
73 : G4ParticleHPFinalState() <<
74 { <<
75 QI.resize(51); <<
76 LR.resize(51); <<
77 for (G4int i = 0; i < 51; ++i) { <<
78 hasXsec = true; <<
79 theXsection[i] = nullptr; <<
80 theEnergyDistribution[i] = nullptr; <<
81 theAngularDistribution[i] = nullptr; <<
82 theEnergyAngData[i] = nullptr; <<
83 theFinalStatePhotons[i] = nullptr; <<
84 QI[i] = 0.0; <<
85 LR[i] = 0; <<
86 } <<
87 } <<
88 <<
89 G4ParticleHPInelasticCompFS::~G4ParticleHPInel <<
90 { <<
91 for (G4int i = 0; i < 51; ++i) { <<
92 if (theXsection[i] != nullptr) delete theX <<
93 if (theEnergyDistribution[i] != nullptr) d <<
94 if (theAngularDistribution[i] != nullptr) <<
95 if (theEnergyAngData[i] != nullptr) delete <<
96 if (theFinalStatePhotons[i] != nullptr) de <<
97 } <<
98 } <<
99 <<
100 void G4ParticleHPInelasticCompFS::InitDistribu <<
101 G4ReactionPr <<
102 { <<
103 if (theAngularDistribution[it] != nullptr) { <<
104 theAngularDistribution[it]->SetTarget(aTar <<
105 theAngularDistribution[it]->SetProjectileR <<
106 } <<
107 <<
108 if (theEnergyAngData[it] != nullptr) { <<
109 theEnergyAngData[it]->SetTarget(aTarget); <<
110 theEnergyAngData[it]->SetProjectileRP(inPa <<
111 } <<
112 } <<
113 58
114 void G4ParticleHPInelasticCompFS::InitGammas(G 59 void G4ParticleHPInelasticCompFS::InitGammas(G4double AR, G4double ZR)
115 { 60 {
116 G4int Z = G4lrint(ZR); << 61 // char the[100] = {""};
117 G4int A = G4lrint(AR); << 62 // std::ostrstream ost(the, 100, std::ios::out);
118 std::ostringstream ost; << 63 // ost <<gammaPath<<"z"<<ZR<<".a"<<AR;
119 ost << gammaPath << "z" << Z << ".a" << A; << 64 // G4String * aName = new G4String(the);
120 G4String aName = ost.str(); << 65 // std::ifstream from(*aName, std::ios::in);
121 std::ifstream from(aName, std::ios::in); << 66
>> 67 std::ostringstream ost;
>> 68 ost <<gammaPath<<"z"<<ZR<<".a"<<AR;
>> 69 G4String aName = ost.str();
>> 70 std::ifstream from(aName, std::ios::in);
>> 71
>> 72 if(!from) return; // no data found for this isotope
>> 73 // std::ifstream theGammaData(*aName, std::ios::in);
>> 74 std::ifstream theGammaData(aName, std::ios::in);
>> 75
>> 76 theGammas.Init(theGammaData);
>> 77 // delete aName;
122 78
123 if (!from) return; // no data found for thi <<
124 std::ifstream theGammaData(aName, std::ios:: <<
125 <<
126 theGammas.Init(theGammaData); <<
127 } 79 }
128 80
129 void G4ParticleHPInelasticCompFS::Init(G4doubl << 81 void G4ParticleHPInelasticCompFS::Init (G4double A, G4double Z, G4int M, G4String & dirName, G4String & aFSType, G4ParticleDefinition*)
130 const G <<
131 { 82 {
132 gammaPath = fManager->GetNeutronHPPath() + " << 83 gammaPath = "/Inelastic/Gammas/"; //only in neutron data base
133 const G4String& tString = dirName; << 84 if(!getenv("G4NEUTRONHPDATA"))
134 SetA_Z(A, Z, M); << 85 throw G4HadronicException(__FILE__, __LINE__, "Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files where Inelastic/Gammas data is found.");
>> 86 G4String tBase = getenv("G4NEUTRONHPDATA");
>> 87 gammaPath = tBase+gammaPath;
>> 88 G4String tString = dirName;
135 G4bool dbool; 89 G4bool dbool;
136 const G4ParticleHPDataUsed& aFile = << 90 G4ParticleHPDataUsed aFile = theNames.GetName(static_cast<G4int>(A), static_cast<G4int>(Z), M, tString, aFSType, dbool);
137 theNames.GetName(theBaseA, theBaseZ, M, tS << 91 G4String filename = aFile.GetName();
138 SetAZMs(aFile); << 92 #ifdef G4PHPDEBUG
139 const G4String& filename = aFile.GetName(); << 93 if( getenv("G4ParticleHPDebug") ) G4cout << " G4ParticleHPInelasticCompFS::Init FILE " << filename << G4endl;
140 #ifdef G4VERBOSE <<
141 if (fManager->GetDEBUG()) <<
142 G4cout << " G4ParticleHPInelasticCompFS::I <<
143 #endif 94 #endif
144 95
145 SetAZMs(A, Z, M, aFile); << 96 SetAZMs( A, Z, M, aFile );
146 << 97 //theBaseA = aFile.GetA();
147 if (!dbool || (theBaseZ <= 2 && (theNDLDataZ << 98 //theBaseZ = aFile.GetZ();
>> 99 //theNDLDataA = (int)aFile.GetA();
>> 100 //theNDLDataZ = aFile.GetZ();
>> 101 //if(!dbool || ( Z<2.5 && ( std::abs(theBaseZ - Z)>0.0001 || std::abs(theBaseA - A)>0.0001)))
>> 102 if ( !dbool || ( Z<2.5 && ( std::abs(theNDLDataZ - Z)>0.0001 || std::abs(theNDLDataA - A)>0.0001)) )
148 { 103 {
149 #ifdef G4VERBOSE << 104 #ifdef G4PHPDEBUG
150 if (fManager->GetDEBUG()) << 105 if(getenv("G4ParticleHPDebug_NamesLogging")) G4cout << "Skipped = "<< filename <<" "<<A<<" "<<Z<<G4endl;
151 G4cout << "Skipped = " << filename << " <<
152 #endif 106 #endif
153 hasAnyData = false; 107 hasAnyData = false;
154 hasFSData = false; << 108 hasFSData = false;
155 hasXsec = false; 109 hasXsec = false;
156 return; 110 return;
157 } 111 }
158 std::istringstream theData(std::ios::in); << 112 // theBaseA = A;
159 fManager->GetDataStream(filename, theData); << 113 // theBaseZ = G4int(Z+.5);
160 if (!theData) //"!" is a operator of ios << 114 //std::ifstream theData(filename, std::ios::in);
>> 115 std::istringstream theData(std::ios::in);
>> 116 G4ParticleHPManager::GetInstance()->GetDataStream(filename,theData);
>> 117 if(!theData) //"!" is a operator of ios
161 { 118 {
162 hasAnyData = false; 119 hasAnyData = false;
163 hasFSData = false; << 120 hasFSData = false;
164 hasXsec = false; 121 hasXsec = false;
>> 122 // theData.close();
165 return; 123 return;
166 } 124 }
167 // here we go 125 // here we go
168 G4int infoType, dataType, dummy; 126 G4int infoType, dataType, dummy;
169 G4int sfType, it; 127 G4int sfType, it;
170 hasFSData = false; << 128 hasFSData = false;
171 while (theData >> infoType) // Loop checkin << 129 while (theData >> infoType) // Loop checking, 11.05.2015, T. Koi
172 { 130 {
173 hasFSData = true; << 131 hasFSData = true;
174 theData >> dataType; 132 theData >> dataType;
175 theData >> sfType >> dummy; 133 theData >> sfType >> dummy;
176 it = 50; 134 it = 50;
177 if (sfType >= 600 || (sfType < 100 && sfTy << 135 if(sfType>=600||(sfType<100&&sfType>=50)) it = sfType%50;
178 it = sfType % 50; << 136 if(dataType==3)
179 if (dataType == 3) <<
180 { 137 {
>> 138 //theData >> dummy >> dummy;
>> 139 //TK110430
>> 140 // QI and LR introudced since G4NDL3.15
181 G4double dqi; 141 G4double dqi;
182 G4int ilr; 142 G4int ilr;
183 theData >> dqi >> ilr; 143 theData >> dqi >> ilr;
184 144
185 QI[it] = dqi * CLHEP::eV; << 145 QI[ it ] = dqi*CLHEP::eV;
186 LR[it] = ilr; << 146 LR[ it ] = ilr;
187 theXsection[it] = new G4ParticleHPVector 147 theXsection[it] = new G4ParticleHPVector;
188 G4int total; 148 G4int total;
189 theData >> total; 149 theData >> total;
190 theXsection[it]->Init(theData, total, CL 150 theXsection[it]->Init(theData, total, CLHEP::eV);
>> 151 //std::cout << theXsection[it]->GetXsec(1*MeV) << std::endl;
191 } 152 }
192 else if (dataType == 4) { << 153 else if(dataType==4)
>> 154 {
193 theAngularDistribution[it] = new G4Parti 155 theAngularDistribution[it] = new G4ParticleHPAngular;
194 theAngularDistribution[it]->Init(theData 156 theAngularDistribution[it]->Init(theData);
195 } 157 }
196 else if (dataType == 5) { << 158 else if(dataType==5)
>> 159 {
197 theEnergyDistribution[it] = new G4Partic 160 theEnergyDistribution[it] = new G4ParticleHPEnergyDistribution;
198 theEnergyDistribution[it]->Init(theData) << 161 theEnergyDistribution[it]->Init(theData);
199 } 162 }
200 else if (dataType == 6) { << 163 else if(dataType==6)
>> 164 {
201 theEnergyAngData[it] = new G4ParticleHPE 165 theEnergyAngData[it] = new G4ParticleHPEnAngCorrelation(theProjectile);
202 // G4cout << this << " CompFS theEn << 166 // G4cout << this << " CompFS theEnergyAngData " << it << theEnergyAngData[it] << G4endl; //GDEB
203 theEnergyAngData[it]->Init(theData); 167 theEnergyAngData[it]->Init(theData);
204 } 168 }
205 else if (dataType == 12) { << 169 else if(dataType==12)
>> 170 {
206 theFinalStatePhotons[it] = new G4Particl 171 theFinalStatePhotons[it] = new G4ParticleHPPhotonDist;
207 theFinalStatePhotons[it]->InitMean(theDa 172 theFinalStatePhotons[it]->InitMean(theData);
208 } 173 }
209 else if (dataType == 13) { << 174 else if(dataType==13)
>> 175 {
210 theFinalStatePhotons[it] = new G4Particl 176 theFinalStatePhotons[it] = new G4ParticleHPPhotonDist;
211 theFinalStatePhotons[it]->InitPartials(t << 177 theFinalStatePhotons[it]->InitPartials(theData);
212 } 178 }
213 else if (dataType == 14) { << 179 else if(dataType==14)
>> 180 {
214 theFinalStatePhotons[it]->InitAngular(th 181 theFinalStatePhotons[it]->InitAngular(theData);
215 } 182 }
216 else if (dataType == 15) { << 183 else if(dataType==15)
>> 184 {
217 theFinalStatePhotons[it]->InitEnergies(t 185 theFinalStatePhotons[it]->InitEnergies(theData);
218 } 186 }
219 else { << 187 else
220 G4ExceptionDescription ed; << 188 {
221 ed << "Z=" << theBaseZ << " A=" << theBa << 189 throw G4HadronicException(__FILE__, __LINE__, "Data-type unknown to G4ParticleHPInelasticCompFS");
222 << " projectile: " << theProjectile->GetPar <<
223 G4Exception("G4ParticleHPInelasticCompFS <<
224 ed, "Data-type unknown"); <<
225 } 190 }
226 } 191 }
>> 192 // theData.close();
227 } 193 }
228 194
229 G4int G4ParticleHPInelasticCompFS::SelectExitC 195 G4int G4ParticleHPInelasticCompFS::SelectExitChannel(G4double eKinetic)
230 { 196 {
>> 197
>> 198 // it = 0 has without Photon
231 G4double running[50]; 199 G4double running[50];
232 running[0] = 0; 200 running[0] = 0;
233 G4int i; << 201 unsigned int i;
234 for (i = 0; i < 50; ++i) { << 202 for(i=0; i<50; i++)
235 if (i != 0) running[i] = running[i - 1]; << 203 {
236 if (theXsection[i] != nullptr) { << 204 if(i!=0) running[i]=running[i-1];
>> 205 if(theXsection[i] != 0)
>> 206 {
237 running[i] += std::max(0., theXsection[i 207 running[i] += std::max(0., theXsection[i]->GetXsec(eKinetic));
238 } 208 }
239 } 209 }
240 G4double random = G4UniformRand(); 210 G4double random = G4UniformRand();
241 G4double sum = running[49]; 211 G4double sum = running[49];
242 G4int it = 50; 212 G4int it = 50;
243 if (0 != sum) { << 213 if(0!=sum)
>> 214 {
244 G4int i0; 215 G4int i0;
245 for (i0 = 0; i0 < 50; ++i0) { << 216 for(i0=0; i0<50; i0++)
>> 217 {
246 it = i0; 218 it = i0;
247 if (random < running[i0] / sum) break; << 219 // G4cout << " SelectExitChannel " << it << " " << random << " " << running[i0]/sum << " " << running[i0] << G4endl; //GDEB
>> 220 if(random < running[i0]/sum) break;
248 } 221 }
249 } 222 }
>> 223 //debug: it = 1;
>> 224 // G4cout << " SelectExitChannel " << it << " " << sum << G4endl; //GDEB
250 return it; 225 return it;
251 } 226 }
252 227
253 // n,p,d,t,he3,a << 228
254 void G4ParticleHPInelasticCompFS::CompositeApp << 229 //n,p,d,t,he3,a
255 << 230 void G4ParticleHPInelasticCompFS::CompositeApply(const G4HadProjectile & theTrack, G4ParticleDefinition * aDefinition)
256 { 231 {
257 // prepare neutron <<
258 if (theResult.Get() == nullptr) theResult.Pu <<
259 theResult.Get()->Clear(); <<
260 G4double eKinetic = theTrack.GetKineticEnerg <<
261 const G4HadProjectile* hadProjectile = &theT <<
262 G4ReactionProduct incidReactionProduct(hadPr <<
263 incidReactionProduct.SetMomentum(hadProjecti <<
264 incidReactionProduct.SetKineticEnergy(eKinet <<
265 <<
266 // prepare target <<
267 for (G4int i = 0; i < 50; ++i) { <<
268 if (theXsection[i] != nullptr) { <<
269 break; <<
270 } <<
271 } <<
272 232
273 G4double targetMass = G4NucleiProperties::Ge << 233 // prepare neutron
274 #ifdef G4VERBOSE << 234 if ( theResult.Get() == NULL ) theResult.Put( new G4HadFinalState );
275 if (fManager->GetDEBUG()) << 235 theResult.Get()->Clear();
276 G4cout << "G4ParticleHPInelasticCompFS::Co << 236 G4double eKinetic = theTrack.GetKineticEnergy();
277 << theBaseZ << " incident " << hadProject << 237 const G4HadProjectile *hadProjectile = &theTrack;
278 << G4endl; << 238 G4ReactionProduct incidReactionProduct( const_cast<G4ParticleDefinition *>(hadProjectile->GetDefinition()) ); // incidReactionProduct
>> 239 incidReactionProduct.SetMomentum( hadProjectile->Get4Momentum().vect() );
>> 240 incidReactionProduct.SetKineticEnergy( eKinetic );
>> 241
>> 242 // prepare target
>> 243 G4int i;
>> 244 for(i=0; i<50; i++)
>> 245 { if(theXsection[i] != 0) { break; } }
>> 246
>> 247 G4double targetMass=0;
>> 248 G4double eps = 0.0001;
>> 249 targetMass = ( G4NucleiProperties::GetNuclearMass(static_cast<G4int>(theBaseA+eps), static_cast<G4int>(theBaseZ+eps))) /
>> 250 theProjectile->GetPDGMass();
>> 251 #ifdef G4PHPDEBUG
>> 252 if( getenv("G4ParticleHPDebug")) G4cout <<this <<" G4ParticleHPInelasticCompFS::CompositeApply A " <<theBaseA <<" Z " <<theBaseZ <<" incident " <<hadProjectile->GetDefinition()->GetParticleName() <<G4endl;
279 #endif 253 #endif
280 G4ReactionProduct theTarget; << 254 // if(theEnergyAngData[i]!=0)
281 G4Nucleus aNucleus; << 255 // targetMass = theEnergyAngData[i]->GetTargetMass();
282 // G4ThreeVector neuVelo = << 256 // else if(theAngularDistribution[i]!=0)
283 // (1./hadProjectile->GetDefinition()->GetPD << 257 // targetMass = theAngularDistribution[i]->GetTargetMass();
284 // = aNucleus.GetBiasedThermalNucleus( targe << 258 // else if(theFinalStatePhotons[50]!=0)
285 // neuVelo, theTrack.GetMaterial()->GetTempe << 259 // targetMass = theFinalStatePhotons[50]->GetTargetMass();
286 // normalization of mass and velocity in neu << 260 G4ReactionProduct theTarget;
287 G4ThreeVector neuVelo = incidReactionProduct << 261 G4Nucleus aNucleus;
288 theTarget = aNucleus.GetBiasedThermalNucleus << 262 G4ThreeVector neuVelo = (1./hadProjectile->GetDefinition()->GetPDGMass())*incidReactionProduct.GetMomentum();
289 << 263 theTarget = aNucleus.GetBiasedThermalNucleus( targetMass, neuVelo, theTrack.GetMaterial()->GetTemperature());
290 << 264 theTarget.SetDefinition( G4IonTable::GetIonTable()->GetIon( G4int(theBaseZ), G4int(theBaseA) , 0.0 ) ); //XX
291 theTarget.SetDefinition(G4IonTable::GetIonTa << 265
292 << 266 // prepare the residual mass
293 // prepare the residual mass << 267 G4double residualMass=0;
294 G4double residualMass = 0; << 268 G4double residualZ = theBaseZ - aDefinition->GetPDGCharge();
295 G4int residualZ = theBaseZ + << 269 G4double residualA = theBaseA - aDefinition->GetBaryonNumber()+1;
296 G4lrint((theProjectile->GetPDGCharge() - a << 270 residualMass = ( G4NucleiProperties::GetNuclearMass(static_cast<G4int>(residualA+eps), static_cast<G4int>(residualZ+eps)) ) /
297 G4int residualA = theBaseA + theProjectile-> << 271 theProjectile->GetPDGMass();
298 residualMass = G4NucleiProperties::GetNuclea << 272
299 << 273 // prepare energy in target rest frame
300 // prepare energy in target rest frame << 274 G4ReactionProduct boosted;
301 G4ReactionProduct boosted; << 275 boosted.Lorentz(incidReactionProduct, theTarget);
302 boosted.Lorentz(incidReactionProduct, theTar << 276 eKinetic = boosted.GetKineticEnergy();
303 eKinetic = boosted.GetKineticEnergy(); << 277 // G4double momentumInCMS = boosted.GetTotalMomentum();
304 << 278
305 // select exit channel for composite FS clas << 279 // select exit channel for composite FS class.
306 G4int it = SelectExitChannel(eKinetic); << 280 G4int it = SelectExitChannel( eKinetic );
307 << 281
308 // E. Mendoza (2018) -- to use JENDL/AN-2005 << 282 // set target and neutron in the relevant exit channel
309 if (theEnergyDistribution[it] == nullptr && << 283 InitDistributionInitialState(incidReactionProduct, theTarget, it);
310 && theEnergyAngData[it] == nullptr) << 284
311 { << 285 G4ReactionProductVector * thePhotons = 0;
312 if (theEnergyDistribution[50] != nullptr | << 286 G4ReactionProductVector * theParticles = 0;
313 || theEnergyAngData[50] != nullptr) << 287 G4ReactionProduct aHadron;
>> 288 aHadron.SetDefinition(aDefinition); // what if only cross-sections exist ==> Na 23 11 @@@@
>> 289 G4double availableEnergy = incidReactionProduct.GetKineticEnergy() + incidReactionProduct.GetMass() - aHadron.GetMass() +
>> 290 (targetMass - residualMass)*theProjectile->GetPDGMass();
>> 291 //080730c
>> 292 if ( availableEnergy < 0 )
314 { 293 {
315 it = 50; << 294 //G4cout << "080730c Adjust availavleEnergy " << G4endl;
>> 295 availableEnergy = 0;
316 } 296 }
317 } << 297 G4int nothingWasKnownOnHadron = 0;
>> 298 G4int dummy;
>> 299 G4double eGamm = 0;
>> 300 G4int iLevel=it-1;
>> 301
>> 302 // TK without photon has it = 0
>> 303 if( 50 == it )
>> 304 {
>> 305
>> 306 // TK Excitation level is not determined
>> 307 iLevel=-1;
>> 308 aHadron.SetKineticEnergy(availableEnergy*residualMass*theProjectile->GetPDGMass()/
>> 309 (aHadron.GetMass()+residualMass*theProjectile->GetPDGMass()));
>> 310
>> 311 //aHadron.SetMomentum(incidReactionProduct.GetMomentum()*(1./incidReactionProduct.GetTotalMomentum())*
>> 312 // std::sqrt(aHadron.GetTotalEnergy()*aHadron.GetTotalEnergy()-
>> 313 // aHadron.GetMass()*aHadron.GetMass()));
318 314
319 // set target and neutron in the relevant ex << 315 //TK add safty 100909
320 InitDistributionInitialState(incidReactionPr << 316 G4double p2 = ( aHadron.GetTotalEnergy()*aHadron.GetTotalEnergy() - aHadron.GetMass()*aHadron.GetMass() );
>> 317 G4double p = 0.0;
>> 318 if ( p2 > 0.0 ) p = std::sqrt( p );
321 319
322 //------------------------------------------ << 320 aHadron.SetMomentum(incidReactionProduct.GetMomentum()*(1./incidReactionProduct.GetTotalMomentum())*p );
323 // Hook for NRESP71MODEL << 321
324 if (fManager->GetUseNRESP71Model() && eKinet << 322 }
325 if (theBaseZ == 6) // If the reaction is << 323 else
326 { 324 {
327 if (theProjectile == G4Neutron::Definiti << 325 while ( iLevel!=-1 && theGammas.GetLevel(iLevel) == 0 ) { iLevel--; } // Loop checking, 11.05.2015, T. Koi
328 if (use_nresp71_model(aDefinition, it, <<
329 } <<
330 } 326 }
331 } <<
332 //------------------------------------------ <<
333 327
334 G4ReactionProductVector* thePhotons = nullpt <<
335 G4ReactionProductVector* theParticles = null <<
336 G4ReactionProduct aHadron; <<
337 aHadron.SetDefinition(aDefinition); // what <<
338 G4double availableEnergy = incidReactionProd <<
339 + incidReactionPr <<
340 + (targetMass - r <<
341 328
342 if (availableEnergy < 0) { << 329 if ( theAngularDistribution[it] != 0 ) // MF4
343 availableEnergy = 0; <<
344 } <<
345 G4int nothingWasKnownOnHadron = 0; <<
346 G4double eGamm = 0; <<
347 G4int iLevel = -1; <<
348 // max gamma energy and index <<
349 G4int imaxEx = theGammas.GetNumberOfLevels() <<
350 <<
351 // without photon has it = 0 <<
352 if (50 == it) { <<
353 // Excitation level is not determined <<
354 aHadron.SetKineticEnergy(availableEnergy * <<
355 <<
356 // TK add safty 100909 <<
357 G4double p2 = <<
358 (aHadron.GetTotalEnergy() * aHadron.GetT <<
359 G4double p = (p2 > 0.0) ? std::sqrt(p2) : <<
360 aHadron.SetMomentum(p * incidReactionProdu <<
361 incidReactionProduct.GetTotalMomentum()) <<
362 } <<
363 else { <<
364 iLevel = imaxEx; <<
365 } <<
366 <<
367 if (theAngularDistribution[it] != nullptr) <<
368 { <<
369 if (theEnergyDistribution[it] != nullptr) <<
370 { 330 {
371 //************************************** << 331 if(theEnergyDistribution[it]!=0) // MF5
372 /* << 332 {
373 aHadron.SetKineticEnergy(theEnergy << 333 //************************************************************
374 G4double eSecN = aHadron.GetKineti << 334 /*
375 */ << 335 aHadron.SetKineticEnergy(theEnergyDistribution[it]->Sample(eKinetic, dummy));
376 //************************************** << 336 G4double eSecN = aHadron.GetKineticEnergy();
377 // EMendoza --> maximum allowable energy << 337 */
378 G4double dqi = 0.0; << 338 //************************************************************
379 if (QI[it] < 0 || 849 < QI[it]) << 339 //EMendoza --> maximum allowable energy should be taken into account.
380 dqi = QI[it]; // For backword compati << 340 G4double dqi = 0.0;
381 G4double MaxEne = eKinetic + dqi; << 341 if ( QI[it] < 0 || 849 < QI[it] ) dqi = QI[it]; //For backword compatibility QI introduced since G4NDL3.15
382 G4double eSecN = 0.; << 342 G4double MaxEne=eKinetic+dqi;
383 << 343 G4double eSecN;
384 G4int icounter = 0; << 344
385 G4int icounter_max = 1024; << 345 G4int icounter=0;
386 G4int dummy = 0; << 346 G4int icounter_max=1024;
387 do { << 347 do {
388 ++icounter; << 348 icounter++;
389 if (icounter > icounter_max) { << 349 if ( icounter > icounter_max ) {
390 G4cout << "Loop-counter exceeded the << 350 G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl;
391 << __FILE__ << "." << G4endl; << 351 break;
392 break; << 352 }
>> 353 eSecN=theEnergyDistribution[it]->Sample(eKinetic, dummy);
>> 354 }while(eSecN>MaxEne); // Loop checking, 11.05.2015, T. Koi
>> 355 aHadron.SetKineticEnergy(eSecN);
>> 356 //************************************************************
>> 357 eGamm = eKinetic-eSecN;
>> 358 for(iLevel=theGammas.GetNumberOfLevels()-1; iLevel>=0; iLevel--)
>> 359 {
>> 360 if(theGammas.GetLevelEnergy(iLevel)<eGamm) break;
393 } 361 }
394 eSecN = theEnergyDistribution[it]->Sam << 362 G4double random = 2*G4UniformRand();
395 } while (eSecN > MaxEne); // Loop check << 363 iLevel+=G4int(random);
396 aHadron.SetKineticEnergy(eSecN); << 364 if(iLevel>theGammas.GetNumberOfLevels()-1)iLevel = theGammas.GetNumberOfLevels()-1;
397 //************************************** <<
398 eGamm = eKinetic - eSecN; <<
399 for (iLevel = imaxEx; iLevel >= 0; --iLe <<
400 if (theGammas.GetLevelEnergy(iLevel) < <<
401 } 365 }
402 if (iLevel < imaxEx && iLevel >= 0) { << 366 else
403 if (G4UniformRand() > 0.5) { << 367 {
404 ++iLevel; << 368 G4double eExcitation = 0;
>> 369 if(iLevel>=0) eExcitation = theGammas.GetLevel(iLevel)->GetLevelEnergy();
>> 370 while (eKinetic-eExcitation < 0 && iLevel>0) // Loop checking, 11.05.2015, T. Koi
>> 371 {
>> 372 iLevel--;
>> 373 eExcitation = theGammas.GetLevel(iLevel)->GetLevelEnergy();
>> 374 }
>> 375 //110610TK BEGIN
>> 376 //Use QI value for calculating excitation energy of residual.
>> 377 G4bool useQI=false;
>> 378 G4double dqi = QI[it];
>> 379 if ( dqi < 0 || 849 < dqi ) useQI = true; //Former libraies does not have values of this range
>> 380
>> 381 if ( useQI )
>> 382 {
>> 383 // QI introudced since G4NDL3.15
>> 384 eExcitation = -QI[it];
>> 385 //Re-evluate iLevel based on this eExcitation
>> 386 iLevel = 0;
>> 387 G4bool find = false;
>> 388 G4int imaxEx = 0;
>> 389 while( theGammas.GetLevel(iLevel+1) != 0 ) // Loop checking, 11.05.2015, T. Koi
>> 390 {
>> 391 G4double maxEx = 0.0;
>> 392 if ( maxEx < theGammas.GetLevel(iLevel)->GetLevelEnergy() )
>> 393 {
>> 394 maxEx = theGammas.GetLevel(iLevel)->GetLevelEnergy();
>> 395 imaxEx = iLevel;
>> 396 }
>> 397 if ( eExcitation < theGammas.GetLevel(iLevel)->GetLevelEnergy() )
>> 398 {
>> 399 find = true;
>> 400 iLevel--;
>> 401 // very small eExcitation, iLevel becomes -1, this is protected below.
>> 402 if ( iLevel == -1 ) iLevel = 0; // But cause energy trouble.
>> 403 break;
>> 404 }
>> 405 iLevel++;
>> 406 }
>> 407 // In case, cannot find proper level, then use the maximum level.
>> 408 if ( !find ) iLevel = imaxEx;
405 } 409 }
>> 410 //110610TK END
>> 411
>> 412 if(getenv("G4ParticleHPDebug") && eKinetic-eExcitation < 0)
>> 413 {
>> 414 throw G4HadronicException(__FILE__, __LINE__, "SEVERE: InelasticCompFS: Consistency of data not good enough, please file report");
>> 415 }
>> 416 if(eKinetic-eExcitation < 0) eExcitation = 0;
>> 417 if(iLevel!= -1) aHadron.SetKineticEnergy(eKinetic - eExcitation);
>> 418
406 } 419 }
407 } << 420 theAngularDistribution[it]->SampleAndUpdate(aHadron);
408 else { << 421
409 G4double eExcitation = 0; << 422 if( theFinalStatePhotons[it] == 0 )
410 for (iLevel = imaxEx; iLevel >= 0; --iLe << 423 {
411 if (theGammas.GetLevelEnergy(iLevel) < << 424 //G4cout << "110610 USE Gamma Level" << G4endl;
>> 425 // TK comment Most n,n* eneter to this
>> 426 thePhotons = theGammas.GetDecayGammas(iLevel);
>> 427 eGamm -= theGammas.GetLevelEnergy(iLevel);
>> 428 if(eGamm>0) // @ ok for now, but really needs an efficient way of correllated sampling @
>> 429 {
>> 430 G4ReactionProduct * theRestEnergy = new G4ReactionProduct;
>> 431 theRestEnergy->SetDefinition(G4Gamma::Gamma());
>> 432 theRestEnergy->SetKineticEnergy(eGamm);
>> 433 G4double costh = 2.*G4UniformRand()-1.;
>> 434 G4double phi = CLHEP::twopi*G4UniformRand();
>> 435 theRestEnergy->SetMomentum(eGamm*std::sin(std::acos(costh))*std::cos(phi),
>> 436 eGamm*std::sin(std::acos(costh))*std::sin(phi),
>> 437 eGamm*costh);
>> 438 if(thePhotons == 0) { thePhotons = new G4ReactionProductVector; }
>> 439 thePhotons->push_back(theRestEnergy);
>> 440 }
412 } 441 }
>> 442 }
>> 443 else if(theEnergyAngData[it] != 0) // MF6
>> 444 {
413 445
414 // Use QI value for calculating excitati << 446 theParticles = theEnergyAngData[it]->Sample(eKinetic);
415 G4bool useQI = false; <<
416 G4double dqi = QI[it]; <<
417 if (dqi < 0 || 849 < dqi) useQI = true; <<
418 <<
419 if (useQI) { <<
420 eExcitation = std::max(0., QI[0] - QI[ <<
421 <<
422 // Re-evaluate iLevel based on this eE <<
423 iLevel = 0; <<
424 G4bool find = false; <<
425 const G4double level_tolerance = 1.0 * <<
426 <<
427 // VI: the first level is ground <<
428 if (0 < imaxEx) { <<
429 for (iLevel = 1; iLevel <= imaxEx; + <<
430 G4double elevel = theGammas.GetLev <<
431 if (std::abs(eExcitation - elevel) <<
432 find = true; <<
433 break; <<
434 } <<
435 if (eExcitation < elevel) { <<
436 find = true; <<
437 iLevel = std::max(iLevel - 1, 0) <<
438 break; <<
439 } <<
440 } <<
441 447
442 // If proper level cannot be found, << 448 //141017 Fix BEGIN
443 if (!find) iLevel = imaxEx; << 449 //Adjust A and Z in the case of miss much between selected data and target nucleus
444 } << 450 if ( theParticles != NULL ) {
>> 451 G4int sumA = 0;
>> 452 G4int sumZ = 0;
>> 453 G4int maxA = 0;
>> 454 G4int jAtMaxA = 0;
>> 455 for ( G4int j = 0 ; j != (G4int)theParticles->size() ; j++ ) {
>> 456 if ( theParticles->at(j)->GetDefinition()->GetBaryonNumber() > maxA ) {
>> 457 maxA = theParticles->at(j)->GetDefinition()->GetBaryonNumber();
>> 458 jAtMaxA = j;
>> 459 }
>> 460 sumA += theParticles->at(j)->GetDefinition()->GetBaryonNumber();
>> 461 sumZ += G4int( theParticles->at(j)->GetDefinition()->GetPDGCharge() + eps );
>> 462 }
>> 463 G4int dA = (G4int)theBaseA + hadProjectile->GetDefinition()->GetBaryonNumber() - sumA;
>> 464 G4int dZ = (G4int)theBaseZ + G4int( hadProjectile->GetDefinition()->GetPDGCharge() + eps ) - sumZ;
>> 465 if ( dA < 0 || dZ < 0 ) {
>> 466 G4int newA = theParticles->at(jAtMaxA)->GetDefinition()->GetBaryonNumber() + dA ;
>> 467 G4int newZ = G4int( theParticles->at(jAtMaxA)->GetDefinition()->GetPDGCharge() + eps ) + dZ;
>> 468 G4ParticleDefinition* pd = G4IonTable::GetIonTable()->GetIon ( newZ , newA );
>> 469 theParticles->at( jAtMaxA )->SetDefinition( pd );
>> 470 }
445 } 471 }
>> 472 //141017 Fix END
446 473
447 if (fManager->GetDEBUG() && eKinetic - e <<
448 throw G4HadronicException( <<
449 __FILE__, __LINE__, <<
450 "SEVERE: InelasticCompFS: Consistenc <<
451 } <<
452 if (eKinetic - eExcitation < 0) eExcitat <<
453 if (iLevel != -1) aHadron.SetKineticEner <<
454 } 474 }
455 theAngularDistribution[it]->SampleAndUpdat << 475 else
456 << 476 {
457 if (theFinalStatePhotons[it] == nullptr) { << 477 // @@@ what to do, if we have photon data, but no info on the hadron itself
458 thePhotons = theGammas.GetDecayGammas(iL << 478 nothingWasKnownOnHadron = 1;
459 eGamm -= theGammas.GetLevelEnergy(iLevel <<
460 } 479 }
461 } <<
462 else if (theEnergyAngData[it] != nullptr) / <<
463 { <<
464 theParticles = theEnergyAngData[it]->Sampl <<
465 480
466 // Adjust A and Z in the case of miss much << 481 //G4cout << "theFinalStatePhotons it " << it << G4endl;
467 if (theParticles != nullptr) { << 482 //G4cout << "theFinalStatePhotons[it] " << theFinalStatePhotons[it] << G4endl;
468 G4int sumA = 0; << 483 //G4cout << "theFinalStatePhotons it " << it << G4endl;
469 G4int sumZ = 0; << 484 //G4cout << "theFinalStatePhotons[it] " << theFinalStatePhotons[it] << G4endl;
470 G4int maxA = 0; << 485 //G4cout << "thePhotons " << thePhotons << G4endl;
471 G4int jAtMaxA = 0; <<
472 for (G4int j = 0; j != (G4int)theParticl <<
473 auto ptr = theParticles->at(j); <<
474 G4int barnum = ptr->GetDefinition()->GetBary <<
475 if (barnum > maxA) { <<
476 maxA = barnum; <<
477 jAtMaxA = j; <<
478 } <<
479 sumA += barnum; <<
480 sumZ += G4lrint(ptr->GetDefinition()-> <<
481 } <<
482 G4int dA = theBaseA + hadProjectile->Get <<
483 G4int dZ = theBaseZ + <<
484 G4lrint(hadProjectile->GetDefinition()->GetP <<
485 if (dA < 0 || dZ < 0) { <<
486 G4int newA = theParticles->at(jAtMaxA) <<
487 G4int newZ = <<
488 G4lrint(theParticles->at(jAtMaxA)->GetDefi <<
489 G4ParticleDefinition* pd = ionTable->G <<
490 theParticles->at(jAtMaxA)->SetDefiniti <<
491 } <<
492 } <<
493 } <<
494 else { <<
495 // @@@ what to do, if we have photon data, <<
496 nothingWasKnownOnHadron = 1; <<
497 } <<
498 486
499 if (theFinalStatePhotons[it] != nullptr) { << 487 if ( theFinalStatePhotons[it] != 0 )
500 // the photon distributions are in the Nuc << 488 {
501 // TK residual rest frame << 489 // the photon distributions are in the Nucleus rest frame.
502 G4ReactionProduct boosted_tmp; << 490 // TK residual rest frame
503 boosted_tmp.Lorentz(incidReactionProduct, << 491 G4ReactionProduct boosted_tmp;
504 G4double anEnergy = boosted_tmp.GetKinetic << 492 boosted_tmp.Lorentz(incidReactionProduct, theTarget);
505 thePhotons = theFinalStatePhotons[it]->Get << 493 G4double anEnergy = boosted_tmp.GetKineticEnergy();
506 G4double aBaseEnergy = theFinalStatePhoton << 494 thePhotons = theFinalStatePhotons[it]->GetPhotons(anEnergy);
507 G4double testEnergy = 0; << 495 G4double aBaseEnergy = theFinalStatePhotons[it]->GetLevelEnergy();
508 if (thePhotons != nullptr && !thePhotons-> << 496 G4double testEnergy = 0;
509 aBaseEnergy -= (*thePhotons)[0]->GetTota << 497 if(thePhotons!=0 && thePhotons->size()!=0)
510 } << 498 { aBaseEnergy-=thePhotons->operator[](0)->GetTotalEnergy(); }
511 if (theFinalStatePhotons[it]->NeedsCascade << 499 if(theFinalStatePhotons[it]->NeedsCascade())
512 while (aBaseEnergy > 0.01 * CLHEP::keV) << 500 {
513 { << 501 while(aBaseEnergy>0.01*CLHEP::keV) // Loop checking, 11.05.2015, T. Koi
514 // cascade down the levels << 502 {
515 G4bool foundMatchingLevel = false; << 503 // cascade down the levels
516 G4int closest = 2; << 504 G4bool foundMatchingLevel = false;
517 G4double deltaEold = -1; << 505 G4int closest = 2;
518 for (G4int j = 1; j < it; ++j) { << 506 G4double deltaEold = -1;
519 if (theFinalStatePhotons[j] != nullp << 507 for(G4int j=1; j<it; j++)
520 testEnergy = theFinalStatePhotons[ << 508 {
521 } << 509 if(theFinalStatePhotons[j]!=0)
522 else { << 510 {
523 testEnergy = 0; << 511 testEnergy = theFinalStatePhotons[j]->GetLevelEnergy();
524 } << 512 }
525 G4double deltaE = std::abs(testEnerg << 513 else
526 if (deltaE < 0.1 * CLHEP::keV) { << 514 {
527 G4ReactionProductVector* theNext = << 515 testEnergy = 0;
528 if (thePhotons != nullptr) thePhot << 516 }
529 aBaseEnergy = testEnergy - theNext << 517 G4double deltaE = std::abs(testEnergy-aBaseEnergy);
>> 518 if(deltaE<0.1*CLHEP::keV)
>> 519 {
>> 520 G4ReactionProductVector * theNext =
>> 521 theFinalStatePhotons[j]->GetPhotons(anEnergy);
>> 522 thePhotons->push_back(theNext->operator[](0));
>> 523 aBaseEnergy = testEnergy-theNext->operator[](0)->GetTotalEnergy();
>> 524 delete theNext;
>> 525 foundMatchingLevel = true;
>> 526 break; // ===>
>> 527 }
>> 528 if(theFinalStatePhotons[j]!=0 && ( deltaE<deltaEold||deltaEold<0.) )
>> 529 {
>> 530 closest = j;
>> 531 deltaEold = deltaE;
>> 532 }
>> 533 } // <=== the break goes here.
>> 534 if(!foundMatchingLevel)
>> 535 {
>> 536 G4ReactionProductVector * theNext =
>> 537 theFinalStatePhotons[closest]->GetPhotons(anEnergy);
>> 538 thePhotons->push_back(theNext->operator[](0));
>> 539 aBaseEnergy = aBaseEnergy-theNext->operator[](0)->GetTotalEnergy();
530 delete theNext; 540 delete theNext;
531 foundMatchingLevel = true; << 541 }
532 break; // ===> << 542 }
533 } <<
534 if (theFinalStatePhotons[j] != nullp <<
535 closest = j; <<
536 deltaEold = deltaE; <<
537 } <<
538 } // <=== the break goes here. <<
539 if (!foundMatchingLevel) { <<
540 G4ReactionProductVector* theNext = t <<
541 if (thePhotons != nullptr) thePhoton <<
542 aBaseEnergy = aBaseEnergy - theNext- <<
543 delete theNext; <<
544 } <<
545 } 543 }
546 } 544 }
547 } << 545 unsigned int i0;
548 << 546 if(thePhotons!=0)
549 if (thePhotons != nullptr) { << 547 {
550 for (auto const & p : *thePhotons) { << 548 for(i0=0; i0<thePhotons->size(); i0++)
551 // back to lab << 549 {
552 p->Lorentz(*p, -1. * theTarget); << 550 // back to lab
553 } << 551 thePhotons->operator[](i0)->Lorentz(*(thePhotons->operator[](i0)), -1.*theTarget);
554 } <<
555 if (nothingWasKnownOnHadron != 0) { <<
556 // In this case, hadron should be isotropi <<
557 // Next 12 lines are Emilio's replacement <<
558 // G4double QM=(incidReactionProduct.GetMa <<
559 // G4double eExcitation = QM-QI[it]; <<
560 // G4double eExcitation = QI[0] - QI[it]; <<
561 // if(eExcitation<20*CLHEP::keV){eExcitati <<
562 <<
563 G4double eExcitation = std::max(0., QI[0] <<
564 <<
565 two_body_reaction(&incidReactionProduct, & <<
566 if (thePhotons == nullptr && eExcitation > <<
567 for (iLevel = imaxEx; iLevel >= 0; --iLe <<
568 if (theGammas.GetLevelEnergy(iLevel) < <<
569 } 552 }
570 thePhotons = theGammas.GetDecayGammas(iL <<
571 } 553 }
572 } << 554 //G4cout << "nothingWasKnownOnHadron " << nothingWasKnownOnHadron << G4endl;
>> 555 if(nothingWasKnownOnHadron)
>> 556 {
>> 557 // TKDB 100405
>> 558 // In this case, hadron should be isotropic in CM
>> 559 // mu and p should be correlated
>> 560 //
>> 561 G4double totalPhotonEnergy = 0.0;
>> 562 if ( thePhotons != 0 )
>> 563 {
>> 564 unsigned int nPhotons = thePhotons->size();
>> 565 unsigned int ii0;
>> 566 for ( ii0=0; ii0<nPhotons; ii0++)
>> 567 {
>> 568 //thePhotons has energies at LAB system
>> 569 totalPhotonEnergy += thePhotons->operator[](ii0)->GetTotalEnergy();
>> 570 }
>> 571 }
>> 572
>> 573 //isotropic distribution in CM
>> 574 G4double mu = 1.0 - 2 * G4UniformRand();
>> 575
>> 576 // need momentums in target rest frame;
>> 577 G4LorentzVector target_in_LAB ( theTarget.GetMomentum() , theTarget.GetTotalEnergy() );
>> 578 G4ThreeVector boostToTargetRest = -target_in_LAB.boostVector();
>> 579 G4LorentzVector proj_in_LAB = hadProjectile->Get4Momentum();
>> 580
>> 581 G4DynamicParticle* proj = new G4DynamicParticle( theProjectile , proj_in_LAB.boost( boostToTargetRest ) );
>> 582 G4DynamicParticle* targ = new G4DynamicParticle( G4IonTable::GetIonTable()->GetIon ( (G4int)theBaseZ , (G4int)theBaseA , totalPhotonEnergy ) , G4ThreeVector(0) );
>> 583 G4DynamicParticle* hadron = new G4DynamicParticle( aHadron.GetDefinition() , G4ThreeVector(0) ); // will be fill momentum
>> 584
>> 585 two_body_reaction ( proj , targ , hadron , mu );
>> 586
>> 587 G4LorentzVector hadron_in_trag_rest = hadron->Get4Momentum();
>> 588 G4LorentzVector hadron_in_LAB = hadron_in_trag_rest.boost ( -boostToTargetRest );
>> 589 aHadron.SetMomentum( hadron_in_LAB.v() );
>> 590 aHadron.SetKineticEnergy ( hadron_in_LAB.e() - hadron_in_LAB.m() );
>> 591
>> 592 delete proj;
>> 593 delete targ;
>> 594 delete hadron;
>> 595
>> 596 //TKDB 100405
>> 597 /*
>> 598 G4double totalPhotonEnergy = 0;
>> 599 if(thePhotons!=0)
>> 600 {
>> 601 unsigned int nPhotons = thePhotons->size();
>> 602 unsigned int i0;
>> 603 for(i0=0; i0<nPhotons; i0++)
>> 604 {
>> 605 totalPhotonEnergy += thePhotons->operator[](i0)->GetTotalEnergy();
>> 606 }
>> 607 }
>> 608 availableEnergy -= totalPhotonEnergy;
>> 609 residualMass += totalPhotonEnergy/theProjectile->GetPDGMass();
>> 610 aHadron.SetKineticEnergy(availableEnergy*residualMass*theProjectile->GetPDGMass()/
>> 611 (aHadron.GetMass()+residualMass*theProjectile->GetPDGMass()));
>> 612 G4double CosTheta = 1.0 - 2.0*G4UniformRand();
>> 613 G4double SinTheta = std::sqrt(1.0 - CosTheta*CosTheta);
>> 614 G4double Phi = twopi*G4UniformRand();
>> 615 G4ThreeVector Vector(std::cos(Phi)*SinTheta, std::sin(Phi)*SinTheta, CosTheta);
>> 616 //aHadron.SetMomentum(Vector* std::sqrt(aHadron.GetTotalEnergy()*aHadron.GetTotalEnergy()-
>> 617 // aHadron.GetMass()*aHadron.GetMass()));
>> 618 G4double p2 = aHadron.GetTotalEnergy()*aHadron.GetTotalEnergy()- aHadron.GetMass()*aHadron.GetMass();
>> 619
>> 620 G4double p = 0.0;
>> 621 if ( p2 > 0.0 )
>> 622 p = std::sqrt ( p2 );
>> 623
>> 624 aHadron.SetMomentum( Vector*p );
>> 625 */
573 626
574 // fill the result <<
575 // Beware - the recoil is not necessarily in <<
576 // Can be calculated from momentum conservat <<
577 // The idea is that the particles ar emitted <<
578 // recoil is on the residual; assumption is <<
579 // the recoil. <<
580 // This needs more design @@@ <<
581 <<
582 G4bool needsSeparateRecoil = false; <<
583 G4int totalBaryonNumber = 0; <<
584 G4int totalCharge = 0; <<
585 G4ThreeVector totalMomentum(0); <<
586 if (theParticles != nullptr) { <<
587 const G4ParticleDefinition* aDef; <<
588 for (std::size_t ii0 = 0; ii0 < theParticl <<
589 aDef = (*theParticles)[ii0]->GetDefiniti <<
590 totalBaryonNumber += aDef->GetBaryonNumb <<
591 totalCharge += G4lrint(aDef->GetPDGCharg <<
592 totalMomentum += (*theParticles)[ii0]->G <<
593 } <<
594 if (totalBaryonNumber <<
595 != theBaseA + hadProjectile->GetDefin <<
596 { <<
597 needsSeparateRecoil = true; <<
598 residualA = theBaseA + hadProjectile->Ge <<
599 - totalBaryonNumber; <<
600 residualZ = theBaseZ + <<
601 G4lrint((hadProjectile->GetDefinition()->Get <<
602 } 627 }
603 } <<
604 628
605 std::size_t nPhotons = 0; << 629 // fill the result
606 if (thePhotons != nullptr) { << 630 // Beware - the recoil is not necessarily in the particles...
607 nPhotons = thePhotons->size(); << 631 // Can be calculated from momentum conservation?
608 } << 632 // The idea is that the particles ar emitted forst, and the gammas only once the
609 << 633 // recoil is on the residual; assumption is that gammas do not contribute to
610 G4DynamicParticle* theSec; << 634 // the recoil.
611 << 635 // This needs more design @@@
612 if (theParticles == nullptr) { << 636
613 theSec = new G4DynamicParticle; << 637 G4int nSecondaries = 2; // the hadron and the recoil
614 theSec->SetDefinition(aHadron.GetDefinitio << 638 G4bool needsSeparateRecoil = false;
615 theSec->SetMomentum(aHadron.GetMomentum()) << 639 G4int totalBaryonNumber = 0;
616 theResult.Get()->AddSecondary(theSec, secI << 640 G4int totalCharge = 0;
617 #ifdef G4VERBOSE << 641 G4ThreeVector totalMomentum(0);
618 if (fManager->GetDEBUG()) << 642 if(theParticles != 0)
619 G4cout << " G4ParticleHPInelasticCompFS: << 643 {
620 << theSec->GetParticleDefinition( << 644 nSecondaries = theParticles->size();
621 << " E= " << theSec->GetKineticEn << 645 const G4ParticleDefinition * aDef;
622 << theResult.Get()->GetNumberOfSe << 646 unsigned int ii0;
>> 647 for(ii0=0; ii0<theParticles->size(); ii0++)
>> 648 {
>> 649 aDef = theParticles->operator[](ii0)->GetDefinition();
>> 650 totalBaryonNumber+=aDef->GetBaryonNumber();
>> 651 totalCharge+=G4int(aDef->GetPDGCharge()+eps);
>> 652 totalMomentum += theParticles->operator[](ii0)->GetMomentum();
>> 653 }
>> 654 if(totalBaryonNumber!=G4int(theBaseA+eps+hadProjectile->GetDefinition()->GetBaryonNumber()))
>> 655 {
>> 656 needsSeparateRecoil = true;
>> 657 nSecondaries++;
>> 658 residualA = G4int(theBaseA+eps+hadProjectile->GetDefinition()->GetBaryonNumber()
>> 659 -totalBaryonNumber);
>> 660 residualZ = G4int(theBaseZ+eps+hadProjectile->GetDefinition()->GetPDGCharge()
>> 661 -totalCharge);
>> 662 }
>> 663 }
>> 664
>> 665 G4int nPhotons = 0;
>> 666 if(thePhotons!=0) { nPhotons = thePhotons->size(); }
>> 667 nSecondaries += nPhotons;
>> 668
>> 669 G4DynamicParticle * theSec;
>> 670
>> 671 if( theParticles==0 )
>> 672 {
>> 673 theSec = new G4DynamicParticle;
>> 674 theSec->SetDefinition(aHadron.GetDefinition());
>> 675 theSec->SetMomentum(aHadron.GetMomentum());
>> 676 theResult.Get()->AddSecondary(theSec);
>> 677 #ifdef G4PHPDEBUG
>> 678 if( getenv("G4ParticleHPDebug")) G4cout << this << " G4ParticleHPInelasticCompFS::BaseApply add secondary1 " << theSec->GetParticleDefinition()->GetParticleName() << " E= " << theSec->GetKineticEnergy() << " NSECO " << theResult.Get()->GetNumberOfSecondaries() << G4endl;
>> 679 #endif
>> 680
>> 681 aHadron.Lorentz(aHadron, theTarget);
>> 682 G4ReactionProduct theResidual;
>> 683 theResidual.SetDefinition(G4IonTable::GetIonTable()
>> 684 ->GetIon(static_cast<G4int>(residualZ), static_cast<G4int>(residualA), 0));
>> 685 theResidual.SetKineticEnergy(aHadron.GetKineticEnergy()*aHadron.GetMass()/theResidual.GetMass());
>> 686
>> 687 //080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #6
>> 688 //theResidual.SetMomentum(-1.*aHadron.GetMomentum());
>> 689 G4ThreeVector incidentNeutronMomentum = incidReactionProduct.GetMomentum();
>> 690 theResidual.SetMomentum(incidentNeutronMomentum - aHadron.GetMomentum());
>> 691
>> 692 theResidual.Lorentz(theResidual, -1.*theTarget);
>> 693 G4ThreeVector totalPhotonMomentum(0,0,0);
>> 694 if(thePhotons!=0)
>> 695 {
>> 696 for(i=0; i<nPhotons; i++)
>> 697 {
>> 698 totalPhotonMomentum += thePhotons->operator[](i)->GetMomentum();
>> 699 }
>> 700 }
>> 701 theSec = new G4DynamicParticle;
>> 702 theSec->SetDefinition(theResidual.GetDefinition());
>> 703 theSec->SetMomentum(theResidual.GetMomentum()-totalPhotonMomentum);
>> 704 theResult.Get()->AddSecondary(theSec);
>> 705 #ifdef G4PHPDEBUG
>> 706 if( getenv("G4ParticleHPDebug")) G4cout << this << " G4ParticleHPInelasticCompFS::BaseApply add secondary2 " << theSec->GetParticleDefinition()->GetParticleName() << " E= " << theSec->GetKineticEnergy() << " NSECO " << theResult.Get()->GetNumberOfSecondaries() << G4endl;
623 #endif 707 #endif
624 <<
625 aHadron.Lorentz(aHadron, theTarget); <<
626 G4ReactionProduct theResidual; <<
627 theResidual.SetDefinition(ionTable->GetIon <<
628 theResidual.SetKineticEnergy(aHadron.GetKi <<
629 / theResidual <<
630 <<
631 // 080612TK contribution from Benoit Pirar <<
632 // theResidual.SetMomentum(-1.*aHadron.Get <<
633 G4ThreeVector incidentNeutronMomentum = in <<
634 theResidual.SetMomentum(incidentNeutronMom <<
635 <<
636 theResidual.Lorentz(theResidual, -1. * the <<
637 G4ThreeVector totalPhotonMomentum(0, 0, 0) <<
638 if (thePhotons != nullptr) { <<
639 for (std::size_t i = 0; i < nPhotons; ++ <<
640 totalPhotonMomentum += (*thePhotons)[i <<
641 } <<
642 } 708 }
643 theSec = new G4DynamicParticle; << 709 else
644 theSec->SetDefinition(theResidual.GetDefin << 710 {
645 theSec->SetMomentum(theResidual.GetMomentu << 711 for(i0=0; i0<theParticles->size(); i0++)
646 theResult.Get()->AddSecondary(theSec, secI << 712 {
647 #ifdef G4VERBOSE << 713 theSec = new G4DynamicParticle;
648 if (fManager->GetDEBUG()) << 714 theSec->SetDefinition(theParticles->operator[](i0)->GetDefinition());
649 G4cout << this << " G4ParticleHPInelasti << 715 theSec->SetMomentum(theParticles->operator[](i0)->GetMomentum());
650 << theSec->GetParticleDefinition( << 716 theResult.Get()->AddSecondary(theSec);
651 << " E= " << theSec->GetKineticEn << 717 #ifdef G4PHPDEBUG
652 << theResult.Get()->GetNumberOfSe << 718 if( getenv("G4ParticleHPDebug")) G4cout << this << " G4ParticleHPInelasticCompFS::BaseApply add secondary3 " << theSec->GetParticleDefinition()->GetParticleName() << " E= " << theSec->GetKineticEnergy() << " NSECO " << theResult.Get()->GetNumberOfSecondaries() << G4endl;
653 #endif 719 #endif
654 } << 720 delete theParticles->operator[](i0);
655 else { << 721 }
656 for (std::size_t i0 = 0; i0 < theParticles << 722 delete theParticles;
657 theSec = new G4DynamicParticle; << 723 if(needsSeparateRecoil && residualZ!=0)
658 theSec->SetDefinition((*theParticles)[i0 << 724 {
659 theSec->SetMomentum((*theParticles)[i0]- << 725 G4ReactionProduct theResidual;
660 theResult.Get()->AddSecondary(theSec, se << 726 theResidual.SetDefinition(G4IonTable::GetIonTable()
661 #ifdef G4VERBOSE << 727 ->GetIon(static_cast<G4int>(residualZ), static_cast<G4int>(residualA), 0));
662 if (fManager->GetDEBUG()) << 728 G4double resiualKineticEnergy = theResidual.GetMass()*theResidual.GetMass();
663 G4cout << " G4ParticleHPInelasticCompF << 729 resiualKineticEnergy += totalMomentum*totalMomentum;
664 << theSec->GetParticleDefinitio << 730 resiualKineticEnergy = std::sqrt(resiualKineticEnergy) - theResidual.GetMass();
665 << " E= " << theSec->GetKinetic << 731 // cout << "Kinetic energy of the residual = "<<resiualKineticEnergy<<endl;
666 << theResult.Get()->GetNumberOf << 732 theResidual.SetKineticEnergy(resiualKineticEnergy);
>> 733
>> 734 //080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #4
>> 735 //theResidual.SetMomentum(-1.*totalMomentum);
>> 736 //G4ThreeVector incidentNeutronMomentum = incidReactionProduct.GetMomentum();
>> 737 //theResidual.SetMomentum(incidentNeutronMomentum - aHadron.GetMomentum());
>> 738 //080717 TK Comment still do NOT include photon's mometum which produce by thePhotons
>> 739 theResidual.SetMomentum( incidReactionProduct.GetMomentum() + theTarget.GetMomentum() - totalMomentum );
>> 740
>> 741 theSec = new G4DynamicParticle;
>> 742 theSec->SetDefinition(theResidual.GetDefinition());
>> 743 theSec->SetMomentum(theResidual.GetMomentum());
>> 744 theResult.Get()->AddSecondary(theSec);
>> 745 #ifdef G4PHPDEBUG
>> 746 if( getenv("G4ParticleHPDebug")) G4cout << this << " G4ParticleHPInelasticCompFS::BaseApply add secondary4 " << theSec->GetParticleDefinition()->GetParticleName() << " E= " << theSec->GetKineticEnergy() << " NSECO " << theResult.Get()->GetNumberOfSecondaries() << G4endl;
667 #endif 747 #endif
668 delete (*theParticles)[i0]; <<
669 } <<
670 delete theParticles; <<
671 if (needsSeparateRecoil && residualZ != 0) <<
672 G4ReactionProduct theResidual; <<
673 theResidual.SetDefinition(ionTable->GetI <<
674 G4double resiualKineticEnergy = theResid <<
675 resiualKineticEnergy += totalMomentum * <<
676 resiualKineticEnergy = std::sqrt(resiual <<
677 theResidual.SetKineticEnergy(resiualKine <<
678 <<
679 // 080612TK contribution from Benoit Pir <<
680 // theResidual.SetMomentum(-1.*totalMome <<
681 // G4ThreeVector incidentNeutronMomentum <<
682 // theResidual.SetMomentum(incidentNeutr <<
683 // 080717 TK Comment still do NOT includ <<
684 theResidual.SetMomentum(incidReactionPro <<
685 - totalMomentum) <<
686 748
687 theSec = new G4DynamicParticle; << 749 }
688 theSec->SetDefinition(theResidual.GetDef <<
689 theSec->SetMomentum(theResidual.GetMomen <<
690 theResult.Get()->AddSecondary(theSec, se <<
691 #ifdef G4VERBOSE <<
692 if (fManager->GetDEBUG()) <<
693 G4cout << " G4ParticleHPInelasticCompF <<
694 << theSec->GetParticleDefinitio <<
695 << " E= " << theSec->GetKinetic <<
696 << theResult.Get()->GetNumberOf <<
697 #endif <<
698 } 750 }
699 } << 751 if(thePhotons!=0)
700 if (thePhotons != nullptr) { << 752 {
701 for (std::size_t i = 0; i < nPhotons; ++i) << 753 for(i=0; i<nPhotons; i++)
702 theSec = new G4DynamicParticle; << 754 {
703 // Bug reported Chao Zhang (Chao.Zhang@u << 755 theSec = new G4DynamicParticle;
704 // 2009 theSec->SetDefinition(G4Gamma::G << 756 //Bug reported Chao Zhang (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) Feb. 25, 2009
705 theSec->SetDefinition((*thePhotons)[i]-> << 757 //theSec->SetDefinition(G4Gamma::Gamma());
706 // But never cause real effect at least << 758 theSec->SetDefinition( thePhotons->operator[](i)->GetDefinition() );
707 theSec->SetMomentum((*thePhotons)[i]->Ge << 759 //But never cause real effect at least with G4NDL3.13 TK
708 theResult.Get()->AddSecondary(theSec, se << 760 theSec->SetMomentum(thePhotons->operator[](i)->GetMomentum());
709 #ifdef G4VERBOSE << 761 theResult.Get()->AddSecondary(theSec);
710 if (fManager->GetDEBUG()) << 762 #ifdef G4PHPDEBUG
711 G4cout << " G4ParticleHPInelasticCompF << 763 if( getenv("G4ParticleHPDebug")) G4cout << this << " G4ParticleHPInelasticCompFS::BaseApply add secondary5 " << theSec->GetParticleDefinition()->GetParticleName() << " E= " << theSec->GetKineticEnergy() << " NSECO " << theResult.Get()->GetNumberOfSecondaries() << G4endl;
712 << theSec->GetParticleDefinitio <<
713 << " E= " << theSec->GetKinetic <<
714 << theResult.Get()->GetNumberOf <<
715 #endif 764 #endif
716 765
717 delete thePhotons->operator[](i); << 766 delete thePhotons->operator[](i);
>> 767 }
>> 768 // some garbage collection
>> 769 delete thePhotons;
718 } 770 }
719 // some garbage collection <<
720 delete thePhotons; <<
721 } <<
722 771
723 G4ParticleDefinition* targ_pd = ionTable->Ge << 772 //080721
724 G4LorentzVector targ_4p_lab( << 773 G4ParticleDefinition* targ_pd = G4IonTable::GetIonTable()->GetIon ( (G4int)theBaseZ , (G4int)theBaseA , 0.0 );
725 theTarget.GetMomentum(), << 774 G4LorentzVector targ_4p_lab ( theTarget.GetMomentum() , std::sqrt( targ_pd->GetPDGMass()*targ_pd->GetPDGMass() + theTarget.GetMomentum().mag2() ) );
726 std::sqrt(targ_pd->GetPDGMass() * targ_pd- << 775 G4LorentzVector proj_4p_lab = theTrack.Get4Momentum();
727 G4LorentzVector proj_4p_lab = theTrack.Get4M << 776 G4LorentzVector init_4p_lab = proj_4p_lab + targ_4p_lab;
728 G4LorentzVector init_4p_lab = proj_4p_lab + << 777 adjust_final_state ( init_4p_lab );
729 adjust_final_state(init_4p_lab); <<
730 778
731 // clean up the primary neutron << 779 // clean up the primary neutron
732 theResult.Get()->SetStatusChange(stopAndKill << 780 theResult.Get()->SetStatusChange(stopAndKill);
733 } 781 }
734 782
735 // Re-implemented by E. Mendoza (2019). Isotro <<
736 // proj: projectile in target-rest-frame (inp <<
737 // targ: target in target-rest-frame (input) <<
738 // product: secondary particle in target-rest <<
739 // resExcitationEnergy: excitation energy of <<
740 // <<
741 void G4ParticleHPInelasticCompFS::two_body_rea <<
742 <<
743 <<
744 <<
745 { <<
746 // CMS system: <<
747 G4ReactionProduct theCMS = *proj + *targ; <<
748 783
749 // Residual definition: <<
750 G4int resZ = G4lrint((proj->GetDefinition()- <<
751 - product->GetDefinition()->GetPDGCharge <<
752 G4int resA = proj->GetDefinition()->GetBaryo <<
753 - product->GetDefinition()->Get <<
754 G4ReactionProduct theResidual; <<
755 theResidual.SetDefinition(ionTable->GetIon(r <<
756 <<
757 // CMS system: <<
758 G4ReactionProduct theCMSproj; <<
759 G4ReactionProduct theCMStarg; <<
760 theCMSproj.Lorentz(*proj, theCMS); <<
761 theCMStarg.Lorentz(*targ, theCMS); <<
762 // final Momentum in the CMS: <<
763 G4double totE = std::sqrt(theCMSproj.GetMass <<
764 + theCMSproj.GetTo <<
765 + std::sqrt(theCMStarg.GetMa <<
766 + theCMStarg.Get <<
767 G4double prodmass = product->GetMass(); <<
768 G4double resmass = theResidual.GetMass() + r <<
769 G4double fmomsquared = (totE * totE - (prodm <<
770 (totE * totE - (prodmass + resmass) * (pro <<
771 G4double fmom = (fmomsquared > 0) ? std::sqr <<
772 <<
773 // random (isotropic direction): <<
774 product->SetMomentum(fmom * G4RandomDirectio <<
775 product->SetTotalEnergy(std::sqrt(prodmass * <<
776 // Back to the LAB system: <<
777 product->Lorentz(*product, -1. * theCMS); <<
778 } <<
779 784
780 G4bool G4ParticleHPInelasticCompFS::use_nresp7 << 785 #include "G4RotationMatrix.hh"
781 << 786 void G4ParticleHPInelasticCompFS::two_body_reaction ( G4DynamicParticle* proj, G4DynamicParticle* targ, G4DynamicParticle* hadron, G4double mu )
782 <<
783 <<
784 { 787 {
785 if (aDefinition == G4Neutron::Definition()) <<
786 { <<
787 // LR: flag LR in ENDF. It indicates wheth <<
788 // it: exit channel (index of the carbon e <<
789 <<
790 // Added by A. R. Garcia (CIEMAT) to inclu <<
791 <<
792 if (LR[itt] > 0) // If there is breakup of <<
793 // MT=52-91 (it=MT-50)). <<
794 { <<
795 // Defining carbon as the target in the <<
796 G4ReactionProduct theCarbon(theTarget); <<
797 <<
798 theCarbon.SetMomentum(G4ThreeVector()); <<
799 theCarbon.SetKineticEnergy(0.); <<
800 <<
801 // Creating four reaction products. <<
802 G4ReactionProduct theProds[4]; <<
803 <<
804 // Applying C(N,N'3A) reaction mechanism <<
805 if (itt == 41) { <<
806 // QI=QM=-7.275 MeV for C-0(N,N')C-C(3 <<
807 // This is not the value of the QI of <<
808 // to the model. So we don't take it. <<
809 // we have calculated: QI=(mn+m12C)-(m <<
810 nresp71_model.ApplyMechanismI_NBeA2A(b <<
811 // N+C --> A[0]+9BE* | 9BE* --> N[1]+8 <<
812 } <<
813 else { <<
814 nresp71_model.ApplyMechanismII_ACN2A(b <<
815 // N+C --> N'[0]+C* | C* --> A[1]+8BE <<
816 } <<
817 <<
818 // Returning to the reference frame wher <<
819 for (auto& theProd : theProds) { <<
820 theProd.Lorentz(theProd, -1. * theTarg <<
821 theResult.Get()->AddSecondary( <<
822 new G4DynamicParticle(theProd.GetDef <<
823 } <<
824 788
825 // Killing the primary neutron. << 789 // Target rest flame
826 theResult.Get()->SetStatusChange(stopAnd << 790 // 4vector in targ rest frame;
>> 791 // targ could have excitation energy (photon energy will be emiited) tricky but,,,
>> 792
>> 793 G4LorentzVector before = proj->Get4Momentum() + targ->Get4Momentum();
>> 794
>> 795 G4ThreeVector p3_proj = proj->GetMomentum();
>> 796 G4ThreeVector d = p3_proj.unit();
>> 797 G4RotationMatrix rot;
>> 798 G4RotationMatrix rot1;
>> 799 rot1.setPhi( CLHEP::pi/2 + d.phi() );
>> 800 G4RotationMatrix rot2;
>> 801 rot2.setTheta( d.theta() );
>> 802 rot=rot2*rot1;
>> 803 proj->SetMomentum( rot*p3_proj );
>> 804
>> 805 // Now proj only has pz component;
>> 806
>> 807 // mu in CM system
>> 808
>> 809 //Valid only for neutron incidence
>> 810 G4DynamicParticle* residual = new G4DynamicParticle ( G4IonTable::GetIonTable()->GetIon ( (G4int)( targ->GetDefinition()->GetPDGCharge() - hadron->GetDefinition()->GetPDGCharge() ) , (G4int)(targ->GetDefinition()->GetBaryonNumber() - hadron->GetDefinition()->GetBaryonNumber()+1) , 0 ) , G4ThreeVector(0) );
>> 811
>> 812 G4double Q = proj->GetDefinition()->GetPDGMass() + targ->GetDefinition()->GetPDGMass()
>> 813 - ( hadron->GetDefinition()->GetPDGMass() + residual->GetDefinition()->GetPDGMass() );
>> 814
>> 815 // Non Relativistic Case
>> 816 G4double A = targ->GetDefinition()->GetPDGMass() / proj->GetDefinition()->GetPDGMass();
>> 817 G4double AA = hadron->GetDefinition()->GetPDGMass() / proj->GetDefinition()->GetPDGMass();
>> 818 G4double E1 = proj->GetKineticEnergy();
>> 819
>> 820 // 101111
>> 821 // In _nat_ data (Q+E1) could become negative value, following line is safty for this case.
>> 822 //if ( (Q+E1) < 0 )
>> 823 if ( ( 1 + (1+A)/A*Q/E1 ) < 0 )
>> 824 {
>> 825 // 1.0e-6 eV is additional safty for numeric precision
>> 826 Q = -( A/(1+A)*E1 ) + 1.0e-6*CLHEP::eV;
>> 827 }
>> 828
>> 829 G4double beta = std::sqrt ( A*(A+1-AA)/AA*( 1 + (1+A)/A*Q/E1 ) );
>> 830 G4double gamma = AA/(A+1-AA)*beta;
>> 831 G4double E3 = AA/G4Pow::GetInstance()->powN((1+A),2)*(beta*beta+1+2*beta*mu)*E1;
>> 832 G4double omega3 = (1+beta*mu)/std::sqrt(beta*beta+1+2*beta*mu);
>> 833 if ( omega3 > 1.0 ) omega3 = 1.0;
>> 834
>> 835 G4double E4 = (A+1-AA)/G4Pow::GetInstance()->powN((1+A),2)*(gamma*gamma+1-2*gamma*mu)*E1;
>> 836 G4double omega4 = (1-gamma*mu)/std::sqrt(gamma*gamma+1-2*gamma*mu);
>> 837 if ( omega4 > 1.0 ) omega4 = 1.0;
>> 838
>> 839 hadron->SetKineticEnergy ( E3 );
>> 840
>> 841 G4double M = hadron->GetDefinition()->GetPDGMass();
>> 842 G4double pmag = std::sqrt ((E3+M)*(E3+M)-M*M) ;
>> 843 G4ThreeVector p ( 0 , pmag*std::sqrt(1-omega3*omega3), pmag*omega3 );
>> 844
>> 845 G4double M4 = residual->GetDefinition()->GetPDGMass();
>> 846 G4double pmag4 = std::sqrt ((E4+M4)*(E4+M4)-M4*M4) ;
>> 847 G4ThreeVector p4 ( 0 , -pmag4*std::sqrt(1-omega4*omega4), pmag4*omega4 );
>> 848
>> 849 // Rotate to orginal target rest flame.
>> 850 p *= rot.inverse();
>> 851 hadron->SetMomentum( p );
>> 852 // Now hadron had 4 momentum in target rest flame
>> 853
>> 854 // TypeA
>> 855 p4 *= rot.inverse();
>> 856 residual->SetMomentum ( p4 );
>> 857
>> 858 //TypeB1
>> 859 //residual->Set4Momentum ( p4_residual );
>> 860 //TypeB2
>> 861 //residual->SetMomentum ( p4_residual.v() );
>> 862
>> 863 // Type A make difference in Momenutum
>> 864 // Type B1 make difference in Mass of residual
>> 865 // Type B2 make difference in total energy.
827 866
828 return true; << 867 delete residual;
829 } <<
830 } <<
831 else if (aDefinition == G4Alpha::Definition( <<
832 { <<
833 // Added by A. R. Garcia (CIEMAT) to inclu <<
834 868
835 if (LR[itt] == 0) // If Z=6, an alpha part <<
836 // residual nucleus LR(f <<
837 { <<
838 // Defining carbon as the target in the <<
839 G4ReactionProduct theCarbon(theTarget); <<
840 theCarbon.SetMomentum(G4ThreeVector()); <<
841 theCarbon.SetKineticEnergy(0.); <<
842 <<
843 // Creating four reaction products. <<
844 G4ReactionProduct theProds[2]; <<
845 <<
846 // Applying C(N,A)9BE reaction mechanism <<
847 nresp71_model.ApplyMechanismABE(boosted, <<
848 // N+C --> A[0]+9BE[1]. <<
849 <<
850 for (auto& theProd : theProds) { <<
851 // Returning to the system of referenc <<
852 theProd.Lorentz(theProd, -1. * theTarg <<
853 theResult.Get()->AddSecondary( <<
854 new G4DynamicParticle(theProd.GetDef <<
855 } <<
856 <<
857 // Killing the primary neutron. <<
858 theResult.Get()->SetStatusChange(stopAnd <<
859 return true; <<
860 } <<
861 G4Exception("G4ParticleHPInelasticCompFS:: <<
862 FatalException, "Alpha product <<
863 } <<
864 return false; <<
865 } 869 }
866 870