Attachment #755 for problem #2468




#include "G4IonTable.hh"
#include "G4ParticleHPManager.hh"
#include <set>
#include <vector>
 
G4ParticleHPContAngularPar::G4ParticleHPContAngularPar( G4ParticleDefinition* projectile )
{  

         for ( G4int j = 0 ; j < nDiscreteEnergies ; j++ ) 
         {
            G4double delta = 0.0;
            if ( theAngular[j].GetLabel() <= fCache.Get()->remaining_energy ) delta = theAngular[j].GetValue(0);
            running[j+1] = running[j] + delta;
         }
         G4double tot_prob_DIS = running[ nDiscreteEnergies ];

        }

         //TK080711
	 // SM 30 March 2022
	 // the remaining energy needs to be lowered by the photon
	 // energy in *any* case. Otherwise additional photons with
	 // too high energy will be produced - therefore the
	 // adjustResult condition has been removed
	 //
	 //if( adjustResult )  fCache.Get()->remaining_energy -= fsEnergy;
	 fCache.Get()->remaining_energy -= fsEnergy;
         //TK080711

         //080801b

         delete [] running;

         //080714
	 // SM 30 March 2022
	 // the remaining energy needs to be lowered by the photon
	 // energy in *any* case. Otherwise additional photons with
	 // too high energy will be produced - therefore the
	 // adjustResult condition has been removed
	 // 
	 //if( adjustResult )  fCache.Get()->remaining_energy -= fsEnergy;
	 fCache.Get()->remaining_energy -= fsEnergy;
         //080714
      }
   }

							G4ParticleHPContAngularPar & angpar2) 
{
  G4int ie,ie1,ie2, ie1Prev, ie2Prev;
  // SM 12.Feb 2022
  //
  // only re-build the interpolation table if we have a new interaction. For several subsequent
  // samplings of final state particles in the same interaction we should use the existing table.
  if ( fCache.Get()->fresh != true ) return;
  nAngularParameters = angpar1.nAngularParameters;
  theManager = angpar1.theManager;
  theEnergy = anEnergy;
  theMinEner = std::min(angpar1.theMinEner,angpar2.theMinEner);
  theMaxEner = std::max(angpar1.theMaxEner,angpar2.theMaxEner);
  // SM 11. Feb 2022
  //
  // re-wrote most of this. The two discrete sets have to be merged. A vector holds the temporary
  // data to be copied to the array in the end. Since the G4ParticleHPList class contains pointers
  // one can not simply assign elements of this type. Each member needs to call the explicit Set()
  // method instead.
  //
  // First average probabilities for those lines that are in both sets
  std::map<G4double,G4int> discEnerOwn1 = angpar1.GetDiscreteEnergiesOwn();
  std::map<G4double,G4int> discEnerOwn2 = angpar2.GetDiscreteEnergiesOwn();
  std::map<G4double,G4int>::const_iterator itedeo1;
  std::map<G4double,G4int>::const_iterator itedeo2;
  std::vector<G4ParticleHPList> vAngular(discEnerOwn1.size());
  for( itedeo1 = discEnerOwn1.begin(); itedeo1 != discEnerOwn1.end(); itedeo1++) {
    G4double discEner = itedeo1->first;
    ie1 = itedeo1->second;
    itedeo2 = discEnerOwn2.find(discEner);
    if( itedeo2 == discEnerOwn2.end() ) {
      ie1 = -1;
    }
    itedeo = discEnerOwn2.find(discEner);
    if( itedeo == discEnerOwn2.end() ) {
      ie2 = 0;
    } else {
      ie2 = -1;
    } else {
      ie2 = itedeo2->second;
    }
    vAngular[ie1].SetLabel(angpar1.theAngular[ie1].GetLabel());
    theAngular[ie].SetLabel(discEner);
    G4double val1, val2;
    for(G4int ip=0; ip<nAngularParameters; ip++) {
      if( ie1 != -1 ) {
	val1 = angpar1.theAngular[ie1].GetValue(ip);
      } else {
	val1 = 0.;
      }
      if( ie2 != -1 ) {
	val2 = angpar2.theAngular[ie2].GetValue(ip);
      } else {
	val2 = 0.;
      }
      G4double value = theInt.Interpolate(aScheme, anEnergy, 
					  angpar1.theEnergy, angpar2.theEnergy,
					  val1,
					  val2);
      if( std::getenv("G4PHPTEST2") ) G4cout << ie1 << " " << ip << " G4ParticleHPContAngularPar::Merge DiscreteEnergies  val1 " << val1 << " val2 " << val2 << " value " << value << G4endl; //GDEB
      
      vAngular[ie1].SetValue(ip, value);
    }
  }
  // add the ones in set2 but not in set1
  std::vector<G4ParticleHPList>::iterator itv;
  for( itedeo2 = discEnerOwn2.begin(); itedeo2 != discEnerOwn2.end(); itedeo2++) {
    G4double discEner = itedeo2->first;
    ie2=itedeo2->second;
    itedeo1 = discEnerOwn1.find(discEner);
    if( itedeo1 == discEnerOwn1.end() ) {
      // not yet in list
      // find position to insert
      bool isInserted=false;
      ie=0;
      for(itv=vAngular.begin();itv!=vAngular.end();itv++,ie++) {
	if ( discEner > itv->GetLabel() ) {
	  G4ParticleHPList newEntry;
	  newEntry.SetLabel(angpar2.theAngular[ie2].GetLabel());
	  for(int ip=0;ip<nAngularParameters;ip++) {
	    newEntry.SetValue(ip,angpar2.theAngular[ie2].GetValue(ip));
	  }
	  vAngular.insert(itv,newEntry);
	  isInserted=true;
	  break;
	}
      }
      if (!isInserted) {
	ie=vAngular.size();
	G4ParticleHPList newEntry;
	newEntry.SetLabel(angpar2.theAngular[ie2].GetLabel());
	for(int ip=0;ip<nAngularParameters;ip++) {
	  newEntry.SetValue(ip,angpar2.theAngular[ie2].GetValue(ip));
	}
	vAngular.push_back(newEntry);
	isInserted=true;
      }
      G4double val1, val2;
      for(G4int ip=0; ip<nAngularParameters; ip++) {
	val1 = 0;
	val2 = angpar2.theAngular[ie2].GetValue(ip);
      G4double value = theInt.Interpolate(aScheme, anEnergy, 
					  angpar1.theEnergy, angpar2.theEnergy,
					  val1,
					  val2);
      if( std::getenv("G4PHPTEST2") ) G4cout << ie << " " << ip << " G4ParticleHPContAngularPar::Merge DiscreteEnergies  val1 " << val1 << " val2 " << val2 << " value " << value << G4endl; //GDEB
      
	vAngular[ie].SetValue(ip, value);
      }
    }
  }
  
  // store new discrete list 
  nDiscreteEnergies = vAngular.size();
  if(theAngular != 0) delete [] theAngular;
  theAngular = new G4ParticleHPList [nDiscreteEnergies];
  theDiscreteEnergiesOwn.clear();
  theDiscreteEnergies.clear();
  for(ie=0; ie<nDiscreteEnergies; ie++) {
    theAngular[ie].SetLabel(vAngular[ie].GetLabel());
    for(int ip=0; ip<nAngularParameters; ip++) {
      theAngular[ie].SetValue(ip,vAngular[ie].GetValue(ip));
    }
    theDiscreteEnergiesOwn[theAngular[ie].GetLabel()] = ie;
    theDiscreteEnergies.insert(theAngular[ie].GetLabel());
  }
  
  nEnergies = nDiscreteEnergies + angpar2.GetNEnergiesTransformed();
  if ( nEnergies > nDiscreteEnergies ) {
    G4ParticleHPList * theNewAngular = new G4ParticleHPList [nEnergies];
    if (theAngular !=0 ) {
      for(ie=0; ie<nDiscreteEnergies; ie++) { 
	theNewAngular[ie].SetLabel(theAngular[ie].GetLabel());
	for(int ip=0; ip<nAngularParameters; ip++) {
	  theNewAngular[ie].SetValue(ip,theAngular[ie].GetValue(ip));
	}
      }
      delete [] theAngular;
    }
    theAngular = theNewAngular;
  }

  //---- Get minimum and maximum energy interpolating
  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

Lines 59-64 Link Here

(-)source-orig/processes/hadronic/models/particle_hp/src/G4ParticleHPContAngularPar.cc (-28 / +118 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 197-203 Link Here
197	for ( G4int j = 0 ; j < nDiscreteEnergies ; j++ )	198	for ( G4int j = 0 ; j < nDiscreteEnergies ; j++ )
198	{	199	{
199	G4double delta = 0.0;	200	G4double delta = 0.0;
200	if ( theAngular[j].GetLabel() <= fCache.Get()->remaining_energy ) delta = theAngular[i].GetValue(0);	201	if ( theAngular[j].GetLabel() <= fCache.Get()->remaining_energy ) delta = theAngular[j].GetValue(0);
201	running[j+1] = running[j] + delta;	202	running[j+1] = running[j] + delta;
202	}	203	}
203	G4double tot_prob_DIS = running[ nDiscreteEnergies ];	204	G4double tot_prob_DIS = running[ nDiscreteEnergies ];
Lines 327-333 Link Here
327	}	328	}
328		329
329	//TK080711	330	//TK080711
330	if( adjustResult ) fCache.Get()->remaining_energy -= fsEnergy;	331	// SM 30 March 2022
		332	// the remaining energy needs to be lowered by the photon
		333	// energy in any case. Otherwise additional photons with
		334	// too high energy will be produced - therefore the
		335	// adjustResult condition has been removed
		336	//
		337	//if( adjustResult ) fCache.Get()->remaining_energy -= fsEnergy;
		338	fCache.Get()->remaining_energy -= fsEnergy;
331	//TK080711	339	//TK080711
332		340
333	//080801b	341	//080801b
Lines 462-468 Link Here
462	delete [] running;	470	delete [] running;
463		471
464	//080714	472	//080714
465	if( adjustResult ) fCache.Get()->remaining_energy -= fsEnergy;	473	// SM 30 March 2022
		474	// the remaining energy needs to be lowered by the photon
		475	// energy in any case. Otherwise additional photons with
		476	// too high energy will be produced - therefore the
		477	// adjustResult condition has been removed
		478	//
		479	//if( adjustResult ) fCache.Get()->remaining_energy -= fsEnergy;
		480	fCache.Get()->remaining_energy -= fsEnergy;
466	//080714	481	//080714
467	}	482	}
468	}	483	}
Lines 742-799 Link Here
742	G4ParticleHPContAngularPar & angpar2)	757	G4ParticleHPContAngularPar & angpar2)
743	{	758	{
744	G4int ie,ie1,ie2, ie1Prev, ie2Prev;	759	G4int ie,ie1,ie2, ie1Prev, ie2Prev;
		760	// SM 12.Feb 2022
		761	//
		762	// only re-build the interpolation table if we have a new interaction. For several subsequent
		763	// samplings of final state particles in the same interaction we should use the existing table.
		764	if ( fCache.Get()->fresh != true ) return;
745	nAngularParameters = angpar1.nAngularParameters;	765	nAngularParameters = angpar1.nAngularParameters;
746	theManager = angpar1.theManager;	766	theManager = angpar1.theManager;
747	theEnergy = anEnergy;	767	theEnergy = anEnergy;
748		768	theMinEner = std::min(angpar1.theMinEner,angpar2.theMinEner);
749	nDiscreteEnergies = theDiscreteEnergies.size();	769	theMaxEner = std::max(angpar1.theMaxEner,angpar2.theMaxEner);
750	std::set<G4double>::const_iterator itede;	770	// SM 11. Feb 2022
		771	//
		772	// re-wrote most of this. The two discrete sets have to be merged. A vector holds the temporary
		773	// data to be copied to the array in the end. Since the G4ParticleHPList class contains pointers
		774	// one can not simply assign elements of this type. Each member needs to call the explicit Set()
		775	// method instead.
		776	//
		777	// First average probabilities for those lines that are in both sets
751	std::map<G4double,G4int> discEnerOwn1 = angpar1.GetDiscreteEnergiesOwn();	778	std::map<G4double,G4int> discEnerOwn1 = angpar1.GetDiscreteEnergiesOwn();
752	std::map<G4double,G4int> discEnerOwn2 = angpar2.GetDiscreteEnergiesOwn();	779	std::map<G4double,G4int> discEnerOwn2 = angpar2.GetDiscreteEnergiesOwn();
753	std::map<G4double,G4int>::const_iterator itedeo;	780	std::map<G4double,G4int>::const_iterator itedeo1;
754	ie = 0;	781	std::map<G4double,G4int>::const_iterator itedeo2;
755	for( itede = theDiscreteEnergies.begin(); itede != theDiscreteEnergies.end(); itede++, ie++ ) {	782	std::vector<G4ParticleHPList> vAngular(discEnerOwn1.size());
756	G4double discEner = *itede;	783	for( itedeo1 = discEnerOwn1.begin(); itedeo1 != discEnerOwn1.end(); itedeo1++) {
757	itedeo = discEnerOwn1.find(discEner);	784	G4double discEner = itedeo1->first;
758	if( itedeo == discEnerOwn1.end() ) {	785	ie1 = itedeo1->second;
759	ie1 = 0;	786	itedeo2 = discEnerOwn2.find(discEner);
760	} else {	787	if( itedeo2 == discEnerOwn2.end() ) {
761	ie1 = -1;
762	}
763	itedeo = discEnerOwn2.find(discEner);
764	if( itedeo == discEnerOwn2.end() ) {
765	ie2 = 0;
766	} else {
767	ie2 = -1;	788	ie2 = -1;
		789	} else {
		790	ie2 = itedeo2->second;
768	}	791	}
769		792	vAngular[ie1].SetLabel(angpar1.theAngular[ie1].GetLabel());
770	theAngular[ie].SetLabel(discEner);
771	G4double val1, val2;	793	G4double val1, val2;
772	for(G4int ip=0; ip<nAngularParameters; ip++) {	794	for(G4int ip=0; ip<nAngularParameters; ip++) {
773	if( ie1 != -1 ) {
774	val1 = angpar1.theAngular[ie1].GetValue(ip);	795	val1 = angpar1.theAngular[ie1].GetValue(ip);
775	} else {
776	val1 = 0.;
777	}
778	if( ie2 != -1 ) {	796	if( ie2 != -1 ) {
779	val2 = angpar2.theAngular[ie2].GetValue(ip);	797	val2 = angpar2.theAngular[ie2].GetValue(ip);
780	} else {	798	} else {
781	val2 = 0.;	799	val2 = 0.;
782	}	800	}
		801	G4double value = theInt.Interpolate(aScheme, anEnergy,
		802	angpar1.theEnergy, angpar2.theEnergy,
		803	val1,
		804	val2);
		805	if( std::getenv("G4PHPTEST2") ) G4cout << ie1 << " " << ip << " G4ParticleHPContAngularPar::Merge DiscreteEnergies val1 " << val1 << " val2 " << val2 << " value " << value << G4endl; //GDEB
783		806
		807	vAngular[ie1].SetValue(ip, value);
		808	}
		809	}
		810	// add the ones in set2 but not in set1
		811	std::vector<G4ParticleHPList>::iterator itv;
		812	for( itedeo2 = discEnerOwn2.begin(); itedeo2 != discEnerOwn2.end(); itedeo2++) {
		813	G4double discEner = itedeo2->first;
		814	ie2=itedeo2->second;
		815	itedeo1 = discEnerOwn1.find(discEner);
		816	if( itedeo1 == discEnerOwn1.end() ) {
		817	// not yet in list
		818	// find position to insert
		819	bool isInserted=false;
		820	ie=0;
		821	for(itv=vAngular.begin();itv!=vAngular.end();itv++,ie++) {
		822	if ( discEner > itv->GetLabel() ) {
		823	G4ParticleHPList newEntry;
		824	newEntry.SetLabel(angpar2.theAngular[ie2].GetLabel());
		825	for(int ip=0;ip<nAngularParameters;ip++) {
		826	newEntry.SetValue(ip,angpar2.theAngular[ie2].GetValue(ip));
		827	}
		828	vAngular.insert(itv,newEntry);
		829	isInserted=true;
		830	break;
		831	}
		832	}
		833	if (!isInserted) {
		834	ie=vAngular.size();
		835	G4ParticleHPList newEntry;
		836	newEntry.SetLabel(angpar2.theAngular[ie2].GetLabel());
		837	for(int ip=0;ip<nAngularParameters;ip++) {
		838	newEntry.SetValue(ip,angpar2.theAngular[ie2].GetValue(ip));
		839	}
		840	vAngular.push_back(newEntry);
		841	isInserted=true;
		842	}
		843	G4double val1, val2;
		844	for(G4int ip=0; ip<nAngularParameters; ip++) {
		845	val1 = 0;
		846	val2 = angpar2.theAngular[ie2].GetValue(ip);
784	G4double value = theInt.Interpolate(aScheme, anEnergy,	847	G4double value = theInt.Interpolate(aScheme, anEnergy,
785	angpar1.theEnergy, angpar2.theEnergy,	848	angpar1.theEnergy, angpar2.theEnergy,
786	val1,	849	val1,
787	val2);	850	val2);
788	if( std::getenv("G4PHPTEST2") ) G4cout << ie << " " << ip << " G4ParticleHPContAngularPar::Merge DiscreteEnergies val1 " << val1 << " val2 " << val2 << " value " << value << G4endl; //GDEB	851	if( std::getenv("G4PHPTEST2") ) G4cout << ie << " " << ip << " G4ParticleHPContAngularPar::Merge DiscreteEnergies val1 " << val1 << " val2 " << val2 << " value " << value << G4endl; //GDEB
789		852
790	theAngular[ie].SetValue(ip, value);	853	vAngular[ie].SetValue(ip, value);
		854	}
791	}	855	}
792	}	856	}
793		857
		858	// store new discrete list
		859	nDiscreteEnergies = vAngular.size();
794	if(theAngular != 0) delete [] theAngular;	860	if(theAngular != 0) delete [] theAngular;
		861	theAngular = new G4ParticleHPList [nDiscreteEnergies];
		862	theDiscreteEnergiesOwn.clear();
		863	theDiscreteEnergies.clear();
		864	for(ie=0; ie<nDiscreteEnergies; ie++) {
		865	theAngular[ie].SetLabel(vAngular[ie].GetLabel());
		866	for(int ip=0; ip<nAngularParameters; ip++) {
		867	theAngular[ie].SetValue(ip,vAngular[ie].GetValue(ip));
		868	}
		869	theDiscreteEnergiesOwn[theAngular[ie].GetLabel()] = ie;
		870	theDiscreteEnergies.insert(theAngular[ie].GetLabel());
		871	}
		872
795	nEnergies = nDiscreteEnergies + angpar2.GetNEnergiesTransformed();	873	nEnergies = nDiscreteEnergies + angpar2.GetNEnergiesTransformed();
796	theAngular = new G4ParticleHPList [nEnergies];	874	if ( nEnergies > nDiscreteEnergies ) {
		875	G4ParticleHPList * theNewAngular = new G4ParticleHPList [nEnergies];
		876	if (theAngular !=0 ) {
		877	for(ie=0; ie<nDiscreteEnergies; ie++) {
		878	theNewAngular[ie].SetLabel(theAngular[ie].GetLabel());
		879	for(int ip=0; ip<nAngularParameters; ip++) {
		880	theNewAngular[ie].SetValue(ip,theAngular[ie].GetValue(ip));
		881	}
		882	}
		883	delete [] theAngular;
		884	}
		885	theAngular = theNewAngular;
		886	}
797		887
798	//---- Get minimum and maximum energy interpolating	888	//---- Get minimum and maximum energy interpolating
799	theMinEner = angpar1.GetMinEner() + (theEnergy-angpar1.GetEnergy()) * (angpar2.GetMinEner()-angpar1.GetMinEner())/(angpar2.GetEnergy()-angpar1.GetEnergy());	889	theMinEner = angpar1.GetMinEner() + (theEnergy-angpar1.GetEnergy()) * (angpar2.GetMinEner()-angpar1.GetMinEner())/(angpar2.GetEnergy()-angpar1.GetEnergy());