Hello everyone,
This is going to be a long post, so in short:
- There is a normalization bug in TestEm5 example which makes simulation results far from experimental data (histogram 12, “(transmit, charged) : space angle dN/dOmega”, and most probably histograms 22, 32, and 42).
- There are bugs in ROOT macro files in TestEm5 which result in SegFault because names of histograms saved by G4 have “h” prefix while macros don’t. From history log seems that latter “h” has been added in front of histogram names (numbers) in G4, but this has not been changed in macro files (26-06-14 mma).
- Hanson data ( A. O. Hanson, L. H. Lanzl, E. M. Lyman, and M. B. Scott Phys. Rev. 84 , 634, 1951) included with TestEm5 are wrongly digitized from the original paper, especially in the low \theta range.
- I’ve found a couple of physics validation comparisons with these wrong data, surprisingly matching it well. Currently, it seems that some models match the wrong data, and some overlap with correct data. I would appreciate G4 developers comment on that, as it make simulation of MSC process questionable in some physics models.
- I will be happy to correct these bugs, so I would appreciate some advise how to do that (report it through Bugzilla? submit a commit through GitHub? this post is enough?)
Ok, so starting from the beginning:
./TestEm5 hanson.mac
cd ./hanson
root hanson.C -l
gives the following error:
root hanson.C -l
root [0]
Processing hanson.C…
Error in TFile::TFile: file ./19um.opt0.root does not exist
as the “hanson.mac” does not specify “.root” extension for the output file. It requires to change analysis output file name in “hanson.mac” (or renaming 19um.X output files in ./hanson dir) from
/analysis/setFileName old-hanson/19um.local
changed to
/analysis/setFileName old-hanson/19um.local.root
After adding “.root” to output file name in “hanson.mac” we get another error from hanson.C:
root hanson.C -l
root [0]
Processing hanson.C…
#0 0x00007fd7f610d687 in __GI___waitpid (pid=13228, stat_loc=stat_loc
entry=0x7ffdd62c0ac8, options=options
entry=0) at …/sysdeps/unix/sysv/linux/waitpid.c:30
#1 0x00007fd7f6078067 in do_system (line=) at …/sysdeps/posix/system.c:149
#2 0x00007fd7f6d1ce63 in TUnixSystem::StackTrace() () from /home/daq/code/root/lib/libCore.so.6.18
#3 0x00007fd7f19870f5 in cling::MultiplexInterpreterCallbacks::PrintStackTrace() () from /home/daq/code/root/lib/libCling.so
#4 0x00007fd7f1986afb in cling_runtime_internal_throwIfInvalidPointer () from /home/daq/code/root/lib/libCling.so
#5 0x00007fd7f74f35d3 in ?? ()
#6 0x00007fd7f74f3430 in ?? ()
#7 0x0d8c27c4e1aaf000 in ?? ()
#8 0x000055be0fb6fc40 in ?? ()
#9 0x000055be106759e8 in ?? ()
#10 0x000055be0fb6fc40 in ?? ()
#11 0x00007fd7f7503358 in ?? ()
#12 0x0000000000000000 in ?? ()
terminate called after throwing an instance of ‘cling::InvalidDerefException’
what(): Trying to dereference null pointer or trying to call routine taking non-null arguments
That’s not very helpful, but I’ve found that it is just a wrong histogram name. G4 analysis output ROOT files have histograms called “h12” while “hanson.C” calls just for just “12”:
// Draw histograms fill by Geant4 TestEm11 simulation
TFile f1(“./19um.opt0.root”);
TH1D* h1 = (TH1D*) f1.Get(“12”);
So it requires to change it into:
TH1D* h1 = (TH1D*) f1.Get(“h12”);
and the following two histograms
TH1D* h1 = (TH1D*) f1.Get(“h12”);
and
TH1D* h3 = (TH1D*) f3.Get(“h12”);
After running ./TestEm5 hanson.mac for three times (to generate opt0, opt3, and local models), we get the following result:
There is an obvious problem with normalization. First I tried a couple of standard tricks like:
Double_t scale = h1->GetXaxis()->GetBinWidth(1)/(h1->Integral());
h1->Scale(scale);
and other combinations of normalizing \theta histograms with integral, bin width etc. but nothing really worked.
Then, I looked into TestEm5 source code, especially “TrackingAction.cc” where the normalization (weight) is calculated for histograms 12, 22, 32 and 42:
> //space angle distribution at exit : dN/dOmega
> //
> G4ThreeVector direction = aTrack->GetMomentumDirection();
> id = 0;
> if (transmit && charged) id = 12;
> else if (transmit && neutral) id = 22;
> else if (reflect && charged) id = 32;
> else if (reflect && neutral) id = 42;
>
> if (id>0) {
> G4double theta = std::acos(direction.x());
> if (theta > 0.0) {
> G4double dteta = analysisManager->GetH1Width(id);
> G4double unit = analysisManager->GetH1Unit(id);
> G4double weight = (unit*unit)/(twopi*std::sin(theta)*dteta);
> /*
> G4cout << "theta, dteta, unit, weight: "
> << theta << " "
> << dteta << " "
> << unit << " "
> << weight << G4endl;
> */
> analysisManager->FillH1(id,theta,weight);
> }
> }
“unit” is (radian/degree) = 0.01745(…), while “dteta” and “theta” are in (degree). Simple dimensional analysis shows that this weight has wrong units (this formula will give (radian/degree^3) unit which is nonsensical, should be 1/degree^2).
So the weight should be calculated as follows:
G4double weight = (unit)/(twopi*std::sin(theta)*dteta);
which will give the correct unit of (1/degree^2), as in the original Hanson et al. article. After that correction, recompiling, and running ./TestEm5 hanson.mac three times more we get:
(notice linear scale)
Pretty good, but I was a bit puzzled with two things: low angle scattering seems to deviate from ‘local’ physics list, and these data which are shipped with TestEm5 have data points at different \theta angle than the original article (notice that data points at 15, 20, 25 degrees are shifted ~2 degrees to the right in the original Hanson et al. data). There are also much more data points in the Hanson’s article:
This prompted me to digitize Hanson et al. data by myself. Then, I noticed that not only 20, 25 degrees point were approximated, but low angle points are quite off:
“Hanson old data” are data points from the ASCII file distributed together with TestEm5, “Hanson new data” are manually digitized by me. What is a bit worrying is that models actually fit better the old wrongly digitized data…
I found some old G4 physics validation plots from DESY, and it seems that these wrong data have been used previously for physics validation. You can easily spot that by looking at 20 and 25 degree points:
Of course, when you plot everything on a log scale, the difference is negligeble:
Interestingly, Penelope physics list deviates the most from the correct data.
Ok, so this is the end of the story. It seems that all models overestimate slightly number of not deflected electrons (\theta<0.25 degree). Is it something what people are aware off? Also, this physics validation:
http://geant4.web.cern.ch/sites/geant4.web.cern.ch/files/geant4/collaboration/working_groups/electromagnetic/verification/angle2.gif
compares against wrong data (with overestimated number of low scattering events) and shows good agreement (!). Are there any more recent validation efforts which verify multiple scattering process? Who shall I contact to repair bugs in TestEm5 example and update Hanson data?
Best,
Weronika
