Merge supports analysis files

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<Long_t> fMergeSeeds;
+
+    /// \brief If the simulation is produced as a merge of multiple files, this stored the number primaries
+    /// for the individual files
+    std::vector<Long_t> 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;

macros/REST_Geant4_MergeRestG4Files.C

@@ -1,3 +1,5 @@
+#include <TFileMerger.h>
+
 #include <string>
 #include <vector>
 
@@ -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<Int_t> 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<Int_t, Int_t>
@@ -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

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<const TRestGeant4Metadata*>(&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;
 }