Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/particle_hp/src/G4ParticleHPContAngularPar.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/hadronic/models/particle_hp/src/G4ParticleHPContAngularPar.cc (Version 11.3.0) and /processes/hadronic/models/particle_hp/src/G4ParticleHPContAngularPar.cc (Version 10.4.p3)


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