| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // neutron_hp -- source file
27 // J.P. Wellisch, Nov-1996
28 // A prototype of the low energy neutron transport model.
29 //
30 // 09-May-06 fix in Sample by T. Koi
31 // 080318 Fix Compilation warnings - gcc-4.3.0 by T. Koi
32 // (This fix has a real effect to the code.)
33 // 080409 Fix div0 error with G4FPE by T. Koi
34 // 080612 Fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #1
35 // 080714 Limiting the sum of energy of secondary particles by T. Koi
36 // 080801 Fix div0 error wiht G4FPE and memory leak by T. Koi
37 // 081024 G4NucleiPropertiesTable:: to G4NucleiProperties::
38 //
39 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
40 //
41 // June-2019 - E. Mendoza --> redefinition of the residual mass to consider incident particles different than neutrons.
42
43 #include "G4ParticleHPContAngularPar.hh"
44 #include "G4PhysicalConstants.hh"
45 #include "G4SystemOfUnits.hh"
46 #include "G4ParticleHPLegendreStore.hh"
47 #include "G4Gamma.hh"
48 #include "G4Electron.hh"
49 #include "G4Positron.hh"
50 #include "G4Neutron.hh"
51 #include "G4Proton.hh"
52 #include "G4Deuteron.hh"
53 #include "G4Triton.hh"
54 #include "G4He3.hh"
55 #include "G4Alpha.hh"
56 #include "G4ParticleHPVector.hh"
57 #include "G4NucleiProperties.hh"
58 #include "G4ParticleHPKallbachMannSyst.hh"
59 #include "G4IonTable.hh"
60 #include <set>
61
62 G4ParticleHPContAngularPar::G4ParticleHPContAngularPar( G4ParticleDefinition* projectile )
63 {
64 theAngular = 0;
65 if ( fCache.Get() == 0 ) cacheInit();
66 fCache.Get()->currentMeanEnergy = -2;
67 fCache.Get()->fresh = true;
68 adjustResult = true;
69 if ( std::getenv( "G4PHP_DO_NOT_ADJUST_FINAL_STATE" ) ) adjustResult = false;
70
71 theMinEner = DBL_MAX;
72 theMaxEner = -DBL_MAX;
73 theProjectile = projectile;
74
75 theEnergy = 0.0;
76 nEnergies = 0;
77 nDiscreteEnergies = 0;
78 nAngularParameters = 0;
79 }
80
81 void G4ParticleHPContAngularPar::Init(std::istream & aDataFile, G4ParticleDefinition* projectile)
82 {
83 adjustResult = true;
84 if ( std::getenv( "G4PHP_DO_NOT_ADJUST_FINAL_STATE" ) ) adjustResult = false;
85
86 theProjectile = projectile;
87
88 aDataFile >> theEnergy >> nEnergies >> nDiscreteEnergies >> nAngularParameters;
89 /*if( std::getenv("G4PHPTEST") )*/
90 theEnergy *= eV;
91 theAngular = new G4ParticleHPList [nEnergies];
92 for(G4int i=0; i<nEnergies; i++)
93 {
94 G4double sEnergy;
95 aDataFile >> sEnergy;
96 sEnergy*=eV;
97 theAngular[i].SetLabel(sEnergy);
98 theAngular[i].Init(aDataFile, nAngularParameters, 1.);
99 theMinEner = std::min(theMinEner,sEnergy);
100 theMaxEner = std::max(theMaxEner,sEnergy);
101 }
102 }
103
104 G4ReactionProduct *
105 G4ParticleHPContAngularPar::Sample(G4double anEnergy, G4double massCode, G4double /*targetMass*/,
106 G4int angularRep, G4int /*interpolE*/ )
107 {
108 if( std::getenv("G4PHPTEST") ) G4cout << " G4ParticleHPContAngularPar::Sample " << anEnergy << " " << massCode << " " << angularRep << G4endl; //GDEB
109 if ( fCache.Get() == 0 ) cacheInit();
110 G4ReactionProduct * result = new G4ReactionProduct;
111 G4int Z = static_cast<G4int>(massCode/1000);
112 G4int A = static_cast<G4int>(massCode-1000*Z);
113 if(massCode==0)
114 {
115 result->SetDefinition(G4Gamma::Gamma());
116 }
117 else if(A==0)
118 {
119 result->SetDefinition(G4Electron::Electron());
120 if(Z==1) result->SetDefinition(G4Positron::Positron());
121 }
122 else if(A==1)
123 {
124 result->SetDefinition(G4Neutron::Neutron());
125 if(Z==1) result->SetDefinition(G4Proton::Proton());
126 }
127 else if(A==2)
128 {
129 result->SetDefinition(G4Deuteron::Deuteron());
130 }
131 else if(A==3)
132 {
133 result->SetDefinition(G4Triton::Triton());
134 if(Z==2) result->SetDefinition(G4He3::He3());
135 }
136 else if(A==4)
137 {
138 result->SetDefinition(G4Alpha::Alpha());
139 if(Z!=2) throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPContAngularPar: Unknown ion case 1");
140 }
141 else
142 {
143 result->SetDefinition(G4IonTable::GetIonTable()->GetIon(Z,A,0));
144 }
145 G4int i(0);
146 G4int it(0);
147 G4double fsEnergy(0);
148 G4double cosTh(0);
149
150 if( angularRep == 1 )
151 {
152 // 080612 Fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #1
153 //if (interpolE == 2)
154 //110609 above was wrong interupition, pointed out by E.Mendoza and D.Cano (CIMAT)
155 //Following are reviesd version written by T.Koi (SLAC)
156 if ( nDiscreteEnergies != 0 )
157 {
158
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();
189
190 G4double * running = new G4double[nEnergies+1];
191 running[0] = 0.0;
192
193 for ( G4int j = 0 ; j < nDiscreteEnergies ; j++ )
194 {
195 G4double delta = 0.0;
196 if ( theAngular[j].GetLabel() <= fCache.Get()->remaining_energy ) delta = theAngular[i].GetValue(0);
197 running[j+1] = running[j] + delta;
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 }
249 }
250 */
251 // Normalize random
252 random *= (tot_prob_DIS + tot_prob_CON);
253 //2nd Judge Discrete or not This shoudl be relatively close to 1 For safty
254 if ( random <= ( tot_prob_DIS / ( tot_prob_DIS + tot_prob_CON ) ) || nDiscreteEnergies == nEnergies )
255 {
256 // Discrete Emission
257 for ( G4int j = 0 ; j < nDiscreteEnergies ; j++ )
258 {
259 //Here we should use i+1
260 if ( random < running[ j+1 ] )
261 {
262 it = j;
263 break;
264 }
265 }
266 fsEnergy = theAngular[ it ].GetLabel();
267
268 G4ParticleHPLegendreStore theStore(1);
269 theStore.Init(0,fsEnergy,nAngularParameters);
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
334 {
335 // Only continue, TK will clean up
336
337 //080714
338 if ( fCache.Get()->fresh == true )
339 {
340 fCache.Get()->remaining_energy = theAngular[ nEnergies-1 ].GetLabel();
341 fCache.Get()->fresh = false;
342 }
343 //080714
344 G4double random = G4UniformRand();
345 G4double * running = new G4double[nEnergies];
346 running[0]=0;
347 G4double weighted = 0;
348 for(i=1; i<nEnergies; i++)
349 {
350 /*
351 if(i!=0)
352 {
353 running[i]=running[i-1];
354 }
355 running[i] += theInt.GetBinIntegral(theManager.GetScheme(i-1),
356 theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
357 theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
358 weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(i-1),
359 theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
360 theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
361 */
362
363 running[i]=running[i-1];
364 if ( fCache.Get()->remaining_energy >= theAngular[i].GetLabel() )
365 {
366 running[i] += theInt.GetBinIntegral(theManager.GetScheme(i-1),
367 theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
368 theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
369 weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(i-1),
370 theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
371 theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
372 }
373 }
374 // cash the mean energy in this distribution
375 //080409 TKDB
376 if ( nEnergies == 1 || running[nEnergies-1] == 0 )
377 fCache.Get()->currentMeanEnergy = 0.0;
378 else
379 {
380 fCache.Get()->currentMeanEnergy = weighted/running[nEnergies-1];
381 }
382
383 //080409 TKDB
384 if ( nEnergies == 1 ) it = 0;
385
386 //080729
387 if ( running[nEnergies-1] != 0 )
388 {
389 for ( i = 1 ; i < nEnergies ; i++ )
390 {
391 it = i;
392 if ( random < running [ i ] / running [ nEnergies-1 ] ) break;
393 }
394 }
395
396 //080714
397 if ( running [ nEnergies-1 ] == 0 ) it = 0;
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;
501 }
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 {
554 G4double random = G4UniformRand();
555 G4double * running = new G4double[nEnergies];
556 running[0]=0;
557 G4double weighted = 0;
558 for(i=1; i<nEnergies; i++)
559 {
560 if(i!=0) running[i]=running[i-1];
561 running[i] += theInt.GetBinIntegral(theManager.GetScheme(i-1),
562 theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
563 theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
564 weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(i-1),
565 theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
566 theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
567 }
568 // cash the mean energy in this distribution
569 //currentMeanEnergy = weighted/running[nEnergies-1];
570 if ( nEnergies == 1 )
571 fCache.Get()->currentMeanEnergy = 0.0;
572 else
573 fCache.Get()->currentMeanEnergy = weighted/running[nEnergies-1];
574
575 //080409 TKDB
576 if ( nEnergies == 1 ) it = 0;
577 //for(i=0; i<nEnergies; i++)
578 for(i=1; i<nEnergies; i++)
579 {
580 it = i;
581 if(random<running[i]/running[nEnergies-1]) break;
582 }
583 if(it<nDiscreteEnergies||it==0)
584 {
585 if(it==0)
586 {
587 fsEnergy = theAngular[0].GetLabel();
588 G4ParticleHPVector theStore;
589 G4int aCounter = 0;
590 for(G4int j=1; j<nAngularParameters; j+=2)
591 {
592 theStore.SetX(aCounter, theAngular[0].GetValue(j));
593 theStore.SetY(aCounter, theAngular[0].GetValue(j+1));
594 aCounter++;
595 }
596 G4InterpolationManager aMan;
597 aMan.Init(angularRep-10, nAngularParameters-1);
598 theStore.SetInterpolationManager(aMan);
599 cosTh = theStore.Sample();
600 }
601 else
602 {
603 fsEnergy = theAngular[it].GetLabel();
604 G4ParticleHPVector theStore;
605 G4InterpolationManager aMan;
606 aMan.Init(angularRep-10, nAngularParameters-1);
607 theStore.SetInterpolationManager(aMan); // Store interpolates f(costh)
608 G4InterpolationScheme currentScheme = theManager.GetInverseScheme(it);
609 G4int aCounter = 0;
610 for(G4int j=1; j<nAngularParameters; j+=2)
611 {
612 theStore.SetX(aCounter, theAngular[it].GetValue(j));
613 theStore.SetY(aCounter, theInt.Interpolate(currentScheme,
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++;
620 }
621 cosTh = theStore.Sample();
622 }
623 }
624 else
625 {
626 G4double x1 = running[it-1]/running[nEnergies-1];
627 G4double x2 = running[it]/running[nEnergies-1];
628 G4double y1 = theAngular[it-1].GetLabel();
629 G4double y2 = theAngular[it].GetLabel();
630 fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it),
631 random,x1,x2,y1,y2);
632 G4ParticleHPVector theBuff1;
633 G4ParticleHPVector theBuff2;
634 G4InterpolationManager aMan;
635 aMan.Init(angularRep-10, nAngularParameters-1);
636 // theBuff1.SetInterpolationManager(aMan); // Store interpolates f(costh)
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));
655 }
656 }
657 G4ParticleHPVector theStore;
658 theStore.SetInterpolationManager(aMan); // Store interpolates f(costh)
659 x1 = y1;
660 x2 = y2;
661 G4double x, y;
662 //for(i=0;i<theBuff1.GetVectorLength(); i++);
663 for(i=0;i<theBuff1.GetVectorLength(); i++)
664 {
665 x = theBuff1.GetX(i); // costh binning identical
666 y1 = theBuff1.GetY(i);
667 y2 = theBuff2.GetY(i);
668 y = theInt.Interpolate(theManager.GetScheme(it),
669 fsEnergy, theAngular[it-1].GetLabel(),
670 theAngular[it].GetLabel(), y1, y2);
671 theStore.SetX(i, x);
672 theStore.SetY(i, y);
673 }
674 cosTh = theStore.Sample();
675 }
676 delete [] running;
677 }
678 else
679 {
680 throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPContAngularPar::Sample: Unknown angular representation");
681 }
682 result->SetKineticEnergy(fsEnergy);
683 G4double phi = twopi*G4UniformRand();
684 G4double theta = std::acos(cosTh);
685 G4double sinth = std::sin(theta);
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;
691 }
692
693
694 #define MERGE_NEW
695
696 void G4ParticleHPContAngularPar::PrepareTableInterpolation(const G4ParticleHPContAngularPar* angParPrev)
697 {
698
699 //----- Discrete energies: store own energies in a map for faster searching
700 G4int ie;
701 for(ie=0; ie<nDiscreteEnergies; ie++) {
702 theDiscreteEnergiesOwn[theAngular[ie].GetLabel()] = ie;
703 }
704 if( !angParPrev ) return;
705
706 //----- Discrete energies: use energies that appear in one or another
707 for(ie=0; ie<nDiscreteEnergies; ie++) {
708 theDiscreteEnergies.insert(theAngular[ie].GetLabel());
709 }
710 G4int nDiscreteEnergiesPrev = angParPrev->GetNDiscreteEnergies();
711 for(ie=0; ie<nDiscreteEnergiesPrev; ie++) {
712 theDiscreteEnergies.insert(angParPrev->theAngular[ie].GetLabel());
713 }
714
715 //--- Get the values for which interpolation will be done : all energies of this and previous ContAngularPar
716 for(ie=nDiscreteEnergies; ie<nEnergies; ie++) {
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
730 }
731 // add the maximum energy
732 theEnergiesTransformed.insert(1.);
733
734 }
735
736 void G4ParticleHPContAngularPar::BuildByInterpolation(G4double anEnergy, G4InterpolationScheme aScheme,
737 G4ParticleHPContAngularPar & angpar1,
738 G4ParticleHPContAngularPar & angpar2)
739 {
740 G4int ie,ie1,ie2, ie1Prev, ie2Prev;
741 nAngularParameters = angpar1.nAngularParameters;
742 theManager = angpar1.theManager;
743 theEnergy = anEnergy;
744
745 nDiscreteEnergies = theDiscreteEnergies.size();
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 }
765
766 theAngular[ie].SetLabel(discEner);
767 G4double val1, val2;
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++) {
812 G4double eT = (*iteet);
813
814 //--- Use eT1 = eT: Get energy and parameters of angpar1 for this eT
815 G4double e1 = (maxEner1-minEner1) * eT + minEner1;
816 //----- Get parameter value corresponding to this e1
817 for(ie1=nDiscreteEnergies1; ie1<nEnergies1; ie1++) {
818 if( (angpar1.theAngular[ie1].GetLabel() - e1) > 1.E-10*e1 ) break;
819 }
820 ie1Prev = ie1 - 1;
821 if( ie1 == 0 ) ie1Prev++;
822 if( ie1 == nEnergies1 ) {
823 ie1--;
824 ie1Prev = ie1;
825 }
826 //--- Use eT2 = eT: Get energy and parameters of angpar2 for this eT
827 G4double e2 = (maxEner2-minEner2) * eT + minEner2;
828 //----- Get parameter value corresponding to this e2
829 for(ie2=nDiscreteEnergies2; ie2<nEnergies2; ie2++) {
830 // G4cout << " GET IE2 " << ie2 << " - " << angpar2.theAngular[ie2].GetLabel() - e2 << " " << angpar2.theAngular[ie2].GetLabel() << " " << e2 <<G4endl;
831 if( (angpar2.theAngular[ie2].GetLabel() - e2) > 1.E-10*e2 ) break;
832 }
833 ie2Prev = ie2 - 1;
834 if( ie2 == 0 ) ie2Prev++;
835 if( ie2 == nEnergies2 ) {
836 ie2--;
837 ie2Prev = ie2;
838 }
839
840 //---- Energy corresponding to energy transformed
841 G4double eN = (theMaxEner-theMinEner) * eT + theMinEner;
842 if( std::getenv("G4PHPTEST2") ) G4cout << ie << " " << ie1 << " " << ie2 << " G4ParticleHPContAngularPar::loop eT " << eT << " -> eN " << eN << " e1 " << e1 << " e2 " << e2 << G4endl; //GDEB
843
844 theAngular[ie].SetLabel(eN);
845
846 for(G4int ip=0; ip<nAngularParameters; ip++) {
847 G4double val1 = theInt.Interpolate2(theManager.GetScheme(ie),
848 e1,
849 angpar1.theAngular[ie1Prev].GetLabel(),
850 angpar1.theAngular[ie1].GetLabel(),
851 angpar1.theAngular[ie1Prev].GetValue(ip),
852 angpar1.theAngular[ie1].GetValue(ip)) * (maxEner1-minEner1);
853 G4double val2 = theInt.Interpolate2(theManager.GetScheme(ie),
854 e2,
855 angpar2.theAngular[ie2Prev].GetLabel(),
856 angpar2.theAngular[ie2].GetLabel(),
857 angpar2.theAngular[ie2Prev].GetValue(ip),
858 angpar2.theAngular[ie2].GetValue(ip)) * (maxEner2-minEner2);
859
860 G4double value = theInt.Interpolate(aScheme, anEnergy,
861 angpar1.theEnergy, angpar2.theEnergy,
862 val1,
863 val2);
864 //value /= (theMaxEner-theMinEner);
865 if ( theMaxEner != theMinEner ) {
866 value /= (theMaxEner-theMinEner);
867 } else if ( value != 0 ) {
868 throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPContAngularPar::PrepareTableInterpolation theMaxEner == theMinEner and value != 0.");
869 }
870 if( std::getenv("G4PHPTEST2") ) G4cout << ie << " " << ip << " G4ParticleHPContAngularPar::Merge val1 " << val1 << " val2 " << val2 << " value " << value << G4endl; //GDEB
871 //- val1 = angpar1.theAngular[ie1-1].GetValue(ip) * (maxEner1-minEner1);
872 //- val2 = angpar2.theAngular[ie2-1].GetValue(ip) * (maxEner2-minEner2);
873 //- if( getenv("G4PHPTEST2") ) G4cout << ie << " " << ip << " G4ParticleHPContAngularPar::MergeOLD val1 " << val1 << " val2 " << val2 << " value " << value << G4endl; //GDEB
874
875 theAngular[ie].SetValue(ip, value);
876 }
877 }
878
879 if( std::getenv("G4PHPTEST2") ) {
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 }
887 }
888
889 void G4ParticleHPContAngularPar::Dump()
890 {
891 G4cout << theEnergy << " " << nEnergies << " " << nDiscreteEnergies << " " << nAngularParameters << G4endl;
892
893 for(G4int ii=0; ii<nEnergies; ii++) {
894 theAngular[ii].Dump();
895 }
896
897 }
898