diff --git a/FCCSW_ecal/neighbours_theta.py b/FCCSW_ecal/neighbours_theta.py new file mode 100644 index 0000000..c58ce82 --- /dev/null +++ b/FCCSW_ecal/neighbours_theta.py @@ -0,0 +1,44 @@ +import os +from Gaudi.Configuration import * + +# Detector geometry +from Configurables import GeoSvc +geoservice = GeoSvc("GeoSvc") +# if FCC_DETECTORS is empty, this should use relative path to working directory +path_to_detector = os.environ.get("FCCDETECTORS", "") +print(path_to_detector) +detectors_to_use=[ + 'Detector/DetFCCeeIDEA-LAr/compact/FCCee_DectMaster_thetamodulemerged.xml', + #'Detector/DetFCCeeIDEA-LAr/compact/FCCee_DectMaster.xml', +] + +# prefix all xmls with path_to_detector +geoservice.detectors = [os.path.join(path_to_detector, _det) for _det in detectors_to_use] +geoservice.OutputLevel = INFO + +# Geant4 service +# Configures the Geant simulation: geometry, physics list and user actions +from Configurables import CreateFCCeeCaloNeighbours +neighbours = CreateFCCeeCaloNeighbours("neighbours", + outputFileName = "neighbours_map_barrel_thetamodulemerged.root", + readoutNamesModuleTheta = ["ECalBarrelModuleThetaMerged"], +# readoutNamesModuleTheta = ["ECalBarrelModuleThetaMerged2"], + systemNamesModuleTheta = ["system"], + systemValuesModuleTheta = [4], + activeFieldNamesModuleTheta = ["layer"], + activeVolumesNumbers = [12], + #activeVolumesTheta = [1.2524, 1.2234, 1.1956, 1.1561, 1.1189, 1.0839, 1.0509, 0.9999, 0.9534, 0.91072], + includeDiagonalCells = False, + readoutNamesVolumes = [], + connectBarrels = False, + OutputLevel = DEBUG) + +# ApplicationMgr +from Configurables import ApplicationMgr +ApplicationMgr( TopAlg = [], + EvtSel = 'NONE', + EvtMax = 1, + # order is important, as GeoSvc is needed by G4SimSvc + ExtSvc = [geoservice, neighbours], + OutputLevel=INFO +) diff --git a/FCCSW_ecal/noise_map_theta.py b/FCCSW_ecal/noise_map_theta.py new file mode 100644 index 0000000..8f7638d --- /dev/null +++ b/FCCSW_ecal/noise_map_theta.py @@ -0,0 +1,75 @@ +from Gaudi.Configuration import * +# Detector geometry +from Configurables import GeoSvc +geoservice = GeoSvc("GeoSvc") +# if FCC_DETECTORS is empty, this should use relative path to working directory +import os +path_to_detector = os.environ.get("FCCDETECTORS", "") +print(path_to_detector) +detectors_to_use=[ + 'Detector/DetFCCeeIDEA-LAr/compact/FCCee_DectMaster_thetamodulemerged.xml' + ] +# prefix all xmls with path_to_detector +geoservice.detectors = [os.path.join(path_to_detector, _det) for _det in detectors_to_use] +geoservice.OutputLevel = INFO + +ecalBarrelReadoutName = "ECalBarrelModuleThetaMerged" +#hcalBarrelReadoutName = "ECalBarrelPhiEta" +hcalBarrelReadoutName = "HCalBarrelReadout" +BarrelNoisePath = os.environ['FCCBASEDIR']+"/LAr_scripts/data/elecNoise_ecalBarrelFCCee_theta.root" +ecalBarrelNoiseHistName = "h_elecNoise_fcc_" + +from Configurables import CellPositionsECalBarrelModuleThetaSegTool +ECalBcells = CellPositionsECalBarrelModuleThetaSegTool("CellPositionsECalBarrel", + readoutName = ecalBarrelReadoutName) +# OutputLevel = DEBUG) +#print(ECalBcells) + +from Configurables import CreateFCCeeCaloNoiseLevelMap, ReadNoiseFromFileTool, ReadNoiseVsThetaFromFileTool +ECalNoiseTool = ReadNoiseVsThetaFromFileTool("ReadNoiseFromFileToolECal", + useSegmentation = False, + cellPositionsTool = ECalBcells, + readoutName = ecalBarrelReadoutName, + noiseFileName = BarrelNoisePath, + elecNoiseHistoName = ecalBarrelNoiseHistName, + setNoiseOffset = False, + activeFieldName = "layer", + addPileup = False, + numRadialLayers = 12, + scaleFactor = 1/1000., #MeV to GeV + OutputLevel = INFO) + +HCalNoiseTool = ReadNoiseFromFileTool("ReadNoiseFromFileToolHCal", + readoutName = hcalBarrelReadoutName, + noiseFileName = BarrelNoisePath, + elecNoiseHistoName = ecalBarrelNoiseHistName, + setNoiseOffset = False, + activeFieldName = "layer", + addPileup = False, + numRadialLayers = 12, + scaleFactor = 1/1000., #MeV to GeV + OutputLevel = INFO) + +noisePerCell = CreateFCCeeCaloNoiseLevelMap("noisePerCell", + ECalBarrelNoiseTool = ECalNoiseTool, + ecalBarrelSysId = 4, + HCalBarrelNoiseTool = HCalNoiseTool, + readoutNamesModuleTheta=[ecalBarrelReadoutName], + systemNamesModuleTheta=["system"], + systemValuesModuleTheta=[4], + activeFieldNamesModuleTheta=["layer"], + activeVolumesNumbers = [12], + #activeVolumesEta = [1.2524, 1.2234, 1.1956, 1.1561, 1.1189, 1.0839, 1.0509, 0.9999, 0.9534, 0.91072], + readoutNamesVolumes = [], + outputFileName = "cellNoise_map_electronicsNoiseLevel_thetamodulemerged.root", + OutputLevel = DEBUG) + +# ApplicationMgr +from Configurables import ApplicationMgr +ApplicationMgr( TopAlg = [], + EvtSel = 'NONE', + EvtMax = 1, + # order is important, as GeoSvc is needed by G4SimSvc + ExtSvc = [geoservice, noisePerCell], + OutputLevel=INFO + ) diff --git a/FCCSW_ecal/printNeighbours.C b/FCCSW_ecal/printNeighbours.C new file mode 100644 index 0000000..0023e22 --- /dev/null +++ b/FCCSW_ecal/printNeighbours.C @@ -0,0 +1,75 @@ +TTree* T = nullptr; +const std::string filename = "neighbours_map_barrel_thetamodulemerged.root"; +const std::string treename = "neighbours"; +ULong64_t cID; +std::vector *neighbours=0; + +// HELPER FUNCTIONS + +// extract layer number from cellID +ULong_t Layer(ULong_t cellID) { + const ULong_t mask = (1<<8) -1; + return (cellID >> 11) & mask; +} + +// extract module number from cellID +ULong_t Module(ULong_t cellID) { + const ULong_t mask = (1<<11) -1; + return (cellID >> 19) & mask; +} + +// extract theta bin from cellID +ULong_t ThetaBin(ULong_t cellID) { + const ULong_t mask = (1<<10) -1; + return (cellID >> 30) & mask; +} + +void printCell(ULong_t cellID) { + cout << "cellID: " << cellID << endl; + cout << "Layer: " << Layer(cellID) << endl; + cout << "Theta bin: " << ThetaBin(cellID) << endl; + cout << "Module: " << Module(cellID) << endl; + cout << endl; +} + +void LoadNeighboursMap() { + if (T==nullptr) { + TFile* f = TFile::Open(filename.c_str(),"READ"); + T = (TTree*) f->Get(treename.c_str()); + T->SetBranchAddress("cellId", &cID); + T->SetBranchAddress("neighbours", &neighbours); + } +} + +void printCellAndNeighbours(ULong64_t iEntry) { + T->GetEntry(iEntry); + cout << "=================================================" << endl; + cout << endl; + printCell(cID); + cout << "Neighbours: " << endl << endl; + for (unsigned int i=0; isize(); i++) { + printCell(neighbours->at(i)); + } + cout << "=================================================" << endl; +} + +void printNeighboursOfCell(ULong_t cellID) { + LoadNeighboursMap(); + for (ULong64_t iEntry=0; iEntryGetEntries(); iEntry++) { + T->GetEntry(iEntry); + if (cID == cellID) { + printCellAndNeighbours(iEntry); + return; + } + } + cout << "CellID not found" << endl; +} + + +void printNeighbours(int n=10) { + LoadNeighboursMap(); + for (int i=0; iUniform(T->GetEntries()); + printCellAndNeighbours(entry); + } +} diff --git a/FCCSW_ecal/runClueAndTopoAndSlidingWindowAndCaloSim.py b/FCCSW_ecal/runClueAndTopoAndSlidingWindowAndCaloSim.py index 31dd4cc..5af9853 100644 --- a/FCCSW_ecal/runClueAndTopoAndSlidingWindowAndCaloSim.py +++ b/FCCSW_ecal/runClueAndTopoAndSlidingWindowAndCaloSim.py @@ -373,6 +373,7 @@ #positionsHECTool = HECcells, #positionsEMFwdTool = ECalFwdcells, #positionsHFwdTool = HCalFwdcells, + readoutName = ecalBarrelReadoutNamePhiEta, seedSigma = 4, neighbourSigma = 2, lastNeighbourSigma = 0, diff --git a/FCCSW_ecal/runTopoAndSlidingWindowAndCaloSim.py b/FCCSW_ecal/runTopoAndSlidingWindowAndCaloSim.py index 043d66e..a4f3ba2 100644 --- a/FCCSW_ecal/runTopoAndSlidingWindowAndCaloSim.py +++ b/FCCSW_ecal/runTopoAndSlidingWindowAndCaloSim.py @@ -399,6 +399,7 @@ #positionsHECTool = HECcells, #positionsEMFwdTool = ECalFwdcells, #positionsHFwdTool = HCalFwdcells, + readoutName = ecalBarrelReadoutNamePhiEta, seedSigma = 4, neighbourSigma = 2, lastNeighbourSigma = 0, diff --git a/FCCSW_ecal/run_thetamodulemerged.py b/FCCSW_ecal/run_thetamodulemerged.py new file mode 100644 index 0000000..914cb9d --- /dev/null +++ b/FCCSW_ecal/run_thetamodulemerged.py @@ -0,0 +1,576 @@ +import os +import copy + +from GaudiKernel.SystemOfUnits import MeV, GeV, tesla + +use_pythia = False +addNoise = False +dumpGDML = False + +# 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 = 50 +#Nevts = 1 +Nevts=1000 +#momentum = 100 # in GeV +momentum = 50 # in GeV +#momentum = 10 # in GeV +_pi = 3.14159 +#thetaMin = 40 # degrees +#thetaMax = 140 # degrees +thetaMin = 90 # degrees +thetaMax = 90 # degrees +phiMin = _pi/2. +phiMax = _pi/2. +#phiMin = 0 +#phiMax = 2 * _pi + +# particle type: 11 electron, 13 muon, 22 photon, 111 pi0, 211 pi+ +pdgCode = 11 +#pdgCode = 22 +#pdgCode = 111 + +# 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 + + +from Gaudi.Configuration import * + +from Configurables import FCCDataSvc +podioevent = FCCDataSvc("EventDataSvc") + +################## Particle gun setup + +from Configurables import GenAlg +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 + # 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") + 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" + +from Configurables import HepMCToEDMConverter +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 +from Configurables import GeoSvc +geoservice = GeoSvc("GeoSvc") +# if FCC_DETECTORS is empty, this should use relative path to working directory +path_to_detector = os.environ.get("FCCDETECTORS", "") +print(path_to_detector) +detectors_to_use=[ + 'Detector/DetFCCeeIDEA-LAr/compact/FCCee_DectMaster_thetamodulemerged.xml', +] +# prefix all xmls with path_to_detector +geoservice.detectors = [ + os.path.join(path_to_detector, _det) for _det in detectors_to_use +] +geoservice.OutputLevel = INFO + +# Geant4 service +# Configures the Geant simulation: geometry, physics list and user actions +from Configurables import ( + SimG4FullSimActions, SimG4Alg, + SimG4PrimariesFromEdmTool, SimG4SaveParticleHistory +) +actions = SimG4FullSimActions() + +if saveG4Hist: + actions.enableHistory=True + actions.energyCut = 1.0 * GeV + saveHistTool = SimG4SaveParticleHistory("saveHistory") + +from Configurables import SimG4Svc +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 +from Configurables import SimG4ConstantMagneticFieldTool +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 = "ECalEndcapPhiEta" +# HCAL +hcalBarrelReadoutName = "HCalBarrelReadout" +hcalEndcapReadoutName = "HCalEndcapReadout" + +# Configure saving of calorimeter hits +ecalBarrelHitsName = "ECalBarrelPositionedHits" +from Configurables import SimG4SaveCalHits +saveECalBarrelTool = SimG4SaveCalHits( + "saveECalBarrelHits", + readoutName = ecalBarrelReadoutName, + OutputLevel = INFO +) +saveECalBarrelTool.CaloHits.Path = ecalBarrelHitsName + +saveECalEndcapTool = SimG4SaveCalHits( + "saveECalEndcapHits", + readoutName = ecalEndcapReadoutName +) +saveECalEndcapTool.CaloHits.Path = "ECalEndcapHits" + +saveHCalTool = SimG4SaveCalHits( + "saveHCalBarrelHits", + readoutName = hcalBarrelReadoutName +) +saveHCalTool.CaloHits.Path = "HCalBarrelPositionedHits" + +#saveHCalEndcapTool = SimG4SaveCalHits( +# "saveHCalEndcapHits", +# readoutName = hcalEndcapReadoutName +#) +#saveHCalEndcapTool.CaloHits.Path = "HCalEndcapHits" + +# next, create the G4 algorithm, giving the list of names of tools ("XX/YY") +from Configurables import SimG4PrimariesFromEdmTool +particle_converter = SimG4PrimariesFromEdmTool("EdmConverter") +particle_converter.GenParticles.Path = genParticlesOutputName + +from Configurables import SimG4Alg +outputTools = [ + saveECalBarrelTool, + saveHCalTool, + saveECalEndcapTool, + #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) +from Configurables import CalibrateInLayersTool +calibEcalBarrel = CalibrateInLayersTool("CalibrateECalBarrel", + 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") + +from Configurables import CalibrateCaloHitsTool +calibHcells = CalibrateCaloHitsTool("CalibrateHCal", invSamplingFraction="31.4") +calibEcalEndcap = CalibrateCaloHitsTool("CalibrateECalEndcap", invSamplingFraction="4.27") +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 + +# Step 1: merge hits into cells according to initial segmentation +from Configurables import CreateCaloCells +EcalBarrelCellsName = "ECalBarrelCells" +createEcalBarrelCells = CreateCaloCells("CreateECalBarrelCells", + doCellCalibration=True, + calibTool = calibEcalBarrel, + 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) +from Configurables import RedoSegmentation +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) +from Configurables import CellPositionsECalBarrelModuleThetaSegTool +from Configurables import CreateCaloCellPositionsFCCee + +cellPositionEcalBarrelTool = CellPositionsECalBarrelModuleThetaSegTool( + "CellPositionsECalBarrel", + readoutName = ecalBarrelReadoutName, + OutputLevel = INFO +) +createEcalBarrelPositionedCells = CreateCaloCellPositionsFCCee( + "ECalBarrelPositionedCells", + OutputLevel = INFO +) +createEcalBarrelPositionedCells.positionsECalBarrelTool = cellPositionEcalBarrelTool +createEcalBarrelPositionedCells.hits.Path = EcalBarrelCellsName +createEcalBarrelPositionedCells.positionedHits.Path = "ECalBarrelPositionedCells" + +cellPositionEcalBarrelTool2= CellPositionsECalBarrelModuleThetaSegTool( + "CellPositionsECalBarrel2", + readoutName = ecalBarrelReadoutName2, + OutputLevel = INFO +) +createEcalBarrelPositionedCells2 = CreateCaloCellPositionsFCCee( + "ECalBarrelPositionedCells2", + OutputLevel = INFO +) +createEcalBarrelPositionedCells2.positionsECalBarrelTool = cellPositionEcalBarrelTool2 +createEcalBarrelPositionedCells2.hits.Path = EcalBarrelCellsName2 +createEcalBarrelPositionedCells2.positionedHits.Path = "ECalBarrelPositionedCells2" + +# Create cells in HCal +# 1. step - merge hits into cells with the default readout +createHcalBarrelCells = CreateCaloCells("CreateHCaloCells", + doCellCalibration=True, + calibTool=calibHcells, + addCellNoise = False, + filterCellNoise = False, + OutputLevel = INFO, + hits="HCalBarrelHits", + cells="HCalBarrelCells") + +# 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" + +#createHcalEndcapCells = CreateCaloCells("CreateHcalEndcapCaloCells", +# doCellCalibration=True, +# calibTool=calibHcalEndcap, +# addCellNoise=False, +# filterCellNoise=False, +# OutputLevel=INFO) +#createHcalEndcapCells.hits.Path="HCalEndcapHits" +#createHcalEndcapCells.cells.Path="HCalEndcapCells" + + +#Empty cells for parts of calorimeter not implemented yet +from Configurables import CreateEmptyCaloCellsCollection +createemptycells = CreateEmptyCaloCellsCollection("CreateEmptyCaloCells") +createemptycells.cells.Path = "emptyCaloCells" + +# Produce sliding window clusters +from Configurables import CaloTowerTool +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.hcalExtBarrelCells.Path = "emptyCaloCells" +towers.hcalEndcapCells.Path = "emptyCaloCells" +towers.hcalFwdCells.Path = "emptyCaloCells" + +# Cluster variables +windE = 9 +windP = 17 +posE = 5 +posP = 11 +dupE = 7 +dupP = 13 +finE = 9 +finP = 17 +# Minimal energy to create a cluster in GeV (FCC-ee detectors have to reconstruct low energy particles) +threshold = 0.040 + +from Configurables import CreateCaloClustersSlidingWindow +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.positionsECalBarrelTool = cellPositionEcalBarrelTool +createEcalBarrelPositionedCaloClusterCells.hits.Path = "CaloClusterCells" +createEcalBarrelPositionedCaloClusterCells.positionedHits.Path = "PositionedCaloClusterCells" + +from Configurables import CorrectCaloClusters +correctCaloClusters = CorrectCaloClusters("correctCaloClusters", + inClusters = createClusters.clusters.Path, + outClusters = "Corrected"+createClusters.clusters.Path, + numLayers = [12], + firstLayerIDs = [0], + lastLayerIDs = [11], + ## 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 + ) + +# TOPO CLUSTERS PRODUCTION +from Configurables import ( + CaloTopoClusterInputTool, CaloTopoClusterFCCee, + TopoCaloNeighbours, TopoCaloNoisyCells +) +createTopoInput = CaloTopoClusterInputTool("CreateTopoInput", + ecalBarrelReadoutName = ecalBarrelReadoutName, + #ecalBarrelReadoutName = ecalBarrelReadoutName2, + ecalEndcapReadoutName = "", + ecalFwdReadoutName = "", + hcalBarrelReadoutName = "", + hcalExtBarrelReadoutName = "", + hcalEndcapReadoutName = "", + hcalFwdReadoutName = "", + OutputLevel = INFO) + +createTopoInput.ecalBarrelCells.Path = "ECalBarrelPositionedCells" +#createTopoInput.ecalBarrelCells.Path = "ECalBarrelPositionedCells2" +createTopoInput.ecalEndcapCells.Path = "emptyCaloCells" +createTopoInput.ecalFwdCells.Path = "emptyCaloCells" +createTopoInput.hcalBarrelCells.Path = "emptyCaloCells" +createTopoInput.hcalExtBarrelCells.Path = "emptyCaloCells" +createTopoInput.hcalEndcapCells.Path = "emptyCaloCells" +createTopoInput.hcalFwdCells.Path = "emptyCaloCells" + +readNeighboursMap = TopoCaloNeighbours("ReadNeighboursMap", + fileName = os.environ['FCCBASEDIR']+"/LAr_scripts/data/neighbours_map_barrel_thetamodulemerged.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, + #positionsECalBarrelTool = cellPositionEcalBarrelTool2, + #positionsHCalBarrelTool = HCalBcells, + #positionsHCalExtBarrelTool = HCalExtBcells, + #positionsEMECTool = EMECcells, + #positionsHECTool = HECcells, + #positionsEMFwdTool = ECalFwdcells, + #positionsHFwdTool = HCalFwdcells, + seedSigma = 4, + neighbourSigma = 2, + lastNeighbourSigma = 0, + OutputLevel = INFO) +createTopoClusters.clusters.Path ="CaloTopoClusters" +createTopoClusters.clusterCells.Path = "CaloTopoClusterCells" + +createEcalBarrelPositionedCaloTopoClusterCells = CreateCaloCellPositionsFCCee( + "ECalBarrelPositionedCaloTopoClusterCells", + OutputLevel = INFO +) +#createEcalBarrelPositionedCaloTopoClusterCells.positionsECalBarrelTool = cellPositionEcalBarrelTool2 +createEcalBarrelPositionedCaloTopoClusterCells.positionsECalBarrelTool = cellPositionEcalBarrelTool +createEcalBarrelPositionedCaloTopoClusterCells.hits.Path = "CaloTopoClusterCells" +createEcalBarrelPositionedCaloTopoClusterCells.positionedHits.Path = "PositionedCaloTopoClusterCells" + +correctCaloTopoClusters = CorrectCaloClusters( + "correctCaloTopoClusters", + inClusters = createTopoClusters.clusters.Path, + outClusters = "Corrected"+createTopoClusters.clusters.Path, + numLayers = [12], + firstLayerIDs = [0], + lastLayerIDs = [11], + ## 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 +from Configurables import PodioOutput +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"] +out.outputCommands = ["keep *", "drop HCal*", "drop emptyCaloCells"] + +import uuid +# 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 +from Configurables import AuditorSvc, ChronoAuditor +chra = ChronoAuditor() +audsvc = AuditorSvc() +audsvc.Auditors = [chra] +genAlg.AuditExecute = True +hepmc_converter.AuditExecute = True +geantsim.AuditExecute = True +createEcalBarrelCells.AuditExecute = True +#createHcalBarrelCells.AuditExecute = True +createEcalBarrelPositionedCells.AuditExecute = True +createTopoClusters.AuditExecute = True +out.AuditExecute = True + +from Configurables import EventCounter +event_counter = EventCounter('event_counter') +event_counter.Frequency = 10 + +ExtSvc = [geoservice, podioevent, geantservice, audsvc] +if dumpGDML: + ExtSvc += [gdmldumpservice] + +from Configurables import ApplicationMgr +ApplicationMgr( + TopAlg = [ + event_counter, + genAlg, + hepmc_converter, + geantsim, + createEcalBarrelCells, + createEcalBarrelPositionedCells, + resegmentEcalBarrel, + createEcalBarrelCells2, + createEcalBarrelPositionedCells2, + createEcalEndcapCells, + ## createHcalBarrelCells, + ## createHcalEndcapCells, + createemptycells, + ## createClusters, + ## createEcalBarrelPositionedCaloClusterCells, + ## correctCaloClusters, + createTopoClusters, + createEcalBarrelPositionedCaloTopoClusterCells, + correctCaloTopoClusters, + out + ], + EvtSel = 'NONE', + EvtMax = Nevts, + ExtSvc = ExtSvc, + StopOnSignal = True, + ) diff --git a/data/elecNoise_ecalBarrelFCCee_theta.root b/data/elecNoise_ecalBarrelFCCee_theta.root new file mode 100644 index 0000000..5279adc Binary files /dev/null and b/data/elecNoise_ecalBarrelFCCee_theta.root differ diff --git a/geometry/README.md b/geometry/README.md index 3a8c500..aa63e40 100644 --- a/geometry/README.md +++ b/geometry/README.md @@ -5,4 +5,9 @@ Upon detector geometry change, the noise values, noise maps and neighbor maps (u - `create_noise_file_chargePreAmp.py` uses the output of `create_capacitance_file.py` to create a TFile with noise vs eta for the different layers. The noise is derived based on [this](https://indico.cern.ch/event/1066234/contributions/4708987/attachments/2387716/4080914/20220209_Brieuc_Francois_Noble_Liquid_Calorimetry_forFCCee_FCCworkshop2022.pdf#page=7). - you then have to go to `../FCCSW_ecal` to run `neighbours.py` and `noise_map.py` (make sure you modify `BarrelNoisePath` in the latter script to point to your new noise values) +If using the theta-module segmentation (with inclined modules and merging of cells along theta and/or module directions) then the following scripts should be used instead: +- `create_capacitance_file_theta.py` +- `create_noise_file_chargePreAmp_theta.py` +- you then have to go to `../FCCSW_ecal` to run `neighbours_theta.py` and `noise_map_theta.py` + NB: At the moment those files are provided from `/eos/user/b/brfranco/rootfile_storage/` (which has restricted writing permissions) as explained in the main README. This should probably be updated to take them from e.g. a public cern box link but meanwhile you can just point `TopoCaloNeighbours` and `TopoCaloNoisyCells` to the local updated files or contact Brieuc so that he updates the files directly there. diff --git a/geometry/create_capacitance_file_theta.py b/geometry/create_capacitance_file_theta.py new file mode 100644 index 0000000..2888eeb --- /dev/null +++ b/geometry/create_capacitance_file_theta.py @@ -0,0 +1,351 @@ +# calculate and save histograms of capacitance per source vs theta +# output is saved in ROOT file with given filename +# +# execute script with +# python create_capacitance_file_theta.py +# + +from ROOT import TH1F, TF1, TF2, TCanvas, TLegend, TFile, gStyle +import ROOT +from math import ceil, sin, cos, atan, exp, log, tan, pi, sqrt, asin, degrees, radians + +ROOT.gROOT.SetBatch(ROOT.kTRUE) + +gStyle.SetPadTickY(1) + +debug = False +verbose = True + +#Dimensions +### FCCee, latest version, with LAr gap size adjusted to obtain 1536 modules (default gives 1545) +# Debug: Number of layers: 12 total thickness 40 +# Info: ECAL cryostat: front: rmin (cm) = 215.4 rmax (cm) = 216.78 dz (cm) = 310 +# Info: ECAL cryostat: back: rmin (cm) = 261.33 rmax (cm) = 271.6 dz (cm) = 310 +# Info: ECAL cryostat: side: rmin (cm) = 216.78 rmax (cm) = 261.33 dz (cm) = 3.38 +# Info: ECAL services: front: rmin (cm) = 216.78 rmax (cm) = 217.28 dz (cm) = 306.62 +# Info: ECAL services: back: rmin (cm) = 257.33 rmax (cm) = 261.33 dz (cm) = 306.62 +# Info: ECAL bath: material = LAr rmin (cm) = 217.28 rmax (cm) = 257.33 thickness in front of ECal (cm) = 0.5 thickness behind ECal (cm) = 4 +# Info: ECAL calorimeter volume rmin (cm) = 217.28 rmax (cm) = 257.33 +# Info: passive inner material = Lead +# and outer material = lArCaloSteel +# thickness of inner part at inner radius (cm) = 0.18 +# thickness of inner part at outer radius (cm) = 0.18 +# thickness of outer part (cm) = 0.01 +# thickness of total (cm) = 0.2 +# rotation angle = 0.872665 +# Info: number of passive plates = 1536 azim. angle difference = 0.00409062 +# Info: distance at inner radius (cm) = 0.888809 +# distance at outer radius (cm) = 1.05264 +# Info: readout material = PCB +# thickness of readout planes (cm) = 0.12 +# number of readout layers = 12 +# Info: thickness of calorimeter (cm) = 40.05 +# length of passive or readout planes (cm) = 56.586 + +filename = "capacitances_perSource_ecalBarrelFCCee_theta.root" +#filename = "capacitances_ecalBarrelFCCee_nLayer_%d_fromAnalytical.root"%numLayers + +# layer 2 require special care as it is separated in several cells and that the shield run beneath the etch: cell 2 signal pad top capa: 0.68 + 0.20 = 0.88, cell 2 signal pad bot: 0.56 + 0.21 = 0.77, cell 3: 0.34 + 2.4 = 2.74, cell 4: 1 + 0.25 = 1.25, cell 5: 1.85 + 0.28 = 2.13 + +# 1*15mm + 11*35mm = 400 +activeTotal = 400.5 +inclinedTotal = 565.86 +tracesPerLayer = [0, 1, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14] # only one trace for strip layer because 4 cells instead of one. Version where we extract all channels from the back +nMergedThetaCells = [4, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4] +nMergedModules = [2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2] +# careful, this is not really the radial spacing, it is, after dilution, the spacing in the parallel direction --> radial depth spacing will not be constant +readoutLayerRadialLengths = [1.500000] * 1 + [3.500000] * 11 +# Detector +rmin = 2172.8 +Nplanes = 1536 +inclination_degree = 50 +angle = inclination_degree / 180. * pi #inclination angle inn radian +passiveThickness = 2.0 #mm +#Segmentation +deltaTheta = 0.009817477/4. +minTheta = 0.58905 +maxTheta = pi-minTheta +numTheta = int(ceil((maxTheta-minTheta)/deltaTheta)) +if verbose: + print("minTheta = ", minTheta) + print("maxTheta = ", maxTheta) + print("numTheta = ", numTheta) + print("Nplanes = ", Nplanes) + print("rmin = ", rmin) + print("activeTotal = ", activeTotal) + print("readoutLayerRadialLengths = ", readoutLayerRadialLengths) + print("inclination (deg) = " , inclination_degree) + print("inclinedTotal = ", inclinedTotal) + print("passiveThickness = ", passiveThickness) + +#PCB dimensions [mm] +hhv = 0.1 +hs = 0.17 +t = 0.035 +w = 0.127 +ws = 0.250 +#hm = 0.250 +hm = 0.2075 +pcbThickness = 7 * t + 2 * hhv + 2 * hs + 2 * hm #mm +print("pcbThickness: %f"%pcbThickness) +#constants: +# distance from signal trace to shield (HS) - from impedance vs. trace width vs. distance to ground layer 2D plot (Z = 50 Ohm) +# trace width (W) - min value +# trace thickness (T) - min value +# distance from shield to the edge of PCB +#http://www.analog.com/media/en/training-seminars/design-handbooks/Basic-Linear-Design/Chapter12.pdf, page 40 +#signal trace +epsilonR = 4.4 # PCB +#conversion factor: 1 inch = 25.4 mm +inch2mm = 25.4 +# capa per length from maxwel1 (pF/mm) +# strip layer has smaller capacitance due to traces running beneath the anti-etch +capa_per_mm = [0.123, 0.062, 0.123, 0.123, 0.123, 0.123, 0.123, 0.123, 0.123, 0.123, 0.123, 0.123] +# multiplicative factors +# for the trace, factor 2 because we have two HV plate / absorber capa per cell +nmultTrace = 2 +# for the shield, where we use maxwell, the extra factor 2 (two signal pad / shield capa) is already accounted for +nmultShield = 1 +# dielectric constants +epsilonRLAr = 1.5 # LAr at 88 K +epsilon0 = 8.854/1000. #pF/mm + + +# Fill the layer length, trace length, etc +readoutLayerParallelLengths = [] +real_radial_separation = [rmin] +real_radial_depth = [] +inclinations_wrt_radial_dir_at_middleRadialDepth = [] +trace_length = [] +numLayers = len(readoutLayerRadialLengths) +dilution_factor = inclinedTotal / activeTotal +trace_length_inner = 0 +trace_length_outer = 0 +outer = False +current_electrode_length = 0 +for idx in range(numLayers):# first pass to get all length parallel to the readout, real radial separation, inclination at the middle of the layer + readoutLayerRadialLengths[idx] *= 10 + parallel_length = readoutLayerRadialLengths[idx] * dilution_factor + # Tricky point: in the xml geo, you define 'radial'segmentation, but these depths will be the one parallel to the plates after scaling by the dilution factor --> even when setting constant radial depth, the geometry builder will make constant parallel length step, not constant radial steps + readoutLayerParallelLengths.append(parallel_length) + if outer: # prepare the starting trace length when starting to extract by the back of the PCB + trace_length_outer += parallel_length + if tracesPerLayer[idx] == 0 and tracesPerLayer[idx - 1] == 0: + outer = True + # sqrt(r**2+(L1+i*L2)**2+2*r*(L1+i*L2)*cos(alpha)) where L1 = 2.68, L2=12.09, r=192, alpha=50) + current_electrode_length += parallel_length + real_radial_separation.append(sqrt(rmin * rmin + current_electrode_length * current_electrode_length + 2 * rmin * current_electrode_length * cos(angle))) + real_radial_depth.append(real_radial_separation[idx+1] - real_radial_separation[idx]) + # treating the fact that radial angle decreases when radial depth increase + # angle comprise by lines from 1) Interaction point to inner right edge of a cell, 2) Interaction point to outer left edge of the considered cell (useful to get the plate angle with radial direction that changes with increasing R) + # based on scalene triangle sine law A/sin(a) = B/sin(b) = C/sin(c) (outer left edge aligned on the Y axis) + inclinations_wrt_radial_dir_at_middleRadialDepth.append(asin(rmin * sin(angle) / (real_radial_separation[idx] + ((real_radial_separation[idx+1] - real_radial_separation[idx]) / 2)))) + + +# second pass to get trace lengths +outer = False +for idx in range(numLayers): + if tracesPerLayer[idx] == 0 and tracesPerLayer[idx - 1] == 0: # we change direction + outer = True + if outer: + trace_length.append(trace_length_outer) + if idx == numLayers - 1: + trace_length_outer == 0 + continue + trace_length_outer -= readoutLayerParallelLengths[idx+1] + else: + trace_length.append(trace_length_inner) + trace_length_inner += readoutLayerParallelLengths[idx] + +print('Readout radial lengths originally asked: ', readoutLayerRadialLengths) +print('Readout parallel lengths: ', readoutLayerParallelLengths) +print("Real radial separation: ", real_radial_separation) +print("Real radial depth: ", real_radial_depth) +print("inclinations_wrt_radial_dir_at_middleRadialDepth: ", [degrees(inclinations) for inclinations in inclinations_wrt_radial_dir_at_middleRadialDepth]) +print("Signal trace length per layer: ", trace_length) + + +gStyle.SetOptStat(0) + +cImpedance = TCanvas("cImpedance","",600,800) +cImpedance.Divide(1,2) +cImpedance.cd(1) +fImpedance = TF2("fImpedance","60/sqrt([0])*log(1.9*(2*x+[1])/(0.8*y+[1]))",0.04,0.2,0.04,0.2) +fImpedance.SetTitle("Impedance vs trace width and distance to ground") +fImpedance.SetParameters(epsilonR, t) +fImpedance.Draw("colz") +fImpedance.GetXaxis().SetTitle("Distance to ground [mm]") +fImpedance.GetYaxis().SetTitle("Trace width [mm]") +cImpedance.cd(2) +fImpedance1D = TF1("fImpedance1D","60/sqrt([0])*log(1.9*(2*x+[1])/(0.8*[2]+[1]))",0.04,0.2) +fImpedance1D.SetTitle("Impedance vs distance to ground") +fImpedance1D.SetParameters(epsilonR, t, w) +fImpedance1D.Draw() +fImpedance1D.GetXaxis().SetTitle("Distance to ground [mm]") +fImpedance1D.GetYaxis().SetTitle("Impedance [#Omega]") + +#prepare the TH1 +hCapTrace = [] +hCapShield = [] +hCapDetector = [] +line_color_number = 1 +line_style_number = 1 +for i in range (0, len(readoutLayerRadialLengths)): + if line_color_number == 8: + line_color_number = 22 + if line_style_number == 8: + line_style_number = 1 + #traces + hCapTrace.append(TH1F()) + hCapTrace[i].SetBins(numTheta, minTheta, maxTheta) + hCapTrace[i].SetLineColor(line_color_number) + hCapTrace[i].SetLineStyle(line_style_number) + hCapTrace[i].SetLineWidth(2) + hCapTrace[i].SetTitle("Stripline capacitance; #theta; Capacitance [pF]") + hCapTrace[i].SetName("hCapacitance_traces"+str(i)) + #shields + hCapShield.append(TH1F()) + hCapShield[i].SetBins(numTheta, minTheta, maxTheta) + hCapShield[i].SetLineColor(line_color_number) + hCapShield[i].SetLineStyle(line_style_number) + hCapShield[i].SetLineWidth(2) + hCapShield[i].SetTitle("Signal pads - ground shields capacitance; #theta; Capacitance [pF]") + hCapShield[i].SetName("hCapacitance_shields"+str(i)) + #area + hCapDetector.append(TH1F()) + hCapDetector[i].SetBins(numTheta, minTheta, maxTheta) + hCapDetector[i].SetLineColor(line_color_number) + hCapDetector[i].SetLineStyle(line_style_number) + hCapDetector[i].SetLineWidth(2) + hCapDetector[i].SetTitle("Signal pad - absorber capacitance; #theta; Capacitance [pF]") + hCapDetector[i].SetName("hCapacitance_detector"+str(i)) + if line_color_number > 8: + line_color_number += 10 + else: + line_color_number += 1 + line_style_number += 1 + +cTrace = TCanvas("cTrace","",600,400) +cShield = TCanvas("cShield","",600,400) +cDetector = TCanvas("cDetector","",600,400) + +legend = TLegend(0.1,0.693,0.8,0.9) +#legend = TLegend(0.135,0.693,0.8,0.892) +legend.SetHeader("Longitudinal layers") +legend.SetNColumns(4) +capa_shield_max = 0 +capa_det_max = 0 +for i in range (0, len(readoutLayerParallelLengths)): + print("--------------") + for index in range(0, numTheta): + theta = minTheta + index * deltaTheta + thetaNext = minTheta + (index+1) * deltaTheta + eta = -log(tan(theta/2.0)) + deltaEta = abs( -log(tan((minTheta + (index+1) * deltaTheta)/2.0)) - eta) + if (debug): + print("theta = ", theta) + print("eta = ", eta) + print("delta eta = ", deltaEta) + + # take into account the inclination in theta + traceLength = trace_length[i] / sin(theta) + #print("Layer %d trace length %f"%(i+1, traceLength)) + #Trace capacitance (stripline) - not used since already accounted for elsewhere + logStripline = log(3.1 * hs / (0.8 * w + t)) + # analytical formula + capacitanceTrace = nmultTrace * 1 / inch2mm * 1.41 * epsilonR / logStripline * traceLength + hCapTrace[i].SetBinContent(index+1, capacitanceTrace) + + #Shield capacitance (microstrip) + cellLength = readoutLayerParallelLengths[i] / sin(theta) + logMicrostrip = log(5.98 * hm / (0.8 * ws + t)) + # analytical formula (nmultShield = 2) + #capacitanceShield = nmultShield * nMergedModules[i] * cellLength * tracesPerLayer[i] * 1 / inch2mm * 0.67 * (epsilonR + 1.41) / logMicrostrip + # from maxwell (nmultShield = 1) + # dont multiply by nMergedThetaCells: the shield/pad capa is reasonably independent of the cell size and the fact that there is some merging + # done for theta cells is already taken into account by the tracesPerLayer + capacitanceShield = nmultShield * nMergedModules[i] * cellLength * tracesPerLayer[i] * capa_per_mm[i] + if capacitanceShield > capa_shield_max: + capa_shield_max = capacitanceShield + hCapShield[i].SetBinContent(index+1, capacitanceShield) + + ##Detector area (C = epsilon*A/d) + #area = ( radius[i] * ( 1 / (tan(2. * atan(exp(- (index + 1) * deltaEta)))) - 1 / (tan(2. * atan(exp(- index * deltaEta))) ) ) + # + radius[i + 1] * ( 1 / (tan(2. * atan(exp(- (index + 1) * deltaEta)))) - 1 / (tan(2. * atan(exp(- index * deltaEta))) ) ) + # ) / 2. * (radius[i+1] - radius[i]) + #distance = (radius[i+1] + radius[i]) / 2. * pi / Nplanes * cos (angle) - pcbThickness / 2. - passiveThickness / 2. + + #Detector area (C = epsilon*A/d) + area = abs( real_radial_separation[i] * ( 1 / tan(thetaNext) - 1 / tan(theta) ) + + real_radial_separation[i + 1] * ( 1 / tan(thetaNext) - 1 / tan(theta)) + ) / 2. * (real_radial_separation[i+1] - real_radial_separation[i]) + # get the cell size perpendicular to the plate direction from the cell size on the circle at given radius and the inclination w.r.t. radial dir, then remove the PCB and lead thickness (no need for any factor here because we are perpendicular to the PCB and lead plates) --> gives the LAr gap size perpendicular + distance = (2 * pi * (real_radial_separation[i+1] + real_radial_separation[i]) / 2. / Nplanes * cos (inclinations_wrt_radial_dir_at_middleRadialDepth[i]) - pcbThickness - passiveThickness) / 2. # divided by two because two lar gap per cell + distance += hhv #the capa is between signal plate and absorber --> need to add distance between HV plate and signal pad + distance += t #the capa is between signal plate and absorber --> need to add distance between HV plate and signal pad + if (abs(theta - pi/2.)<1e-4): + print("LAr gap size (perpendicular) + hhv + t: %f mm"%distance) + capacitanceDetector = nMergedModules[i] * nMergedThetaCells[i] * 2 * epsilon0 * epsilonRLAr * area / distance # factor 2 is because there are 2 LAr gaps for each cell + hCapDetector[i].SetBinContent(index+1, capacitanceDetector) + if capacitanceDetector > capa_det_max: + capa_det_max = capacitanceDetector + if (abs(theta - pi/2.)<1e-4): + print("layer %d" %(i+1), "theta=%f" %theta, ": capacitanceTrace %.0f pF," %capacitanceTrace, "capacitanceShield %.0f pF" %capacitanceShield, "capacitanceDetector %.0f pF" %capacitanceDetector) + #, "distance %.1f mm" %distance + + #Draw + cTrace.cd() + if i == 0: + hCapTrace[i].Draw() + else: + hCapTrace[i].Draw("same") + legend.AddEntry(hCapTrace[i],"layer %d"%(i+1),"l") + cShield.cd() + if i == 0: + hCapShield[i].Draw() + else: + hCapShield[i].Draw("same") + cDetector.cd() + if i == 0: + hCapDetector[i].Draw() + else: + hCapDetector[i].Draw("same") + +maximum = capa_shield_max + +plots = TFile(filename,"RECREATE") + +for i in range (0, len(readoutLayerParallelLengths)): + hCapTrace[i].SetMinimum(0) + hCapTrace[i].SetMaximum(maximum*1.8) + hCapTrace[i].Write() + hCapShield[i].SetMinimum(0) + hCapShield[i].SetMaximum(capa_shield_max*1.5) + hCapShield[i].Write() + hCapDetector[i].SetMinimum(0) + hCapDetector[i].SetMaximum(capa_det_max*1.5) + hCapDetector[i].Write() + +cTrace.cd() +legend.Draw() +cTrace.Update() +cTrace.Write() +cTrace.Print("capa_trace.png") +cTrace.Print("capa_trace.pdf") +cShield.cd() +legend.Draw() +cShield.Update() +cShield.Write() +cShield.Print("capa_shield.png") +cShield.Print("capa_shield.pdf") +cDetector.cd() +legend.Draw() +cDetector.Update() +cDetector.Write() +cDetector.Print("capa_detector.png") +cDetector.Print("capa_detector.pdf") + +fImpedance.Write() +fImpedance1D.Write() + +#closeInput = raw_input("Press ENTER to exit") diff --git a/geometry/create_noise_file_chargePreAmp_theta.py b/geometry/create_noise_file_chargePreAmp_theta.py new file mode 100644 index 0000000..f79307e --- /dev/null +++ b/geometry/create_noise_file_chargePreAmp_theta.py @@ -0,0 +1,348 @@ +# uses the output of create_capacitance_file_theta.py to create +# a TFile with noise vs theta for the different layers. +# The noise is derived based on https://indico.cern.ch/event/1066234/contributions/4708987/attachments/2387716/4080914/20220209_Brieuc_Francois_Noble_Liquid_Calorimetry_forFCCee_FCCworkshop2022.pdf#page=7 + +# execute script with +# python create_noise_file_chargePreAmp_theta.py + +# the input file name is configured with capa_filename +# noise and capacitances per layer are saved in root files +# capacitances_ecalBarrelFCCee_theta.root +# elecNoise_ecalBarrelFCCee_theta.root +# in folder noise_capa_ + +from ROOT import TH1F, TCanvas, TLegend, TFile, gStyle, gPad +import ROOT +import itertools +from datetime import date +import os, sys +from numpy import ones,vstack +from numpy.linalg import lstsq +from math import floor + +ROOT.gROOT.SetBatch(ROOT.kTRUE) +ROOT.gStyle.SetPadTickY(1) +ROOT.gStyle.SetOptTitle(0) + +def rescaleaxis(g, scale = 50e-12/1e-9): + """This function rescales the x-axis on a TGraph.""" + N = g.GetN() + x = g.GetX() + for i in range(N): + x[i] *= scale + g.GetHistogram().Delete() + g.SetHistogram(0) + return + +# you know the noise charge rms from Martin (as a function of capacitance): https://indico.cern.ch/event/1066234/contributions/4708987/attachments/2387716/4080914/20220209_Brieuc_Francois_Noble_Liquid_Calorimetry_forFCCee_FCCworkshop2022.pdf#page=7 +# you need to further get, for each layer, the collected charge corresponding to an energy deposit of 1 MeV in the cell (cell considered including the energy in absorber and PCB), the cell merging strategy does not matter yet here to first approximation because the 1 MeV current equivalent will be the same in any merging scenarios (the current will be shared in more gaps when merging many cells, but all these current will then be 'summed' before to reach the readout). Merging many cells will pay back later, when more signal will be collected per read out channel compared to merging less cells, for the same noise values + +# Retrieving a capa dependent function for the noise charge rms in terms of number of electrons: https://indico.cern.ch/event/1066234/contributions/4708987/attachments/2387716/4080914/20220209_Brieuc_Francois_Noble_Liquid_Calorimetry_forFCCee_FCCworkshop2022.pdf#page=7 +points_capa_noise = [(100, 4375), (500, 6750)] +x_coords, y_coords = zip(*points_capa_noise) +A = vstack([x_coords, ones(len(x_coords))]).T +m, c = lstsq(A, y_coords,rcond=-1)[0] +#print("Line Solution is y = {m}x + {c}".format(m=m, c=c)) +def get_noise_charge_rms(capacitance): + return m * capacitance + c # number of electrons + +# Get the equivalent of 1 MeV energy deposit in a cell (absorber + Lar) in terms of number of electrons in the charge pre-amplifier + shaper +r_recomb = 0.04 +w_lar = 23.6 # eV needed to create a ion/electron pair +def get_ref_charge(SF, E_dep = 1 * pow(10, 6)): #E_dep en eV, choose 1 MeV + return E_dep * SF * (1 - r_recomb) / (2 * w_lar) # nA, the factor 2 comes from: Q_tot = I_0 * t_drift * 1/2 (rectangle --> triangle), t_drift cancels out from the formula to get I_0 which has v_drift/d_gap (Ramo Shockley). + # Assumption: shaping time is similar or bigger to drift time + +# Sampling fraction in each layer +# numbers updated for latest model with 1536 modules +SFfcc = [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 +nLayers = len(SFfcc) + +print("##########: ", get_ref_charge(0.16)) + +SF_rounded_forPrint = [] +for SF in SFfcc: + SF_rounded_forPrint.append(round(SF,2)) +print('SF:', SF_rounded_forPrint) + +#output_folder = "noise_capa_" + date.today().strftime("%y%m%d") +output_folder = "noise_vs_capa_chargePreAmp" + +#filename = "ecalBarrelFCCee_"+flagImpedance+"Ohm_"+flagTraces+"_"+str(flagsShieldsWidth)+"shieldWidth" +capa_filename = "capacitances_perSource_ecalBarrelFCCee_theta.root" +if not os.path.exists(capa_filename): + print("Error: capacitance file does not exist, please run first python create_capacitance_file.py") + sys.exit(1) +fIn = TFile(capa_filename, "r") + +output_folder = "noise_capa_" + date.today().strftime("%y%m%d") +if not os.path.isdir(output_folder): + os.mkdir(output_folder) +fSaveAll = TFile(os.path.join(output_folder, "capacitances_ecalBarrelFCCee_theta.root"),"RECREATE") +fSave = TFile(os.path.join(output_folder, "elecNoise_ecalBarrelFCCee_theta.root"),"RECREATE") + +gStyle.SetOptStat(0) + +#TH1 with capacitances +hCapShield = [] +hCapTrace = [] +hCapDetector = [] +hCapTotal = [] +# electronic noise histograms +h_elecNoise_fcc = [] # default total noise shield + detector capacitance (without trace capacitance) -> to be used in FCCSW as noise estimation +h_elecNoise_all = [] # total noise shield + trace + detector capacitance +h_elecNoise_withTraceCap = [] # total noise without trace capacitance (as in ATLAS - trace cap. can be neglected): from shield + detector capacitance +h_elecNoise_shield = [] +h_elecNoise_trace = [] +h_elecNoise_detector = [] + +#Read graphs from files +for i in range (0, nLayers): + nameShield = "hCapacitance_shields"+str(i) + hCapShield.append(fIn.Get(nameShield)) + nameTrace = "hCapacitance_traces"+str(i) + hCapTrace.append(fIn.Get(nameTrace)) + nameDetector = "hCapacitance_detector"+str(i) + hCapDetector.append(fIn.Get(nameDetector)) + +index = 0 +nbins = hCapShield[index].GetNbinsX() +thetaMin = hCapShield[index].GetXaxis().GetBinLowEdge(1) +thetaMax = hCapShield[index].GetXaxis().GetBinUpEdge(nbins) +print("number of bins ", nbins, ", thetaMin ", thetaMin, ", thetaMax ", thetaMax) + +maximumCap = 0. +maximumNoise = 0. +maximumNoiseWithTrace = 0. + +line_color_number = 1 +line_style_number = 1 +for i in range (0, nLayers): + if line_color_number == 10: + line_color_number = 28 + if line_style_number > 10: + line_style_number = 1 + #Prepare electronic noise histograms + h_elecNoise_fcc.append( TH1F() ) + h_elecNoise_fcc[i].SetLineWidth(3) + h_elecNoise_fcc[i].SetLineColor(line_color_number) + h_elecNoise_fcc[i].SetLineStyle(line_style_number) + h_elecNoise_fcc[i].SetBins(nbins, thetaMin, thetaMax) + h_elecNoise_fcc[i].SetTitle("Default electronic noise: shield + detector capacitance; #theta; Electronic noise [MeV]") + h_elecNoise_fcc[i].SetName("h_elecNoise_fcc_"+str(i+1)) + + h_elecNoise_all.append( TH1F() ) + h_elecNoise_all[i].SetLineWidth(3) + h_elecNoise_all[i].SetLineColor(line_color_number) + h_elecNoise_all[i].SetLineStyle(line_style_number) + h_elecNoise_all[i].SetBins(nbins, thetaMin, thetaMax) + h_elecNoise_all[i].SetTitle("Total electronic noise: shield + trace + detector capacitance; #theta; Electronic noise [MeV]") + h_elecNoise_all[i].SetName("h_elecNoise_all_"+str(i+1)) + + h_elecNoise_withTraceCap.append( TH1F() ) + h_elecNoise_withTraceCap[i].SetLineWidth(3) + h_elecNoise_withTraceCap[i].SetLineColor(line_color_number) + h_elecNoise_withTraceCap[i].SetLineStyle(line_style_number) + h_elecNoise_withTraceCap[i].SetBins(nbins, thetaMin, thetaMax) + h_elecNoise_withTraceCap[i].SetTitle("Electronic noise with trace capacitance; #theta; Electronic noise [MeV]") + h_elecNoise_withTraceCap[i].SetName("h_elecNoise_withTraceCap_"+str(i+1)) + + h_elecNoise_shield.append( TH1F() ) + h_elecNoise_shield[i].SetLineWidth(3) + h_elecNoise_shield[i].SetLineColor(line_color_number) + h_elecNoise_shield[i].SetLineStyle(line_style_number) + h_elecNoise_shield[i].SetBins(nbins, thetaMin, thetaMax) + h_elecNoise_shield[i].SetTitle("Electronic noise - shields; #theta; Electronic noise [MeV]") + h_elecNoise_shield[i].SetName("h_elecNoise_shield_"+str(i+1)) + + h_elecNoise_trace.append( TH1F() ) + h_elecNoise_trace[i].SetLineWidth(3) + h_elecNoise_trace[i].SetLineColor(line_color_number) + h_elecNoise_trace[i].SetLineStyle(line_style_number) + h_elecNoise_trace[i].SetBins(nbins, thetaMin, thetaMax) + h_elecNoise_trace[i].SetTitle("Electronic noise -traces; #theta; Electronic noise [MeV]") + h_elecNoise_trace[i].SetName("h_elecNoise_trace_"+str(i+1)) + + h_elecNoise_detector.append( TH1F() ) + h_elecNoise_detector[i].SetLineWidth(3) + h_elecNoise_detector[i].SetLineColor(line_color_number) + h_elecNoise_detector[i].SetLineStyle(line_style_number) + h_elecNoise_detector[i].SetBins(nbins, thetaMin, thetaMax) + h_elecNoise_detector[i].SetTitle("Electronic noise - detector; #theta; Electronic noise [MeV]") + h_elecNoise_detector[i].SetName("h_elecNoise_detector_"+str(i+1)) + + #Total capacitance plot (shield + trace + detector) + hCapTotal.append( TH1F() ) + hCapTotal[i].SetBins(nbins, thetaMin, thetaMax) + hCapTotal[i].SetLineColor(line_color_number) + hCapTotal[i].SetLineStyle(line_style_number) + hCapTotal[i].SetLineWidth(3) + hCapTotal[i].SetTitle("Total capacitance; #theta; Capacitance [pF]") + hCapTotal[i].SetName("hCapacitance"+str(i)) + + ref_charge_1mev = get_ref_charge(SFfcc[i]) + if i != 0: + print(noise_charge_rms) + print(i, " ", ref_charge_1mev) + for ibin in range(0, nbins+1): + capShield = hCapShield[i].GetBinContent(ibin) + capTrace = hCapTrace[i].GetBinContent(ibin) + capDetector = hCapDetector[i].GetBinContent(ibin) + #total capacitance + #hCapTotal[i].SetBinContent( ibin, capShield + capTrace + capDetector ) + hCapTotal[i].SetBinContent( ibin, capShield + capDetector ) + # Compute the NOISE + noise_charge_rms = get_noise_charge_rms(capShield + capDetector) + noise = noise_charge_rms / (ref_charge_1mev) + noiseWithTrace = get_noise_charge_rms(capShield + capDetector + capTrace) / ref_charge_1mev + noiseShield = get_noise_charge_rms(capShield) / ref_charge_1mev + noiseTrace = get_noise_charge_rms(capTrace) / ref_charge_1mev + noiseDetector = get_noise_charge_rms(capDetector) / ref_charge_1mev + #find maximum for drawing of histograms + if noise > maximumNoise: + maximumNoise = noise + if noiseWithTrace > maximumNoiseWithTrace: + maximumNoiseWithTrace = noiseWithTrace + if (capShield + capDetector)>maximumCap: + maximumCap = (capShield + capDetector) + + #fill histogram + #default: without traces + h_elecNoise_fcc[i].SetBinContent(ibin, noise) + h_elecNoise_all[i].SetBinContent(ibin, noise) + h_elecNoise_withTraceCap[i].SetBinContent(ibin, noiseWithTrace) + h_elecNoise_shield[i].SetBinContent(ibin, noiseShield) + h_elecNoise_trace[i].SetBinContent(ibin, noiseTrace) + h_elecNoise_detector[i].SetBinContent(ibin, noiseDetector) + if ibin==floor((nbins+1)/2): + print("layer %d" %(i+1), "eta==0: capacitance %.0f pF," %( capShield + capTrace + capDetector ), "total elec. noise %.4f MeV" %noise, "elec. noise without trace cap. %.4f MeV" %noiseWithTrace) + line_color_number += 1 + line_style_number += 1 + +#print maximumCap, maximumNoise, maximumNoiseWithTrace + +cCapacitance = TCanvas("cCapacitance","Capacitance per cell",800,600) +cCapacitance.cd() +#legend = TLegend(0.135,0.573,0.466,0.872) +legend = TLegend(0.135,0.693,0.8,0.892) +legend.SetBorderSize(0) +legend.SetLineColor(0) +legend.SetLineStyle(0) +legend.SetLineWidth(0) +legend.SetFillColor(0) +legend.SetFillStyle(0) +legend.SetHeader("Longitudinal layers") +legend.SetNColumns(4) +for i, h in enumerate(hCapTotal): + h.SetMinimum(0) + h.SetMaximum(maximumCap*1.5) + h.GetYaxis().SetTitleOffset(1.4) + if i == 0: + h.Draw("") + else: + h.Draw("same") + legend.AddEntry(h,"Layer " + str(i+1),"l") + + +# Prepare "nice" plots & save all capacitances + noise +fSaveAll.cd() + +for h in hCapTotal: + h.SetMinimum(0.) + #h.SetMaximum(maximumCap*1.3) + h.GetYaxis().SetTitleOffset(1.4) + h.Write() + +legend.Draw() +cCapacitance.Update() +cCapacitance.Write() +cCapacitance.Print(os.path.join(output_folder, "cCapacitance.png")) + +#maximumCap = 1200. +#maximumNoise = 0.04 +#maximumNoiseWithTrace = 0.015 + +for h in itertools.chain(h_elecNoise_fcc, h_elecNoise_all): + h.SetMinimum(0.) + h.SetMaximum(maximumNoise*1.5) + h.GetYaxis().SetTitleOffset(1.4) + h.Write() + +for h in itertools.chain(h_elecNoise_withTraceCap, h_elecNoise_shield, h_elecNoise_trace, h_elecNoise_detector): + h.SetMinimum(0.) + h.SetMaximum(maximumNoiseWithTrace*1.5) + h.GetYaxis().SetTitleOffset(1.4) + h.Write() + +cNoise = TCanvas("cNoise","Electronic noise per cell",800,600) +cNoise.cd() +for i, h in enumerate(h_elecNoise_fcc): + if i == 0: + h.Draw("") + else: + h.Draw("same") + +legend.Draw() +cNoise.Update() +cNoise.Write() +cNoise.Print(os.path.join(output_folder, "cNoise.png")) + +cNoiseWithTrace = TCanvas("cNoiseWithTrace","Electronic noise with trace cap. per cell",800,600) +cNoiseWithTrace.cd() +for i, h in enumerate(h_elecNoise_withTraceCap): + if i == 0: + h.Draw("") + else: + h.Draw("same") + +legend.Draw() +cNoiseWithTrace.Update() +cNoiseWithTrace.Write() +cNoiseWithTrace.Print(os.path.join(output_folder, "cNoiseWithTrace.png")) + +legendP = TLegend(0.1,0.6,0.43,0.9) +legendP.SetHeader("Capacitance") + +cCapParts = TCanvas("cCapParts","",1200,1000) +cCapParts.Divide(3,4) +for i in range (0, nLayers): + cCapParts.cd(i+1) + if i < 7: + hCapTotal[i].SetMaximum(100) + else: + hCapTotal[i].SetMaximum(200) + + hCapTotal[i].SetTitle("Layer "+str(i+1)) + hCapTotal[i].GetXaxis().SetTitleSize(0.045) + hCapTotal[i].GetYaxis().SetTitleSize(0.045) + hCapTotal[i].GetYaxis().SetTitleOffset(1.15) + hCapTotal[i].SetLineColor(1) + hCapTotal[i].SetLineStyle(1) + hCapShield[i].SetLineColor(2) + hCapShield[i].SetLineStyle(2) + hCapTrace[i].SetLineColor(3) + hCapTrace[i].SetLineStyle(3) + hCapDetector[i].SetLineColor(4) + hCapDetector[i].SetLineStyle(4) + hCapTotal[i].Draw() + hCapShield[i].Draw("same") + #hCapTrace[i].Draw("same") + hCapDetector[i].Draw("same") + gPad.Update() + if i==0: + legendP.AddEntry(hCapTotal[i],"total cap.","l") + legendP.AddEntry(hCapShield[i],"shield cap.","l") + #legendP.AddEntry(hCapTrace[i],"trace cap.","l") + legendP.AddEntry(hCapDetector[i],"detector cap.","l") + legendP.Draw() + +cCapParts.Write() +cCapParts.Print(os.path.join(output_folder, "cCapParts.png")) + +#Save final noise plot (to be used in FCCSW) +fSave.cd() +for h in h_elecNoise_fcc: + h.Write() + +#closeInput = raw_input("Press ENTER to exit") +