debugging energy conservation

.gitignore

@@ -2,3 +2,4 @@ sim/build/
 *__pycache__*
 */log/*
 *root
+*.ipynb*

sim/exampleB4b.cc

@@ -47,6 +47,7 @@
 #include "G4UIExecutive.hh"
 
 #include "CaloTree.h"
+#include "B4bPhysicsList.hh"
 
 #
 using namespace std;
@@ -128,11 +129,14 @@ int main(int argc, char **argv)
 
   //   optical physics from examples/extended/optical/OpNovice2
   //  auto physicsList = new FTFP_BERT;
-  auto physicsList = new QGSP_BERT;
+  // auto physicsList = new QGSP_BERT;
+  auto physicsList = new CustomPhysicsList();
 
   physicsList->ReplacePhysics(new G4EmStandardPhysics_option4());
   G4OpticalPhysics *opticalPhysics = new G4OpticalPhysics();
   physicsList->RegisterPhysics(opticalPhysics);
+  physicsList->RemovePhysics("photonNuclear");
+  physicsList->RemovePhysics("positronNuclear");
 
   // G4Cerenkov* theCerenkovProcess=new G4Cerenkov("Cerenkov");
   // theCerenkovProcess->SetTrackSecondariesFirst(true);
@@ -186,7 +190,6 @@ int main(int argc, char **argv)
 
     // string evtmax="100";
     command = "/run/beamOn " + histo->getParamS("numberOfEvents");
-    ;
     cout << "command: " << command << endl;
     UImanager->ApplyCommand(command);
   }

sim/include/B4bPhysicsList.hh

@@ -0,0 +1,33 @@
+#include "QGSP_BERT.hh"
+#include "G4Gamma.hh"
+#include "G4ProcessManager.hh"
+
+class CustomPhysicsList : public QGSP_BERT
+{
+public:
+    CustomPhysicsList() : QGSP_BERT() {}
+    ~CustomPhysicsList() {}
+
+    void ConstructProcess() override
+    {
+        QGSP_BERT::ConstructProcess(); // Call base class to set up standard processes
+
+        // Now, remove the photonNuclear process
+        G4ParticleDefinition *gamma = G4Gamma::GammaDefinition();
+        G4ProcessManager *pmanager = gamma->GetProcessManager();
+        if (pmanager)
+        {
+            G4int nProcesses = pmanager->GetProcessListLength();
+            for (G4int i = 0; i < nProcesses; i++)
+            {
+                G4VProcess *process = (*pmanager->GetProcessList())[i];
+                if (process->GetProcessName() == "photonNuclear")
+                {
+                    pmanager->RemoveProcess(process);
+                    G4cout << "Removed process: " << process->GetProcessName() << G4endl;
+                    break;
+                }
+            }
+        }
+    }
+};

sim/include/CaloTree.h

@@ -39,6 +39,7 @@ class CaloTree
   //  called fro SteppingAction...
   void accumulateHits(CaloHit aHit);
   void accumulateEnergy(double eleak, int type);
+  void accumulateDeposits(double edep, int step);
   void saveBeamXYZE(string, int, float, float, float, float);
 
   // for histogrming...
@@ -157,6 +158,8 @@ class CaloTree
   double m_eCentruth;
   double m_eScintruth;
 
+  std::vector<double> m_eDEPs;
+
   // scintillation hit variables
   int m_nhits3dSS;
   vector<int> m_id3dSS; // channel ID   xxxyyyzzz

sim/paramBatch03_single.mac

@@ -45,5 +45,5 @@
 /run/initialize
 #
 # /tracking/verbose 1
-/tracking/verbose 0
+/tracking/verbose 1
 #

sim/src/B4bSteppingAction.cc

@@ -39,6 +39,7 @@
 
 #include "G4PhysicalConstants.hh"
 #include "G4SystemOfUnits.hh"
+#include "G4Positron.hh"
 
 #include "CaloID.h"  // including CaloID, CaloHit, CaloTree
 #include "CaloHit.h" // including CaloID, CaloHit, CaloTree
@@ -91,7 +92,6 @@ void B4bSteppingAction::UserSteppingAction(const G4Step *step)
 
   if (particleDef == opticalphoton)
   {
-    // cout<<"this is an optical photon.  Need to stop it"<<endl;
     track->SetTrackStatus(fStopAndKill);
     return;
   }
@@ -123,10 +123,10 @@ void B4bSteppingAction::UserSteppingAction(const G4Step *step)
   // if(tkstatus==fStopAndKill && absPdgCode>100 && track->GetTrackID()==1) {
   // if(tkstatus==fStopAndKill && absPdgCode>100) {
   // save secondaries from hadronic interactions and photon/lepton DIS...
-  if (tkstatus == fStopAndKill)
-  {
-    // fEventAction->FillSecondaries(step);
-  }
+  // if (tkstatus == fStopAndKill)
+  //{
+  //  // fEventAction->FillSecondaries(step);
+  //}
 
   // if (charge == 0.0)
   //{
@@ -159,9 +159,45 @@ void B4bSteppingAction::UserSteppingAction(const G4Step *step)
   {
     // std::cout<<"Stepping Action:  track goes outside the world volume"<<std::endl;
     double eLeak = step->GetPostStepPoint()->GetKineticEnergy();
+    if (particle == G4Positron::Positron())
+    {
+      eLeak += 2 * electron_mass_c2;
+    }
     hh->accumulateEnergy(eLeak / GeV, -99);
   }
 
+  double e_lost = 0.0;
+  bool addPositronMass = false;
+  const std::vector<const G4Track *> *secondaries = step->GetSecondaryInCurrentStep();
+  for (auto sec : *secondaries)
+  {
+    e_lost += sec->GetKineticEnergy();
+    if (particle != G4Positron::Positron() && sec->GetParticleDefinition() == G4Positron::Positron())
+    {
+      e_lost += 2 * electron_mass_c2;
+      addPositronMass = true;
+    }
+  }
+  double e_change = step->GetPreStepPoint()->GetKineticEnergy() - step->GetPostStepPoint()->GetKineticEnergy() - e_lost - edep;
+  hh->accumulateDeposits(edep / GeV, step->GetTrack()->GetCurrentStepNumber());
+  if (fabs(e_change) >= 3 * electron_mass_c2)
+  {
+    std::cout << "Step number: " << step->GetTrack()->GetCurrentStepNumber()
+              << ", Particle name: " << particleName
+              << ", Number of secondaries: " << secondaries->size()
+              << ", Energy change: MeV " << e_change / MeV
+              << ", Energy lost to secondaries: " << e_lost
+              << ", Energy deposited: " << edep
+              << ", Delta energy: " << step->GetPreStepPoint()->GetKineticEnergy() - step->GetPostStepPoint()->GetKineticEnergy()
+              << ", Pre-step position: (" << step->GetPreStepPoint()->GetPosition().x() << ", " << step->GetPreStepPoint()->GetPosition().y() << ", " << step->GetPreStepPoint()->GetPosition().z() << ")"
+              << ", Post-step position: (" << step->GetPostStepPoint()->GetPosition().x() << ", " << step->GetPostStepPoint()->GetPosition().y() << ", " << step->GetPostStepPoint()->GetPosition().z() << ")"
+              << ", Pre-step time: " << step->GetPreStepPoint()->GetGlobalTime()
+              << ", Post-step time: " << step->GetPostStepPoint()->GetGlobalTime()
+              << ", Pre-step momentum: (" << step->GetPreStepPoint()->GetMomentum().x() << ", " << step->GetPreStepPoint()->GetMomentum().y() << ", " << step->GetPreStepPoint()->GetMomentum().z() << ")"
+              << ", Post-step momentum: (" << step->GetPostStepPoint()->GetMomentum().x() << ", " << step->GetPostStepPoint()->GetMomentum().y() << ", " << step->GetPostStepPoint()->GetMomentum().z() << ")"
+              << std::endl;
+  }
+
   if (thisName.compare(0, 5, "World") == 0)
   {
     // outside the volume

sim/src/CaloTree.cc

@@ -21,6 +21,7 @@
 #include "TPaveText.h"
 #include "TText.h"
 #include "TTree.h"
+#include <numeric>
 
 #include "CaloHit.h"
 #include "CaloID.h"
@@ -343,6 +344,13 @@ void CaloTree::EndEvent()
     m_nhitstruth = m_pidtruth.size();
     //
     tree->Fill();
+    std::cout << "Look into energy deposition in the calorimeter..." << std::endl;
+    std::cout << " total absolute energy deposits " << std::accumulate(m_eDEPs.begin(), m_eDEPs.end(), 0.0) << std::endl;
+    for (int i = 0; i < 200; i++)
+    {
+      if (m_eDEPs[i] > 0.0)
+        std::cout << "  i=" << i << "  eDEP=" << m_eDEPs[i] << std::endl;
+    }
   } //  end of if((eventCounts-1)<getParamI("eventsInNtupe"))
 
   //   analyze this event.
@@ -413,6 +421,12 @@ void CaloTree::clearCaloTree()
   m_eCentruth = 0.0;
   m_eScintruth = 0.0;
 
+  m_eDEPs.clear();
+  for (unsigned int i = 0; i < 200; i++)
+  {
+    m_eDEPs.push_back(0.0);
+  }
+
   m_nhits3dSS = 0;
   m_id3dSS.clear();
   m_type3dSS.clear();
@@ -527,6 +541,15 @@ void CaloTree::accumulateEnergy(double edep, int type = 0)
     m_eCentruth += edep;
 }
 
+void CaloTree::accumulateDeposits(double edep, int step)
+{
+  if (step < 0)
+    step = 0;
+  if (step >= 200)
+    step = 199;
+  m_eDEPs[step] += edep;
+}
+
 // ########################################################################
 void CaloTree::analyze()
 {