Generate minus num of photons @ G4Cerenkov

cxwx · 2023-07-19T08:46:06.380Z

Hi all,

I set water with RINDEX

                        const int n = 4;
                        double e[n] = {2.034*eV,2.82*eV,4.136*eV, 4.58*eV};
                        double rind[n] = {1.3435, 1.353, 1.3608, 1.3610};
                        fMPT_water->AddProperty("RINDEX",e, rind,n);

and use G4Cerenkov process
Found that it will cause Infinite loop@
G4Cerenkov.cc at special situation

    do
    {
      rand            = G4UniformRand();
      NumberOfPhotons = MeanNumberOfPhotons1 -
                        rand * (MeanNumberOfPhotons1 - MeanNumberOfPhotons2);
      N =
        G4UniformRand() * std::max(MeanNumberOfPhotons1, MeanNumberOfPhotons2);
       
      // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
    } while(N > NumberOfPhotons);

which the MeanNumberOfPhotons1 <0 that is wrong

Geant4 version 11.1.1, 11.1.2, on debian11,
detail value are

dp: 5.77021e-07 
nMin: 1.3435 
nMax: 1.361 
Pmin: 4.00298e-06 
PMax: 4.58e-06 
CaiMin: 1.07494e-06 
CaiMax:1.38692e-06 
ge: 3.11984e-07
NumOfPhotons: -0.00127874

cxwx · 2023-07-20T08:13:20.257Z

there should be something wrong with

  // If n(Pmin) < 1/Beta, and n(Pmax) >= 1/Beta, then we need to find a P such
  // that the value of n(P) == 1/Beta. Interpolation is performed by the
  // GetEnergy() and Value() methods of the G4MaterialPropertiesTable and
  // the Value() method of G4PhysicsVector.
  else
  {
    Pmin = Rindex->GetEnergy(BetaInverse);
    dp   = Pmax - Pmin;

    G4double CAImin = CerenkovAngleIntegrals->Value(Pmin);
    ge              = CAImax - CAImin;

    if(verboseLevel > 1)
    {
      G4cout << "CAImin = " << CAImin << G4endl << "ge = " << ge << G4endl;
    }
  }
}

that it use CerenkovAngleIntegrals->Value(Pmin);
cause (dp - ge * BetaInverse * BetaInverse) <0
which is generate by G4.
~~however I cannot the reproduce same bug on difference system, not sure why~~

cxwx · 2023-07-20T11:27:23.675Z

I’ve tried to reproduce the problem with the below code, I thought it is a bug, did I miss something?
It should be zero but return minus(-)

#include <iostream>
using namespace std;
#include <G4PhysicsTable.hh>
#include <G4PhysicsFreeVector.hh>
#include <G4MaterialPropertiesTable.hh>

void getNumOfCerenPhotons()
{
	double eV = 1.e-6;
	G4PhysicsTable* thePhysicsTable;

	thePhysicsTable = new G4PhysicsTable(1);

	const int n = 4;
	double e[n] = {2.034*eV,2.82*eV,4.136*eV, 4.58*eV};
	//double rind[n] = {1.3435, 1.353, 1.3608, 1.360};
	double rind[n] = {1.3435, 1.353, 1.3608, 1.361};
	G4PhysicsFreeVector* cerenkovIntegral = nullptr;
	G4MaterialPropertiesTable * atable = new G4MaterialPropertiesTable();
	atable->AddProperty("RINDEX", e, rind, n);

	G4MaterialPropertyVector* refractiveIndex = atable->GetProperty(kRINDEX);
	cerenkovIntegral = new G4PhysicsFreeVector();
	double currentRI = (*refractiveIndex)[0];
	if(currentRI > 1.0)
	{
		G4double currentPM  = refractiveIndex->Energy(0);
		G4double currentCAI = 0.0;

		cerenkovIntegral->InsertValues(currentPM, currentCAI);

		G4double prevPM  = currentPM;
		G4double prevCAI = currentCAI;
		G4double prevRI  = currentRI;

		for(std::size_t ii = 1; ii < refractiveIndex->GetVectorLength(); ++ii)
		{
			currentRI  = (*refractiveIndex)[ii];
			currentPM  = refractiveIndex->Energy(ii);
			currentCAI = prevCAI + (currentPM - prevPM) * 0.5 *
				(1.0 / (prevRI * prevRI) +
				 1.0 / (currentRI * currentRI));

			cerenkovIntegral->InsertValues(currentPM, currentCAI);

			prevPM  = currentPM;
			prevCAI = currentCAI;
			prevRI  = currentRI;
		}
	}
	thePhysicsTable->insertAt(0, cerenkovIntegral);
	(*thePhysicsTable)[0]->DumpValues();

	double beta = 0.735288;
	//double beta = 0.734;
	G4double BetaInverse = 1. / beta;

	G4PhysicsVector* CerenkovAngleIntegrals = ((*thePhysicsTable)(0));

	int length = CerenkovAngleIntegrals->GetVectorLength();
	cout<<length<<endl;
	auto Rindex = atable->GetProperty(kRINDEX);
	G4double Pmin = Rindex->Energy(0);
	G4double Pmax = Rindex->GetMaxEnergy();

	G4double nMin = Rindex->GetMinValue();
	G4double nMax = Rindex->GetMaxValue();

	G4double CAImax = (*CerenkovAngleIntegrals)[length - 1];

	G4double dp, ge;
	if(nMax < BetaInverse)
	{
		dp = 0.0;
		ge = 0.0;
	}
	else if(nMin > BetaInverse)
	{
		dp = Pmax - Pmin;
		ge = CAImax;
	}
	else
	{
		Pmin = Rindex->GetEnergy(BetaInverse);
		dp   = Pmax - Pmin;

		G4double CAImin = CerenkovAngleIntegrals->Value(Pmin);
		ge              = CAImax - CAImin;
		cout.precision(7);
		cout
			<<dp<<" "
			<<nMin<<" "
			<<nMax<<" "
			<<Pmin<<" "
			<<Pmax<<" "
			<<CAImin<<" "
			<<CAImax<<" "
			<<ge<<" "
			<<endl;
	}

	//cause <0 numbers

	cout<<(dp - ge * BetaInverse * BetaInverse)<<endl;

}

at my situation, the MeanNumberOfPhotons base on (beta(postPoint)+beta(prePoint)) >0
bt MeanNumberOfPhotons1(prePoint) <0 and MeanNumberOfPhotons1(postPoint) =0,

beta = 0.7350 <0
beta = 0.7349 >0 ***
beta = 0.7348 <0
beta = 0.7347 =0

dsawkey · 2023-07-24T17:08:26.237Z

Hi, are you able to reproduce this with a Geant4 example? Could you please try OpNovice2 with a macro similar to below. I haven’t been able to reproduce this, but of course there are many parameters that may be relevant.

/control/verbose 2
/tracking/verbose 1
/run/verbose 1
/run/numberOfThreads 1
/process/optical/verbose 1
/control/cout/ignoreThreadsExcept 0

/opnovice2/boxProperty RINDEX    0.000002034 1.3435  0.00000282 1.353 0.000004136 1.3608 0.00000458 1.361 
/opnovice2/boxProperty ABSLENGTH 0.000002 10000 0.000005 2000 0.000008 3000

/opnovice2/worldProperty RINDEX    0.000002 1.01    0.000008 1.01
/opnovice2/worldProperty ABSLENGTH 0.000002 1000000 0.000005 2000000 0.000008 3000000

/run/initialize
#
/gun/particle e-
/gun/energy 3 MeV
/gun/position 0 0 0 cm
# will be normalized
/gun/direction .5 .44 0
/opnovice2/gun/optPhotonPolar
/run/initialize
/run/printProgress 1
/run/beamOn 1000

cxwx · 2023-07-26T13:44:24.377Z

It depends on the random seed,
I use G4EmStandardPhysics without option.
I can repoduce it on a server (debian 11), with beamOn(INT_MAX) but I haven’t repoduce on my mac os (ver11.3) yet.
not sure why, even if I use the same random seed.

beta = 0.7350 <0
beta = 0.7349 >0 ***
beta = 0.7348 <0
beta = 0.7347 =0

I’ve use the gdb to pointing the errors
the problem is the geant4 Cerenkov check the >0 with beta = 0.7349>0 (prestep point+ post step point),
so it checked successfully. see

  G4double beta = (pPreStepPoint->GetBeta() + pPostStepPoint->GetBeta()) * 0.5;

  G4double MeanNumberOfPhotons =
    GetAverageNumberOfPhotons(charge, beta, aMaterial, Rindex);

  if(MeanNumberOfPhotons <= 0.0)
  {
    // return unchanged particle and no secondaries
    aParticleChange.SetNumberOfSecondaries(0);
    return pParticleChange;
  }

but after that, it use (beta=0.7350)<0(presteppoint) - (beta=0.7347) = 0(post steppoint)

  G4double beta1 = pPreStepPoint->GetBeta();
  G4double beta2 = pPostStepPoint->GetBeta();

  G4double MeanNumberOfPhotons1 =
    GetAverageNumberOfPhotons(charge, beta1, aMaterial, Rindex);
  G4double MeanNumberOfPhotons2 =
    GetAverageNumberOfPhotons(charge, beta2, aMaterial, Rindex);

so the totally num of photons is <0, which will cause endless loop for my program.

    do
    {
      rand            = G4UniformRand();
      NumberOfPhotons = MeanNumberOfPhotons1 -
                        rand * (MeanNumberOfPhotons1 - MeanNumberOfPhotons2);
      N =
        G4UniformRand() * std::max(MeanNumberOfPhotons1, MeanNumberOfPhotons2);
      // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
    } while(N > NumberOfPhotons);

so if the script return the

beta = 0.7350 <0
beta = 0.7349 >0 ***
beta = 0.7348 <0
beta = 0.7347 =0

It should cause the same problem at some situation.
like average beta ~ 0.7349.. and prebeta ~ 0.7350.. post beta ~ 0.7347..

may be I miss something that cause difference on mac server, but I thought
the beta should not appear like <0 >0 =0

dsawkey · 2023-07-26T15:43:45.121Z

cxwx:

(beta=0.7350)<0(presteppoint) - (beta=0.737) = 0(post steppoint)

Do you have an electric field? How does beta increase?

cxwx · 2023-07-26T15:46:20.487Z

0.7350 > 0.7347
not 0.737,
no field,
sorry for that

dsawkey · 2023-07-26T16:41:59.292Z

Thanks. I’m still struggling to understand what is going on. Could you please put a G4cout (or gdb) in the do-while loop that is going into an infinite loop. What are the values of MeanNumberOfPhotons1 and MeanNumberOfPhotons2, and beta1 and beta2.

cxwx · 2023-07-26T17:18:25.465Z

OK, it takes sometime to recompile g4 with DEBUG mode.
bug at

(gdb) where
#0  CLHEP::MixMaxRng::flat (this=0x555555586800)
    at /home/chenxu/geant4-v11.1.2/source/externals/clhep/include/CLHEP/Random/MixMaxRng.h:66
#1  0x00007ffff4897a82 in G4Cerenkov::PostStepDoIt (this=0x555557141980, aTrack=...,
    aStep=...)
    at /home/chenxu/geant4-v11.1.2/source/processes/electromagnetic/xrays/src/G4Cerenkov.cc:350
#2  0x00007ffff5b18260 in G4SteppingManager::InvokePSDIP (this=this@entry=0x5555559128b0,
    np=np@entry=4)
    ......

#1 endless cycling.

N:0 NumberOfPhotons:-0.00047203 beta1:0.735288 beta2:0.734536 MeanNumberOfPhotons:3.49966e-06 MeanNumberOfPhotons1:-0.00127874 MeanNumberOfPhotons2:0
N:0 NumberOfPhotons:-0.000271159 beta1:0.735288 beta2:0.734536 MeanNumberOfPhotons:3.49966e-06 MeanNumberOfPhotons1:-0.00127874 MeanNumberOfPhotons2:0
N:0 NumberOfPhotons:-0.000853648 beta1:0.735288 beta2:0.734536 MeanNumberOfPhotons:3.49966e-06 MeanNumberOfPhotons1:-0.00127874 MeanNumberOfPhotons2:0
N:0 NumberOfPhotons:-0.00103103 beta1:0.735288 beta2:0.734536 MeanNumberOfPhotons:3.49966e-06 MeanNumberOfPhotons1:-0.00127874 MeanNumberOfPhotons2:0
N:0 NumberOfPhotons:-0.000365598 beta1:0.735288 beta2:0.734536 MeanNumberOfPhotons:3.49966e-06 MeanNumberOfPhotons1:-0.00127874 MeanNumberOfPhotons2:0
N:0 NumberOfPhotons:-0.000302929 beta1:0.735288 beta2:0.734536 MeanNumberOfPhotons:3.49966e-06 MeanNumberOfPhotons1:-0.00127874 MeanNumberOfPhotons2:0
N:0 NumberOfPhotons:-0.000564136 beta1:0.735288 beta2:0.734536 MeanNumberOfPhotons:3.49966e-06 MeanNumberOfPhotons1:-0.00127874 MeanNumberOfPhotons2:0
N:0 NumberOfPhotons:-0.000573361 beta1:0.735288 beta2:0.734536 MeanNumberOfPhotons:3.49966e-06 MeanNumberOfPhotons1:-0.00127874 MeanNumberOfPhotons2:0
N:0 NumberOfPhotons:-0.00101918 beta1:0.735288 beta2:0.734536 MeanNumberOfPhotons:3.49966e-06 MeanNumberOfPhotons1:-0.00127874 MeanNumberOfPhotons2:0
N:0 NumberOfPhotons:-0.000675364 beta1:0.735288 beta2:0.734536 MeanNumberOfPhotons:3.49966e-06 MeanNumberOfPhotons1:-0.00127874 MeanNumberOfPhotons2:0
N:0 NumberOfPhotons:-0.00112915 beta1:0.735288 beta2:0.734536 MeanNumberOfPhotons:3.49966e-06 MeanNumberOfPhotons1:-0.00127874 MeanNumberOfPhotons2:0
N:0 NumberOfPhotons:-0.000483529 beta1:0.735288 beta2:0.734536 MeanNumberOfPhotons:3.49966e-06 MeanNumberOfPhotons1:-0.00127874 MeanNumberOfPhotons2:0
N:0 NumberOfPhotons:-1.02927e-05 beta1:0.735288 beta2:0.734536 MeanNumberOfPhotons:3.49966e-06 MeanNumberOfPhotons1:-0.00127874 MeanNumberOfPhotons2:0
N:0 NumberOfPhotons:-0.000220846 beta1:0.735288 beta2:0.734536 MeanNumberOfPhotons:3.49966e-06 MeanNumberOfPhotons1:-0.00127874 MeanNumberOfPhotons2:0
N:0 NumberOfPhotons:-0.00119925 beta1:0.735288 beta2:0.734536 MeanNumberOfPhotons:3.49966e-06 MeanNumberOfPhotons1:-0.00127874 MeanNumberOfPhotons2:0

while (N=0) always > (MeanNumberOfPhotons1<0)-(MeanNumberOfPhotons2=0)
and MeanNumberOfPhotons is >0

dsawkey · 2023-07-26T18:34:29.270Z

Thanks for doing this! I’d like to try to understand what is going on a bit more. (Writing down my questions so they are documented, not that you have to answer them!)

How does (G4int)G4Poisson(3.5e-6) return 1 or more? There should be protection against returning 0:

260   fNumPhotons          = (G4int) G4Poisson(MeanNumberOfPhotons);
261 
262   if(fNumPhotons <= 0 || !fStackingFlag)

How can MeanNumberOfPhotons1 be less than MeanNumberOfPhotons2, since beta1 > beta2?
Do we need all of MeanNumberOfPhotons, MeanNumberOfPhotons1, MeanNumberOfPhotons2.
Geant4/processes/electromagnetic/xrays/src/G4Cerenkov.cc

cxwx · 2023-07-26T20:37:13.938Z

some info may be needed

step_length: 0.00252619
MeanNumberOfPhotons: 3.4996567673462968312e-06
fNumPhotons: 1
(fNumPhotons<=0): false

however G4Poisson is a random number may return depends on difference system or compiler.

MeanNumberOfPhotons2 = 0
because in the G4Cerenkov::GetAverageNumberOfPhotons

  if(nMax < BetaInverse)
  {
    dp = 0.0;
    ge = 0.0;
  }

but MeanNumberOfPhotons1 gets <0 which I thought unreasonable

dsawkey · 2023-08-09T17:41:54.595Z

Thanks again for reporting this. I was able to reproduce it. As you describe, when it occurs MeanNumberOfPhotons2 = 0 and MeanNumberOfPhotons1 < 0 (but << 1). This leads to the infinite loop you correctly diagnose.

Responding to my own question, G4Poisson(x) has to occasionally return 1 for x>0. G4Poisson returns a value drawn from the Poisson distribution with mean x.

I will patch the code shortly.