From eb0e0771bf8a39e2985dde1e162595138aad9605 Mon Sep 17 00:00:00 2001 From: Erich Varnes Date: Wed, 7 Aug 2024 06:41:43 -0700 Subject: [PATCH] update test scripts --- ...el_ReconstructionTopoClusters_crosstalk.py | 445 ++++++ .../tests/options/run_thetamodulemerged.py | 448 +++--- .../run_thetamodulemerged_SW_benchmarkCorr.py | 1354 +++++++++-------- 3 files changed, 1331 insertions(+), 916 deletions(-) create mode 100644 RecFCCeeCalorimeter/tests/options/runEcalBarrel_ReconstructionTopoClusters_crosstalk.py diff --git a/RecFCCeeCalorimeter/tests/options/runEcalBarrel_ReconstructionTopoClusters_crosstalk.py b/RecFCCeeCalorimeter/tests/options/runEcalBarrel_ReconstructionTopoClusters_crosstalk.py new file mode 100644 index 00000000..80056b9d --- /dev/null +++ b/RecFCCeeCalorimeter/tests/options/runEcalBarrel_ReconstructionTopoClusters_crosstalk.py @@ -0,0 +1,445 @@ +# +# IMPORTS +# +from Configurables import ApplicationMgr +#from Configurables import EventCounter +from Configurables import AuditorSvc, ChronoAuditor +# Input/output +from Configurables import k4DataSvc, PodioInput +from Configurables import PodioOutput +# Geometry +from Configurables import GeoSvc +# Create cells +from Configurables import CreateCaloCells +from Configurables import CreateEmptyCaloCellsCollection +# Cell positioning tools +from Configurables import CreateCaloCellPositionsFCCee +from Configurables import CellPositionsECalBarrelModuleThetaSegTool +from Configurables import CellPositionsECalEndcapTurbineSegTool +# Redo segmentation for ECAL +from Configurables import RedoSegmentation +# Change HCAL segmentation +from Configurables import RewriteBitfield +# 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 +# Decorate clusters with shower shape parameters +from Configurables import AugmentClustersFCCee +# Read crosstalk map +from Configurables import ReadCaloCrosstalkMap +# 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 +# +inputfile = "https://fccsw.web.cern.ch/fccsw/filesForSimDigiReco/ALLEGRO/ALLEGRO_o1_v03/ALLEGRO_sim.root" +# note - this file probably contains the old ecal endcap segmentation so we disable the endcap digitisation later +Nevts = 50 # -1 means all events +dumpGDML = False + +# - 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 = False +# saveCells = False +saveHits = True +saveCells = True +saveClusterCells = True +doCrosstalk = True # switch on/off the crosstalk + +# ECAL barrel parameters for digitisation +samplingFraction=[0.37586625991994105] * 1 + [0.13459486704309379] * 1 + [0.142660085165352] * 1 + [0.14768106642302886] * 1 + [0.15205230356024715] * 1 + [0.15593671843591686] * 1 + [0.15969313426201745] * 1 + [0.16334257010426537] * 1 + [0.16546584993953908] * 1 + [0.16930439771304764] * 1 + [0.1725913708958098] * 1 +upstreamParameters = [[0.025582045561310333, -0.9524128168665387, -53.10089405478649, 1.283851527438571, -295.30650178662637, -284.8945817377308]] +downstreamParameters = [[0.0018280333929494054, 0.004932212590963076, 0.8409676097173655, -1.2676690014715288, 0.005347798049886769, 4.161741293789687]] + +ecalBarrelLayers = len(samplingFraction) +resegmentECalBarrel = False + +# - parameters for clustering +# +doSWClustering = True +doTopoClustering = True + +# calculate cluster energy and barycenter per layer and save it as extra parameters +addShapeParameters = True +ecalBarrelThetaWeights = [-1, 3.0, 3.0, 3.0, 4.25, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0] # to be recalculated for V03, separately for topo and calo clusters... + +# +# ALGORITHMS AND SERVICES SETUP +# + +# Input: load the output of the SIM step +evtsvc = k4DataSvc('EventDataSvc') +evtsvc.input = inputfile +input_reader = PodioInput('InputReader') +podioevent = k4DataSvc("EventDataSvc") + +# Detector geometry +# prefix all xmls with path_to_detector +# if K4GEO is empty, this should use relative path to working directory +geoservice = GeoSvc("GeoSvc") +path_to_detector = os.environ.get("K4GEO", "") +detectors_to_use = [ + 'FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ALLEGRO_o1_v03.xml' +] +geoservice.detectors = [ + os.path.join(path_to_detector, _det) for _det in detectors_to_use +] +geoservice.OutputLevel = INFO + +# GDML dump of detector model +if dumpGDML: + from Configurables import GeoToGdmlDumpSvc + gdmldumpservice = GeoToGdmlDumpSvc("GeoToGdmlDumpSvc") + +# Digitisation (merging hits into cells, EM scale calibration via sampling fractions) + +# - ECAL readouts +ecalBarrelReadoutName = "ECalBarrelModuleThetaMerged" +ecalBarrelReadoutName2 = "ECalBarrelModuleThetaMerged2" +ecalEndcapReadoutName = "ECalEndcapTurbine" + +hcalBarrelReadoutName = "" +hcalBarrelReadoutName2 = "" +hcalEndcapReadoutName = "" + +# - EM scale calibration (sampling fraction) +# * ECAL barrel +calibEcalBarrel = CalibrateInLayersTool("CalibrateECalBarrel", + samplingFraction=samplingFraction, + readoutName=ecalBarrelReadoutName, + layerFieldName="layer") +# * ECAL endcap +calibEcalEndcap = CalibrateCaloHitsTool( + "CalibrateECalEndcap", invSamplingFraction="4.27") # FIXME: to be updated for ddsim + +# Create cells in ECal barrel (needed if one wants to apply cell calibration, +# which is not performed by ddsim) + +# read the crosstalk map +readCrosstalkMap = ReadCaloCrosstalkMap("ReadCrosstalkMap", + fileName="https://fccsw.web.cern.ch/fccsw/filesForSimDigiReco/ALLEGRO/ALLEGRO_o1_v03/xtalk_neighbours_map_ecalB_thetamodulemerged.root", + OutputLevel=INFO) + +# - merge hits into cells according to initial segmentation +ecalBarrelCellsName = "ECalBarrelCells" +createEcalBarrelCells = CreateCaloCells("CreateECalBarrelCells", + doCellCalibration=True, + calibTool=calibEcalBarrel, + crosstalksTool=readCrosstalkMap, + addCrosstalk=doCrosstalk, + addCellNoise=False, + filterCellNoise=False, + addPosition=True, + OutputLevel=INFO, + hits=ecalBarrelReadoutName, + cells=ecalBarrelCellsName) + +# - add to Ecal barrel cells the position information +# (good for physics, all coordinates set properly) +cellPositionEcalBarrelTool = CellPositionsECalBarrelModuleThetaSegTool( + "CellPositionsECalBarrel", + readoutName=ecalBarrelReadoutName, + OutputLevel=INFO +) +ecalBarrelPositionedCellsName = ecalBarrelReadoutName + "Positioned" +createEcalBarrelPositionedCells = CreateCaloCellPositionsFCCee( + "CreateECalBarrelPositionedCells", + OutputLevel=INFO +) +createEcalBarrelPositionedCells.positionsTool = cellPositionEcalBarrelTool +createEcalBarrelPositionedCells.hits.Path = ecalBarrelCellsName +createEcalBarrelPositionedCells.positionedHits.Path = ecalBarrelPositionedCellsName + +# Create cells in ECal endcap (needed if one wants to apply cell calibration, +# which is not performed by ddsim) +#ecalEndcapCellsName = "ECalEndcapCells" +#createEcalEndcapCells = CreateCaloCells("CreateEcalEndcapCaloCells", +# doCellCalibration=True, +# calibTool=calibEcalEndcap, +# addCellNoise=False, +# filterCellNoise=False, +# OutputLevel=INFO, +# hits=ecalEndcapReadoutName, +# cells=ecalEndcapCellsName) + +# Add to Ecal endcap cells the position information +# (good for physics, all coordinates set properly) +#cellPositionEcalEndcapTool = CellPositionsECalEndcapTurbineSegTool( +# "CellPositionsECalEndcap", +# readoutName=ecalEndcapReadoutName, +# OutputLevel=INFO +#) +#ecalEndcapPositionedCellsName = "ECalEndcapPositionedCells" +#createEcalEndcapPositionedCells = CreateCaloCellPositionsFCCee( +# "CreateECalEndcapPositionedCells", +# OutputLevel=INFO +#) +#createEcalEndcapPositionedCells.positionsTool = cellPositionEcalEndcapTool +#createEcalEndcapPositionedCells.hits.Path = ecalEndcapCellsName +#createEcalEndcapPositionedCells.positionedHits.Path = ecalEndcapPositionedCellsName + +hcalBarrelCellsName = "emptyCaloCells" +hcalBarrelPositionedCellsName = "emptyCaloCells" +hcalBarrelCellsName2 = "emptyCaloCells" +hcalBarrelPositionedCellsName2 = "emptyCaloCells" +cellPositionHcalBarrelTool = None +cellPositionHcalBarrelTool2 = None + +# Empty cells for parts of calorimeter not implemented yet +createemptycells = CreateEmptyCaloCellsCollection("CreateEmptyCaloCells") +createemptycells.cells.Path = "emptyCaloCells" + +# Produce sliding window clusters +if doSWClustering: + towers = CaloTowerToolFCCee("towers", + deltaThetaTower=4 * 0.009817477 / 4, deltaPhiTower=2 * 2 * pi / 1536., + ecalBarrelReadoutName=ecalBarrelReadoutName, + #ecalEndcapReadoutName=ecalEndcapReadoutName, + ecalEndcapReadoutName="", + ecalFwdReadoutName="", + hcalBarrelReadoutName=hcalBarrelReadoutName2, + hcalExtBarrelReadoutName="", + hcalEndcapReadoutName="", + hcalFwdReadoutName="", + OutputLevel=INFO) + towers.ecalBarrelCells.Path = ecalBarrelPositionedCellsName + #towers.ecalEndcapCells.Path = ecalEndcapPositionedCellsName + towers.ecalEndcapCells.Path = "emptyCaloCells" + towers.ecalFwdCells.Path = "emptyCaloCells" + towers.hcalBarrelCells.Path = hcalBarrelPositionedCellsName2 + towers.hcalExtBarrelCells.Path = "emptyCaloCells" + towers.hcalEndcapCells.Path = "emptyCaloCells" + towers.hcalFwdCells.Path = "emptyCaloCells" + + # Cluster variables + windT = 9 + windP = 17 + posT = 5 + posP = 11 + dupT = 7 + dupP = 13 + finT = 9 + finP = 17 + # Minimal energy to create a cluster in GeV (FCC-ee detectors have to reconstruct low energy particles) + threshold = 0.040 + + createClusters = CreateCaloClustersSlidingWindowFCCee("CreateClusters", + towerTool=towers, + nThetaWindow=windT, nPhiWindow=windP, + nThetaPosition=posT, nPhiPosition=posP, + nThetaDuplicates=dupT, nPhiDuplicates=dupP, + nThetaFinal=finT, nPhiFinal=finP, + energyThreshold=threshold, + energySharingCorrection=False, + attachCells=True, + OutputLevel=INFO + ) + createClusters.clusters.Path = "CaloClusters" + createClusters.clusterCells.Path = "CaloClusterCells" + + if addShapeParameters: + augmentCaloClusters = AugmentClustersFCCee("augmentCaloClusters", + inClusters=createClusters.clusters.Path, + outClusters="Augmented" + createClusters.clusters.Path, + systemIDs=[4], + systemNames=["EMB"], + numLayers=[ecalBarrelLayers], + readoutNames=[ecalBarrelReadoutName], + layerFieldNames=["layer"], + thetaRecalcWeights=[ecalBarrelThetaWeights], + OutputLevel=INFO + ) + +if doTopoClustering: + # Produce topoclusters (ECAL only or ECAL+HCAL) + createTopoInput = CaloTopoClusterInputTool("CreateTopoInput", + ecalBarrelReadoutName=ecalBarrelReadoutName, + ecalEndcapReadoutName="", + ecalFwdReadoutName="", + hcalBarrelReadoutName=hcalBarrelReadoutName2, + hcalExtBarrelReadoutName="", + hcalEndcapReadoutName="", + hcalFwdReadoutName="", + OutputLevel=INFO) + + createTopoInput.ecalBarrelCells.Path = ecalBarrelPositionedCellsName + createTopoInput.ecalEndcapCells.Path = "emptyCaloCells" + createTopoInput.ecalFwdCells.Path = "emptyCaloCells" + createTopoInput.hcalBarrelCells.Path = hcalBarrelPositionedCellsName2 + createTopoInput.hcalExtBarrelCells.Path = "emptyCaloCells" + createTopoInput.hcalEndcapCells.Path = "emptyCaloCells" + createTopoInput.hcalFwdCells.Path = "emptyCaloCells" + cellPositionHcalBarrelNoSegTool = None + cellPositionHcalExtBarrelTool = None + + neighboursMap = "https://fccsw.web.cern.ch/fccsw/filesForSimDigiReco/ALLEGRO/ALLEGRO_o1_v03/neighbours_map_ecalB_thetamodulemerged.root" + noiseMap = "https://fccsw.web.cern.ch/fccsw/filesForSimDigiReco/ALLEGRO/ALLEGRO_o1_v03/cellNoise_map_electronicsNoiseLevel_ecalB_thetamodulemerged.root" + + readNeighboursMap = TopoCaloNeighbours("ReadNeighboursMap", + fileName=neighboursMap, + OutputLevel=INFO) + + # Noise levels per cell + readNoisyCellsMap = TopoCaloNoisyCells("ReadNoisyCellsMap", + fileName=noiseMap, + OutputLevel=INFO) + + createTopoClusters = CaloTopoClusterFCCee("CreateTopoClusters", + TopoClusterInput=createTopoInput, + # expects neighbours map from cellid->vec < neighbourIds > + neigboursTool=readNeighboursMap, + # tool to get noise level per cellid + noiseTool=readNoisyCellsMap, + # cell positions tools for all sub - systems + positionsECalBarrelTool=cellPositionEcalBarrelTool, + positionsHCalBarrelTool=cellPositionHcalBarrelTool2, + # positionsHCalBarrelNoSegTool=cellPositionHcalBarrelNoSegTool, + # positionsHCalExtBarrelTool=cellPositionHcalExtBarrelTool, + # positionsHCalExtBarrelTool = HCalExtBcells, + # positionsEMECTool = EMECcells, + # positionsHECTool = HECcells, + # positionsEMFwdTool = ECalFwdcells, + # positionsHFwdTool = HCalFwdcells, + noSegmentationHCal=False, + seedSigma=4, + neighbourSigma=2, + lastNeighbourSigma=0, + OutputLevel=INFO) + createTopoClusters.clusters.Path = "CaloTopoClusters" + createTopoClusters.clusterCells.Path = "CaloTopoClusterCells" + + + # Correction below is for EM-only clusters + # Need something different for EM+HCAL + if addShapeParameters: + augmentCaloTopoClusters = AugmentClustersFCCee("augmentCaloTopoClusters", + inClusters=createTopoClusters.clusters.Path, + outClusters="Augmented" + createTopoClusters.clusters.Path, + systemIDs=[4], + systemNames=["EMB"], + numLayers=[ecalBarrelLayers], + readoutNames=[ecalBarrelReadoutName], + layerFieldNames=["layer"], + thetaRecalcWeights=[ecalBarrelThetaWeights], + OutputLevel=INFO) +# Output +out = PodioOutput("out", + OutputLevel=INFO) +out.filename = "ALLEGRO_sim_digi_reco.root" + +# drop the unpositioned ECal barrel cells +out.outputCommands = ["keep *", "drop HCal*", "drop emptyCaloCells", "drop ECalBarrelCells*"] +out.outputCommands.append("drop %s" % ecalBarrelReadoutName) +out.outputCommands.append("drop %s" % ecalBarrelReadoutName2) +out.outputCommands.append("drop ECalBarrelCellsMerged") + +# drop lumi, vertex, DCH, Muons (unless want to keep for event display) +out.outputCommands.append("drop Lumi*") +# out.outputCommands.append("drop Vertex*") +# out.outputCommands.append("drop DriftChamber_simHits*") +out.outputCommands.append("drop MuonTagger*") + +# drop hits/positioned cells/cluster cells if desired +if not saveHits: + out.outputCommands.append("drop %s_contributions" % ecalBarrelReadoutName) + out.outputCommands.append("drop %s_contributions" % ecalBarrelReadoutName2) +if not saveCells: + out.outputCommands.append("drop %s" % ecalBarrelPositionedCellsName) + out.outputCommands.append("drop %s" % ecalBarrelPositionedCellsName2) + +if not saveClusterCells: + out.outputCommands.append("drop Calo*ClusterCells*") + +# CPU information +chra = ChronoAuditor() +audsvc = AuditorSvc() +audsvc.Auditors = [chra] +out.AuditExecute = True + +# Event counter +#event_counter = EventCounter('event_counter') +#event_counter.Frequency = 10 + +# Configure list of external services +ExtSvc = [geoservice, podioevent, audsvc] +if dumpGDML: + ExtSvc += [gdmldumpservice] + +# Setup alg sequence +TopAlg = [ +# event_counter, + input_reader, + createEcalBarrelCells, + createEcalBarrelPositionedCells, +# createEcalEndcapCells, +# createEcalEndcapPositionedCells +] +createEcalBarrelCells.AuditExecute = True +createEcalBarrelPositionedCells.AuditExecute = True +#createEcalEndcapCells.AuditExecute = True + +if resegmentECalBarrel: + TopAlg += [ + resegmentEcalBarrelTool, + createEcalBarrelCells2, + createEcalBarrelPositionedCells2, + ] + resegmentEcalBarrelTool.AuditExecute = True + createEcalBarrelCells2.AuditExecute = True + createEcalBarrelPositionedCells2.AuditExecute = True + +if doSWClustering or doTopoClustering: + TopAlg += [createemptycells] + createemptycells.AuditExecute = True + + if doSWClustering: + TopAlg += [createClusters] + createClusters.AuditExecute = True + + if addShapeParameters: + TopAlg += [augmentCaloClusters] + augmentCaloClusters.AuditExecute = True + + if doTopoClustering: + TopAlg += [createTopoClusters] + createTopoClusters.AuditExecute = True + + if addShapeParameters: + TopAlg += [augmentCaloTopoClusters] + augmentCaloTopoClusters.AuditExecute = True + +TopAlg += [ + out +] + +ApplicationMgr( + TopAlg=TopAlg, + EvtSel='NONE', + EvtMax=Nevts, + ExtSvc=ExtSvc, + StopOnSignal=True, +) + diff --git a/RecFCCeeCalorimeter/tests/options/run_thetamodulemerged.py b/RecFCCeeCalorimeter/tests/options/run_thetamodulemerged.py index 7b09cf78..66141c20 100644 --- a/RecFCCeeCalorimeter/tests/options/run_thetamodulemerged.py +++ b/RecFCCeeCalorimeter/tests/options/run_thetamodulemerged.py @@ -1,16 +1,13 @@ from Configurables import ApplicationMgr -from Configurables import EventCounter +#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 CreateCaloClustersSlidingWindow +from Configurables import CaloTowerTool from Configurables import CreateEmptyCaloCellsCollection from Configurables import CreateCaloCellPositionsFCCee from Configurables import CellPositionsECalBarrelModuleThetaSegTool -from Configurables import CellPositionsECalEndcapTurbineSegTool from Configurables import RedoSegmentation from Configurables import CreateCaloCells from Configurables import CalibrateCaloHitsTool @@ -26,36 +23,20 @@ from Configurables import HepMCToEDMConverter from Configurables import GenAlg from Configurables import FCCDataSvc -from Configurables import RewriteBitfield -from Configurables import ReadCaloCrosstalkMap -from Gaudi.Configuration import INFO +from Configurables import CaloTopoClusterInputTool +from Configurables import TopoCaloNeighbours +from Configurables import TopoCaloNoisyCells +from Configurables import CaloTopoClusterFCCee +from Gaudi.Configuration import * import os -from GaudiKernel.SystemOfUnits import GeV, tesla, mm -from GaudiKernel.PhysicalConstants import pi, halfpi, twopi -from math import cos, sin, tan +from GaudiKernel.SystemOfUnits import GeV, tesla use_pythia = False addNoise = False dumpGDML = False runHCal = False -# 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 = False -saveCells = False -saveClusterCells = False - -# clustering -doClustering = True -# NOTE: since topoclustering requires root files with noise and neighbours that are -# not in the release, topoclusters are disabled -# cluster energy corrections -# simple parametrisations of up/downstream losses -applyUpDownstreamCorrections = True -# calculate cluster energy and barycenter per layer and save it as extra parameters -addShapeParameters = True # Input for simulations (momentum is expected in GeV!) # Parameters for the particle gun simulations, dummy if use_pythia is set @@ -65,30 +46,23 @@ # (in strips: 0.5625/4=0.14) # Nevts = 20000 -Nevts = 2 +Nevts = 100 # Nevts = 1 # Nevts=1000 - -# particle momentum and direction -# momentum = 100 # in GeV -momentum = 10 # in GeV -# momentum = 10 # in GeV -thetaMin = 45 # degrees -thetaMax = 135 # degrees +# momentum = 100 # in GeV +# momentum = 50 # in GeV +momentum = 10 # in GeV +_pi = 3.14159 +thetaMin = 40 # degrees +thetaMax = 140 # degrees # thetaMin = 89 # thetaMax = 91 -# thetaMin = 90 # degrees -# thetaMax = 90 # degrees -# phiMin = halfpi -# phiMax = halfpi +# thetaMin = 90 # degrees +# thetaMax = 90 # degrees +# phiMin = _pi/2. +# phiMax = _pi/2. phiMin = 0 -phiMax = twopi - -# particle origin -# origR = 1000.0*mm -origR = 0.0 * mm -origTheta = halfpi -origPhi = 0.0 +phiMax = 2 * _pi # particle type: 11 electron, 13 muon, 22 photon, 111 pi0, 211 pi+ pdgCode = 11 @@ -121,6 +95,14 @@ pythia8gentool.printPythiaStatistics = False pythia8gentool.pythiaExtraSettings = [""] genAlg.SignalProvider = pythia8gentool + # to smear the primary vertex position: + # from Configurables import GaussSmearVertex + # smeartool = GaussSmearVertex() + # smeartool.xVertexSigma = 0.5*units.mm + # smeartool.yVertexSigma = 0.5*units.mm + # smeartool.zVertexSigma = 40.0*units.mm + # smeartool.tVertexSigma = 180.0*units.picosecond + # genAlg.VertexSmearingTool = smeartool else: from Configurables import MomentumRangeParticleGun pgun = MomentumRangeParticleGun("ParticleGun") @@ -129,30 +111,12 @@ pgun.MomentumMax = momentum * GeV pgun.PhiMin = phiMin pgun.PhiMax = phiMax - pgun.ThetaMin = thetaMin * pi / 180. - pgun.ThetaMax = thetaMax * pi / 180. + pgun.ThetaMin = thetaMin * _pi / 180. + pgun.ThetaMax = thetaMax * _pi / 180. genAlg.SignalProvider = pgun genAlg.hepmc.Path = "hepmc" -# smear/shift vertex -if origR > 0.0: - origX = origR * cos(origPhi) - origY = origR * sin(origPhi) - origZ = origR / tan(origTheta) - print("Particle gun will be moved to %f %f %f" % (origX, origY, origZ)) - from Configurables import GaussSmearVertex - vertexSmearAndShiftTool = GaussSmearVertex() - vertexSmearAndShiftTool.xVertexSigma = 0. - vertexSmearAndShiftTool.yVertexSigma = 0. - vertexSmearAndShiftTool.tVertexSigma = 0. - vertexSmearAndShiftTool.xVertexMean = origX - vertexSmearAndShiftTool.yVertexMean = origY - vertexSmearAndShiftTool.zVertexMean = origZ - vertexSmearAndShiftTool.tVertexMean = 0. - genAlg.VertexSmearingTool = vertexSmearAndShiftTool - -# hepMC -> EDM converter hepmc_converter = HepMCToEDMConverter() hepmc_converter.hepmc.Path = "hepmc" genParticlesOutputName = "genParticles" @@ -167,7 +131,8 @@ path_to_detector = os.environ.get("K4GEO", "") print(path_to_detector) detectors_to_use = [ - 'FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ALLEGRO_o1_v03.xml' + # '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 = [ @@ -210,7 +175,7 @@ if magneticField == 1: field = SimG4ConstantMagneticFieldTool( "SimG4ConstantMagneticFieldTool", - FieldComponentZ=-2 * tesla, + FieldComponentZ=-2*tesla, FieldOn=True, IntegratorStepper="ClassicalRK4" ) @@ -228,15 +193,13 @@ # ECAL ecalBarrelReadoutName = "ECalBarrelModuleThetaMerged" ecalBarrelReadoutName2 = "ECalBarrelModuleThetaMerged2" -ecalEndcapReadoutName = "ECalEndcapTurbine" +ecalEndcapReadoutName = "ECalEndcapPhiEta" # HCAL if runHCal: hcalBarrelReadoutName = "HCalBarrelReadout" - hcalBarrelReadoutName2 = "BarHCal_Readout_phitheta" hcalEndcapReadoutName = "HCalEndcapReadout" else: hcalBarrelReadoutName = "" - hcalBarrelReadoutName2 = "" hcalEndcapReadoutName = "" # Configure saving of calorimeter hits @@ -293,17 +256,13 @@ # Digitization (Merging hits into cells, EM scale calibration) # EM scale calibration (sampling fraction) calibEcalBarrel = CalibrateInLayersTool("CalibrateECalBarrel", - samplingFraction=[0.3775596654349802] * 1 + [0.13400227700041234] * 1 + [0.14390509963164044] * 1 + [0.14998482026270935] * 1 + [0.15457673722531148] * 1 + [0.15928098152159675] * 1 + [0.1635367867767212] * 1 + [0.16801070646031507] * 1 + [0.1713409944779989] * 1 + [0.17580195406064622] * 1 + [0.17966699467772812] * 1, + samplingFraction=[0.36599110182660616] * 1 + [0.1366222373338866] * 1 + [0.1452035173747207] * 1 + [0.1504319190969367] * 1 + [0.15512713637727382] * 1 + [0.1592916726494782] * 1 + [ + 0.16363478857307595] * 1 + [0.1674697333180323] * 1 + [0.16998205747422343] * 1 + [0.1739146363733975] * 1 + [0.17624609543603845] * 1 + [0.1768613530850488] * 1, readoutName=ecalBarrelReadoutName, layerFieldName="layer") -#calibEcalEndcap = CalibrateCaloHitsTool( -# "CalibrateECalEndcap", invSamplingFraction="4.27") -calibEcalEndcap = CalibrateInLayersTool("CalibrateECalEndcap", - samplingFraction = [0.16419] * 1 + [0.192898] * 1 + [0.18783] * 1 + [0.193203] * 1 + [0.193928] * 1 + [0.192286] * 1 + [0.199959] * 1 + [0.200153] * 1 + [0.212635] * 1 + [0.180345] * 1 + [0.18488] * 1 + [0.194762] * 1 + [0.197775] * 1 + [0.200504] * 1 + [0.205555] * 1 + [0.203601] * 1 + [0.210877] * 1 + [0.208376] * 1 + [0.216345] * 1 + [0.201452] * 1 + [0.202134] * 1 + [0.207566] * 1 + [0.208152] * 1 + [0.209889] * 1 + [0.211743] * 1 + [0.213188] * 1 + [0.215864] * 1 + [0.22972] * 1 + [0.192515] * 1 + [0.0103233] * 1, - readoutName=ecalEndcapReadoutName, - layerFieldName="layer") - +calibEcalEndcap = CalibrateCaloHitsTool( + "CalibrateECalEndcap", invSamplingFraction="4.27") if runHCal: calibHcells = CalibrateCaloHitsTool( "CalibrateHCal", invSamplingFraction="31.4") @@ -317,18 +276,11 @@ # (merging several modules and severla theta readout cells). # Add noise at this step if you derived the noise already assuming merged cells -# read the crosstalk map -readCrosstalkMap = ReadCaloCrosstalkMap("ReadCrosstalkMap", - fileName="https://fccsw.web.cern.ch/fccsw/filesForSimDigiReco/ALLEGRO/ALLEGRO_o1_v03/xtalk_neighbours_map_ecalB_thetamodulemerged.root", - OutputLevel=INFO) - # Step 1: merge hits into cells according to initial segmentation ecalBarrelCellsName = "ECalBarrelCells" createEcalBarrelCells = CreateCaloCells("CreateECalBarrelCells", doCellCalibration=True, calibTool=calibEcalBarrel, - crosstalksTool=readCrosstalkMap, - addCrosstalk=False, addCellNoise=False, filterCellNoise=False, addPosition=True, @@ -394,35 +346,18 @@ # Create cells in ECal endcap -ecalEndcapCellsName = "ECalEndcapCells" createEcalEndcapCells = CreateCaloCells("CreateEcalEndcapCaloCells", doCellCalibration=True, calibTool=calibEcalEndcap, addCellNoise=False, filterCellNoise=False, - OutputLevel=INFO, - hits=ecalEndcapReadoutName, - cells=ecalEndcapCellsName) - -# Add to Ecal endcap cells the position information -# (good for physics, all coordinates set properly) -cellPositionEcalEndcapTool = CellPositionsECalEndcapTurbineSegTool( - "CellPositionsECalEndcap", - readoutName=ecalEndcapReadoutName, - OutputLevel=INFO -) -ecalEndcapPositionedCellsName = "ECalEndcapPositionedCells" -createEcalEndcapPositionedCells = CreateCaloCellPositionsFCCee( - "CreateECalEndcapPositionedCells", - OutputLevel=INFO -) -createEcalEndcapPositionedCells.positionsTool = cellPositionEcalEndcapTool -createEcalEndcapPositionedCells.hits.Path = ecalEndcapCellsName -createEcalEndcapPositionedCells.positionedHits.Path = ecalEndcapPositionedCellsName + OutputLevel=INFO) +createEcalEndcapCells.hits.Path = "ECalEndcapHits" +createEcalEndcapCells.cells.Path = "ECalEndcapCells" if runHCal: # Create cells in HCal - # 1 - merge hits into cells with the default readout + # 1. step - merge hits into cells with the default readout hcalBarrelCellsName = "HCalBarrelCells" createHcalBarrelCells = CreateCaloCells("CreateHCalBarrelCells", doCellCalibration=True, @@ -434,14 +369,22 @@ cells=hcalBarrelCellsName, OutputLevel=INFO) - # 2 - attach positions to the cells + # createHcalEndcapCells = CreateCaloCells("CreateHcalEndcapCaloCells", + # doCellCalibration=True, + # calibTool=calibHcalEndcap, + # addCellNoise=False, + # filterCellNoise=False, + # OutputLevel=INFO) + # createHcalEndcapCells.hits.Path="HCalEndcapHits" + # createHcalEndcapCells.cells.Path="HCalEndcapCells" + from Configurables import CellPositionsHCalBarrelPhiThetaSegTool cellPositionHcalBarrelTool = CellPositionsHCalBarrelPhiThetaSegTool( "CellPositionsHCalBarrel", readoutName=hcalBarrelReadoutName, OutputLevel=INFO ) - hcalBarrelPositionedCellsName = "HCalBarrelPositionedCells" + hcalBarrelPositionedCellsName = "HCalBarrelPositionedCells" createHcalBarrelPositionedCells = CreateCaloCellPositionsFCCee( "CreateHcalBarrelPositionedCells", OutputLevel=INFO @@ -449,158 +392,196 @@ createHcalBarrelPositionedCells.positionsTool = cellPositionHcalBarrelTool createHcalBarrelPositionedCells.hits.Path = hcalBarrelCellsName createHcalBarrelPositionedCells.positionedHits.Path = hcalBarrelPositionedCellsName - - # 3 - compute new cellID of cells based on new readout - removing row information - hcalBarrelCellsName2 = "HCalBarrelCells2" - rewriteHCalBarrel = RewriteBitfield("RewriteHCalBarrel", - # old bitfield (readout) - oldReadoutName=hcalBarrelReadoutName, - # specify which fields are going to be deleted - removeIds=["row"], - # new bitfield (readout), with new segmentation - newReadoutName=hcalBarrelReadoutName2, - debugPrint=10, - OutputLevel=INFO) - # clusters are needed, with deposit position and cellID in bits - rewriteHCalBarrel.inhits.Path = hcalBarrelCellsName - rewriteHCalBarrel.outhits.Path = hcalBarrelCellsName2 - - # 4 - attach positions to the new cells - from Configurables import CellPositionsHCalBarrelPhiThetaSegTool - hcalBarrelPositionedCellsName2 = "HCalBarrelPositionedCells2" - cellPositionHcalBarrelTool2 = CellPositionsHCalBarrelPhiThetaSegTool( - "CellPositionsHCalBarrel2", - readoutName=hcalBarrelReadoutName2, - OutputLevel=INFO - ) - createHcalBarrelPositionedCells2 = CreateCaloCellPositionsFCCee( - "CreateHCalBarrelPositionedCells2", - OutputLevel=INFO - ) - createHcalBarrelPositionedCells2.positionsTool = cellPositionHcalBarrelTool2 - createHcalBarrelPositionedCells2.hits.Path = hcalBarrelCellsName2 - createHcalBarrelPositionedCells2.positionedHits.Path = hcalBarrelPositionedCellsName2 - - # createHcalEndcapCells = CreateCaloCells("CreateHcalEndcapCaloCells", - # doCellCalibration=True, - # calibTool=calibHcalEndcap, - # addCellNoise=False, - # filterCellNoise=False, - # OutputLevel=INFO) - # createHcalEndcapCells.hits.Path="HCalEndcapHits" - # createHcalEndcapCells.cells.Path="HCalEndcapCells" - else: hcalBarrelCellsName = "emptyCaloCells" hcalBarrelPositionedCellsName = "emptyCaloCells" - hcalBarrelCellsName2 = "emptyCaloCells" - hcalBarrelPositionedCellsName2 = "emptyCaloCells" cellPositionHcalBarrelTool = None - cellPositionHcalBarrelTool2 = None - + # Empty cells for parts of calorimeter not implemented yet createemptycells = CreateEmptyCaloCellsCollection("CreateEmptyCaloCells") createemptycells.cells.Path = "emptyCaloCells" -# Produce sliding window clusters (ECAL only) -towers = CaloTowerToolFCCee("towers", - deltaThetaTower=4 * 0.009817477/4, deltaPhiTower=2 * 2 * pi / 1536., - ecalBarrelReadoutName=ecalBarrelReadoutName, - ecalEndcapReadoutName=ecalEndcapReadoutName, - ecalFwdReadoutName="", - hcalBarrelReadoutName="", - hcalExtBarrelReadoutName="", - hcalEndcapReadoutName="", - hcalFwdReadoutName="", - OutputLevel=INFO) -towers.ecalBarrelCells.Path = ecalBarrelPositionedCellsName -towers.ecalEndcapCells.Path = ecalEndcapPositionedCellsName +# Produce sliding window clusters +towers = CaloTowerTool("towers", + deltaEtaTower=0.01, deltaPhiTower=2*_pi/768, + radiusForPosition=2160 + 40 / 2.0, + # ecalBarrelReadoutName = ecalBarrelReadoutNamePhiEta, + ecalBarrelReadoutName=ecalBarrelReadoutName, + ecalEndcapReadoutName=ecalEndcapReadoutName, + ecalFwdReadoutName="", + hcalBarrelReadoutName=hcalBarrelReadoutName, + hcalExtBarrelReadoutName="", + hcalEndcapReadoutName="", + hcalFwdReadoutName="", + OutputLevel=INFO) +towers.ecalBarrelCells.Path = ecalBarrelCellsName +towers.ecalEndcapCells.Path = "ECalEndcapCells" towers.ecalFwdCells.Path = "emptyCaloCells" -towers.hcalBarrelCells.Path = "emptyCaloCells" +towers.hcalBarrelCells.Path = hcalBarrelCellsName towers.hcalExtBarrelCells.Path = "emptyCaloCells" towers.hcalEndcapCells.Path = "emptyCaloCells" towers.hcalFwdCells.Path = "emptyCaloCells" # Cluster variables -windT = 9 +windE = 9 windP = 17 -posT = 5 +posE = 5 posP = 11 -dupT = 7 +dupE = 7 dupP = 13 -finT = 9 +finE = 9 finP = 17 # Minimal energy to create a cluster in GeV (FCC-ee detectors have to reconstruct low energy particles) threshold = 0.040 -createClusters = CreateCaloClustersSlidingWindowFCCee("CreateClusters", - towerTool=towers, - nThetaWindow=windT, nPhiWindow=windP, - nThetaPosition=posT, nPhiPosition=posP, - nThetaDuplicates=dupT, nPhiDuplicates=dupP, - nThetaFinal=finT, nPhiFinal=finP, - energyThreshold=threshold, - energySharingCorrection=False, - attachCells=True, - OutputLevel=INFO - ) +createClusters = CreateCaloClustersSlidingWindow("CreateClusters", + towerTool=towers, + nEtaWindow=windE, nPhiWindow=windP, + nEtaPosition=posE, nPhiPosition=posP, + nEtaDuplicates=dupE, nPhiDuplicates=dupP, + nEtaFinal=finE, nPhiFinal=finP, + energyThreshold=threshold, + energySharingCorrection=False, + attachCells=True, + OutputLevel=INFO + ) createClusters.clusters.Path = "CaloClusters" createClusters.clusterCells.Path = "CaloClusterCells" +createEcalBarrelPositionedCaloClusterCells = CreateCaloCellPositionsFCCee( + "ECalBarrelPositionedCaloClusterCells", + OutputLevel=INFO +) +createEcalBarrelPositionedCaloClusterCells.positionsTool = cellPositionEcalBarrelTool +createEcalBarrelPositionedCaloClusterCells.hits.Path = "CaloClusterCells" +createEcalBarrelPositionedCaloClusterCells.positionedHits.Path = "PositionedCaloClusterCells" + correctCaloClusters = CorrectCaloClusters("correctCaloClusters", inClusters=createClusters.clusters.Path, - outClusters="Corrected" + createClusters.clusters.Path, - systemIDs=[4], - numLayers=[11], + outClusters="Corrected"+createClusters.clusters.Path, + numLayers=[0], firstLayerIDs=[0], - lastLayerIDs=[10], + lastLayerIDs=[-1], + # readoutNames = [ecalBarrelReadoutNamePhiEta], 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]], + # upstreamParameters = [[0.02729094887360858, -1.378665489864182, -68.40424543618059, 3.6930827214130053, -5528.714729126099, -1630.7911298009794]], + upstreamParameters=[ + [0.02729094887360858, -1.378665489864182, -68.40424543618059, 3.6930827214130053, -5528.714729126099, -1630.7911298009794]], upstreamFormulas=[ ['[0]+[1]/(x-[2])', '[0]+[1]/(x-[2])']], - # do not split the following line or it will break scripts that update the values of the corrections - downstreamParameters = [[0.0018280333929494054, 0.004932212590963076, 0.8409676097173655, -1.2676690014715288, 0.005347798049886769, 4.161741293789687]], + # downstreamParameters = [[-0.0032351643028483354, 0.006597484738888312, 0.8972024981692965, -1.0207168610322181, 0.017878133854084398, 9.108099243443101]], + downstreamParameters=[ + [-0.0032351643028483354, 0.006597484738888312, 0.8972024981692965, -1.0207168610322181, 0.017878133854084398, 9.108099243443101]], downstreamFormulas=[ ['[0]+[1]*x', '[0]+[1]/sqrt(x)', '[0]+[1]/x']], OutputLevel=INFO ) -augmentCaloClusters = AugmentClustersFCCee("augmentCaloClusters", - inClusters=createClusters.clusters.Path, - outClusters="Augmented" + createClusters.clusters.Path, - systemIDs=[4], - systemNames=["EMB"], - numLayers=[11], - readoutNames=[ecalBarrelReadoutName], - layerFieldNames=["layer"], - thetaRecalcWeights=[ - # [-1, 3.0, 3.0, 3.0, 4.25, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0] # -1 : use linear weights - [-1, 3.0, 3.0, 3.0, 4.25, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0] # -1 : use linear weights - ], - OutputLevel=INFO - ) +# TOPO CLUSTERS PRODUCTION +createTopoInput = CaloTopoClusterInputTool("CreateTopoInput", + ecalBarrelReadoutName=ecalBarrelReadoutName, + ecalEndcapReadoutName="", + ecalFwdReadoutName="", + hcalBarrelReadoutName=hcalBarrelReadoutName, + hcalExtBarrelReadoutName="", + hcalEndcapReadoutName="", + hcalFwdReadoutName="", + OutputLevel=INFO) + +createTopoInput.ecalBarrelCells.Path = ecalBarrelPositionedCellsName +# createTopoInput.ecalBarrelCells.Path = "ECalBarrelPositionedCells2" +createTopoInput.ecalEndcapCells.Path = "emptyCaloCells" +createTopoInput.ecalFwdCells.Path = "emptyCaloCells" +createTopoInput.hcalBarrelCells.Path = hcalBarrelPositionedCellsName +createTopoInput.hcalExtBarrelCells.Path = "emptyCaloCells" +createTopoInput.hcalEndcapCells.Path = "emptyCaloCells" +createTopoInput.hcalFwdCells.Path = "emptyCaloCells" +cellPositionHcalBarrelNoSegTool = None +cellPositionHcalExtBarrelTool = None + +readNeighboursMap = TopoCaloNeighbours("ReadNeighboursMap", + fileName=os.environ['FCCBASEDIR'] + + "/LAr_scripts/data/neighbours_map_barrel_thetamodulemerged.root", + #"/LAr_scripts/data/neighbours_map_HCalBarrel.root", + OutputLevel=INFO) + +# Noise levels per cell +readNoisyCellsMap = TopoCaloNoisyCells("ReadNoisyCellsMap", + fileName=os.environ['FCCBASEDIR'] + + "/LAr_scripts/data/cellNoise_map_electronicsNoiseLevel_thetamodulemerged.root", + OutputLevel=INFO) +createTopoClusters = CaloTopoClusterFCCee("CreateTopoClusters", + TopoClusterInput=createTopoInput, + # expects neighbours map from cellid->vec < neighbourIds > + neigboursTool=readNeighboursMap, + # tool to get noise level per cellid + noiseTool=readNoisyCellsMap, + # cell positions tools for all sub - systems + positionsECalBarrelTool=cellPositionEcalBarrelTool, + positionsHCalBarrelTool=cellPositionHcalBarrelTool, + positionsHCalBarrelNoSegTool=cellPositionHcalBarrelNoSegTool, + positionsHCalExtBarrelTool=cellPositionHcalExtBarrelTool, + # positionsHCalExtBarrelTool = HCalExtBcells, + # positionsEMECTool = EMECcells, + # positionsHECTool = HECcells, + # positionsEMFwdTool = ECalFwdcells, + # positionsHFwdTool = HCalFwdcells, + noSegmentationHCal=False, + seedSigma=4, + neighbourSigma=2, + lastNeighbourSigma=0, + OutputLevel=INFO) +createTopoClusters.clusters.Path = "CaloTopoClusters" +createTopoClusters.clusterCells.Path = "CaloTopoClusterCells" + +createEcalBarrelPositionedCaloTopoClusterCells = CreateCaloCellPositionsFCCee( + "ECalBarrelPositionedCaloTopoClusterCells", + OutputLevel=INFO +) +# createEcalBarrelPositionedCaloTopoClusterCells.positionsTool = cellPositionEcalBarrelTool2 +createEcalBarrelPositionedCaloTopoClusterCells.positionsTool = cellPositionEcalBarrelTool +createEcalBarrelPositionedCaloTopoClusterCells.hits.Path = "CaloTopoClusterCells" +createEcalBarrelPositionedCaloTopoClusterCells.positionedHits.Path = "PositionedCaloTopoClusterCells" + +correctCaloTopoClusters = CorrectCaloClusters( + "correctCaloTopoClusters", + inClusters=createTopoClusters.clusters.Path, + outClusters="Corrected"+createTopoClusters.clusters.Path, + numLayers=[0], + firstLayerIDs=[0], + lastLayerIDs=[-1], + # readoutNames = [ecalBarrelReadoutNamePhiEta], + readoutNames=[ecalBarrelReadoutName], + # upstreamParameters = [[0.02729094887360858, -1.378665489864182, -68.40424543618059, 3.6930827214130053, -5528.714729126099, -1630.7911298009794]], + upstreamParameters=[[0.02729094887360858, -1.378665489864182, -68.40424543618059, + 3.6930827214130053, -5528.714729126099, -1630.7911298009794]], + upstreamFormulas=[['[0]+[1]/(x-[2])', '[0]+[1]/(x-[2])']], + # downstreamParameters = [[-0.0032351643028483354, 0.006597484738888312, 0.8972024981692965, -1.0207168610322181, 0.017878133854084398, 9.108099243443101]], + downstreamParameters=[[-0.0032351643028483354, 0.006597484738888312, + 0.8972024981692965, -1.0207168610322181, 0.017878133854084398, 9.108099243443101]], + downstreamFormulas=[['[0]+[1]*x', '[0]+[1]/sqrt(x)', '[0]+[1]/x']], + OutputLevel=INFO +) # Output out = PodioOutput("out", OutputLevel=INFO) +# out.outputCommands = ["keep *"] +# out.outputCommands = ["keep *", "drop ECalBarrelHits", "drop HCal*", "drop ECalBarrelCellsStep*", "drop ECalBarrelPositionedHits", "drop emptyCaloCells", "drop CaloClusterCells"] +# out.outputCommands = ["keep *", "drop ECalBarrelHits", "drop HCal*", "drop ECalBarrelCellsStep*", "drop ECalBarrelPositionedHits", "drop emptyCaloCells", "drop CaloClusterCells", "drop %s"%ecalBarrelCellsName, "drop %s"%createEcalBarrelPositionedCells.positionedHits.Path] +# out.outputCommands = ["keep *", "drop ECalBarrelHits", "drop HCal*", "drop *ells*", "drop ECalBarrelPositionedHits", "drop emptyCaloCells"] +# out.outputCommands = ["keep *", "drop HCal*", "drop ECalBarrel*", "drop emptyCaloCells"] if runHCal: out.outputCommands = ["keep *", "drop emptyCaloCells"] else: out.outputCommands = ["keep *", "drop HCal*", "drop emptyCaloCells"] -if not saveCells: - out.outputCommands.append("drop ECal*Cells*") -if not saveClusterCells: - out.outputCommands.append("drop *ClusterCells*") -if not saveHits: - out.outputCommands.append("drop ECal*Hits*") - -out.filename = "./output_evts_" + str(Nevts) + "_pdg_" + str(pdgCode) + "_" + str(momentum) + "_GeV" + "_ThetaMinMax_" + str(thetaMin) + "_" + str( - thetaMax) + "_PhiMinMax_" + str(phiMin) + "_" + str(phiMax) + "_MagneticField_" + str(magneticField) + "_Noise" + str(addNoise) + ".root" +# out.filename = "root/output_fullCalo_SimAndDigi_withTopoCluster_MagneticField_"+str(magneticField)+"_pMin_"+str(momentum*1000)+"_MeV"+"_ThetaMinMax_"+str(thetaMin)+"_"+str(thetaMax)+"_pdgId_"+str(pdgCode)+"_pythia"+str(use_pythia)+"_Noise"+str(addNoise)+".root" +out.filename = "./root_merge/output_evts_"+str(Nevts)+"_pdg_"+str(pdgCode)+"_"+str(momentum)+"_GeV"+"_ThetaMinMax_"+str(thetaMin)+"_"+str( + thetaMax)+"_PhiMinMax_"+str(phiMin)+"_"+str(phiMax)+"_MagneticField_"+str(magneticField)+"_Noise"+str(addNoise)+".root" # CPU information chra = ChronoAuditor() @@ -613,17 +594,18 @@ createEcalBarrelPositionedCells.AuditExecute = True if runHCal: createHcalBarrelCells.AuditExecute = True +createTopoClusters.AuditExecute = True out.AuditExecute = True -event_counter = EventCounter('event_counter') -event_counter.Frequency = 10 +#event_counter = EventCounter('event_counter') +#event_counter.Frequency = 10 ExtSvc = [geoservice, podioevent, geantservice, audsvc] if dumpGDML: ExtSvc += [gdmldumpservice] TopAlg = [ - event_counter, +# event_counter, genAlg, hepmc_converter, geantsim, @@ -632,36 +614,22 @@ resegmentEcalBarrel, createEcalBarrelCells2, createEcalBarrelPositionedCells2, - createEcalEndcapCells, - createEcalEndcapPositionedCells, + createEcalEndcapCells ] - if runHCal: TopAlg += [ createHcalBarrelCells, createHcalBarrelPositionedCells, - rewriteHCalBarrel, - createHcalBarrelPositionedCells2, # createHcalEndcapCells ] - -if doClustering: - TopAlg += [ - createemptycells, - createClusters, - ] - - if applyUpDownstreamCorrections: - TopAlg += [ - correctCaloClusters, - ] - - if addShapeParameters: - TopAlg += [ - augmentCaloClusters, - ] - TopAlg += [ + createemptycells, + # createClusters, + # createEcalBarrelPositionedCaloClusterCells, + # correctCaloClusters, + createTopoClusters, + createEcalBarrelPositionedCaloTopoClusterCells, + correctCaloTopoClusters, out ] diff --git a/RecFCCeeCalorimeter/tests/options/run_thetamodulemerged_SW_benchmarkCorr.py b/RecFCCeeCalorimeter/tests/options/run_thetamodulemerged_SW_benchmarkCorr.py index 45186ee7..b12d58d1 100644 --- a/RecFCCeeCalorimeter/tests/options/run_thetamodulemerged_SW_benchmarkCorr.py +++ b/RecFCCeeCalorimeter/tests/options/run_thetamodulemerged_SW_benchmarkCorr.py @@ -1,676 +1,678 @@ -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 CreateEmptyCaloCellsCollection -from Configurables import CreateCaloCellPositionsFCCee -from Configurables import CellPositionsECalBarrelModuleThetaSegTool -from Configurables import CellPositionsECalEndcapTurbineSegTool -from Configurables import RedoSegmentation -from Configurables import CreateCaloCells -from Configurables import CalibrateCaloHitsTool -from Configurables import CalibrateInLayersTool -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 -from Configurables import GeoSvc -from Configurables import HepMCToEDMConverter -from Configurables import GenAlg -from Configurables import FCCDataSvc -from Configurables import ReadCaloCrosstalkMap -from Gaudi.Configuration import INFO -# from Gaudi.Configuration import * - -import os - -from GaudiKernel.SystemOfUnits import GeV, tesla, mm -from GaudiKernel.PhysicalConstants import pi, halfpi, twopi -from math import cos, sin, tan - -use_pythia = False -addNoise = False -dumpGDML = False -runHCal = True -# 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 = False -saveCells = False -saveClusterCells = False - -# cluster energy corrections -# simple parametrisations of up/downstream losses -applyUpDownstreamBenchmarkCorrections = True - -# Input for simulations (momentum is expected in GeV!) -# Parameters for the particle gun simulations, dummy if use_pythia is set -# to True -# theta from 80 to 100 degrees corresponds to -0.17 < eta < 0.17 -# reminder: cell granularity in theta = 0.5625 degrees -# (in strips: 0.5625/4=0.14) - -# Nevts = 20000 -Nevts = 2 -# Nevts = 1 -# Nevts=1000 - -# particle momentum and direction -# momentum = 100 # in GeV -momentum = 10. # in GeV -# momentum = 10 # in GeV -thetaMin = 69. # degrees -thetaMax = 69. # degrees -# thetaMin = 89 -# thetaMax = 91 -# thetaMin = 90 # degrees -# thetaMax = 90 # degrees -# phiMin = halfpi -# phiMax = halfpi -phiMin = 0 -phiMax = twopi - -# particle origin -# origR = 1000.0*mm -origR = 0.0 * mm -origTheta = halfpi -origPhi = 0.0 - -# particle type: 11 electron, 13 muon, 22 photon, 111 pi0, 211 pi+ -pdgCode = 211 -# pdgCode = 22 -# pdgCode = 111 -# pdgCode = 211 - -# Set to true if history from Geant4 decays is needed (e.g. to get the -# photons from pi0) -saveG4Hist = False -if (pdgCode == 111): - saveG4Hist = True - -magneticField = False - - -podioevent = FCCDataSvc("EventDataSvc") - -# Particle gun setup - -genAlg = GenAlg() -if use_pythia: - from Configurables import PythiaInterface - pythia8gentool = PythiaInterface() - pythia8gentool.pythiacard = os.path.join( - os.environ.get('PWD', ''), - "MCGeneration/ee_Zgamma_inclusive.cmd" - ) - # pythia8gentool.pythiacard = "MCGeneration/ee_Z_ee.cmd" - pythia8gentool.printPythiaStatistics = False - pythia8gentool.pythiaExtraSettings = [""] - genAlg.SignalProvider = pythia8gentool -else: - from Configurables import MomentumRangeParticleGun - pgun = MomentumRangeParticleGun("ParticleGun") - pgun.PdgCodes = [pdgCode] - pgun.MomentumMin = momentum * GeV - pgun.MomentumMax = momentum * GeV - pgun.PhiMin = phiMin - pgun.PhiMax = phiMax - pgun.ThetaMin = thetaMin * pi / 180. - pgun.ThetaMax = thetaMax * pi / 180. - genAlg.SignalProvider = pgun - -genAlg.hepmc.Path = "hepmc" - -# smear/shift vertex -if origR > 0.0: - origX = origR * cos(origPhi) - origY = origR * sin(origPhi) - origZ = origR / tan(origTheta) - print("Particle gun will be moved to %f %f %f" % (origX, origY, origZ)) - from Configurables import GaussSmearVertex - vertexSmearAndShiftTool = GaussSmearVertex() - vertexSmearAndShiftTool.xVertexSigma = 0. - vertexSmearAndShiftTool.yVertexSigma = 0. - vertexSmearAndShiftTool.tVertexSigma = 0. - vertexSmearAndShiftTool.xVertexMean = origX - vertexSmearAndShiftTool.yVertexMean = origY - vertexSmearAndShiftTool.zVertexMean = origZ - vertexSmearAndShiftTool.tVertexMean = 0. - genAlg.VertexSmearingTool = vertexSmearAndShiftTool - -# hepMC -> EDM converter -hepmc_converter = HepMCToEDMConverter() -hepmc_converter.hepmc.Path = "hepmc" -genParticlesOutputName = "genParticles" -hepmc_converter.GenParticles.Path = genParticlesOutputName -hepmc_converter.hepmcStatusList = [] -hepmc_converter.OutputLevel = INFO - -# Simulation setup -# Detector geometry -geoservice = GeoSvc("GeoSvc") -# 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 = [ - 'FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ALLEGRO_o1_v03.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 -actions = SimG4FullSimActions() - -if saveG4Hist: - actions.enableHistory = True - actions.energyCut = 1.0 * GeV - saveHistTool = SimG4SaveParticleHistory("saveHistory") - -geantservice = SimG4Svc( - "SimG4Svc", - detector='SimG4DD4hepDetector', - physicslist="SimG4FtfpBert", - actions=actions -) - -# from Configurables import GeoToGdmlDumpSvc -if dumpGDML: - from Configurables import GeoToGdmlDumpSvc - gdmldumpservice = GeoToGdmlDumpSvc("GeoToGdmlDumpSvc") - -# Fixed seed to have reproducible results, change it for each job if you -# split one production into several jobs -# Mind that if you leave Gaudi handle random seed and some job start within -# the same second (very likely) you will have duplicates -geantservice.randomNumbersFromGaudi = False -geantservice.seedValue = 4242 - -# Range cut -geantservice.g4PreInitCommands += ["/run/setCut 0.1 mm"] - -# Magnetic field -if magneticField == 1: - field = SimG4ConstantMagneticFieldTool( - "SimG4ConstantMagneticFieldTool", - FieldComponentZ=-2 * tesla, - FieldOn=True, - IntegratorStepper="ClassicalRK4" - ) -else: - field = SimG4ConstantMagneticFieldTool( - "SimG4ConstantMagneticFieldTool", - FieldOn=False - ) - -# Geant4 algorithm -# Translates EDM to G4Event, passes the event to G4, writes out outputs -# via tools and a tool that saves the calorimeter hits - -# Detector readouts -# ECAL -ecalBarrelReadoutName = "ECalBarrelModuleThetaMerged" -ecalBarrelReadoutName2 = "ECalBarrelModuleThetaMerged2" -ecalEndcapReadoutName = "ECalEndcapTurbine" -# HCAL -if runHCal: - hcalBarrelReadoutName = "HCalBarrelReadout" - hcalBarrelReadoutName2 = "BarHCal_Readout_phitheta" - hcalEndcapReadoutName = "HCalEndcapReadout" -else: - hcalBarrelReadoutName = "" - hcalBarrelReadoutName2 = "" - hcalEndcapReadoutName = "" - -# Configure saving of calorimeter hits -ecalBarrelHitsName = "ECalBarrelPositionedHits" -saveECalBarrelTool = SimG4SaveCalHits( - "saveECalBarrelHits", - readoutName=ecalBarrelReadoutName, - OutputLevel=INFO -) -saveECalBarrelTool.CaloHits.Path = ecalBarrelHitsName - -saveECalEndcapTool = SimG4SaveCalHits( - "saveECalEndcapHits", - readoutName=ecalEndcapReadoutName -) -saveECalEndcapTool.CaloHits.Path = "ECalEndcapHits" - -if runHCal: - hcalBarrelHitsName = "HCalBarrelPositionedHits" - saveHCalTool = SimG4SaveCalHits( - "saveHCalBarrelHits", - readoutName=hcalBarrelReadoutName - ) - saveHCalTool.CaloHits.Path = hcalBarrelHitsName - - # saveHCalEndcapTool = SimG4SaveCalHits( - # "saveHCalEndcapHits", - # readoutName = hcalEndcapReadoutName - # ) - # saveHCalEndcapTool.CaloHits.Path = "HCalEndcapHits" - -# next, create the G4 algorithm, giving the list of names of tools ("XX/YY") -particle_converter = SimG4PrimariesFromEdmTool("EdmConverter") -particle_converter.GenParticles.Path = genParticlesOutputName - -outputTools = [ - saveECalBarrelTool, - saveECalEndcapTool -] -if runHCal: - outputTools += [ - saveHCalTool, - # saveHCalEndcapTool - ] - -if saveG4Hist: - outputTools += [saveHistTool] - -geantsim = SimG4Alg("SimG4Alg", - outputs=outputTools, - eventProvider=particle_converter, - OutputLevel=INFO) - -# Digitization (Merging hits into cells, EM scale calibration) -# EM scale calibration (sampling fraction) -calibEcalBarrel = CalibrateInLayersTool("CalibrateECalBarrel", - samplingFraction=[0.3864252122990472] * 1 + [0.13597644835735828] * 1 + [0.14520427829645913] * 1 + [0.1510076084632846] * 1 + [0.1552347580991012] * 1 + [0.159694330729184] * 1 + [0.1632954482794191] * 1 + [0.16720711037339814] * 1 + [0.17047749048884808] * 1 + [0.17461698117974286] * 1 + [0.1798984163980135] * 1 + [0.17920355117405806] * 1, - readoutName=ecalBarrelReadoutName, - layerFieldName="layer") - -#calibEcalEndcap = CalibrateCaloHitsTool( -# "CalibrateECalEndcap", invSamplingFraction="4.27") -calibEcalEndcap = CalibrateInLayersTool("CalibrateECalEndcap", - samplingFraction = [0.16419] * 1 + [0.192898] * 1 + [0.18783] * 1 + [0.193203] * 1 + [0.193928] * 1 + [0.192286] * 1 + [0.199959] * 1 + [0.200153] * 1 + [0.212635] * 1 + [0.180345] * 1 + [0.18488] * 1 + [0.194762] * 1 + [0.197775] * 1 + [0.200504] * 1 + [0.205555] * 1 + [0.203601] * 1 + [0.210877] * 1 + [0.208376] * 1 + [0.216345] * 1 + [0.201452] * 1 + [0.202134] * 1 + [0.207566] * 1 + [0.208152] * 1 + [0.209889] * 1 + [0.211743] * 1 + [0.213188] * 1 + [0.215864] * 1 + [0.22972] * 1 + [0.192515] * 1 + [0.0103233] * 1, - readoutName=ecalEndcapReadoutName, - layerFieldName="layer") - - -if runHCal: - calibHcells = CalibrateCaloHitsTool( - "CalibrateHCal", invSamplingFraction="30.4") - calibHcalEndcap = CalibrateCaloHitsTool( - "CalibrateHCalEndcap", invSamplingFraction="31.7") - -# Create cells in ECal barrel -# 1. step - merge hits into cells with theta and module segmentation -# (module is a 'physical' cell i.e. lead + LAr + PCB + LAr +lead) -# 2. step - rewrite the cellId using the merged theta-module segmentation -# (merging several modules and severla theta readout cells). -# Add noise at this step if you derived the noise already assuming merged cells - -# read the crosstalk map -readCrosstalkMap = ReadCaloCrosstalkMap("ReadCrosstalkMap", - fileName="https://fccsw.web.cern.ch/fccsw/filesForSimDigiReco/ALLEGRO/ALLEGRO_o1_v03/xtalk_neighbours_map_ecalB_thetamodulemerged.root", - OutputLevel=INFO) - -# Step 1: merge hits into cells according to initial segmentation -ecalBarrelCellsName = "ECalBarrelCells" -createEcalBarrelCells = CreateCaloCells("CreateECalBarrelCells", - doCellCalibration=True, - calibTool=calibEcalBarrel, - crosstalksTool=readCrosstalkMap, - addCrosstalk=False, - addCellNoise=False, - filterCellNoise=False, - addPosition=True, - OutputLevel=INFO, - hits=ecalBarrelHitsName, - cells=ecalBarrelCellsName) - -# Step 2a: compute new cellID of cells based on new readout -# (merged module-theta segmentation with variable merging vs layer) -resegmentEcalBarrel = RedoSegmentation("ReSegmentationEcal", - # old bitfield (readout) - oldReadoutName=ecalBarrelReadoutName, - # specify which fields are going to be altered (deleted/rewritten) - oldSegmentationIds=["module", "theta"], - # new bitfield (readout), with new segmentation (merged modules and theta cells) - newReadoutName=ecalBarrelReadoutName2, - OutputLevel=INFO, - debugPrint=200, - inhits=ecalBarrelCellsName, - outhits="ECalBarrelCellsMerged") - -# Step 2b: merge new cells with same cellID together -# do not apply cell calibration again since cells were already -# calibrated in Step 1 -ecalBarrelCellsName2 = "ECalBarrelCells2" -createEcalBarrelCells2 = CreateCaloCells("CreateECalBarrelCells2", - doCellCalibration=False, - addCellNoise=False, - filterCellNoise=False, - OutputLevel=INFO, - hits="ECalBarrelCellsMerged", - cells=ecalBarrelCellsName2) - -# Add to Ecal barrel cells the position information -# (good for physics, all coordinates set properly) - -cellPositionEcalBarrelTool = CellPositionsECalBarrelModuleThetaSegTool( - "CellPositionsECalBarrel", - readoutName=ecalBarrelReadoutName, - OutputLevel=INFO -) -ecalBarrelPositionedCellsName = "ECalBarrelPositionedCells" -createEcalBarrelPositionedCells = CreateCaloCellPositionsFCCee( - "CreateECalBarrelPositionedCells", - OutputLevel=INFO -) -createEcalBarrelPositionedCells.positionsTool = cellPositionEcalBarrelTool -createEcalBarrelPositionedCells.hits.Path = ecalBarrelCellsName -createEcalBarrelPositionedCells.positionedHits.Path = ecalBarrelPositionedCellsName - -cellPositionEcalBarrelTool2 = CellPositionsECalBarrelModuleThetaSegTool( - "CellPositionsECalBarrel2", - readoutName=ecalBarrelReadoutName2, - OutputLevel=INFO -) -createEcalBarrelPositionedCells2 = CreateCaloCellPositionsFCCee( - "CreateECalBarrelPositionedCells2", - OutputLevel=INFO -) -createEcalBarrelPositionedCells2.positionsTool = cellPositionEcalBarrelTool2 -createEcalBarrelPositionedCells2.hits.Path = ecalBarrelCellsName2 -createEcalBarrelPositionedCells2.positionedHits.Path = "ECalBarrelPositionedCells2" - - -# Create cells in ECal endcap -createEcalEndcapCells = CreateCaloCells("CreateEcalEndcapCaloCells", - doCellCalibration=True, - calibTool=calibEcalEndcap, - addCellNoise=False, - filterCellNoise=False, - OutputLevel=INFO) -createEcalEndcapCells.hits.Path = "ECalEndcapHits" -createEcalEndcapCells.cells.Path = "ECalEndcapCells" - -# Add to Ecal endcap cells the position information -# (good for physics, all coordinates set properly) -# not yet merged!! - cellPositionEcalEndcapTool = CellPositionsECalEndcapTurbineSegTool( - "CellPositionsECalEndcap", - readoutName=ecalEndcapReadoutName, - OutputLevel=INFO -) -ecalEndcapPositionedCellsName = "ECalEndcapPositionedCells" -createEcalEndcapPositionedCells = CreateCaloCellPositionsFCCee( - "CreateECalEndcapPositionedCells", - OutputLevel=INFO -) -createEcalEndcapPositionedCells.positionsTool = cellPositionEcalEndcapTool -createEcalEndcapPositionedCells.hits.Path = ecalEndcapCellsName -createEcalEndcapPositionedCells.positionedHits.Path = ecalEndcapPositionedCellsName - -if runHCal: - # Create cells in HCal - # 1 - merge hits into cells with the default readout - hcalBarrelCellsName = "HCalBarrelCells" - createHcalBarrelCells = CreateCaloCells("CreateHCalBarrelCells", - doCellCalibration=True, - calibTool=calibHcells, - addCellNoise=False, - filterCellNoise=False, - addPosition=True, - hits=hcalBarrelHitsName, - cells=hcalBarrelCellsName, - OutputLevel=INFO) - - # 2 - attach positions to the cells (cell positions needed for RedoSegmentation!) - from Configurables import CellPositionsHCalBarrelPhiThetaSegTool - cellPositionHcalBarrelTool = CellPositionsHCalBarrelPhiThetaSegTool( - "CellPositionsHCalBarrel", - readoutName=hcalBarrelReadoutName, - OutputLevel=INFO - ) - hcalBarrelPositionedCellsName = "HCalBarrelPositionedCells" - createHcalBarrelPositionedCells = CreateCaloCellPositionsFCCee( - "CreateHcalBarrelPositionedCells", - OutputLevel=INFO - ) - createHcalBarrelPositionedCells.positionsTool = cellPositionHcalBarrelTool - createHcalBarrelPositionedCells.hits.Path = hcalBarrelCellsName - createHcalBarrelPositionedCells.positionedHits.Path = hcalBarrelPositionedCellsName - - # 3 - compute new cellID of cells based on new readout - removing row information - # We use a RedoSegmentation. Using a RewriteBitField with removeIds=["row"], - # wont work because there are tiles with same layer/theta/phi but different row - # as a consequence there will be multiple cells with same cellID in the output collection - # and this will screw up the SW clustering - hcalBarrelCellsName2 = "HCalBarrelCells2" - - # first we create new hits with the readout without the row information - # and then merge them into new cells - rewriteHCalBarrel = RedoSegmentation("ReSegmentationHcal", - # old bitfield (readout) - 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) - newReadoutName=hcalBarrelReadoutName2, - OutputLevel=INFO, - debugPrint=200, - inhits=hcalBarrelPositionedCellsName, - outhits="HCalBarrelCellsWithoutRow") - - createHcalBarrelCells2 = CreateCaloCells("CreateHCalBarrelCells2", - doCellCalibration=False, - addCellNoise=False, - filterCellNoise=False, - OutputLevel=INFO, - hits=rewriteHCalBarrel.outhits.Path, - cells=hcalBarrelCellsName2) - - # 4 - attach positions to the new cells - from Configurables import CellPositionsHCalBarrelPhiThetaSegTool - hcalBarrelPositionedCellsName2 = "HCalBarrelPositionedCells2" - cellPositionHcalBarrelTool2 = CellPositionsHCalBarrelPhiThetaSegTool( - "CellPositionsHCalBarrel2", - readoutName=hcalBarrelReadoutName2, - OutputLevel=INFO - ) - createHcalBarrelPositionedCells2 = CreateCaloCellPositionsFCCee( - "CreateHCalBarrelPositionedCells2", - OutputLevel=INFO - ) - createHcalBarrelPositionedCells2.positionsTool = cellPositionHcalBarrelTool2 - createHcalBarrelPositionedCells2.hits.Path = hcalBarrelCellsName2 - createHcalBarrelPositionedCells2.positionedHits.Path = hcalBarrelPositionedCellsName2 - - - # createHcalEndcapCells = CreateCaloCells("CreateHcalEndcapCaloCells", - # doCellCalibration=True, - # calibTool=calibHcalEndcap, - # addCellNoise=False, - # filterCellNoise=False, - # OutputLevel=INFO) - # createHcalEndcapCells.hits.Path="HCalEndcapHits" - # createHcalEndcapCells.cells.Path="HCalEndcapCells" - -else: - hcalBarrelCellsName = "emptyCaloCells" - hcalBarrelPositionedCellsName = "emptyCaloCells" - hcalBarrelCellsName2 = "emptyCaloCells" - hcalBarrelPositionedCellsName2 = "emptyCaloCells" - cellPositionHcalBarrelTool = None - cellPositionHcalBarrelTool2 = None - -# Empty cells for parts of calorimeter not implemented yet -createemptycells = CreateEmptyCaloCellsCollection("CreateEmptyCaloCells") -createemptycells.cells.Path = "emptyCaloCells" - -# Produce sliding window clusters (ECAL+HCAL) -towers = CaloTowerToolFCCee("towers", - #deltaThetaTower=4 * 0.009817477 / 4, deltaPhiTower=2 * 2 * pi / 1536., - deltaThetaTower= 0.022180, - deltaPhiTower=2 * pi / 256., - max_layer=25, - ecalBarrelReadoutName=ecalBarrelReadoutName, - ecalEndcapReadoutName=ecalEndcapReadoutName, - ecalFwdReadoutName="", - hcalBarrelReadoutName=hcalBarrelReadoutName2, - hcalExtBarrelReadoutName="", - hcalEndcapReadoutName="", - hcalFwdReadoutName="", - OutputLevel=INFO) -towers.ecalBarrelCells.Path = ecalBarrelPositionedCellsName -towers.ecalEndcapCells.Path = ecalEndcapPositionedCellsName -towers.ecalFwdCells.Path = "emptyCaloCells" - -towers.hcalBarrelCells.Path = hcalBarrelPositionedCellsName2 -towers.hcalExtBarrelCells.Path = "emptyCaloCells" -towers.hcalEndcapCells.Path = "emptyCaloCells" -towers.hcalFwdCells.Path = "emptyCaloCells" - -# Cluster variables (not optimized) -windT = 18 -windP = 34 -posT = 10 -posP = 22 -dupT = 14 -dupP = 26 -finT = 18 -finP = 34 -# Minimal energy to create a cluster in GeV (FCC-ee detectors have to reconstruct low energy particles) -threshold = 0.5 - -createClusters = CreateCaloClustersSlidingWindowFCCee("CreateClusters", - towerTool=towers, - nThetaWindow=windT, nPhiWindow=windP, - nThetaPosition=posT, nPhiPosition=posP, - nThetaDuplicates=dupT, nPhiDuplicates=dupP, - nThetaFinal=finT, nPhiFinal=finP, - energyThreshold=threshold, - energySharingCorrection=False, - attachCells=True, - OutputLevel=INFO - ) -createClusters.clusters.Path = "CaloClusters" -createClusters.clusterCells.Path = "CaloClusterCells" - -correctCaloClusters = CorrectCaloClusters("correctCaloClusters", - inClusters=createClusters.clusters.Path, - outClusters="Corrected" + createClusters.clusters.Path, - systemIDs=[4,8], - numLayers=[12,13], - firstLayerIDs=[0,0], - lastLayerIDs=[11,12], - readoutNames=[ecalBarrelReadoutName,hcalBarrelReadoutName2], - # do not split the following line or it will break scripts that update the values of the corrections - upstreamParameters = [[0.03900891447361534, -4.322941016402328, -139.1811369546787, 0.498342628339746, -3.3545078429754813, -13.99996971344221],[]], - upstreamFormulas=[['[0]+[1]/(x-[2])', '[0]+[1]/(x-[2])'],[]], - # do not split the following line or it will break scripts that update the values of the corrections - downstreamParameters = [[-0.0027661744480442195, 0.006059143775380306, 0.9788596364251927, -1.4951749409378743, -0.08491999337012696, 16.017621428757778],[]], - downstreamFormulas=[['[0]+[1]*x', '[0]+[1]/sqrt(x)', '[0]+[1]/x'],[]], - ## Approximate parameters for the first energy estimate obtained from the benchmark calib for 100 GeV pion - benchmarkParamsApprox = [1.22, 1, 1.04, -0.0019, 25.69, 0.], - ## below parameters and formulas for two functions with ene switch at 9~GeV, the constant term is set to zero - #benchmarkParametrization = [[1.0109, 2.4768, 1., 1.0837, -5.278, -1.9938, 0.0006, -0.2715, 45.29, -5674.31, 126.04, 0.],[2.0432, 0.2902, -0.0709, 0.004, 1., 1.2099, -432.2, 13.73, -0.0023, -0.2572, 1.99, -0.968, 0.259, -0.0152, 0.]], - #benchmarkFormulas = [['[0]+[1]/sqrt(x)', '[0]', '[0]+[1]/(x+[2])', '[0]+[1]/x', '[0]+[1]/(x+[2])', '[0]'],['[0]+[1]*x+[2]*x*x+[3]*x*x*x', '[0]', '[0]+[1]/(x+[2])**2', '[0]+[1]/x', '[0]+[1]*x+[2]*x*x+[3]*x*x*x', '[0]']], - ## below parameters and formulas for one function (redireved in Nov 2023) - benchmarkParametrization = [[2.04, 1.0626, 1., -3.5, 1.0182, -0.26, 0.0004, 45.9, -5906.11, -129.27, 0.],[]], - benchmarkFormulas = [['[0]/sqrt(x)+[1]', '[0]', '[0]/x+[1]', '[0]/x+[1]', '[0]+[1]/(x-[2])', '[0]'],[]], - benchmarkEneSwitch = -1., - upstreamCorr = False, - downstreamCorr = False, - benchmarkCorr = True, - OutputLevel=INFO - ) - -# Output -out = PodioOutput("out", - OutputLevel=INFO) - -if runHCal: - out.outputCommands = ["keep *", "drop emptyCaloCells"] -else: - out.outputCommands = ["keep *", "drop HCal*", "drop emptyCaloCells"] - -if not saveCells: - out.outputCommands.append("drop ECal*Cells*") -if not saveClusterCells: - out.outputCommands.append("drop *ClusterCells*") -if not saveHits: - out.outputCommands.append("drop ECal*Hits*") - out.outputCommands.append("drop HCal*Hits*") - -out.filename = "./output_evts_" + str(Nevts) + "_pdg_" + str(pdgCode) + "_pMin_" + str(momentum) + "_GeV" + "_ThetaMinMax_" + str(thetaMin) + "_" + str( - thetaMax) + "_PhiMinMax_" + str(phiMin) + "_" + str(phiMax) + "_MagneticField_" + str(magneticField) + "_Noise" + str(addNoise) + ".root" - -# CPU information -chra = ChronoAuditor() -audsvc = AuditorSvc() -audsvc.Auditors = [chra] -genAlg.AuditExecute = True -hepmc_converter.AuditExecute = True -geantsim.AuditExecute = True -createEcalBarrelCells.AuditExecute = True -createEcalBarrelPositionedCells.AuditExecute = True -if runHCal: - createHcalBarrelCells.AuditExecute = True -out.AuditExecute = True - -event_counter = EventCounter('event_counter') -event_counter.Frequency = 10 - -ExtSvc = [geoservice, podioevent, geantservice, audsvc] -if dumpGDML: - ExtSvc += [gdmldumpservice] - -TopAlg = [ - event_counter, - genAlg, - hepmc_converter, - geantsim, - createEcalBarrelCells, - createEcalBarrelPositionedCells, - resegmentEcalBarrel, - createEcalBarrelCells2, - createEcalBarrelPositionedCells2, - createEcalEndcapCells, - createEcalEndcapPositionedCells, -] - -if runHCal: - TopAlg += [ - createHcalBarrelCells, - createHcalBarrelPositionedCells, - rewriteHCalBarrel, - createHcalBarrelCells2, - createHcalBarrelPositionedCells2, - # createHcalEndcapCells - ] - -TopAlg += [ - createemptycells, - createClusters, -] - -if applyUpDownstreamBenchmarkCorrections: - TopAlg += [ - correctCaloClusters, - ] - -TopAlg += [ - out -] - -ApplicationMgr( - TopAlg=TopAlg, - EvtSel='NONE', - EvtMax=Nevts, - ExtSvc=ExtSvc, - StopOnSignal=True, -) +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 CreateEmptyCaloCellsCollection +from Configurables import CreateCaloCellPositionsFCCee +from Configurables import CellPositionsECalBarrelModuleThetaSegTool +from Configurables import CellPositionsECalEndcapTurbineSegTool +from Configurables import RedoSegmentation +from Configurables import CreateCaloCells +from Configurables import CalibrateCaloHitsTool +from Configurables import CalibrateInLayersTool +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 +from Configurables import GeoSvc +from Configurables import HepMCToEDMConverter +from Configurables import GenAlg +from Configurables import FCCDataSvc +from Configurables import ReadCaloCrosstalkMap +from Gaudi.Configuration import INFO +# from Gaudi.Configuration import * + +import os + +from GaudiKernel.SystemOfUnits import GeV, tesla, mm +from GaudiKernel.PhysicalConstants import pi, halfpi, twopi +from math import cos, sin, tan + +use_pythia = False +addNoise = False +dumpGDML = False +runHCal = True +# 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 = False +saveCells = False +saveClusterCells = False + +# cluster energy corrections +# simple parametrisations of up/downstream losses +applyUpDownstreamBenchmarkCorrections = True + +# Input for simulations (momentum is expected in GeV!) +# Parameters for the particle gun simulations, dummy if use_pythia is set +# to True +# theta from 80 to 100 degrees corresponds to -0.17 < eta < 0.17 +# reminder: cell granularity in theta = 0.5625 degrees +# (in strips: 0.5625/4=0.14) + +# Nevts = 20000 +Nevts = 2 +# Nevts = 1 +# Nevts=1000 + +# particle momentum and direction +# momentum = 100 # in GeV +momentum = 10. # in GeV +# momentum = 10 # in GeV +thetaMin = 69. # degrees +thetaMax = 69. # degrees +# thetaMin = 89 +# thetaMax = 91 +# thetaMin = 90 # degrees +# thetaMax = 90 # degrees +# phiMin = halfpi +# phiMax = halfpi +phiMin = 0 +phiMax = twopi + +# particle origin +# origR = 1000.0*mm +origR = 0.0 * mm +origTheta = halfpi +origPhi = 0.0 + +# particle type: 11 electron, 13 muon, 22 photon, 111 pi0, 211 pi+ +pdgCode = 211 +# pdgCode = 22 +# pdgCode = 111 +# pdgCode = 211 + +# Set to true if history from Geant4 decays is needed (e.g. to get the +# photons from pi0) +saveG4Hist = False +if (pdgCode == 111): + saveG4Hist = True + +magneticField = False + + +podioevent = FCCDataSvc("EventDataSvc") + +# Particle gun setup + +genAlg = GenAlg() +if use_pythia: + from Configurables import PythiaInterface + pythia8gentool = PythiaInterface() + pythia8gentool.pythiacard = os.path.join( + os.environ.get('PWD', ''), + "MCGeneration/ee_Zgamma_inclusive.cmd" + ) + # pythia8gentool.pythiacard = "MCGeneration/ee_Z_ee.cmd" + pythia8gentool.printPythiaStatistics = False + pythia8gentool.pythiaExtraSettings = [""] + genAlg.SignalProvider = pythia8gentool +else: + from Configurables import MomentumRangeParticleGun + pgun = MomentumRangeParticleGun("ParticleGun") + pgun.PdgCodes = [pdgCode] + pgun.MomentumMin = momentum * GeV + pgun.MomentumMax = momentum * GeV + pgun.PhiMin = phiMin + pgun.PhiMax = phiMax + pgun.ThetaMin = thetaMin * pi / 180. + pgun.ThetaMax = thetaMax * pi / 180. + genAlg.SignalProvider = pgun + +genAlg.hepmc.Path = "hepmc" + +# smear/shift vertex +if origR > 0.0: + origX = origR * cos(origPhi) + origY = origR * sin(origPhi) + origZ = origR / tan(origTheta) + print("Particle gun will be moved to %f %f %f" % (origX, origY, origZ)) + from Configurables import GaussSmearVertex + vertexSmearAndShiftTool = GaussSmearVertex() + vertexSmearAndShiftTool.xVertexSigma = 0. + vertexSmearAndShiftTool.yVertexSigma = 0. + vertexSmearAndShiftTool.tVertexSigma = 0. + vertexSmearAndShiftTool.xVertexMean = origX + vertexSmearAndShiftTool.yVertexMean = origY + vertexSmearAndShiftTool.zVertexMean = origZ + vertexSmearAndShiftTool.tVertexMean = 0. + genAlg.VertexSmearingTool = vertexSmearAndShiftTool + +# hepMC -> EDM converter +hepmc_converter = HepMCToEDMConverter() +hepmc_converter.hepmc.Path = "hepmc" +genParticlesOutputName = "genParticles" +hepmc_converter.GenParticles.Path = genParticlesOutputName +hepmc_converter.hepmcStatusList = [] +hepmc_converter.OutputLevel = INFO + +# Simulation setup +# Detector geometry +geoservice = GeoSvc("GeoSvc") +# 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 = [ + 'FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ALLEGRO_o1_v03.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 +actions = SimG4FullSimActions() + +if saveG4Hist: + actions.enableHistory = True + actions.energyCut = 1.0 * GeV + saveHistTool = SimG4SaveParticleHistory("saveHistory") + +geantservice = SimG4Svc( + "SimG4Svc", + detector='SimG4DD4hepDetector', + physicslist="SimG4FtfpBert", + actions=actions +) + +# from Configurables import GeoToGdmlDumpSvc +if dumpGDML: + from Configurables import GeoToGdmlDumpSvc + gdmldumpservice = GeoToGdmlDumpSvc("GeoToGdmlDumpSvc") + +# Fixed seed to have reproducible results, change it for each job if you +# split one production into several jobs +# Mind that if you leave Gaudi handle random seed and some job start within +# the same second (very likely) you will have duplicates +geantservice.randomNumbersFromGaudi = False +geantservice.seedValue = 4242 + +# Range cut +geantservice.g4PreInitCommands += ["/run/setCut 0.1 mm"] + +# Magnetic field +if magneticField == 1: + field = SimG4ConstantMagneticFieldTool( + "SimG4ConstantMagneticFieldTool", + FieldComponentZ=-2 * tesla, + FieldOn=True, + IntegratorStepper="ClassicalRK4" + ) +else: + field = SimG4ConstantMagneticFieldTool( + "SimG4ConstantMagneticFieldTool", + FieldOn=False + ) + +# Geant4 algorithm +# Translates EDM to G4Event, passes the event to G4, writes out outputs +# via tools and a tool that saves the calorimeter hits + +# Detector readouts +# ECAL +ecalBarrelReadoutName = "ECalBarrelModuleThetaMerged" +ecalBarrelReadoutName2 = "ECalBarrelModuleThetaMerged2" +ecalEndcapReadoutName = "ECalEndcapTurbine" +# HCAL +if runHCal: + hcalBarrelReadoutName = "HCalBarrelReadout" + hcalBarrelReadoutName2 = "BarHCal_Readout_phitheta" + hcalEndcapReadoutName = "HCalEndcapReadout" +else: + hcalBarrelReadoutName = "" + hcalBarrelReadoutName2 = "" + hcalEndcapReadoutName = "" + +# Configure saving of calorimeter hits +ecalBarrelHitsName = "ECalBarrelPositionedHits" +saveECalBarrelTool = SimG4SaveCalHits( + "saveECalBarrelHits", + readoutName=ecalBarrelReadoutName, + OutputLevel=INFO +) +saveECalBarrelTool.CaloHits.Path = ecalBarrelHitsName + +ecalEndcapHitsName = "ECalEndcapPositionedHits" +saveECalEndcapTool = SimG4SaveCalHits( + "saveECalEndcapHits", + readoutName=ecalEndcapReadoutName +) +saveECalEndcapTool.CaloHits.Path = ecalEndcapHitsName + +if runHCal: + hcalBarrelHitsName = "HCalBarrelPositionedHits" + saveHCalTool = SimG4SaveCalHits( + "saveHCalBarrelHits", + readoutName=hcalBarrelReadoutName + ) + saveHCalTool.CaloHits.Path = hcalBarrelHitsName + + # saveHCalEndcapTool = SimG4SaveCalHits( + # "saveHCalEndcapHits", + # readoutName = hcalEndcapReadoutName + # ) + # saveHCalEndcapTool.CaloHits.Path = "HCalEndcapHits" + +# next, create the G4 algorithm, giving the list of names of tools ("XX/YY") +particle_converter = SimG4PrimariesFromEdmTool("EdmConverter") +particle_converter.GenParticles.Path = genParticlesOutputName + +outputTools = [ + saveECalBarrelTool, + saveECalEndcapTool +] +if runHCal: + outputTools += [ + saveHCalTool, + # saveHCalEndcapTool + ] + +if saveG4Hist: + outputTools += [saveHistTool] + +geantsim = SimG4Alg("SimG4Alg", + outputs=outputTools, + eventProvider=particle_converter, + OutputLevel=INFO) + +# Digitization (Merging hits into cells, EM scale calibration) +# EM scale calibration (sampling fraction) +calibEcalBarrel = CalibrateInLayersTool("CalibrateECalBarrel", + samplingFraction=[0.3864252122990472] * 1 + [0.13597644835735828] * 1 + [0.14520427829645913] * 1 + [0.1510076084632846] * 1 + [0.1552347580991012] * 1 + [0.159694330729184] * 1 + [0.1632954482794191] * 1 + [0.16720711037339814] * 1 + [0.17047749048884808] * 1 + [0.17461698117974286] * 1 + [0.1798984163980135] * 1 + [0.17920355117405806] * 1, + readoutName=ecalBarrelReadoutName, + layerFieldName="layer") + +#calibEcalEndcap = CalibrateCaloHitsTool( +# "CalibrateECalEndcap", invSamplingFraction="4.27") +calibEcalEndcap = CalibrateInLayersTool("CalibrateECalEndcap", + samplingFraction = [0.16419] * 1 + [0.192898] * 1 + [0.18783] * 1 + [0.193203] * 1 + [0.193928] * 1 + [0.192286] * 1 + [0.199959] * 1 + [0.200153] * 1 + [0.212635] * 1 + [0.180345] * 1 + [0.18488] * 1 + [0.194762] * 1 + [0.197775] * 1 + [0.200504] * 1 + [0.205555] * 1 + [0.203601] * 1 + [0.210877] * 1 + [0.208376] * 1 + [0.216345] * 1 + [0.201452] * 1 + [0.202134] * 1 + [0.207566] * 1 + [0.208152] * 1 + [0.209889] * 1 + [0.211743] * 1 + [0.213188] * 1 + [0.215864] * 1 + [0.22972] * 1 + [0.192515] * 1 + [0.0103233] * 1, + readoutName=ecalEndcapReadoutName, + layerFieldName="layer") + + +if runHCal: + calibHcells = CalibrateCaloHitsTool( + "CalibrateHCal", invSamplingFraction="30.4") + calibHcalEndcap = CalibrateCaloHitsTool( + "CalibrateHCalEndcap", invSamplingFraction="31.7") + +# Create cells in ECal barrel +# 1. step - merge hits into cells with theta and module segmentation +# (module is a 'physical' cell i.e. lead + LAr + PCB + LAr +lead) +# 2. step - rewrite the cellId using the merged theta-module segmentation +# (merging several modules and severla theta readout cells). +# Add noise at this step if you derived the noise already assuming merged cells + +# read the crosstalk map +readCrosstalkMap = ReadCaloCrosstalkMap("ReadCrosstalkMap", + fileName="https://fccsw.web.cern.ch/fccsw/filesForSimDigiReco/ALLEGRO/ALLEGRO_o1_v03/xtalk_neighbours_map_ecalB_thetamodulemerged.root", + OutputLevel=INFO) + +# Step 1: merge hits into cells according to initial segmentation +ecalBarrelCellsName = "ECalBarrelCells" +createEcalBarrelCells = CreateCaloCells("CreateECalBarrelCells", + doCellCalibration=True, + calibTool=calibEcalBarrel, + crosstalksTool=readCrosstalkMap, + addCrosstalk=False, + addCellNoise=False, + filterCellNoise=False, + addPosition=True, + OutputLevel=INFO, + hits=ecalBarrelHitsName, + cells=ecalBarrelCellsName) + +# Step 2a: compute new cellID of cells based on new readout +# (merged module-theta segmentation with variable merging vs layer) +resegmentEcalBarrel = RedoSegmentation("ReSegmentationEcal", + # old bitfield (readout) + oldReadoutName=ecalBarrelReadoutName, + # specify which fields are going to be altered (deleted/rewritten) + oldSegmentationIds=["module", "theta"], + # new bitfield (readout), with new segmentation (merged modules and theta cells) + newReadoutName=ecalBarrelReadoutName2, + OutputLevel=INFO, + debugPrint=200, + inhits=ecalBarrelCellsName, + outhits="ECalBarrelCellsMerged") + +# Step 2b: merge new cells with same cellID together +# do not apply cell calibration again since cells were already +# calibrated in Step 1 +ecalBarrelCellsName2 = "ECalBarrelCells2" +createEcalBarrelCells2 = CreateCaloCells("CreateECalBarrelCells2", + doCellCalibration=False, + addCellNoise=False, + filterCellNoise=False, + OutputLevel=INFO, + hits="ECalBarrelCellsMerged", + cells=ecalBarrelCellsName2) + +# Add to Ecal barrel cells the position information +# (good for physics, all coordinates set properly) + +cellPositionEcalBarrelTool = CellPositionsECalBarrelModuleThetaSegTool( + "CellPositionsECalBarrel", + readoutName=ecalBarrelReadoutName, + OutputLevel=INFO +) +ecalBarrelPositionedCellsName = "ECalBarrelPositionedCells" +createEcalBarrelPositionedCells = CreateCaloCellPositionsFCCee( + "CreateECalBarrelPositionedCells", + OutputLevel=INFO +) +createEcalBarrelPositionedCells.positionsTool = cellPositionEcalBarrelTool +createEcalBarrelPositionedCells.hits.Path = ecalBarrelCellsName +createEcalBarrelPositionedCells.positionedHits.Path = ecalBarrelPositionedCellsName + +cellPositionEcalBarrelTool2 = CellPositionsECalBarrelModuleThetaSegTool( + "CellPositionsECalBarrel2", + readoutName=ecalBarrelReadoutName2, + OutputLevel=INFO +) +createEcalBarrelPositionedCells2 = CreateCaloCellPositionsFCCee( + "CreateECalBarrelPositionedCells2", + OutputLevel=INFO +) +createEcalBarrelPositionedCells2.positionsTool = cellPositionEcalBarrelTool2 +createEcalBarrelPositionedCells2.hits.Path = ecalBarrelCellsName2 +createEcalBarrelPositionedCells2.positionedHits.Path = "ECalBarrelPositionedCells2" + + +# Create cells in ECal endcap +ecalEndcapCellsName = "ECalEndcapCells" +createEcalEndcapCells = CreateCaloCells("CreateEcalEndcapCaloCells", + doCellCalibration=True, + calibTool=calibEcalEndcap, + addCellNoise=False, + filterCellNoise=False, + OutputLevel=INFO, + hits=ecalEndcapHitsName, + cells=ecalEndcapCellsName) + +# Add to Ecal endcap cells the position information +# (good for physics, all coordinates set properly) +# not yet merged!! +cellPositionEcalEndcapTool = CellPositionsECalEndcapTurbineSegTool( + "CellPositionsECalEndcap", + readoutName=ecalEndcapReadoutName, + OutputLevel=INFO +) +ecalEndcapPositionedCellsName = "ECalEndcapPositionedCells" +createEcalEndcapPositionedCells = CreateCaloCellPositionsFCCee( + "CreateECalEndcapPositionedCells", + OutputLevel=INFO +) +createEcalEndcapPositionedCells.positionsTool = cellPositionEcalEndcapTool +createEcalEndcapPositionedCells.hits.Path = ecalEndcapCellsName +createEcalEndcapPositionedCells.positionedHits.Path = ecalEndcapPositionedCellsName + +if runHCal: + # Create cells in HCal + # 1 - merge hits into cells with the default readout + hcalBarrelCellsName = "HCalBarrelCells" + createHcalBarrelCells = CreateCaloCells("CreateHCalBarrelCells", + doCellCalibration=True, + calibTool=calibHcells, + addCellNoise=False, + filterCellNoise=False, + addPosition=True, + hits=hcalBarrelHitsName, + cells=hcalBarrelCellsName, + OutputLevel=INFO) + + # 2 - attach positions to the cells (cell positions needed for RedoSegmentation!) + from Configurables import CellPositionsHCalBarrelPhiThetaSegTool + cellPositionHcalBarrelTool = CellPositionsHCalBarrelPhiThetaSegTool( + "CellPositionsHCalBarrel", + readoutName=hcalBarrelReadoutName, + OutputLevel=INFO + ) + hcalBarrelPositionedCellsName = "HCalBarrelPositionedCells" + createHcalBarrelPositionedCells = CreateCaloCellPositionsFCCee( + "CreateHcalBarrelPositionedCells", + OutputLevel=INFO + ) + createHcalBarrelPositionedCells.positionsTool = cellPositionHcalBarrelTool + createHcalBarrelPositionedCells.hits.Path = hcalBarrelCellsName + createHcalBarrelPositionedCells.positionedHits.Path = hcalBarrelPositionedCellsName + + # 3 - compute new cellID of cells based on new readout - removing row information + # We use a RedoSegmentation. Using a RewriteBitField with removeIds=["row"], + # wont work because there are tiles with same layer/theta/phi but different row + # as a consequence there will be multiple cells with same cellID in the output collection + # and this will screw up the SW clustering + hcalBarrelCellsName2 = "HCalBarrelCells2" + + # first we create new hits with the readout without the row information + # and then merge them into new cells + rewriteHCalBarrel = RedoSegmentation("ReSegmentationHcal", + # old bitfield (readout) + 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) + newReadoutName=hcalBarrelReadoutName2, + OutputLevel=INFO, + debugPrint=200, + inhits=hcalBarrelPositionedCellsName, + outhits="HCalBarrelCellsWithoutRow") + + createHcalBarrelCells2 = CreateCaloCells("CreateHCalBarrelCells2", + doCellCalibration=False, + addCellNoise=False, + filterCellNoise=False, + OutputLevel=INFO, + hits=rewriteHCalBarrel.outhits.Path, + cells=hcalBarrelCellsName2) + + # 4 - attach positions to the new cells + from Configurables import CellPositionsHCalBarrelPhiThetaSegTool + hcalBarrelPositionedCellsName2 = "HCalBarrelPositionedCells2" + cellPositionHcalBarrelTool2 = CellPositionsHCalBarrelPhiThetaSegTool( + "CellPositionsHCalBarrel2", + readoutName=hcalBarrelReadoutName2, + OutputLevel=INFO + ) + createHcalBarrelPositionedCells2 = CreateCaloCellPositionsFCCee( + "CreateHCalBarrelPositionedCells2", + OutputLevel=INFO + ) + createHcalBarrelPositionedCells2.positionsTool = cellPositionHcalBarrelTool2 + createHcalBarrelPositionedCells2.hits.Path = hcalBarrelCellsName2 + createHcalBarrelPositionedCells2.positionedHits.Path = hcalBarrelPositionedCellsName2 + + + # createHcalEndcapCells = CreateCaloCells("CreateHcalEndcapCaloCells", + # doCellCalibration=True, + # calibTool=calibHcalEndcap, + # addCellNoise=False, + # filterCellNoise=False, + # OutputLevel=INFO) + # createHcalEndcapCells.hits.Path="HCalEndcapHits" + # createHcalEndcapCells.cells.Path="HCalEndcapCells" + +else: + hcalBarrelCellsName = "emptyCaloCells" + hcalBarrelPositionedCellsName = "emptyCaloCells" + hcalBarrelCellsName2 = "emptyCaloCells" + hcalBarrelPositionedCellsName2 = "emptyCaloCells" + cellPositionHcalBarrelTool = None + cellPositionHcalBarrelTool2 = None + +# Empty cells for parts of calorimeter not implemented yet +createemptycells = CreateEmptyCaloCellsCollection("CreateEmptyCaloCells") +createemptycells.cells.Path = "emptyCaloCells" + +# Produce sliding window clusters (ECAL+HCAL) +towers = CaloTowerToolFCCee("towers", + #deltaThetaTower=4 * 0.009817477 / 4, deltaPhiTower=2 * 2 * pi / 1536., + deltaThetaTower= 0.022180, + deltaPhiTower=2 * pi / 256., + max_layer=25, + ecalBarrelReadoutName=ecalBarrelReadoutName, + ecalEndcapReadoutName=ecalEndcapReadoutName, + ecalFwdReadoutName="", + hcalBarrelReadoutName=hcalBarrelReadoutName2, + hcalExtBarrelReadoutName="", + hcalEndcapReadoutName="", + hcalFwdReadoutName="", + OutputLevel=INFO) +towers.ecalBarrelCells.Path = ecalBarrelPositionedCellsName +towers.ecalEndcapCells.Path = ecalEndcapPositionedCellsName +towers.ecalFwdCells.Path = "emptyCaloCells" + +towers.hcalBarrelCells.Path = hcalBarrelPositionedCellsName2 +towers.hcalExtBarrelCells.Path = "emptyCaloCells" +towers.hcalEndcapCells.Path = "emptyCaloCells" +towers.hcalFwdCells.Path = "emptyCaloCells" + +# Cluster variables (not optimized) +windT = 18 +windP = 34 +posT = 10 +posP = 22 +dupT = 14 +dupP = 26 +finT = 18 +finP = 34 +# Minimal energy to create a cluster in GeV (FCC-ee detectors have to reconstruct low energy particles) +threshold = 0.5 + +createClusters = CreateCaloClustersSlidingWindowFCCee("CreateClusters", + towerTool=towers, + nThetaWindow=windT, nPhiWindow=windP, + nThetaPosition=posT, nPhiPosition=posP, + nThetaDuplicates=dupT, nPhiDuplicates=dupP, + nThetaFinal=finT, nPhiFinal=finP, + energyThreshold=threshold, + energySharingCorrection=False, + attachCells=True, + OutputLevel=INFO + ) +createClusters.clusters.Path = "CaloClusters" +createClusters.clusterCells.Path = "CaloClusterCells" + +correctCaloClusters = CorrectCaloClusters("correctCaloClusters", + inClusters=createClusters.clusters.Path, + outClusters="Corrected" + createClusters.clusters.Path, + systemIDs=[4,8], + numLayers=[12,13], + firstLayerIDs=[0,0], + lastLayerIDs=[11,12], + readoutNames=[ecalBarrelReadoutName,hcalBarrelReadoutName2], + # do not split the following line or it will break scripts that update the values of the corrections + upstreamParameters = [[0.03900891447361534, -4.322941016402328, -139.1811369546787, 0.498342628339746, -3.3545078429754813, -13.99996971344221],[]], + upstreamFormulas=[['[0]+[1]/(x-[2])', '[0]+[1]/(x-[2])'],[]], + # do not split the following line or it will break scripts that update the values of the corrections + downstreamParameters = [[-0.0027661744480442195, 0.006059143775380306, 0.9788596364251927, -1.4951749409378743, -0.08491999337012696, 16.017621428757778],[]], + downstreamFormulas=[['[0]+[1]*x', '[0]+[1]/sqrt(x)', '[0]+[1]/x'],[]], + ## Approximate parameters for the first energy estimate obtained from the benchmark calib for 100 GeV pion + benchmarkParamsApprox = [1.22, 1, 1.04, -0.0019, 25.69, 0.], + ## below parameters and formulas for two functions with ene switch at 9~GeV, the constant term is set to zero + #benchmarkParametrization = [[1.0109, 2.4768, 1., 1.0837, -5.278, -1.9938, 0.0006, -0.2715, 45.29, -5674.31, 126.04, 0.],[2.0432, 0.2902, -0.0709, 0.004, 1., 1.2099, -432.2, 13.73, -0.0023, -0.2572, 1.99, -0.968, 0.259, -0.0152, 0.]], + #benchmarkFormulas = [['[0]+[1]/sqrt(x)', '[0]', '[0]+[1]/(x+[2])', '[0]+[1]/x', '[0]+[1]/(x+[2])', '[0]'],['[0]+[1]*x+[2]*x*x+[3]*x*x*x', '[0]', '[0]+[1]/(x+[2])**2', '[0]+[1]/x', '[0]+[1]*x+[2]*x*x+[3]*x*x*x', '[0]']], + ## below parameters and formulas for one function (redireved in Nov 2023) + benchmarkParametrization = [[2.04, 1.0626, 1., -3.5, 1.0182, -0.26, 0.0004, 45.9, -5906.11, -129.27, 0.],[]], + benchmarkFormulas = [['[0]/sqrt(x)+[1]', '[0]', '[0]/x+[1]', '[0]/x+[1]', '[0]+[1]/(x-[2])', '[0]'],[]], + benchmarkEneSwitch = -1., + upstreamCorr = False, + downstreamCorr = False, + benchmarkCorr = True, + OutputLevel=INFO + ) + +# Output +out = PodioOutput("out", + OutputLevel=INFO) + +if runHCal: + out.outputCommands = ["keep *", "drop emptyCaloCells"] +else: + out.outputCommands = ["keep *", "drop HCal*", "drop emptyCaloCells"] + +if not saveCells: + out.outputCommands.append("drop ECal*Cells*") +if not saveClusterCells: + out.outputCommands.append("drop *ClusterCells*") +if not saveHits: + out.outputCommands.append("drop ECal*Hits*") + out.outputCommands.append("drop HCal*Hits*") + +out.filename = "./output_evts_" + str(Nevts) + "_pdg_" + str(pdgCode) + "_pMin_" + str(momentum) + "_GeV" + "_ThetaMinMax_" + str(thetaMin) + "_" + str( + thetaMax) + "_PhiMinMax_" + str(phiMin) + "_" + str(phiMax) + "_MagneticField_" + str(magneticField) + "_Noise" + str(addNoise) + ".root" + +# CPU information +chra = ChronoAuditor() +audsvc = AuditorSvc() +audsvc.Auditors = [chra] +genAlg.AuditExecute = True +hepmc_converter.AuditExecute = True +geantsim.AuditExecute = True +createEcalBarrelCells.AuditExecute = True +createEcalBarrelPositionedCells.AuditExecute = True +if runHCal: + createHcalBarrelCells.AuditExecute = True +out.AuditExecute = True + +#event_counter = EventCounter('event_counter') +#event_counter.Frequency = 10 + +ExtSvc = [geoservice, podioevent, geantservice, audsvc] +if dumpGDML: + ExtSvc += [gdmldumpservice] + +TopAlg = [ +# event_counter, + genAlg, + hepmc_converter, + geantsim, + createEcalBarrelCells, + createEcalBarrelPositionedCells, + resegmentEcalBarrel, + createEcalBarrelCells2, + createEcalBarrelPositionedCells2, + createEcalEndcapCells, + createEcalEndcapPositionedCells, +] + +if runHCal: + TopAlg += [ + createHcalBarrelCells, + createHcalBarrelPositionedCells, + rewriteHCalBarrel, + createHcalBarrelCells2, + createHcalBarrelPositionedCells2, + # createHcalEndcapCells + ] + +TopAlg += [ + createemptycells, + createClusters, +] + +if applyUpDownstreamBenchmarkCorrections: + TopAlg += [ + correctCaloClusters, + ] + +TopAlg += [ + out +] + +ApplicationMgr( + TopAlg=TopAlg, + EvtSel='NONE', + EvtMax=Nevts, + ExtSvc=ExtSvc, + StopOnSignal=True, +)