/**
* @file SolidRunAction.cc
* @author: (modified by) Ibrahin Pinera
* @date 2016 SoLid - University of Antwerp
*/
#include "SolidRunAction.hh"
#include "SolidRunMessenger.hh"
#include "SolidSteppingAction.hh"
#include "SolidPrimaryGeneratorAction.hh"
//#ifdef G4MULTITHREADED
//#include "G4MTRunManager.hh"
//#else
#include "G4RunManager.hh"
//#endif
#include "G4Run.hh"
#include "G4UImanager.hh"
#ifdef G4VIS_USE
#include "G4VVisManager.hh"
#include "G4VisExecutive.hh"
#endif
#include "G4ios.hh"
#include "Randomize.hh"
#include "TFile.h"
#include "TBranch.h"
#include "TLeaf.h"
#include "TTree.h"
#include "TH2I.h"
#include "TH1D.h"
#include "TAxis.h"
#include "TError.h"
#include "TDirectory.h"
#include "TROOT.h"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
SolidRunAction::SolidRunAction(TFile* outp_, SolidPrimaryGeneratorAction* pga)
{
//
outp = outp_;
genAction = pga;
G4RunManager::GetRunManager()->SetPrintProgress(1000);
fSeed = 123456789;
m_storeT2 = false;
m_storeT6 = false;
m_storeT7 = false;
m_storeT8 = false;
m_storeT9 = true;
m_storeLightT = true;
m_storeT9cut = 1;
m_raMess = new SolidRunMessenger(this);
lightOutput = 0;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
SolidRunAction::~SolidRunAction()
{;}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
void SolidRunAction::BeginOfRunAction(const G4Run* aRun)
{
//start the timer clock to calculate run times
start = time(NULL);
//
if (fSeed == -1) {
//set the random seed to the CPU clock
fSeed = start;
}
G4Random::setTheSeed(fSeed);
G4cout << "Random generator seed value="
<< G4Random::getTheSeed() << G4endl;
// set random seed for rand48() in SolidIBDKinematics.cc
srand48(fSeed);
// Go to output file
outp->cd();
// gROOT->ProcessLine("include <vector>"); //required for the vectors in T3 Branch
//create the ROOT tree for event data
T0 = new TTree("T0","Simulation setup information");
T0->Branch( "Event", &Stats.Evt, "Evt/I");
T0->Branch( "NumberCubes", &Stats.NCubes, "NCubes/I");
T0->Branch( "Seed", &fSeed, "Seed/L");
T0->Branch( "RealTime", &Stats.real_time, "RealTime/D");
T1 = new TTree("T1","Primary vertex");
T1->Branch( "pdg", &primary_pdgid);
T1->Branch( "energy", &primary_energy);
T1->Branch( "posx", &primary_posx);
T1->Branch( "posy", &primary_posy);
T1->Branch( "posz", &primary_posz);
T1->Branch( "momx", &primary_momx);
T1->Branch( "momy", &primary_momy);
T1->Branch( "momz", &primary_momz);
T1->Branch( "volid", &primary_vol);
//These 4 entries below should only be active for IBD SoLO simulations
if (genAction->IsIBDSoLoOn())
{
T1->Branch( "ltrue", &primary_ltrue);
T1->Branch( "rxtrx", &primary_rxtrx);
T1->Branch( "rxtry", &primary_rxtry);
T1->Branch( "rxtrz", &primary_rxtrz);
}
if( m_storeT2 )
{
T2 = new TTree("T2","ibd information");
T2->Branch( "LiCapture", &Stats.LiCapture, "LiCapture/I");
T2->Branch( "PVTCapture", &Stats.PVTCapture, "PVTCapture/I");
T2->Branch( "CaptureCube", &Stats.CapCube, "CaptureCube/I");
//T2->Branch( "DepositEnergy", &Stats.Deposit, "DepositE/F");
T2->Branch( "LiCaptureEnergy", &CapEneLi, "LiCaptureEnergy/D");
T2->Branch( "PvtCaptureEnergy",&CapEnePVT, "PvtCaptureEnergy/D");
T2->Branch( "CapPosX", &Stats.CapPos[0], "InitPosX/D");
T2->Branch( "CapPosY", &Stats.CapPos[1], "InitPosY/D");
T2->Branch( "CapPosZ", &Stats.CapPos[2], "InitPosZ/D");
// T2->Branch( "ThermalPosX", &Stats.ThermPos[0], "ThermalPosX/D");
// T2->Branch( "ThermalPosY", &Stats.ThermPos[1], "ThermalPosY/D");
// T2->Branch( "ThermalPosZ", &Stats.ThermPos[2], "ThermalPosZ/D");
// T2->Branch( "EneCube", &Stats.EneCube, "EneCube[NCubes]/F");
// T2->Branch( "CapE", &Stats.CapEne, "CapE/F");
// T2->Branch( "Detected", &Stats.Detected, "Detected/I");
// T2->Branch( "Generated1", &Stats.Generated1, "Generated1/I");
// T2->Branch( "Generated2", &Stats.Generated2, "Generated2/I");
// T2->Branch( "WLS", &Stats.WLS, "WLS/I");
T2->Branch( "ModerationTime", &ModTime, "ModerationTime/D");
T2->Branch( "ThermalTime", &ThermalTime, "ThermalTime/D");
T2->Branch( "CaptureTime", &CaptureTime, "CaptureTime/D");
T2->Branch( "ThermalRadius", &ThermalRadius, "ThermalRadius/D");
T2->Branch( "ModDistance", &ModDistance, "ModDistance/D");
T2->Branch( "CapDistance", &CaptureDistance, "CaptureDistance/D");
T2->Branch( "NScatter", &NScatter, "NScatter/I");
// T2->Branch( "LiTotCross", &LiTotCross, "LiTotCross/I");
// T2->Branch( "LiFastCross", &LiFastCross, "LiFastCross/I");
// T2->Branch( "LiThermalCross", &LiThermalCross, "LiThermalCross/I");
// T2->Branch( "FurthestDistance",&Furthest, "FurthestDistance/D");
// T2->Branch( "HDPEReflection", &NeutronReflect, "HDPEReflection/I");
// T2->Branch( "NeutronEscape", &NeutronEscape, "NeutronEscape/I");
T2->Branch( "CaptureVol", &CaptureVol, "CaptureVol/I");
T2->Branch( "CaptureCube", &CaptureCube, "CaptureCube/I");
T2->Branch( "NeutronDeltaR",&NeutronDeltaR, "NeutronDeltaR/D");
T2->Branch( "PositronEndX", &PositronEndX, "PositronEndX/D");
T2->Branch( "PositronEndY", &PositronEndY, "PositronEndY/D");
T2->Branch( "PositronEndZ", &PositronEndZ, "PositronEndZ/D");
T2->Branch( "PositronDistance",&PositronDistance, "PositronDistance/D");
T2->Branch( "AnnihilationTime",&AnnihilationTime, "AnnihilationTime/D");
T2->Branch( "AnnihilationVol", &AnnihilationVol, "AnnihilationVol/I");
T2->Branch( "AnnihilationCube",&AnnihilationCube, "AnnihilationCube/I");
T2->Branch( "PositronDeltaR", &PositronDeltaR, "PositronDeltaR/D");
T2->Branch( "ibdDeltaT",&ibdDeltaT, "ibdDeltaT/D");
T2->Branch( "ibdDeltaR",&ibdDeltaR, "ibdDeltaR/D");
}
T5 = new TTree("T5","Cube-level information from sensitive detectors");
//T5->Branch( "cubeID", &cubenumber);
T5->Branch( "pdg", &pdg);
T5->Branch( "trackid", &partID);
T5->Branch( "parentid", &mothID);
//cube position, format: zzxxyy
T5->Branch( "volid", &cubenumber);
T5->Branch( "edep_tot", &edepTot);
T5->Branch( "edep_pvt", &edepPvt);
T5->Branch( "edep_zns", &edepLi);
// T5->Branch( "time", &timeOfBirth);
T5->Branch( "time_last", &timeOfDemise);
T5->Branch( "tracklength", &distanceTravelled);
T5->Branch( "EdepPosX", &EdepPosX);
T5->Branch( "EdepPosY", &EdepPosY);
T5->Branch( "EdepPosZ", &EdepPosZ);
if( m_storeT6 )
{
T6 = new TTree("T6","particles tracking info");
T6->Branch( "pdg", &partPDG);
T6->Branch( "trackid", &partTrackID);
T6->Branch( "parentid", &partParentID);
T6->Branch( "energy", &partEnergy);
T6->Branch( "deltaE", &partDeltaE);
T6->Branch( "Edepos", &partDeposited);
T6->Branch( "posx", &partX);
T6->Branch( "posy", &partY);
T6->Branch( "posz", &partZ);
T6->Branch( "time", &partTime);
T6->Branch( "volid", &partVolID);
T6->Branch( "cubeid", &partCubeID);
//debug
// T6->Branch( "process", &partProcess);
}
if( m_storeT7 )
{
T7 = new TTree("T7","particles entering the detector volume");
T7->Branch( "pdg", &enteringPartID);
T7->Branch( "trackid", &enteringTrackID);
T7->Branch( "parentid", &enteringParentID);
T7->Branch( "energy", &enteringEnergy);
// T7->Branch( "deltaE", &enteringDeposited);
T7->Branch( "posx", &enteringX);
T7->Branch( "posy", &enteringY);
T7->Branch( "posz", &enteringZ);
T7->Branch( "momx", &enteringXmom);
T7->Branch( "momy", &enteringYmom);
T7->Branch( "momz", &enteringZmom);
T7->Branch( "time", &enteringTime);
T7->Branch( "prevolid", &enteringPreVolID);
T7->Branch( "postvolid", &enteringPostVolID);
}
if( m_storeT8 )
{
T8 = new TTree("T8","spallation neutrons from muons");
T8->Branch( "trackid", &nGenerationTrackID);
T8->Branch( "parentid", &nGenerationParentID);
T8->Branch( "energy", &nGenerationEnergy);
T8->Branch( "posx", &nGenerationPosX);
T8->Branch( "posy", &nGenerationPosY);
T8->Branch( "posz", &nGenerationPosZ);
T8->Branch( "time", &nGenerationTime);
T8->Branch( "cubeid", &nGenerationCube);
T8->Branch( "volid", &nGenerationVol);
T8->Branch( "process", &nGenerationProcess);
}
if( m_storeT9 ) {
T9 = new TTree("T9","statistical process info from MC");
T9->Branch( "pdg", &summaryPartPDG);
T9->Branch( "trackid", &summaryPartID);
T9->Branch( "parentid", &summaryParentID);
T9->Branch( "parentPDG", &summaryParentPDG);
T9->Branch( "initialEkin", &summaryInitialEkin);
T9->Branch( "nSteps", &summaryNsteps);
T9->Branch( "nCubes", &summaryNcubes);
T9->Branch( "firstVolID", &summaryFirstVol);
T9->Branch( "lastVolID", &summaryLastVol);
T9->Branch( "CreatorProcessID", &summaryCreatorProcessID);
T9->Branch( "LastProcessID", &summaryLastProcessID);
T9->Branch( "posx", &summaryPosX);
T9->Branch( "posy", &summaryPosY);
T9->Branch( "posz", &summaryPosZ);
T9->Branch( "PxInPVT", &summaryMomX);
T9->Branch( "PyInPVT", &summaryMomY);
T9->Branch( "PzInPVT", &summaryMomZ);
/* -----
the branches below are kept here commented, while doing a validation of the output.
They can be restored if needed.
*/
// T9->Branch( "nIonis", &summaryNIonisation);
// T9->Branch( "nCompton", &summaryNCompton);
// T9->Branch( "nBrem", &summaryNBrem);
// T9->Branch( "CreatorProcess", &summaryCreatorProcess);
// T9->Branch( "LastProcess", &summaryLastProcess);
}
G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
// G4UImanager* UI = G4UImanager::GetUIpointer();
// UI->ApplyCommand("/tracking/storeTrajectory 1");
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
void SolidRunAction::EndOfRunAction(const G4Run*)
{
// Go to output file
outp->cd();
Stats.Evt++;
T0->Fill();
T0->Write("",TObject::kOverwrite);
T1->Write("",TObject::kOverwrite);
if( m_storeT2 )
T2->Write("",TObject::kOverwrite);
// T3->Write("",TObject::kOverwrite);
// T4->Write("",TObject::kOverwrite);
T5->Write("",TObject::kOverwrite);
if( m_storeT6 )
T6->Write("",TObject::kOverwrite);
if( m_storeT7 )
T7->Write("",TObject::kOverwrite);
if( m_storeT8 )
T8->Write("",TObject::kOverwrite);
if( m_storeT9 )
T9->Write("",TObject::kOverwrite);
outp->Close();
//display run time and write to file Runtime.out
time_t end = time(NULL);
G4int elapsed = end - start;
G4cout<<"Run Completed in "<<elapsed/3600<<":"<<(elapsed%3600)/60<<":"<<((elapsed%3600)%60)<<G4endl;
//#ifdef G4VIS_USE
// if (G4VVisManager::GetConcreteInstance()) {
// G4VisManager *visMan = (G4VisManager *) G4VVisManager::GetConcreteInstance();
// // Apply command only when current viewer is available
// if(visMan->GetCurrentViewer()!=0) {
// G4UImanager::GetUIpointer()->ApplyCommand("/vis/viewer/update");
// }
// }
//#endif
}
void SolidRunAction::UpdateTracking(void){
// T3->Fill();
}
void SolidRunAction::UpdateStatistics(RunTally aRunTally){
Stats = aRunTally;
lightOutput = 1;
if(m_storeLightT && edepTot.size()<=0 ) lightOutput = 0;
if (lightOutput) {
T1->Fill();
if( m_storeT2)
T2->Fill();
// T3->Fill();
// T4->Fill();
T5->Fill();
if( m_storeT6)
T6->Fill();
if( m_storeT7)
T7->Fill();
if( m_storeT8)
T8->Fill();
if( m_storeT9)
T9->Fill();
}
// printf("Moderation Length = %e mm, \t Moderation Time = %1.2f ns , \n", ModDistance, ModTime);
// printf("Thermal Radius = %e mm, \t Thermal Neutron Time = %1.2f ns ,\n", ThermalRadius, ThermalTime);
// printf("Capture Distance = %e, \t Capture Time = %1.2f ns \n", CaptureDistance, CaptureTime);
// printf("Crosses before moderated = %d, \t \t Thermal Crosses = %d, \t Total = %d \n", LiFastCross, LiThermalCross, LiTotCross);
}
void SolidRunAction::SaveFile(G4int){
Stats.Evt++;
T0->Fill();
T0->Write("",TObject::kOverwrite);
T1->Write("",TObject::kOverwrite);
if( m_storeT2 )
T2->Write("",TObject::kOverwrite);
// T3->Write("",TObject::kOverwrite);
// T4->Write("",TObject::kOverwrite);
T5->Write("",TObject::kOverwrite);
if( m_storeT6 )
T6->Write("",TObject::kOverwrite);
if( m_storeT7 )
T7->Write("",TObject::kOverwrite);
if( m_storeT8 )
T8->Write("",TObject::kOverwrite);
if( m_storeT9 )
T9->Write("",TObject::kOverwrite);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....