Attachment #749 for problem #2468

Lines 59-64 Link Here

(-)../source-orig/processes/hadronic/models/particle_hp/src/G4ParticleHPContAngularPar.cc (-26 / +102 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 194-200 Link Here
194	for ( G4int j = 0 ; j < nDiscreteEnergies ; j++ )	195	for ( G4int j = 0 ; j < nDiscreteEnergies ; j++ )
195	{	196	{
196	G4double delta = 0.0;	197	G4double delta = 0.0;
197	if ( theAngular[j].GetLabel() <= fCache.Get()->remaining_energy ) delta = theAngular[i].GetValue(0);	198	if ( theAngular[j].GetLabel() <= fCache.Get()->remaining_energy ) delta = theAngular[j].GetValue(0);
198	running[j+1] = running[j] + delta;	199	running[j+1] = running[j] + delta;
199	}	200	}
200	G4double tot_prob_DIS = running[ nDiscreteEnergies ];	201	G4double tot_prob_DIS = running[ nDiscreteEnergies ];
Lines 739-796 Link Here
739	G4ParticleHPContAngularPar & angpar2)	740	G4ParticleHPContAngularPar & angpar2)
740	{	741	{
741	G4int ie,ie1,ie2, ie1Prev, ie2Prev;	742	G4int ie,ie1,ie2, ie1Prev, ie2Prev;
		743	// SM 12.Feb 2022
		744	//
		745	// only re-build the interpolation table if we have a new interaction. For several subsequent
		746	// samplings of final state particles in the same interaction we should use the existing table.
		747	if ( fCache.Get()->fresh != true ) return;
742	nAngularParameters = angpar1.nAngularParameters;	748	nAngularParameters = angpar1.nAngularParameters;
743	theManager = angpar1.theManager;	749	theManager = angpar1.theManager;
744	theEnergy = anEnergy;	750	theEnergy = anEnergy;
745		751	theMinEner = std::min(angpar1.theMinEner,angpar2.theMinEner);
746	nDiscreteEnergies = theDiscreteEnergies.size();	752	theMaxEner = std::max(angpar1.theMaxEner,angpar2.theMaxEner);
747	std::set<G4double>::const_iterator itede;	753	// SM 11. Feb 2022
		754	//
		755	// re-wrote most of this. The two discrete sets have to be merged. A vector holds the temporary
		756	// data to be copied to the array in the end. Since the G4ParticleHPList class contains pointers
		757	// one can not simply assign elements of this type. Each member needs to call the explicit Set()
		758	// method instead.
		759	//
		760	// First average probabilities for those lines that are in both sets
748	std::map<G4double,G4int> discEnerOwn1 = angpar1.GetDiscreteEnergiesOwn();	761	std::map<G4double,G4int> discEnerOwn1 = angpar1.GetDiscreteEnergiesOwn();
749	std::map<G4double,G4int> discEnerOwn2 = angpar2.GetDiscreteEnergiesOwn();	762	std::map<G4double,G4int> discEnerOwn2 = angpar2.GetDiscreteEnergiesOwn();
750	std::map<G4double,G4int>::const_iterator itedeo;	763	std::map<G4double,G4int>::const_iterator itedeo1;
751	ie = 0;	764	std::map<G4double,G4int>::const_iterator itedeo2;
752	for( itede = theDiscreteEnergies.begin(); itede != theDiscreteEnergies.end(); itede++, ie++ ) {	765	std::vector<G4ParticleHPList> vAngular(discEnerOwn1.size());
753	G4double discEner = *itede;	766	for( itedeo1 = discEnerOwn1.begin(); itedeo1 != discEnerOwn1.end(); itedeo1++) {
754	itedeo = discEnerOwn1.find(discEner);	767	G4double discEner = itedeo1->first;
755	if( itedeo == discEnerOwn1.end() ) {	768	ie1 = itedeo1->second;
756	ie1 = 0;	769	itedeo2 = discEnerOwn2.find(discEner);
757	} else {	770	if( itedeo2 == discEnerOwn2.end() ) {
758	ie1 = -1;
759	}
760	itedeo = discEnerOwn2.find(discEner);
761	if( itedeo == discEnerOwn2.end() ) {
762	ie2 = 0;
763	} else {
764	ie2 = -1;	771	ie2 = -1;
		772	} else {
		773	ie2 = itedeo2->second;
765	}	774	}
766		775	vAngular[ie1].SetLabel(angpar1.theAngular[ie1].GetLabel());
767	theAngular[ie].SetLabel(discEner);
768	G4double val1, val2;	776	G4double val1, val2;
769	for(G4int ip=0; ip<nAngularParameters; ip++) {	777	for(G4int ip=0; ip<nAngularParameters; ip++) {
770	if( ie1 != -1 ) {
771	val1 = angpar1.theAngular[ie1].GetValue(ip);	778	val1 = angpar1.theAngular[ie1].GetValue(ip);
772	} else {
773	val1 = 0.;
774	}
775	if( ie2 != -1 ) {	779	if( ie2 != -1 ) {
776	val2 = angpar2.theAngular[ie2].GetValue(ip);	780	val2 = angpar2.theAngular[ie2].GetValue(ip);
777	} else {	781	} else {
778	val2 = 0.;	782	val2 = 0.;
779	}	783	}
		784	G4double value = theInt.Interpolate(aScheme, anEnergy,
		785	angpar1.theEnergy, angpar2.theEnergy,
		786	val1,
		787	val2);
		788	if( std::getenv("G4PHPTEST2") ) G4cout << ie1 << " " << ip << " G4ParticleHPContAngularPar::Merge DiscreteEnergies val1 " << val1 << " val2 " << val2 << " value " << value << G4endl; //GDEB
780		789
		790	vAngular[ie1].SetValue(ip, value);
		791	}
		792	}
		793	// add the ones in set2 but not in set1
		794	std::vector<G4ParticleHPList>::iterator itv;
		795	for( itedeo2 = discEnerOwn2.begin(); itedeo2 != discEnerOwn2.end(); itedeo2++) {
		796	G4double discEner = itedeo2->first;
		797	ie2=itedeo2->second;
		798	itedeo1 = discEnerOwn1.find(discEner);
		799	if( itedeo1 == discEnerOwn1.end() ) {
		800	// not yet in list
		801	// find position to insert
		802	bool isInserted=false;
		803	ie=0;
		804	for(itv=vAngular.begin();itv!=vAngular.end();itv++,ie++) {
		805	if ( discEner > itv->GetLabel() ) {
		806	G4ParticleHPList newEntry;
		807	newEntry.SetLabel(angpar2.theAngular[ie2].GetLabel());
		808	for(int ip=0;ip<nAngularParameters;ip++) {
		809	newEntry.SetValue(ip,angpar2.theAngular[ie2].GetValue(ip));
		810	}
		811	vAngular.insert(itv,newEntry);
		812	isInserted=true;
		813	break;
		814	}
		815	}
		816	if (!isInserted) {
		817	ie=vAngular.size();
		818	G4ParticleHPList newEntry;
		819	newEntry.SetLabel(angpar2.theAngular[ie2].GetLabel());
		820	for(int ip=0;ip<nAngularParameters;ip++) {
		821	newEntry.SetValue(ip,angpar2.theAngular[ie2].GetValue(ip));
		822	}
		823	vAngular.push_back(newEntry);
		824	isInserted=true;
		825	}
		826	G4double val1, val2;
		827	for(G4int ip=0; ip<nAngularParameters; ip++) {
		828	val1 = 0;
		829	val2 = angpar2.theAngular[ie2].GetValue(ip);
781	G4double value = theInt.Interpolate(aScheme, anEnergy,	830	G4double value = theInt.Interpolate(aScheme, anEnergy,
782	angpar1.theEnergy, angpar2.theEnergy,	831	angpar1.theEnergy, angpar2.theEnergy,
783	val1,	832	val1,
784	val2);	833	val2);
785	if( std::getenv("G4PHPTEST2") ) G4cout << ie << " " << ip << " G4ParticleHPContAngularPar::Merge DiscreteEnergies val1 " << val1 << " val2 " << val2 << " value " << value << G4endl; //GDEB	834	if( std::getenv("G4PHPTEST2") ) G4cout << ie << " " << ip << " G4ParticleHPContAngularPar::Merge DiscreteEnergies val1 " << val1 << " val2 " << val2 << " value " << value << G4endl; //GDEB
786		835
787	theAngular[ie].SetValue(ip, value);	836	vAngular[ie].SetValue(ip, value);
		837	}
788	}	838	}
789	}	839	}
790		840
		841	// store new discrete list
		842	nDiscreteEnergies = vAngular.size();
791	if(theAngular != 0) delete [] theAngular;	843	if(theAngular != 0) delete [] theAngular;
		844	theAngular = new G4ParticleHPList [nDiscreteEnergies];
		845	theDiscreteEnergiesOwn.clear();
		846	theDiscreteEnergies.clear();
		847	for(ie=0; ie<nDiscreteEnergies; ie++) {
		848	theAngular[ie].SetLabel(vAngular[ie].GetLabel());
		849	for(int ip=0; ip<nAngularParameters; ip++) {
		850	theAngular[ie].SetValue(ip,vAngular[ie].GetValue(ip));
		851	}
		852	theDiscreteEnergiesOwn[theAngular[ie].GetLabel()] = ie;
		853	theDiscreteEnergies.insert(theAngular[ie].GetLabel());
		854	}
		855
792	nEnergies = nDiscreteEnergies + angpar2.GetNEnergiesTransformed();	856	nEnergies = nDiscreteEnergies + angpar2.GetNEnergiesTransformed();
793	theAngular = new G4ParticleHPList [nEnergies];	857	if ( nEnergies > nDiscreteEnergies ) {
		858	G4ParticleHPList * theNewAngular = new G4ParticleHPList [nEnergies];
		859	if (theAngular !=0 ) {
		860	for(ie=0; ie<nDiscreteEnergies; ie++) {
		861	theNewAngular[ie].SetLabel(theAngular[ie].GetLabel());
		862	for(int ip=0; ip<nAngularParameters; ip++) {
		863	theNewAngular[ie].SetValue(ip,theAngular[ie].GetValue(ip));
		864	}
		865	}
		866	delete [] theAngular;
		867	}
		868	theAngular = theNewAngular;
		869	}
794		870
795	//---- Get minimum and maximum energy interpolating	871	//---- Get minimum and maximum energy interpolating
796	theMinEner = angpar1.GetMinEner() + (theEnergy-angpar1.GetEnergy()) * (angpar2.GetMinEner()-angpar1.GetMinEner())/(angpar2.GetEnergy()-angpar1.GetEnergy());	872	theMinEner = angpar1.GetMinEner() + (theEnergy-angpar1.GetEnergy()) * (angpar2.GetMinEner()-angpar1.GetMinEner())/(angpar2.GetEnergy()-angpar1.GetEnergy());

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