Introduce examine_material.C

Macro which compares the names of materials as generated from a simulation
to those obtained directly from opening the gemometry binaries directory.
User specifies the setup he wishes to examine. Readme includes suggested
commands.

macro/geometry/README.md

@@ -203,3 +203,21 @@ Allows to run the `check_radlen` script on sets of geometries present in the sub
   Example for a version: `./macro/geometry/check_radlen_bulk.sh v22a` \
   Special tags which cannot be used in this way: `git`, `all`, `main`
 
+
+---
+
+## examine_materials.C
+
+Root macro which compares the materials in the simulated TGeoManager to the materials as extracted from the geometry
+binaries for a specified setup. Compares only the names.
+
+Sugggested commands
+
+root -q examine_materials.C
+
+root -q 'examine_materials.C("sis100_muon_jpsi_DEC21")'
+
+ for file in `ls ../../geometry/setup/setup_mcbm_beam_2022_*`; do tag=`echo $file | awk -F'/' '{print $5}' | sed 's|setup_|"|'  | sed 's_.C$_"_'`; echo $tag; root -q './examine_material.C('$tag')'; done > MCBM_2022_EXAMINE
+
+
+

macro/geometry/examine_material.C

+/* Copyright (C) 2023 Facility for Antiproton and Ion Research in Europe, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Eoin Clerkin [committer] */
+
+#define MAX_MATERIALS 200
+
+#include <cstdlib>
+#include <fstream>
+#include <iostream>
+
+#include <libgen.h>
+
+#define NODE_DEPTH 25
+
+/*  The materials in the transported TGeoManager will be compared to the materials in the subsystem binaries.
+    The general idea pre-writing will be to run transport for a user specified setup generating the TGeoManager. 
+    The number of volumes for each material will be recorded in a standard map. Each subsystem will be extracted
+    from the CBMSETUP, and the materials will be extraced from the respective binary and then substracted from 
+    the from the TGeoManager version original. After going through all subsystems, if there is no volumes in the
+    original standard, the two TGeoManager and the material in the subsystem match, otherwise there is a material
+    mismatch detected. Output to the screen will help diagonsis. 					       */
+
+char* padd(char* source_string, unsigned int size)
+{
+
+  char* padded_string = (char*) malloc(sizeof(char) * (size + 1));
+  strncpy(padded_string, source_string, size + 1);
+
+  unsigned int i = strlen(source_string);
+  if (i < size)
+    for (; i < size; i++)
+      padded_string[i] = ' ';
+  padded_string[size] = '\0';
+
+  return padded_string;
+};
+
+typedef struct MAT_VOL mat_vol_t;
+
+struct MAT_VOL {
+  int count;  // Number of volumes
+  char name[100];
+};
+
+int REG_MATERIALS = 0;  // counter for the currently number of registered materials.
+
+int get_register_count_material(const char* MATERIAL_NAME, mat_vol_t* material_array)
+{
+  int i = 0;
+
+  while (i <= MAX_MATERIALS && i <= REG_MATERIALS) {
+    // Register
+    if (REG_MATERIALS == i) {
+      strcpy((material_array[i]).name, MATERIAL_NAME);
+      REG_MATERIALS++;
+      material_array[i].count = 0;
+      break;
+    };
+    // Found
+    if (!strcmp((material_array[i]).name, MATERIAL_NAME)) { break; };
+    i++;
+  };
+  material_array[i].count++;
+
+  return 0;
+}
+
+int extract_mat_vol(char* fileName, mat_vol_t* material_array)
+{
+  TGeoNode* node[NODE_DEPTH];
+  TGeoMedium* med;
+  TGeoManager* gman = nullptr;
+  TGeoVolume* top   = nullptr;
+
+  TFile file(fileName, "READ");
+  if (!file.IsOpen()) {
+    std::cerr << "Was a root geo file supplied?" << endl;
+    exit(1);
+  };
+
+  file.GetListOfKeys()->Print();
+  file.GetSize();
+  char geo_tag[100];
+  strcpy(geo_tag, fileName);
+  char* ptr = strstr(geo_tag, ".geo.root");
+  *(ptr)    = '\0';
+
+  gman = (TGeoManager*) file.Get("FAIRGeom");
+  if (gman != nullptr) { top = gman->GetTopVolume(); };
+  if (top == nullptr) top = (TGeoVolume*) file.Get(geo_tag);
+  if (top == nullptr) top = (TGeoVolume*) file.Get("top");
+  if (top == nullptr) top = (TGeoVolume*) file.Get("TOP");
+  if (top == nullptr) top = (TGeoVolume*) file.Get("Top");
+  if (top == nullptr) top = (TGeoVolume*) file.Get("geometryFromGDML");
+  if (top == nullptr) top = (TGeoVolume*) file.Get(file.GetListOfKeys()->First()->GetName());
+  if (top == nullptr) {
+    std::cerr << "No Top Volume found. Is the TGeoManager or Top Volume unusally named" << std::endl;
+    exit(1);
+  };
+
+  TObjArray* nodes = top->GetNodes();
+
+  int i_array[NODE_DEPTH], num_array[NODE_DEPTH], j;
+  for (int i = 0; i < NODE_DEPTH; i++)
+    i_array[i] = 0;
+  for (int i = 0; i < NODE_DEPTH; i++)
+    num_array[i] = 0;
+
+  j            = 0;
+  i_array[0]   = 0;
+  num_array[0] = nodes->GetEntries();
+
+  while (num_array[0] != 0) {
+    if (j == 0) { node[0] = static_cast<TGeoNode*>(nodes->At(i_array[0])); }
+    else {
+      node[j] = static_cast<TGeoNode*>(node[j - 1]->GetDaughter(i_array[j]));
+    };
+    i_array[j]++;  // Update number.
+    med              = node[j]->GetMedium();
+    num_array[j + 1] = node[j]->GetNdaughters();
+
+    get_register_count_material(med->GetMaterial()->GetName(), material_array);
+
+    if (num_array[j + 1] > 0) {
+      j++;
+      num_array[j + 2] = 0;
+    };
+
+    while (i_array[j] == num_array[j]) {
+      num_array[j] = 0;
+      i_array[j]   = 0;  // Reset the counter before falling back.
+      j--;
+    };
+
+    if (j >= NODE_DEPTH) {
+      std::cerr << "Maximum depth of " << NODE_DEPTH << " exceeded. Increase NODE_DEPTH in header of macro," << endl;
+      exit(NODE_DEPTH);
+    };
+  };
+
+  file.Close();
+  return 0;
+}
+
+int examine_material(const char* setup = "sis100_electron")
+{
+  TString srcDir    = gSystem->Getenv("VMCWORKDIR");
+  ECbmEngine engine = kGeant3;
+
+  CbmTransport run;
+  run.SetEngine(engine);
+  run.AddInput(srcDir + "/input/urqmd.auau.10gev.centr.root", kUnigen);
+  run.SetOutFileName("examine.tra.root", kTRUE);
+  run.SetParFileName("examaine.par.root");
+  run.SetGeoFileName("examine.geo.root");
+  run.LoadSetup(setup);
+  run.SetBeamPosition(0, 0, 0.1, 0.1);
+  run.SetRandomSeed(1234);  // I fix this number so everything is deterministic
+  run.Run(1);
+
+  mat_vol_t* transport_materials = (mat_vol_t*) malloc(sizeof(mat_vol_t) * MAX_MATERIALS);
+  for (int i = 0; i < MAX_MATERIALS; i++) {
+    transport_materials[i].count = 0;
+    strcpy(transport_materials[i].name, "");
+  };
+
+  extract_mat_vol("examine.geo.root", transport_materials);
+
+  // Detector and passive subsystems
+  for (int i = 0; i < 120; i++)
+    printf("#");
+  printf("\n");
+
+  TString fileName;
+  char fileNameC[100];
+  CbmSetup* cs                  = run.GetSetup();
+  mat_vol_t* detector_materials = (mat_vol_t*) malloc(sizeof(mat_vol_t) * 500);
+  for (int i = 0; i < MAX_MATERIALS; i++) {
+    detector_materials[i].count = 0;
+    strcpy(detector_materials[i].name, "");
+  };
+
+  REG_MATERIALS = 0;
+  for (ECbmModuleId subsystem = ECbmModuleId::kMvd; subsystem < ECbmModuleId::kCave; ++subsystem) {
+    if (cs->IsActive(subsystem)) {
+
+      cs->GetGeoFileName(subsystem, fileName);  // returns true
+
+
+      // Handling double beampipe
+      if ((subsystem == ECbmModuleId::kPipe) && strchr(fileName.Data(), ':')) {
+
+        fileName = srcDir + "/geometry/" + fileName;  // Get the right directory
+        printf("fileName is %s\n", fileName.Data());
+        strcpy(fileNameC, fileName.Data());
+        char* change = strchr(fileNameC, ':');
+        *change      = '\0';
+        printf("FIRST PIPE fileName now is %s\n", fileNameC);
+        extract_mat_vol(fileNameC, detector_materials);
+        strcpy(fileNameC, srcDir + "/geometry/" + (change + 1));
+        printf("SECOND PIPE fileName now is %s\n", fileNameC);
+        extract_mat_vol(fileNameC, detector_materials);
+
+        // Other subystems
+      }
+      else {
+
+        fileName = srcDir + "/geometry/" + fileName;  // Get the right directory
+        printf("fileName is %s\n", fileName.Data());
+        strcpy(fileNameC, fileName.Data());
+        extract_mat_vol(fileNameC, detector_materials);
+      };
+    };
+  }
+
+  // Quick Sort Technique
+  mat_vol_t temp_mat;
+  for (int i = 0; i < (MAX_MATERIALS - 1); ++i) {
+    for (int j = 0; j < (MAX_MATERIALS - 1 - i); ++j) {
+      if ((detector_materials[j].count) < (detector_materials[j + 1].count)) {
+        temp_mat                  = detector_materials[j + 1];
+        detector_materials[j + 1] = detector_materials[j];
+        detector_materials[j]     = temp_mat;
+      };
+    };
+  };
+  for (int i = 0; i < (MAX_MATERIALS - 1); ++i) {
+    for (int j = 0; j < (MAX_MATERIALS - 1 - i); ++j) {
+      if ((transport_materials[j].count) < (transport_materials[j + 1].count)) {
+        temp_mat                   = transport_materials[j + 1];
+        transport_materials[j + 1] = transport_materials[j];
+        transport_materials[j]     = temp_mat;
+      };
+    };
+  };
+  int SUMT = 0, SUMD = 0;
+  for (int j = 0; j < MAX_MATERIALS; j++) {
+    if ((transport_materials + j)->count > 0) SUMT += (transport_materials + j)->count;
+    if ((detector_materials + j)->count > 0) SUMD += (detector_materials + j)->count;
+  };
+  (void) printf("TOTAL NUMBER OF VOLUMES\tTRANSPORTED: %d\t DETECTORS: %d\n", SUMT, SUMD);
+
+  printf("%s\t\tTransported\t\tSubsystems\t\tDifference \n", padd("Material", 25));
+  int FAIL = 0;  // Success
+  int DIFF = 0;
+
+  for (int i = 0; i < MAX_MATERIALS; i++) {
+    if (0 < (transport_materials + i)->count) {
+      printf("%s\t\t%d", padd((transport_materials + i)->name, 25), (transport_materials + i)->count);
+      for (int j = 0; j < MAX_MATERIALS; j++) {
+        if (!strcmp((transport_materials[i].name), (detector_materials[j].name)))
+          printf("\t\t%d", (detector_materials + i)->count);
+      };
+      DIFF = ((transport_materials + i)->count) - ((detector_materials + i)->count);
+      printf("\t\t%d\n", DIFF);
+    };
+    if ((strcmp(transport_materials[i].name, "dummy") != 0) && (strcmp(transport_materials[i].name, "air") != 0))
+      FAIL += abs(DIFF);
+  };
+
+  if (0 == FAIL) { (void) printf("[SUCCESS] Excluding air and dummy materials, no material errors were detected. \n"); }
+  else {
+    (void) printf("[FAILED] Excluding air and dummy materials, %d issues with materials were detected. \n", FAIL);
+  };
+
+  return FAIL;  // Success (when FAIL=0)
+}