ch5/src/RunAction.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // gpaterno, October 2025
0027 //
0028 /// \file RunAction.cc
0029 /// \brief Implementation of the RunAction class
0030 //
0031 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0032
0033 #include "RunAction.hh"
0034 #include "RunActionMessenger.hh"
0035 #include "DetectorConstruction.hh"
0036 #include "Run.hh"
0037
0038 #include "G4AnalysisManager.hh"
0039 #include "G4RunManager.hh"
0040 #include "G4Run.hh"
0041 #include "G4RegionStore.hh"
0042 #include "G4AccumulableManager.hh"
0043 #include "G4LogicalVolumeStore.hh"
0044 #include "G4LogicalVolume.hh"
0045 #include "G4UnitsTable.hh"
0046 #include "G4SystemOfUnits.hh"
0047
0048 #include <iostream>
0049 #include <fstream>
0050 #include <string>
0051
0052 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0053
0054 RunAction::RunAction():G4UserRunAction()
0055 {
0056     G4RunManager::GetRunManager()->SetPrintProgress(1000);
0057
0058     fMessenger = new RunActionMessenger(this);
0059
0060     //Create the analysis manager
0061     fAnalysisManager = G4AnalysisManager::Instance();
0062     fAnalysisManager->SetDefaultFileType("root");
0063     fAnalysisManager->SetFileName(fFileName);
0064 #ifdef G4MULTITHREADED
0065     fAnalysisManager->SetNtupleMerging(true);
0066 #else
0067     fAnalysisManager->SetNtupleMerging(false);
0068 #endif
0069     fAnalysisManager->SetVerboseLevel(0);
0070
0071     //Creating the ntuple to score the particles impinging on the virtual screens
0072     fAnalysisManager->CreateNtuple("scoring_ntuple","virtual scoring screens");
0073     fAnalysisManager->CreateNtupleIColumn("screenID");
0074     fAnalysisManager->CreateNtupleSColumn("particle");
0075     fAnalysisManager->CreateNtupleDColumn("x");
0076     fAnalysisManager->CreateNtupleDColumn("y");
0077     fAnalysisManager->CreateNtupleDColumn("px");
0078     fAnalysisManager->CreateNtupleDColumn("py");
0079     fAnalysisManager->CreateNtupleDColumn("pz");
0080     fAnalysisManager->CreateNtupleDColumn("t");
0081     fAnalysisManager->CreateNtupleIColumn("eventID");
0082     fAnalysisManager->FinishNtuple();
0083
0084     //Creating the ntuple to score Edep in the Radiatior
0085     fAnalysisManager->CreateNtuple("edep_rad","radiator");
0086     fAnalysisManager->CreateNtupleDColumn("edep");
0087     fAnalysisManager->CreateNtupleIColumn("eventID");
0088     fAnalysisManager->FinishNtuple();
0089
0090     //Creating the ntuple to score Edep in the Converter
0091     fAnalysisManager->CreateNtuple("edep_conv","converter");
0092     fAnalysisManager->CreateNtupleDColumn("edep");
0093     fAnalysisManager->CreateNtupleIColumn("eventID");
0094     fAnalysisManager->FinishNtuple();
0095
0096     //Creating the ntuple to score Edep in the spheres of the Granular Target
0097     fAnalysisManager->CreateNtuple("edep_spheres","granular target");
0098     fAnalysisManager->CreateNtupleIColumn("volumeID");
0099     fAnalysisManager->CreateNtupleDColumn("edep");
0100     fAnalysisManager->CreateNtupleIColumn("eventID");
0101     fAnalysisManager->FinishNtuple();
0102
0103     //Creating the ntuple to score the particles exiting the crystals
0104     fAnalysisManager->CreateNtuple("scoring_ntuple2","particle leaving crystals");
0105     fAnalysisManager->CreateNtupleSColumn("particle");
0106     fAnalysisManager->CreateNtupleDColumn("x");
0107     fAnalysisManager->CreateNtupleDColumn("y");
0108     fAnalysisManager->CreateNtupleDColumn("z");
0109     fAnalysisManager->CreateNtupleDColumn("px");
0110     fAnalysisManager->CreateNtupleDColumn("py");
0111     fAnalysisManager->CreateNtupleDColumn("pz");
0112     fAnalysisManager->CreateNtupleDColumn("t");
0113     fAnalysisManager->CreateNtupleIColumn("eventID");
0114     fAnalysisManager->CreateNtupleIColumn("trackID");
0115     fAnalysisManager->FinishNtuple();
0116 }
0117
0118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0119
0120 G4Run* RunAction::GenerateRun() {return new Run;}
0121
0122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0123
0124 void RunAction::BeginOfRunAction(const G4Run* run)
0125 {
0126     //merging is only for root
0127     if (fFileName.find(".csv") != std::string::npos) {
0128         fAnalysisManager->SetNtupleMerging(false);
0129     }
0130
0131     //open the output file
0132     if (!fIsFileOpened) {
0133         fAnalysisManager->OpenFile(fFileName);
0134         fIsFileOpened = true;
0135     }
0136
0137     //print begin message
0138     if (IsMaster()) {
0139         G4int NumberOfEventToBeProcessed = run->GetNumberOfEventToBeProcessed();
0140
0141         G4cout
0142         << G4endl
0143         << "--------------------Begin of Global Run-----------------------"
0144         << G4endl
0145         << "Number of events to be processed: " << NumberOfEventToBeProcessed
0146         << G4endl
0147         << "--------------------------------------------------------------"
0148         << G4endl;
0149     }
0150 }
0151
0152 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0153
0154 void RunAction::EndOfRunAction(const G4Run* run)
0155 {
0156     //Write results on the root file
0157     if (fIsFileOpened) {
0158         fAnalysisManager->Write();
0159         fAnalysisManager->CloseFile();
0160
0161         G4int threadID = G4Threading::G4GetThreadId();
0162         if (threadID > 0) {
0163             G4cout << "writing results of thread " << G4Threading::G4GetThreadId()
0164                    << " to root file: " << fFileName << G4endl;
0165         }
0166     }
0167
0168     //Get the number of events
0169     G4int nofEvents = run->GetNumberOfEvent();
0170     if (nofEvents == 0) return;
0171
0172     //Get results from the Edep (in the Radiator Crystal) scorer
0173     const Run* aRun = static_cast<const Run*>(run);
0174     G4double Edep  = aRun->GetEdep();
0175     G4double Edep2 = aRun->GetEdep2();
0176     G4int nGoodEvents = aRun->GoodEvents();
0177     G4double stdEdep = std::sqrt(Edep2 - Edep*Edep/nGoodEvents);
0178
0179     //Print and write results to a text file
0180     if (IsMaster()) {
0181         //detector construction instance
0182         const DetectorConstruction* detectorConstruction =
0183             static_cast<const DetectorConstruction*>(G4RunManager::GetRunManager()
0184                 ->GetUserDetectorConstruction());
0185
0186         G4String matName =
0187             detectorConstruction->GetCrystalVolume()->GetMaterial()->GetName();
0188
0189         //set the significant digits to print
0190         G4cout.precision(8);
0191
0192         //print the results on screen
0193         G4cout << G4endl
0194         << "--------------------End of Global Run-----------------------"
0195         << G4endl
0196         << " The run had " << nofEvents << " events";
0197         G4cout << G4endl
0198         << " Edep in " << "the Radiator Crystal" << " (made of "
0199         << matName << "): " << Edep/MeV
0200         << " +/- " << stdEdep/MeV << " MeV" << G4endl
0201         << "------------------------------------------------------------"
0202         << G4endl << G4endl;
0203     }
0204
0205     //end message
0206     G4cout
0207         << G4endl
0208         << " The run consisted of " << nofEvents << " particles" << G4endl
0209         << "------------------------------------------------------------"
0210         << G4endl << G4endl;
0211 }
0212
0213 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0214
0215 void RunAction::SetFileName(G4String filename)
0216 {
0217     if (filename != "") fFileName = filename;
0218 }
0219
0220 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0221