Merge pull request #13 from giovannimarchiori/gmarchio-main-2023-09-0…

…8-newclustering

Make topoclustering work with the new theta-module readout for the FCC-ee calorimeter

FCCSW_ecal/neighbours_theta.py

@@ -0,0 +1,44 @@
+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", "")
+print(path_to_detector)
+detectors_to_use=[
+    'Detector/DetFCCeeIDEA-LAr/compact/FCCee_DectMaster_thetamodulemerged.xml',
+    #'Detector/DetFCCeeIDEA-LAr/compact/FCCee_DectMaster.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
+
+# Geant4 service
+# Configures the Geant simulation: geometry, physics list and user actions
+from Configurables import CreateFCCeeCaloNeighbours
+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)
+
+# 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
+)

FCCSW_ecal/noise_map_theta.py

@@ -0,0 +1,75 @@
+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
+import os
+path_to_detector = os.environ.get("FCCDETECTORS", "")
+print(path_to_detector)
+detectors_to_use=[
+                    'Detector/DetFCCeeIDEA-LAr/compact/FCCee_DectMaster_thetamodulemerged.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"
+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",
+                                             useSegmentation = False,
+                                             cellPositionsTool = ECalBcells,
+                                             readoutName = ecalBarrelReadoutName,
+                                             noiseFileName = BarrelNoisePath,
+                                             elecNoiseHistoName = ecalBarrelNoiseHistName,
+                                             setNoiseOffset = False,
+                                             activeFieldName = "layer",
+                                             addPileup = False,
+                                             numRadialLayers = 12,
+                                             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)
+
+noisePerCell = CreateFCCeeCaloNoiseLevelMap("noisePerCell", 
+                                            ECalBarrelNoiseTool = ECalNoiseTool, 
+                                            ecalBarrelSysId = 4,
+                                            HCalBarrelNoiseTool = HCalNoiseTool,
+                                            readoutNamesModuleTheta=[ecalBarrelReadoutName],
+                                            systemNamesModuleTheta=["system"],
+                                            systemValuesModuleTheta=[4],
+                                            activeFieldNamesModuleTheta=["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)
+
+# ApplicationMgr
+from Configurables import ApplicationMgr
+ApplicationMgr( TopAlg = [],
+                EvtSel = 'NONE',
+                EvtMax   = 1,
+                # order is important, as GeoSvc is needed by G4SimSvc
+                ExtSvc = [geoservice, noisePerCell],
+                OutputLevel=INFO
+               )

FCCSW_ecal/printNeighbours.C

@@ -0,0 +1,75 @@
+TTree* T = nullptr;
+const std::string filename = "neighbours_map_barrel_thetamodulemerged.root";
+const std::string treename = "neighbours";
+ULong64_t cID;
+std::vector<unsigned long> *neighbours=0;
+
+// HELPER FUNCTIONS
+
+// extract layer number from cellID
+ULong_t Layer(ULong_t cellID) {
+  const ULong_t mask = (1<<8) -1;
+  return (cellID >> 11) & mask;
+}
+
+// extract module number from cellID
+ULong_t Module(ULong_t cellID) {
+  const ULong_t mask = (1<<11) -1;
+  return (cellID >> 19) & mask;
+}
+
+// extract theta bin from cellID
+ULong_t ThetaBin(ULong_t cellID) {
+  const ULong_t mask = (1<<10) -1;
+  return (cellID >> 30) & mask;
+}
+
+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 << endl;
+}
+
+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) {
+  T->GetEntry(iEntry);
+  cout << "=================================================" << endl;
+  cout << endl;
+  printCell(cID);
+  cout << "Neighbours: " << endl << endl;
+  for (unsigned int i=0; i<neighbours->size(); i++) {
+    printCell(neighbours->at(i));
+  }
+  cout << "=================================================" << endl;
+}
+
+void printNeighboursOfCell(ULong_t cellID) {
+  LoadNeighboursMap();
+  for (ULong64_t iEntry=0; iEntry<T->GetEntries(); iEntry++) {
+    T->GetEntry(iEntry);
+    if (cID == cellID) {
+      printCellAndNeighbours(iEntry);
+      return;
+    }
+  }
+  cout << "CellID not found" << endl;
+}
+
+
+void printNeighbours(int n=10) {
+  LoadNeighboursMap();
+  for (int i=0; i<n; i++) {
+    int entry = (int) gRandom->Uniform(T->GetEntries());
+    printCellAndNeighbours(entry);
+  }
+}

FCCSW_ecal/runClueAndTopoAndSlidingWindowAndCaloSim.py

@@ -373,6 +373,7 @@
                                      #positionsHECTool = HECcells,
                                      #positionsEMFwdTool = ECalFwdcells,
                                      #positionsHFwdTool = HCalFwdcells,
+                                     readoutName = ecalBarrelReadoutNamePhiEta,
                                      seedSigma = 4,
                                      neighbourSigma = 2,
                                      lastNeighbourSigma = 0,

FCCSW_ecal/runTopoAndSlidingWindowAndCaloSim.py

@@ -399,6 +399,7 @@
                                      #positionsHECTool = HECcells,
                                      #positionsEMFwdTool = ECalFwdcells,
                                      #positionsHFwdTool = HCalFwdcells,
+                                     readoutName = ecalBarrelReadoutNamePhiEta,
                                      seedSigma = 4,
                                      neighbourSigma = 2,
                                      lastNeighbourSigma = 0,