Dear Geant4 Users
I am working with the clustering example of the DNA folder (Geant4-DNA examples).
This example provides multiple output (singleSSB, complexSSB, …) in a histogram form. I would like to extract this information as a root file in form a tuple.
I included some lines to create and fill tuples but the output does not agree with the histograms, I would like to have some help to identify the error that I am making.
This is the RunAccion.cc file:
The main difference is the CreateNtuple() function
//===================================
/// \file RunAction.cc
/// \brief Implementation of the RunAction class
#include “Analysis.hh”
#include “RunAction.hh”
#include “RunInitObserver.hh”
#include “RunActionMessenger.hh”
#include “G4Run.hh”
#include “G4RunManager.hh”
//…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…
RunAction::RunAction()
: G4UserRunAction()
{
fFileName=“clusters_output”;
fpRunMessenger = new RunActionMessenger(this);
CreateHistogram();
CreateNtuple();
}
//…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…
RunAction::~RunAction()
{
delete fpRunMessenger;
//delete G4AnalysisManager::Instance();
}
//…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…
void RunAction::BeginOfRunAction(const G4Run*)
{
//
RunInitManager::Instance()->Initialize();
// Get analysis manager
G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
// Open an output file
analysisManager->OpenFile(fFileName);
G4cout << "\n----> Histogram file is opened in " <<
fFileName << “.” << analysisManager->GetFileType() << G4endl;
}
//…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…
void RunAction::CreateHistogram()
{
G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
//analysisManager->SetHistoDirectoryName(“histograms”);
analysisManager->SetVerboseLevel(3);
analysisManager->SetFirstHistoId(1);
analysisManager->CreateH1(“1”,“simpleSSB”,75,0.,75);
analysisManager->CreateH1(“2”,“complexSSB”,75,0.,75);
analysisManager->CreateH1(“3”,“DSB”,75,0.,75);
analysisManager->CreateH1(“4”,“cluster size”,20,1,20);
analysisManager->CreateH1(“5”,“edep”,1,0.,1E6);
}
void RunAction::CreateNtuple()
{
G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
//analysisManager->SetNtupleDirectoryName(“ntuple”);
//analysisManager->SetNtupleMerging(true);
analysisManager->SetVerboseLevel(3);
analysisManager->SetFirstNtupleId(1);
analysisManager->CreateNtuple(“Clustering”, “Scoring”);
analysisManager->CreateNtupleIColumn(1,“singleSSB”);
analysisManager->CreateNtupleIColumn(1,“complexSSB”);
analysisManager->CreateNtupleIColumn(1,“DSB”);
analysisManager->CreateNtupleIColumn(1,“Dose”);
analysisManager->FinishNtuple(1);
}
//…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…
void RunAction::EndOfRunAction(const G4Run*)
{
//G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
//analysisManager->Write();
//analysisManager->CloseFile();
WriteHT();
}
//…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…
void RunAction::WriteHT()
{
G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
analysisManager->Write();
analysisManager->CloseFile();
analysisManager->SetVerboseLevel(3);
}
This is the EventAction .cc file:
The diference with the example is the funtion FillNtuple(), which is incharge to fill the tuple.
//====================================
#include “EventAction.hh”
#include “Analysis.hh”
#include “ClusteringAlgo.hh”
#include “G4Event.hh”
#include “G4LogicalVolume.hh”
#include “G4LogicalVolumeStore.hh”
#include “G4SystemOfUnits.hh”
#include “Randomize.hh”
//#include “G4VAnalysisManager.hh”
//…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…
EventAction::EventAction():G4UserEventAction()
{
//default parameter values
fEdep=0.;
// Create clustering algorithm
// These default values have been tuned for the Physics List G4EmDNAPhysics
// to reproduce data published by:
// Francis et al. 2011 Comput. Meth. Programs. Biomed. 2011 101(3)
fpClustering = new ClusteringAlgo(3.3nanometer,2,0.2,5eV,37.5*eV);
}
//…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…
EventAction::~EventAction()
{
delete fpClustering;
}
//…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…
void EventAction::BeginOfEventAction( const G4Event*)
{
fEdep=0.;
fpClustering->Purge();
}
//…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…
void EventAction::EndOfEventAction( const G4Event*)
{
FillNtuple();
std::map<G4int,G4int> sizeDistribution = fpClustering->RunClustering();
G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
//analysisManager->SetVerboseLevel(3);
analysisManager->FillH1(1,fpClustering->GetSSB());
analysisManager->FillH1(2,fpClustering->GetComplexSSB());
analysisManager->FillH1(3,fpClustering->GetDSB());
while ( !sizeDistribution.empty() )
{
analysisManager->FillH1(4,
sizeDistribution.begin()->first,
sizeDistribution.begin()->second);
sizeDistribution.erase(sizeDistribution.begin());
}
analysisManager->FillH1(5,
(fEdep/joule)/
(G4LogicalVolumeStore::GetInstance()->
GetVolume(“Target”)->GetMass()/kg));
}
void EventAction::FillNtuple()
{
singleSSB = fpClustering->GetSSB();
complexSSB = fpClustering->GetComplexSSB();
DSB = fpClustering->GetDSB();
Dose = (fEdep/joule)/(G4LogicalVolumeStore::GetInstance()->GetVolume(“Target”)->GetMass()/kg);
G4cout << "\n----> SingleSSB " << singleSSB << G4endl;
G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
analysisManager->SetVerboseLevel(3);
analysisManager->FillNtupleIColumn(1,singleSSB);
analysisManager->FillNtupleIColumn(1,complexSSB);
analysisManager->FillNtupleIColumn(1,DSB);
analysisManager->FillNtupleIColumn(1,Dose);
analysisManager-> AddNtupleRow(1);
}
//=======================================