Hello,
My output consists in a single root file containing histograms and ntuples. I use both SteppingAction and SensitiveDetector to fill them. There a lot of Ntuples so I will omit some to keep the thread readable.
in the constructor of RunAction:
RunAction::RunAction(const string& configFileName, EventAction\* eventAction)
: fEventAction(eventAction)
{
DefineCommands();
// Create analysis manager
// The choice of the output format is done via the specified
// file extension.
auto analysisManager = G4AnalysisManager::Instance();
analysisManager->SetVerboseLevel(1);
analysisManager->SetNtupleMerging(true);
analysisManager->SetActivation(true);
// Book histograms, ntuples
//
// Creating histograms
analysisManager->CreateH1(“Ecalo”,“Edep in crystal”, 100, Ecalo_min\*MeV, Ecalo_max\*MeV); //0
analysisManager->CreateH1(“Lcalo”,“trackL in crystal”, 100, 0., 50\*m); //1
analysisManager->CreateH1(“Elyso”,“Edep in lyso”, 100, 0.\*MeV, 50.\*MeV); //2
analysisManager->CreateH1(“Eps”,“Edep in preshower”, 100, 0.\*MeV, 1000\*MeV); //3
analysisManager->CreateH1(“Ecalo_SA”,“Edep in crystal (SA)”, 100, Ecalo_min\*MeV, Ecalo_max\*MeV); //4
analysisManager->CreateH1(“N_scint_CH1”,“N photons scint in S”, 100, 0, 15000); //5
analysisManager->CreateH1(“N_cer_CH1”,“N photons cer in S”, 10, 0, 30); //6
analysisManager->CreateH1(“N_scint_CH2”,“N photons scint in C”, 100, 0, 100000); //7
analysisManager->CreateH1(“N_cer_CH2”,“N photons cer in C”, 100, 0, 2000); //8
// Create ntuples, first beam information and energy deposit information
analysisManager->CreateNtuple(“primaryEdep”, “Edep and TrackL”);
analysisManager->CreateNtupleDColumn(“primaryEnergy”); //0
analysisManager->CreateNtupleDColumn(“primaryX”); //1
analysisManager->CreateNtupleDColumn(“primaryY”); //2
analysisManager->CreateNtupleDColumn(“primaryZ”); //3
analysisManager->CreateNtupleDColumn(“primaryPx”); //4
analysisManager->CreateNtupleDColumn(“primaryPy”); //5
analysisManager->CreateNtupleDColumn(“primaryPz”); //6
analysisManager->CreateNtupleDColumn(“Ecalo”); //7
analysisManager->CreateNtupleDColumn(“Elyso”); //8
analysisManager->CreateNtupleDColumn(“Epreshower”); //9
analysisManager->CreateNtupleDColumn(“Lcalo”); //10
analysisManager->CreateNtupleDColumn(“Ecalo_SA”); //11
analysisManager->CreateNtupleIColumn(“eventID”); //12
analysisManager->FinishNtuple();
//ntuple for information about arrival of photons
analysisManager->CreateNtuple(“photoDet”, “photons arrival”);
analysisManager->CreateNtupleIColumn(“N_scint_CH1”); //0
analysisManager->CreateNtupleIColumn(“N_cer_CH1”); //1
analysisManager->CreateNtupleIColumn(“N_scint_CH2”); //2
analysisManager->CreateNtupleIColumn(“N_cer_CH2”); //3
analysisManager->CreateNtupleIColumn(“N_scint_calo”); //4
analysisManager->CreateNtupleIColumn(“N_cer_calo”); //5
//vector branches don’t need to be filled in EventAction
//time vectors
analysisManager->CreateNtupleFColumn(“time_scint_CH1”, fEventAction->m_time_scint_CH1);
analysisManager->CreateNtupleFColumn(“time_cer_CH1”, fEventAction->m_time_cer_CH1);
analysisManager->CreateNtupleFColumn(“time_scint_CH2”, fEventAction->m_time_scint_CH2);
analysisManager->CreateNtupleFColumn(“time_cer_CH2”, fEventAction->m_time_cer_CH2);
//more branches following
analysisManager->FinishNtuple();
-BeginOfRunAction()
void RunAction::BeginOfRunAction(const G4Run\* /\*run\*/)
{
//inform the runManager to save random number seed
//G4RunManager::GetRunManager()->SetRandomNumberStore(true);
// Get analysis manager
auto analysisManager = G4AnalysisManager::Instance();
// Open an output file
//
if (!fname.empty()) fname = fname;
G4String fileName = fpath + fname + “.root”;
analysisManager->OpenFile(fileName);
G4cout << "Using " << analysisManager->GetType() << G4endl;
}
-EndOfRunAction()
void RunAction::EndOfRunAction(const G4Run\* /\*run\*/)
{
auto analysisManager = G4AnalysisManager::Instance();
analysisManager->Write();
analysisManager->CloseFile();
}
Some information is saved by SteppingAction()
void SteppingAction::UserSteppingAction(const G4Step\* step)
{
// get volume of the current step
auto volume = step->GetPreStepPoint()->GetTouchableHandle()->GetVolume();
auto theTrack = step->GetTrack();
auto theTrackInfo = dynamic_cast<TrackInformation\*>(theTrack->GetUserInformation());
//Beam information
if(G4StrUtil::contains(volume->GetName(), “World”)){
if ( theTrack->GetParentID() == 0 ) // select only the primary particle
{
G4StepPoint\* thePrePoint = step->GetPreStepPoint();
const G4ThreeVector& thePrePosition = thePrePoint->GetPosition();
G4VPhysicalVolume\* thePrePV = thePrePoint->GetPhysicalVolume();
G4String thePrePVName = "";
if (thePrePV) thePrePVName = thePrePV->GetName();
}
}
SensitiveDetectors are defined for SiPM volumes, which is not the volume where I defined the UserLimits:
G4bool SiPMSD::ProcessHits(G4Step\* aStep, G4TouchableHistory\*)
{
if (!aStep) return false;
G4Track\* theTrack = aStep->GetTrack();
//check if it's an optical photon, SiPM only measure optical photons
if (theTrack->GetDefinition() != G4OpticalPhoton::OpticalPhotonDefinition()) {
return false;
}
//get step points of the step
G4StepPoint\* thePrePoint = aStep->GetPreStepPoint();
G4StepPoint\* thePostPoint = aStep->GetPostStepPoint();
G4String process = theTrack->GetCreatorProcess()->GetProcessName();
auto theTrackInfo = dynamic_cast<TrackInformation\*>(theTrack->GetUserInformation());
G4int parentID = theTrack->GetParentID();
G4double photonEnergy = theTrack->GetTotalEnergy();
G4double photonTime = theTrack->GetGlobalTime();
G4int trackID = theTrack->GetTrackID();
G4ThreeVector photonPosition = thePrePoint->GetPosition(); //a step belongs to the volume of its preStepPoint
G4ThreeVector photonMomentum = thePrePoint->GetMomentumDirection(); //a step belongs to the volume of its preStepPoint
G4String parentName = theTrackInfo->GetParentName();
G4double prodEnergy = theTrackInfo->GetParticleEnergy();
G4double prodTime = theTrackInfo->GetParticleTime();
G4ThreeVector prodPosition = theTrackInfo->GetParticlePosition();
G4ThreeVector prodMomentum = theTrackInfo->GetParticleMomentumDirection();
//define a hit and set its members
SiPMHit\* hit = new SiPMHit();
hit->SetCreatorProcess(process);
hit->SetEnergy(photonEnergy);
hit->SetTime(photonTime);
hit->SetTrackID(trackID);
hit->SetPos(photonPosition);
hit->SetMomentum(photonMomentum);
hit->SetParentName(parentName);
hit->SetProdEnergy(prodEnergy);
hit->SetProdTime(prodTime);
hit->SetProdPosition(prodPosition);
hit->SetProdMomentum(prodMomentum);
//add the hit to the hit collection
fSiPMHitCollection->insert(hit);
return true;
}
Thank you for your time