Attachment #774 for problem #2468




    adjustResult = true;
  }

  G4ParticleHPContAngularPar(G4ParticleHPContAngularPar & val)
  {
    theEnergy         = val.theEnergy;
    nEnergies         = val.nEnergies;
    nDiscreteEnergies = val.nDiscreteEnergies;
    nAngularParameters= val.nAngularParameters;
    theProjectile     = val.theProjectile;
    theManager        = val.theManager;
    theInt            = val.theInt;
    adjustResult      = val.adjustResult;
    theMinEner        = val.theMinEner;
    theMaxEner        = val.theMaxEner;
    theEnergiesTransformed = val.theEnergiesTransformed;
    theDiscreteEnergies = val.theDiscreteEnergies;
    theDiscreteEnergiesOwn = val.theDiscreteEnergiesOwn;
    fCache.Put(0);
    theAngular        = new G4ParticleHPList[nEnergies];
    for(int ie=0;ie<nEnergies;++ie) {
      theAngular[ie].SetLabel(val.theAngular[ie].GetLabel());
      for(int ip=0;ip<nAngularParameters;++ip) {
	theAngular[ie].SetValue(ip,val.theAngular[ie].GetValue(ip));
      }
    }
  }

  G4ParticleHPContAngularPar(G4ParticleDefinition* projectile); 

  ~G4ParticleHPContAngularPar()

  G4ReactionProduct* Sample(G4double anEnergy, G4double massCode, G4double mass, 
                            G4int angularRep, G4int interpol);
  
  G4double GetEnergy() const
  { 
    if( std::getenv("G4PHPTEST") )
      G4cout << this << " G4ParticleHPContAngularPar::GetEnergy "

             G4ParticleHPContAngularPar & store2);
    // NOTE: this interpolates legendre coefficients

  void PrepareTableInterpolation(); // SM 26 Apr 2022 const G4ParticleHPContAngularPar* angularPrev);
  
  G4double MeanEnergyOfThisInteraction()
  {

    fCache.Get()->fresh = true;
  }

  void Dump() const;

private:
  

Lines 78-86 Link Here

(-)source-orig/processes/hadronic/models/particle_hp/include/G4ParticleHPContEnergyAngular.hh (-2 / +2 lines)
78	theAngular[i].SetInterpolation(theInterpolation);	78	theAngular[i].SetInterpolation(theInterpolation);
79	#ifndef PHP_AS_HP	79	#ifndef PHP_AS_HP
80	if( i != 0 ) {	80	if( i != 0 ) {
81	theAngular[i].PrepareTableInterpolation(&(theAngular[i-1]));	81	theAngular[i].PrepareTableInterpolation(); // SM 26 Apr 2022 &(theAngular[i-1]));
82	} else {	82	} else {
83	theAngular[i].PrepareTableInterpolation((G4ParticleHPContAngularPar*)0);	83	theAngular[i].PrepareTableInterpolation(); // SM 26 Apr 2022 ((G4ParticleHPContAngularPar*)0);
84	}	84	}
85	#endif	85	#endif
86	}	86	}

Lines 59-64 Link Here

(-)source-orig/processes/hadronic/models/particle_hp/src/G4ParticleHPContAngularPar.cc (-82 / +310 lines)
59	#include "G4IonTable.hh"	59	#include "G4IonTable.hh"
60	#include "G4ParticleHPManager.hh"	60	#include "G4ParticleHPManager.hh"
61	#include <set>	61	#include <set>
		62	#include <vector>
62		63
63	G4ParticleHPContAngularPar::G4ParticleHPContAngularPar( G4ParticleDefinition* projectile )	64	G4ParticleHPContAngularPar::G4ParticleHPContAngularPar( G4ParticleDefinition* projectile )
64	{	65	{
Lines 168-173 Link Here
168	//Continues Emssions, low to high LAST	169	//Continues Emssions, low to high LAST
169	fCache.Get()->remaining_energy = std::max ( theAngular[0].GetLabel() , theAngular[nEnergies-1].GetLabel() );	170	fCache.Get()->remaining_energy = std::max ( theAngular[0].GetLabel() , theAngular[nEnergies-1].GetLabel() );
170	fCache.Get()->fresh = false;	171	fCache.Get()->fresh = false;
		172	if( std::getenv("G4PHPTEST") ) G4cout << " G4ParticleHPContAngularPar::Sample fresh remaining_energy " << fCache.Get()->remaining_energy << G4endl;
		173
171	}	174	}
172		175
173	//Cheating for small remaining_energy	176	//Cheating for small remaining_energy
Lines 175-180 Link Here
175	if ( nDiscreteEnergies == nEnergies )	178	if ( nDiscreteEnergies == nEnergies )
176	{	179	{
177	fCache.Get()->remaining_energy = std::max ( fCache.Get()->remaining_energy , theAngular[nDiscreteEnergies-1].GetLabel() ); //Minimum Line	180	fCache.Get()->remaining_energy = std::max ( fCache.Get()->remaining_energy , theAngular[nDiscreteEnergies-1].GetLabel() ); //Minimum Line
		181	if( std::getenv("G4PHPTEST") ) G4cout << " G4ParticleHPContAngularPar::Sample cheat case for small remaining_energy " << fCache.Get()->remaining_energy << G4endl;
178	}	182	}
179	else	183	else
180	{	184	{
Lines 187-192 Link Here
187	if ( theAngular[j].GetValue(0) != 0.0 ) break;	191	if ( theAngular[j].GetValue(0) != 0.0 ) break;
188	}	192	}
189	fCache.Get()->remaining_energy = std::max ( fCache.Get()->remaining_energy , std::min ( theAngular[nDiscreteEnergies-1].GetLabel() , cont_min ) ); //Minimum Line or grid	193	fCache.Get()->remaining_energy = std::max ( fCache.Get()->remaining_energy , std::min ( theAngular[nDiscreteEnergies-1].GetLabel() , cont_min ) ); //Minimum Line or grid
		194	if( std::getenv("G4PHPTEST") ) G4cout << " G4ParticleHPContAngularPar::Sample max line or grid remaining_energy " << fCache.Get()->remaining_energy << G4endl;
190	}	195	}
191	//	196	//
192	G4double random = G4UniformRand();	197	G4double random = G4UniformRand();
Lines 197-203 Link Here
197	for ( G4int j = 0 ; j < nDiscreteEnergies ; j++ )	202	for ( G4int j = 0 ; j < nDiscreteEnergies ; j++ )
198	{	203	{
199	G4double delta = 0.0;	204	G4double delta = 0.0;
200	if ( theAngular[j].GetLabel() <= fCache.Get()->remaining_energy ) delta = theAngular[i].GetValue(0);	205	if ( theAngular[j].GetLabel() <= fCache.Get()->remaining_energy ) delta = theAngular[j].GetValue(0);
201	running[j+1] = running[j] + delta;	206	running[j+1] = running[j] + delta;
202	}	207	}
203	G4double tot_prob_DIS = running[ nDiscreteEnergies ];	208	G4double tot_prob_DIS = running[ nDiscreteEnergies ];
Lines 327-333 Link Here
327	}	332	}
328		333
329	//TK080711	334	//TK080711
330	if( adjustResult ) fCache.Get()->remaining_energy -= fsEnergy;	335	// SM 30 March 2022
		336	// the remaining energy needs to be lowered by the photon
		337	// energy in any case. Otherwise additional photons with
		338	// too high energy will be produced - therefore the
		339	// adjustResult condition has been removed
		340	//
		341	//if( adjustResult ) fCache.Get()->remaining_energy -= fsEnergy;
		342	if( std::getenv("G4PHPTEST") ) G4cout << " G4ParticleHPContAngularPar::Sample remaining_energy - fsEnergy " << fCache.Get()->remaining_energy << " - " << fsEnergy << G4endl;
		343	fCache.Get()->remaining_energy -= fsEnergy;
331	//TK080711	344	//TK080711
332		345
333	//080801b	346	//080801b
Lines 343-348 Link Here
343	{	356	{
344	fCache.Get()->remaining_energy = theAngular[ nEnergies-1 ].GetLabel();	357	fCache.Get()->remaining_energy = theAngular[ nEnergies-1 ].GetLabel();
345	fCache.Get()->fresh = false;	358	fCache.Get()->fresh = false;
		359	if( std::getenv("G4PHPTEST") ) G4cout << " G4ParticleHPContAngularPar::Sample 2nd fresh remaining_energy " << fCache.Get()->remaining_energy << G4endl;
346	}	360	}
347	//080714	361	//080714
348	G4double random = G4UniformRand();	362	G4double random = G4UniformRand();
Lines 462-468 Link Here
462	delete [] running;	476	delete [] running;
463		477
464	//080714	478	//080714
465	if( adjustResult ) fCache.Get()->remaining_energy -= fsEnergy;	479	// SM 30 March 2022
		480	// the remaining energy needs to be lowered by the photon
		481	// energy in any case. Otherwise additional photons with
		482	// too high energy will be produced - therefore the
		483	// adjustResult condition has been removed
		484	//
		485	//if( adjustResult ) fCache.Get()->remaining_energy -= fsEnergy;
		486	if( std::getenv("G4PHPTEST") ) G4cout << " G4ParticleHPContAngularPar::Sample 2nd remaining_energy - fsEnergy " << fCache.Get()->remaining_energy << " - " << fsEnergy << G4endl;
		487	fCache.Get()->remaining_energy -= fsEnergy;
466	//080714	488	//080714
467	}	489	}
468	}	490	}
Lines 697-710 Link Here
697		719
698	#define MERGE_NEW	720	#define MERGE_NEW
699		721
700	void G4ParticleHPContAngularPar::PrepareTableInterpolation(const G4ParticleHPContAngularPar* angParPrev)	722	void G4ParticleHPContAngularPar::PrepareTableInterpolation() // SM 26 Apr 2022 const G4ParticleHPContAngularPar* angParPrev)
701	{	723	{
702		724
703	//----- Discrete energies: store own energies in a map for faster searching	725	//----- Discrete energies: store own energies in a map for faster searching
		726	/*
		727	* SM 27 Apr 2022
		728	*
		729	* the data files sometimes have identical discrete energies (likely typos)
		730	* which would lead to overwriting the already existing index and hence
		731	* creating a hole in the lookup table.
		732	* No attempt is made here to correct for the energies - rather an epsilon
		733	* is subtracted from the energy in order to uniquely identify the line
		734	*/
704	G4int ie;	735	G4int ie;
705	for(ie=0; ie<nDiscreteEnergies; ie++) {	736	for(ie=0; ie<nDiscreteEnergies; ie++) {
706	theDiscreteEnergiesOwn[theAngular[ie].GetLabel()] = ie;	737	// check if energy is already present and subtract epsilon if that's the case
		738	double myE = theAngular[ie].GetLabel();
		739	while ( theDiscreteEnergiesOwn.find(myE) != theDiscreteEnergiesOwn.end() ) {
		740	myE -= 1e-6;
		741	}
		742	theDiscreteEnergiesOwn[myE] = ie;
707	}	743	}
		744	/*
		745	* SM 26 Apr 2022
		746	*
		747	* the approach here makes no sense. It would work only for two sets that
		748	* have identical min and max energy. If the 2 sets differ in min, max or
		749	* both, the energy inserted would be normalized to its original set but
		750	* interpreted with the new - which is not correct.
		751	*
		752	* Disable the code for now and simply return ...
		753	*/
		754
		755	return;
		756
		757	/* SM 26 Apr 2022
		758	*
		759
708	if( !angParPrev ) return;	760	if( !angParPrev ) return;
709		761
710	//----- Discrete energies: use energies that appear in one or another	762	//----- Discrete energies: use energies that appear in one or another
Lines 721-727 Link Here
721	G4double ener = theAngular[ie].GetLabel();	773	G4double ener = theAngular[ie].GetLabel();
722	G4double enerT = (ener-theMinEner)/(theMaxEner-theMinEner);	774	G4double enerT = (ener-theMinEner)/(theMaxEner-theMinEner);
723	theEnergiesTransformed.insert(enerT);	775	theEnergiesTransformed.insert(enerT);
724	//- if( getenv("G4PHPTEST2") ) G4cout <<this << " G4ParticleHPContAngularPar::PrepareTableInterpolation theEnergiesTransformed1 " << enerT << G4endl; //GDEB	776	if( getenv("G4PHPTEST2") ) G4cout <<this << " G4ParticleHPContAngularPar::PrepareTableInterpolation theEnergiesTransformed1 " << enerT << G4endl; //GDEB
725	}	777	}
726	G4int nEnergiesPrev = angParPrev->GetNEnergies();	778	G4int nEnergiesPrev = angParPrev->GetNEnergies();
727	G4double minEnerPrev = angParPrev->GetMinEner();	779	G4double minEnerPrev = angParPrev->GetMinEner();
Lines 730-740 Link Here
730	G4double ener = angParPrev->theAngular[ie].GetLabel();	782	G4double ener = angParPrev->theAngular[ie].GetLabel();
731	G4double enerT = (ener-minEnerPrev)/(maxEnerPrev-minEnerPrev);	783	G4double enerT = (ener-minEnerPrev)/(maxEnerPrev-minEnerPrev);
732	theEnergiesTransformed.insert(enerT);	784	theEnergiesTransformed.insert(enerT);
733	//- if( getenv("G4PHPTEST2") ) G4cout << this << " G4ParticleHPContAngularPar::PrepareTableInterpolation theEnergiesTransformed2 " << enerT << G4endl; //GDEB	785	if( getenv("G4PHPTEST2") ) G4cout << this << " G4ParticleHPContAngularPar::PrepareTableInterpolation theEnergiesTransformed2 " << enerT << G4endl; //GDEB
734	}	786	}
735	// add the maximum energy	787	// add the maximum energy
736	theEnergiesTransformed.insert(1.);	788	//theEnergiesTransformed.insert(1.);
737		789
		790	*
		791	*/ // SM 26 Apr 2022
738	}	792	}
739		793
740	void G4ParticleHPContAngularPar::BuildByInterpolation(G4double anEnergy, G4InterpolationScheme aScheme,	794	void G4ParticleHPContAngularPar::BuildByInterpolation(G4double anEnergy, G4InterpolationScheme aScheme,
Lines 742-825 Link Here
742	G4ParticleHPContAngularPar & angpar2)	796	G4ParticleHPContAngularPar & angpar2)
743	{	797	{
744	G4int ie,ie1,ie2, ie1Prev, ie2Prev;	798	G4int ie,ie1,ie2, ie1Prev, ie2Prev;
745	nAngularParameters = angpar1.nAngularParameters;	799	// SM 12.Feb 2022
746	theManager = angpar1.theManager;	800	//
		801	// only re-build the interpolation table if we have a new interaction. For several subsequent
		802	// samplings of final state particles in the same interaction we should use the existing table.
		803	if ( fCache.Get()->fresh != true ) return;
		804
		805	// SM 10. Apr 2022
		806	//
		807	// make copies of angpar1 and angpar2. Since these are given by reference it can not
		808	// be excluded that one of them is this. Hence this code uses potentially the old this
		809	// for creating the new this - which leads to memory corruption if the old is not stored
		810	// as separarte object for lookup
		811	const G4ParticleHPContAngularPar copyAngpar1(angpar1),copyAngpar2(angpar2);
		812
		813	nAngularParameters = copyAngpar1.nAngularParameters;
		814	theManager = copyAngpar1.theManager;
747	theEnergy = anEnergy;	815	theEnergy = anEnergy;
748		816	theMinEner = DBL_MAX; // SM 27 Apr 2022 min and max will be re-calculated after interpolation
749	nDiscreteEnergies = theDiscreteEnergies.size();	817	theMaxEner = -DBL_MAX;
750	std::set<G4double>::const_iterator itede;	818	// SM 11. Feb 2022
751	std::map<G4double,G4int> discEnerOwn1 = angpar1.GetDiscreteEnergiesOwn();	819	//
752	std::map<G4double,G4int> discEnerOwn2 = angpar2.GetDiscreteEnergiesOwn();	820	// re-wrote most of this. The two discrete sets have to be merged. A vector holds the temporary
753	std::map<G4double,G4int>::const_iterator itedeo;	821	// data to be copied to the array in the end. Since the G4ParticleHPList class contains pointers
754	ie = 0;	822	// one can not simply assign elements of this type. Each member needs to call the explicit Set()
755	for( itede = theDiscreteEnergies.begin(); itede != theDiscreteEnergies.end(); itede++, ie++ ) {	823	// method instead.
756	G4double discEner = *itede;	824	//
757	itedeo = discEnerOwn1.find(discEner);	825	// First average probabilities for those lines that are in both sets
758	if( itedeo == discEnerOwn1.end() ) {	826	const std::map<G4double,G4int> discEnerOwn1 = copyAngpar1.GetDiscreteEnergiesOwn();
759	ie1 = 0;	827	const std::map<G4double,G4int> discEnerOwn2 = copyAngpar2.GetDiscreteEnergiesOwn();
760	} else {	828	std::map<G4double,G4int>::const_iterator itedeo1;
761	ie1 = -1;	829	std::map<G4double,G4int>::const_iterator itedeo2;
		830	std::vector<G4ParticleHPList *> vAngular(discEnerOwn1.size());
		831	for( itedeo1 = discEnerOwn1.begin(); itedeo1 != discEnerOwn1.end(); ++itedeo1) {
		832	G4double discEner = itedeo1->first;
		833	if ( discEner < theMinEner ) {
		834	theMinEner = discEner;
762	}	835	}
763	itedeo = discEnerOwn2.find(discEner);	836	if ( discEner > theMaxEner ) {
764	if( itedeo == discEnerOwn2.end() ) {	837	theMaxEner = discEner;
765	ie2 = 0;	838	}
766	} else {	839	ie1 = itedeo1->second;
		840	itedeo2 = discEnerOwn2.find(discEner);
		841	if( itedeo2 == discEnerOwn2.end() ) {
767	ie2 = -1;	842	ie2 = -1;
		843	} else {
		844	ie2 = itedeo2->second;
768	}	845	}
769		846	vAngular[ie1] = new G4ParticleHPList();
770	theAngular[ie].SetLabel(discEner);	847	vAngular[ie1]->SetLabel(copyAngpar1.theAngular[ie1].GetLabel());
771	G4double val1, val2;	848	G4double val1, val2;
772	for(G4int ip=0; ip<nAngularParameters; ip++) {	849	for(G4int ip=0; ip<nAngularParameters; ++ip) {
773	if( ie1 != -1 ) {	850	val1 = copyAngpar1.theAngular[ie1].GetValue(ip);
774	val1 = angpar1.theAngular[ie1].GetValue(ip);
775	} else {
776	val1 = 0.;
777	}
778	if( ie2 != -1 ) {	851	if( ie2 != -1 ) {
779	val2 = angpar2.theAngular[ie2].GetValue(ip);	852	val2 = copyAngpar2.theAngular[ie2].GetValue(ip);
780	} else {	853	} else {
781	val2 = 0.;	854	val2 = 0.;
782	}	855	}
		856	G4double value = theInt.Interpolate(aScheme, anEnergy,
		857	copyAngpar1.theEnergy, copyAngpar2.theEnergy,
		858	val1,
		859	val2);
		860	if( std::getenv("G4PHPTEST2") ) G4cout << ie1 << " " << ip << " G4ParticleHPContAngularPar::Merge DiscreteEnergies val1 " << val1 << " val2 " << val2 << " value " << value << G4endl; //GDEB
783		861
		862	vAngular[ie1]->SetValue(ip, value);
		863	}
		864	}
		865	// add the ones in set2 but not in set1
		866	std::vector<G4ParticleHPList *>::iterator itv;
		867	for( itedeo2 = discEnerOwn2.begin(); itedeo2 != discEnerOwn2.end(); ++itedeo2) {
		868	G4double discEner = itedeo2->first;
		869	ie2=itedeo2->second;
		870	bool notFound = true;
		871	itedeo1 = discEnerOwn1.find(discEner);
		872	if( itedeo1 != discEnerOwn1.end() ) {
		873	notFound = false;
		874	}
		875	if( notFound ) {
		876	// not yet in list
		877	if ( discEner < theMinEner ) {
		878	theMinEner = discEner;
		879	}
		880	if ( discEner > theMaxEner ) {
		881	theMaxEner = discEner;
		882	}
		883	// find position to insert
		884	bool isInserted=false;
		885	ie=0;
		886	for(itv=vAngular.begin();itv!=vAngular.end();++itv,++ie) {
		887	if ( discEner > (*itv)->GetLabel() ) {
		888	itv=vAngular.insert(itv,new G4ParticleHPList);
		889	(*itv)->SetLabel(copyAngpar2.theAngular[ie2].GetLabel());
		890	isInserted=true;
		891	break;
		892	}
		893	}
		894	if (!isInserted) {
		895	ie=vAngular.size();
		896	vAngular.push_back(new G4ParticleHPList);
		897	vAngular[ie]->SetLabel(copyAngpar2.theAngular[ie2].GetLabel());
		898	isInserted=true;
		899	}
		900	G4double val1, val2;
		901	for(G4int ip=0; ip<nAngularParameters; ++ip) {
		902	val1 = 0;
		903	val2 = copyAngpar2.theAngular[ie2].GetValue(ip);
784	G4double value = theInt.Interpolate(aScheme, anEnergy,	904	G4double value = theInt.Interpolate(aScheme, anEnergy,
785	angpar1.theEnergy, angpar2.theEnergy,	905	copyAngpar1.theEnergy, copyAngpar2.theEnergy,
786	val1,	906	val1,
787	val2);	907	val2);
788	if( std::getenv("G4PHPTEST2") ) G4cout << ie << " " << ip << " G4ParticleHPContAngularPar::Merge DiscreteEnergies val1 " << val1 << " val2 " << val2 << " value " << value << G4endl; //GDEB	908	if( std::getenv("G4PHPTEST2") ) G4cout << ie << " " << ip << " G4ParticleHPContAngularPar::Merge DiscreteEnergies val1 " << val1 << " val2 " << val2 << " value " << value << G4endl; //GDEB
789		909
790	theAngular[ie].SetValue(ip, value);	910	vAngular[ie]->SetValue(ip, value);
		911	}
791	}	912	}
792	}	913	}
793		914
		915	// store new discrete list
		916	nDiscreteEnergies = vAngular.size();
794	if(theAngular != 0) delete [] theAngular;	917	if(theAngular != 0) delete [] theAngular;
795	nEnergies = nDiscreteEnergies + angpar2.GetNEnergiesTransformed();	918	theAngular = 0;
796	theAngular = new G4ParticleHPList [nEnergies];	919	if ( nDiscreteEnergies > 0 ) {
		920	theAngular = new G4ParticleHPList [nDiscreteEnergies];
		921	}
		922	theDiscreteEnergiesOwn.clear();
		923	theDiscreteEnergies.clear();
		924	for(ie=0; ie<nDiscreteEnergies; ++ie) {
		925	theAngular[ie].SetLabel(vAngular[ie]->GetLabel());
		926	for(int ip=0; ip<nAngularParameters; ++ip) {
		927	theAngular[ie].SetValue(ip,vAngular[ie]->GetValue(ip));
		928	}
		929	theDiscreteEnergiesOwn[theAngular[ie].GetLabel()] = ie;
		930	theDiscreteEnergies.insert(theAngular[ie].GetLabel());
		931	}
		932
		933	/*
		934	* SM 26 Apr 2022
		935	*
		936	* the continuous energies need to be made from scratch like the discrete
		937	* ones. Therefore the re-assignemnt of theAngular needs to be done
		938	* after the continuous energy set is also finalized. Only then the
		939	* total number of nEnergies is known and the array can be allocated.
		940	*
		941	* Hence the code here is commented out.
		942
		943	nEnergies = nDiscreteEnergies + copyAngpar2.GetNEnergiesTransformed();
		944	if ( nEnergies > nDiscreteEnergies ) {
		945	G4ParticleHPList * theNewAngular = new G4ParticleHPList [nEnergies];
		946	if (theAngular !=0 ) {
		947	for(ie=0; ie<nDiscreteEnergies; ++ie) {
		948	theNewAngular[ie].SetLabel(theAngular[ie].GetLabel());
		949	for(int ip=0; ip<nAngularParameters; ++ip) {
		950	theNewAngular[ie].SetValue(ip,theAngular[ie].GetValue(ip));
		951	}
		952	}
		953	delete [] theAngular;
		954	}
		955	theAngular = theNewAngular;
		956	}
		957
		958	*
		959	*/ // SM 26 Apr 2022
797		960
798	//---- Get minimum and maximum energy interpolating	961	//---- Get minimum and maximum energy interpolating
799	theMinEner = angpar1.GetMinEner() + (theEnergy-angpar1.GetEnergy()) * (angpar2.GetMinEner()-angpar1.GetMinEner())/(angpar2.GetEnergy()-angpar1.GetEnergy());	962	// SM 27 Apr 2022 don't use theMinEner or theMaxEner here, since the transformed
800	theMaxEner = angpar1.GetMaxEner() + (theEnergy-angpar1.GetEnergy()) * (angpar2.GetMaxEner()-angpar1.GetMaxEner())/(angpar2.GetEnergy()-angpar1.GetEnergy());	963	// energies need the interpolated range from the original Angpar
		964	G4double interMinEner = copyAngpar1.GetMinEner() + (theEnergy-copyAngpar1.GetEnergy()) * (copyAngpar2.GetMinEner()-copyAngpar1.GetMinEner())/(copyAngpar2.GetEnergy()-copyAngpar1.GetEnergy());
		965	G4double interMaxEner = copyAngpar1.GetMaxEner() + (theEnergy-copyAngpar1.GetEnergy()) * (copyAngpar2.GetMaxEner()-copyAngpar1.GetMaxEner())/(copyAngpar2.GetEnergy()-copyAngpar1.GetEnergy());
801		966
802	if( std::getenv("G4PHPTEST2") ) G4cout << " G4ParticleHPContAngularPar::Merge E " << anEnergy << " minmax " << theMinEner << " " << theMaxEner << G4endl; //GDEB	967	if( std::getenv("G4PHPTEST2") ) G4cout << " G4ParticleHPContAngularPar::Merge E " << anEnergy << " minmax " << interMinEner << " " << interMaxEner << G4endl; //GDEB
803		968
804	//--- Loop to energies of new set	969	//--- Loop to energies of new set
805	std::set<G4double> energiesTransformed = angpar2.GetEnergiesTransformed();	970	//SM 26 Apr 2022 std::set<G4double> energiesTransformed = copyAngpar2.GetEnergiesTransformed();
806	std::set<G4double>::const_iterator iteet = energiesTransformed.begin();	971	//SM 26 Apr 2022 std::set<G4double>::const_iterator iteet = energiesTransformed.begin();
807	G4int nEnergies1 = angpar1.GetNEnergies();	972	theEnergiesTransformed.clear();
808	G4int nDiscreteEnergies1 = angpar1.GetNDiscreteEnergies();	973
809	G4double minEner1 = angpar1.GetMinEner();	974	G4int nEnergies1 = copyAngpar1.GetNEnergies();
810	G4double maxEner1 = angpar1.GetMaxEner();	975	G4int nDiscreteEnergies1 = copyAngpar1.GetNDiscreteEnergies();
811	G4int nEnergies2 = angpar2.GetNEnergies();	976	G4double minEner1 = copyAngpar1.GetMinEner();
812	G4int nDiscreteEnergies2 = angpar2.GetNDiscreteEnergies();	977	G4double maxEner1 = copyAngpar1.GetMaxEner();
813	G4double minEner2 = angpar2.GetMinEner();	978	G4int nEnergies2 = copyAngpar2.GetNEnergies();
814	G4double maxEner2 = angpar2.GetMaxEner();	979	G4int nDiscreteEnergies2 = copyAngpar2.GetNDiscreteEnergies();
815	for(ie=nDiscreteEnergies; ie<nEnergies; ie++,iteet++) {	980	G4double minEner2 = copyAngpar2.GetMinEner();
		981	G4double maxEner2 = copyAngpar2.GetMaxEner();
		982
		983	// SM 26 Apr 2022: first build the list of transformed energies normalized
		984	// to the new min max by assuming that the min-max range of
		985	// each set would be scalable to the new, interpolated min
		986	// max range
		987
		988	for(ie1=nDiscreteEnergies1; ie1<nEnergies1; ++ie1) {
		989	G4double e1 = copyAngpar1.theAngular[ie1].GetLabel();
		990	G4double eT = (e1-minEner1);
		991	if ( maxEner1 != minEner1) {
		992	eT /= (maxEner1-minEner1);
		993	}
		994	if ( eT >= 0 && eT <= 1 ) {
		995	theEnergiesTransformed.insert(eT);
		996	}
		997	}
		998	for(ie2=nDiscreteEnergies2; ie2<nEnergies2; ++ie2) {
		999	G4double e2 = copyAngpar2.theAngular[ie2].GetLabel();
		1000	G4double eT = (e2-minEner2);
		1001	if ( maxEner2 != minEner2) {
		1002	eT /= (maxEner2-minEner2);
		1003	}
		1004	if ( eT >= 0 && eT <= 1 ) {
		1005	theEnergiesTransformed.insert(eT);
		1006	}
		1007	}
		1008
		1009	// now the list of energies is complete
		1010
		1011	nEnergies = nDiscreteEnergies+theEnergiesTransformed.size();
		1012
		1013	// create final array of angular parameters
		1014
		1015	G4ParticleHPList * theNewAngular = new G4ParticleHPList [nEnergies];
		1016
		1017	// copy the discrete energies and interpolated parameters to new array
		1018
		1019	if (theAngular !=0 ) {
		1020	for(ie=0; ie<nDiscreteEnergies; ++ie) {
		1021	theNewAngular[ie].SetLabel(theAngular[ie].GetLabel());
		1022	for(int ip=0; ip<nAngularParameters; ++ip) {
		1023	theNewAngular[ie].SetValue(ip,theAngular[ie].GetValue(ip));
		1024	}
		1025	}
		1026	delete [] theAngular;
		1027	}
		1028	theAngular = theNewAngular;
		1029
		1030	// interpolate the continuous energies for new array
		1031	std::set<G4double>::const_iterator iteet = theEnergiesTransformed.begin();
		1032
		1033	for(ie=nDiscreteEnergies; ie<nEnergies; ++ie,++iteet) {
816	G4double eT = (*iteet);	1034	G4double eT = (*iteet);
817		1035
818	//--- Use eT1 = eT: Get energy and parameters of angpar1 for this eT	1036	//--- Use eT1 = eT: Get energy and parameters of copyAngpar1 for this eT
819	G4double e1 = (maxEner1-minEner1) * eT + minEner1;	1037	G4double e1 = (maxEner1-minEner1) * eT + minEner1;
820	//----- Get parameter value corresponding to this e1	1038	//----- Get parameter value corresponding to this e1
821	for(ie1=nDiscreteEnergies1; ie1<nEnergies1; ie1++) {	1039	for(ie1=nDiscreteEnergies1; ie1<nEnergies1; ++ie1) {
822	if( (angpar1.theAngular[ie1].GetLabel() - e1) > 1.E-10*e1 ) break;	1040	if( (copyAngpar1.theAngular[ie1].GetLabel() - e1) > 1.E-10*e1 ) break;
823	}	1041	}
824	ie1Prev = ie1 - 1;	1042	ie1Prev = ie1 - 1;
825	if( ie1 == 0 ) ie1Prev++;	1043	if( ie1 == 0 ) ie1Prev++;
Lines 827-838 Link Here
827	ie1--;	1045	ie1--;
828	ie1Prev = ie1;	1046	ie1Prev = ie1;
829	}	1047	}
830	//--- Use eT2 = eT: Get energy and parameters of angpar2 for this eT	1048	//--- Use eT2 = eT: Get energy and parameters of copyAngpar2 for this eT
831	G4double e2 = (maxEner2-minEner2) * eT + minEner2;	1049	G4double e2 = (maxEner2-minEner2) * eT + minEner2;
832	//----- Get parameter value corresponding to this e2	1050	//----- Get parameter value corresponding to this e2
833	for(ie2=nDiscreteEnergies2; ie2<nEnergies2; ie2++) {	1051	for(ie2=nDiscreteEnergies2; ie2<nEnergies2; ++ie2) {
834	// G4cout << " GET IE2 " << ie2 << " - " << angpar2.theAngular[ie2].GetLabel() - e2 << " " << angpar2.theAngular[ie2].GetLabel() << " " << e2 <<G4endl;	1052	// G4cout << " GET IE2 " << ie2 << " - " << copyAngpar2.theAngular[ie2].GetLabel() - e2 << " " << copyAngpar2.theAngular[ie2].GetLabel() << " " << e2 <<G4endl;
835	if( (angpar2.theAngular[ie2].GetLabel() - e2) > 1.E-10*e2 ) break;	1053	if( (copyAngpar2.theAngular[ie2].GetLabel() - e2) > 1.E-10*e2 ) break;
836	}	1054	}
837	ie2Prev = ie2 - 1;	1055	ie2Prev = ie2 - 1;
838	if( ie2 == 0 ) ie2Prev++;	1056	if( ie2 == 0 ) ie2Prev++;
Lines 842-896 Link Here
842	}	1060	}
843		1061
844	//---- Energy corresponding to energy transformed	1062	//---- Energy corresponding to energy transformed
845	G4double eN = (theMaxEner-theMinEner) * eT + theMinEner;	1063	G4double eN = (interMaxEner-interMinEner) * eT + interMinEner;
846	if( std::getenv("G4PHPTEST2") ) G4cout << ie << " " << ie1 << " " << ie2 << " G4ParticleHPContAngularPar::loop eT " << eT << " -> eN " << eN << " e1 " << e1 << " e2 " << e2 << G4endl; //GDEB	1064	if( std::getenv("G4PHPTEST2") ) G4cout << ie << " " << ie1 << " " << ie2 << " G4ParticleHPContAngularPar::loop eT " << eT << " -> eN " << eN << " e1 " << e1 << " e2 " << e2 << G4endl; //GDEB
847		1065
848	theAngular[ie].SetLabel(eN);	1066	theAngular[ie].SetLabel(eN);
		1067	if ( eN < theMinEner ) {
		1068	theMinEner = eN;
		1069	}
		1070	if ( eN > theMaxEner ) {
		1071	theMaxEner = eN;
		1072	}
849		1073
850	for(G4int ip=0; ip<nAngularParameters; ip++) {	1074	for(G4int ip=0; ip<nAngularParameters; ++ip) {
851	G4double val1 = theInt.Interpolate2(theManager.GetScheme(ie),	1075	G4double val1 = theInt.Interpolate2(theManager.GetScheme(ie),
852	e1,	1076	e1,
853	angpar1.theAngular[ie1Prev].GetLabel(),	1077	copyAngpar1.theAngular[ie1Prev].GetLabel(),
854	angpar1.theAngular[ie1].GetLabel(),	1078	copyAngpar1.theAngular[ie1].GetLabel(),
855	angpar1.theAngular[ie1Prev].GetValue(ip),	1079	copyAngpar1.theAngular[ie1Prev].GetValue(ip),
856	angpar1.theAngular[ie1].GetValue(ip)) * (maxEner1-minEner1);	1080	copyAngpar1.theAngular[ie1].GetValue(ip)) * (maxEner1-minEner1);
857	G4double val2 = theInt.Interpolate2(theManager.GetScheme(ie),	1081	G4double val2 = theInt.Interpolate2(theManager.GetScheme(ie),
858	e2,	1082	e2,
859	angpar2.theAngular[ie2Prev].GetLabel(),	1083	copyAngpar2.theAngular[ie2Prev].GetLabel(),
860	angpar2.theAngular[ie2].GetLabel(),	1084	copyAngpar2.theAngular[ie2].GetLabel(),
861	angpar2.theAngular[ie2Prev].GetValue(ip),	1085	copyAngpar2.theAngular[ie2Prev].GetValue(ip),
862	angpar2.theAngular[ie2].GetValue(ip)) * (maxEner2-minEner2);	1086	copyAngpar2.theAngular[ie2].GetValue(ip)) * (maxEner2-minEner2);
863		1087
864	G4double value = theInt.Interpolate(aScheme, anEnergy,	1088	G4double value = theInt.Interpolate(aScheme, anEnergy,
865	angpar1.theEnergy, angpar2.theEnergy,	1089	copyAngpar1.theEnergy, copyAngpar2.theEnergy,
866	val1,	1090	val1,
867	val2);	1091	val2);
868	//value /= (theMaxEner-theMinEner);	1092	//value /= (interMaxEner-interMinEner);
869	if ( theMaxEner != theMinEner ) {	1093	if ( interMaxEner != interMinEner ) {
870	value /= (theMaxEner-theMinEner);	1094	value /= (interMaxEner-interMinEner);
871	} else if ( value != 0 ) {	1095	} else if ( value != 0 ) {
872	throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPContAngularPar::PrepareTableInterpolation theMaxEner == theMinEner and value != 0.");	1096	throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPContAngularPar::PrepareTableInterpolation interMaxEner == interMinEner and value != 0.");
873	}	1097	}
874	if( std::getenv("G4PHPTEST2") ) G4cout << ie << " " << ip << " G4ParticleHPContAngularPar::Merge val1 " << val1 << " val2 " << val2 << " value " << value << G4endl; //GDEB	1098	if( std::getenv("G4PHPTEST2") ) G4cout << ie << " " << ip << " G4ParticleHPContAngularPar::Merge val1 " << val1 << " val2 " << val2 << " value " << value << G4endl; //GDEB
875	//- val1 = angpar1.theAngular[ie1-1].GetValue(ip) * (maxEner1-minEner1);	1099	//- val1 = copyAngpar1.theAngular[ie1-1].GetValue(ip) * (maxEner1-minEner1);
876	//- val2 = angpar2.theAngular[ie2-1].GetValue(ip) * (maxEner2-minEner2);	1100	//- val2 = copyAngpar2.theAngular[ie2-1].GetValue(ip) * (maxEner2-minEner2);
877	//- if( getenv("G4PHPTEST2") ) G4cout << ie << " " << ip << " G4ParticleHPContAngularPar::MergeOLD val1 " << val1 << " val2 " << val2 << " value " << value << G4endl; //GDEB	1101	//- if( getenv("G4PHPTEST2") ) G4cout << ie << " " << ip << " G4ParticleHPContAngularPar::MergeOLD val1 " << val1 << " val2 " << val2 << " value " << value << G4endl; //GDEB
878		1102
879	theAngular[ie].SetValue(ip, value);	1103	theAngular[ie].SetValue(ip, value);
880	}	1104	}
881	}	1105	}
882		1106
		1107	for(itv = vAngular.begin(); itv != vAngular.end(); ++itv) {
		1108	delete (*itv);
		1109	}
		1110
883	if( std::getenv("G4PHPTEST2") ) {	1111	if( std::getenv("G4PHPTEST2") ) {
884	G4cout << " G4ParticleHPContAngularPar::Merge ANGPAR1 " << G4endl; //GDEB	1112	G4cout << " G4ParticleHPContAngularPar::Merge ANGPAR1 " << G4endl; //GDEB
885	angpar1.Dump();	1113	copyAngpar1.Dump();
886	G4cout << " G4ParticleHPContAngularPar::Merge ANGPAR2 " << G4endl;	1114	G4cout << " G4ParticleHPContAngularPar::Merge ANGPAR2 " << G4endl;
887	angpar2.Dump();	1115	copyAngpar2.Dump();
888	G4cout << " G4ParticleHPContAngularPar::Merge ANGPARNEW " << G4endl;	1116	G4cout << " G4ParticleHPContAngularPar::Merge ANGPARNEW " << G4endl;
889	Dump();	1117	Dump();
890	}	1118	}
891	}	1119	}
892		1120
893	void G4ParticleHPContAngularPar::Dump()	1121	void G4ParticleHPContAngularPar::Dump() const
894	{	1122	{
895	G4cout << theEnergy << " " << nEnergies << " " << nDiscreteEnergies << " " << nAngularParameters << G4endl;	1123	G4cout << theEnergy << " " << nEnergies << " " << nDiscreteEnergies << " " << nAngularParameters << G4endl;
896		1124

Lines 91-100 Link Here

(-)source-orig/processes/hadronic/models/particle_hp/src/G4ParticleHPContEnergyAngular.cc (-1 / +18 lines)
91	fCacheAngular.Get()->SetTarget(GetTarget());	91	fCacheAngular.Get()->SetTarget(GetTarget());
92	fCacheAngular.Get()->SetTargetCode(theTargetCode);	92	fCacheAngular.Get()->SetTargetCode(theTargetCode);
93	fCacheAngular.Get()->SetPrimary(GetProjectileRP());	93	fCacheAngular.Get()->SetPrimary(GetProjectileRP());
		94
94	result = fCacheAngular.Get()->Sample(anEnergy, massCode, targetMass,	95	result = fCacheAngular.Get()->Sample(anEnergy, massCode, targetMass,
95	theAngularRep, theInterpolation);	96	theAngularRep, theInterpolation);
96	currentMeanEnergy.Put( fCacheAngular.Get()->MeanEnergyOfThisInteraction() );	97	currentMeanEnergy.Put( fCacheAngular.Get()->MeanEnergyOfThisInteraction() );
97	fCacheAngular.Get()->ClearHistories();	98	// SM 12.Feb 2022
		99	//
		100	// do not call ClearHistories() here because multiple photons need to be produced and they need
		101	// the history to get the remaining energy correct
		102	//
		103	// fCacheAngular.Get()->ClearHistories();
		104	//
		105	// SM 12.Feb 2022
98		106
99	// delete fAngular; //fix end	107	// delete fAngular; //fix end
100		108
Lines 130-133 Link Here
130	for ( G4int i = 0 ; i< nEnergy ; i++ )	138	for ( G4int i = 0 ; i< nEnergy ; i++ )
131	theAngular[i].ClearHistories();	139	theAngular[i].ClearHistories();
132	}	140	}
		141	// SM 12.Feb 2022
		142	//
		143	// added fCacheAngular ClearHistories() - this is the one actually used!
		144	// Maybe theAngular does not even need ClearHistories()?
		145	//
		146	if (fCacheAngular.Get() != 0) {
		147	fCacheAngular.Get()->ClearHistories();
		148	}
		149	// SM 12.Feb 2022
133	}	150	}

Return to problem 2468

Lines 79-84 Link Here

(-)source-orig/processes/hadronic/models/particle_hp/include/G4ParticleHPContAngularPar.hh (-3 / +28 lines)
79	adjustResult = true;	79	adjustResult = true;
80	}	80	}
81		81
		82	G4ParticleHPContAngularPar(G4ParticleHPContAngularPar & val)
		83	{
		84	theEnergy = val.theEnergy;
		85	nEnergies = val.nEnergies;
		86	nDiscreteEnergies = val.nDiscreteEnergies;
		87	nAngularParameters= val.nAngularParameters;
		88	theProjectile = val.theProjectile;
		89	theManager = val.theManager;
		90	theInt = val.theInt;
		91	adjustResult = val.adjustResult;
		92	theMinEner = val.theMinEner;
		93	theMaxEner = val.theMaxEner;
		94	theEnergiesTransformed = val.theEnergiesTransformed;
		95	theDiscreteEnergies = val.theDiscreteEnergies;
		96	theDiscreteEnergiesOwn = val.theDiscreteEnergiesOwn;
		97	fCache.Put(0);
		98	theAngular = new G4ParticleHPList[nEnergies];
		99	for(int ie=0;ie<nEnergies;++ie) {
		100	theAngular[ie].SetLabel(val.theAngular[ie].GetLabel());
		101	for(int ip=0;ip<nAngularParameters;++ip) {
		102	theAngular[ie].SetValue(ip,val.theAngular[ie].GetValue(ip));
		103	}
		104	}
		105	}
		106
82	G4ParticleHPContAngularPar(G4ParticleDefinition* projectile);	107	G4ParticleHPContAngularPar(G4ParticleDefinition* projectile);
83		108
84	~G4ParticleHPContAngularPar()	109	~G4ParticleHPContAngularPar()
Lines 92-98 Link Here
92	G4ReactionProduct* Sample(G4double anEnergy, G4double massCode, G4double mass,	117	G4ReactionProduct* Sample(G4double anEnergy, G4double massCode, G4double mass,
93	G4int angularRep, G4int interpol);	118	G4int angularRep, G4int interpol);
94		119
95	G4double GetEnergy()	120	G4double GetEnergy() const
96	{	121	{
97	if( std::getenv("G4PHPTEST") )	122	if( std::getenv("G4PHPTEST") )
98	G4cout << this << " G4ParticleHPContAngularPar::GetEnergy "	123	G4cout << this << " G4ParticleHPContAngularPar::GetEnergy "
Lines 125-131 Link Here
125	G4ParticleHPContAngularPar & store2);	150	G4ParticleHPContAngularPar & store2);
126	// NOTE: this interpolates legendre coefficients	151	// NOTE: this interpolates legendre coefficients
127		152
128	void PrepareTableInterpolation(const G4ParticleHPContAngularPar* angularPrev);	153	void PrepareTableInterpolation(); // SM 26 Apr 2022 const G4ParticleHPContAngularPar* angularPrev);
129		154
130	G4double MeanEnergyOfThisInteraction()	155	G4double MeanEnergyOfThisInteraction()
131	{	156	{
Lines 182-188 Link Here
182	fCache.Get()->fresh = true;	207	fCache.Get()->fresh = true;
183	}	208	}
184		209
185	void Dump();	210	void Dump() const;
186		211
187	private:	212	private:
188		213