| 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 // 09-May-06 fix in Sample by T. Koi 30 // 09-May-06 fix in Sample by T. Koi
31 // 080318 Fix Compilation warnings - gcc-4.3.0 31 // 080318 Fix Compilation warnings - gcc-4.3.0 by T. Koi
32 // (This fix has a real effect to the c << 32 // (This fix has a real effect to the code.)
33 // 080409 Fix div0 error with G4FPE by T. Koi 33 // 080409 Fix div0 error with G4FPE by T. Koi
34 // 080612 Fix contribution from Benoit Pirard 34 // 080612 Fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #1
35 // 080714 Limiting the sum of energy of second 35 // 080714 Limiting the sum of energy of secondary particles by T. Koi
36 // 080801 Fix div0 error wiht G4FPE and memory 36 // 080801 Fix div0 error wiht G4FPE and memory leak by T. Koi
37 // 081024 G4NucleiPropertiesTable:: to G4Nucle 37 // 081024 G4NucleiPropertiesTable:: to G4NucleiProperties::
38 // 38 //
39 // P. Arce, June-2014 Conversion neutron_hp to 39 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
40 // 40 //
41 // June-2019 - E. Mendoza --> redefinition of << 41 // June-2019 - E. Mendoza --> redefinition of the residual mass to consider incident particles different than neutrons.
42 // different than neutrons. <<
43 // <<
44 // V. Ivanchenko, July-2023 Basic revision of <<
45 // <<
46 42
47 #include "G4ParticleHPContAngularPar.hh" 43 #include "G4ParticleHPContAngularPar.hh"
48 <<
49 #include "G4ParticleDefinition.hh" <<
50 #include "G4Alpha.hh" <<
51 #include "G4Deuteron.hh" <<
52 #include "G4Electron.hh" <<
53 #include "G4Gamma.hh" <<
54 #include "G4He3.hh" <<
55 #include "G4IonTable.hh" <<
56 #include "G4Neutron.hh" <<
57 #include "G4NucleiProperties.hh" <<
58 #include "G4ParticleHPKallbachMannSyst.hh" <<
59 #include "G4ParticleHPLegendreStore.hh" <<
60 #include "G4ParticleHPManager.hh" <<
61 #include "G4ParticleHPVector.hh" <<
62 #include "G4PhysicalConstants.hh" 44 #include "G4PhysicalConstants.hh"
>> 45 #include "G4SystemOfUnits.hh"
>> 46 #include "G4ParticleHPLegendreStore.hh"
>> 47 #include "G4Gamma.hh"
>> 48 #include "G4Electron.hh"
63 #include "G4Positron.hh" 49 #include "G4Positron.hh"
>> 50 #include "G4Neutron.hh"
64 #include "G4Proton.hh" 51 #include "G4Proton.hh"
65 #include "G4SystemOfUnits.hh" << 52 #include "G4Deuteron.hh"
66 #include "G4Triton.hh" 53 #include "G4Triton.hh"
67 << 54 #include "G4He3.hh"
>> 55 #include "G4Alpha.hh"
>> 56 #include "G4ParticleHPVector.hh"
>> 57 #include "G4NucleiProperties.hh"
>> 58 #include "G4ParticleHPKallbachMannSyst.hh"
>> 59 #include "G4IonTable.hh"
68 #include <set> 60 #include <set>
69 #include <vector> << 61
70 << 62 G4ParticleHPContAngularPar::G4ParticleHPContAngularPar( G4ParticleDefinition* projectile )
71 G4ParticleHPContAngularPar::G4ParticleHPContAn << 63 {
72 { << 64 theAngular = 0;
73 theProjectile = (nullptr == p) ? G4Neutron:: << 65 if ( fCache.Get() == 0 ) cacheInit();
74 toBeCached v; << 66 fCache.Get()->currentMeanEnergy = -2;
75 fCache.Put(v); << 67 fCache.Get()->fresh = true;
76 if (G4ParticleHPManager::GetInstance()->GetD <<
77 } <<
78 <<
79 G4ParticleHPContAngularPar::G4ParticleHPContAn <<
80 { <<
81 theEnergy = val.theEnergy; <<
82 nEnergies = val.nEnergies; <<
83 nDiscreteEnergies = val.nDiscreteEnergies; <<
84 nAngularParameters = val.nAngularParameters; <<
85 theProjectile = val.theProjectile; <<
86 theManager = val.theManager; <<
87 theInt = val.theInt; <<
88 adjustResult = val.adjustResult; <<
89 theMinEner = val.theMinEner; <<
90 theMaxEner = val.theMaxEner; <<
91 theEnergiesTransformed = val.theEnergiesTran <<
92 theDiscreteEnergies = val.theDiscreteEnergie <<
93 theDiscreteEnergiesOwn = val.theDiscreteEner <<
94 toBeCached v; <<
95 fCache.Put(v); <<
96 const std::size_t esize = nEnergies > 0 ? nE <<
97 theAngular = new G4ParticleHPList[esize]; <<
98 for (G4int ie = 0; ie < nEnergies; ++ie) { <<
99 theAngular[ie].SetLabel(val.theAngular[ie] <<
100 for (G4int ip = 0; ip < nAngularParameters <<
101 theAngular[ie].SetValue(ip, val.theAngul <<
102 } <<
103 } <<
104 } <<
105 <<
106 G4ParticleHPContAngularPar::~G4ParticleHPContA <<
107 { <<
108 delete[] theAngular; <<
109 } <<
110 <<
111 void G4ParticleHPContAngularPar::Init(std::ist <<
112 { <<
113 adjustResult = true; 68 adjustResult = true;
114 if (G4ParticleHPManager::GetInstance()->GetD << 69 if ( std::getenv( "G4PHP_DO_NOT_ADJUST_FINAL_STATE" ) ) adjustResult = false;
115 70
116 theProjectile = (nullptr == p) ? G4Neutron:: << 71 theMinEner = DBL_MAX;
>> 72 theMaxEner = -DBL_MAX;
>> 73 theProjectile = projectile;
117 74
118 aDataFile >> theEnergy >> nEnergies >> nDisc << 75 theEnergy = 0.0;
119 theEnergy *= eV; << 76 nEnergies = 0;
120 const std::size_t esize = nEnergies > 0 ? nE << 77 nDiscreteEnergies = 0;
121 theAngular = new G4ParticleHPList[esize]; << 78 nAngularParameters = 0;
122 G4double sEnergy; <<
123 for (G4int i = 0; i < nEnergies; ++i) { <<
124 aDataFile >> sEnergy; <<
125 sEnergy *= eV; <<
126 theAngular[i].SetLabel(sEnergy); <<
127 theAngular[i].Init(aDataFile, nAngularPara <<
128 theMinEner = std::min(theMinEner, sEnergy) <<
129 theMaxEner = std::max(theMaxEner, sEnergy) <<
130 } <<
131 } 79 }
132 80
133 G4ReactionProduct* G4ParticleHPContAngularPar: << 81 void G4ParticleHPContAngularPar::Init(std::istream & aDataFile, G4ParticleDefinition* projectile)
134 << 82 {
135 << 83 adjustResult = true;
136 { << 84 if ( std::getenv( "G4PHP_DO_NOT_ADJUST_FINAL_STATE" ) ) adjustResult = false;
137 // The following line is needed because it m << 85
138 adjustResult = true; << 86 theProjectile = projectile;
139 if (G4ParticleHPManager::GetInstance()->GetD << 87
140 << 88 aDataFile >> theEnergy >> nEnergies >> nDiscreteEnergies >> nAngularParameters;
141 auto result = new G4ReactionProduct; << 89 /*if( std::getenv("G4PHPTEST") )*/
142 auto Z = static_cast<G4int>(massCode / 1000) << 90 theEnergy *= eV;
143 auto A = static_cast<G4int>(massCode - 1000 << 91 theAngular = new G4ParticleHPList [nEnergies];
144 if (massCode == 0) { << 92 for(G4int i=0; i<nEnergies; i++)
145 result->SetDefinition(G4Gamma::Gamma()); << 93 {
146 } << 94 G4double sEnergy;
147 else if (A == 0) { << 95 aDataFile >> sEnergy;
148 result->SetDefinition(G4Electron::Electron << 96 sEnergy*=eV;
149 if (Z == 1) result->SetDefinition(G4Positr << 97 theAngular[i].SetLabel(sEnergy);
150 } << 98 theAngular[i].Init(aDataFile, nAngularParameters, 1.);
151 else if (A == 1) { << 99 theMinEner = std::min(theMinEner,sEnergy);
152 result->SetDefinition(G4Neutron::Neutron() << 100 theMaxEner = std::max(theMaxEner,sEnergy);
153 if (Z == 1) result->SetDefinition(G4Proton << 101 }
154 } <<
155 else if (A == 2) { <<
156 result->SetDefinition(G4Deuteron::Deuteron <<
157 } <<
158 else if (A == 3) { <<
159 result->SetDefinition(G4Triton::Triton()); <<
160 if (Z == 2) result->SetDefinition(G4He3::H <<
161 } <<
162 else if (A == 4) { <<
163 result->SetDefinition(G4Alpha::Alpha()); <<
164 if (Z != 2) <<
165 throw G4HadronicException(__FILE__, __LI <<
166 "G4ParticleHPC <<
167 } <<
168 else { <<
169 result->SetDefinition(G4IonTable::GetIonTa <<
170 } 102 }
171 103
172 G4int i(0); << 104 G4ReactionProduct *
173 G4int it(0); << 105 G4ParticleHPContAngularPar::Sample(G4double anEnergy, G4double massCode, G4double /*targetMass*/,
174 G4double fsEnergy(0); << 106 G4int angularRep, G4int /*interpolE*/ )
175 G4double cosTh(0); << 107 {
176 /* << 108 if( std::getenv("G4PHPTEST") ) G4cout << " G4ParticleHPContAngularPar::Sample " << anEnergy << " " << massCode << " " << angularRep << G4endl; //GDEB
177 G4cout << "G4ParticleHPContAngularPar::Sampl << 109 if ( fCache.Get() == 0 ) cacheInit();
178 << " angularRep=" << angularRep << " << 110 G4ReactionProduct * result = new G4ReactionProduct;
179 << " Ne=" << nEnergies << G4endl; << 111 G4int Z = static_cast<G4int>(massCode/1000);
180 */ << 112 G4int A = static_cast<G4int>(massCode-1000*Z);
181 if (angularRep == 1) { << 113 if(massCode==0)
182 if (nDiscreteEnergies != 0) { << 114 {
183 // 1st check remaining_energy << 115 result->SetDefinition(G4Gamma::Gamma());
184 // if this is the first set it. (How?) << 116 }
185 if (fCache.Get().fresh) { << 117 else if(A==0)
186 // Discrete Lines, larger energies com << 118 {
187 // Continues Emssions, low to high << 119 result->SetDefinition(G4Electron::Electron());
188 fCache.Get().remaining_energy = << 120 if(Z==1) result->SetDefinition(G4Positron::Positron());
189 std::max(theAngular[0].GetLabel(), t << 121 }
190 fCache.Get().fresh = false; << 122 else if(A==1)
191 } << 123 {
192 << 124 result->SetDefinition(G4Neutron::Neutron());
193 // Cheating for small remaining_energy << 125 if(Z==1) result->SetDefinition(G4Proton::Proton());
194 // Temporary solution << 126 }
195 if (nDiscreteEnergies == nEnergies) { << 127 else if(A==2)
196 fCache.Get().remaining_energy = << 128 {
197 std::max(fCache.Get().remaining_ener << 129 result->SetDefinition(G4Deuteron::Deuteron());
198 theAngular[nDiscreteEnergie << 130 }
199 } << 131 else if(A==3)
200 else { << 132 {
201 G4double cont_min = 0.0; << 133 result->SetDefinition(G4Triton::Triton());
202 for (G4int j = nDiscreteEnergies; j < << 134 if(Z==2) result->SetDefinition(G4He3::He3());
203 cont_min = theAngular[j].GetLabel(); << 135 }
204 if (theAngular[j].GetValue(0) != 0.0 << 136 else if(A==4)
205 } << 137 {
206 fCache.Get().remaining_energy = std::m << 138 result->SetDefinition(G4Alpha::Alpha());
207 fCache.Get().remaining_energy, std:: << 139 if(Z!=2) throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPContAngularPar: Unknown ion case 1");
208 << 140 }
209 } << 141 else
210 << 142 {
211 G4double random = G4UniformRand(); << 143 result->SetDefinition(G4IonTable::GetIonTable()->GetIon(Z,A,0));
212 auto running = new G4double[nEnergies + << 144 }
213 running[0] = 0.0; << 145 G4int i(0);
>> 146 G4int it(0);
>> 147 G4double fsEnergy(0);
>> 148 G4double cosTh(0);
214 149
215 G4double delta; << 150 if( angularRep == 1 )
216 for (G4int j = 0; j < nDiscreteEnergies; << 151 {
217 delta = 0.0; << 152 // 080612 Fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #1
218 if (theAngular[j].GetLabel() <= fCache << 153 //if (interpolE == 2)
219 delta = theAngular[j].GetValue(0); << 154 //110609 above was wrong interupition, pointed out by E.Mendoza and D.Cano (CIMAT)
220 running[j + 1] = running[j] + delta; << 155 //Following are reviesd version written by T.Koi (SLAC)
221 } << 156 if ( nDiscreteEnergies != 0 )
>> 157 {
222 158
223 G4double tot_prob_DIS = std::max(running << 159 //1st check remaining_energy
>> 160 // if this is the first set it. (How?)
>> 161 if ( fCache.Get()->fresh == true )
>> 162 {
>> 163 //Discrete Lines, larger energies come first
>> 164 //Continues Emssions, low to high LAST
>> 165 fCache.Get()->remaining_energy = std::max ( theAngular[0].GetLabel() , theAngular[nEnergies-1].GetLabel() );
>> 166 fCache.Get()->fresh = false;
>> 167 }
>> 168
>> 169 //Cheating for small remaining_energy
>> 170 //TEMPORAL SOLUTION
>> 171 if ( nDiscreteEnergies == nEnergies )
>> 172 {
>> 173 fCache.Get()->remaining_energy = std::max ( fCache.Get()->remaining_energy , theAngular[nDiscreteEnergies-1].GetLabel() ); //Minimum Line
>> 174 }
>> 175 else
>> 176 {
>> 177 //G4double cont_min = theAngular[nDiscreteEnergies].GetLabel();
>> 178 //if ( theAngular[nDiscreteEnergies].GetLabel() == 0.0 ) cont_min = theAngular[nDiscreteEnergies+1].GetLabel();
>> 179 G4double cont_min=0.0;
>> 180 for ( G4int j = nDiscreteEnergies ; j < nEnergies ; j++ )
>> 181 {
>> 182 cont_min = theAngular[j].GetLabel();
>> 183 if ( theAngular[j].GetValue(0) != 0.0 ) break;
>> 184 }
>> 185 fCache.Get()->remaining_energy = std::max ( fCache.Get()->remaining_energy , std::min ( theAngular[nDiscreteEnergies-1].GetLabel() , cont_min ) ); //Minimum Line or grid
>> 186 }
>> 187 //
>> 188 G4double random = G4UniformRand();
224 189
225 G4double delta1; << 190 G4double * running = new G4double[nEnergies+1];
226 for (G4int j = nDiscreteEnergies; j < nE << 191 running[0] = 0.0;
227 delta1 = 0.0; <<
228 G4double e_low = 0.0; <<
229 G4double e_high = 0.0; <<
230 if (theAngular[j].GetLabel() <= fCache <<
231 delta1 = theAngular[j].GetValue(0); <<
232 <<
233 // To calculate Prob. e_low and e_high <<
234 // There are two cases: <<
235 // 1: theAngular[nDiscreteEnergies].Ge <<
236 // delta1 should be used between j- <<
237 // At j = nDiscreteEnergies (the fi <<
238 if (theAngular[j].GetLabel() != 0) { <<
239 if (j == nDiscreteEnergies) { <<
240 e_low = 0.0 / eV; <<
241 } <<
242 else { <<
243 if (j < 1) j = 1; // Protection a <<
244 e_low = theAngular[j - 1].GetLabel <<
245 } <<
246 e_high = theAngular[j].GetLabel() / <<
247 } <<
248 192
249 // 2: theAngular[nDiscreteEnergies].Ge << 193 for ( G4int j = 0 ; j < nDiscreteEnergies ; j++ )
250 // delta1 should be used between j << 194 {
251 if (theAngular[j].GetLabel() == 0.0) { << 195 G4double delta = 0.0;
252 e_low = theAngular[j].GetLabel() / e << 196 if ( theAngular[j].GetLabel() <= fCache.Get()->remaining_energy ) delta = theAngular[i].GetValue(0);
253 if (j != nEnergies - 1) { << 197 running[j+1] = running[j] + delta;
254 e_high = theAngular[j + 1].GetLabe << 198 }
>> 199 G4double tot_prob_DIS = running[ nDiscreteEnergies ];
>> 200
>> 201 for ( G4int j = nDiscreteEnergies ; j < nEnergies ; j++ )
>> 202 {
>> 203 G4double delta = 0.0;
>> 204 G4double e_low = 0.0;
>> 205 G4double e_high = 0.0;
>> 206 if ( theAngular[j].GetLabel() <= fCache.Get()->remaining_energy ) delta = theAngular[j].GetValue(0);
>> 207
>> 208 //To calculate Prob. e_low and e_high should be in eV
>> 209 //There are two case
>> 210 //1:theAngular[nDiscreteEnergies].GetLabel() != 0.0
>> 211 // delta should be used between j-1 and j
>> 212 // At j = nDiscreteEnergies (the first) e_low should be set explicitly
>> 213 if ( theAngular[j].GetLabel() != 0 )
>> 214 {
>> 215 if ( j == nDiscreteEnergies ) {
>> 216 e_low = 0.0/eV;
>> 217 } else {
>> 218 e_low = theAngular[j-1].GetLabel()/eV;
>> 219 }
>> 220 e_high = theAngular[j].GetLabel()/eV;
>> 221 }
>> 222 //2:theAngular[nDiscreteEnergies].GetLabel() == 0.0
>> 223 // delta should be used between j and j+1
>> 224 if ( theAngular[j].GetLabel() == 0.0 ) {
>> 225 e_low = theAngular[j].GetLabel()/eV;
>> 226 if ( j != nEnergies-1 ) {
>> 227 e_high = theAngular[j+1].GetLabel()/eV;
>> 228 } else {
>> 229 e_high = theAngular[j].GetLabel()/eV;
>> 230 if ( theAngular[j].GetValue(0) != 0.0 ) {
>> 231 throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPContAngularPar: Unexpected non zero value of theAngular[nEnergies-1].GetValue(0)");
>> 232 }
>> 233 }
>> 234 }
>> 235
>> 236 running[j+1] = running[j] + ( ( e_high - e_low ) * delta );
>> 237 }
>> 238 G4double tot_prob_CON = running[ nEnergies ] - running[ nDiscreteEnergies ];
>> 239
>> 240 /*
>> 241 For FPE debugging
>> 242 if (tot_prob_DIS + tot_prob_CON == 0 ) {
>> 243 G4cout << "TKDB tot_prob_DIS + tot_prob_CON " << tot_prob_DIS + tot_prob_CON << G4endl;
>> 244 G4cout << "massCode " << massCode << G4endl;
>> 245 G4cout << "nDiscreteEnergies " << nDiscreteEnergies << " nEnergies " << nEnergies << G4endl;
>> 246 for ( int j = nDiscreteEnergies ; j < nEnergies ; j++ ) {
>> 247 G4cout << j << " " << theAngular[j].GetLabel() << " " << theAngular[j].GetValue(0) << G4endl;
>> 248 }
255 } 249 }
256 else { << 250 */
257 e_high = theAngular[j].GetLabel() << 251 // Normalize random
258 } << 252 random *= (tot_prob_DIS + tot_prob_CON);
259 } << 253 //2nd Judge Discrete or not This shoudl be relatively close to 1 For safty
260 << 254 if ( random <= ( tot_prob_DIS / ( tot_prob_DIS + tot_prob_CON ) ) || nDiscreteEnergies == nEnergies )
261 running[j + 1] = running[j] + ((e_high << 255 {
262 } << 256 // Discrete Emission
263 G4double tot_prob_CON = std::max(running << 257 for ( G4int j = 0 ; j < nDiscreteEnergies ; j++ )
264 << 258 {
265 // Give up in the pathological case of n << 259 //Here we should use i+1
266 if (tot_prob_DIS == 0.0 && tot_prob_CON << 260 if ( random < running[ j+1 ] )
267 delete[] running; << 261 {
268 return result; << 262 it = j;
269 } << 263 break;
270 // Normalize random << 264 }
271 random *= (tot_prob_DIS + tot_prob_CON); << 265 }
272 // 2nd Judge Discrete or not << 266 fsEnergy = theAngular[ it ].GetLabel();
273 << 267
274 // This should be relatively close to 1 << 268 G4ParticleHPLegendreStore theStore(1);
275 if (random <= (tot_prob_DIS / (tot_prob_ << 269 theStore.Init(0,fsEnergy,nAngularParameters);
276 || nDiscreteEnergies == nEnergies) << 270 for (G4int j=0;j<nAngularParameters;j++)
>> 271 {
>> 272 theStore.SetCoeff(0,j,theAngular[it].GetValue(j));
>> 273 }
>> 274 // use it to sample.
>> 275 cosTh = theStore.SampleMax(fsEnergy);
>> 276 //Done
>> 277 }
>> 278 else
>> 279 {
>> 280 // Continuous Emission
>> 281 for ( G4int j = nDiscreteEnergies ; j < nEnergies ; j++ )
>> 282 {
>> 283 //Here we should use i
>> 284 if ( random < running[ j ] )
>> 285 {
>> 286 it = j;
>> 287 break;
>> 288 }
>> 289 }
>> 290
>> 291 G4double x1 = running[it-1];
>> 292 G4double x2 = running[it];
>> 293
>> 294 G4double y1 = 0.0;
>> 295 if ( it != nDiscreteEnergies )
>> 296 y1 = theAngular[it-1].GetLabel();
>> 297 G4double y2 = theAngular[it].GetLabel();
>> 298
>> 299 fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it),
>> 300 random,x1,x2,y1,y2);
>> 301
>> 302 G4ParticleHPLegendreStore theStore(2);
>> 303 theStore.Init(0,y1,nAngularParameters);
>> 304 theStore.Init(1,y2,nAngularParameters);
>> 305 theStore.SetManager(theManager);
>> 306 for (G4int j=0;j<nAngularParameters;j++)
>> 307 {
>> 308 G4int itt = it;
>> 309 if ( it == nDiscreteEnergies ) itt = it+1; //"This case "it-1" has data for Discrete, so we will use an extrpolate values it and it+1
>> 310 if ( it == 0 )
>> 311 {
>> 312 //Safty for unexpected it = 0;
>> 313 //G4cout << "110611 G4ParticleHPContAngularPar::Sample it = 0; invetigation required " << G4endl;
>> 314 itt = it+1;
>> 315 }
>> 316 theStore.SetCoeff(0,j,theAngular[itt-1].GetValue(j));
>> 317 theStore.SetCoeff(1,j,theAngular[itt].GetValue(j));
>> 318 }
>> 319 // use it to sample.
>> 320 cosTh = theStore.SampleMax(fsEnergy);
>> 321
>> 322 //Done
>> 323 }
>> 324
>> 325 //TK080711
>> 326 if( adjustResult ) fCache.Get()->remaining_energy -= fsEnergy;
>> 327 //TK080711
>> 328
>> 329 //080801b
>> 330 delete[] running;
>> 331 //080801b
>> 332 }
>> 333 else
277 { 334 {
278 // Discrete Emission << 335 // Only continue, TK will clean up
279 for (G4int j = 0; j < nDiscreteEnergie <<
280 // Here we should use i+1 <<
281 if (random < running[j + 1]) { <<
282 it = j; <<
283 break; <<
284 } <<
285 } <<
286 fsEnergy = theAngular[it].GetLabel(); <<
287 336
288 G4ParticleHPLegendreStore theStore(1); << 337 //080714
289 theStore.Init(0, fsEnergy, nAngularPar << 338 if ( fCache.Get()->fresh == true )
290 for (G4int j = 0; j < nAngularParamete << 339 {
291 theStore.SetCoeff(0, j, theAngular[i << 340 fCache.Get()->remaining_energy = theAngular[ nEnergies-1 ].GetLabel();
292 } << 341 fCache.Get()->fresh = false;
293 // use it to sample. << 342 }
294 cosTh = theStore.SampleMax(fsEnergy); << 343 //080714
295 // Done << 344 G4double random = G4UniformRand();
296 } << 345 G4double * running = new G4double[nEnergies];
297 else { << 346 running[0]=0;
298 // Continuous emission << 347 G4double weighted = 0;
299 for (G4int j = nDiscreteEnergies; j < << 348 for(i=1; i<nEnergies; i++)
300 // Here we should use i << 349 {
301 if (random < running[j]) { << 350 /*
302 it = j; << 351 if(i!=0)
303 break; << 352 {
304 } << 353 running[i]=running[i-1];
305 } << 354 }
306 << 355 running[i] += theInt.GetBinIntegral(theManager.GetScheme(i-1),
307 if (it < 1) it = 1; // Protection aga << 356 theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
308 << 357 theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
309 G4double x1 = running[it - 1]; << 358 weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(i-1),
310 G4double x2 = running[it]; << 359 theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
311 << 360 theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
312 G4double y1 = 0.0; << 361 */
313 if (it != nDiscreteEnergies) y1 = theA << 362
314 G4double y2 = theAngular[it].GetLabel( << 363 running[i]=running[i-1];
315 << 364 if ( fCache.Get()->remaining_energy >= theAngular[i].GetLabel() )
316 fsEnergy = theInt.Interpolate(theManag << 365 {
317 << 366 running[i] += theInt.GetBinIntegral(theManager.GetScheme(i-1),
318 G4ParticleHPLegendreStore theStore(2); << 367 theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
319 theStore.Init(0, y1, nAngularParameter << 368 theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
320 theStore.Init(1, y2, nAngularParameter << 369 weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(i-1),
321 theStore.SetManager(theManager); << 370 theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
322 G4int itt; << 371 theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
323 for (G4int j = 0; j < nAngularParamete << 372 }
324 itt = it; << 373 }
325 if (it == nDiscreteEnergies) itt = i << 374 // cash the mean energy in this distribution
326 // "This case "it-1" has data for Di << 375 //080409 TKDB
327 // it+1 << 376 if ( nEnergies == 1 || running[nEnergies-1] == 0 )
328 theStore.SetCoeff(0, j, theAngular[i << 377 fCache.Get()->currentMeanEnergy = 0.0;
329 theStore.SetCoeff(1, j, theAngular[i << 378 else
330 } << 379 {
331 // use it to sample. << 380 fCache.Get()->currentMeanEnergy = weighted/running[nEnergies-1];
332 cosTh = theStore.SampleMax(fsEnergy); << 381 }
333 << 382
334 // Done << 383 //080409 TKDB
335 } << 384 if ( nEnergies == 1 ) it = 0;
336 << 385
337 // The remaining energy needs to be lowe << 386 //080729
338 // Otherwise additional photons with too << 387 if ( running[nEnergies-1] != 0 )
339 // adjustResult condition has been remov << 388 {
340 fCache.Get().remaining_energy -= fsEnerg << 389 for ( i = 1 ; i < nEnergies ; i++ )
341 delete[] running; << 390 {
342 << 391 it = i;
343 // end (nDiscreteEnergies != 0) branch << 392 if ( random < running [ i ] / running [ nEnergies-1 ] ) break;
344 } << 393 }
345 else { << 394 }
346 // Only continue, TK will clean up << 395
347 if (fCache.Get().fresh) { << 396 //080714
348 fCache.Get().remaining_energy = theAng << 397 if ( running [ nEnergies-1 ] == 0 ) it = 0;
349 fCache.Get().fresh = false; << 398 //080714
>> 399
>> 400 if (it<nDiscreteEnergies||it==0)
>> 401 {
>> 402 if(it == 0)
>> 403 {
>> 404 fsEnergy = theAngular[it].GetLabel();
>> 405 G4ParticleHPLegendreStore theStore(1);
>> 406 theStore.Init(0,fsEnergy,nAngularParameters);
>> 407 for(i=0;i<nAngularParameters;i++)
>> 408 {
>> 409 theStore.SetCoeff(0,i,theAngular[it].GetValue(i));
>> 410 }
>> 411 // use it to sample.
>> 412 cosTh = theStore.SampleMax(fsEnergy);
>> 413 }
>> 414 else
>> 415 {
>> 416 G4double e1, e2;
>> 417 e1 = theAngular[it-1].GetLabel();
>> 418 e2 = theAngular[it].GetLabel();
>> 419 fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it),
>> 420 random,
>> 421 running[it-1]/running[nEnergies-1],
>> 422 running[it]/running[nEnergies-1],
>> 423 e1, e2);
>> 424 // fill a Legendrestore
>> 425 G4ParticleHPLegendreStore theStore(2);
>> 426 theStore.Init(0,e1,nAngularParameters);
>> 427 theStore.Init(1,e2,nAngularParameters);
>> 428 for(i=0;i<nAngularParameters;i++)
>> 429 {
>> 430 theStore.SetCoeff(0,i,theAngular[it-1].GetValue(i));
>> 431 theStore.SetCoeff(1,i,theAngular[it].GetValue(i));
>> 432 }
>> 433 // use it to sample.
>> 434 theStore.SetManager(theManager);
>> 435 cosTh = theStore.SampleMax(fsEnergy);
>> 436 }
>> 437 }
>> 438 else // continuum contribution
>> 439 {
>> 440 G4double x1 = running[it-1]/running[nEnergies-1];
>> 441 G4double x2 = running[it]/running[nEnergies-1];
>> 442 G4double y1 = theAngular[it-1].GetLabel();
>> 443 G4double y2 = theAngular[it].GetLabel();
>> 444 fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it),
>> 445 random,x1,x2,y1,y2);
>> 446 G4ParticleHPLegendreStore theStore(2);
>> 447 theStore.Init(0,y1,nAngularParameters);
>> 448 theStore.Init(1,y2,nAngularParameters);
>> 449 theStore.SetManager(theManager);
>> 450 for(i=0;i<nAngularParameters;i++)
>> 451 {
>> 452 theStore.SetCoeff(0,i,theAngular[it-1].GetValue(i));
>> 453 theStore.SetCoeff(1,i,theAngular[it].GetValue(i));
>> 454 }
>> 455 // use it to sample.
>> 456 cosTh = theStore.SampleMax(fsEnergy);
>> 457 }
>> 458 delete [] running;
>> 459
>> 460 //080714
>> 461 if( adjustResult ) fCache.Get()->remaining_energy -= fsEnergy;
>> 462 //080714
>> 463 }
>> 464 }
>> 465 else if(angularRep==2)
>> 466 {
>> 467 // first get the energy (already the right for this incoming energy)
>> 468 G4int j;
>> 469 G4double * running = new G4double[nEnergies];
>> 470 running[0]=0;
>> 471 G4double weighted = 0;
>> 472 if( std::getenv("G4PHPTEST") ) G4cout << " G4ParticleHPContAngularPar::Sample nEnergies " << nEnergies << G4endl;
>> 473 for(j=1; j<nEnergies; j++)
>> 474 {
>> 475 if(j!=0) running[j]=running[j-1];
>> 476 running[j] += theInt.GetBinIntegral(theManager.GetScheme(j-1),
>> 477 theAngular[j-1].GetLabel(), theAngular[j].GetLabel(),
>> 478 theAngular[j-1].GetValue(0), theAngular[j].GetValue(0));
>> 479 weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(j-1),
>> 480 theAngular[j-1].GetLabel(), theAngular[j].GetLabel(),
>> 481 theAngular[j-1].GetValue(0), theAngular[j].GetValue(0));
>> 482 if( std::getenv("G4PHPTEST") ) G4cout << " G4ParticleHPContAngularPar::Sample " << j << " running " << running[j]
>> 483 << " " << theManager.GetScheme(j-1) << " " << theAngular[j-1].GetLabel() << " " << theAngular[j].GetLabel() << " " << theAngular[j-1].GetValue(0) << " " << theAngular[j].GetValue(0) << G4endl; //GDEB
>> 484 }
>> 485 // cash the mean energy in this distribution
>> 486 //080409 TKDB
>> 487 //currentMeanEnergy = weighted/running[nEnergies-1];
>> 488 if ( nEnergies == 1 )
>> 489 fCache.Get()->currentMeanEnergy = 0.0;
>> 490 else
>> 491 fCache.Get()->currentMeanEnergy = weighted/running[nEnergies-1];
>> 492
>> 493 G4int itt(0);
>> 494 G4double randkal = G4UniformRand();
>> 495 //080409 TKDB
>> 496 //for(i=0; i<nEnergies; i++)
>> 497 for(j=1; j<nEnergies; j++)
>> 498 {
>> 499 itt = j;
>> 500 if(randkal<running[j]/running[nEnergies-1]) break;
350 } 501 }
351 << 502
>> 503 // interpolate the secondary energy.
>> 504 G4double x, x1,x2,y1,y2;
>> 505 if(itt==0) itt=1;
>> 506 x = randkal*running[nEnergies-1];
>> 507 x1 = running[itt-1];
>> 508 x2 = running[itt];
>> 509 G4double compoundFraction;
>> 510 // interpolate energy
>> 511 y1 = theAngular[itt-1].GetLabel();
>> 512 y2 = theAngular[itt].GetLabel();
>> 513 fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(itt-1),
>> 514 x, x1,x2,y1,y2);
>> 515 if( std::getenv("G4PHPTEST") ) G4cout << itt << " G4particleHPContAngularPar fsEnergy " << fsEnergy << " " << theManager.GetInverseScheme(itt-1) << " x " << x << " " << x1 << " " << x2 << " y " << y1 << " " << y2 << G4endl; //GDEB
>> 516 // for theta interpolate the compoundFractions
>> 517 G4double cLow = theAngular[itt-1].GetValue(1);
>> 518 G4double cHigh = theAngular[itt].GetValue(1);
>> 519 compoundFraction = theInt.Interpolate(theManager.GetScheme(itt),
>> 520 fsEnergy, y1, y2, cLow,cHigh);
>> 521
>> 522 if ( compoundFraction > 1.0 ) compoundFraction = 1.0; // Protection against unphysical interpolation
>> 523
>> 524 if( std::getenv("G4PHPTEST") ) G4cout << itt << " G4particleHPContAngularPar compoundFraction " << compoundFraction << " E " << fsEnergy << " " << theManager.GetScheme(itt) << " ener " << fsEnergy << " y " << y1 << " " << y2 << " cLH " << cLow << " " << cHigh << G4endl; //GDEB
>> 525 delete [] running;
>> 526
>> 527 // get cosTh
>> 528 G4double incidentEnergy = anEnergy;
>> 529 G4double incidentMass = theProjectile->GetPDGMass();
>> 530 G4double productEnergy = fsEnergy;
>> 531 G4double productMass = result->GetMass();
>> 532 G4int targetZ = G4int(fCache.Get()->theTargetCode/1000);
>> 533 G4int targetA = G4int(fCache.Get()->theTargetCode-1000*targetZ);
>> 534 // To correspond to natural composition (-nat-) data files.
>> 535 if ( targetA == 0 )
>> 536 targetA = G4int ( fCache.Get()->theTarget->GetMass()/amu_c2 + 0.5 );
>> 537 G4double targetMass = fCache.Get()->theTarget->GetMass();
>> 538 G4int incidentA=G4int(incidentMass/amu_c2 + 0.5 );
>> 539 G4int incidentZ=G4int(theProjectile->GetPDGCharge()+ 0.5 );
>> 540 G4int residualA = targetA+incidentA-A;
>> 541 G4int residualZ = targetZ+incidentZ-Z;
>> 542 G4double residualMass =G4NucleiProperties::GetNuclearMass( residualA , residualZ );
>> 543 G4ParticleHPKallbachMannSyst theKallbach(compoundFraction,
>> 544 incidentEnergy, incidentMass,
>> 545 productEnergy, productMass,
>> 546 residualMass, residualA, residualZ,
>> 547 targetMass, targetA, targetZ,
>> 548 incidentA,incidentZ,A,Z);
>> 549 cosTh = theKallbach.Sample(anEnergy);
>> 550 if( std::getenv("G4PHPTEST") ) G4cout << " G4ParticleHPKallbachMannSyst::Sample resulttest " << cosTh << G4endl; //GDEB
>> 551 }
>> 552 else if(angularRep>10&&angularRep<16)
>> 553 {
352 G4double random = G4UniformRand(); 554 G4double random = G4UniformRand();
353 auto running = new G4double[nEnergies]; << 555 G4double * running = new G4double[nEnergies];
354 running[0] = 0; << 556 running[0]=0;
355 G4double weighted = 0; 557 G4double weighted = 0;
356 for (i = 1; i < nEnergies; i++) { << 558 for(i=1; i<nEnergies; i++)
357 running[i] = running[i - 1]; << 559 {
358 if (fCache.Get().remaining_energy >= t << 560 if(i!=0) running[i]=running[i-1];
359 running[i] += theInt.GetBinIntegral( << 561 running[i] += theInt.GetBinIntegral(theManager.GetScheme(i-1),
360 theManager.GetScheme(i - 1), theAn << 562 theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
361 theAngular[i - 1].GetValue(0), the << 563 theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
362 weighted += theInt.GetWeightedBinInt << 564 weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(i-1),
363 theManager.GetScheme(i - 1), theAn << 565 theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
364 theAngular[i - 1].GetValue(0), the << 566 theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
365 } << 567 }
366 } << 568 // cash the mean energy in this distribution
367 << 569 //currentMeanEnergy = weighted/running[nEnergies-1];
368 // Cache the mean energy in this distrib << 570 if ( nEnergies == 1 )
369 if (nEnergies == 1 || running[nEnergies << 571 fCache.Get()->currentMeanEnergy = 0.0;
370 fCache.Get().currentMeanEnergy = 0.0; << 572 else
371 } << 573 fCache.Get()->currentMeanEnergy = weighted/running[nEnergies-1];
372 else { << 574
373 fCache.Get().currentMeanEnergy = weigh << 575 //080409 TKDB
374 } << 576 if ( nEnergies == 1 ) it = 0;
375 << 577 //for(i=0; i<nEnergies; i++)
376 if (nEnergies == 1) it = 0; << 578 for(i=1; i<nEnergies; i++)
377 if (running[nEnergies - 1] != 0) { << 579 {
378 for (i = 1; i < nEnergies; i++) { << 580 it = i;
379 it = i; << 581 if(random<running[i]/running[nEnergies-1]) break;
380 if (random < running[i] / running[nE <<
381 } <<
382 } 582 }
383 << 583 if(it<nDiscreteEnergies||it==0)
384 if (running[nEnergies - 1] == 0) it = 0; << 584 {
385 if (it < nDiscreteEnergies || it == 0) { << 585 if(it==0)
386 if (it == 0) { << 586 {
387 fsEnergy = theAngular[it].GetLabel() << 587 fsEnergy = theAngular[0].GetLabel();
388 G4ParticleHPLegendreStore theStore(1 << 588 G4ParticleHPVector theStore;
389 theStore.Init(0, fsEnergy, nAngularP << 589 G4int aCounter = 0;
390 for (i = 0; i < nAngularParameters; << 590 for(G4int j=1; j<nAngularParameters; j+=2)
391 theStore.SetCoeff(0, i, theAngular << 591 {
>> 592 theStore.SetX(aCounter, theAngular[0].GetValue(j));
>> 593 theStore.SetY(aCounter, theAngular[0].GetValue(j+1));
>> 594 aCounter++;
392 } 595 }
393 // use it to sample. << 596 G4InterpolationManager aMan;
394 cosTh = theStore.SampleMax(fsEnergy) << 597 aMan.Init(angularRep-10, nAngularParameters-1);
>> 598 theStore.SetInterpolationManager(aMan);
>> 599 cosTh = theStore.Sample();
395 } 600 }
396 else { << 601 else
397 G4double e1, e2; << 602 {
398 e1 = theAngular[it - 1].GetLabel(); << 603 fsEnergy = theAngular[it].GetLabel();
399 e2 = theAngular[it].GetLabel(); << 604 G4ParticleHPVector theStore;
400 fsEnergy = theInt.Interpolate(theMan << 605 G4InterpolationManager aMan;
401 runnin << 606 aMan.Init(angularRep-10, nAngularParameters-1);
402 runnin << 607 theStore.SetInterpolationManager(aMan); // Store interpolates f(costh)
403 // fill a Legendrestore << 608 G4InterpolationScheme currentScheme = theManager.GetInverseScheme(it);
404 G4ParticleHPLegendreStore theStore(2 << 609 G4int aCounter = 0;
405 theStore.Init(0, e1, nAngularParamet << 610 for(G4int j=1; j<nAngularParameters; j+=2)
406 theStore.Init(1, e2, nAngularParamet << 611 {
407 for (i = 0; i < nAngularParameters; << 612 theStore.SetX(aCounter, theAngular[it].GetValue(j));
408 theStore.SetCoeff(0, i, theAngular << 613 theStore.SetY(aCounter, theInt.Interpolate(currentScheme,
409 theStore.SetCoeff(1, i, theAngular << 614 random,
>> 615 running[it-1]/running[nEnergies-1],
>> 616 running[it]/running[nEnergies-1],
>> 617 theAngular[it-1].GetValue(j+1),
>> 618 theAngular[it].GetValue(j+1)));
>> 619 aCounter++;
410 } 620 }
411 // use it to sample. << 621 cosTh = theStore.Sample();
412 theStore.SetManager(theManager); <<
413 cosTh = theStore.SampleMax(fsEnergy) <<
414 } 622 }
415 } 623 }
416 else { // continuum contribution << 624 else
417 G4double x1 = running[it - 1] / runnin << 625 {
418 G4double x2 = running[it] / running[nE << 626 G4double x1 = running[it-1]/running[nEnergies-1];
419 G4double y1 = theAngular[it - 1].GetLa << 627 G4double x2 = running[it]/running[nEnergies-1];
>> 628 G4double y1 = theAngular[it-1].GetLabel();
420 G4double y2 = theAngular[it].GetLabel( 629 G4double y2 = theAngular[it].GetLabel();
421 fsEnergy = theInt.Interpolate(theManag << 630 fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it),
422 G4ParticleHPLegendreStore theStore(2); << 631 random,x1,x2,y1,y2);
423 theStore.Init(0, y1, nAngularParameter << 632 G4ParticleHPVector theBuff1;
424 theStore.Init(1, y2, nAngularParameter << 633 G4ParticleHPVector theBuff2;
425 theStore.SetManager(theManager); << 634 G4InterpolationManager aMan;
426 for (i = 0; i < nAngularParameters; i+ << 635 aMan.Init(angularRep-10, nAngularParameters-1);
427 theStore.SetCoeff(0, i, theAngular[i << 636 // theBuff1.SetInterpolationManager(aMan); // Store interpolates f(costh)
428 theStore.SetCoeff(1, i, theAngular[i << 637 // theBuff2.SetInterpolationManager(aMan); // Store interpolates f(costh)
>> 638 // Bug Report #1366 from L. Russell
>> 639 //for(i=0; i<nAngularParameters; i++) // i=1 ist wichtig!
>> 640 //{
>> 641 // theBuff1.SetX(i, theAngular[it-1].GetValue(i));
>> 642 // theBuff1.SetY(i, theAngular[it-1].GetValue(i+1));
>> 643 // theBuff2.SetX(i, theAngular[it].GetValue(i));
>> 644 // theBuff2.SetY(i, theAngular[it].GetValue(i+1));
>> 645 // i++;
>> 646 //}
>> 647 {
>> 648 G4int j;
>> 649 for(i=0,j=1; i<nAngularParameters; i++,j+=2)
>> 650 {
>> 651 theBuff1.SetX(i, theAngular[it-1].GetValue(j));
>> 652 theBuff1.SetY(i, theAngular[it-1].GetValue(j+1));
>> 653 theBuff2.SetX(i, theAngular[it].GetValue(j));
>> 654 theBuff2.SetY(i, theAngular[it].GetValue(j+1));
429 } 655 }
430 // use it to sample. <<
431 cosTh = theStore.SampleMax(fsEnergy); <<
432 } <<
433 delete[] running; <<
434 <<
435 // The remaining energy needs to be lowe <<
436 // *any* case. Otherwise additional pho <<
437 // produced - therefore the adjustResul <<
438 <<
439 fCache.Get().remaining_energy -= fsEnerg <<
440 // end if (nDiscreteEnergies != 0) <<
441 } <<
442 // end of (angularRep == 1) branch <<
443 } <<
444 else if (angularRep == 2) { <<
445 // first get the energy (already the right <<
446 G4int j; <<
447 auto running = new G4double[nEnergies]; <<
448 running[0] = 0; <<
449 G4double weighted = 0; <<
450 for (j = 1; j < nEnergies; ++j) { <<
451 if (j != 0) running[j] = running[j - 1]; <<
452 running[j] += theInt.GetBinIntegral(theM <<
453 theA <<
454 theA <<
455 weighted += theInt.GetWeightedBinIntegra <<
456 theManager.GetScheme(j - 1), theAngula <<
457 theAngular[j - 1].GetValue(0), theAngu <<
458 } <<
459 <<
460 // Cache the mean energy in this distribut <<
461 if (nEnergies == 1) <<
462 fCache.Get().currentMeanEnergy = 0.0; <<
463 else <<
464 fCache.Get().currentMeanEnergy = weighte <<
465 <<
466 G4int itt(0); <<
467 G4double randkal = G4UniformRand(); <<
468 for (j = 1; j < nEnergies; ++j) { <<
469 itt = j; <<
470 if (randkal*running[nEnergies - 1] < run <<
471 } <<
472 <<
473 // Interpolate the secondary energy <<
474 G4double x, x1, x2, y1, y2; <<
475 if (itt == 0) itt = 1; <<
476 x = randkal * running[nEnergies - 1]; <<
477 x1 = running[itt - 1]; <<
478 x2 = running[itt]; <<
479 G4double compoundFraction; <<
480 // interpolate energy <<
481 y1 = theAngular[itt - 1].GetLabel(); <<
482 y2 = theAngular[itt].GetLabel(); <<
483 fsEnergy = theInt.Interpolate(theManager.G <<
484 <<
485 // For theta, interpolate the compoundFrac <<
486 G4double cLow = theAngular[itt - 1].GetVal <<
487 G4double cHigh = theAngular[itt].GetValue( <<
488 compoundFraction = theInt.Interpolate(theM <<
489 <<
490 if (compoundFraction > 1.0) <<
491 compoundFraction = 1.0; // Protection a <<
492 <<
493 delete[] running; <<
494 <<
495 // get cosTh <<
496 G4double incidentEnergy = anEnergy; <<
497 G4double incidentMass = theProjectile->Get <<
498 G4double productEnergy = fsEnergy; <<
499 G4double productMass = result->GetMass(); <<
500 auto targetZ = G4int(fCache.Get().theTarge <<
501 auto targetA = G4int(fCache.Get().theTarge <<
502 <<
503 // To correspond to natural composition (- <<
504 if (targetA == 0) targetA = G4int(fCache.G <<
505 G4double targetMass = fCache.Get().theTarg <<
506 auto incidentA = G4int(incidentMass / amu_ <<
507 auto incidentZ = G4int(theProjectile->GetP <<
508 G4int residualA = targetA + incidentA - A; <<
509 G4int residualZ = targetZ + incidentZ - Z; <<
510 G4double residualMass = G4NucleiProperties <<
511 <<
512 G4ParticleHPKallbachMannSyst theKallbach( <<
513 compoundFraction, incidentEnergy, incide <<
514 residualA, residualZ, targetMass, target <<
515 cosTh = theKallbach.Sample(anEnergy); <<
516 // end (angularRep == 2) branch <<
517 } <<
518 else if (angularRep > 10 && angularRep < 16) <<
519 G4double random = G4UniformRand(); <<
520 auto running = new G4double[nEnergies]; <<
521 running[0] = 0; <<
522 G4double weighted = 0; <<
523 for (i = 1; i < nEnergies; ++i) { <<
524 if (i != 0) running[i] = running[i - 1]; <<
525 running[i] += theInt.GetBinIntegral(theM <<
526 theA <<
527 theA <<
528 weighted += theInt.GetWeightedBinIntegra <<
529 theManager.GetScheme(i - 1), theAngula <<
530 theAngular[i - 1].GetValue(0), theAngu <<
531 } <<
532 <<
533 // Cache the mean energy in this distribut <<
534 if (nEnergies == 1) <<
535 fCache.Get().currentMeanEnergy = 0.0; <<
536 else <<
537 fCache.Get().currentMeanEnergy = weighte <<
538 <<
539 if (nEnergies == 1) it = 0; <<
540 for (i = 1; i < nEnergies; i++) { <<
541 it = i; <<
542 if (random < running[i] / running[nEnerg <<
543 } <<
544 <<
545 if (it < nDiscreteEnergies || it == 0) { <<
546 if (it == 0) { <<
547 fsEnergy = theAngular[0].GetLabel(); <<
548 G4ParticleHPVector theStore; <<
549 G4int aCounter = 0; <<
550 for (G4int j = 1; j < nAngularParamete <<
551 theStore.SetX(aCounter, theAngular[0 <<
552 theStore.SetY(aCounter, theAngular[0 <<
553 aCounter++; <<
554 } 656 }
555 G4InterpolationManager aMan; <<
556 aMan.Init(angularRep - 10, nAngularPar <<
557 theStore.SetInterpolationManager(aMan) <<
558 cosTh = theStore.Sample(); <<
559 } <<
560 else { <<
561 fsEnergy = theAngular[it].GetLabel(); <<
562 G4ParticleHPVector theStore; 657 G4ParticleHPVector theStore;
563 G4InterpolationManager aMan; << 658 theStore.SetInterpolationManager(aMan); // Store interpolates f(costh)
564 aMan.Init(angularRep - 10, nAngularPar << 659 x1 = y1;
565 theStore.SetInterpolationManager(aMan) << 660 x2 = y2;
566 G4InterpolationScheme currentScheme = << 661 G4double x, y;
567 G4int aCounter = 0; << 662 //for(i=0;i<theBuff1.GetVectorLength(); i++);
568 for (G4int j = 1; j < nAngularParamete << 663 for(i=0;i<theBuff1.GetVectorLength(); i++)
569 theStore.SetX(aCounter, theAngular[i << 664 {
570 theStore.SetY(aCounter, theInt.Inter << 665 x = theBuff1.GetX(i); // costh binning identical
571 << 666 y1 = theBuff1.GetY(i);
572 << 667 y2 = theBuff2.GetY(i);
573 << 668 y = theInt.Interpolate(theManager.GetScheme(it),
574 << 669 fsEnergy, theAngular[it-1].GetLabel(),
575 ++aCounter; << 670 theAngular[it].GetLabel(), y1, y2);
>> 671 theStore.SetX(i, x);
>> 672 theStore.SetY(i, y);
576 } 673 }
577 cosTh = theStore.Sample(); 674 cosTh = theStore.Sample();
578 } 675 }
>> 676 delete [] running;
579 } 677 }
580 else { << 678 else
581 G4double x1 = running[it - 1] / running[ << 679 {
582 G4double x2 = running[it] / running[nEne << 680 throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPContAngularPar::Sample: Unknown angular representation");
583 G4double y1 = theAngular[it - 1].GetLabe <<
584 G4double y2 = theAngular[it].GetLabel(); <<
585 fsEnergy = theInt.Interpolate(theManager <<
586 G4ParticleHPVector theBuff1; <<
587 G4ParticleHPVector theBuff2; <<
588 G4InterpolationManager aMan; <<
589 aMan.Init(angularRep - 10, nAngularParam <<
590 <<
591 G4int j; <<
592 for (i = 0, j = 1; i < nAngularParameter <<
593 theBuff1.SetX(i, theAngular[it - 1].Ge <<
594 theBuff1.SetY(i, theAngular[it - 1].Ge <<
595 theBuff2.SetX(i, theAngular[it].GetVal <<
596 theBuff2.SetY(i, theAngular[it].GetVal <<
597 } <<
598 <<
599 G4ParticleHPVector theStore; <<
600 theStore.SetInterpolationManager(aMan); <<
601 x1 = y1; <<
602 x2 = y2; <<
603 G4double x, y; <<
604 for (i = 0; i < theBuff1.GetVectorLength <<
605 x = theBuff1.GetX(i); // costh binnin <<
606 y1 = theBuff1.GetY(i); <<
607 y2 = theBuff2.GetY(i); <<
608 y = theInt.Interpolate(theManager.GetS <<
609 theAngular[it]. <<
610 theStore.SetX(i, x); <<
611 theStore.SetY(i, y); <<
612 } <<
613 cosTh = theStore.Sample(); <<
614 } 681 }
615 delete[] running; << 682 result->SetKineticEnergy(fsEnergy);
616 } << 683 G4double phi = twopi*G4UniformRand();
617 else { << 684 G4double theta = std::acos(cosTh);
618 throw G4HadronicException(__FILE__, __LINE << 685 G4double sinth = std::sin(theta);
619 "G4ParticleHPCon << 686 G4double mtot = result->GetTotalMomentum();
>> 687 G4ThreeVector tempVector(mtot*sinth*std::cos(phi), mtot*sinth*std::sin(phi), mtot*std::cos(theta) );
>> 688 result->SetMomentum(tempVector);
>> 689 // return the result.
>> 690 return result;
620 } 691 }
621 //G4cout << " Efin=" << fsEnergy << G4endl; <<
622 result->SetKineticEnergy(fsEnergy); <<
623 692
624 G4double phi = twopi * G4UniformRand(); <<
625 if(cosTh > 1.0) { cosTh = 1.0; } <<
626 else if (cosTh < -1.0) { cosTh = -1.0; } <<
627 G4double sinth = std::sqrt((1.0 - cosTh)*(1. <<
628 G4double mtot = result->GetTotalMomentum(); <<
629 G4ThreeVector tempVector(mtot * sinth * std: <<
630 result->SetMomentum(tempVector); <<
631 return result; <<
632 } <<
633 693
634 void G4ParticleHPContAngularPar::PrepareTableI << 694 #define MERGE_NEW
635 { <<
636 // Discrete energies: store own energies in <<
637 // <<
638 // The data files sometimes have identical d <<
639 // which would lead to overwriting the alrea <<
640 // creating a hole in the lookup table. <<
641 // No attempt is made here to correct for th <<
642 // is subtracted from the energy in order to <<
643 <<
644 for (G4int ie = 0; ie < nDiscreteEnergies; i <<
645 // check if energy is already present and <<
646 G4double myE = theAngular[ie].GetLabel(); <<
647 while (theDiscreteEnergiesOwn.find(myE) != <<
648 myE -= 1e-6; <<
649 } <<
650 theDiscreteEnergiesOwn[myE] = ie; <<
651 } <<
652 return; <<
653 } <<
654 695
655 void G4ParticleHPContAngularPar::BuildByInterp << 696 void G4ParticleHPContAngularPar::PrepareTableInterpolation(const G4ParticleHPContAngularPar* angParPrev)
656 <<
657 <<
658 <<
659 { 697 {
660 G4int ie, ie1, ie2, ie1Prev, ie2Prev; <<
661 // Only rebuild the interpolation table if t <<
662 // For several subsequent samplings of final <<
663 // interaction the existing table should be <<
664 if (!fCache.Get().fresh) return; <<
665 <<
666 // Make copies of angpar1 and angpar2. Since <<
667 // it can not be excluded that one of them i <<
668 // potentially the old "this" for creating t <<
669 // memory corruption if the old is not store <<
670 const G4ParticleHPContAngularPar copyAngpar1 <<
671 <<
672 nAngularParameters = copyAngpar1.nAngularPar <<
673 theManager = copyAngpar1.theManager; <<
674 theEnergy = anEnergy; <<
675 theMinEner = DBL_MAX; // min and max will b <<
676 theMaxEner = -DBL_MAX; <<
677 698
678 // The two discrete sets must be merged. A v << 699 //----- Discrete energies: store own energies in a map for faster searching
679 // be copied to the array in the end. Since << 700 G4int ie;
680 // contains pointers, can't simply assign el << 701 for(ie=0; ie<nDiscreteEnergies; ie++) {
681 // needs to call the explicit Set() method i <<
682 <<
683 // First, average probabilities for those li <<
684 const std::map<G4double, G4int> discEnerOwn1 <<
685 const std::map<G4double, G4int> discEnerOwn2 <<
686 std::map<G4double, G4int>::const_iterator it <<
687 std::map<G4double, G4int>::const_iterator it <<
688 std::vector<G4ParticleHPList*> vAngular(disc <<
689 G4double discEner1; <<
690 for (itedeo1 = discEnerOwn1.cbegin(); itedeo <<
691 discEner1 = itedeo1->first; <<
692 if (discEner1 < theMinEner) { <<
693 theMinEner = discEner1; <<
694 } <<
695 if (discEner1 > theMaxEner) { <<
696 theMaxEner = discEner1; <<
697 } <<
698 ie1 = itedeo1->second; <<
699 itedeo2 = discEnerOwn2.find(discEner1); <<
700 if (itedeo2 == discEnerOwn2.cend()) { <<
701 ie2 = -1; <<
702 } <<
703 else { <<
704 ie2 = itedeo2->second; <<
705 } <<
706 vAngular[ie1] = new G4ParticleHPList(); <<
707 vAngular[ie1]->SetLabel(copyAngpar1.theAng <<
708 G4double val1, val2; <<
709 for (G4int ip = 0; ip < nAngularParameters <<
710 val1 = copyAngpar1.theAngular[ie1].GetVa <<
711 if (ie2 != -1) { <<
712 val2 = copyAngpar2.theAngular[ie2].Get <<
713 } <<
714 else { <<
715 val2 = 0.; <<
716 } <<
717 G4double value = theInt.Interpolate(aSch <<
718 copy <<
719 vAngular[ie1]->SetValue(ip, value); <<
720 } <<
721 } // itedeo1 loop <<
722 <<
723 // Add the ones in set2 but not in set1 <<
724 std::vector<G4ParticleHPList*>::const_iterat <<
725 G4double discEner2; <<
726 for (itedeo2 = discEnerOwn2.cbegin(); itedeo <<
727 discEner2 = itedeo2->first; <<
728 ie2 = itedeo2->second; <<
729 G4bool notFound = true; <<
730 itedeo1 = discEnerOwn1.find(discEner2); <<
731 if (itedeo1 != discEnerOwn1.cend()) { <<
732 notFound = false; <<
733 } <<
734 if (notFound) { <<
735 // not yet in list <<
736 if (discEner2 < theMinEner) { <<
737 theMinEner = discEner2; <<
738 } <<
739 if (discEner2 > theMaxEner) { <<
740 theMaxEner = discEner2; <<
741 } <<
742 // find position to insert <<
743 G4bool isInserted = false; <<
744 ie = 0; <<
745 for (itv = vAngular.cbegin(); itv != vAn <<
746 if (discEner2 > (*itv)->GetLabel()) { <<
747 itv = vAngular.insert(itv, new G4Par <<
748 (*itv)->SetLabel(copyAngpar2.theAngu <<
749 isInserted = true; <<
750 break; <<
751 } <<
752 } <<
753 if (!isInserted) { <<
754 ie = (G4int)vAngular.size(); <<
755 vAngular.push_back(new G4ParticleHPLis <<
756 vAngular[ie]->SetLabel(copyAngpar2.the <<
757 isInserted = true; <<
758 } <<
759 <<
760 G4double val1, val2; <<
761 for (G4int ip = 0; ip < nAngularParamete <<
762 val1 = 0; <<
763 val2 = copyAngpar2.theAngular[ie2].Get <<
764 G4double value = theInt.Interpolate(aS <<
765 co <<
766 vAngular[ie]->SetValue(ip, value); <<
767 } <<
768 } // end if(notFound) <<
769 } // end loop on itedeo2 <<
770 <<
771 // Store new discrete list <<
772 nDiscreteEnergies = (G4int)vAngular.size(); <<
773 delete[] theAngular; <<
774 theAngular = nullptr; <<
775 if (nDiscreteEnergies > 0) { <<
776 theAngular = new G4ParticleHPList[nDiscret <<
777 } <<
778 theDiscreteEnergiesOwn.clear(); <<
779 theDiscreteEnergies.clear(); <<
780 for (ie = 0; ie < nDiscreteEnergies; ++ie) { <<
781 theAngular[ie].SetLabel(vAngular[ie]->GetL <<
782 for (G4int ip = 0; ip < nAngularParameters <<
783 theAngular[ie].SetValue(ip, vAngular[ie] <<
784 } <<
785 theDiscreteEnergiesOwn[theAngular[ie].GetL 702 theDiscreteEnergiesOwn[theAngular[ie].GetLabel()] = ie;
786 theDiscreteEnergies.insert(theAngular[ie]. <<
787 } 703 }
>> 704 if( !angParPrev ) return;
788 705
789 // The continuous energies need to be made f << 706 //----- Discrete energies: use energies that appear in one or another
790 // ones. Therefore the re-assignemnt of theA << 707 for(ie=0; ie<nDiscreteEnergies; ie++) {
791 // after the continuous energy set is also f << 708 theDiscreteEnergies.insert(theAngular[ie].GetLabel());
792 // total number of nEnergies is known and th <<
793 <<
794 // Get minimum and maximum energy interpolat <<
795 // Don't use theMinEner or theMaxEner here, <<
796 // need the interpolated range from the orig <<
797 G4double interMinEner = copyAngpar1.GetMinEn <<
798 + (theEnergy - copyA <<
799 * (copyAngpar2.G <<
800 / (copyAngpar2.G <<
801 G4double interMaxEner = copyAngpar1.GetMaxEn <<
802 + (theEnergy - copyA <<
803 * (copyAngpar2.G <<
804 / (copyAngpar2.G <<
805 <<
806 // Loop to energies of new set <<
807 theEnergiesTransformed.clear(); <<
808 <<
809 G4int nEnergies1 = copyAngpar1.GetNEnergies( <<
810 G4int nDiscreteEnergies1 = copyAngpar1.GetND <<
811 G4double minEner1 = copyAngpar1.GetMinEner() <<
812 G4double maxEner1 = copyAngpar1.GetMaxEner() <<
813 G4int nEnergies2 = copyAngpar2.GetNEnergies( <<
814 G4int nDiscreteEnergies2 = copyAngpar2.GetND <<
815 G4double minEner2 = copyAngpar2.GetMinEner() <<
816 G4double maxEner2 = copyAngpar2.GetMaxEner() <<
817 <<
818 // First build the list of transformed energ <<
819 // to the new min max by assuming that the m <<
820 // each set would be scalable to the new, in <<
821 // max range <<
822 <<
823 G4double e1(0.); <<
824 G4double eTNorm1(0.); <<
825 for (ie1 = nDiscreteEnergies1; ie1 < nEnergi <<
826 e1 = copyAngpar1.theAngular[ie1].GetLabel( <<
827 eTNorm1 = (e1 - minEner1); <<
828 if (maxEner1 != minEner1) eTNorm1 /= (maxE <<
829 if (eTNorm1 >= 0 && eTNorm1 <= 1) theEnerg <<
830 } 709 }
831 << 710 G4int nDiscreteEnergiesPrev = angParPrev->GetNDiscreteEnergies();
832 G4double e2(0.); << 711 for(ie=0; ie<nDiscreteEnergiesPrev; ie++) {
833 G4double eTNorm2(0.); << 712 theDiscreteEnergies.insert(angParPrev->theAngular[ie].GetLabel());
834 for (ie2 = nDiscreteEnergies2; ie2 < nEnergi << 713 }
835 e2 = copyAngpar2.theAngular[ie2].GetLabel( << 714
836 eTNorm2 = (e2 - minEner2); << 715 //--- Get the values for which interpolation will be done : all energies of this and previous ContAngularPar
837 if (maxEner2 != minEner2) eTNorm2 /= (maxE << 716 for(ie=nDiscreteEnergies; ie<nEnergies; ie++) {
838 if (eTNorm2 >= 0 && eTNorm2 <= 1) theEnerg << 717 G4double ener = theAngular[ie].GetLabel();
>> 718 G4double enerT = (ener-theMinEner)/(theMaxEner-theMinEner);
>> 719 theEnergiesTransformed.insert(enerT);
>> 720 //- if( getenv("G4PHPTEST2") ) G4cout <<this << " G4ParticleHPContAngularPar::PrepareTableInterpolation theEnergiesTransformed1 " << enerT << G4endl; //GDEB
>> 721 }
>> 722 G4int nEnergiesPrev = angParPrev->GetNEnergies();
>> 723 G4double minEnerPrev = angParPrev->GetMinEner();
>> 724 G4double maxEnerPrev = angParPrev->GetMaxEner();
>> 725 for(ie=nDiscreteEnergiesPrev; ie<nEnergiesPrev; ie++) {
>> 726 G4double ener = angParPrev->theAngular[ie].GetLabel();
>> 727 G4double enerT = (ener-minEnerPrev)/(maxEnerPrev-minEnerPrev);
>> 728 theEnergiesTransformed.insert(enerT);
>> 729 //- if( getenv("G4PHPTEST2") ) G4cout << this << " G4ParticleHPContAngularPar::PrepareTableInterpolation theEnergiesTransformed2 " << enerT << G4endl; //GDEB
839 } 730 }
>> 731 // add the maximum energy
>> 732 theEnergiesTransformed.insert(1.);
840 733
841 // Now the list of energies is complete << 734 }
842 nEnergies = nDiscreteEnergies + (G4int)theEn <<
843 735
844 // Create final array of angular parameters << 736 void G4ParticleHPContAngularPar::BuildByInterpolation(G4double anEnergy, G4InterpolationScheme aScheme,
845 const std::size_t esize = nEnergies > 0 ? nE << 737 G4ParticleHPContAngularPar & angpar1,
846 auto theNewAngular = new G4ParticleHPList[es << 738 G4ParticleHPContAngularPar & angpar2)
847 << 739 {
848 // Copy discrete energies and interpolated p << 740 G4int ie,ie1,ie2, ie1Prev, ie2Prev;
849 << 741 nAngularParameters = angpar1.nAngularParameters;
850 if (theAngular != nullptr) { << 742 theManager = angpar1.theManager;
851 for (ie = 0; ie < nDiscreteEnergies; ++ie) << 743 theEnergy = anEnergy;
852 theNewAngular[ie].SetLabel(theAngular[ie <<
853 for (G4int ip = 0; ip < nAngularParamete <<
854 theNewAngular[ie].SetValue(ip, theAngu <<
855 } <<
856 } <<
857 delete[] theAngular; <<
858 } <<
859 theAngular = theNewAngular; <<
860 744
861 // Interpolate the continuous energies for n << 745 nDiscreteEnergies = theDiscreteEnergies.size();
862 auto iteet = theEnergiesTransformed.begin(); << 746 std::set<G4double>::const_iterator itede;
>> 747 std::map<G4double,G4int> discEnerOwn1 = angpar1.GetDiscreteEnergiesOwn();
>> 748 std::map<G4double,G4int> discEnerOwn2 = angpar2.GetDiscreteEnergiesOwn();
>> 749 std::map<G4double,G4int>::const_iterator itedeo;
>> 750 ie = 0;
>> 751 for( itede = theDiscreteEnergies.begin(); itede != theDiscreteEnergies.end(); itede++, ie++ ) {
>> 752 G4double discEner = *itede;
>> 753 itedeo = discEnerOwn1.find(discEner);
>> 754 if( itedeo == discEnerOwn1.end() ) {
>> 755 ie1 = 0;
>> 756 } else {
>> 757 ie1 = -1;
>> 758 }
>> 759 itedeo = discEnerOwn2.find(discEner);
>> 760 if( itedeo == discEnerOwn2.end() ) {
>> 761 ie2 = 0;
>> 762 } else {
>> 763 ie2 = -1;
>> 764 }
863 765
864 G4double e1Interp(0.); << 766 theAngular[ie].SetLabel(discEner);
865 G4double e2Interp(0.); << 767 G4double val1, val2;
866 for (ie = nDiscreteEnergies; ie < nEnergies; << 768 for(G4int ip=0; ip<nAngularParameters; ip++) {
>> 769 if( ie1 != -1 ) {
>> 770 val1 = angpar1.theAngular[ie1].GetValue(ip);
>> 771 } else {
>> 772 val1 = 0.;
>> 773 }
>> 774 if( ie2 != -1 ) {
>> 775 val2 = angpar2.theAngular[ie2].GetValue(ip);
>> 776 } else {
>> 777 val2 = 0.;
>> 778 }
>> 779
>> 780 G4double value = theInt.Interpolate(aScheme, anEnergy,
>> 781 angpar1.theEnergy, angpar2.theEnergy,
>> 782 val1,
>> 783 val2);
>> 784 if( std::getenv("G4PHPTEST2") ) G4cout << ie << " " << ip << " G4ParticleHPContAngularPar::Merge DiscreteEnergies val1 " << val1 << " val2 " << val2 << " value " << value << G4endl; //GDEB
>> 785
>> 786 theAngular[ie].SetValue(ip, value);
>> 787 }
>> 788 }
>> 789
>> 790 if(theAngular != 0) delete [] theAngular;
>> 791 nEnergies = nDiscreteEnergies + angpar2.GetNEnergiesTransformed();
>> 792 theAngular = new G4ParticleHPList [nEnergies];
>> 793
>> 794 //---- Get minimum and maximum energy interpolating
>> 795 theMinEner = angpar1.GetMinEner() + (theEnergy-angpar1.GetEnergy()) * (angpar2.GetMinEner()-angpar1.GetMinEner())/(angpar2.GetEnergy()-angpar1.GetEnergy());
>> 796 theMaxEner = angpar1.GetMaxEner() + (theEnergy-angpar1.GetEnergy()) * (angpar2.GetMaxEner()-angpar1.GetMaxEner())/(angpar2.GetEnergy()-angpar1.GetEnergy());
>> 797
>> 798 if( std::getenv("G4PHPTEST2") ) G4cout << " G4ParticleHPContAngularPar::Merge E " << anEnergy << " minmax " << theMinEner << " " << theMaxEner << G4endl; //GDEB
>> 799
>> 800 //--- Loop to energies of new set
>> 801 std::set<G4double> energiesTransformed = angpar2.GetEnergiesTransformed();
>> 802 std::set<G4double>::const_iterator iteet = energiesTransformed.begin();
>> 803 G4int nEnergies1 = angpar1.GetNEnergies();
>> 804 G4int nDiscreteEnergies1 = angpar1.GetNDiscreteEnergies();
>> 805 G4double minEner1 = angpar1.GetMinEner();
>> 806 G4double maxEner1 = angpar1.GetMaxEner();
>> 807 G4int nEnergies2 = angpar2.GetNEnergies();
>> 808 G4int nDiscreteEnergies2 = angpar2.GetNDiscreteEnergies();
>> 809 G4double minEner2 = angpar2.GetMinEner();
>> 810 G4double maxEner2 = angpar2.GetMaxEner();
>> 811 for(ie=nDiscreteEnergies; ie<nEnergies; ie++,iteet++) {
867 G4double eT = (*iteet); 812 G4double eT = (*iteet);
868 813
869 //--- Use eT1 = eT: Get energy and paramet << 814 //--- Use eT1 = eT: Get energy and parameters of angpar1 for this eT
870 e1Interp = (maxEner1 - minEner1) * eT + mi << 815 G4double e1 = (maxEner1-minEner1) * eT + minEner1;
871 //----- Get parameter value corresponding << 816 //----- Get parameter value corresponding to this e1
872 for (ie1 = nDiscreteEnergies1; ie1 < nEner << 817 for(ie1=nDiscreteEnergies1; ie1<nEnergies1; ie1++) {
873 if ((copyAngpar1.theAngular[ie1].GetLabe << 818 if( (angpar1.theAngular[ie1].GetLabel() - e1) > 1.E-10*e1 ) break;
874 } 819 }
875 ie1Prev = ie1 - 1; 820 ie1Prev = ie1 - 1;
876 if (ie1 == 0) ++ie1Prev; << 821 if( ie1 == 0 ) ie1Prev++;
877 if (ie1 == nEnergies1) { << 822 if( ie1 == nEnergies1 ) {
878 ie1--; 823 ie1--;
879 ie1Prev = ie1; 824 ie1Prev = ie1;
880 } 825 }
881 << 826 //--- Use eT2 = eT: Get energy and parameters of angpar2 for this eT
882 //--- Use eT2 = eT: Get energy and paramet << 827 G4double e2 = (maxEner2-minEner2) * eT + minEner2;
883 e2Interp = (maxEner2 - minEner2) * eT + mi << 828 //----- Get parameter value corresponding to this e2
884 //----- Get parameter value corresponding << 829 for(ie2=nDiscreteEnergies2; ie2<nEnergies2; ie2++) {
885 for (ie2 = nDiscreteEnergies2; ie2 < nEner << 830 // G4cout << " GET IE2 " << ie2 << " - " << angpar2.theAngular[ie2].GetLabel() - e2 << " " << angpar2.theAngular[ie2].GetLabel() << " " << e2 <<G4endl;
886 if ((copyAngpar2.theAngular[ie2].GetLabe << 831 if( (angpar2.theAngular[ie2].GetLabel() - e2) > 1.E-10*e2 ) break;
887 } 832 }
888 ie2Prev = ie2 - 1; 833 ie2Prev = ie2 - 1;
889 if (ie2 == 0) ++ie2Prev; << 834 if( ie2 == 0 ) ie2Prev++;
890 if (ie2 == nEnergies2) { << 835 if( ie2 == nEnergies2 ) {
891 ie2--; 836 ie2--;
892 ie2Prev = ie2; 837 ie2Prev = ie2;
893 } 838 }
894 839
895 //---- Energy corresponding to energy tran << 840 //---- Energy corresponding to energy transformed
896 G4double eN = (interMaxEner - interMinEner << 841 G4double eN = (theMaxEner-theMinEner) * eT + theMinEner;
897 << 842 if( std::getenv("G4PHPTEST2") ) G4cout << ie << " " << ie1 << " " << ie2 << " G4ParticleHPContAngularPar::loop eT " << eT << " -> eN " << eN << " e1 " << e1 << " e2 " << e2 << G4endl; //GDEB
>> 843
898 theAngular[ie].SetLabel(eN); 844 theAngular[ie].SetLabel(eN);
899 if (eN < theMinEner) { << 845
900 theMinEner = eN; << 846 for(G4int ip=0; ip<nAngularParameters; ip++) {
901 } << 847 G4double val1 = theInt.Interpolate2(theManager.GetScheme(ie),
902 if (eN > theMaxEner) { << 848 e1,
903 theMaxEner = eN; << 849 angpar1.theAngular[ie1Prev].GetLabel(),
904 } << 850 angpar1.theAngular[ie1].GetLabel(),
905 << 851 angpar1.theAngular[ie1Prev].GetValue(ip),
906 G4double val1(0.); << 852 angpar1.theAngular[ie1].GetValue(ip)) * (maxEner1-minEner1);
907 G4double val2(0.); << 853 G4double val2 = theInt.Interpolate2(theManager.GetScheme(ie),
908 G4double value(0.); << 854 e2,
909 for (G4int ip = 0; ip < nAngularParameters << 855 angpar2.theAngular[ie2Prev].GetLabel(),
910 val1 = theInt.Interpolate2( << 856 angpar2.theAngular[ie2].GetLabel(),
911 theManager.GetScheme(ie), e1Int << 857 angpar2.theAngular[ie2Prev].GetValue(ip),
912 copyAngpar1.theAngular[ie1].Get << 858 angpar2.theAngular[ie2].GetValue(ip)) * (maxEner2-minEner2);
913 copyAngpar1.theAngular[ie1].Get << 859
914 * (maxEner1 - minEner1); << 860 G4double value = theInt.Interpolate(aScheme, anEnergy,
915 val2 = theInt.Interpolate2( << 861 angpar1.theEnergy, angpar2.theEnergy,
916 theManager.GetScheme(ie), e2Int << 862 val1,
917 copyAngpar2.theAngular[ie2].Get << 863 val2);
918 copyAngpar2.theAngular[ie2].Get << 864 //value /= (theMaxEner-theMinEner);
919 * (maxEner2 - minEner2); << 865 if ( theMaxEner != theMinEner ) {
920 << 866 value /= (theMaxEner-theMinEner);
921 value = theInt.Interpolate(aScheme, anEn << 867 } else if ( value != 0 ) {
922 val1, val2); << 868 throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPContAngularPar::PrepareTableInterpolation theMaxEner == theMinEner and value != 0.");
923 if (interMaxEner != interMinEner) { << 869 }
924 value /= (interMaxEner - interMinEner) << 870 if( std::getenv("G4PHPTEST2") ) G4cout << ie << " " << ip << " G4ParticleHPContAngularPar::Merge val1 " << val1 << " val2 " << val2 << " value " << value << G4endl; //GDEB
925 } << 871 //- val1 = angpar1.theAngular[ie1-1].GetValue(ip) * (maxEner1-minEner1);
926 else if (value != 0) { << 872 //- val2 = angpar2.theAngular[ie2-1].GetValue(ip) * (maxEner2-minEner2);
927 throw G4HadronicException(__FILE__, __ << 873 //- if( getenv("G4PHPTEST2") ) G4cout << ie << " " << ip << " G4ParticleHPContAngularPar::MergeOLD val1 " << val1 << " val2 " << val2 << " value " << value << G4endl; //GDEB
928 "G4ParticleH << 874
929 "interMaxEne <<
930 } <<
931 theAngular[ie].SetValue(ip, value); 875 theAngular[ie].SetValue(ip, value);
932 } 876 }
933 } // end loop on nDiscreteEnergies << 877 }
934 878
935 for (itv = vAngular.cbegin(); itv != vAngula << 879 if( std::getenv("G4PHPTEST2") ) {
936 delete (*itv); << 880 G4cout << " G4ParticleHPContAngularPar::Merge ANGPAR1 " << G4endl; //GDEB
>> 881 angpar1.Dump();
>> 882 G4cout << " G4ParticleHPContAngularPar::Merge ANGPAR2 " << G4endl;
>> 883 angpar2.Dump();
>> 884 G4cout << " G4ParticleHPContAngularPar::Merge ANGPARNEW " << G4endl;
>> 885 Dump();
>> 886 }
937 } 887 }
938 888
939 void G4ParticleHPContAngularPar::Dump() const << 889 void G4ParticleHPContAngularPar::Dump()
940 { 890 {
941 G4cout << theEnergy << " " << nEnergies << " << 891 G4cout << theEnergy << " " << nEnergies << " " << nDiscreteEnergies << " " << nAngularParameters << G4endl;
942 << G4endl; <<
943 892
944 for (G4int ii = 0; ii < nEnergies; ++ii) << 893 for(G4int ii=0; ii<nEnergies; ii++) {
945 theAngular[ii].Dump(); 894 theAngular[ii].Dump();
>> 895 }
>> 896
946 } 897 }
947 898