Attachment #764 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 "

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

  void Dump() const;

private:
  

Lines 59-64 Link Here

(-)source-orig/processes/hadronic/models/particle_hp/src/G4ParticleHPContAngularPar.cc (-72 / +181 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 733-739 Link Here
733	//- if( getenv("G4PHPTEST2") ) G4cout << this << " G4ParticleHPContAngularPar::PrepareTableInterpolation theEnergiesTransformed2 " << enerT << G4endl; //GDEB	755	//- if( getenv("G4PHPTEST2") ) G4cout << this << " G4ParticleHPContAngularPar::PrepareTableInterpolation theEnergiesTransformed2 " << enerT << G4endl; //GDEB
734	}	756	}
735	// add the maximum energy	757	// add the maximum energy
736	theEnergiesTransformed.insert(1.);	758	//theEnergiesTransformed.insert(1.);
737		759
738	}	760	}
739		761
Lines 742-825 Link Here
742	G4ParticleHPContAngularPar & angpar2)	764	G4ParticleHPContAngularPar & angpar2)
743	{	765	{
744	G4int ie,ie1,ie2, ie1Prev, ie2Prev;	766	G4int ie,ie1,ie2, ie1Prev, ie2Prev;
745	nAngularParameters = angpar1.nAngularParameters;	767	// SM 12.Feb 2022
746	theManager = angpar1.theManager;	768	//
		769	// only re-build the interpolation table if we have a new interaction. For several subsequent
		770	// samplings of final state particles in the same interaction we should use the existing table.
		771	if ( fCache.Get()->fresh != true ) return;
		772
		773	// SM 10. Apr 2022
		774	//
		775	// make copies of angpar1 and angpar2. Since these are given by reference it can not
		776	// be excluded that one of them is this. Hence this code uses potentially the old this
		777	// for creating the new this - which leads to memory corruption if the old is not stored
		778	// as separarte object for lookup
		779	const G4ParticleHPContAngularPar copyAngpar1(angpar1),copyAngpar2(angpar2);
		780
		781	nAngularParameters = copyAngpar1.nAngularParameters;
		782	theManager = copyAngpar1.theManager;
747	theEnergy = anEnergy;	783	theEnergy = anEnergy;
748		784	theMinEner = std::min(copyAngpar1.theMinEner,copyAngpar2.theMinEner);
749	nDiscreteEnergies = theDiscreteEnergies.size();	785	theMaxEner = std::max(copyAngpar1.theMaxEner,copyAngpar2.theMaxEner);
750	std::set<G4double>::const_iterator itede;	786	// SM 11. Feb 2022
751	std::map<G4double,G4int> discEnerOwn1 = angpar1.GetDiscreteEnergiesOwn();	787	//
752	std::map<G4double,G4int> discEnerOwn2 = angpar2.GetDiscreteEnergiesOwn();	788	// 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;	789	// data to be copied to the array in the end. Since the G4ParticleHPList class contains pointers
754	ie = 0;	790	// 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++ ) {	791	// method instead.
756	G4double discEner = *itede;	792	//
757	itedeo = discEnerOwn1.find(discEner);	793	// First average probabilities for those lines that are in both sets
758	if( itedeo == discEnerOwn1.end() ) {	794	const std::map<G4double,G4int> discEnerOwn1 = copyAngpar1.GetDiscreteEnergiesOwn();
759	ie1 = 0;	795	const std::map<G4double,G4int> discEnerOwn2 = copyAngpar2.GetDiscreteEnergiesOwn();
760	} else {	796	std::map<G4double,G4int>::const_iterator itedeo1;
761	ie1 = -1;	797	std::map<G4double,G4int>::const_iterator itedeo2;
762	}	798	std::vector<G4ParticleHPList *> vAngular(discEnerOwn1.size());
763	itedeo = discEnerOwn2.find(discEner);	799	for( itedeo1 = discEnerOwn1.begin(); itedeo1 != discEnerOwn1.end(); ++itedeo1) {
764	if( itedeo == discEnerOwn2.end() ) {	800	G4double discEner = itedeo1->first;
765	ie2 = 0;	801	ie1 = itedeo1->second;
766	} else {	802	itedeo2 = discEnerOwn2.find(discEner);
		803	if( itedeo2 == discEnerOwn2.end() ) {
767	ie2 = -1;	804	ie2 = -1;
		805	} else {
		806	ie2 = itedeo2->second;
768	}	807	}
769		808	vAngular[ie1] = new G4ParticleHPList();
770	theAngular[ie].SetLabel(discEner);	809	vAngular[ie1]->SetLabel(copyAngpar1.theAngular[ie1].GetLabel());
771	G4double val1, val2;	810	G4double val1, val2;
772	for(G4int ip=0; ip<nAngularParameters; ip++) {	811	for(G4int ip=0; ip<nAngularParameters; ++ip) {
773	if( ie1 != -1 ) {	812	val1 = copyAngpar1.theAngular[ie1].GetValue(ip);
774	val1 = angpar1.theAngular[ie1].GetValue(ip);
775	} else {
776	val1 = 0.;
777	}
778	if( ie2 != -1 ) {	813	if( ie2 != -1 ) {
779	val2 = angpar2.theAngular[ie2].GetValue(ip);	814	val2 = copyAngpar2.theAngular[ie2].GetValue(ip);
780	} else {	815	} else {
781	val2 = 0.;	816	val2 = 0.;
782	}	817	}
		818	G4double value = theInt.Interpolate(aScheme, anEnergy,
		819	copyAngpar1.theEnergy, copyAngpar2.theEnergy,
		820	val1,
		821	val2);
		822	if( std::getenv("G4PHPTEST2") ) G4cout << ie1 << " " << ip << " G4ParticleHPContAngularPar::Merge DiscreteEnergies val1 " << val1 << " val2 " << val2 << " value " << value << G4endl; //GDEB
783		823
		824	vAngular[ie1]->SetValue(ip, value);
		825	}
		826	}
		827	// add the ones in set2 but not in set1
		828	std::vector<G4ParticleHPList *>::iterator itv;
		829	for( itedeo2 = discEnerOwn2.begin(); itedeo2 != discEnerOwn2.end(); ++itedeo2) {
		830	G4double discEner = itedeo2->first;
		831	ie2=itedeo2->second;
		832	bool notFound = true;
		833	for(itv = vAngular.begin(); itv != vAngular.end(); ++itv) {
		834	if ( discEner == (*itv)->GetLabel() ) {
		835	notFound = false;
		836	break;
		837	}
		838	}
		839	if( notFound ) {
		840	// not yet in list
		841	// find position to insert
		842	bool isInserted=false;
		843	ie=0;
		844	for(itv=vAngular.begin();itv!=vAngular.end();++itv,++ie) {
		845	if ( discEner > (*itv)->GetLabel() ) {
		846	itv=vAngular.insert(itv,new G4ParticleHPList);
		847	(*itv)->SetLabel(copyAngpar2.theAngular[ie2].GetLabel());
		848	isInserted=true;
		849	break;
		850	}
		851	}
		852	if (!isInserted) {
		853	ie=vAngular.size();
		854	vAngular.push_back(new G4ParticleHPList);
		855	vAngular[ie]->SetLabel(copyAngpar2.theAngular[ie2].GetLabel());
		856	isInserted=true;
		857	}
		858	G4double val1, val2;
		859	for(G4int ip=0; ip<nAngularParameters; ++ip) {
		860	val1 = 0;
		861	val2 = copyAngpar2.theAngular[ie2].GetValue(ip);
784	G4double value = theInt.Interpolate(aScheme, anEnergy,	862	G4double value = theInt.Interpolate(aScheme, anEnergy,
785	angpar1.theEnergy, angpar2.theEnergy,	863	copyAngpar1.theEnergy, copyAngpar2.theEnergy,
786	val1,	864	val1,
787	val2);	865	val2);
788	if( std::getenv("G4PHPTEST2") ) G4cout << ie << " " << ip << " G4ParticleHPContAngularPar::Merge DiscreteEnergies val1 " << val1 << " val2 " << val2 << " value " << value << G4endl; //GDEB	866	if( std::getenv("G4PHPTEST2") ) G4cout << ie << " " << ip << " G4ParticleHPContAngularPar::Merge DiscreteEnergies val1 " << val1 << " val2 " << val2 << " value " << value << G4endl; //GDEB
789		867
790	theAngular[ie].SetValue(ip, value);	868	vAngular[ie]->SetValue(ip, value);
		869	}
791	}	870	}
792	}	871	}
793		872
		873	// store new discrete list
		874	nDiscreteEnergies = vAngular.size();
794	if(theAngular != 0) delete [] theAngular;	875	if(theAngular != 0) delete [] theAngular;
795	nEnergies = nDiscreteEnergies + angpar2.GetNEnergiesTransformed();	876	theAngular = new G4ParticleHPList [nDiscreteEnergies];
796	theAngular = new G4ParticleHPList [nEnergies];	877	theDiscreteEnergiesOwn.clear();
		878	theDiscreteEnergies.clear();
		879	for(ie=0; ie<nDiscreteEnergies; ++ie) {
		880	theAngular[ie].SetLabel(vAngular[ie]->GetLabel());
		881	for(int ip=0; ip<nAngularParameters; ++ip) {
		882	theAngular[ie].SetValue(ip,vAngular[ie]->GetValue(ip));
		883	}
		884	theDiscreteEnergiesOwn[theAngular[ie].GetLabel()] = ie;
		885	theDiscreteEnergies.insert(theAngular[ie].GetLabel());
		886	}
		887
		888	nEnergies = nDiscreteEnergies + copyAngpar2.GetNEnergiesTransformed();
		889	if ( nEnergies > nDiscreteEnergies ) {
		890	G4ParticleHPList * theNewAngular = new G4ParticleHPList [nEnergies];
		891	if (theAngular !=0 ) {
		892	for(ie=0; ie<nDiscreteEnergies; ++ie) {
		893	theNewAngular[ie].SetLabel(theAngular[ie].GetLabel());
		894	for(int ip=0; ip<nAngularParameters; ++ip) {
		895	theNewAngular[ie].SetValue(ip,theAngular[ie].GetValue(ip));
		896	}
		897	}
		898	delete [] theAngular;
		899	}
		900	theAngular = theNewAngular;
		901	}
797		902
798	//---- Get minimum and maximum energy interpolating	903	//---- Get minimum and maximum energy interpolating
799	theMinEner = angpar1.GetMinEner() + (theEnergy-angpar1.GetEnergy()) * (angpar2.GetMinEner()-angpar1.GetMinEner())/(angpar2.GetEnergy()-angpar1.GetEnergy());	904	theMinEner = copyAngpar1.GetMinEner() + (theEnergy-copyAngpar1.GetEnergy()) * (copyAngpar2.GetMinEner()-copyAngpar1.GetMinEner())/(copyAngpar2.GetEnergy()-copyAngpar1.GetEnergy());
800	theMaxEner = angpar1.GetMaxEner() + (theEnergy-angpar1.GetEnergy()) * (angpar2.GetMaxEner()-angpar1.GetMaxEner())/(angpar2.GetEnergy()-angpar1.GetEnergy());	905	theMaxEner = copyAngpar1.GetMaxEner() + (theEnergy-copyAngpar1.GetEnergy()) * (copyAngpar2.GetMaxEner()-copyAngpar1.GetMaxEner())/(copyAngpar2.GetEnergy()-copyAngpar1.GetEnergy());
801		906
802	if( std::getenv("G4PHPTEST2") ) G4cout << " G4ParticleHPContAngularPar::Merge E " << anEnergy << " minmax " << theMinEner << " " << theMaxEner << G4endl; //GDEB	907	if( std::getenv("G4PHPTEST2") ) G4cout << " G4ParticleHPContAngularPar::Merge E " << anEnergy << " minmax " << theMinEner << " " << theMaxEner << G4endl; //GDEB
803		908
804	//--- Loop to energies of new set	909	//--- Loop to energies of new set
805	std::set<G4double> energiesTransformed = angpar2.GetEnergiesTransformed();	910	std::set<G4double> energiesTransformed = copyAngpar2.GetEnergiesTransformed();
806	std::set<G4double>::const_iterator iteet = energiesTransformed.begin();	911	std::set<G4double>::const_iterator iteet = energiesTransformed.begin();
807	G4int nEnergies1 = angpar1.GetNEnergies();	912	G4int nEnergies1 = copyAngpar1.GetNEnergies();
808	G4int nDiscreteEnergies1 = angpar1.GetNDiscreteEnergies();	913	G4int nDiscreteEnergies1 = copyAngpar1.GetNDiscreteEnergies();
809	G4double minEner1 = angpar1.GetMinEner();	914	G4double minEner1 = copyAngpar1.GetMinEner();
810	G4double maxEner1 = angpar1.GetMaxEner();	915	G4double maxEner1 = copyAngpar1.GetMaxEner();
811	G4int nEnergies2 = angpar2.GetNEnergies();	916	G4int nEnergies2 = copyAngpar2.GetNEnergies();
812	G4int nDiscreteEnergies2 = angpar2.GetNDiscreteEnergies();	917	G4int nDiscreteEnergies2 = copyAngpar2.GetNDiscreteEnergies();
813	G4double minEner2 = angpar2.GetMinEner();	918	G4double minEner2 = copyAngpar2.GetMinEner();
814	G4double maxEner2 = angpar2.GetMaxEner();	919	G4double maxEner2 = copyAngpar2.GetMaxEner();
815	for(ie=nDiscreteEnergies; ie<nEnergies; ie++,iteet++) {	920	for(ie=nDiscreteEnergies; ie<nEnergies; ++ie,++iteet) {
816	G4double eT = (*iteet);	921	G4double eT = (*iteet);
817		922
818	//--- Use eT1 = eT: Get energy and parameters of angpar1 for this eT	923	//--- Use eT1 = eT: Get energy and parameters of copyAngpar1 for this eT
819	G4double e1 = (maxEner1-minEner1) * eT + minEner1;	924	G4double e1 = (maxEner1-minEner1) * eT + minEner1;
820	//----- Get parameter value corresponding to this e1	925	//----- Get parameter value corresponding to this e1
821	for(ie1=nDiscreteEnergies1; ie1<nEnergies1; ie1++) {	926	for(ie1=nDiscreteEnergies1; ie1<nEnergies1; ++ie1) {
822	if( (angpar1.theAngular[ie1].GetLabel() - e1) > 1.E-10*e1 ) break;	927	if( (copyAngpar1.theAngular[ie1].GetLabel() - e1) > 1.E-10*e1 ) break;
823	}	928	}
824	ie1Prev = ie1 - 1;	929	ie1Prev = ie1 - 1;
825	if( ie1 == 0 ) ie1Prev++;	930	if( ie1 == 0 ) ie1Prev++;
Lines 827-838 Link Here
827	ie1--;	932	ie1--;
828	ie1Prev = ie1;	933	ie1Prev = ie1;
829	}	934	}
830	//--- Use eT2 = eT: Get energy and parameters of angpar2 for this eT	935	//--- Use eT2 = eT: Get energy and parameters of copyAngpar2 for this eT
831	G4double e2 = (maxEner2-minEner2) * eT + minEner2;	936	G4double e2 = (maxEner2-minEner2) * eT + minEner2;
832	//----- Get parameter value corresponding to this e2	937	//----- Get parameter value corresponding to this e2
833	for(ie2=nDiscreteEnergies2; ie2<nEnergies2; ie2++) {	938	for(ie2=nDiscreteEnergies2; ie2<nEnergies2; ++ie2) {
834	// G4cout << " GET IE2 " << ie2 << " - " << angpar2.theAngular[ie2].GetLabel() - e2 << " " << angpar2.theAngular[ie2].GetLabel() << " " << e2 <<G4endl;	939	// 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;	940	if( (copyAngpar2.theAngular[ie2].GetLabel() - e2) > 1.E-10*e2 ) break;
836	}	941	}
837	ie2Prev = ie2 - 1;	942	ie2Prev = ie2 - 1;
838	if( ie2 == 0 ) ie2Prev++;	943	if( ie2 == 0 ) ie2Prev++;
Lines 847-868 Link Here
847		952
848	theAngular[ie].SetLabel(eN);	953	theAngular[ie].SetLabel(eN);
849		954
850	for(G4int ip=0; ip<nAngularParameters; ip++) {	955	for(G4int ip=0; ip<nAngularParameters; ++ip) {
851	G4double val1 = theInt.Interpolate2(theManager.GetScheme(ie),	956	G4double val1 = theInt.Interpolate2(theManager.GetScheme(ie),
852	e1,	957	e1,
853	angpar1.theAngular[ie1Prev].GetLabel(),	958	copyAngpar1.theAngular[ie1Prev].GetLabel(),
854	angpar1.theAngular[ie1].GetLabel(),	959	copyAngpar1.theAngular[ie1].GetLabel(),
855	angpar1.theAngular[ie1Prev].GetValue(ip),	960	copyAngpar1.theAngular[ie1Prev].GetValue(ip),
856	angpar1.theAngular[ie1].GetValue(ip)) * (maxEner1-minEner1);	961	copyAngpar1.theAngular[ie1].GetValue(ip)) * (maxEner1-minEner1);
857	G4double val2 = theInt.Interpolate2(theManager.GetScheme(ie),	962	G4double val2 = theInt.Interpolate2(theManager.GetScheme(ie),
858	e2,	963	e2,
859	angpar2.theAngular[ie2Prev].GetLabel(),	964	copyAngpar2.theAngular[ie2Prev].GetLabel(),
860	angpar2.theAngular[ie2].GetLabel(),	965	copyAngpar2.theAngular[ie2].GetLabel(),
861	angpar2.theAngular[ie2Prev].GetValue(ip),	966	copyAngpar2.theAngular[ie2Prev].GetValue(ip),
862	angpar2.theAngular[ie2].GetValue(ip)) * (maxEner2-minEner2);	967	copyAngpar2.theAngular[ie2].GetValue(ip)) * (maxEner2-minEner2);
863		968
864	G4double value = theInt.Interpolate(aScheme, anEnergy,	969	G4double value = theInt.Interpolate(aScheme, anEnergy,
865	angpar1.theEnergy, angpar2.theEnergy,	970	copyAngpar1.theEnergy, copyAngpar2.theEnergy,
866	val1,	971	val1,
867	val2);	972	val2);
868	//value /= (theMaxEner-theMinEner);	973	//value /= (theMaxEner-theMinEner);
Lines 872-896 Link Here
872	throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPContAngularPar::PrepareTableInterpolation theMaxEner == theMinEner and value != 0.");	977	throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPContAngularPar::PrepareTableInterpolation theMaxEner == theMinEner and value != 0.");
873	}	978	}
874	if( std::getenv("G4PHPTEST2") ) G4cout << ie << " " << ip << " G4ParticleHPContAngularPar::Merge val1 " << val1 << " val2 " << val2 << " value " << value << G4endl; //GDEB	979	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);	980	//- val1 = copyAngpar1.theAngular[ie1-1].GetValue(ip) * (maxEner1-minEner1);
876	//- val2 = angpar2.theAngular[ie2-1].GetValue(ip) * (maxEner2-minEner2);	981	//- 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	982	//- if( getenv("G4PHPTEST2") ) G4cout << ie << " " << ip << " G4ParticleHPContAngularPar::MergeOLD val1 " << val1 << " val2 " << val2 << " value " << value << G4endl; //GDEB
878		983
879	theAngular[ie].SetValue(ip, value);	984	theAngular[ie].SetValue(ip, value);
880	}	985	}
881	}	986	}
882		987
		988	for(itv = vAngular.begin(); itv != vAngular.end(); ++itv) {
		989	delete (*itv);
		990	}
		991
883	if( std::getenv("G4PHPTEST2") ) {	992	if( std::getenv("G4PHPTEST2") ) {
884	G4cout << " G4ParticleHPContAngularPar::Merge ANGPAR1 " << G4endl; //GDEB	993	G4cout << " G4ParticleHPContAngularPar::Merge ANGPAR1 " << G4endl; //GDEB
885	angpar1.Dump();	994	copyAngpar1.Dump();
886	G4cout << " G4ParticleHPContAngularPar::Merge ANGPAR2 " << G4endl;	995	G4cout << " G4ParticleHPContAngularPar::Merge ANGPAR2 " << G4endl;
887	angpar2.Dump();	996	copyAngpar2.Dump();
888	G4cout << " G4ParticleHPContAngularPar::Merge ANGPARNEW " << G4endl;	997	G4cout << " G4ParticleHPContAngularPar::Merge ANGPARNEW " << G4endl;
889	Dump();	998	Dump();
890	}	999	}
891	}	1000	}
892		1001
893	void G4ParticleHPContAngularPar::Dump()	1002	void G4ParticleHPContAngularPar::Dump() const
894	{	1003	{
895	G4cout << theEnergy << " " << nEnergies << " " << nDiscreteEnergies << " " << nAngularParameters << G4endl;	1004	G4cout << theEnergy << " " << nEnergies << " " << nDiscreteEnergies << " " << nAngularParameters << G4endl;
896		1005

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 (-2 / +27 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 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