Merge pull request #15 from giovannimarchiori/gmarchio-main-20231103-…

…k4geomigration-FCCDetPR63

Migration to k4geo

.gitignore

@@ -139,5 +139,5 @@ dmypy.json
 
 # my stuff 
 *.sub
-*.root
+
 

FCCSW_ecal/neighbours_theta.py

@@ -1,44 +1,41 @@
+from Configurables import ApplicationMgr
+from Configurables import CreateFCCeeCaloNeighbours
 import os
 from Gaudi.Configuration import *
 
 # Detector geometry
 from Configurables import GeoSvc
 geoservice = GeoSvc("GeoSvc")
-# if FCC_DETECTORS is empty, this should use relative path to working directory
-path_to_detector = os.environ.get("FCCDETECTORS", "")
+# if K4GEO is empty, this should use relative path to working directory
+path_to_detector = os.environ.get("K4GEO", "")
 print(path_to_detector)
-detectors_to_use=[
-    'Detector/DetFCCeeIDEA-LAr/compact/FCCee_DectMaster_thetamodulemerged.xml',
-    #'Detector/DetFCCeeIDEA-LAr/compact/FCCee_DectMaster.xml',
+detectors_to_use = [
+    'FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/ALLEGRO_o1_v02.xml'
 ]
 
 # prefix all xmls with path_to_detector
-geoservice.detectors = [os.path.join(path_to_detector, _det) for _det in detectors_to_use]
+geoservice.detectors = [os.path.join(
+    path_to_detector, _det) for _det in detectors_to_use]
 geoservice.OutputLevel = INFO
 
 # Geant4 service
 # Configures the Geant simulation: geometry, physics list and user actions
-from Configurables import CreateFCCeeCaloNeighbours
-neighbours = CreateFCCeeCaloNeighbours("neighbours", 
+neighbours = CreateFCCeeCaloNeighbours("neighbours",
                                        outputFileName = "neighbours_map_barrel_thetamodulemerged.root",
-                                       readoutNamesModuleTheta = ["ECalBarrelModuleThetaMerged"],
-#                                       readoutNamesModuleTheta = ["ECalBarrelModuleThetaMerged2"],
-                                       systemNamesModuleTheta = ["system"],
-                                       systemValuesModuleTheta = [4],
-                                       activeFieldNamesModuleTheta = ["layer"],
-                                       activeVolumesNumbers = [12],
-                                       #activeVolumesTheta = [1.2524, 1.2234, 1.1956, 1.1561, 1.1189, 1.0839, 1.0509, 0.9999, 0.9534, 0.91072],
-                                       includeDiagonalCells = False,
-                                       readoutNamesVolumes = [],
-                                       connectBarrels = False, 
-                                       OutputLevel = DEBUG)
+                                       readoutNames=["ECalBarrelModuleThetaMerged"],
+                                       systemNames=["system"],
+                                       systemValues=[4],
+                                       activeFieldNames=["layer"],
+                                       activeVolumesNumbers=[12],
+                                       includeDiagonalCells=False,
+                                       connectBarrels=False,
+                                       OutputLevel=DEBUG)
 
 # ApplicationMgr
-from Configurables import ApplicationMgr
-ApplicationMgr( TopAlg = [],
-                EvtSel = 'NONE',
-                EvtMax   = 1,
-                # order is important, as GeoSvc is needed by G4SimSvc
-                ExtSvc = [geoservice, neighbours],
-                OutputLevel=INFO
-)
+ApplicationMgr(TopAlg=[],
+               EvtSel='NONE',
+               EvtMax=1,
+               # order is important, as GeoSvc is needed by G4SimSvc
+               ExtSvc=[geoservice, neighbours],
+               OutputLevel=INFO
+               )

FCCSW_ecal/noise_map_theta.py

@@ -2,28 +2,25 @@
 # Detector geometry
 from Configurables import GeoSvc
 geoservice = GeoSvc("GeoSvc")
-# if FCC_DETECTORS is empty, this should use relative path to working directory
 import os
-path_to_detector = os.environ.get("FCCDETECTORS", "")
+# if K4GEO is empty, this should use relative path to working directory
+path_to_detector = os.environ.get("K4GEO", "")
 print(path_to_detector)
 detectors_to_use=[
-                    'Detector/DetFCCeeIDEA-LAr/compact/FCCee_DectMaster_thetamodulemerged.xml'
-                  ]
+    'FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/ALLEGRO_o1_v02.xml'
+]
 # prefix all xmls with path_to_detector
 geoservice.detectors = [os.path.join(path_to_detector, _det) for _det in detectors_to_use]
 geoservice.OutputLevel = INFO
 
 ecalBarrelReadoutName = "ECalBarrelModuleThetaMerged"
-#hcalBarrelReadoutName = "ECalBarrelPhiEta"
-hcalBarrelReadoutName = "HCalBarrelReadout"
+hcalBarrelReadoutName = "BarHCal_Readout_phitheta"
 BarrelNoisePath = os.environ['FCCBASEDIR']+"/LAr_scripts/data/elecNoise_ecalBarrelFCCee_theta.root"
 ecalBarrelNoiseHistName = "h_elecNoise_fcc_"
 
 from Configurables import CellPositionsECalBarrelModuleThetaSegTool
 ECalBcells = CellPositionsECalBarrelModuleThetaSegTool("CellPositionsECalBarrel",
                                                        readoutName = ecalBarrelReadoutName)
-#                                                       OutputLevel = DEBUG)
-#print(ECalBcells)
 
 from Configurables import CreateFCCeeCaloNoiseLevelMap, ReadNoiseFromFileTool, ReadNoiseVsThetaFromFileTool
 ECalNoiseTool = ReadNoiseVsThetaFromFileTool("ReadNoiseFromFileToolECal",
@@ -39,28 +36,30 @@
                                              scaleFactor = 1/1000., #MeV to GeV
                                              OutputLevel = INFO)
 
-HCalNoiseTool = ReadNoiseFromFileTool("ReadNoiseFromFileToolHCal",
-                                      readoutName = hcalBarrelReadoutName,
-                                      noiseFileName = BarrelNoisePath,
-                                      elecNoiseHistoName = ecalBarrelNoiseHistName,
-                                      setNoiseOffset = False,
-                                      activeFieldName = "layer",
-                                      addPileup = False,
-                                      numRadialLayers = 12,
-                                      scaleFactor = 1/1000., #MeV to GeV
-                                      OutputLevel = INFO)
+# HCAL noise file has yet to be created/implemented
+# HCalNoiseTool = ReadNoiseFromFileTool("ReadNoiseFromFileToolHCal",
+#                                       readoutName = hcalBarrelReadoutName,
+#                                       noiseFileName = BarrelNoisePath,
+#                                       elecNoiseHistoName = ecalBarrelNoiseHistName,
+#                                       setNoiseOffset = False,
+#                                       activeFieldName = "layer",
+#                                       addPileup = False,
+#                                       numRadialLayers = 12,
+#                                       scaleFactor = 1/1000., #MeV to GeV
+#                                       OutputLevel = INFO)
+HCalNoiseTool = None
 
+# when HCal noise is ready, add also hcal readout info
 noisePerCell = CreateFCCeeCaloNoiseLevelMap("noisePerCell", 
                                             ECalBarrelNoiseTool = ECalNoiseTool, 
                                             ecalBarrelSysId = 4,
                                             HCalBarrelNoiseTool = HCalNoiseTool,
-                                            readoutNamesModuleTheta=[ecalBarrelReadoutName],
-                                            systemNamesModuleTheta=["system"],
-                                            systemValuesModuleTheta=[4],
-                                            activeFieldNamesModuleTheta=["layer"],
+                                            readoutNames=[ecalBarrelReadoutName],
+                                            systemNames=["system"],
+                                            systemValues=[4],
+                                            activeFieldNames=["layer"],
                                             activeVolumesNumbers = [12],
                                             #activeVolumesEta = [1.2524, 1.2234, 1.1956, 1.1561, 1.1189, 1.0839, 1.0509, 0.9999, 0.9534, 0.91072],
-                                            readoutNamesVolumes = [],
                                             outputFileName = "cellNoise_map_electronicsNoiseLevel_thetamodulemerged.root",
                                             OutputLevel = DEBUG)
 

FCCSW_ecal/printNeighbours.C

@@ -1,75 +1,149 @@
-TTree* T = nullptr;
-const std::string filename = "neighbours_map_barrel_thetamodulemerged.root";
+TTree *T = nullptr;
+// const std::string filename = "neighbours_map_barrel_thetamodulemerged.root";
+const std::string filename = "neighbours_map_HCalBarrel.root";
 const std::string treename = "neighbours";
 ULong64_t cID;
-std::vector<unsigned long> *neighbours=0;
+std::vector<unsigned long> *neighbours = 0;
+
+bool useHCalReadoutWithRows = false;
 
 // HELPER FUNCTIONS
+ULong_t ReadNbitsAtPositionMFromCellID(int n, int m, ULong_t cellID)
+{
+  const ULong_t mask = (1<<n) - 1;
+  return (cellID >> m) & mask;
+}
 
+// ECalBarrel: system:4,cryo:1,type:3,subtype:3,layer:8,module:11,theta:10
+// HCalBarrel with row: system:4,layer:5,row:9,theta:9,phi:10
+// HCalBarrel without row: system:4,layer:5,theta:9,phi:10
+
+// extract system ID from cellID
+ULong_t SystemID(ULong_t cellID)
+{
+  return ReadNbitsAtPositionMFromCellID(4, 0, cellID);
+}
 // extract layer number from cellID
-ULong_t Layer(ULong_t cellID) {
-  const ULong_t mask = (1<<8) -1;
-  return (cellID >> 11) & mask;
+ULong_t Layer(ULong_t cellID)
+{
+  if (SystemID(cellID)==4)
+    return ReadNbitsAtPositionMFromCellID(8, 11, cellID);
+  if (SystemID(cellID)==8)
+    return ReadNbitsAtPositionMFromCellID(5, 4, cellID);
+  return 999999999;
 }
 
 // extract module number from cellID
-ULong_t Module(ULong_t cellID) {
-  const ULong_t mask = (1<<11) -1;
-  return (cellID >> 19) & mask;
+ULong_t ECalBarrelModule(ULong_t cellID)
+{
+  if (SystemID(cellID)==4)
+    return ReadNbitsAtPositionMFromCellID(11, 19, cellID);
+  return 999999999;
 }
 
 // extract theta bin from cellID
-ULong_t ThetaBin(ULong_t cellID) {
-  const ULong_t mask = (1<<10) -1;
-  return (cellID >> 30) & mask;
+ULong_t ThetaBin(ULong_t cellID)
+{
+  if (SystemID(cellID)==4)
+    return ReadNbitsAtPositionMFromCellID(10, 30, cellID);
+  if (SystemID(cellID)==8)
+  {
+    if (useHCalReadoutWithRows)
+      return ReadNbitsAtPositionMFromCellID(9, 18, cellID);
+    else
+      return ReadNbitsAtPositionMFromCellID(9, 9, cellID);
+  }
+  return 999999999;
+}
+
+// extract row number from cellID
+ULong_t HCalBarrelRow(ULong_t cellID)
+{
+  if (SystemID(cellID)==8)
+  {
+    if (useHCalReadoutWithRows)
+      return ReadNbitsAtPositionMFromCellID(9, 9, cellID);
+  }
+  return 999999999;
+}
+
+// extract phi number from cellID
+ULong_t HCalBarrelPhiBin(ULong_t cellID)
+{
+  if (SystemID(cellID)==8) {
+    if (useHCalReadoutWithRows)
+      return ReadNbitsAtPositionMFromCellID(10, 27, cellID);
+    else
+      return ReadNbitsAtPositionMFromCellID(10, 18, cellID);
+  }
+  return 999999999;   
 }
 
-void printCell(ULong_t cellID) {
+void printCell(ULong_t cellID)
+{
   cout << "cellID: " << cellID << endl;
-  cout << "Layer: " << Layer(cellID) << endl;
-  cout << "Theta bin: " << ThetaBin(cellID) << endl;
-  cout << "Module: " << Module(cellID) << endl;
+  cout << "System: " << SystemID(cellID) << endl;
+  if (SystemID(cellID)==4) {  
+    cout << "Layer: " << Layer(cellID) << endl;
+    cout << "Theta bin: " << ThetaBin(cellID) << endl;
+    cout << "Module: " << ECalBarrelModule(cellID) << endl;
+  }
+  if (SystemID(cellID)==8) {  
+    cout << "Layer: " << Layer(cellID) << endl;
+    cout << "Row: " << HCalBarrelRow(cellID) << endl;
+    cout << "Theta bin: " << ThetaBin(cellID) << endl;
+    cout << "Phi bin: " << HCalBarrelPhiBin(cellID) << endl;
+  }
   cout << endl;
 }
 
-void LoadNeighboursMap() {
-  if (T==nullptr) {
-    TFile* f = TFile::Open(filename.c_str(),"READ");
-    T = (TTree*) f->Get(treename.c_str());
+void LoadNeighboursMap()
+{
+  if (T == nullptr)
+  {
+    TFile *f = TFile::Open(filename.c_str(), "READ");
+    T = (TTree *)f->Get(treename.c_str());
     T->SetBranchAddress("cellId", &cID);
     T->SetBranchAddress("neighbours", &neighbours);
   }
 }
 
-void printCellAndNeighbours(ULong64_t iEntry) {
+void printCellAndNeighbours(ULong64_t iEntry)
+{
   T->GetEntry(iEntry);
   cout << "=================================================" << endl;
   cout << endl;
   printCell(cID);
-  cout << "Neighbours: " << endl << endl;
-  for (unsigned int i=0; i<neighbours->size(); i++) {
+  cout << "Neighbours: " << endl
+       << endl;
+  for (unsigned int i = 0; i < neighbours->size(); i++)
+  {
     printCell(neighbours->at(i));
   }
   cout << "=================================================" << endl;
 }
 
-void printNeighboursOfCell(ULong_t cellID) {
+void printNeighboursOfCell(ULong_t cellID)
+{
   LoadNeighboursMap();
-  for (ULong64_t iEntry=0; iEntry<T->GetEntries(); iEntry++) {
+  for (ULong64_t iEntry = 0; iEntry < T->GetEntries(); iEntry++)
+  {
     T->GetEntry(iEntry);
-    if (cID == cellID) {
+    if (cID == cellID)
+    {
       printCellAndNeighbours(iEntry);
       return;
     }
   }
   cout << "CellID not found" << endl;
 }
 
-
-void printNeighbours(int n=10) {
+void printNeighbours(int n = 10)
+{
   LoadNeighboursMap();
-  for (int i=0; i<n; i++) {
-    int entry = (int) gRandom->Uniform(T->GetEntries());
+  for (int i = 0; i < n; i++)
+  {
+    int entry = (int)gRandom->Uniform(T->GetEntries());
     printCellAndNeighbours(entry);
   }
 }