Going On with the Check of Geant4

Nov 14 2013

Author: Andrey Karpov

Summary of the previous article
New suspicious code fragments
Conclusion

This is the correct article about the results of checking the Geant4 project, which I have written after the previous incorrect one. Let me remind you the whole story. I have recently checked an old version of the Geant4 library and reported the results in the article "Copy-Paste and Muons". Why old version? Nobody is perfect, and we finally made a mistake ourselves. To find out which exactly, see the previous article. This time I offer you a brief report about checking Geant4 of the version 10.0-beta.

Summary of the previous article

In the article "Copy-Paste and muons" I was reasoning about how useful the static analysis methodology is and about the PVS-Studio analyzer's diagnostic abilities. I checked an ancient version of the Geant4 project (version 4.9.4), found a number of supposedly incorrect code fragments and described them in the article.

Geant4 (for GEometry ANd Tracking) is a platform for "the simulation of the passage of particles through matter," using Monte Carlo methods. It is the successor of the GEANT series of software toolkits developed by CERN, and the first to use object oriented programming (in C++). Its development, maintenance and user support are taken care by the international Geant4 Collaboration. Application areas include high energy physics and nuclear experiments, medical, accelerator and space physics studies. The software is used by a number of research projects around the world.

The project website: http://geant4.org.

In the previous article, I described at least 16 suspicious code fragments. The recent check of the new version has revealed only 10 of those. The rest are either fixed or discarded together with the pieces of code. I won't discuss those defects once again in this article; if you want to see them, please refer to the previous article (there are comments below each sample telling if a particular bug is fixed in the new version or not.)

Sorry for such a strange format of my investigation, but I hope that it won't in any way prevent the developers from fixing some defects in the Geant4 project and that the PVS-Studio will attract their attention.

If I'm correct, the previous library version dates back to 2011. Much has changed since then, and it's no wonder that some new strange code fragments were found. Let's see if there is anything new or anything I missed during the previous check.

New suspicious code fragments

A complete list of all the suspicious fragments that caught my eye is saved in the file geant4_new.txt. But please don't rely solely on this list; the developers should check the project themselves and study all the warnings. We can grant a free registration key to the Geant4 developers for some time so that they could check their code: see the feedback page.

Identical functions

G4double G4CsvAnalysisManager::GetH2Xmin(G4int /*id*/) const
{
  ExceptionForHistograms("GetH2Xmin");
  return 0;
}

G4double G4CsvAnalysisManager::GetH2Xmax(G4int /*id*/) const
{
  ExceptionForHistograms("GetH2Xmin");
  return 0;
}

PVS-Studio's diagnostic message: V524 It is odd that the body of 'GetH2Xmax' function is fully equivalent to the body of 'GetH2Xmin' function. _G4analysis-archive g4csvanalysismanager.cc 933

The GetH2Xmax() function should probably call the ExceptionForHistograms() function with a different parameter:

ExceptionForHistograms("GetH2Xmax");

This bug doesn't look serious. As far as I understand, this is a construct to handle an exception. However, I decided to mention this Copy-Paste bug anyway.

Zero energy

The function CalculateTotalEnergy () sums up values in the 'Etot' variable which suddenly gets zeroed. That's a very strange thing.

G4double G4RKFieldIntegrator::CalculateTotalEnergy(const
  G4KineticTrackVector& Barions)
{
  G4double Etot = 0;
  ....
  for(G4int c2 = c1 + 1; c2 < nBarion; c2++)
  {  
    ....
    //  Esk2
    Etot += t1*std::pow(Alpha/pi, 3/2)*
            std::exp(-Alpha*r12*r12);

    // Eyuk
    Etot += ....;

    // Ecoul
    Etot += 1.44*p1->GetDefinition()->GetPDGCharge()*
            p2->GetDefinition()->GetPDGCharge()/r12*
            Erf(std::sqrt(Alpha)*r12);

    // Epaul
    Etot = 0;
    ....
  }
  ....
}

PVS-Studio's diagnostic message: V519 The 'Etot' variable is assigned values twice successively. Perhaps this is a mistake. Check lines: 80, 83. _G4processes-archive g4rkfieldintegrator.cc 83

Strange logic

G4double G4EmBiasingManager::ApplySecondaryBiasing(....)
{
  ....
  if(0 == nsplit) { 
    if(safety > fSafetyMin) ....
  } if(1 == nsplit) { 
    weight = ApplyRussianRoulette(vd, index);
  } else {
    G4double tmpEnergy = pPartChange->GetProposedKineticEnergy();
    G4ThreeVector tmpMomDir = ....
    weight = ApplySplitting(vd, track, currentModel, index, tcut);
    pPartChange->SetProposedKineticEnergy(tmpEnergy);
    pPartChange->ProposeMomentumDirection(tmpMomDir);
  }
  ....
}

PVS-Studio's diagnostic message: V646 Consider inspecting the application's logic. It's possible that 'else' keyword is missing. _G4processes-archive g4embiasingmanager.cc 299

The code formatting suggests that the programmer was using the "else if" construct. But I can't see any here. If we format the code properly, we'll get the following:

if(0 == nsplit) { 
  if(safety > fSafetyMin) ....
}

if(1 == nsplit) { 
  weight = ApplyRussianRoulette(vd, index);
} else {
  G4double tmpEnergy = pPartChange->GetProposedKineticEnergy();
  G4ThreeVector tmpMomDir = ....
  weight = ApplySplitting(vd, track, currentModel, index, tcut);
  pPartChange->SetProposedKineticEnergy(tmpEnergy);
  pPartChange->ProposeMomentumDirection(tmpMomDir);
}

Notice that the block referring to the 'else' branch is executed every time whenever the "1 != nsplit" condition is true. It suspect that the programmer wanted the program logic to be somewhat different.

A similar issue can be found in the following fragment:

V646 Consider inspecting the application's logic. It's possible that 'else' keyword is missing. _G4processes-archive g4embiasingmanager.cc 347

Incomplete code?

void G4MolecularDecayTable::AddExcitedState(const G4String& label)
{
  channelsMap::iterator channelsIter =
    fDecayChannelsMap.find(label);
  if(channelsIter != fDecayChannelsMap.end())
  {
    G4String errMsg = "Excited state" + label +
                      " already registered in the decay table.";
    G4Exception("G4MolecularDecayTable::AddExcitedState",
                "G4MolecularDecayTable003",
                FatalErrorInArgument, errMsg);
     return;
  }
  fDecayChannelsMap[label] ;
}

PVS-Studio's diagnostic message: V607 Ownerless expression 'fDecayChannelsMap[label]'. _G4processes-archive g4moleculardecaytable.cc 140

The end of the function is very strange:

fDecayChannelsMap[label] ;

What's it? Something is missing? What did the programmer intend to do with an array cell?

Kinematics

The next sample is rather lengthy. I shortened it as much as I could, but unfortunately it is still pretty large. Read it and note the values the 'id' variable takes.

void G4QMDCollision::CalKinematicsOfBinaryCollisions(
  G4double dt)
{
  ....
  G4int id = 0;
  ....
  if ( secs )
  {
    ....
    id++;
    ....
  }
  if ( std::abs ( eini - efin ) < fepse*10 ) 
    ....
  else
  {  
    ....             
    for ( G4int i0i = 0 ; i0i < id-1 ; i0i++ )
    {
      theSystem->DeleteParticipant( i0i+n0 );
    }
    ....
  }
  ....
}

PVS-Studio's diagnostic message: V621 Consider inspecting the 'for' operator. It's possible that the loop will be executed incorrectly or won't be executed at all. _G4processes-archive g4qmdcollision.cc 228

If the "if ( secs )" condition is false, the variable 'id' will remain equal to zero. In this case we might get the following loop:

for ( G4int i0i = 0 ; i0i < -1 ; i0i++ )

And that will be a very strange loop indeed. I guess something is wrong with the CalKinematicsOfBinaryCollisions() function's logic.

Miscellaneous

There are a few more warnings I didn't describe in the previous article. Neither will I in this one. Just one sample:

class G4HadronicException : public std::exception
{
  ....
};

inline G4double G4GeneralPhaseSpaceDecay::Pmx(
  G4double e, G4double p1, G4double p2)
{
   if (e-p1-p2 < 0 )
   {  
    G4HadronicException(__FILE__, __LINE__,
      "G4GeneralPhaseSpaceDecay::Pmx "
      "energy in cms > mass1+mass2");
   }
  ....
}

PVS-Studio's diagnostic message. V596 The object was created but it is not being used. The 'throw' keyword could be missing: throw G4HadronicException(FOO); _G4processes-archive g4generalphasespacedecay.hh 116

There are some mistakes with missing 'throw' operators. They result in an object of the 'G4HadronicException' type being created and at once destroyed, the program going on to work with incorrect data.

For other examples of such typos, see the file geant4_new.txt. There you will also find some warnings related to microoptimizations.