Project

General

Profile

Bug #15071

Nutools does not recognize nuclear PDG codes

Added by Steven Gardiner over 2 years ago. Updated over 2 years ago.

Status:
Resolved
Priority:
Normal
Start date:
01/10/2017
Due date:
% Done:

0%

Estimated time:
Duration:

Description

Currently, when nutools encounters PDG codes for nuclei, it often fails to recognize them. Running the LArSoft FHiCL script prodmarley.fcl (in the larsim repository) will always give an error message (created by the ConvertPrimaryToGeant4 class in nutools) similar to the following example:

"%MSG-w ConvertPrimaryToGeant4: RootOutput:out1@EndJob 10-Jan-2017 15:37:15 CST PostEndRun
The following unknown PDG codes were present in the simb::MCTruth input.
They were not processed by Geant4.
Unknown PDG code = 1000180390, seen 2 times.
Unknown PDG code = 1000190390, seen 11 times.
Unknown PDG code = 1000190400, seen 87 times.
%MSG
Art has completed and will exit with status 0."

History

#1 Updated by Robert Hatcher over 2 years ago

... when nutools encounters PDG codes for nuclei, it often fails to recognize them ... It's not so much "often" as "under some circumstances".

What's happening here is that the part of NuTools that converts a MCTruth object to fill Geant4 primary objects is relying on ROOT's class TDatabasePDG. By default that class when invoked and is uninitialized, it reads the file $ROOTSYS/etc/pdg_table.txt for its information. That file does not have a complete list of PDG codes for most nuclei. But when GENIE has been invoked via GENIEHelper it forces TDatabasePDG to read from $GENIE/data/evgen/catalogues/pdg/genie_pdg_table.txt before it can be accidentally invoked otherwise (and thus read the default ROOT version). The GENIE file does indeed have the desired PDG codes:

$ egrep '1000180390|000190390|1000190400' $GENIE/data/evgen/catalogues/pdg/genie_pdg_table.txt
  999 Ar39          1000180390 0 100 Ion        54 36.2858295   0.00000e+00 -100  -1 -100 -1 -1 0
  999 K39           1000190390 0 100 Ion        57 36.2847535   0.00000e+00 -100  -1 -100 -1 -1 0
  999 K40           1000190400 0 100 Ion        57 37.2165193   0.00000e+00 -100  -1 -100 -1 -1 0

So, it's a matter of getting the right file read, once and only once, in all cases. Part of the problem is that the mechanism for directing TDatabasePDG's "favorite" data file is ... not very flexible. And another confounding issue is that while the class is meant to be a "singleton" isn't not very good at ensuring that it only reads one file and one can accidentally have it read two different files, or the same file twice (which screws up the branching factor info, as well as generating error messages if I remember correctly). Which is why GENIE wraps it in a class PDGLibrary:

One could use (in the module constructor):

  #include "PDG/PDGLibrary.h" 
  [...]
  PDGLibrary* pdglib = PDGLibrary::Instance(); // get message out of the way

Though that has two (smallish?) drawbacks:

1. It requires linkage to GENIE libraries in places otherwise not needed.
2. Invoking this causes the GENIE banner to be output

Ideally the "best" way to deal with this would to make improvements to TDatabasePDG and get them accepted/adopted by the ROOT team. And/or one could approach ROOT for adopting the GENIE file (with it's nuclei codes, but also some GENIE specific codes).

Alas, either of those is a "long game" approach, so no one ever seems to do it. NuTools could also do an approach that wrapped TDatabasePDG (and pointed it at the GENIE file) but would have to ensure that it played "nice" with the GENIE version, lest it introduce problems when running the GENIE event generator.

All that said, one should take note that this:

"%MSG-w ConvertPrimaryToGeant4: RootOutput:out1@EndJob 10-Jan-2017 15:37:15 CST PostEndRun
The following unknown PDG codes were present in the simb::MCTruth input.
They were not processed by Geant4.
Unknown PDG code = 1000180390, seen 2 times.
Unknown PDG code = 1000190390, seen 11 times.
Unknown PDG code = 1000190400, seen 87 times.
%MSG
Art has completed and will exit with status 0." 

is "PostEndRun". The job didn't fail, what it did do was not hand those nuclei to Geant4 for transport. Since I'm going to guess that those are remnant recoil nuclei with little kinetic energy, I'm not sure how much this affects anything.

#2 Updated by Steven Gardiner over 2 years ago

Very enlightening, thanks! My proposed fix doesn't require anything other than Geant4 itself. It queries the G4ParticleTable first, then asks the G4IonTable to get the definition (which is created on the fly by the table if it doesn't exist yet) if G4ParticleTable didn't recognize the PDG code:

// Get Geant4's definition of the particle.
G4ParticleDefinition* particleDefinition;

if(pdgCode==0){
  particleDefinition = fParticleTable->FindParticle("opticalphoton");
}
else
  particleDefinition = fParticleTable->FindParticle(pdgCode);

if ( pdgCode > 1000000000) { // If the particle is a nucleus
  LOG_DEBUG("ConvertPrimaryToGeant4") << ": %%% Nuclear PDG code = " << pdgCode
                                      << " (x,y,z,t)=(" << x
                                      << "," << y
                                      << "," << z
                                      << "," << t << ")" 
                                      << " P=" << momentum.P()
                                      << ", E=" << momentum.E();
  // If the particle table doesn't have a definition yet, ask the ion
  // table for one. This will create a new ion definition as needed.
  if (!particleDefinition) {
    int Z = (pdgCode % 10000000) / 10000; // atomic number
    int A = (pdgCode % 10000) / 10; // mass number
    particleDefinition = fParticleTable->GetIonTable()->GetIon(Z, A, 0.);
  }
}

Would this kind of thing be an acceptable workaround until one of the "long game" approaches emerges?

And you're right, this is typically only a minor issue.

#3 Updated by Lynn Garren over 2 years ago

  • Status changed from New to Resolved


Also available in: Atom PDF