Hi all,
I am using Ntuples to store events detected by a pixelated CZT detector and events which escape the world volume. For this I am creating two Ntuples where the ID 0 Ntuple stores the CZT events while ID 1 Ntuple stores the events escaping the world volume. They are stored in a .root
file which at the EndOfRunAction
is read using G4AnalysisReader and the info is converted to a .fits
file.
When I run this with single energy it works fine, but when I use \control\foreach
to specify multiple energy runs, the analysis manager creates Ntuples with higher IDs than just 0 and 1. Because I am using a Ntuple ID to store the info to fits
files only the events corresponding to the first energy are stored.
So my question is how can I reset the Ntuple IDs at the start of every run?
Thanks in advance,
Sujay
Here is the relevant part of my RunAction
:
void DetBoxRunAction::BeginOfRunAction(const G4Run* aRun)
{
// Total events
fnofEvents = aRun->GetNumberOfEventToBeProcessed();
// Set time dependent seed
fseed = time(NULL);
CLHEP::HepRandom::setTheSeed(fseed);
// Get analysis manager
auto analysisManager = G4AnalysisManager::Instance();
G4cout << "Using " << analysisManager->GetType() << G4endl;
// Set verbose and merging
analysisManager->SetVerboseLevel(1);
//analysisManager->SetNtupleMerging(true);
// Open file to save ntuples
fRootfile = "CZTRunData_" + std::to_string(fRunCnt);
analysisManager->OpenFile(fRootfile);
// Creating ntuples
// Ntuple ID 0
analysisManager->CreateNtuple("CZTEvents", "Event data");
analysisManager->CreateNtupleIColumn("EventID");
analysisManager->CreateNtupleIColumn("EvtOrderID");
analysisManager->CreateNtupleIColumn("PixIDX");
analysisManager->CreateNtupleDColumn("Energy");
analysisManager->CreateNtupleIColumn("Mult");
analysisManager->FinishNtuple();
// Ntuple ID 1
analysisManager->CreateNtuple("OutOfWorldEvents", "Event pos and ene");
analysisManager->CreateNtupleIColumn("EventID");
analysisManager->CreateNtupleDColumn("Theta");
analysisManager->CreateNtupleDColumn("Phi");
analysisManager->CreateNtupleDColumn("Energy");
analysisManager->FinishNtuple();
}
//***********************************************/
void DetBoxRunAction::EndOfRunAction(const G4Run* aRun)
{
// get primary energy
fInpEne = fPrimaryGeneratorAction->fgps->GetParticleEnergy() / keV;
// Define detected events file name
std::stringstream evtfitsfname;
evtfitsfname << "CZTEvents_" << std::setw(4) << std::setfill('0') << std::setprecision(0)
<< std::fixed << fInpEne << "_keV";
fEvtFitsfile = evtfitsfname.str();
// Define OutOfWorld events file name
std::stringstream OWfitsfname;
OWfitsfname << "OutOfWorldEvents_" << std::setw(4) << std::setfill('0') << std::setprecision(0)
<< std::fixed << fInpEne << "_keV";
fOWFitsfile = OWfitsfname.str();
if (isMaster)
{
auto analysisManager = G4AnalysisManager::Instance();
// save histograms & ntuple
analysisManager->Write();
analysisManager->CloseFile();
G4String rootfname = analysisManager->GetFileName() + ".root";
readroot(rootfname); // Read the root file using G4AnalysisReader (function given below)
writeEvtfits(fEvtFitsfile); // write CZT events fits file
writeOWfits(fOWFitsfile); // write Out of world events fits file
remove(rootfname);
fRunCnt++;
}
}
//***********************************************/
void DetBoxRunAction::readroot(G4String rootfname)
{
// Code to create (or get) analysis reader
auto analysisReader = G4AnalysisReader::Instance();
// Read Events ntuple
analysisReader->SetFileName(rootfname);
G4int ntupleId = analysisReader->GetNtuple("CZTEvents");
if ( ntupleId == 0 )
{
G4int EventID, EvtOrder, boxID, moduleID, pixID, mult, pixIDX;
G4double energy;
analysisReader->SetNtupleIColumn(0, "EventID", EventID);
analysisReader->SetNtupleIColumn(0, "EvtOrderID", EvtOrder);
analysisReader->SetNtupleIColumn(0, "PixIDX", pixIDX);
analysisReader->SetNtupleDColumn(0, "Energy", energy);
analysisReader->SetNtupleIColumn(0, "Mult", mult);
while ( analysisReader->GetNtupleRow(ntupleId) )
{
boxID = pixIDX / 5120;
moduleID = (pixIDX % 5120) / 256;
pixID = (pixIDX % 5120) % 256;
// populate vectors to pass to fits.
feventID.push_back(EventID);
fevtorderID.push_back(EvtOrder);
fpixIDX.push_back(pixIDX);
fboxID.push_back(boxID);
fmoduleID.push_back(moduleID);
fpixID.push_back(pixID);
fEvtenergy.push_back(energy);
fmult.push_back(mult);
}
}
// Read Events ntuple
ntupleId = analysisReader->GetNtuple("OutOfWorldEvents");
if (ntupleId == 1)
{
G4int OWEventID;
G4double OWtheta, OWphi, OWenergy;
analysisReader->SetNtupleIColumn(1, "EventID", OWEventID);
analysisReader->SetNtupleDColumn(1, "Theta", OWtheta);
analysisReader->SetNtupleDColumn(1, "Phi", OWphi);
analysisReader->SetNtupleDColumn(1, "Energy", OWenergy);
while ( analysisReader->GetNtupleRow(ntupleId) )
{
fOWeventID.push_back(OWEventID);
fOWtheta.push_back(OWtheta);
fOWphi.push_back(OWphi);
fOWenergy.push_back(OWenergy);
}
}
}