ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ParticleHPContAngularPar.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ParticleHPContAngularPar.cc
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 
44 #include "G4PhysicalConstants.hh"
45 #include "G4SystemOfUnits.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"
59 #include "G4IonTable.hh"
60 #include <set>
61 
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 
73  theProjectile = projectile;
74 
75  theEnergy = 0.0;
76  nEnergies = 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 
89  /*if( std::getenv("G4PHPTEST") )*/
90  theEnergy *= eV;
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 
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  {
120  if(Z==1) result->SetDefinition(G4Positron::Positron());
121  }
122  else if(A==1)
123  {
125  if(Z==1) result->SetDefinition(G4Proton::Proton());
126  }
127  else if(A==2)
128  {
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();
298 
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));
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();
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();
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));
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));
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  }
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;
606  aMan.Init(angularRep-10, nAngularParameters-1);
607  theStore.SetInterpolationManager(aMan); // Store interpolates f(costh)
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();
631  random,x1,x2,y1,y2);
632  G4ParticleHPVector theBuff1;
633  G4ParticleHPVector theBuff2;
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);
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);
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 
697 {
698 
699  //----- Discrete energies: store own energies in a map for faster searching
700  G4int ie;
701  for(ie=0; ie<nDiscreteEnergies; 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 
737  G4ParticleHPContAngularPar & angpar1,
738  G4ParticleHPContAngularPar & angpar2)
739 {
740  G4int ie,ie1,ie2, ie1Prev, ie2Prev;
742  theManager = angpar1.theManager;
743  theEnergy = anEnergy;
744 
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;
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
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++) {
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);
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 
890 {
891  G4cout << theEnergy << " " << nEnergies << " " << nDiscreteEnergies << " " << nAngularParameters << G4endl;
892 
893  for(G4int ii=0; ii<nEnergies; ii++) {
894  theAngular[ii].Dump();
895  }
896 
897 }