Multifunctional detector not agreeing with mesh scorer

I think I am having trouble reading results from my multifunctional detector properly. I set up a general particle source which varies in intensity which I can confirm from viewing it in opengl

However when I try to score dose from this using a multifunctional detector all the sources appear to deposit dose a nearly uniform intensity, giving a line profile down them like

The peaks appear to get more noisy as the flux should decrease which makes me feel like they are being normalised according to the number of photons travelling through the region, which is pretty non-intuitive to me.

When I instead run a scoring mesh using the macro commands I get a line profile like

which follows the trend I expect. However I only get it when I use the sum_w() method
MeshScoreMap scMap = fScoringMesh->GetScoreMap();
MeshScoreMap::const_iterator msMapItr = scMap.find(psName);
std::map<G4int, G4StatDouble*> * score = msMapItr->second->GetMap();
std::map<G4int, G4StatDouble*>::iterator value = score->find(idx);
double totED = value->second->sum_w() / unitValue;
Which I am unable to use with my multifunctional detector. My understanding is that the sum_w method seems to multiply the entry by some track weight but I have no idea how to access this.

Is anyone able to help out with why a multifunctional detector is not returning the correct scoring values?
At the moment I am accessing the multifunctional detector values through
if (IsMaster())
Run* reRun = (Run*) run;
G4SDManager* SDman = G4SDManager::GetSDMpointer();
G4MultiFunctionalDetector* mfd = (G4MultiFunctionalDetector*)(SDman->FindSensitiveDetector(“phantomSD”));
G4String Unit = mfd->GetPrimitive(0)->GetUnit();
G4double UnitValue = mfd->GetPrimitive(0)->GetUnitValue();

        G4THitsMap<G4double>* Dose = reRun->GetHitsMap("phantomSD/DoseDeposit");
        for (int Index = 0; Index < 100 * 100 * 3; Index++)
            G4double* totED = (*Dose)[Index];
            if (totED)
                std::cout << "Dose at index " << Index << " is " << *totED << std::endl;


Eternal gratitude to anyone who can help out with this

how do you set up the gps? does this involve track weights / biasing?

Hey Weller
I use a bias with the macro
/gps/hist/type biasx
and following that set a histogram.

These are applied as macro commands after running /run/initialize and before running beamOn. Do I need to somehow convey this information to the sensitive detector?

The biasing is automatically applied to the tracks as a weight (see G4Track::GetWeight()), and it’ll be propagated automatically to secondaries.

If you’re using biasing, then you should include the weighting when you fill histograms. I think the scorer is taking care of that for you automatically (at least, because you’ve requested the sum_w() (sum of weights) method. In your SD, you need to make sure that you either store the weight in your N-tuple for offline access, or included when filling your histograms.

I think in this case the weighting should not be included, as biasx is used to create different intensities in the four different regions. as shown, accounting for the weighting (command based scoring with the mesh automatically includes the weights) compensates for the lower number of primaries and all four peaks are correctly but unintentionally scaled to the same height.

can you try to
use 4 separate gps sources with different intensities instead? this should leave the track weight at 1…

Hello Weller and mkelsey,

I made a slight mistake in the first post sorry, in order to get dose accumulating in the scoring mesh macro I actually need to use the sum_wx method which returns a dose reading, however it ends up giving the same issue as my sensitive detector that the values appear to be normalised.
Is this all expected behaviour that intensity biases are incompatible with sensitive detectors/ scoring meshes, or is there a way to make them work? The end result I want is to output a 3D array of the dose measured at each position

Unfortunately dividing the GPS into 4 is just a simplification for troubleshooting, my intent is to use the biasing to create gaussian-like shapes in x and y to match beam profiles I measure during experiments. I would need to use potentially millions of particle sources to achieve this at the same resolution which was not achievable when I first tried it.

If I am not mistaken, you could use the sigma_r (or sigma_x/y) parameter instead of bias:

Although the profiles can sometimes appear gaussian I don’t think sigma_x/y will work for them since they are really arbitrary depending on the experiment. Often they can be off centre, offset, asymmetry and decay faster than an actual gaussian.

Would this same behaviour also be seen with angular and energy biasing?
Is there be any way I could correct it during a stepping action or something similar?
Thanks for the help,

I don’t know, but you can easily check the track weight in your user stepping action:

G4cout << "TrackWeight: " << aStep->GetPreStepPoint()->GetWeight() << G4endl;

Would it be feasible to manipulate the gps at code level instead of macro configuration?

void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) {
    G4double x0 = ... // code to calculate...

Okay, after some reflection I think I have found the source of my confusion. I had thought the purpose of bias would be to set a non constant intensity across a source, now I’m thinking it is actually used to vary the sampling rate across a source?
Does this mean it is impossible to create a source with arbitrary intensity?

Would it be feasible to manipulate the gps at code level instead of macro configuration?

The idea being that this would potentially define the sources faster than macro commands?
I think I was going to try this before finding the bias commands but found only had instructions in my manual for the macro versions.
Do you know of any other documentation?

And would I be able to change the detector setup between runs?

no. It may just not be possible with macro commands. Angular distributions however are mentioned in the manual:

The user can define histograms for several reasons: angular distributions in either θ
or ϕ; energy distributions; energy per nucleon distributions; or biasing of x, y, z, θ, ϕ, or energy

There is also nothing like that in the exgps example.

The idea being that this enables you to define arbitrary intensity distributions :slight_smile:

You can make your simulation read the source intensity distribution from a file or some command line argument, and then pass the content to your G4VUserPrimaryGeneratorAction class via G4VUserActionInitialization.

Thank you Weller.
My current plan moving forward is to try defining all the sources individually at the code level and hope this evaluates in a reasonable time unlike when I tried doing it using macro commands.
I’ve figured out how to do all the same definitions I otherwise would have in macro, my main confusion at the moment is where to define the sources. My understanding is that generateprimaries is called at the beginning of every event meaning I would be consistently repeating the definitions which could cause issues.
I have veered pretty far away from the original topic of this post so if necessary I will start a new one.
Thanks again for your time,

edit; on second thoughts from examples involving G4ParticleGun it may be that I need to do the random number generation for each event myself and position the source appropriately then

that’s what I had in mind here

You can also use the particle gun for that. Might be even better suited as you don’t need all the additional/predefined features of gps…

(I would not worry about the performance differences, as long as you don’t make severe mistakes such as reading in a distribution from file for every single event.)

1 Like