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.2)


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