From c99d6d17fb4e0fa3926bd48a8aebab0f96c8b6d6 Mon Sep 17 00:00:00 2001
From: Chad Lantz <clantz2@illinois.edu>
Date: Wed, 10 Jul 2024 22:34:15 -0500
Subject: [PATCH] Update to Geant4.11.2.2 by removing legacy G4AnalysisManager
 and made things thread safe

---
 include/ASCIIPrimaryGenerator.hh |  1 -
 include/DetectorConstruction.hh  |  1 +
 include/RunAction.hh             |  2 +-
 include/lgAnalysis.hh            | 36 -------------------
 src/ASCIIPrimaryGenerator.cc     |  4 ---
 src/DetectorConstruction.cc      | 21 +++++------
 src/DetectorMessenger.cc         |  4 +--
 src/EventAction.cc               | 22 ++++++------
 src/PhysicsList.cc               |  2 +-
 src/PrimaryGeneratorAction.cc    | 14 +++-----
 src/RunAction.cc                 | 62 ++++++++++++++++----------------
 11 files changed, 64 insertions(+), 105 deletions(-)
 delete mode 100644 include/lgAnalysis.hh

diff --git a/include/ASCIIPrimaryGenerator.hh b/include/ASCIIPrimaryGenerator.hh
index 4bb8a9b..82b9a20 100644
--- a/include/ASCIIPrimaryGenerator.hh
+++ b/include/ASCIIPrimaryGenerator.hh
@@ -35,7 +35,6 @@
 #define ASCIIPrimaryGenerator_h 1
 
 #include "G4VPrimaryGenerator.hh"
-#include "lgAnalysis.hh"
 
 #include <fstream>
 #include <vector>
diff --git a/include/DetectorConstruction.hh b/include/DetectorConstruction.hh
index 2072b14..7fd2b23 100644
--- a/include/DetectorConstruction.hh
+++ b/include/DetectorConstruction.hh
@@ -53,6 +53,7 @@ class DetectorConstruction : public G4VUserDetectorConstruction
 
   public:
     virtual G4VPhysicalVolume* Construct();
+    void ConstructSDandField();
 
     void BuildWorld           ();
     void BuildTrapezoidLG     ();
diff --git a/include/RunAction.hh b/include/RunAction.hh
index e48e114..d606d2c 100644
--- a/include/RunAction.hh
+++ b/include/RunAction.hh
@@ -55,7 +55,7 @@ class RunAction : public G4UserRunAction
     virtual void EndOfRunAction(const G4Run* aRun);
 
     std::vector< std::vector<double>* > GetVectors(){ return fPtrVec; }
-    inline  void ClearVectors(){ for(auto vec : fPtrVec) vec->clear(); }
+    inline void ClearVectors(){ for(auto vec : fPtrVec) vec->clear(); }
 
   private:
     G4Timer* fTimer;
diff --git a/include/lgAnalysis.hh b/include/lgAnalysis.hh
deleted file mode 100644
index 7644179..0000000
--- a/include/lgAnalysis.hh
+++ /dev/null
@@ -1,36 +0,0 @@
-// ********************************************************************
-// * License and Disclaimer                                           *
-// *                                                                  *
-// * The  Geant4 software  is  copyright of the Copyright Holders  of *
-// * the Geant4 Collaboration.  It is provided  under  the terms  and *
-// * conditions of the Geant4 Software License,  included in the file *
-// * LICENSE and available at  http://cern.ch/geant4/license .  These *
-// * include a list of copyright holders.                             *
-// *                                                                  *
-// * Neither the authors of this software system, nor their employing *
-// * institutes,nor the agencies providing financial support for this *
-// * work  make  any representation or  warranty, express or implied, *
-// * regarding  this  software system or assume any liability for its *
-// * use.  Please see the license in the file  LICENSE  and URL above *
-// * for the full disclaimer and the limitation of liability.         *
-// *                                                                  *
-// * This  code  implementation is the result of  the  scientific and *
-// * technical work of the GEANT4 collaboration.                      *
-// * By using,  copying,  modifying or  distributing the software (or *
-// * any work based  on the software)  you  agree  to acknowledge its *
-// * use  in  resulting  scientific  publications,  and indicate your *
-// * acceptance of all terms of the Geant4 Software license.          *
-// ********************************************************************
-//
-//
-/// \file lgAnalysis.hh
-/// \brief Selection of the analysis technology
-
-#ifndef lgAnalysis_h
-#define lgAnalysis_h 1
-
-#include "g4root.hh"
-//#include "g4cvs.hh"
-//#include "g4xml.hh"
-
-#endif
diff --git a/src/ASCIIPrimaryGenerator.cc b/src/ASCIIPrimaryGenerator.cc
index c1adf2a..087604a 100644
--- a/src/ASCIIPrimaryGenerator.cc
+++ b/src/ASCIIPrimaryGenerator.cc
@@ -79,10 +79,6 @@ void ASCIIPrimaryGenerator::GeneratePrimaryVertex(G4Event* event)
 
   GetNextEvent();
 
-  // Pass the original event number through to the ouput
-  auto analysisManager = G4AnalysisManager::Instance();
-  analysisManager->FillNtupleIColumn(7, fEventNo );
-
   //For each photon read in, add a new particle and vertex to the event
   int nPhotons = fPositionVec->size();
   for(int i = 0; i < nPhotons; i++){
diff --git a/src/DetectorConstruction.cc b/src/DetectorConstruction.cc
index 8cb9a39..586c785 100644
--- a/src/DetectorConstruction.cc
+++ b/src/DetectorConstruction.cc
@@ -38,7 +38,6 @@
 #include <math.h>
 
 // GEANT4 //
-#include "G4Element.hh"
 #include "G4SDManager.hh"
 #include "G4RunManager.hh"
 #include "G4GeometryManager.hh"
@@ -164,7 +163,8 @@ void DetectorConstruction::BuildWorld(){
                       0,                  //copy number
                       checkOverlaps);     //overlaps checking
 
-  G4VisAttributes* boxVisAtt_world= new G4VisAttributes(G4VisAttributes::Invisible);
+  G4VisAttributes* boxVisAtt_world= new G4VisAttributes();
+  boxVisAtt_world ->SetVisibility(false);
 
 	m_logicWorld ->SetVisAttributes(boxVisAtt_world);
 
@@ -249,10 +249,6 @@ void DetectorConstruction::BuildTrapezoidLG( ){
  */
 void DetectorConstruction::BuildPMT(){
 
-  G4SDManager* SDman = G4SDManager::GetSDMpointer();
-  PMTSD* PMT = new PMTSD("MyPMT");
-  SDman->AddNewDetector( PMT );
-
   m_solidPMT =
     new G4Tubs("PMT",       //name
               0.0*mm,       //Inner radius
@@ -268,10 +264,17 @@ void DetectorConstruction::BuildPMT(){
 
   G4VisAttributes* VisAtt_PMT = new G4VisAttributes(G4Colour(1.0,1.0,0.6,0.7));
   m_logicPMT->SetVisAttributes(VisAtt_PMT);
-  m_logicPMT->SetSensitiveDetector( PMT );
+  
 
 }
 
+void DetectorConstruction::ConstructSDandField(){
+  G4SDManager* SDman = G4SDManager::GetSDMpointer();
+  PMTSD* PMT = new PMTSD("MyPMT");
+  SDman->AddNewDetector( PMT );
+  m_logicPMT->SetSensitiveDetector( PMT );
+}
+
 /*
 *
 */
@@ -528,7 +531,7 @@ void DetectorConstruction::UseCADModel(G4String fileName){
   m_Surfvec.clear();
 
 
-  G4String fileType = fileName.substr( fileName.last('.') + 1, fileName.size() - fileName.last('.'));
+  G4String fileType = fileName.substr( fileName.find_last_of('.') + 1, fileName.size() - fileName.find_last_of('.'));
 
   if(fileType == "stl"){
     auto mesh = CADMesh::TessellatedMesh::FromSTL((char*) fileName.c_str());
@@ -668,7 +671,6 @@ void DetectorConstruction::SetSurfaceSigmaAlpha(G4double v){
  * of the light guide
  */
 void DetectorConstruction::AddSurfaceMPV(const char* c, G4MaterialPropertyVector* mpv){
-  mpv->SetSpline(true);
   materials->GetMPTArray().at(1)->AddProperty(c, mpv);
   materials->AlSurface->SetMaterialPropertiesTable(materials->GetMPTArray().at(1));
   G4cout << "The MPT for the surface is now: " << G4endl;
@@ -681,7 +683,6 @@ void DetectorConstruction::AddSurfaceMPV(const char* c, G4MaterialPropertyVector
  * which contains the light guide
  */
 void DetectorConstruction::AddGasMPV(const char* c, G4MaterialPropertyVector* mpv){
-  mpv->SetSpline(true);
   m_GasMPT->AddProperty(c, mpv);
   G4cout << "The MPT for the gas is now: " << G4endl;
   m_GasMPT->DumpTable();
diff --git a/src/DetectorMessenger.cc b/src/DetectorMessenger.cc
index f80d887..e22e10c 100644
--- a/src/DetectorMessenger.cc
+++ b/src/DetectorMessenger.cc
@@ -419,7 +419,7 @@ void DetectorMessenger::SetNewValue(G4UIcommand* command,G4String newValue)
     // Convert string to physics vector
     // string format is property name, then pairs of energy, value
     // space delimited
-    G4MaterialPropertyVector* mpv = new G4MaterialPropertyVector();
+    G4MaterialPropertyVector* mpv = new G4MaterialPropertyVector(true);
     G4cout << newValue << G4endl;
     std::istringstream instring(newValue);
     G4String prop;
@@ -466,7 +466,7 @@ void DetectorMessenger::SetNewValue(G4UIcommand* command,G4String newValue)
     // Convert string to physics vector
     // string format is property name, then pairs of energy, value
     // space delimited
-    G4MaterialPropertyVector* mpv = new G4MaterialPropertyVector();
+    G4MaterialPropertyVector* mpv = new G4MaterialPropertyVector(true);
     G4cout << newValue << G4endl;
     std::istringstream instring(newValue);
     G4String prop;
diff --git a/src/EventAction.cc b/src/EventAction.cc
index 45c22c1..35837dc 100644
--- a/src/EventAction.cc
+++ b/src/EventAction.cc
@@ -25,18 +25,11 @@
 // @Author Chad Lantz
 #include "EventAction.hh"
 
-#include "G4VProcess.hh"
-
-#include "G4ParticleDefinition.hh"
-#include "G4ParticleTypes.hh"
-#include "G4SDManager.hh"
+#include "G4AnalysisManager.hh"
 #include "G4RunManager.hh"
 #include "G4Event.hh"
-#include "G4Track.hh"
 #include "G4ios.hh"
 
-#include "PMTHit.hh"
-#include "lgAnalysis.hh"
 #include "RunAction.hh"
 
 /*
@@ -68,9 +61,18 @@ void EventAction::BeginOfEventAction(const G4Event* event){
  */
 void EventAction::EndOfEventAction(const G4Event* event){
   // fill ntuple  //
+
   G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
+
+  RunAction *runAction = ((RunAction*)G4RunManager::GetRunManager()->GetUserRunAction());
+  std::vector< std::vector<double>* > ptrVec = runAction->GetVectors();
+
+  for(unsigned int i = 0; i < ptrVec[0]->size(); i++){
+    analysisManager->FillH2(0, ptrVec[0]->at(i), ptrVec[1]->at(i));
+  }
+
   analysisManager->FillNtupleIColumn( 0, event->GetEventID());
   analysisManager->AddNtupleRow();
-
-  ((RunAction*)G4RunManager::GetRunManager()->GetUserRunAction())->ClearVectors();
+  
+  runAction->ClearVectors();
 }
diff --git a/src/PhysicsList.cc b/src/PhysicsList.cc
index 4aceecb..30e3016 100644
--- a/src/PhysicsList.cc
+++ b/src/PhysicsList.cc
@@ -209,7 +209,7 @@ void PhysicsList::ConstructOp(){
   fCerenkovProcess->SetMaxBetaChangePerStep(10.0);
   fCerenkovProcess->SetTrackSecondariesFirst(true);
   fScintillationProcess = new G4Scintillation("Scintillation");
-  fScintillationProcess->SetScintillationYieldFactor(1.);
+  // fScintillationProcess->SetScintillationYieldFactor(1.);
   fScintillationProcess->SetTrackSecondariesFirst(true);
   fAbsorptionProcess = new G4OpAbsorption();
   fRayleighScatteringProcess = new G4OpRayleigh();
diff --git a/src/PrimaryGeneratorAction.cc b/src/PrimaryGeneratorAction.cc
index 749f98a..140ab36 100644
--- a/src/PrimaryGeneratorAction.cc
+++ b/src/PrimaryGeneratorAction.cc
@@ -25,14 +25,12 @@
 /// \file /src/PrimaryGeneratorAction.cc
 /// \brief Implementation of the PrimaryGeneratorAction class
 #include "PrimaryGeneratorAction.hh"
-#include "lgAnalysis.hh"
-
-#include "Randomize.hh"
 
 #include "G4Event.hh"
 #include "G4ParticleTable.hh"
 #include "G4ParticleDefinition.hh"
 #include "G4SystemOfUnits.hh"
+#include <G4String.hh>
 
 
 
@@ -84,10 +82,6 @@ void PrimaryGeneratorAction::GeneratePrimariesFromRootFile(G4Event* anEvent){
 
   fInputTree->GetEntry(evNo);
 
-  //Pass the original event number to the output in case they aren't sequential
-  auto analysisManager = G4AnalysisManager::Instance();
-  analysisManager->FillNtupleIColumn(7, anEvent->GetEventID() );
-
   G4int nPhotons = x->size();
   for(G4int i = 0; i < nPhotons; i++){
      G4PrimaryParticle* particle = new G4PrimaryParticle(particleDefinition);
@@ -106,8 +100,8 @@ void PrimaryGeneratorAction::SetInputFile(G4String _name){
     return;
   }
   G4String check = _name;
-  check.toLower();
-  if(check.contains(".root")){
+  check = G4StrUtil::to_lower_copy(check);
+  if(G4StrUtil::contains(check,".root")){
     fInputFile = new TFile( _name.c_str(), "READ");
 
     if(!fInputFile->IsOpen()){
@@ -131,7 +125,7 @@ void PrimaryGeneratorAction::SetInputFile(G4String _name){
     fInputTree->SetBranchAddress("Energy",&Energy);
     fInputTree->SetBranchAddress("Time",&time);
 
-  }else if(check.contains(".txt")){
+  }else if(G4StrUtil::contains(check,".txt")){
     fUseASCIIInput = true;
 
     G4cout << "Using " << _name << " as source" << G4endl;
diff --git a/src/RunAction.cc b/src/RunAction.cc
index aa87e78..b19eca6 100644
--- a/src/RunAction.cc
+++ b/src/RunAction.cc
@@ -28,7 +28,6 @@
 // @Author Chad Lantz
 
 
-// Make this appear first!
 #include "G4Timer.hh"
 
 #include "G4RunManager.hh"
@@ -36,9 +35,8 @@
 
 #include "RunAction.hh"
 #include "EventAction.hh"
-#include "lgAnalysis.hh"
 
-#include "G4VAnalysisManager.hh"
+#include "G4AnalysisManager.hh"
 #include "G4SystemOfUnits.hh"
 
 
@@ -52,38 +50,24 @@ RunAction::RunAction(G4String fileName)
 
   G4RunManager::GetRunManager()->SetPrintProgress(1);
 
-}
-
-/*
- *
- */
-RunAction::~RunAction(){
-  delete fTimer;
-  delete G4AnalysisManager::Instance();
-  for(auto vec : fPtrVec ) delete vec;
-}
-
-/*
- *
- */
-void RunAction::BeginOfRunAction(const G4Run* aRun){
-  G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
-  fTimer->Start();
-
-  // Create analysis manager. The choice of analysis
-  // technology is done via selection of a namespace
-  // in B4Analysis.hh
+  // Get analysis manager and set the type
   auto analysisManager = G4AnalysisManager::Instance();
-  G4cout << "Using " << analysisManager->GetType() << G4endl;
+  analysisManager->SetDefaultFileType("root");
+  analysisManager->SetVerboseLevel(1);
+  // Only merge in MT mode to avoid warning when running in Sequential mode
+#ifdef G4MULTITHREADED
+  analysisManager->SetNtupleMerging(true);
+#endif
 
   // Open an output file
-  if(m_fileName == "") m_fileName = "output";
+  if(m_fileName == "") m_fileName = "output.root";
   analysisManager->OpenFile(m_fileName);
+  
+  analysisManager->SetHistoDirectoryName("histo");
+  analysisManager->SetNtupleDirectoryName("ntuple");
 
+  analysisManager->CreateH2("position", "Emission Position of Detected Photons", 100, -39, 30, 100, -30, 30);
 
-  // Create directories
-  analysisManager->SetVerboseLevel(1);
-  analysisManager->SetNtupleMerging(true);
 
   // Creating ntuple
   // Populate the vector with valid pointers and assign them to branches
@@ -99,6 +83,24 @@ void RunAction::BeginOfRunAction(const G4Run* aRun){
   analysisManager->CreateNtupleDColumn( "theta",   *fPtrVec[5] );
   analysisManager->CreateNtupleDColumn( "energy",  *fPtrVec[6] );
   analysisManager->FinishNtuple();
+
+}
+
+/*
+ *
+ */
+RunAction::~RunAction(){
+  delete fTimer;
+  delete G4AnalysisManager::Instance();
+  for(auto vec : fPtrVec ) delete vec;
+}
+
+/*
+ *
+ */
+void RunAction::BeginOfRunAction(const G4Run* aRun){
+  fTimer->Start();
+
 }
 
 /*
@@ -113,5 +115,5 @@ void RunAction::EndOfRunAction(const G4Run* aRun){
   auto analysisManager = G4AnalysisManager::Instance();
   // save histograms & ntuple
   analysisManager->Write();
-  analysisManager->CloseFile();
+  analysisManager->CloseFile(false);
 }
-- 
GitLab