diff --git a/inc/TRestGeant4Metadata.h b/inc/TRestGeant4Metadata.h index fcbe445..e962ece 100644 --- a/inc/TRestGeant4Metadata.h +++ b/inc/TRestGeant4Metadata.h @@ -137,6 +137,14 @@ class TRestGeant4Metadata : public TRestMetadata { /// If it is zero its value will be assigned using the system timestamp. Long_t fSeed = 0; + /// \brief If the simulation is produced as a merge of multiple files, this stored the seeds of the + /// individual files + std::vector fMergeSeeds; + + /// \brief If the simulation is produced as a merge of multiple files, this stored the number primaries + /// for the individual files + std::vector fMergeNEvents; + /// \brief If this parameter is set to 'true' it will save all events even if they leave no energy in the /// sensitive volume (used for debugging purposes). It is set to 'false' by default. Bool_t fSaveAllEvents = false; @@ -379,11 +387,9 @@ class TRestGeant4Metadata : public TRestMetadata { void SetActiveVolume(const TString& name, Double_t chance, Double_t maxStep = 0); - void SetSimulationTime(Long64_t time) { fSimulationTime = time; } - void PrintMetadata() override; - void Merge(const TRestGeant4Metadata&); + void Merge(const TRestMetadata&) override; TRestGeant4Metadata(); TRestGeant4Metadata(const char* configFilename, const std::string& name = ""); @@ -393,7 +399,7 @@ class TRestGeant4Metadata : public TRestMetadata { TRestGeant4Metadata(const TRestGeant4Metadata& metadata); TRestGeant4Metadata& operator=(const TRestGeant4Metadata& metadata); - ClassDefOverride(TRestGeant4Metadata, 15); + ClassDefOverride(TRestGeant4Metadata, 16); // Allow modification of otherwise inaccessible / immutable members that shouldn't be modified by the user friend class SteppingAction; diff --git a/macros/REST_Geant4_MergeRestG4Files.C b/macros/REST_Geant4_MergeRestG4Files.C index 3cf1b77..e52a524 100644 --- a/macros/REST_Geant4_MergeRestG4Files.C +++ b/macros/REST_Geant4_MergeRestG4Files.C @@ -1,3 +1,5 @@ +#include + #include #include @@ -59,13 +61,26 @@ void REST_Geant4_MergeRestG4Files(const char* outputFilename, const char* inputF TRestGeant4Event* mergeEvent = nullptr; auto mergeEventTree = mergeRun->GetEventTree(); + auto mergeAnalysisTree = mergeRun->GetAnalysisTree(); mergeEventTree->Branch("TRestGeant4EventBranch", "TRestGeant4Event", &mergeEvent); + const string outputTempFilename = + string(outputFilename).substr(0, string(outputFilename).size() - 5) + ".temp.root"; + { + TFileMerger merger; + merger.OutputFile(outputTempFilename.c_str()); + merger.AddObjectNames("EventTree AnalysisTree"); + for (unsigned int i = 0; i < inputFiles.size(); i++) { + merger.AddFile(inputFiles[i].c_str()); + } + merger.PartialMerge(TFileMerger::kAll | TFileMerger::kIncremental | TFileMerger::kOnlyListed); + } + set eventIds; // std::set is sorted from lower to higher automatically long long eventCounter = 0; // iterate over all other files - for (int i = 0; i < inputFiles.size(); i++) { + for (unsigned int i = 0; i < inputFiles.size(); i++) { cout << "Processing file " << i + 1 << "/" << inputFiles.size() << endl; map @@ -81,7 +96,10 @@ void REST_Geant4_MergeRestG4Files(const char* outputFilename, const char* inputF TRestGeant4Event* event = nullptr; auto eventTree = run.GetEventTree(); eventTree->SetBranchAddress("TRestGeant4EventBranch", &event); - for (int j = 0; j < eventTree->GetEntries(); j++) { + for (unsigned int j = 0; j < eventTree->GetEntries(); j++) { + eventCounter++; + continue; // not working at this time (this logic works, but we are using TFileMerger which does + // not work with this) eventTree->GetEntry(j); *mergeEvent = *event; @@ -106,33 +124,76 @@ void REST_Geant4_MergeRestG4Files(const char* outputFilename, const char* inputF eventIds.insert(mergeEvent->GetID()); mergeEventTree->Fill(); - mergeRun->GetAnalysisTree()->Fill(); + // TODO: this just adds empty entries to the analysis tree. It should be merged + mergeAnalysisTree->Fill(); + } - eventCounter++; + // add the remaining metadata of the first file to the merged file + if (i == 0) { + // iterate over all keys of the file to find inheritance of TRestMetadata + TIter nextkey(run.GetInputFile()->GetListOfKeys()); + TKey* key; + while ((key = (TKey*)nextkey())) { + const auto obj = key->ReadObj(); + if (obj->InheritsFrom("TRestMetadata")) { + if (obj->InheritsFrom("TRestGeant4Metadata")) { + // This is merged and added later + continue; + } + const auto metadataKey = (TRestMetadata*)obj; + mergeRun->GetOutputFile()->cd(); + metadataKey->Write(); + } + } } } - cout << "Output filename: " << mergeRun->GetOutputFileName() << endl; - cout << "Output file: " << mergeRun->GetOutputFile() << endl; - mergeRun->GetOutputFile()->cd(); - gGeoManager->Write("Geometry", TObject::kOverwrite); + if (gGeoManager != nullptr) { + gGeoManager->Write("Geometry", TObject::kOverwrite); + } - mergeMetadata.SetName("geant4Metadata"); mergeMetadata.Write(); + mergeRun->UpdateOutputFile(); mergeRun->CloseFile(); - // Open the file again to check the number of events + // At this point we have two files: the "mergeRun" file with the updated metadata and the "temp" file with + // the event and analysis tree (the event tree here does not have updated event ids) Open the file again + + { + auto fileWithMetadata = TFile::Open(mergeRun->GetOutputFileName()); + auto fileWithTrees = TFile::Open(outputTempFilename.c_str(), "UPDATE"); + + // copy all objects from the file with metadata except "EventTree" and "AnalysisTree" + TIter nextkey(fileWithMetadata->GetListOfKeys()); + TKey* key; + while ((key = (TKey*)nextkey())) { + const auto obj = key->ReadObj(); + if (obj->InheritsFrom("TTree")) { + continue; + } + fileWithTrees->cd(); + obj->Write(); + } + + // close files + fileWithMetadata->Close(); + fileWithTrees->Close(); + } + + // replace the "run" file by the temp file + remove(mergeRun->GetOutputFileName()); + rename(outputTempFilename.c_str(), mergeRun->GetOutputFileName()); + + // to check the number of events TRestRun runCheck(outputFilename); if (runCheck.GetEntries() != eventCounter) { cerr << "ERROR: number of events in the output file (" << runCheck.GetEntries() << ") does not match the number of events in the input files (" << eventCounter << ")" << endl; exit(1); } - cout << "Number of events in the output file: " << runCheck.GetEntries() << " matches internal count" - << endl; } #endif diff --git a/src/TRestGeant4Metadata.cxx b/src/TRestGeant4Metadata.cxx index 184bb3b..6125e76 100644 --- a/src/TRestGeant4Metadata.cxx +++ b/src/TRestGeant4Metadata.cxx @@ -1433,7 +1433,15 @@ void TRestGeant4Metadata::PrintMetadata() { TRestMetadata::PrintMetadata(); RESTMetadata << "Geant4 version: " << GetGeant4Version() << RESTendl; - RESTMetadata << "Random seed: " << GetSeed() << RESTendl; + if (!fIsMerge) { + RESTMetadata << "Random seed: " << GetSeed() << RESTendl; + } else { + RESTMetadata << "This Geant4 simulation is a merge of " << fMergeSeeds.size() << " files" << RESTendl; + for (unsigned int i = 0; i < fMergeSeeds.size(); i++) { + RESTMetadata << " - File " << i << " Random seed: " << fMergeSeeds[i] + << " - Number of primaries: " << fMergeNEvents[i] << RESTendl; + } + } RESTMetadata << "GDML geometry: " << GetGdmlReference() << RESTendl; RESTMetadata << "GDML materials reference: " << GetMaterialsReference() << RESTendl; RESTMetadata << "Sub-event time delay: " << GetSubEventTimeDelay() << " us" << RESTendl; @@ -1496,7 +1504,9 @@ void TRestGeant4Metadata::PrintMetadata() { Int_t TRestGeant4Metadata::GetActiveVolumeID(const TString& name) { Int_t id; for (id = 0; id < (Int_t)fActiveVolumes.size(); id++) { - if (fActiveVolumes[id] == name) return id; + if (fActiveVolumes[id] == name) { + return id; + } } return -1; } @@ -1547,7 +1557,9 @@ Bool_t TRestGeant4Metadata::isVolumeStored(const TString& volume) const { Double_t TRestGeant4Metadata::GetStorageChance(const TString& volume) { Int_t id; for (id = 0; id < (Int_t)fActiveVolumes.size(); id++) { - if (fActiveVolumes[id] == volume) return fChance[id]; + if (fActiveVolumes[id] == volume) { + return fChance[id]; + } } RESTWarning << "TRestGeant4Metadata::GetStorageChance. Volume " << volume << " not found" << RESTendl; @@ -1580,19 +1592,23 @@ size_t TRestGeant4Metadata::GetGeant4VersionMajor() const { return std::stoi(majorVersion.Data()); } -void TRestGeant4Metadata::Merge(const TRestGeant4Metadata& metadata) { +void TRestGeant4Metadata::Merge(const TRestMetadata& metadata) { + const auto restMetadata = *dynamic_cast(&metadata); fIsMerge = true; fSeed = 0; // seed makes no sense in a merged file + fMergeSeeds.push_back(restMetadata.fSeed); + fMergeNEvents.push_back(restMetadata.fNEvents); - fNEvents += metadata.fNEvents; - fNRequestedEntries += metadata.fNRequestedEntries; - fSimulationTime += metadata.fSimulationTime; + fNEvents += restMetadata.fNEvents; + fNRequestedEntries += restMetadata.fNRequestedEntries; + fSimulationTime += restMetadata.fSimulationTime; } TRestGeant4Metadata::TRestGeant4Metadata(const TRestGeant4Metadata& metadata) { *this = metadata; } TRestGeant4Metadata& TRestGeant4Metadata::operator=(const TRestGeant4Metadata& metadata) { fIsMerge = metadata.fIsMerge; + fName = metadata.fName; fGeant4GeometryInfo = metadata.fGeant4GeometryInfo; fGeant4PhysicsInfo = metadata.fGeant4PhysicsInfo; fGeant4PrimaryGeneratorInfo = metadata.fGeant4PrimaryGeneratorInfo; @@ -1621,5 +1637,7 @@ TRestGeant4Metadata& TRestGeant4Metadata::operator=(const TRestGeant4Metadata& m fKillVolumes = metadata.fKillVolumes; fRegisterEmptyTracks = metadata.fRegisterEmptyTracks; fMagneticField = metadata.fMagneticField; + fMergeNEvents = metadata.fMergeNEvents; + fMergeSeeds = metadata.fMergeSeeds; return *this; }