Hello everyone,
I created a boxed mesh water phantom to get energy deposition in that. I created the PhantomHit and PhantomSD classes to implement this purpose.
Now, I want to save results to root files. How can I do it?
Thank you in advance.
Does looking into example 4c help? There could be some inspiration…
Thank you for your response.
I tried to create a RunAction and EventAction to process the results.
RunAction.cc
#include "RunAction.hh"
#include "G4UnitsTable.hh"
#include "G4SystemOfUnits.hh"
#include "G4AnalysisManager.hh"
#include "G4RunManager.hh"
#include <fstream>
#include <iostream>
#include "G4ios.hh"
#include <ctime>
using namespace std;
RunAction::RunAction()
{
G4RunManager::GetRunManager()->SetPrintProgress(1);
auto analysisManager = G4AnalysisManager::Instance();
analysisManager->SetVerboseLevel(1);
analysisManager->SetNtupleMerging(true);
// Creating ntuple
//
analysisManager->CreateNtuple("BoxMesh", "Depositerd Energy");
analysisManager->CreateNtupleIColumn("xID");
analysisManager->CreateNtupleIColumn("yID");
analysisManager->CreateNtupleIColumn("zID");
analysisManager->CreateNtupleDColumn("eDep");
analysisManager->FinishNtuple();
}
double RunAction::diffclock(clock_t clock1, clock_t clock2)
{
double diffticks = clock1 - clock2;
double diffms = (diffticks*1000)/CLOCKS_PER_SEC;
return diffms;
}
void RunAction::BeginOfRunAction(const G4Run*)
{
begin = clock();
auto analysisManager = G4AnalysisManager::Instance();
G4String fileName = "data.root";
analysisManager->OpenFile(fileName);
}
void RunAction::EndOfRunAction(const G4Run* run)
{
end = clock();
auto analysisManager = G4AnalysisManager::Instance();
// Print
//
if (IsMaster())
{
TextFile.open("time.txt", std::ios::out);
TextFile << " Elapsed time in seconds: " << (double(diffclock(end,begin)/1000)) << " s" << G4endl;
TextFile.close();
}
// save histograms & ntuple
//
analysisManager->Write();
analysisManager->CloseFile();
}
EventAction.cc
#include "EventAction.hh"
#include "PhantomSD.hh"
#include "DetectorConstruction.hh"
#include "PhantomHit.hh"
#include "G4AnalysisManager.hh"
#include "G4RunManager.hh"
#include "G4Event.hh"
#include "G4SDManager.hh"
#include "G4HCofThisEvent.hh"
#include "G4UnitsTable.hh"
#include "Randomize.hh"
#include <iomanip>
PhantomHitsCollection*
EventAction::GetHitsCollection(G4int hcID,
const G4Event* event) const
{
auto hitsCollection = static_cast<PhantomHitsCollection*>(event->GetHCofThisEvent()->GetHC(hcID));
return hitsCollection;
}
void EventAction::BeginOfEventAction(const G4Event* event)
{
// Get hits collections IDs
if (fHCID == -1)
{
fHCID = G4SDManager::GetSDMpointer()->GetCollectionID("PhantomROIHitsCollection");
}
}
void EventAction::EndOfEventAction(const G4Event* event)
{
// Get hits collections
auto absoHC = GetHitsCollection(fHCID, event);
G4int entries = absoHC->entries();
//G4cout << "Number of entries: " << entries << G4endl;
// Get analysis manager
auto analysisManager = G4AnalysisManager::Instance();
G4RunManager *runManager = G4RunManager::GetRunManager();
//const G4Event *event = runManager->GetCurrentEvent();
G4int Total_Events_To_Be_Processed = runManager->GetNumberOfEventsToBeProcessed();
G4int event_id = event->GetEventID();
int NumberOfThreads = G4Threading::G4GetNumberOfCores();
for (int i = G4Threading::GetNumberOfRunningWorkerThreads(); i >= 1; i--)
{
if(event_id == (i*(Total_Events_To_Be_Processed/NumberOfThreads)-1))
{
for (G4int i=0; i<entries; i++)
{
auto absoHit = (*absoHC)[i];
G4int xID = absoHit->GetXID();
G4int yID = absoHit->GetYID();
G4int zID = absoHit->GetZID();
G4double energyDep = absoHit->GetEdep();
vecX.push_back(xID);
vecY.push_back(yID);
vecZ.push_back(zID);
vecEnDep.push_back(energyDep);
}
G4cout << "Vector size: " << vecEnDep.size() << G4endl;
if(vecEnDep.size()>0)
{
for(int hitCount=0; hitCount<(int)vecX.size(); hitCount++)
{
analysisManager->FillNtupleIColumn(0, vecX[hitCount]);
analysisManager->FillNtupleIColumn(1, vecY[hitCount]);
analysisManager->FillNtupleIColumn(2, vecZ[hitCount]);
analysisManager->FillNtupleDColumn(3, vecEnDep[hitCount]/MeV);
analysisManager->AddNtupleRow();
}
}
}
}
}
But what I want is to back up the correct energy deposition result with the index (xID, yID, zID) of each voxel like an excel sheet of 4 columns. I don’t know if my above method is correct?
You’ve got it – “like an excel sheet” is exactly what an N-tuple (“TTree” in ROOT language) is: each “column” (or “TBranch” in ROOT) is one of the variables, and each row is one entry for all of the correlated values of those variables – an event, or a hit, or whatever.
Whether you’ve got all the details of your specific code right, I leave to you, but you have just the right picture.
Thank you,
I have written the results in the root file. I plotted the results, but I think it’s incorrect.
Because the number of voxels is 8000 (20x20x20), while the number of entries in the root file is 40000.
And there was nothing that appeared on the canvas. I don’t have much experience with the root, can you check it for me?
void ana()
{
TChain chain("BoxMesh");
chain.Add("data.root");
int dxPos;
int dyPos;
int dzPos;
double eDep;
chain.SetBranchAddress("xID", &dxPos);
chain.SetBranchAddress("yID", &dyPos);
chain.SetBranchAddress("zID", &dzPos);
chain.SetBranchAddress("eDep", &eDep);
Int_t nevent = chain.GetEntries();
cout << nevent << endl;
TH3F * hist = new TH3F("hist","hist", 20,-20, 20, 20, -20, 20, 20, -20, 20);
for (Int_t i=0;i<nevent;i++)
{
chain.GetEvent(i);
hist->Fill(dxPos, dyPos, dzPos, eDep);
}
hist->SetFillColor(kCyan);
hist->Draw("iso");
}
You seem to be doing some odd stuff in your EventAction.
First, why do you copy the hits into a bunch of vectors, then loop over the vectors to fill the N-tuple rows? Why not just put the analysisManger->FillNtuple... lines directly into the loop over hit entries, and get rid of the vectors entirely?
Second, why the loop over threads? The EndOfEventAction is called separate for each individual thread, and you don’t have to deal with any of that stuff yourself. G4AnalysisManager will take care of the thread merging, because you set it up that way in your RunAction.
Third, why the conditional check on event ID? Why not just write out an N-tuple row for each event unconditionally?
I intend to write the final results after the last event has finished. It will be a waste of time if I write all hits (each voxel created a hit per event). So I added thread and event conditions to the EventAction class. I attempted to save all hits per event, but the simulation run time was extremely long. Do you have any other ideas?
Can I use a function to accumulate values per event and write to the root file in the EndOfRunAction method?
Are you wanting event-by-event values (for instance, to see how the shape changes per event)? Or do you just want to accumulate final totals?
If you want per-event, I would suggest the following: git rid of the thread stuff and separate vectors, you don’t really need them. In your loop over hits, only record entries where Edep > 0. That should reduce the number of rows from 8000 per event to many fewer. If you want to do it this way, I’d recommend adding a column for eventID, and you can make plots per event, or sum everything in ROOT for the whole job.
If you want to accumulate totals for the run, this isn’t the best way. Instead, create a container or something in your RunAction, and pass a reference to it to your EventAction (some people pass the pointer to RunAction into their EventAction constructor, and save it in a data member). At the end of each event, accumulate the Edep from your voxels into that container.
For example, in your RunAction.hh, create a container with an associated accumulator function:
typedef std::array<G4int,3> voxel;
std::map<voxel, G4double> Etotal;
void AddEdep(G4int xID, G4int yID, G4int zID, G4double Edep) {
Etotal[voxel(xID,yID,zID)] += Edep;
}
Then, in your EventAction.cc, assuming you store RunAction* ra, loop over your hits and do
ra->AddEdep(absoHit->GetXID(), absoHit->GetYID(), absoHit->GetZID(), absoHit->GetEdep();
Finally, in RunAction::EndOfRunAction, you can loop over the entries in the std::map<> Etotal, pull apart each key into (x,y,z) indices, and store those indices and the energy into your N-tuple.
This is done systematically in most of electromagnetic or hadronic examples …
The simplest : TestEm4. Then : TestEm1 … etc …
Thanks, @mkelsey, I have implemented my simulation in both of the above ways.
But I still have a bit of an error related to merging the root file.
In my RunAction, I summed up the results from each event
#include "RunAction.hh"
#include "G4UnitsTable.hh"
#include "G4SystemOfUnits.hh"
#include "G4AnalysisManager.hh"
#include "G4RunManager.hh"
#include <fstream>
#include <iostream>
#include "G4ios.hh"
#include <ctime>
using namespace std;
RunAction::RunAction()
{
G4RunManager::GetRunManager()->SetPrintProgress(1);
auto analysisManager = G4AnalysisManager::Instance();
//analysisManager->SetVerboseLevel(1);
analysisManager->SetNtupleMerging(true);
// Creating histograms
analysisManager->CreateH1("Eabs", "Edep spectrum in phantom", 10000, 0 * MeV, 10 * MeV);
// Creating ntuple
//
analysisManager->CreateNtuple("BoxMesh", "Depositerd Energy");
analysisManager->CreateNtupleIColumn("xID");
analysisManager->CreateNtupleIColumn("yID");
analysisManager->CreateNtupleIColumn("zID");
analysisManager->CreateNtupleDColumn("eDep");
//analysisManager->CreateNtupleIColumn("eventID");
analysisManager->FinishNtuple();
Voxels = new sVoxel**[40];
for (int ix=0; ix<40; ix++)
{
Voxels[ix] = new sVoxel*[40];
for (int iy=0; iy<40; iy++)
{
Voxels[ix][iy] = new sVoxel[40];
for (int iz=0; iz<40; iz++)
{
Voxels[ix][iy][iz].voxelID = 0;
Voxels[ix][iy][iz].eDep = 0.;
Voxels[ix][iy][iz].nEvents = 0;
}
}
}
}
double RunAction::diffclock(clock_t clock1, clock_t clock2)
{
double diffticks = clock1 - clock2;
double diffms = (diffticks*1000)/CLOCKS_PER_SEC;
return diffms;
}
void RunAction::BeginOfRunAction(const G4Run*)
{
begin = clock();
auto analysisManager = G4AnalysisManager::Instance();
G4String fileName = "data.root";
analysisManager->OpenFile(fileName);
}
void RunAction::EndOfRunAction(const G4Run* run)
{
end = clock();
auto analysisManager = G4AnalysisManager::Instance();
// Print
//
if (IsMaster())
{
TextFile.open("time.txt", std::ios::out);
TextFile << " Elapsed time in seconds: " << (double(diffclock(end,begin)/1000)) << " s" << G4endl;
TextFile.close();
}
for(int xID=0; xID<40; xID++)
{
for(int yID=0; yID<40; yID++)
{
for(int zID=0; zID<40; zID++)
{
analysisManager->FillNtupleIColumn(0, xID);
analysisManager->FillNtupleIColumn(1, yID);
analysisManager->FillNtupleIColumn(2, zID);
analysisManager->FillNtupleDColumn(3, Voxels[xID][yID][zID].eDep/MeV);
analysisManager->AddNtupleRow();
}
}
}
// save histograms & ntuple
//
analysisManager->Write();
analysisManager->CloseFile();
}
void RunAction::AddEdep(G4int pXID, G4int pYID, G4int pZID, G4double pEDep)
{
Voxels[pXID][pYID][pZID].voxelID = 0;
Voxels[pXID][pYID][pZID].eDep += pEDep;
Voxels[pXID][pYID][pZID].nEvents += 1;
}
If I used
analysisManager->SetNtupleMerging(true);
There is an error
-------- WWWW ------- G4Exception-START -------- WWWW -------
*** G4Exception : Analysis_W001
issued by : G4TNtupleManager<NT,FT>::FillNtupleTColumn
Ntuple 0 does not exist.
*** This is just a warning message. ***
-------- WWWW -------- G4Exception-END --------- WWWW -------
if I didn’t use the command. There was no error. I can use the “hadd” command to sum up root files manually.