From 29a43948bbb9350beb215f85e15cc94c63a20b7f Mon Sep 17 00:00:00 2001 From: Valentin Volkl Date: Tue, 19 Oct 2021 19:09:59 +0200 Subject: [PATCH] fixes for new nightlies --- .gitignore | 25 + fast-sim-and-analysis/FccFastSimGeneration.md | 2 +- .../starterkit/FccFastSimAnalysis/doPlots.py | 451 ++++++++++++++++++ .../FccCaloPerformance/FccCaloPerformance.md | 41 +- .../FccCaloPerformance/runCaloSim.py | 194 ++++++++ ...loSystem_ReconstructionSW_noiseFromFile.py | 138 ++++++ 6 files changed, 829 insertions(+), 22 deletions(-) create mode 100755 fast-sim-and-analysis/fccanalyses/doc/starterkit/FccFastSimAnalysis/doPlots.py create mode 100644 full-detector-simulations/FccCaloPerformance/runCaloSim.py create mode 100644 full-detector-simulations/FccCaloPerformance/runFullCaloSystem_ReconstructionSW_noiseFromFile.py diff --git a/.gitignore b/.gitignore index 8650c6aa..6519e01f 100644 --- a/.gitignore +++ b/.gitignore @@ -1,9 +1,34 @@ _book + __pycache__ plots .ipynb_checkpoints *.ipynb +*.input +*.output +*-test.md +*.lhe +*.mod +*.cmd* +semaphore +*.sin* +*.la +*.makefile +*.phs +*.vgx +*.evx +iniseed +.libs +*.LHE +*.hst +option.py +*.vg +*.f90 +*.lo +*.gdml +*.tcl* + *.root *.log *.out diff --git a/fast-sim-and-analysis/FccFastSimGeneration.md b/fast-sim-and-analysis/FccFastSimGeneration.md index 753034bb..56163d37 100644 --- a/fast-sim-and-analysis/FccFastSimGeneration.md +++ b/fast-sim-and-analysis/FccFastSimGeneration.md @@ -165,7 +165,7 @@ which mg5_aMC Herwig is available as standalone program: -```bash +``` which Herwig ``` diff --git a/fast-sim-and-analysis/fccanalyses/doc/starterkit/FccFastSimAnalysis/doPlots.py b/fast-sim-and-analysis/fccanalyses/doc/starterkit/FccFastSimAnalysis/doPlots.py new file mode 100755 index 00000000..de28e032 --- /dev/null +++ b/fast-sim-and-analysis/fccanalyses/doc/starterkit/FccFastSimAnalysis/doPlots.py @@ -0,0 +1,451 @@ +#!/usr/bin/env python +import sys, os +import os.path +import ntpath +import importlib +import ROOT +import copy +import re + +#__________________________________________________________ +def removekey(d, key): + r = dict(d) + del r[key] + return r + +#__________________________________________________________ +def mapHistos(var, label, sel, param): + print ('run plots for var:{} label:{} selection:{}'.format(var,label,sel)) + signal=param.plots[label]['signal'] + backgrounds=param.plots[label]['backgrounds'] + + hsignal = {} + for s in signal: + hsignal[s]=[] + for f in signal[s]: + fin=param.inputDir+f+'_'+sel+'_histo.root' + if not os.path.isfile(fin): + print ('file {} does not exist, skip'.format(fin)) + else: + tf=ROOT.TFile(fin) + h=tf.Get(var) + hh = copy.deepcopy(h) + scaleSig=1. + try: + scaleSig=param.scaleSig + except AttributeError: + print ('no scale signal, using 1') + param.scaleSig=scaleSig + print ('scaleSig ',scaleSig) + hh.Scale(param.intLumi*scaleSig) + + if len(hsignal[s])==0: + hsignal[s].append(hh) + else: + hh.Add(hsignal[s][0]) + hsignal[s][0]=hh + + + hbackgrounds = {} + for b in backgrounds: + hbackgrounds[b]=[] + for f in backgrounds[b]: + fin=param.inputDir+f+'_'+sel+'_histo.root' + if not os.path.isfile(fin): + print ('file {} does not exist, skip'.format(fin)) + else: + tf=ROOT.TFile(fin) + h=tf.Get(var) + hh = copy.deepcopy(h) + hh.Scale(param.intLumi) + if len(hbackgrounds[b])==0: + hbackgrounds[b].append(hh) + else: + hh.Add(hbackgrounds[b][0]) + hbackgrounds[b][0]=hh + + for s in hsignal: + if len(hsignal[s])==0: + hsignal=removekey(hsignal,s) + + for b in hbackgrounds: + if len(hbackgrounds[b])==0: + hbackgrounds=removekey(hbackgrounds,b) + + return hsignal,hbackgrounds + +#__________________________________________________________ +def runPlots(var,sel,param,hsignal,hbackgrounds,extralab): + legsize = 0.04*(len(hbackgrounds)+len(hsignal)) + legCoord=[0.68, 0.86-legsize, 0.96, 0.88] + try: + legCoord=param.legendCoord + except AttributeError: + print ('no legCoord, using default one...') + legCoord=[0.68, 0.86-legsize, 0.96, 0.88] + + leg = ROOT.TLegend(legCoord[0],legCoord[1],legCoord[2],legCoord[3]) + leg.SetFillColor(0) + leg.SetFillStyle(0) + leg.SetLineColor(0) + leg.SetShadowColor(10) + leg.SetTextSize(0.035) + leg.SetTextFont(42) + + for s in hsignal: + leg.AddEntry(hsignal[s][0],param.legend[s],"l") + for b in hbackgrounds: + leg.AddEntry(hbackgrounds[b][0],param.legend[b],"f") + + + yields={} + for s in hsignal: + yields[s]=[param.legend[s],hsignal[s][0].Integral(0,-1), hsignal[s][0].GetEntries()] + for b in hbackgrounds: + yields[b]=[param.legend[b],hbackgrounds[b][0].Integral(0,-1), hbackgrounds[b][0].GetEntries()] + + histos=[] + colors=[] + + nsig=len(hsignal) + + for s in hsignal: + histos.append(hsignal[s][0]) + colors.append(param.colors[s]) + + for b in hbackgrounds: + histos.append(hbackgrounds[b][0]) + colors.append(param.colors[b]) + + intLumiab = param.intLumi/1e+06 + + + lt = "FCCAnalyses: FCC-hh Simulation (Delphes)" + rt = "#sqrt{{s}} = {:.1f} TeV, L = {:.0f} ab^{{-1}}".format(param.energy,intLumiab) + + if 'ee' in param.collider: + lt = "FCCAnalyses: FCC-ee Simulation (Delphes)" + rt = "#sqrt{{s}} = {:.1f} GeV, L = {:.0f} ab^{{-1}}".format(param.energy,intLumiab) + + customLabel="" + try: + customLabel=param.customLabel + except AttributeError: + print ('no customLable, using nothing...') + + if 'AAAyields' in var: + drawStack(var, 'events', leg, lt, rt, param.formats, param.outdir+"/"+sel, False , True , histos, colors, param.ana_tex, extralab, param.scaleSig, customLabel, nsig, yields) + return + + if 'stack' in param.stacksig: + if 'lin' in param.yaxis: + drawStack(var+"_stack_lin", 'events', leg, lt, rt, param.formats, param.outdir+"/"+sel, False , True , histos, colors, param.ana_tex, extralab, param.scaleSig, customLabel, nsig) + if 'log' in param.yaxis: + drawStack(var+"_stack_log", 'events', leg, lt, rt, param.formats, param.outdir+"/"+sel, True , True , histos, colors, param.ana_tex, extralab, param.scaleSig, customLabel, nsig) + if 'lin' not in param.yaxis and 'log' not in param.yaxis: + print ('unrecognised option in formats, should be [\'lin\',\'log\']'.format(param.formats)) + + if 'nostack' in param.stacksig: + if 'lin' in param.yaxis: + drawStack(var+"_nostack_lin", 'events', leg, lt, rt, param.formats, param.outdir+"/"+sel, False , False , histos, colors, param.ana_tex, extralab, param.scaleSig, customLabel, nsig) + if 'log' in param.yaxis: + drawStack(var+"_nostack_log", 'events', leg, lt, rt, param.formats, param.outdir+"/"+sel, True , False , histos, colors, param.ana_tex, extralab, param.scaleSig, customLabel, nsig) + if 'lin' not in param.yaxis and 'log' not in param.yaxis: + print ('unrecognised option in formats, should be [\'lin\',\'log\']'.format(param.formats)) + if 'stack' not in param.stacksig and 'nostack' not in param.stacksig: + print ('unrecognised option in stacksig, should be [\'stack\',\'nostack\']'.format(param.formats)) + + +#_____________________________________________________________________________________________________________ +def drawStack(name, ylabel, legend, leftText, rightText, formats, directory, logY, stacksig, histos, colors, ana_tex, extralab, scaleSig, customLabel, nsig, yields=None): + + canvas = ROOT.TCanvas(name, name, 600, 600) + canvas.SetLogy(logY) + canvas.SetTicks(1,1) + canvas.SetLeftMargin(0.14) + canvas.SetRightMargin(0.08) + + + # first retrieve maximum + sumhistos = histos[0].Clone() + iterh = iter(histos) + next(iterh) + + unit = 'GeV' + if 'TeV' in str(histos[1].GetXaxis().GetTitle()): + unit = 'TeV' + + if unit in str(histos[1].GetXaxis().GetTitle()): + bwidth=sumhistos.GetBinWidth(1) + if bwidth.is_integer(): + ylabel+=' / {} {}'.format(int(bwidth), unit) + else: + ylabel+=' / {:.2f} {}'.format(bwidth, unit) + + for h in iterh: + sumhistos.Add(h) + + maxh = sumhistos.GetMaximum() + minh = sumhistos.GetMinimum() + + if logY: + canvas.SetLogy(1) + + # define stacked histo + hStack = ROOT.THStack("hstack","") + hStackBkg = ROOT.THStack("hstackbkg","") + + # first plot backgrounds + histos[nsig].SetLineWidth(0) + histos[nsig].SetLineColor(ROOT.kBlack) + histos[nsig].SetFillColor(colors[nsig]) + hStack.Add(histos[nsig]) + hStackBkg.Add(histos[nsig]) + + # now loop over other background (skipping first) + iterh = iter(histos) + for i in range(nsig): + next(iterh) + next(iterh) + + k = nsig+1 + for h in iterh: + h.SetLineWidth(0) + h.SetLineColor(ROOT.kBlack) + h.SetFillColor(colors[k]) + hStack.Add(h) + hStackBkg.Add(h) + k += 1 + + if not stacksig: + hStack.Draw("hist") + + # define stacked signal histo + hStackSig = ROOT.THStack("hstacksig","") + + # finally add signal on top + for l in range(nsig): + histos[l].SetLineWidth(3) + histos[l].SetLineColor(colors[l]) + if stacksig: + hStack.Add(histos[l]) + else: + hStackSig.Add(histos[l]) + + if stacksig: + hStack.Draw("hist") + + xlabel = histos[1].GetXaxis().GetTitle() + + hStack.GetXaxis().SetTitle(xlabel) + hStack.GetYaxis().SetTitle(ylabel) + + hStack.GetYaxis().SetTitleOffset(1.95) + hStack.GetXaxis().SetTitleOffset(1.40) + + lowY=0. + if logY: + highY=200.*maxh/ROOT.gPad.GetUymax() + threshold=0.5 + bin_width=hStack.GetXaxis().GetBinWidth(1) + lowY=threshold*bin_width + hStack.SetMaximum(highY) + hStack.SetMinimum(lowY) + else: + hStack.SetMaximum(1.3*maxh) + hStack.SetMinimum(0.) + + escape_scale_Xaxis=True + hStacklast = hStack.GetStack().Last() + lowX_is0=True + lowX=hStacklast.GetBinCenter(1)-(hStacklast.GetBinWidth(1)/2.) + highX_ismax=False + highX=hStacklast.GetBinCenter(hStacklast.GetNbinsX())+(hStacklast.GetBinWidth(1)/2.) + + if escape_scale_Xaxis==False: + for i_bin in range( 1, hStacklast.GetNbinsX()+1 ): + bkg_val=hStacklast.GetBinContent(i_bin) + sig_val=histos[0].GetBinContent(i_bin) + if bkg_val/maxh>0.1 and i_bin<15 and lowX_is0==True : + lowX_is0=False + lowX=hStacklast.GetBinCenter(i_bin)-(hStacklast.GetBinWidth(i_bin)/2.) + + val_to_compare=bkg_val + if sig_val>bkg_val : val_to_compare=sig_val + if val_to_compare15 and highX_ismax==False: + highX_ismax=True + highX=hStacklast.GetBinCenter(i_bin)+(hStacklast.GetBinWidth(i_bin)/2.) + highX*=1.1 + # protections + if lowXhStacklast.GetBinCenter(hStacklast.GetNbinsX())+(hStacklast.GetBinWidth(1)/2.) : + highX=hStacklast.GetBinCenter(hStacklast.GetNbinsX())+(hStacklast.GetBinWidth(1)/2.) + if lowX>=highX : + lowX=hStacklast.GetBinCenter(1)-(hStacklast.GetBinWidth(1)/2.) + highX=hStacklast.GetBinCenter(hStacklast.GetNbinsX())+(hStacklast.GetBinWidth(1)/2.) + hStack.GetXaxis().SetLimits(int(lowX),int(highX)) + + if not stacksig: + if 'AAAyields' not in name: hStackSig.Draw("same hist nostack") + + legend.Draw() + + Text = ROOT.TLatex() + Text.SetNDC() + Text.SetTextAlign(31); + Text.SetTextSize(0.04) + + text = '#it{' + leftText +'}' + Text.DrawLatex(0.90, 0.94, text) + + text = '#it{'+customLabel+'}' + Text.SetTextAlign(12); + Text.SetNDC(ROOT.kTRUE) + Text.SetTextSize(0.04) + Text.DrawLatex(0.18, 0.85, text) + + rightText = re.split(",", rightText) + text = '#bf{#it{' + rightText[0] +'}}' + + Text.SetTextAlign(12); + Text.SetNDC(ROOT.kTRUE) + Text.SetTextSize(0.04) + Text.DrawLatex(0.18, 0.81, text) + + rightText[1]=rightText[1].replace(" ","") + text = '#bf{#it{' + rightText[1] +'}}' + Text.SetTextSize(0.035) + Text.DrawLatex(0.18, 0.76, text) + + text = '#bf{#it{' + ana_tex +'}}' + Text.SetTextSize(0.04) + Text.DrawLatex(0.18, 0.71, text) + + text = '#bf{#it{' + extralab +'}}' + Text.SetTextSize(0.025) + Text.DrawLatex(0.18, 0.66, text) + + text = '#bf{#it{' + 'Signal scale=' + str(scaleSig)+'}}' + Text.SetTextSize(0.025) + if scaleSig!=1:Text.DrawLatex(0.18, 0.63, text) + + canvas.RedrawAxis() + canvas.GetFrame().SetBorderSize( 12 ) + canvas.Modified() + canvas.Update() + + if 'AAAyields' in name: + dummyh=ROOT.TH1F("","",1,0,1) + dummyh.SetStats(0) + dummyh.GetXaxis().SetLabelOffset(999) + dummyh.GetXaxis().SetLabelSize(0) + dummyh.GetYaxis().SetLabelOffset(999) + dummyh.GetYaxis().SetLabelSize(0) + dummyh.Draw("AH") + legend.Draw() + + Text.SetNDC() + Text.SetTextAlign(31); + Text.SetTextSize(0.04) + + text = '#it{' + leftText +'}' + Text.DrawLatex(0.90, 0.92, text) + + text = '#bf{#it{' + rightText[0] +'}}' + Text.SetTextAlign(12); + Text.SetNDC(ROOT.kTRUE) + Text.SetTextSize(0.04) + Text.DrawLatex(0.18, 0.83, text) + + text = '#bf{#it{' + rightText[1] +'}}' + Text.SetTextSize(0.035) + Text.DrawLatex(0.18, 0.78, text) + + text = '#bf{#it{' + ana_tex +'}}' + Text.SetTextSize(0.04) + Text.DrawLatex(0.18, 0.73, text) + + text = '#bf{#it{' + extralab +'}}' + Text.SetTextSize(0.025) + Text.DrawLatex(0.18, 0.68, text) + + text = '#bf{#it{' + 'Signal scale=' + str(scaleSig)+'}}' + Text.SetTextSize(0.04) + Text.DrawLatex(0.18, 0.55, text) + + dy=0 + text = '#bf{#it{' + 'Process' +'}}' + Text.SetTextSize(0.035) + Text.DrawLatex(0.18, 0.45, text) + + text = '#bf{#it{' + 'Yields' +'}}' + Text.SetTextSize(0.035) + Text.DrawLatex(0.5, 0.45, text) + + text = '#bf{#it{' + 'Raw MC' +'}}' + Text.SetTextSize(0.035) + Text.DrawLatex(0.75, 0.45, text) + + for y in yields: + text = '#bf{#it{' + yields[y][0] +'}}' + Text.SetTextSize(0.035) + Text.DrawLatex(0.18, 0.4-dy*0.05, text) + + stry=str(yields[y][1]) + stry=stry.split('.')[0] + text = '#bf{#it{' + stry +'}}' + Text.SetTextSize(0.035) + Text.DrawLatex(0.5, 0.4-dy*0.05, text) + + stry=str(yields[y][2]) + stry=stry.split('.')[0] + text = '#bf{#it{' + stry +'}}' + Text.SetTextSize(0.035) + Text.DrawLatex(0.75, 0.4-dy*0.05, text) + + + dy+=1 + #canvas.Modified() + #canvas.Update() + + + printCanvas(canvas, name, formats, directory) + + + + +#____________________________________________________ +def printCanvas(canvas, name, formats, directory): + + if format != "": + if not os.path.exists(directory) : + os.system("mkdir -p "+directory) + for f in formats: + outFile = os.path.join(directory, name) + "." + f + canvas.SaveAs(outFile) + + + +#__________________________________________________________ +if __name__=="__main__": + ROOT.gROOT.SetBatch(True) + ROOT.gErrorIgnoreLevel = ROOT.kWarning + # param file + paramFile = sys.argv[1] + + module_path = os.path.abspath(paramFile) + module_dir = os.path.dirname(module_path) + base_name = os.path.splitext(ntpath.basename(paramFile))[0] + + sys.path.insert(0, module_dir) + param = importlib.import_module(base_name) + + counter=0 + for var in param.variables: + for label, sels in param.selections.items(): + for sel in sels: + hsignal,hbackgrounds=mapHistos(var,label,sel, param) + runPlots(var+"_"+label,sel,param,hsignal,hbackgrounds,param.extralabel[sel]) + if counter==0: runPlots("AAAyields_"+label,sel,param,hsignal,hbackgrounds,param.extralabel[sel]) + counter+=1 diff --git a/full-detector-simulations/FccCaloPerformance/FccCaloPerformance.md b/full-detector-simulations/FccCaloPerformance/FccCaloPerformance.md index 84f2a696..07c4001e 100644 --- a/full-detector-simulations/FccCaloPerformance/FccCaloPerformance.md +++ b/full-detector-simulations/FccCaloPerformance/FccCaloPerformance.md @@ -68,7 +68,7 @@ geantservice = SimG4Svc("SimG4Svc") from Configurables import GeoToGdmlDumpSvc geodumpservice = GeoToGdmlDumpSvc("GeoDump") -geodumpservice.gdml="FCCee_IDEA.gdml" +geodumpservice.gdml="DetFCCeeCLD.gdml" from Configurables import ApplicationMgr ApplicationMgr( TopAlg = [], @@ -84,25 +84,28 @@ ApplicationMgr( TopAlg = [], A job with this configuration can be executed with ```bash -fccrun dumpGeo_fccee.py +# if you want to use your local job options file: +#fccrun dumpGeo_fccee.py + +fccrun $FCCSW/Examples/options/export_detector_gdml.py ``` Note the printout of the GeoSvc and make sure the information is as expected. If there is something unclear or missing make sure to open an [issue](https://github.com/HEP-FCC/FCCSW/issues)! Take a look at the resulting gdml file. Although it is text-based it is not really human-readable for a geometry of this size, but you can check the number of lines and volume names if you are familiar with the geometry. ```bash -tail FCCee_IDEA.gdml +tail DetFCCeeCLD.gdml ``` ```bash # count occurences of "physvol" -grep -c ">h_GenParticles_P") +events.Draw("sqrt(pow(GenParticles.momentum.x,2) + pow(GenParticles.momentum.y,2) +pow(GenParticles.momentum.z,2))>>h_GenParticles_P") c.Draw() ``` @@ -161,7 +163,7 @@ events = f.Get("events") c = ROOT.TCanvas("c_ECalBarrelPositions_xy", "", 700, 600) # draw hits for first five events -events.Draw("ECalBarrelPositions.position.y:ECalBarrelPositions.position.x", "", "", 5, 0) +events.Draw("ECalBarrelPositionedCells.position.y:ECalBarrelPositionedCells.position.x", "", "", 5, 0) c.Draw() ``` @@ -171,11 +173,8 @@ c.Draw() Now that the detector response is simulated, it is time to reconstruct the signals. FCCSW includes another configuration to run a Sliding Window reconstruction: ```bash -fccrun $FCCSWBASEDIR/share/FCCSW/RecFCCeeCalorimeter/options/runFullCaloSystem_ReconstructionSW_noiseFromFile.py \ - -n 100 \ - --input fccee_idea_LAr_pgun.root \ - --noiseFileName http://fccsw.web.cern.ch/fccsw/testsamples/elecNoise_ecalBarrelFCCee_50Ohm_traces1_4shieldWidth.root \ - --filename output_allCalo_reco_noise.root + +k4run $K4RECCALORIMETER/RecFCCeeCalorimeter/tests/options/runFullCaloSystem_ReconstructionSW_noiseFromFile.py -n 50 --input.EventDataSvc fccee_idea_LAr_pgun.root --noiseFileName http://fccsw.web.cern.ch/fccsw/testsamples/elecNoise_ecalBarrelFCCee_50Ohm_traces1_4shieldWidth.root --filename output_allCalo_reco_noise.root ``` This configuration inludes electronics noise especially calculated for this detector geometry. which is overlayed on the branch `ECalBarrelCells` containing information on all cells in the ECal Barrel. @@ -189,7 +188,7 @@ events = f.Get("events") c = ROOT.TCanvas("c_ECalBarrelCellsNoise_energy", "", 700, 600) h = ROOT.TH1F("h_ECalBarrelCells_energy", ";ECal Barrel Cells Energy [GeV]; Cells", 80, -0.2 ,1) -events.Draw("ECalBarrelCells.core.energy >> h_ECalBarrelCells_energy", "", "", 10, 0) +events.Draw("ECalBarrelCells.energy >> h_ECalBarrelCells_energy", "", "", 10, 0) h.GetYaxis().SetRangeUser(0.2, 8*10**6) @@ -201,7 +200,7 @@ c2 = ROOT.TCanvas("c_ECalBarrelCells_energy", "", 700, 600) #h2 = ROOT.TH1F("h_ECalBarrelCellsNoise_energy", ";ECall Barrel Cells Energy with Noise [GeV]; Events", 80, -0.2 ,1) h2 = h.Clone("h_ECalBarrelCellsNoise_energy") h2.SetTitle(";ECal Barrel Cells Energy with Noise [GeV]; Cells") -events.Draw("ECalBarrelCellsNoise.core.energy>> h_ECalBarrelCellsNoise_energy", "", "", 10, 0) +events.Draw("ECalBarrelCellsNoise.energy>> h_ECalBarrelCellsNoise_energy", "", "", 10, 0) h2.GetYaxis().SetRangeUser(0.2, 8*10**6) h2.SetLineColor(ROOT.kBlack) @@ -224,7 +223,7 @@ events = f.Get("events") c = ROOT.TCanvas("c_CaloClusters_energy", "", 700, 600) hEn = ROOT.TH1F("h_CaloClusters_energy", ";ECal Calo Cluster Energy [GeV]; # Clusters", 120, 0 ,8) -events.Draw("CaloClusters.core.energy >> h_CaloClusters_energy") +events.Draw("CaloClustersFinal.energy >> h_CaloClusters_energy") c.Draw() @@ -240,7 +239,7 @@ events = f.Get("events") c = ROOT.TCanvas("c_CaloClusters_energyFit", "", 700, 600) hEn = ROOT.TH1F("h_CaloClusters_energy", ";ECal Calo Cluster Energy [GeV]; # Clusters", 120, 0 ,8) -events.Draw("CaloClusters.core.energy >> h_CaloClusters_energy") +events.Draw("CaloClustersFinal.energy >> h_CaloClusters_energy") fitPre = ROOT.TF1("fitPre","gaus", hEn.GetMean() - 1. * hEn.GetRMS(), hEn.GetMean() + 1. * hEn.GetRMS()) resultPre = hEn.Fit(fitPre, "SQRN") @@ -250,8 +249,8 @@ mean = result.Get().Parameter(1) sigma = result.Get().Parameter(2) dMean = result.Get().Error(1) dSigma = result.Get().Error(2) -print "mean:", round(mean,2), "[GeV]" -print "sigma:", round(sigma ,2), "[GeV]" +print("mean:", round(mean,2), "[GeV]") +print("sigma:", round(sigma ,2), "[GeV]") fit.Draw("SAME") c.Draw() diff --git a/full-detector-simulations/FccCaloPerformance/runCaloSim.py b/full-detector-simulations/FccCaloPerformance/runCaloSim.py new file mode 100644 index 00000000..e1011e69 --- /dev/null +++ b/full-detector-simulations/FccCaloPerformance/runCaloSim.py @@ -0,0 +1,194 @@ +import os + +from GaudiKernel.SystemOfUnits import MeV, GeV, tesla + +# Input for simulations (momentum is expected in GeV!) +momentum = 100 +# theta from 80 to 100 degrees corresponds to -0.17 < eta < 0.17 +thetaMin = 80. +thetaMax = 100. +magneticField = False + +from Gaudi.Configuration import * + +from Configurables import FCCDataSvc +podioevent = FCCDataSvc("EventDataSvc") + +################## Particle gun setup +_pi = 3.14159 + +from Configurables import MomentumRangeParticleGun +pgun = MomentumRangeParticleGun("ParticleGun_Electron") +pgun.PdgCodes = [11] +pgun.MomentumMin = momentum * GeV +pgun.MomentumMax = momentum * GeV +pgun.PhiMin = 0 +pgun.PhiMax = 2 * _pi +pgun.ThetaMin = thetaMin * _pi / 180. +pgun.ThetaMax = thetaMax * _pi / 180. + +from Configurables import GenAlg +genalg_pgun = GenAlg() +genalg_pgun.SignalProvider = pgun +genalg_pgun.hepmc.Path = "hepmc" + +from Configurables import HepMCToEDMConverter +hepmc_converter = HepMCToEDMConverter() +hepmc_converter.hepmc.Path="hepmc" +hepmc_converter.GenParticles.Path="GenParticles" + +################## 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", "") +detectors_to_use=[ + '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 = WARNING + +# Geant4 service +# Configures the Geant simulation: geometry, physics list and user actions +from Configurables import SimG4Svc +geantservice = SimG4Svc("SimG4Svc", detector='SimG4DD4hepDetector', physicslist="SimG4FtfpBert", actions="SimG4FullSimActions") + +# 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 = "ECalBarrelEta" +ecalBarrelReadoutNamePhiEta = "ECalBarrelPhiEta" +# HCAL +hcalReadoutName = "HCalBarrelReadout" +extHcalReadoutName = "HCalExtBarrelReadout" + +# Configure saving of calorimeter hits +from Configurables import SimG4SaveCalHits +saveecalbarreltool = SimG4SaveCalHits("saveECalBarrelHits", readoutNames = [ecalBarrelReadoutName]) +saveecalbarreltool.CaloHits.Path = "ECalBarrelHits" + +savehcaltool = SimG4SaveCalHits("saveHCalBarrelHits",readoutNames = [hcalReadoutName]) +savehcaltool.CaloHits.Path = "HCalBarrelHits" + +# 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 = "GenParticles" + +from Configurables import SimG4Alg +geantsim = SimG4Alg("SimG4Alg", + outputs= ["SimG4SaveCalHits/saveECalBarrelHits", + "SimG4SaveCalHits/saveHCalBarrelHits", + ], + 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", + # sampling fraction obtained using SamplingFractionInLayers from DetStudies package + samplingFraction = [0.24833] * 1 + [0.09482] * 1 + [0.12242] * 1 + [0.14182] * 1 + [0.15667] * 1 + [0.16923] * 1 + [0.17980] * 1 + [0.20085] * 1, + readoutName = ecalBarrelReadoutName, + layerFieldName = "layer") + +from Configurables import CalibrateCaloHitsTool +calibHcells = CalibrateCaloHitsTool("CalibrateHCal", invSamplingFraction="41.66") + +# Create cells in ECal barrel +# 1. step - merge hits into cells with default Eta segmentation +# 2. step - rewrite the cellId using the Phi-Eta segmentation +from Configurables import CreateCaloCells +createEcalBarrelCellsStep1 = CreateCaloCells("CreateECalBarrelCellsStep1", + doCellCalibration=True, + calibTool = calibEcalBarrel, + addCellNoise=False, filterCellNoise=False, + OutputLevel=INFO, + hits="ECalBarrelHits", + cells="ECalBarrelCellsStep1") + +# Use Phi-Eta segmentation in ECal barrel +from Configurables import RedoSegmentation +resegmentEcalBarrel = RedoSegmentation("ReSegmentationEcal", + # old bitfield (readout) + oldReadoutName = ecalBarrelReadoutName, + # specify which fields are going to be altered (deleted/rewritten) + oldSegmentationIds = ["module"], + # new bitfield (readout), with new segmentation + newReadoutName = ecalBarrelReadoutNamePhiEta, + OutputLevel = INFO, + inhits = "ECalBarrelPositions", + outhits = "ECalBarrelCellsStep2") +createEcalBarrelCells = CreateCaloCells("CreateECalBarrelCells", + doCellCalibration=False, + addCellNoise=False, filterCellNoise=False, + OutputLevel=INFO, + hits="ECalBarrelCellsStep2", + cells="ECalBarrelCells") + +# 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") + +################ Output +from Configurables import PodioOutput +out = PodioOutput("out", + OutputLevel=INFO) +#Save information about generated particles & calorimeter cells, drop G4 hits and the intermediate steps +#out.outputCommands = ["drop *", "keep ECalBarrelCells", "keep HCalBarrelCells", "keep GenParticles","keep GenVertices"] +#Save all +out.outputCommands = ["keep *"] + +import uuid +out.filename = "output_fullCalo_SimAndDigi_"+str(momentum)+"GeV_"+uuid.uuid4().hex+".root" + +#CPU information +from Configurables import AuditorSvc, ChronoAuditor +chra = ChronoAuditor() +audsvc = AuditorSvc() +audsvc.Auditors = [chra] +genalg_pgun.AuditExecute = True +hepmc_converter.AuditExecute = True +geantsim.AuditExecute = True +createEcalBarrelCellsStep1.AuditExecute = True +resegmentEcalBarrel.AuditExecute = True +createEcalBarrelCells.AuditExecute = True +createHcalBarrelCells.AuditExecute = True +out.AuditExecute = True + +from Configurables import ApplicationMgr +ApplicationMgr( + TopAlg = [genalg_pgun, + hepmc_converter, + geantsim, + createEcalBarrelCellsStep1, + resegmentEcalBarrel, + createEcalBarrelCells, + createHcalBarrelCells, + out + ], + EvtSel = 'NONE', + EvtMax = 10, + ExtSvc = [geoservice, podioevent, geantservice, audsvc], + ) diff --git a/full-detector-simulations/FccCaloPerformance/runFullCaloSystem_ReconstructionSW_noiseFromFile.py b/full-detector-simulations/FccCaloPerformance/runFullCaloSystem_ReconstructionSW_noiseFromFile.py new file mode 100644 index 00000000..bad3e246 --- /dev/null +++ b/full-detector-simulations/FccCaloPerformance/runFullCaloSystem_ReconstructionSW_noiseFromFile.py @@ -0,0 +1,138 @@ +# Setup +# Names of cells collections +ecalBarrelCellsName = "ECalBarrelCells" +# Readouts +ecalBarrelReadoutName = "ECalBarrelPhiEta" + +# Number of events +num_events = -1 + +from Gaudi.Configuration import * +from Configurables import ApplicationMgr, FCCDataSvc, PodioOutput +import sys + +podioevent = FCCDataSvc("EventDataSvc") +import glob +podioevent.inputs=glob.glob("output_fullCalo_SimAndDigi_*.root") +# reads HepMC text file and write the HepMC::GenEvent to the data service +from Configurables import PodioInput +podioinput = PodioInput("PodioReader", + collections = [ecalBarrelCellsName, + "GenParticles", + ]) + + +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", "") +detectors_to_use=[ + '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 = WARNING + + +ecalBarrelNoisePath = "http://fccsw.web.cern.ch/fccsw/testsamples/elecNoise_ecalBarrelFCCee_50Ohm_traces1_4shieldWidth.root" +ecalBarrelNoiseHistName = "h_elecNoise_fcc_" + +# add noise, create all existing cells in detector +from Configurables import NoiseCaloCellsFromFileTool +noiseBarrel = NoiseCaloCellsFromFileTool("NoiseBarrel", + readoutName = ecalBarrelReadoutName, + noiseFileName = ecalBarrelNoisePath, + elecNoiseHistoName = ecalBarrelNoiseHistName, + activeFieldName = "layer", + addPileup = False, + numRadialLayers = 8) +from Configurables import TubeLayerPhiEtaCaloTool +barrelGeometry = TubeLayerPhiEtaCaloTool("EcalBarrelGeo", + readoutName = ecalBarrelReadoutName, + activeVolumeName = "LAr_sensitive", + activeFieldName = "layer", + fieldNames = ["system"], + fieldValues = [4], + activeVolumesNumber = 8) +from Configurables import CreateCaloCellsNoise +createEcalBarrelCells = CreateCaloCellsNoise("CreateECalBarrelCells", + geometryTool = barrelGeometry, + doCellCalibration=False, # already calibrated + addCellNoise=True, filterCellNoise=False, + noiseTool = noiseBarrel, + hits=ecalBarrelCellsName, + cells=ecalBarrelCellsName+"Noise", + OutputLevel=INFO) + +#Empty cells for parts of calorimeter not implemented yet +from Configurables import CreateEmptyCaloCellsCollection +createemptycells =CreateEmptyCaloCellsCollection("CreateEmptyCaloCells") +createemptycells.cells.Path = "emptyCaloCells" + +#Create calorimeter clusters +from GaudiKernel.PhysicalConstants import pi + +from Configurables import CaloTowerTool +towers = CaloTowerTool("towers", + deltaEtaTower = 0.01, deltaPhiTower = 2*pi/704., + ecalBarrelReadoutName = ecalBarrelReadoutName, + ecalEndcapReadoutName = "", + ecalFwdReadoutName = "", + hcalBarrelReadoutName = "", + hcalExtBarrelReadoutName = "", + hcalEndcapReadoutName = "", + hcalFwdReadoutName = "", + OutputLevel = INFO) +towers.ecalBarrelCells.Path = ecalBarrelCellsName + "Noise" +towers.ecalEndcapCells.Path = "emptyCaloCells" +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 +# approx in GeV: changed from default of 12 in FCC-hh +threshold = 2 + +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) +createClusters.clusters.Path = "CaloClustersFinal" +createClusters.clusterCells.Path = "CaloClusterCellsFinal" + +import uuid +out = PodioOutput("out", filename="output_allCalo_reco_noise_" + uuid.uuid4().hex + ".root") + +#CPU information +from Configurables import AuditorSvc, ChronoAuditor +chra = ChronoAuditor() +audsvc = AuditorSvc() +audsvc.Auditors = [chra] +podioinput.AuditExecute = True +createClusters.AuditExecute = True +out.AuditExecute = True + +ApplicationMgr( + TopAlg = [podioinput, + createemptycells, + createEcalBarrelCells, + createClusters, + out + ], + EvtSel = 'NONE', + EvtMax = num_events, + ExtSvc = [podioevent, geoservice])