Merge pull request #32 from giovannimarchiori/gmarchio-main-20240720

update noise tools

FCCSW_ecal/noise_map_theta.py

@@ -1,14 +1,14 @@
 from Configurables import GeoSvc
 from Configurables import ApplicationMgr
 from Configurables import CreateFCCeeCaloNoiseLevelMap
-# from Configurables import ReadNoiseFromFileTool
-from Configurables import ReadNoiseVsThetaFromFileTool
+#from Configurables import ReadNoiseVsThetaFromFileTool
+from Configurables import NoiseCaloCellsVsThetaFromFileTool
 from Configurables import ConstNoiseTool
 from Configurables import CellPositionsECalBarrelModuleThetaSegTool
 import os
 from Gaudi.Configuration import INFO, DEBUG
 
-doHCal = False
+doHCal = True
 
 # Detector geometry
 geoservice = GeoSvc("GeoSvc")
@@ -29,24 +29,37 @@
 # names of file and histograms with noise per layer vs theta for barrel ECAL
 BarrelNoisePath = os.environ['FCCBASEDIR'] + \
     "/LAr_scripts/data/elecNoise_ecalBarrelFCCee_theta.root"
-ecalBarrelNoiseHistName = "h_elecNoise_fcc_"
+ecalBarrelNoiseRMSHistName = "h_elecNoise_fcc_"
 
 # cell positioning and noise tool for the ecal barrel
 ECalBcells = CellPositionsECalBarrelModuleThetaSegTool("CellPositionsECalBarrel",
                                                        readoutName=ecalBarrelReadoutName)
 
-ECalNoiseTool = ReadNoiseVsThetaFromFileTool("ReadNoiseFromFileToolECal",
-                                             useSegmentation=False,
-                                             cellPositionsTool=ECalBcells,
-                                             readoutName=ecalBarrelReadoutName,
-                                             noiseFileName=BarrelNoisePath,
-                                             elecNoiseHistoName=ecalBarrelNoiseHistName,
-                                             setNoiseOffset=False,
-                                             activeFieldName="layer",
-                                             addPileup=False,
-                                             numRadialLayers=11,
-                                             scaleFactor=1 / 1000.,  # MeV to GeV
-                                             OutputLevel=INFO)
+#ECalNoiseTool = ReadNoiseVsThetaFromFileTool("ReadNoiseFromFileToolECal",
+#                                             useSegmentation=False,
+#                                             cellPositionsTool=ECalBcells,
+#                                             readoutName=ecalBarrelReadoutName,
+#                                             noiseFileName=BarrelNoisePath,
+#                                             elecNoiseHistoName=ecalBarrelNoiseHistName,
+#                                             setNoiseOffset=False,
+#                                             activeFieldName="layer",
+#                                             addPileup=False,
+#                                             numRadialLayers=11,
+#                                             scaleFactor=1 / 1000.,  # MeV to GeV
+#                                             OutputLevel=INFO)
+
+ECalNoiseTool = NoiseCaloCellsVsThetaFromFileTool("NoiseCaloCellsVsThetaFromFileTool",
+                                                  cellPositionsTool=ECalBcells,
+                                                  readoutName=ecalBarrelReadoutName,
+                                                  noiseFileName=BarrelNoisePath,
+                                                  elecNoiseRMSHistoName=ecalBarrelNoiseRMSHistName,
+                                                  setNoiseOffset=False,
+                                                  activeFieldName="layer",
+                                                  addPileup=False,
+                                                  numRadialLayers=11,
+                                                  scaleFactor=1 / 1000.,  # MeV to GeV
+                                                  OutputLevel=DEBUG)
+
 
 if doHCal:
     # noise tool for the HCAL barrel
@@ -63,9 +76,17 @@
     #                                       OutputLevel = INFO)
     # ConstNoiseTool provides constant noise for all calo subsystems
     # here we are going to use it only for hcal barrel
-    HCalNoiseTool = ConstNoiseTool("ConstNoiseTool")
+    HCalNoiseTool = ConstNoiseTool("ConstNoiseTool",
+                                   detectors = [ "HCAL_Barrel"],
+                                   detectorsNoiseRMS = [0.0115/4],
+                                   OutputLevel = DEBUG)
+    #HCalNoiseTool = ConstNoiseTool("ConstNoiseTool",
+    #                               detectors = ["ECAL_Barrel", "ECAL_Endcap", "HCAL_Barrel", "HCAL_Endcap"],
+    #                               detectorsNoiseRMS = [0.0075/4, 0.0075/4, 0.0115/4, 0.0115/4],
+    #                               OutputLevel = DEBUG)
 
     # create the noise file for ECAL+HCAL barrel cells
+    # the tool wants the system IDs, maybe we could have passed the names instead
     noisePerCell = CreateFCCeeCaloNoiseLevelMap("noisePerCell",
                                                 ECalBarrelNoiseTool=ECalNoiseTool,
                                                 ecalBarrelSysId=4,
@@ -79,7 +100,6 @@
                                                 activeFieldNames=[
                                                     "layer", "layer"],
                                                 activeVolumesNumbers=[11, 13],
-                                                # activeVolumesEta = [1.2524, 1.2234, 1.1956, 1.1561, 1.1189, 1.0839, 1.0509, 0.9999, 0.9534, 0.91072],
                                                 activeVolumesTheta=[
                                                     [],
                                                     [

FCCSW_ecal/run_thetamodulemerged.py

@@ -1,59 +1,83 @@
+#
+# IMPORTS
+#
 from Configurables import ApplicationMgr
 from Configurables import EventCounter
 from Configurables import AuditorSvc, ChronoAuditor
-from Configurables import PodioOutput
-from Configurables import CaloTowerToolFCCee
-from Configurables import CreateCaloClustersSlidingWindowFCCee
-from Configurables import CorrectCaloClusters
-from Configurables import CalibrateCaloClusters
-from Configurables import AugmentClustersFCCee
-from Configurables import CreateEmptyCaloCellsCollection
-from Configurables import CreateCaloCellPositionsFCCee
-from Configurables import CellPositionsECalBarrelModuleThetaSegTool
-from Configurables import RedoSegmentation
-from Configurables import CreateCaloCells
-from Configurables import CalibrateCaloHitsTool
-from Configurables import CalibrateInLayersTool
+# Event generation
+from Configurables import GenAlg
+# G4 simulation
 from Configurables import SimG4Alg
 from Configurables import SimG4PrimariesFromEdmTool
 from Configurables import SimG4SaveCalHits
 from Configurables import SimG4ConstantMagneticFieldTool
 from Configurables import SimG4Svc
 from Configurables import SimG4FullSimActions
 from Configurables import SimG4SaveParticleHistory
+# Geometry
 from Configurables import GeoSvc
-from Configurables import HepMCToEDMConverter
-from Configurables import GenAlg
+# Input/output
 from Configurables import FCCDataSvc
+from Configurables import PodioOutput
+from Configurables import HepMCToEDMConverter
+# Create cells
+from Configurables import CreateCaloCells
+from Configurables import CreateEmptyCaloCellsCollection
+# Cell positioning tools
+from Configurables import CreateCaloCellPositionsFCCee
+from Configurables import CellPositionsECalBarrelModuleThetaSegTool
+# Redo segmentation for ECAL and HCAL
+from Configurables import RedoSegmentation
+# Read noise values from file and generate noise in cells
+from Configurables import NoiseCaloCellsVsThetaFromFileTool
+# Apply sampling fraction corrections
+from Configurables import CalibrateCaloHitsTool
+from Configurables import CalibrateInLayersTool
+# Up/down stream correction
+from Configurables import CorrectCaloClusters
+# SW clustering
+from Configurables import CaloTowerToolFCCee
+from Configurables import CreateCaloClustersSlidingWindowFCCee
+# Topo clustering
 from Configurables import CaloTopoClusterInputTool
 from Configurables import TopoCaloNeighbours
 from Configurables import TopoCaloNoisyCells
 from Configurables import CaloTopoClusterFCCee
-from Gaudi.Configuration import INFO
-# , VERBOSE, DEBUG
-# from Gaudi.Configuration import *
-
-import os
-
+# Decorate clusters with shower shape parameters
+from Configurables import AugmentClustersFCCee
+# MVA calibration
+from Configurables import CalibrateCaloClusters
+# Logger
+from Gaudi.Configuration import INFO, VERBOSE, DEBUG
+# Units and physical constants
 from GaudiKernel.SystemOfUnits import GeV, tesla, mm
 from GaudiKernel.PhysicalConstants import pi, halfpi, twopi
+# python libraries
+import os
 from math import cos, sin, tan
 
+#
+# SETTINGS
+#
 
-# general settings
+# - general settings
+#
 use_pythia = False  # use pythia or particle gun
-addNoise = False  # add noise or not to the cell energy
-dumpGDML = False  # create GDML file of detector model
-runHCal = False  # simulate only the ECAL or both ECAL+HCAL
+addNoise = False     # add noise or not to the cell energy
+dumpGDML = False    # create GDML file of detector model
+runHCal = False     # simulate only the ECAL or both ECAL+HCAL
 
+# - what to save in output file
+#
 # for big productions, save significant space removing hits and cells
 # however, hits and cluster cells might be wanted for small productions for detailed event displays
 # also, cluster cells are needed for the MVA training
 saveHits = True
 saveCells = True
 saveClusterCells = True
 
-# clustering
+# - clustering
+#
 doSWClustering = True
 doTopoClustering = True
 
@@ -64,12 +88,13 @@
 
 # BDT regression from total cluster energy and fraction of energy in each layer (after correction for sampling fraction)
 # not to be applied for ECAL+HCAL clustering (MVA trained only on ECAL)
-applyMVAClusterEnergyCalibration = True and not runHCal
+applyMVAClusterEnergyCalibration = False and not runHCal
 
 # calculate cluster energy and barycenter per layer and save it as extra parameters
 addShapeParameters = True and not runHCal
 
-# Input for simulations (momentum is expected in GeV!)
+# - input for simulations (momentum is expected in GeV!)
+#
 
 # Parameters for the particle gun simulations, dummy if use_pythia is set
 # to True
@@ -120,8 +145,11 @@
 magneticField = False
 
 
-podioevent = FCCDataSvc("EventDataSvc")
+#
+# ALGORITHMS AND SERVICES SETUP
+#
 
+podioevent = FCCDataSvc("EventDataSvc")
 
 # Particle gun setup
 genAlg = GenAlg()
@@ -242,7 +270,7 @@
 # ECAL
 ecalBarrelReadoutName = "ECalBarrelModuleThetaMerged"  # barrel, original segmentation (baseline)
 ecalBarrelReadoutName2 = "ECalBarrelModuleThetaMerged2"  # barrel, after re-segmentation (for optimisation studies)
-ecalEndcapReadoutName = "ECalEndcapPhiEta"  # endcap, from FCC-hh sim, to be replaced by proper readout/detector geometry
+ecalEndcapReadoutName = "ECalEndcapTurbine"  # endcap, turbine-like (baseline)
 # HCAL
 if runHCal:
     hcalBarrelReadoutName = "HCalBarrelReadout"  # barrel, layer-row-theta-phi based (can be used to fill various cell collections with different readouts)
@@ -318,7 +346,7 @@
                                         readoutName=ecalBarrelReadoutName,
                                         layerFieldName="layer")
 
-# - ECAL endcap
+# - ECAL endcap (TO BE UPDATED)
 calibEcalEndcap = CalibrateCaloHitsTool(
     "CalibrateECalEndcap", invSamplingFraction="4.27")
 
@@ -415,6 +443,55 @@
 createEcalEndcapCells.hits.Path = "ECalEndcapHits"
 createEcalEndcapCells.cells.Path = "ECalEndcapCells"
 
+if addNoise:
+    ecalBarrelNoisePath = os.environ['FCCBASEDIR']+"/LAr_scripts/data/elecNoise_ecalBarrelFCCee_theta.root"
+    ecalBarrelNoiseRMSHistName = "h_elecNoise_fcc_"
+    from Configurables import NoiseCaloCellsVsThetaFromFileTool
+    noiseBarrel = NoiseCaloCellsVsThetaFromFileTool("NoiseBarrel",
+                                                    cellPositionsTool=cellPositionEcalBarrelTool,
+                                                    readoutName=ecalBarrelReadoutName,
+                                                    noiseFileName=ecalBarrelNoisePath,
+                                                    elecNoiseRMSHistoName=ecalBarrelNoiseRMSHistName,
+                                                    setNoiseOffset=False,
+                                                    activeFieldName="layer",
+                                                    addPileup=False,
+                                                    filterNoiseThreshold=0,
+                                                    numRadialLayers=11,
+                                                    scaleFactor=1 / 1000.,  # MeV to GeV
+                                                    OutputLevel=DEBUG)
+
+    # needs to be migrated!
+    #from Configurables import TubeLayerPhiEtaCaloTool
+    #barrelGeometry = TubeLayerPhiEtaCaloTool("EcalBarrelGeo",
+    #                                         readoutName=ecalBarrelReadoutNamePhiEta,
+    #                                         activeVolumeName="LAr_sensitive",
+    #                                         activeFieldName="layer",
+    #                                         activeVolumesNumber=12,
+    #                                         fieldNames=["system"],
+    #                                         fieldValues=[4])
+
+    # cells with noise not filtered
+    # createEcalBarrelCellsNoise = CreateCaloCells("CreateECalBarrelCellsNoise",
+    #                                              doCellCalibration=False,
+    #                                              addCellNoise=True,
+    #                                              filterCellNoise=False,
+    #                                              OutputLevel=INFO,
+    #                                              hits="ECalBarrelCellsStep2",
+    #                                              noiseTool=noiseBarrel,
+    #                                              geometryTool=barrelGeometry,
+    #                                              cells=EcalBarrelCellsName)
+
+    # cells with noise filtered
+    # createEcalBarrelCellsNoise = CreateCaloCells("CreateECalBarrelCellsNoise_filtered",
+    #                                              doCellCalibration=False,
+    #                                              addCellNoise=True,
+    #                                              filterCellNoise=True,
+    #                                              OutputLevel=INFO,
+    #                                              hits="ECalBarrelCellsStep2",
+    #                                              noiseTool=noiseBarrel,
+    #                                              geometryTool=barrelGeometry,
+    #                                              cells=EcalBarrelCellsName)
+
 if runHCal:
     # Create cells in HCal barrel
     # 1 - merge hits into cells with the default readout
@@ -459,7 +536,7 @@
                                          oldReadoutName=hcalBarrelReadoutName,
                                          # specify which fields are going to be altered (deleted/rewritten)
                                          oldSegmentationIds=["row", "theta", "phi"],
-                                         # new bitfield (readout), with new segmentation (merged modules and theta cells)
+                                         # new bitfield (readout), with new segmentation (theta-phi grid)
                                          newReadoutName=hcalBarrelReadoutName2,
                                          OutputLevel=INFO,
                                          debugPrint=200,
@@ -500,6 +577,8 @@
     # createHcalEndcapCells.hits.Path="HCalEndcapHits"
     # createHcalEndcapCells.cells.Path="HCalEndcapCells"
 
+    # TODO: noise...
+
 else:
     hcalBarrelCellsName = "emptyCaloCells"
     hcalBarrelPositionedCellsName = "emptyCaloCells"
@@ -565,6 +644,7 @@
                                                   numLayers=[11],
                                                   firstLayerIDs=[0],
                                                   lastLayerIDs=[10],
+                                                  systemIDs=[4],
                                                   readoutNames=[ecalBarrelReadoutName],
                                                   # do not split the following line or it will break scripts that update the values of the corrections
                                                   upstreamParameters = [[0.025582045561310333, -0.9524128168665387, -53.10089405478649, 1.283851527438571, -295.30650178662637, -284.8945817377308]],
@@ -683,6 +763,7 @@
             "correctCaloTopoClusters",
             inClusters=createTopoClusters.clusters.Path,
             outClusters="Corrected" + createTopoClusters.clusters.Path,
+            systemIDs=[4],
             numLayers=[11],
             firstLayerIDs=[0],
             lastLayerIDs=[10],