Hi,
I have a 2D histogram to plot how energy is deposited in a material. I want to find the moliere radius of the material, the radius within which 90% of the energy is deposited.
I have tried fitting slices in Y and looking at the parameters that gives, however none of these seems to be sensible values.
I have attached the histogram for tungsten which has a moliere radius of 9.3mm, and ideally I’d like the to be able to find this just from the histogram.
On the X axis is the the depth of the material, i.e. the direction my beam is travelling.
On the Y axis is the radius or distance from the beam line.
Each bin contains the total energy deposited there.
Thanks,
Alex.
Dear @alex.smith3 ,
Perhaps @moneta or @couet can help you.
Cheers,
Vincenzo
Dear @alex.smith3,
Below there is a script to reproduce your problem.
The idea is to project your 2D histogram along the y axis and then start to sum the bin content until you reach the 90%. Obviously the finer is the binning on the y axis the more precise the value that you found.
{
//create a 2D histogram
TH2D *dep_energy = new TH2D("h2_dep_energy","Z vs radius :Z: Radius",320,0.,160.,800,0.,200.);
TRandom3 *rnd= new TRandom3(1);
double z, r;
int i=0;
int bin_90_energy=0;
double sigma =40.;
//fill the histogram uniformly along x axis and with the absolute value of a gaussian on y axis
for(i=0; i<1e6; ++i)
dep_energy->Fill(rnd->Uniform(200.) ,abs(rnd->Gaus(0.,sigma)));
//project the histogram on y axis
TH1D *dep_energy_py= dep_energy->ProjectionY();
double integral=0.;
double total_event=dep_energy_py->GetEntries();
i=1;
//looping over the bin of the new histogram to found at which bin the count is over the 90% of the total count
while(integral<0.90 && i <= dep_energy_py->GetNbinsX())
{
integral+= dep_energy_py->GetBinContent(i)/total_event;
++i;
}
//reduce i by 1 to account for the i++ before exiting the while loop
bin_90_energy = i-1;
// print the center of the bin I found
cout<<"Integral with 90% of content is found at radius = "<<dep_energy_py->GetBinCenter(bin_90_energy)<<endl;
return;
}
Stefano
PS.
You can also double check your result trying to evaluate the radius that contains 90% and 95% of your energy distribution and your second radius should be twice the first one
Dear Stefano,
Thanks for your help!
Your code didn’t seem to get close to 90% energy deposition, stopping at 70% when it ran out of bins, and I’m not sure why. I modified the code which gets a somewhat sensible value for moliere radius - 10.6mm instead of 9.3mm, however this completely fails when set to 95% - getting 17.0mm which is not close to double.
I’ve attached my modified code below for anyone else who may want something similar.
{
gROOT->Reset();
gROOT->SetStyle("Plain");
// Open file filled by Geant4 simulation
TFile f("TungstenOnly.root");
// Get ntuple
TNtuple* ntuple = (TNtuple*)f.Get("B4_SpacialData");
TCanvas* c1 = new TCanvas("c1", "", 20, 20, 800, 800);
gPad->SetLogz(1);
ntuple->Draw("(X^2+Y^2)^(1/2):Z>>hXYE(200,-85,85,200,0,115)","Edep","colz");
TH1D *dep_energy_py = hXYE->ProjectionY();
double integral=0.;
//double total_event=dep_energy_py->GetEntries();
double total_event=0.;
i=1;
while(i <= dep_energy_py->GetNbinsX())
{
total_event+= dep_energy_py->GetBinContent(i);
++i;
}
cout<<"Total events = " <<total_event<<endl;
i=1;
//looping over the bin of the new histogram to found at which bin the count is over the 90% of the total count
while(integral<0.90 && i <= dep_energy_py->GetNbinsX())
{
integral+= dep_energy_py->GetBinContent(i)/total_event;
cout << dep_energy_py->GetBinContent(i)/total_event<<endl;
++i;
}
cout <<"Integral reached = "<< integral <<endl;
//reduce i by 1 to account for the i++ before exiting the while loop
bin_90_energy = i-1;
// print the center of the bin I found
cout<<"Integral with 90% of content is found at radius = "<<dep_energy_py->GetBinCenter(bin_90_energy)<<endl;
return;
}
I think the issue may be with GetEntries compared to GetBinContent, but I’m not sure.
Any more advice would be much appreciated
Alex.
The content of your 2D histogram bins are not integer but double.
Instead of
dep_energy_py->GetEntries();
you should use
dep_energy_py->Integral();
To improve the results try to use 400 or 800 as bin on the y-axis.
Stefano
This appears to lead to the same issue as my code has, producing exactly the same values. I think the issue is potentially in the data generation or collection.
Thanks for your help,
Alex.
