diff --git a/benchmarks/lowq2_reconstruction/Snakefile b/benchmarks/lowq2_reconstruction/Snakefile index dbb5c94d..5bb3bf19 100644 --- a/benchmarks/lowq2_reconstruction/Snakefile +++ b/benchmarks/lowq2_reconstruction/Snakefile @@ -1,6 +1,9 @@ SIMOUTDIR="sim_output/lowq2_reconstruction/" ANALYSISDIR=SIMOUTDIR+"analysis/" +# Common list of plugins to ignore so eicrecon runs faster +LOWQ2_IGNORE_PLUGINS="janatop,LUMISPECCAL,ECTOF,BTOF,FOFFMTRK,RPOTS,B0TRK,MPGD,ECTRK,DRICH,DIRC,pid,tracking,EEMC,BEMC,FEMC,EHCAL,BHCAL,FHCAL,B0ECAL,ZDC,BTRK,BVTX,PFRICH,richgeo,evaluator,pid_lut,reco,rootfile" + ########################################################################################## ### Rules for checking the reconstruction of the electron momentum ########################################################################################## @@ -12,9 +15,9 @@ rule filter_hepmc_for_training: warmup="warmup/epic_ip6_extended.edm4hep.root", script=workflow.source_path("filterHEPMC3.py"), output: - SIMOUTDIR+"Low-Q2_Training_Events.hepmc3.tree.root", + SIMOUTDIR+"Low-Q2_Training_Events_{BEAMENERGY}.hepmc3.tree.root", params: - events="root://dtn-eic.jlab.org//volatile/eic/EPIC/EVGEN/SIDIS/pythia6-eic/1.0.0/18x275/q2_0to1/pythia_ep_noradcor_18x275_q2_0.000000001_1.0_run1.ab.hepmc3.tree.root", + events="root://dtn-eic.jlab.org//volatile/eic/EPIC/EVGEN/SIDIS/pythia6-eic/1.0.0/{BEAMENERGY}/q2_0to1/pythia_ep_noradcor_{BEAMENERGY}_q2_0.000000001_1.0_run1.ab.hepmc3.tree.root", shell: """ python {input.script} --outFile {output} --inFile {params.events} @@ -24,10 +27,11 @@ rule filter_hepmc_for_training: rule lowq2_training_sim: input: warmup="warmup/epic_ip6_extended.edm4hep.root", - events=SIMOUTDIR+"Low-Q2_Training_Events.hepmc3.tree.root", + events=SIMOUTDIR+"Low-Q2_Training_Events_{BEAMENERGY}.hepmc3.tree.root", + xml=lambda wildcards: workflow.source_path(f"epic_ip6_extended_{wildcards.BEAMENERGY}.xml"), geometry_lib=find_epic_libraries(), output: - SIMOUTDIR+"Low-Q2_Training_SimEvents_{CAMPAIGN}.edm4hep.root", + SIMOUTDIR+"Low-Q2_Training_SimEvents_{BEAMENERGY}_{CAMPAIGN}.edm4hep.root", params: DD4HEP_HASH=get_spack_package_hash("dd4hep"), NPSIM_HASH=get_spack_package_hash("npsim"), @@ -39,7 +43,7 @@ rule lowq2_training_sim: --random.seed 1 \ --inputFiles {input.events} \ --numberOfEvents 1000000 \ - --compactFile $DETECTOR_PATH/epic_ip6_extended.xml \ + --compactFile {input.xml} \ --printLevel WARNING \ --outputFile {output} \ --physics.rangecut 100*m @@ -48,20 +52,26 @@ rule lowq2_training_sim: # Run EICrecon to create TaggerTrackerReconstructedParticles rule lowq2_reconstruction_particles_eicrecon: params: - xml=os.getenv("DETECTOR_PATH")+"/epic_ip6_extended.xml", output_collections="TaggerTrackerReconstructedParticles,MCParticles,EventHeader", - ignore_plugins="janatop,LUMISPECCAL,ECTOF,BTOF,FOFFMTRK,RPOTS,B0TRK,MPGD,ECTRK,DRICH,DIRC,pid,tracking,EEMC,BEMC,FEMC,EHCAL,BHCAL,FHCAL,B0ECAL,ZDC,BTRK,BVTX,PFRICH,richgeo,evaluator,pid_lut,reco,rootfile", + ignore_plugins=LOWQ2_IGNORE_PLUGINS, + ebeam=lambda wildcards: wildcards.BEAMENERGY.split('x')[0], EICRECON_HASH=get_spack_package_hash("eicrecon"), input: - data=SIMOUTDIR+"Low-Q2_Training_SimEvents_{CAMPAIGN}.edm4hep.root", + data=SIMOUTDIR+"Low-Q2_Training_SimEvents_{BEAMENERGY}_{CAMPAIGN}.edm4hep.root", + xml=lambda wildcards: workflow.source_path(f"epic_ip6_extended_{wildcards.BEAMENERGY}.xml") output: - eicreconfile=ANALYSISDIR+"Low-Q2_test_Particles_{CAMPAIGN}.eicrecon.edm4hep.root", + eicreconfile=ANALYSISDIR+"Low-Q2_test_Particles_{BEAMENERGY}_{CAMPAIGN}.eicrecon.edm4hep.root", cache: True shell: """ - eicrecon -Ppodio:output_file={output.eicreconfile} -PLOWQ2:TaggerTrackerTransportationPreML:beamE=18.0 -PLOWQ2:TaggerTrackerTransportationPostML:beamE=18.0 -Pdd4hep:xml_files={params.xml} \ - -Pjana:nevents=500000 \ - -Ppodio:output_collections={params.output_collections} -Pplugins_to_ignore={params.ignore_plugins} {input.data} + eicrecon -Ppodio:output_file={output.eicreconfile} \ + -PLOWQ2:TaggerTrackerTransportationPreML:beamE={params.ebeam} \ + -PLOWQ2:TaggerTrackerTransportationPostML:beamE={params.ebeam} \ + -Pdd4hep:xml_files={input.xml} \ + -Pjana:nevents=500000 \ + -Ppodio:output_collections={params.output_collections} \ + -Pplugins_to_ignore={params.ignore_plugins} \ + {input.data} """ @@ -72,32 +82,39 @@ rule lowq2_reconstruction_particles_eicrecon: # Run EICrecon and create input tensors for reconstruction. rule lowq2_reconstruction_tensors_eicrecon: params: - xml=os.getenv("DETECTOR_PATH")+"/epic_ip6_extended.xml", output_collections="BackwardsBeamlineHits,TaggerTrackerFeatureTensor,TaggerTrackerTargetTensor,MCParticles,EventHeader", - ignore_plugins="janatop,LUMISPECCAL,ECTOF,BTOF,FOFFMTRK,RPOTS,B0TRK,MPGD,ECTRK,DRICH,DIRC,pid,tracking,EEMC,BEMC,FEMC,EHCAL,BHCAL,FHCAL,B0ECAL,ZDC,BTRK,BVTX,PFRICH,richgeo,evaluator,pid_lut,reco,rootfile", + ignore_plugins=LOWQ2_IGNORE_PLUGINS, + ebeam=lambda wildcards: wildcards.BEAMENERGY.split('x')[0], EICRECON_HASH=get_spack_package_hash("eicrecon"), input: - data=SIMOUTDIR+"Low-Q2_Training_SimEvents_{CAMPAIGN}.edm4hep.root", + data=SIMOUTDIR+"Low-Q2_Training_SimEvents_{BEAMENERGY}_{CAMPAIGN}.edm4hep.root", + xml=lambda wildcards: workflow.source_path(f"epic_ip6_extended_{wildcards.BEAMENERGY}.xml") output: - eicreconfile=ANALYSISDIR+"Low-Q2_Training_Tensors_{CAMPAIGN}.eicrecon.edm4hep.root", + eicreconfile=ANALYSISDIR+"Low-Q2_Training_Tensors_{BEAMENERGY}_{CAMPAIGN}.eicrecon.edm4hep.root", cache: True shell: """ - eicrecon -Ppodio:output_file={output.eicreconfile} -PLOWQ2:TaggerTrackerTransportationPreML:beamE=18.0 -Pdd4hep:xml_files={params.xml} \ - -Pjana:nevents=500000 -Pjana:nskip=500000 \ - -Ppodio:output_collections={params.output_collections} -Pplugins_to_ignore={params.ignore_plugins} {input.data} + eicrecon -Ppodio:output_file={output.eicreconfile} \ + -PLOWQ2:TaggerTrackerTransportationPreML:beamE={params.ebeam} \ + -Pdd4hep:xml_files={input.xml} \ + -Pjana:nevents=500000 -Pjana:nskip=500000 \ + -Ppodio:output_collections={params.output_collections} \ + -Pplugins_to_ignore={params.ignore_plugins} \ + {input.data} """ # Processes the simulation output data for training rule lowq2_steering_reconstruction_preparation: + params: + ebeam=lambda wildcards: wildcards.BEAMENERGY.split('x')[0] input: script=workflow.source_path("cleanData.C"), - data=ANALYSISDIR+"Low-Q2_Training_Tensors_{CAMPAIGN}.eicrecon.edm4hep.root", + data=ANALYSISDIR+"Low-Q2_Training_Tensors_{BEAMENERGY}_{CAMPAIGN}.eicrecon.edm4hep.root" output: - rootfile=ANALYSISDIR+"Low-Q2_Steering_Reconstruction_Data{CAMPAIGN}.root", + rootfile=ANALYSISDIR+"Low-Q2_Steering_Reconstruction_Data_{BEAMENERGY}_{CAMPAIGN}.root" shell: """ - root -l -b -q '{input.script}("{input.data}", "{output.rootfile}",18.0)' + root -l -b -q '{input.script}("{input.data}", "{output.rootfile}",{params.ebeam})' """ # Trains a regression model to predict the MCParticle momentum from the lowq2 tagger tracker hits. @@ -106,9 +123,9 @@ rule lowq2_steering_reconstruction_training: script=workflow.source_path("SteeringRegression.py"), scriptdata=workflow.source_path("LoadData.py"), scriptmodel=workflow.source_path("RegressionModel.py"), - data=ANALYSISDIR+"Low-Q2_Steering_Reconstruction_Data{CAMPAIGN}.root", + data=ANALYSISDIR+"Low-Q2_Steering_Reconstruction_Data_{BEAMENERGY}_{CAMPAIGN}.root", output: - onnxfile=ANALYSISDIR+"Low-Q2_Steering_Reconstruction_Test_{CAMPAIGN}.onnx", + onnxfile=ANALYSISDIR+"Low-Q2_Steering_Reconstruction_Test_{BEAMENERGY}_{CAMPAIGN}.onnx", shell: """ python {input.script} --dataFiles {input.data} --outModelFile {output.onnxfile} --epochs 1000 @@ -117,20 +134,28 @@ rule lowq2_steering_reconstruction_training: # Run eicrecon with the newly trained neural network rule lowq2_steering_reconstruction_particles_trained_eicrecon: params: - xml=os.getenv("DETECTOR_PATH")+"/epic_ip6_extended.xml", output_collections="TaggerTrackerReconstructedParticles,MCParticles,EventHeader", - ignore_plugins="janatop,LUMISPECCAL,ECTOF,BTOF,FOFFMTRK,RPOTS,B0TRK,MPGD,ECTRK,DRICH,DIRC,pid,tracking,EEMC,BEMC,FEMC,EHCAL,BHCAL,FHCAL,B0ECAL,ZDC,BTRK,BVTX,PFRICH,richgeo,evaluator,pid_lut,reco,rootfile", + ignore_plugins=LOWQ2_IGNORE_PLUGINS, + ebeam=lambda wildcards: wildcards.BEAMENERGY.split('x')[0], EICRECON_HASH=get_spack_package_hash("eicrecon"), input: - data=SIMOUTDIR+"Low-Q2_Training_SimEvents_{CAMPAIGN}.edm4hep.root", - onnxfile=ANALYSISDIR+"Low-Q2_Steering_Reconstruction_Test_{CAMPAIGN}.onnx", + data=SIMOUTDIR+"Low-Q2_Training_SimEvents_{BEAMENERGY}_{CAMPAIGN}.edm4hep.root", + onnxfile=ANALYSISDIR+"Low-Q2_Steering_Reconstruction_Test_{BEAMENERGY}_{CAMPAIGN}.onnx", + xml=lambda wildcards: workflow.source_path(f"epic_ip6_extended_{wildcards.BEAMENERGY}.xml") output: - eicreconfile=ANALYSISDIR+"Low-Q2_retrained_Particles_{CAMPAIGN}.eicrecon.edm4hep.root", + eicreconfile=ANALYSISDIR+"Low-Q2_retrained_Particles_{BEAMENERGY}_{CAMPAIGN}.eicrecon.edm4hep.root" cache: True shell: """ - eicrecon -Ppodio:output_file={output.eicreconfile} -Pjana:nevents=500000 -PLOWQ2:TaggerTrackerTransportationPreML:beamE=18.0 -PLOWQ2:TaggerTrackerTransportationPostML:beamE=18.0 -Pdd4hep:xml_files={params.xml} \ - -PLOWQ2:TaggerTrackerTransportationInference:modelPath={input.onnxfile} -Ppodio:output_collections={params.output_collections} -Pplugins_to_ignore={params.ignore_plugins} {input.data} + eicrecon -Ppodio:output_file={output.eicreconfile} \ + -Pjana:nevents=500000 \ + -PLOWQ2:TaggerTrackerTransportationPreML:beamE={params.ebeam} \ + -PLOWQ2:TaggerTrackerTransportationPostML:beamE={params.ebeam} \ + -Pdd4hep:xml_files={input.xml} \ + -PLOWQ2:TaggerTrackerTransportationInference:modelPath={input.onnxfile} \ + -Ppodio:output_collections={params.output_collections} \ + -Pplugins_to_ignore={params.ignore_plugins} \ + {input.data} """ ########################################################################################## @@ -139,27 +164,52 @@ rule lowq2_steering_reconstruction_particles_trained_eicrecon: # Produce hitsograms of th resolution of TaggerTrackerReconstructedParticles compared to their MCParticles rule lowq2_reconstruction_particles_test: + params: + ebeam=lambda wildcards: wildcards.BEAMENERGY.split('x')[0], + pbeam=lambda wildcards: wildcards.BEAMENERGY.split('x')[1] input: script=workflow.source_path("reconstructionAnalysis.C"), - data=ANALYSISDIR+"Low-Q2_{STAGE}_Particles_{CAMPAIGN}.eicrecon.edm4hep.root", + etpscript=workflow.source_path("format_Energy_Theta_Phi_Resolutions.C"), + eq2script=workflow.source_path("format_Energy_Q2_Acceptance.C"), + q2recoscript=workflow.source_path("format_Q2_Reconstruction.C"), + data=ANALYSISDIR+"Low-Q2_{STAGE}_Particles_{BEAMENERGY}_{CAMPAIGN}.eicrecon.edm4hep.root" + output: + rootfile=ANALYSISDIR+"Low-Q2_{STAGE}_Resolution_Results_{BEAMENERGY}_{CAMPAIGN}.root", + momentumCanvas=ANALYSISDIR+"{STAGE}_momentum_resolution_{BEAMENERGY}_{CAMPAIGN}.png", + energyThetaPhiCanvas=ANALYSISDIR+"{STAGE}_energy_theta_phi_resolution_{BEAMENERGY}_{CAMPAIGN}.png", + relationCanvas=ANALYSISDIR+"{STAGE}_relation_resolution_{BEAMENERGY}_{CAMPAIGN}.png", + resolutionGraphsCanvas=ANALYSISDIR+"{STAGE}_resolution_graphs_{BEAMENERGY}_{CAMPAIGN}.png", + acceptanceCanvas=ANALYSISDIR+"{STAGE}_acceptance_{BEAMENERGY}_{CAMPAIGN}.png", + energyQ2acceptanceCanvas=ANALYSISDIR+"{STAGE}_energy_Q2_acceptance_{BEAMENERGY}_{CAMPAIGN}.png", + xQ2acceptanceCanvas=ANALYSISDIR+"{STAGE}_x_Q2_acceptance_{BEAMENERGY}_{CAMPAIGN}.png", + Q2reconstructionCanvas=ANALYSISDIR+"{STAGE}_Q2_Reconstruction_{BEAMENERGY}_{CAMPAIGN}.png", + shell: + """ + root -l -b -q '{input.script}("{input.data}", {params.ebeam}, {params.pbeam}, "{output.rootfile}", "{output.momentumCanvas}", "{output.relationCanvas}", "{output.resolutionGraphsCanvas}", "{output.acceptanceCanvas}", "{output.energyQ2acceptanceCanvas}", "{output.xQ2acceptanceCanvas}")' + root -l -b -q '{input.etpscript}("{output.rootfile}", "{output.energyThetaPhiCanvas}")' + root -l -b -q '{input.eq2script}("{output.rootfile}", "{output.energyQ2acceptanceCanvas}")' + root -l -b -q '{input.q2recoscript}("{output.rootfile}", "{output.Q2reconstructionCanvas}")' + """ +rule lowq2_reconstruction_x_Q2_plots: + input: + script=workflow.source_path("format_X_Q2_Acceptance.C"), + data5=ANALYSISDIR+"Low-Q2_test_Resolution_Results_5x41_{CAMPAIGN}.root", + data10=ANALYSISDIR+"Low-Q2_test_Resolution_Results_10x100_{CAMPAIGN}.root", + data18=ANALYSISDIR+"Low-Q2_test_Resolution_Results_18x275_{CAMPAIGN}.root" output: - rootfile=ANALYSISDIR+"Low-Q2_{STAGE}_Resolution_Results_{CAMPAIGN}.root", - momentumCanvas=ANALYSISDIR+"{STAGE}_momentum_resolution_{CAMPAIGN}.png", - energyThetaPhiCanvas=ANALYSISDIR+"{STAGE}_energy_theta_phi_resolution_{CAMPAIGN}.png", - relationCanvas=ANALYSISDIR+"{STAGE}_relation_resolution_{CAMPAIGN}.png", - resolutionGraphsCanvas=ANALYSISDIR+"{STAGE}_resolution_graphs_{CAMPAIGN}.png", + xQ2Canvas=ANALYSISDIR+"test_X_Q2_Acceptance_{CAMPAIGN}.png" shell: """ - root -l -b -q '{input.script}("{input.data}", 18, "{output.rootfile}", "{output.momentumCanvas}", "{output.energyThetaPhiCanvas}", "{output.relationCanvas}", "{output.resolutionGraphsCanvas}")' + root -l -b -q '{input.script}("{input.data5}", "{input.data10}", "{input.data18}", "{output.xQ2Canvas}")' """ # Check the resloutions and offsets are within tollerance rule lowq2_reconstruction_particles_check: input: script=workflow.source_path("checkResolutions.C"), - rootfile=ANALYSISDIR+"Low-Q2_{STAGE}_Resolution_Results_{CAMPAIGN}.root" + rootfile=ANALYSISDIR+"Low-Q2_{STAGE}_Resolution_Results_{BEAMENERGY}_{CAMPAIGN}.root" output: - jsonfile=ANALYSISDIR+"Low-Q2_{STAGE}_Resolution_Results_{CAMPAIGN}.json" + jsonfile=ANALYSISDIR+"Low-Q2_{STAGE}_Resolution_Results_{BEAMENERGY}_{CAMPAIGN}.json" shell: """ root -l -b -q '{input.script}("{input.rootfile}", "{output.jsonfile}")' @@ -171,21 +221,24 @@ rule lowq2_reconstruction_particles_check: ########################################################################################## def get_onnx_input(wildcards): if wildcards.STAGE == "retrained": - return ANALYSISDIR + f"Low-Q2_Steering_Reconstruction_Test_{wildcards.CAMPAIGN}.onnx" + return ANALYSISDIR + f"Low-Q2_Steering_Reconstruction_Test_{wildcards.BEAMENERGY}_{wildcards.CAMPAIGN}.onnx" else: return [] rule lowq2_reconstruction: input: - rootfile=ANALYSISDIR+"Low-Q2_{STAGE}_Resolution_Results_{CAMPAIGN}.root", - momentumCanvas=ANALYSISDIR+"{STAGE}_momentum_resolution_{CAMPAIGN}.png", - energyThetaPhiCanvas=ANALYSISDIR+"{STAGE}_energy_theta_phi_resolution_{CAMPAIGN}.png", - relationCanvas=ANALYSISDIR+"{STAGE}_relation_resolution_{CAMPAIGN}.png", - resolutionGraphsCanvas=ANALYSISDIR+"{STAGE}_resolution_graphs_{CAMPAIGN}.png", - jsonfile=ANALYSISDIR+"Low-Q2_{STAGE}_Resolution_Results_{CAMPAIGN}.json", + rootfile=ANALYSISDIR+"Low-Q2_{STAGE}_Resolution_Results_{BEAMENERGY}_{CAMPAIGN}.root", + momentumCanvas=ANALYSISDIR+"{STAGE}_momentum_resolution_{BEAMENERGY}_{CAMPAIGN}.png", + energyThetaPhiCanvas=ANALYSISDIR+"{STAGE}_energy_theta_phi_resolution_{BEAMENERGY}_{CAMPAIGN}.png", + relationCanvas=ANALYSISDIR+"{STAGE}_relation_resolution_{BEAMENERGY}_{CAMPAIGN}.png", + resolutionGraphsCanvas=ANALYSISDIR+"{STAGE}_resolution_graphs_{BEAMENERGY}_{CAMPAIGN}.png", + acceptanceCanvas=ANALYSISDIR+"{STAGE}_acceptance_{BEAMENERGY}_{CAMPAIGN}.png", + energyQ2acceptanceCanvas=ANALYSISDIR+"{STAGE}_energy_Q2_acceptance_{BEAMENERGY}_{CAMPAIGN}.png", + xQ2acceptanceCanvas=ANALYSISDIR+"{STAGE}_x_Q2_acceptance_{BEAMENERGY}_{CAMPAIGN}.png", + jsonfile=ANALYSISDIR+"Low-Q2_{STAGE}_Resolution_Results_{BEAMENERGY}_{CAMPAIGN}.json", onnxfile=get_onnx_input, output: - directory("results/lowq2_reconstruction/{STAGE}_{CAMPAIGN}/") + directory("results/lowq2_reconstruction/{STAGE}_{BEAMENERGY}_{CAMPAIGN}/") shell: """ mkdir {output} @@ -197,8 +250,8 @@ rule lowq2_reconstruction: ########################################################################################## rule lowq2_reconstruction_local: input: - "results/lowq2_reconstruction/test_local/" + "results/lowq2_reconstruction/test_18x275_local/" rule lowq2_reconstruction_retrained: input: - "results/lowq2_reconstruction/retrained_local/" + "results/lowq2_reconstruction/retrained_18x275_local/" diff --git a/benchmarks/lowq2_reconstruction/epic_ip6_extended_10x100.xml b/benchmarks/lowq2_reconstruction/epic_ip6_extended_10x100.xml new file mode 100644 index 00000000..c2b264f4 --- /dev/null +++ b/benchmarks/lowq2_reconstruction/epic_ip6_extended_10x100.xml @@ -0,0 +1,138 @@ + + + + + + + + + + + + + + # EPIC Detector + - https://github.com/eic/epic + - https://github.com/eic/ip6 + + + + + EPIC + + + + + + + + ## Main Constant Definitions + + The ip6 (or other ip) defines should be included first. + These files have only a define tags. + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + ## World Volume + + The world is a simple box, but could be a union of multiple regions. + + + + + + + + + ## Detector Subsystems + + ### IP Subsystems + + The interaction point subsystems are included before the central detector subsystems. + This is becuase the IP subsystems, for example the beampipe, will define paramters + which are subsquently used in the central detector construction -- e.g. the vertex tracker + uses the beampipe OD to help define its placement. + + The IP subsystems include the Far forward and backward regions. The list of subsystem includes: + - Interaction region beampipe + - B0 tracker + - Off-momentum tracker + - Far forward roman pots + - Zero Degree Calorimeter + - Beam line magnets. + - and more... + + + + + ## Central tracking detectors + + + + + + ## Central beam pipe + + + + + ## Far foward detectors + + + + + ## Far backward detectors + + + diff --git a/benchmarks/lowq2_reconstruction/epic_ip6_extended_18x275.xml b/benchmarks/lowq2_reconstruction/epic_ip6_extended_18x275.xml new file mode 100644 index 00000000..b6a6d306 --- /dev/null +++ b/benchmarks/lowq2_reconstruction/epic_ip6_extended_18x275.xml @@ -0,0 +1,138 @@ + + + + + + + + + + + + + + # EPIC Detector + - https://github.com/eic/epic + - https://github.com/eic/ip6 + + + + + EPIC + + + + + + + + ## Main Constant Definitions + + The ip6 (or other ip) defines should be included first. + These files have only a define tags. + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + ## World Volume + + The world is a simple box, but could be a union of multiple regions. + + + + + + + + + ## Detector Subsystems + + ### IP Subsystems + + The interaction point subsystems are included before the central detector subsystems. + This is becuase the IP subsystems, for example the beampipe, will define paramters + which are subsquently used in the central detector construction -- e.g. the vertex tracker + uses the beampipe OD to help define its placement. + + The IP subsystems include the Far forward and backward regions. The list of subsystem includes: + - Interaction region beampipe + - B0 tracker + - Off-momentum tracker + - Far forward roman pots + - Zero Degree Calorimeter + - Beam line magnets. + - and more... + + + + + ## Central tracking detectors + + + + + + ## Central beam pipe + + + + + ## Far foward detectors + + + + + ## Far backward detectors + + + diff --git a/benchmarks/lowq2_reconstruction/epic_ip6_extended_5x41.xml b/benchmarks/lowq2_reconstruction/epic_ip6_extended_5x41.xml new file mode 100644 index 00000000..b2184b95 --- /dev/null +++ b/benchmarks/lowq2_reconstruction/epic_ip6_extended_5x41.xml @@ -0,0 +1,138 @@ + + + + + + + + + + + + + + # EPIC Detector + - https://github.com/eic/epic + - https://github.com/eic/ip6 + + + + + EPIC + + + + + + + + ## Main Constant Definitions + + The ip6 (or other ip) defines should be included first. + These files have only a define tags. + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + ## World Volume + + The world is a simple box, but could be a union of multiple regions. + + + + + + + + + ## Detector Subsystems + + ### IP Subsystems + + The interaction point subsystems are included before the central detector subsystems. + This is becuase the IP subsystems, for example the beampipe, will define paramters + which are subsquently used in the central detector construction -- e.g. the vertex tracker + uses the beampipe OD to help define its placement. + + The IP subsystems include the Far forward and backward regions. The list of subsystem includes: + - Interaction region beampipe + - B0 tracker + - Off-momentum tracker + - Far forward roman pots + - Zero Degree Calorimeter + - Beam line magnets. + - and more... + + + + + ## Central tracking detectors + + + + + + ## Central beam pipe + + + + + ## Far foward detectors + + + + + ## Far backward detectors + + + diff --git a/benchmarks/lowq2_reconstruction/format_Energy_Q2_Acceptance.C b/benchmarks/lowq2_reconstruction/format_Energy_Q2_Acceptance.C new file mode 100644 index 00000000..e4906aa4 --- /dev/null +++ b/benchmarks/lowq2_reconstruction/format_Energy_Q2_Acceptance.C @@ -0,0 +1,91 @@ +// Format the log10(Q2) vs E acceptance histogram from reconstructionAnalysis output +// Usage: root -l 'format_Energy_Q2_Acceptance.C("reconstruction_results.root", "Energy_Q2_Acceptance.png")' + +#include +#include +#include +#include +#include +#include + +void format_Energy_Q2_Acceptance(const char* infile = "reconstruction_results.root", + const char* outname = "Energy_Q2_Acceptance.png") { + // Style + gStyle->SetOptStat(0); + gStyle->SetPalette(kViridis); + gStyle->SetNumberContours(60); + + // Load histogram + TFile* f = TFile::Open(infile, "READ"); + if (!f || f->IsZombie()) { + printf("Cannot open %s\n", infile); + return; + } + + TH2* h = (TH2*)f->Get("hLog10Q2_vs_E_acceptance"); + if (!h) { + printf("Histogram hLog10Q2_vs_E_acceptance not found\n"); + f->Close(); + return; + } + + // Canvas + TCanvas* c = new TCanvas("c", "Energy vs Q2 Acceptance", 900, 700); + c->SetLeftMargin(0.12); + c->SetRightMargin(0.16); + c->SetBottomMargin(0.12); + c->SetTopMargin(0.08); + + // Draw + h->SetTitle(";E_{e^{-}} [GeV];log_{10}(Q^{2}) [GeV^{2}]"); + h->SetMinimum(0.0); + h->SetMaximum(1.0); + h->Draw("COLZ"); + + double lm = gPad->GetLeftMargin(); + double rm = gPad->GetRightMargin(); + double tm = gPad->GetTopMargin(); + + // Title + TLatex title; + title.SetNDC(); + title.SetTextFont(72); + // title.SetTextAlign(33); // right-top + title.SetTextSize(0.035); + title.SetTextAlign(12); + title.DrawLatex(lm, 1.0 - tm + 0.3*tm, "ePIC Simulation 25.10.4"); + + TLatex title2; + title2.SetNDC(); + title2.SetTextFont(72); + title2.SetTextAlign(32); // right-top + title2.SetTextSize(0.035); + double x2 = 1.0 - rm; + double y2 = 1.0 - tm + 0.3*tm; // adjust 0.5*tm to move up/down + title2.DrawLatex(x2, y2, "Low Q^{2} Tagger: Q^{2} vs E Acceptance"); + + c->Update(); + + // Adjust palette (color bar) to match main histogram vertical extent + TPaletteAxis* pal = (TPaletteAxis*)h->GetListOfFunctions()->FindObject("palette"); + if (pal) { + double y1 = c->GetBottomMargin(); // lower bound at bottom margin + double y2 = 1.0 - c->GetTopMargin(); // upper bound at top margin + pal->SetY1NDC(y1); + pal->SetY2NDC(y2); + + // Set a consistent width for the palette axis + double x2 = 1.0 - c->GetRightMargin() + 0.05; // slight outward shift + double x1 = x2 - 0.046; // width of color bar + pal->SetX1NDC(x1); + pal->SetX2NDC(x2); + + pal->SetLabelSize(0.035); + pal->SetTitleSize(0.035); + c->Modified(); + c->Update(); + } + c->SaveAs(outname); + + f->Close(); +} diff --git a/benchmarks/lowq2_reconstruction/format_Energy_Theta_Phi_Resolutions.C b/benchmarks/lowq2_reconstruction/format_Energy_Theta_Phi_Resolutions.C new file mode 100644 index 00000000..c5d0e9d0 --- /dev/null +++ b/benchmarks/lowq2_reconstruction/format_Energy_Theta_Phi_Resolutions.C @@ -0,0 +1,157 @@ +// Macro to build a nicely formatted canvas for energy, scattering angle and filtered phi resolutions +// Reads histograms previously produced by reconstructionAnalysis and saved in a ROOT file. +// Histograms expected: +// 2D: E_vs_E, theta_vs_theta, phi_filtered_vs_phi +// 1D: E_res, theta_diff, phi_filtered_diff +// Usage (ROOT): .x EnergyThetaPhiResolutions.C("reconstruction_results.root", "energy_theta_phi_filtered_nice.png") +// Or compile: root -l -b -q 'EnergyThetaPhiResolutions.C("reconstruction_results.root")' + +#include +#include +#include +#include + +#include "TFile.h" +#include "TH1.h" +#include "TH2.h" +#include "TCanvas.h" +#include "TStyle.h" +#include "TColor.h" +#include "TLatex.h" +#include "TLegend.h" +#include "TSystem.h" + +// Helper to annotate mean and RMS for 1D histograms +void AnnotateStats(TH1* h, const char* units="") { + if(!h) return; + double mean = h->GetMean(); + double rms = h->GetRMS(); + TLatex *lat = new TLatex(); + lat->SetNDC(); + lat->SetTextSize(0.04); + lat->SetTextColor(kGray+2); + lat->DrawLatex(0.15, 0.85, Form("Mean = %.3g %s", mean, units)); + lat->DrawLatex(0.15, 0.78, Form("RMS = %.3g %s", rms, units)); +} + +// Apply common style to 2D histograms +void Style2D(TH2* h, const char* title) { + if(!h) return; + h->SetTitle(title); + h->GetXaxis()->SetTitleOffset(1.1); + h->GetYaxis()->SetTitleOffset(1.1); + h->GetXaxis()->SetTitleSize(0.045); + h->GetYaxis()->SetTitleSize(0.045); + h->GetXaxis()->SetLabelSize(0.04); + h->GetYaxis()->SetLabelSize(0.04); +} + +// Apply common style to 1D histograms +void Style1D(TH1* h, const char* title, Color_t col) { + if(!h) return; + h->SetTitle(title); + h->SetLineColor(col); + h->SetMarkerColor(col); + h->SetMarkerStyle(20); + h->SetLineWidth(2); + h->GetXaxis()->SetTitleOffset(1.1); + h->GetYaxis()->SetTitleOffset(1.1); + h->GetXaxis()->SetTitleSize(0.045); + h->GetYaxis()->SetTitleSize(0.045); + h->GetXaxis()->SetLabelSize(0.04); + h->GetYaxis()->SetLabelSize(0.04); +} + +void format_Energy_Theta_Phi_Resolutions(const char* inFile="reconstruction_results.root", + const char* outPng="energy_theta_phi_filtered_nice.png") { + // Global style tweaks + gStyle->SetOptStat(0); + gStyle->SetNumberContours(60); + gStyle->SetPalette(kViridis); + + std::unique_ptr f(TFile::Open(inFile, "READ")); + if(!f || f->IsZombie()) { + std::cerr << "ERROR: Cannot open input file: " << inFile << std::endl; + return; + } + + // Fetch histograms + TH2* hE_vs_E = dynamic_cast(f->Get("E_vs_E")); + TH2* hTheta_vs_theta = dynamic_cast(f->Get("theta_vs_theta")); + TH2* hPhi_filtered_vs_phi = dynamic_cast(f->Get("phi_filtered_vs_phi")); + TH1* hE_res = dynamic_cast(f->Get("E_res")); + TH1* hTheta_diff = dynamic_cast(f->Get("theta_diff")); + TH1* hPhi_filtered_diff = dynamic_cast(f->Get("phi_filtered_diff")); + + // Basic existence check + if(!hE_vs_E || !hTheta_vs_theta || !hPhi_filtered_vs_phi || !hE_res || !hTheta_diff || !hPhi_filtered_diff) { + std::cerr << "ERROR: Missing one or more required histograms in file: " << inFile << std::endl; + return; + } + + // Scale energy resolution to percentage (fraction -> %) + if(hE_res) hE_res->Scale(100.0); + + // Style histograms + Style2D(hE_vs_E, "Reco vs MC Energy; E_{e} MC [GeV]; E_{e} reco [GeV]"); + Style2D(hTheta_vs_theta, "Reco vs MC Scattering Angle; #theta_{MC} [mrad]; #theta_{reco} [mrad]"); + Style2D(hPhi_filtered_vs_phi, "Reco vs MC #phi ( #theta > 1 mrad ); #phi_{MC} [deg]; #phi_{reco} [deg]"); + + Style1D(hE_res, "Energy Resolution; #DeltaE/E_{MC} [%]", kAzure+2); + Style1D(hTheta_diff, "Scattering Angle Difference; #Delta#theta [mrad]", kOrange+7); + Style1D(hPhi_filtered_diff, "#phi Difference ( #theta > 1 mrad ); #Delta#phi [deg]", kGreen+2); + // Create canvas with extra top space (dedicated title pad) without altering 3x2 plot layout + TCanvas *c = new TCanvas("energy_theta_phi_filtered_nice", "Energy/ScattAngle/Phi (Filtered) Resolutions", 2600, 1600); + TPad *padTitle = new TPad("padTitle","Title",0.0,0.95,1.0,1.0); // ~7% height for title + padTitle->SetFillColor(0); padTitle->SetBorderSize(0); padTitle->SetMargin(0,0,0,0); padTitle->Draw(); + TPad *padMain = new TPad("padMain","Plots",0.0,0.0,1.0,0.95); + padMain->SetFillColor(0); padMain->SetBorderSize(0); padMain->Draw(); + padMain->cd(); + padMain->Divide(3,2,0.005,0.005); + + // Adjust margins for each of the 6 plot pads (keep plotting area proportions) + for(int pad=1; pad<=6; ++pad) { + padMain->cd(pad); + gPad->SetTopMargin(0.04); + gPad->SetLeftMargin(0.12); + gPad->SetBottomMargin(0.13); + gPad->SetRightMargin(0.05); + } + + // 2D plots with log z (suppress main titles, keep axis titles) + padMain->cd(1); hE_vs_E->SetTitle(""); hE_vs_E->Draw("COL"); gPad->SetLogz(); + padMain->cd(2); hTheta_vs_theta->SetTitle(""); hTheta_vs_theta->Draw("COL"); gPad->SetLogz(); + padMain->cd(3); hPhi_filtered_vs_phi->SetTitle(""); hPhi_filtered_vs_phi->Draw("COL"); gPad->SetLogz(); + + // 1D resolution/difference plots (suppress titles) + padMain->cd(4); hE_res->SetTitle(""); hE_res->Draw("HIST"); hE_res->GetYaxis()->SetMaxDigits(3); AnnotateStats(hE_res, "%"); + padMain->cd(5); hTheta_diff->SetTitle(""); hTheta_diff->Draw("HIST"); hTheta_diff->GetYaxis()->SetMaxDigits(3); AnnotateStats(hTheta_diff, "mrad"); + padMain->cd(6); hPhi_filtered_diff->SetTitle(""); hPhi_filtered_diff->Draw("HIST"); hPhi_filtered_diff->GetYaxis()->SetMaxDigits(3); AnnotateStats(hPhi_filtered_diff, "deg"); + + // Overall title in dedicated title pad (top-left) + padTitle->cd(); + + TLatex overall; + overall.SetNDC(); + overall.SetTextFont(72); + overall.SetTextSize(0.6); + overall.SetTextAlign(13); + overall.DrawLatex(0.02, 0.64, "ePIC Simulation 25.10.4"); + + TLatex details; + details.SetNDC(); + details.SetTextFont(72); + details.SetTextSize(0.6); + details.SetTextAlign(33); // right-top alignment + details.DrawLatex(0.98, 0.64, "Low-Q^{2} Tagger Resolutions"); + + // Improve aesthetics + // for(int pad=1; pad<=6; ++pad) { + // c->cd(pad); + // gPad->SetGrid(); + // } + + c->Update(); + c->SaveAs(outPng); + +} diff --git a/benchmarks/lowq2_reconstruction/format_Q2_Reconstruction.C b/benchmarks/lowq2_reconstruction/format_Q2_Reconstruction.C new file mode 100644 index 00000000..f274d46c --- /dev/null +++ b/benchmarks/lowq2_reconstruction/format_Q2_Reconstruction.C @@ -0,0 +1,91 @@ +// Format the log10(Q^2) reconstruction vs MC histogram from reconstructionAnalysis output +// Usage: root -l 'format_Q2_Reconstruction.C("reconstruction_results.root", "log10Q2_Reconstruction.png")' + +#include +#include +#include +#include +#include +#include +#include + +void format_Q2_Reconstruction(const char* infile = "reconstruction_results.root", + const char* outname = "log10Q2_Reconstruction.png") { + // Style + gStyle->SetOptStat(0); + gStyle->SetPalette(kViridis); + gStyle->SetNumberContours(60); + + // Load histogram + TFile* f = TFile::Open(infile, "READ"); + if (!f || f->IsZombie()) { + printf("Cannot open %s\n", infile); + return; + } + + TH2* h = (TH2*)f->Get("log10Q2_vs_log10Q2"); + if (!h) { + printf("Histogram log10Q2_vs_log10Q2 not found\n"); + f->Close(); + return; + } + + // Canvas + TCanvas* c = new TCanvas("c", "Q2 Reconstruction", 1200, 700); + c->SetLeftMargin(0.10); + c->SetRightMargin(0.10); + c->SetBottomMargin(0.12); + c->SetTopMargin(0.08); + c->SetLogz(); + + // Draw + h->SetTitle(";log_{10}(Q^{2})_{MC} [GeV^{2}];log_{10}(Q^{2})_{reco} [GeV^{2}]"); + h->Draw("COLZ"); + + + double lm = gPad->GetLeftMargin(); + double rm = gPad->GetRightMargin(); + double tm = gPad->GetTopMargin(); + + // Title + TLatex title; + title.SetNDC(); + title.SetTextFont(72); + title.SetTextSize(0.035); + title.SetTextAlign(12); + title.DrawLatex(lm, 1.0 - tm + 0.3*tm, "ePIC Simulation 25.10.4"); + + TLatex title2; + title2.SetNDC(); + title2.SetTextFont(72); + title2.SetTextAlign(32); // right-top + title2.SetTextSize(0.035); + double x2 = 1.0 - rm; + double y2 = 1.0 - tm + 0.3*tm; // adjust 0.3*tm to move up/down + title2.DrawLatex(x2, y2, "Low Q^{2} Tagger: Q^{2} Reconstruction"); + + c->Update(); + + // Adjust palette (color bar) to match main histogram vertical extent + TPaletteAxis* pal = (TPaletteAxis*)h->GetListOfFunctions()->FindObject("palette"); + if (pal) { + double y1 = c->GetBottomMargin(); // lower bound at bottom margin + double y2 = 1.0 - c->GetTopMargin(); // upper bound at top margin + pal->SetY1NDC(y1); + pal->SetY2NDC(y2); + + // Set a consistent width for the palette axis + double x2_pal = 1.0 - c->GetRightMargin() + 0.05; // slight outward shift + double x1_pal = x2_pal - 0.046; // width of color bar + pal->SetX1NDC(x1_pal); + pal->SetX2NDC(x2_pal); + + pal->SetLabelSize(0.035); + pal->SetTitleSize(0.035); + c->Modified(); + c->Update(); + } + c->SaveAs(outname); + + f->Close(); +} diff --git a/benchmarks/lowq2_reconstruction/format_X_Q2_Acceptance.C b/benchmarks/lowq2_reconstruction/format_X_Q2_Acceptance.C new file mode 100644 index 00000000..5754f9cb --- /dev/null +++ b/benchmarks/lowq2_reconstruction/format_X_Q2_Acceptance.C @@ -0,0 +1,204 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +void format_X_Q2_Acceptance(TString file5 = "lowq2-5x41/reconstruction_results.root", + TString file10 = "lowq2-10x100/reconstruction_results.root", + TString file18 = "lowq2-18x275/reconstruction_results.root", + TString outname = "X_Q2_Acceptance.png", + bool show_unresolvable_region = true) { + + // Global style tweaks + gStyle->SetOptStat(0); + gStyle->SetPadGridX(kFALSE); // Turn off X grid + gStyle->SetPadGridY(kFALSE); // Turn off Y grid + + bool add_epic_logo = false; + + // Define input files and their labels + std::vector files = { + file5, + file10, + file18 + }; + + std::vector labels = { + "5x41 GeV", + "10x100 GeV", + "18x275 GeV" + }; + + std::vector colors = {kP6Blue, kP6Yellow, kP6Red}; + std::vector fillStyles = {3004, 3005, 3006}; // Different hatched patterns for each energy + + TString histName = "hLog10Q2_vs_log10x_acceptance"; + TString meanHistName = "Q2_rel_res_vs_log10Q2_1"; // FitSlicesY creates this name from the relative resolution histogram + TString resolutionHistName = "Q2_rel_res_vs_log10Q2_2"; // FitSlicesY creates this name from the relative resolution histogram + + // Create canvas + TCanvas *c1 = new TCanvas("c1", "x-Q^{2} Acceptance with Resolution Regions", 1000, 700); + + TLegend *leg = new TLegend(0.16, 0.57, 0.47, 0.88); + leg->SetBorderSize(0); + leg->SetFillStyle(0); + + bool firstDraw = true; + double resolution_threshold = 1.0; + + printf("Looking for histogram: %s\\n", histName.Data()); + printf("Looking for resolution histogram: %s\\n", resolutionHistName.Data()); + + // Loop over files to draw acceptance regions + for (size_t i = 0; i < files.size(); i++) { + TFile *f = TFile::Open(files[i]); + if (!f || f->IsZombie()) { + printf("Cannot open file: %s\n", files[i].Data()); + continue; + } + + // Determine Q2 threshold for this specific file + double q2_threshold = -10.0; // Default fallback + TH1 *h_mean = (TH1*)f->Get(meanHistName); + TH1 *h_res = (TH1*)f->Get(resolutionHistName); + if (h_res && h_mean) { + for (int bin = 1; bin <= h_res->GetNbinsX(); bin++) { + double log10q2 = h_res->GetBinCenter(bin); + double rel_resolution = h_res->GetBinContent(bin); // This is now relative resolution directly + double mean_offset = h_mean->GetBinContent(bin); + + double error = h_res->GetBinError(bin); + + // Skip bins with insufficient data or invalid fits + if (error == 0) continue; // Require minimum entries + if (rel_resolution <= 0) continue; // Skip non-positive resolution + if (rel_resolution > 2.0) continue; // Skip unreasonably large relative resolution (>200%) + + // std::cout << "Bin " << bin << ": log10q2 = " << log10q2 << ", rel_resolution = " << rel_resolution << std::endl; + + // Look for where relative resolution exceeds threshold (poor resolution region) + if (rel_resolution < resolution_threshold && std::abs(mean_offset) < resolution_threshold) { + q2_threshold = log10q2; + break; + } + } + } + + printf("%s: Q2 threshold = log10(Q2) = %.2f (Q2 = %.2e GeV^2)\n", + labels[i].Data(), q2_threshold, std::pow(10.0, q2_threshold)); + + TH2 *h = (TH2*)f->Get(histName); + if (!h) { + printf("Cannot find histogram %s in %s\n", histName.Data(), files[i].Data()); + f->Close(); + continue; + } + + printf("Found histogram %s with %d entries, max value: %f\n", + histName.Data(), (int)h->GetEntries(), h->GetMaximum()); + + h->SetContour(1); + h->SetContourLevel(0, 0.005); + h->SetLineColor(colors[i]); + h->SetLineWidth(2); + + if (firstDraw) { + h->Draw("CONT3"); + h->SetTitle(";log_{10}(x);log_{10}(Q^{2}) [GeV^{2}]"); + firstDraw = false; + } else { + h->Draw("CONT3 SAME"); + } + // Add to legend + leg->AddEntry(h, labels[i], "l"); + // break; + + if (show_unresolvable_region){ + TH2 *cut_hist = (TH2*)h->Clone("cut_hist"); + //Set range to below the Q2 cut + cut_hist->GetYaxis()->SetRangeUser(-10.0, q2_threshold); + + // Set fill properties with unique pattern for each energy + cut_hist->SetFillColorAlpha(colors[i],0.3); + cut_hist->SetFillStyle(fillStyles[i]); // Use different pattern for each energy + cut_hist->SetLineWidth(0); // No outline + + // Draw as filled contour + cut_hist->SetContour(1); + cut_hist->SetContourLevel(0, 0.005); + cut_hist->Draw("CONT0 F SAME"); + + leg->AddEntry(cut_hist, Form("%s Smallest resolvable bin", labels[i].Data()), "f"); + } + + + + // f->Close(); + } + + + // Ensure we're back on the main canvas before drawing titles + // c1->cd(); + + double lm = gPad->GetLeftMargin(); + double rm = gPad->GetRightMargin(); + double tm = gPad->GetTopMargin(); + + TLatex title; + title.SetNDC(); + title.SetTextFont(72); + title.SetTextSize(0.035); + title.SetTextAlign(12); + title.DrawLatex(lm, 1.0 - tm + 0.3*tm, "ePIC Simulation 25.10.4"); + + TLatex title2; + title2.SetNDC(); + title2.SetTextFont(72); + title2.SetTextAlign(32); + title2.SetTextSize(0.035); + double x2 = 1.0 - rm; + double y2 = 1.0 - tm + 0.3*tm; + title2.DrawLatex(x2, y2, "Low Q^{2} Tagger: x-Q^{2} Acceptance"); + + // Add description if showing unresolvable region + // if (show_unresolvable_region) { + // TLatex desc; + // desc.SetNDC(); + // desc.SetTextFont(42); + // desc.SetTextSize(0.025); + // desc.SetTextAlign(12); + // desc.DrawLatex(0.5, 0.85, "Shaded area: Poor Q^{2} Resolution Region"); + // } + + // Ensure we're on the main canvas + // c1->cd(); + + // Draw the legend + leg->Draw(); + + if(add_epic_logo) + { + // ===== Add ePIC logo to the figure if desired ====== + TImage *logo = TImage::Open("EPIC-logo_black_small.png"); + TPad *pad2 = new TPad("pad2", "Pad 2", 0.8, 0.8, 0.93, 0.93); // Create a new pad and then draw the image in it + pad2->Draw(); + pad2->cd(); // Enter the new pad + logo->Draw(); + } + + + c1->Update(); + + c1->SaveAs(outname); +} \ No newline at end of file diff --git a/benchmarks/lowq2_reconstruction/reconstructionAnalysis.C b/benchmarks/lowq2_reconstruction/reconstructionAnalysis.C index 8a0c5aa0..b1b06270 100644 --- a/benchmarks/lowq2_reconstruction/reconstructionAnalysis.C +++ b/benchmarks/lowq2_reconstruction/reconstructionAnalysis.C @@ -8,15 +8,19 @@ #include "edm4eic/ReconstructedParticleCollection.h" #include "TCanvas.h" #include "TStyle.h" - -void reconstructionAnalysis(TString inFile = "/home/simong/EIC/detector_benchmarks_anl/sim_output/lowq2_reconstruction/analysis/Low-Q2_retrained_Particles_new.eicrecon.edm4hep.root", - float beamEnergy = 18.0, - TString outFile = "reconstruction_results.root", - TString momentumCanvasName = "momentum_resolution.png", - TString energyThetaPhiCanvasName = "energy_theta_phi_resolution.png", - TString relationCanvasName = "relation_resolution.png", - TString resolutionGraphsCanvasName = "resolution_graphs.png", - std::string particleCollectionName = "TaggerTrackerReconstructedParticles") { +#include "TMath.h" + +void reconstructionAnalysis( TString inFile = "/home/simong/EIC/detector_benchmarks_anl/sim_output/lowq2_reconstruction/analysis/Low-Q2_retrained_Particles_new.eicrecon.edm4hep.root", + float electronEnergy = 18.0, + float protonEnergy = 275.0, + TString outFile = "reconstruction_results.root", + TString momentumCanvasName = "momentum_resolution.png", + TString relationCanvasName = "relation_resolution.png", + TString resolutionGraphsCanvasName = "resolution_graphs.png", + TString acceptanceCanvasName = "acceptance_canvas.png", + TString xQ2acceptanceCanvasName = "x_Q2_acceptance_canvas.png", + TString Q2xResolutionCanvasName = "Q2_x_resolution.png", + std::string particleCollectionName = "TaggerTrackerReconstructedParticles") { //Set ROOT style gStyle->SetPadLeftMargin(0.1); // Set left margin @@ -33,36 +37,156 @@ void reconstructionAnalysis(TString inFile = "/home/simong/EIC/detec ROOT::RDataFrame d0("events",inFile, {"MCParticles",particleCollectionName}); - auto filterDF = d0.Define("SimParticles", "MCParticles[MCParticles.generatorStatus==1 && MCParticles.PDG==11]") - .Filter("SimParticles.size()==1") - .Filter(particleCollectionName+".size()==1"); + // Base MC scattered electron plus real beam particle definitions + auto baseDF = d0 + // Scattered (final state) electron (generatorStatus==1, PDG==11) + .Define("SimParticles", "MCParticles[MCParticles.generatorStatus==1 && MCParticles.PDG==11]") + .Define("hasElectron", "SimParticles.size()==1") + // Beam particles (incoming) with generatorStatus==4 + .Define("BeamElectrons", "MCParticles[MCParticles.generatorStatus==4 && MCParticles.PDG==11]") + .Define("BeamProtons", "MCParticles[MCParticles.generatorStatus==4 && MCParticles.PDG==2212]") + .Define("hasBeamElectron", "BeamElectrons.size()==1") + .Define("hasBeamProton", "BeamProtons.size()==1") + // Use scattered electron momentum for mc truth k' + .Define("mc_momentum", "hasElectron ? SimParticles[0].momentum : MCParticles[0].momentum") + .Define("px_mc", "mc_momentum.x") + .Define("py_mc", "mc_momentum.y") + .Define("pz_mc", "mc_momentum.z") + .Define("E_mc", "sqrt(px_mc*px_mc + py_mc*py_mc + pz_mc*pz_mc + (hasElectron ? SimParticles[0].mass*SimParticles[0].mass : MCParticles[0].mass*MCParticles[0].mass))") + // Beam electron 4-vector (fallback to ideal if missing) + .Define("beam_e_px", [electronEnergy](const ROOT::RVec& be){ + if(be.size()==1) return be[0].momentum.x; + // ideal head-on along -z + return 0.0; }, {"BeamElectrons"}) + .Define("beam_e_py", [electronEnergy](const ROOT::RVec& be){ + if(be.size()==1) return be[0].momentum.y; return 0.0; }, {"BeamElectrons"}) + .Define("beam_e_pz", [electronEnergy](const ROOT::RVec& be){ + if(be.size()==1) return be[0].momentum.z; // fallback approximate -|p| + const double me = 0.00051099895; const double pe = std::sqrt(std::max(0.0, electronEnergy*electronEnergy - me*me)); return -pe; }, {"BeamElectrons"}) + .Define("beam_e_E", [electronEnergy](const ROOT::RVec& be){ + if(be.size()==1){ + auto p = be[0].momentum; + double m = be[0].mass; + return std::sqrt(p.x*p.x + p.y*p.y + p.z*p.z + m*m); + } + return (double)electronEnergy; }, {"BeamElectrons"}) + // Beam proton 4-vector (fallback to ideal if missing) + .Define("beam_p_px", [] (const ROOT::RVec& bp){ if(bp.size()==1) return bp[0].momentum.x; return 0.0; }, {"BeamProtons"}) + .Define("beam_p_py", [] (const ROOT::RVec& bp){ if(bp.size()==1) return bp[0].momentum.y; return 0.0; }, {"BeamProtons"}) + .Define("beam_p_pz", [protonEnergy](const ROOT::RVec& bp){ + if(bp.size()==1) return bp[0].momentum.z; // fallback +|p| + const double Mp = 0.9382720813; const double pp = std::sqrt(std::max(0.0, protonEnergy*protonEnergy - Mp*Mp)); return pp; }, {"BeamProtons"}) + .Define("beam_p_E", [protonEnergy](const ROOT::RVec& bp){ + if(bp.size()==1){ + auto p = bp[0].momentum; + double m = bp[0].mass; + return std::sqrt(p.x*p.x + p.y*p.y + p.z*p.z + m*m); + } + return (double)protonEnergy; }, {"BeamProtons"}) + // Raw polar angle and derived scattering angle (pi - raw) + .Define("theta_mc_raw", "std::atan2(std::sqrt(px_mc*px_mc + py_mc*py_mc), pz_mc)") + .Define("theta_mc", "TMath::Pi() - theta_mc_raw") + .Define("theta_mc_mrad", "1000.0*theta_mc") + .Define("phi_mc_rad", "std::atan2(py_mc, px_mc)") + .Define("phi_mc", "TMath::RadToDeg()*phi_mc_rad") + // Invariant kinematics using real beam 4-vectors when available: + // k0 (beam e), P (beam p), k' (scattered e) + // q = k0 - k' + // Q2 = -q^2 (>0) + // x = Q2 / (2 P·q) + // W^2 = M_p^2 + 2 P·q - Q2 + .Define("Q2_mc", [](double px, double py, double pz, double E, + double bex,double bey,double bez,double beE){ + const double qE = beE - E; + const double qx = bex - px; + const double qy = bey - py; + const double qz = bez - pz; + const double q2 = qE*qE - (qx*qx + qy*qy + qz*qz); + return -q2; }, {"px_mc","py_mc","pz_mc","E_mc","beam_e_px","beam_e_py","beam_e_pz","beam_e_E"}) + .Define("W_mc", [](double Q2,double px,double py,double pz,double E, + double bex,double bey,double bez,double beE, + double bpx,double bpy,double bpz,double bpE){ + const double qE = beE - E; + const double qx = bex - px; + const double qy = bey - py; + const double qz = bez - pz; + const double Pdotq = bpE*qE - (bpx*qx + bpy*qy + bpz*qz); + const double Mp = 0.9382720813; // GeV + const double W2 = Mp*Mp + 2.0*Pdotq - Q2; + return W2>0 ? std::sqrt(W2) : 0.0; }, {"Q2_mc","px_mc","py_mc","pz_mc","E_mc","beam_e_px","beam_e_py","beam_e_pz","beam_e_E","beam_p_px","beam_p_py","beam_p_pz","beam_p_E"}) + .Define("log10Q2_mc", "Q2_mc>0 ? std::log10(Q2_mc) : -10.0") + .Define("x_bj", [](double Q2,double px,double py,double pz,double E, + double bex,double bey,double bez,double beE, + double bpx,double bpy,double bpz,double bpE){ + const double qE = beE - E; + const double qx = bex - px; + const double qy = bey - py; + const double qz = bez - pz; + const double Pdotq = bpE*qE - (bpx*qx + bpy*qy + bpz*qz); + const double denom = 2.0*Pdotq; + return (Q2>0.0 && denom>0.0) ? Q2/denom : 1e-10; }, {"Q2_mc","px_mc","py_mc","pz_mc","E_mc","beam_e_px","beam_e_py","beam_e_pz","beam_e_E","beam_p_px","beam_p_py","beam_p_pz","beam_p_E"}) + .Define("log10x_mc", "x_bj>0 ? std::log10(x_bj) : -10.0"); + + // Denominator: events with a valid scattered electron + auto denomDF = baseDF.Filter("hasElectron"); + // Numerator: additionally require single reconstructed particle + auto filterDF = denomDF.Filter(particleCollectionName+".size()==1"); // Plot x,y,z momentum resolution as a function of the MCParticle momentum auto momentumDF = filterDF .Define("reco_momentum", particleCollectionName+"[0].momentum") - .Define("mc_momentum", "SimParticles[0].momentum") .Define("px_rec", "reco_momentum.x") .Define("py_rec", "reco_momentum.y") .Define("pz_rec", "reco_momentum.z") - .Define("px_mc", "mc_momentum.x") - .Define("py_mc", "mc_momentum.y") - .Define("pz_mc", "mc_momentum.z") - // Calculate theta and phi for reco and MC - .Define("theta_rec", "std::atan2(std::sqrt(px_rec*px_rec + py_rec*py_rec), pz_rec)") - .Define("theta_mc", "std::atan2(std::sqrt(px_mc*px_mc + py_mc*py_mc), pz_mc)") + // Calculate theta and phi for reco + .Define("theta_rec_raw", "std::atan2(std::sqrt(px_rec*px_rec + py_rec*py_rec), pz_rec)") + .Define("theta_rec", "TMath::Pi() - theta_rec_raw") + .Define("theta_rec_mrad", "1000.0*theta_rec") .Define("phi_rec_rad", "std::atan2(py_rec, px_rec)") - .Define("phi_mc_rad", "std::atan2(py_mc, px_mc)") .Define("phi_rec", "TMath::RadToDeg()*phi_rec_rad") - .Define("phi_mc", "TMath::RadToDeg()*phi_mc_rad") // Calculate resolutions .Define("px_diff", "(px_rec - px_mc)") .Define("py_diff", "(py_rec - py_mc)") .Define("pz_res", "(pz_rec - pz_mc)/pz_mc") .Define("E_rec", particleCollectionName+"[0].energy") - .Define("E_mc", "std::sqrt(px_mc*px_mc + py_mc*py_mc + pz_mc*pz_mc + SimParticles[0].mass*SimParticles[0].mass)") .Define("E_res", "(E_rec - E_mc)/E_mc") - .Define("theta_diff", "(theta_rec - theta_mc)") - .Define("phi_diff", "TMath::RadToDeg()*ROOT::VecOps::DeltaPhi(phi_rec_rad, phi_mc_rad)"); + .Define("theta_diff", "(theta_rec - theta_mc)") + .Define("theta_diff_mrad", "1000.0*theta_diff") + .Define("phi_diff", "TMath::RadToDeg()*ROOT::VecOps::DeltaPhi(phi_rec_rad, phi_mc_rad)") + // Reconstructed Q2 and x using ideal beam kinematics (reconstruction doesn't have access to MC beam particles) + .Define("Q2_rec", [electronEnergy](float px, float py, float pz, float E){ + const double me = 0.00051099895; // GeV + const double Ee = electronEnergy; + const double pe = std::sqrt(std::max(0.0, Ee*Ee - me*me)); + const double dE = Ee - E; + const double dqx = -px; + const double dqy = -py; + const double dqz = -pe - pz; + const double q2 = (dE*dE - (dqx*dqx + dqy*dqy + dqz*dqz)); + return -q2; + }, {"px_rec","py_rec","pz_rec","E_rec"}) + .Define("x_bj_rec", [electronEnergy, protonEnergy](double Q2, float px, float py, float pz, float E){ + const double Mp = 0.9382720813; // GeV + const double me = 0.00051099895; // GeV + const double Ee = electronEnergy; + const double Ep = protonEnergy; + const double pe = std::sqrt(std::max(0.0, Ee*Ee - me*me)); + const double pp = std::sqrt(std::max(0.0, Ep*Ep - Mp*Mp)); + const double qE = Ee - E; + const double qx = -px; + const double qy = -py; + const double qz = -pe - pz; + const double Pdotq = Ep*qE - (0.0*qx + 0.0*qy + pp*qz); + const double denom = 2.0 * Pdotq; + return (Q2 > 0.0 && denom > 0.0) ? (Q2 / denom) : 1e-10; + }, {"Q2_rec","px_rec","py_rec","pz_rec","E_rec"}) + .Define("log10Q2_rec", "Q2_rec>0 ? std::log10(Q2_rec) : -10.0") + .Define("log10x_rec", "x_bj_rec>0 ? std::log10(x_bj_rec) : -10.0") + .Define("Q2_diff", "(Q2_rec - Q2_mc)") + .Define("x_diff", "(x_bj_rec - x_bj)") + .Define("log10Q2_diff", "(log10Q2_rec - log10Q2_mc)") + .Define("log10x_diff", "(log10x_rec - log10x_mc)") + .Define("Q2_rel_resolution", "Q2_mc > 0 ? Q2_diff / Q2_mc : 1.0"); //Print the size of the original DataFrame std::cout << "Original DataFrame size: " << d0.Count().GetValue() << std::endl; @@ -72,7 +196,7 @@ void reconstructionAnalysis(TString inFile = "/home/simong/EIC/detec int momentumBins = 100; float momentumXRange[2] = {-0.1, 0.1}; float momentumYRange[2] = {-0.1, 0.1}; - float momentumZRange[2] = {-beamEnergy, 0}; + float momentumZRange[2] = {-electronEnergy, 0.0f}; int momentumResolutionBins = 100; float momentumDifferenceXRange[2] = {-0.05, 0.05}; @@ -82,19 +206,38 @@ void reconstructionAnalysis(TString inFile = "/home/simong/EIC/detec int energyBins = 100; int thetaBins = 100; int phiBins = 100; - float energyRange[2] = {3.5, beamEnergy}; // GeV - float thetaRange[2] = {3.134, TMath::Pi()}; // radians from 3.1 to pi + float energyRange[2] = {0.0f, electronEnergy}; + float thetaRange[2] = {0.0, 8.0}; // scattering angle range (pi - theta_raw) in mrad (forward angles) float phiRange[2] = {-180, 180}; // degrees from -180 to 180 int resolutionBins = 100; float energyResolutionRange[2] = {-0.1, 0.1}; - float thetaResolutionRange[2] = {-0.003, 0.003}; + float thetaResolutionRange[2] = {-3.0, 3.0}; // mrad float phiResolutionRange[2] = {-90, 90}; // degrees from -90 to 90 + // Q2 and x binning configuration + int Q2Bins = 200; + float Q2Range[2] = {0.0, 1.0}; + int xBins = 200; + float xRange[2] = {0.0, 0.1}; + int log10Q2Bins = 200; + float log10Q2Range[2] = {-9.0, 0.0}; + int log10xBins = 200; + float log10xRange[2] = {-13.0, 0.0}; + int Q2DiffBins = 200; + float Q2DiffRange[2] = {-0.2, 0.2}; + int xDiffBins = 200; + float xDiffRange[2] = {-0.02, 0.02}; + int log10DiffBins = 200; + float log10Q2DiffRange[2] = {-3, 3}; + float log10xDiffRange[2] = {-3, 3}; + int Q2RelResBins = 200; + float Q2RelResRange[2] = {-2.0, 2.0}; // 0 to 200% relative resolution + // Plot reconstructed vs montecarlo momentum components - auto px_Hist = momentumDF.Histo2D({"px_vs_px", "Reconstructed vs MC px; px reconstructed [GeV]; px MC [GeV]", momentumBins, momentumXRange[0], momentumXRange[1], momentumBins, momentumXRange[0], momentumXRange[1]}, "px_rec", "px_mc"); - auto py_Hist = momentumDF.Histo2D({"py_vs_py", "Reconstructed vs MC py; py reconstructed [GeV]; py MC [GeV]", momentumBins, momentumYRange[0], momentumYRange[1], momentumBins, momentumYRange[0], momentumYRange[1]}, "py_rec", "py_mc"); - auto pz_Hist = momentumDF.Histo2D({"pz_vs_pz", "Reconstructed vs MC pz; pz reconstructed [GeV]; pz MC [GeV]", momentumBins, momentumZRange[0], momentumZRange[1], momentumBins, momentumZRange[0], momentumZRange[1]}, "pz_rec", "pz_mc"); + auto px_Hist = momentumDF.Histo2D({"px_vs_px", "Reconstructed vs MC px; px MC [GeV]; px reconstructed [GeV]", momentumBins, momentumXRange[0], momentumXRange[1], momentumBins, momentumXRange[0], momentumXRange[1]}, "px_mc", "px_rec"); + auto py_Hist = momentumDF.Histo2D({"py_vs_py", "Reconstructed vs MC py; py MC [GeV]; py reconstructed [GeV]", momentumBins, momentumYRange[0], momentumYRange[1], momentumBins, momentumYRange[0], momentumYRange[1]}, "py_mc", "py_rec"); + auto pz_Hist = momentumDF.Histo2D({"pz_vs_pz", "Reconstructed vs MC pz; pz MC [GeV]; pz reconstructed [GeV]", momentumBins, momentumZRange[0], momentumZRange[1], momentumBins, momentumZRange[0], momentumZRange[1]}, "pz_mc", "pz_rec"); // Plot individual momentum resolutions auto px_diff_Hist = momentumDF.Histo1D({"px_diff", "px difference; px difference [GeV]; Entries", momentumResolutionBins, momentumDifferenceXRange[0], momentumDifferenceXRange[1]}, "px_diff"); @@ -102,25 +245,160 @@ void reconstructionAnalysis(TString inFile = "/home/simong/EIC/detec auto pz_res_Hist = momentumDF.Histo1D({"pz_res", "pz resolution; pz resolution [GeV]; Entries", momentumResolutionBins, momentumResolutionZRange[0], momentumResolutionZRange[1]}, "pz_res"); // Plot reconstructed vs montecarlo energy, theta and phi - auto E_Hist = momentumDF.Histo2D({"E_vs_E", "Reconstructed vs MC energy; E reconstructed [GeV]; E MC [GeV]", energyBins, energyRange[0], energyRange[1], energyBins, energyRange[0], energyRange[1]}, "E_rec", "E_mc"); - auto theta_Hist = momentumDF.Histo2D({"theta_vs_theta", "Reconstructed vs MC theta; theta reconstructed [rad]; theta MC [rad]", thetaBins, thetaRange[0], thetaRange[1], thetaBins, thetaRange[0], thetaRange[1]}, "theta_rec", "theta_mc"); - auto phi_Hist = momentumDF.Histo2D({"phi_vs_phi", "Reconstructed vs MC phi; phi reconstructed [deg]; phi MC [deg]", phiBins, phiRange[0], phiRange[1], phiBins, phiRange[0], phiRange[1]}, "phi_rec", "phi_mc"); + auto E_Hist = momentumDF.Histo2D({"E_vs_E", "Reconstructed vs MC energy; E MC [GeV]; E reconstructed [GeV]", energyBins, energyRange[0], energyRange[1], energyBins, energyRange[0], energyRange[1]}, "E_mc", "E_rec"); + auto theta_Hist = momentumDF.Histo2D({"theta_vs_theta", "Reconstructed vs MC scattering angle; scattering angle MC [mrad]; scattering angle reco [mrad]", thetaBins, thetaRange[0], thetaRange[1], thetaBins, thetaRange[0], thetaRange[1]}, "theta_mc_mrad", "theta_rec_mrad"); + auto phi_Hist = momentumDF.Histo2D({"phi_vs_phi", "Reconstructed vs MC phi; phi MC [deg]; phi reconstructed [deg]", phiBins, phiRange[0], phiRange[1], phiBins, phiRange[0], phiRange[1]}, "phi_mc", "phi_rec"); auto E_res_Hist = momentumDF.Histo1D({"E_res", "E resolution; E resolution [GeV]; Entries", resolutionBins, energyResolutionRange[0], energyResolutionRange[1]}, "E_res"); - auto theta_diff_Hist = momentumDF.Histo1D({"theta_diff", "theta difference; theta difference [rad]; Entries", resolutionBins, thetaResolutionRange[0], thetaResolutionRange[1]}, "theta_diff"); + auto theta_diff_Hist = momentumDF.Histo1D({"theta_diff", "theta difference; theta difference [mrad]; Entries", resolutionBins, thetaResolutionRange[0], thetaResolutionRange[1]}, "theta_diff_mrad"); auto phi_diff_Hist = momentumDF.Histo1D({"phi_diff", "phi difference; phi difference [deg]; Entries", resolutionBins, phiResolutionRange[0], phiResolutionRange[1]}, "phi_diff"); // Plot Reconstructed energy, theta and phi resolutions as a function of each reconstructed value of energy, thata and phi auto E_res_vs_E_Hist = momentumDF.Histo2D({"E_res_vs_E", "E resolution vs E reconstructed; E reconstructed [GeV]; E resolution [GeV]", energyBins, energyRange[0], energyRange[1], resolutionBins, energyResolutionRange[0], energyResolutionRange[1]}, "E_rec", "E_res"); - auto E_res_vs_theta_Hist = momentumDF.Histo2D({"E_res_vs_theta", "E resolution vs theta reconstructed; theta reconstructed [rad]; E resolution [GeV]", thetaBins, thetaRange[0], thetaRange[1], resolutionBins, energyResolutionRange[0], energyResolutionRange[1]}, "theta_rec", "E_res"); + auto E_res_vs_theta_Hist = momentumDF.Histo2D({"E_res_vs_theta", "E resolution vs theta reconstructed; theta reconstructed [mrad]; E resolution [GeV]", thetaBins, thetaRange[0], thetaRange[1], resolutionBins, energyResolutionRange[0], energyResolutionRange[1]}, "theta_rec_mrad", "E_res"); auto E_res_vs_phi_Hist = momentumDF.Histo2D({"E_res_vs_phi", "E resolution vs phi reconstructed; phi reconstructed [deg]; E resolution [GeV]", phiBins, phiRange[0], phiRange[1], resolutionBins, energyResolutionRange[0], energyResolutionRange[1]}, "phi_rec", "E_res"); - auto theta_diff_vs_E_Hist = momentumDF.Histo2D({"theta_diff_vs_E", "theta difference vs E reconstructed; E reconstructed [GeV]; theta difference [rad]", energyBins, energyRange[0], energyRange[1], resolutionBins, thetaResolutionRange[0], thetaResolutionRange[1]}, "E_rec", "theta_diff"); - auto theta_diff_vs_theta_Hist = momentumDF.Histo2D({"theta_diff_vs_theta", "theta difference vs theta reconstructed; theta reconstructed [rad]; theta difference [rad]", thetaBins, thetaRange[0], thetaRange[1], resolutionBins, thetaResolutionRange[0], thetaResolutionRange[1]}, "theta_rec", "theta_diff"); - auto theta_diff_vs_phi_Hist = momentumDF.Histo2D({"theta_diff_vs_phi", "theta difference vs phi reconstructed; phi reconstructed [deg]; theta difference [rad]", phiBins, phiRange[0], phiRange[1], resolutionBins, thetaResolutionRange[0], thetaResolutionRange[1]}, "phi_rec", "theta_diff"); - auto phi_diff_vs_E_Hist = momentumDF.Histo2D({"phi_diff_vs_E", "phi difference vs E reconstructed; E reconstructed [GeV]; phi difference [rad]", energyBins, energyRange[0], energyRange[1], resolutionBins, phiResolutionRange[0], phiResolutionRange[1]}, "E_rec", "phi_diff"); - auto phi_diff_vs_theta_Hist = momentumDF.Histo2D({"phi_diff_vs_theta", "phi difference vs theta reconstructed; theta reconstructed [rad]; phi difference [deg]", thetaBins, thetaRange[0], thetaRange[1], resolutionBins, phiResolutionRange[0], phiResolutionRange[1]}, "theta_rec", "phi_diff"); + auto theta_diff_vs_E_Hist = momentumDF.Histo2D({"theta_diff_vs_E", "theta difference vs E reconstructed; E reconstructed [GeV]; theta difference [mrad]", energyBins, energyRange[0], energyRange[1], resolutionBins, thetaResolutionRange[0], thetaResolutionRange[1]}, "E_rec", "theta_diff_mrad"); + auto theta_diff_vs_theta_Hist = momentumDF.Histo2D({"theta_diff_vs_theta", "theta difference vs theta reconstructed; theta reconstructed [mrad]; theta difference [mrad]", thetaBins, thetaRange[0], thetaRange[1], resolutionBins, thetaResolutionRange[0], thetaResolutionRange[1]}, "theta_rec_mrad", "theta_diff_mrad"); + auto theta_diff_vs_phi_Hist = momentumDF.Histo2D({"theta_diff_vs_phi", "theta difference vs phi reconstructed; phi reconstructed [deg]; theta difference [mrad]", phiBins, phiRange[0], phiRange[1], resolutionBins, thetaResolutionRange[0], thetaResolutionRange[1]}, "phi_rec", "theta_diff_mrad"); + auto phi_diff_vs_E_Hist = momentumDF.Histo2D({"phi_diff_vs_E", "phi difference vs E reconstructed; E reconstructed [GeV]; phi difference [deg]", energyBins, energyRange[0], energyRange[1], resolutionBins, phiResolutionRange[0], phiResolutionRange[1]}, "E_rec", "phi_diff"); + auto phi_diff_vs_theta_Hist = momentumDF.Histo2D({"phi_diff_vs_theta", "phi difference vs theta reconstructed; theta reconstructed [mrad]; phi difference [deg]", thetaBins, thetaRange[0], thetaRange[1], resolutionBins, phiResolutionRange[0], phiResolutionRange[1]}, "theta_rec_mrad", "phi_diff"); auto phi_diff_vs_phi_Hist = momentumDF.Histo2D({"phi_diff_vs_phi", "phi difference vs phi reconstructed; phi reconstructed [deg]; phi difference [deg]", phiBins, phiRange[0], phiRange[1], resolutionBins, phiResolutionRange[0], phiResolutionRange[1]}, "phi_rec", "phi_diff"); + // Filtered dataframe for scattering angle > 1 mrad (0.001 rad) + auto phiThetaFilteredDF = momentumDF.Filter("theta_mc > 0.001"); + auto phi_filtered_Hist = phiThetaFilteredDF.Histo2D({"phi_filtered_vs_phi", "Reco vs MC phi (theta>1mrad); phi MC [deg]; phi reco [deg]", phiBins, phiRange[0], phiRange[1], phiBins, phiRange[0], phiRange[1]}, "phi_mc", "phi_rec"); + auto phi_filtered_diff_Hist = phiThetaFilteredDF.Histo1D({"phi_filtered_diff", "Phi difference (theta>1mrad); phi difference [deg]; Entries", resolutionBins, phiResolutionRange[0], phiResolutionRange[1]}, "phi_diff"); + + // Q2 and x resolution plots + auto Q2_Hist = momentumDF.Histo2D({"Q2_vs_Q2", "Reconstructed vs MC Q^{2}; Q^{2} MC [GeV^{2}]; Q^{2} reconstructed [GeV^{2}]", Q2Bins, Q2Range[0], Q2Range[1], Q2Bins, Q2Range[0], Q2Range[1]}, "Q2_mc", "Q2_rec"); + auto x_Hist = momentumDF.Histo2D({"x_vs_x", "Reconstructed vs MC x; x MC; x reconstructed", xBins, xRange[0], xRange[1], xBins, xRange[0], xRange[1]}, "x_bj", "x_bj_rec"); + auto log10Q2_comp_Hist = momentumDF.Histo2D({"log10Q2_vs_log10Q2", "Reconstructed vs MC log10(Q^{2}); log10(Q^{2}) MC; log10(Q^{2}) reconstructed", log10Q2Bins, log10Q2Range[0], log10Q2Range[1], log10Q2Bins, log10Q2Range[0], log10Q2Range[1]}, "log10Q2_mc", "log10Q2_rec"); + auto log10x_comp_Hist = momentumDF.Histo2D({"log10x_vs_log10x", "Reconstructed vs MC log10(x); log10(x) MC; log10(x) reconstructed", log10xBins, log10xRange[0], log10xRange[1], log10xBins, log10xRange[0], log10xRange[1]}, "log10x_mc", "log10x_rec"); + + auto Q2_diff_Hist = momentumDF.Histo1D({"Q2_diff", "Q^{2} difference; Q^{2} difference [GeV^{2}]; Entries", Q2DiffBins, Q2DiffRange[0], Q2DiffRange[1]}, "Q2_diff"); + auto x_diff_Hist = momentumDF.Histo1D({"x_diff", "x difference; x difference; Entries", xDiffBins, xDiffRange[0], xDiffRange[1]}, "x_diff"); + auto log10Q2_diff_Hist = momentumDF.Histo1D({"log10Q2_diff", "log10(Q^{2}) difference; log10(Q^{2}) difference; Entries", log10DiffBins, log10Q2DiffRange[0], log10Q2DiffRange[1]}, "log10Q2_diff"); + auto log10x_diff_Hist = momentumDF.Histo1D({"log10x_diff", "log10(x) difference; log10(x) difference; Entries", log10DiffBins, log10xDiffRange[0], log10xDiffRange[1]}, "log10x_diff"); + + auto Q2_diff_vs_Q2_Hist = momentumDF.Histo2D({"Q2_diff_vs_Q2", "Q^{2} difference vs MC Q^{2}; Q^{2} MC [GeV^{2}]; Q^{2} difference [GeV^{2}]", Q2Bins, Q2Range[0], Q2Range[1], Q2DiffBins, Q2DiffRange[0], Q2DiffRange[1]}, "Q2_rec", "Q2_diff"); + auto x_diff_vs_x_Hist = momentumDF.Histo2D({"x_diff_vs_x", "x difference vs MC x; x MC; x difference", xBins, xRange[0], xRange[1], xDiffBins, xDiffRange[0], xDiffRange[1]}, "x_bj", "x_diff"); + auto log10Q2_diff_vs_log10Q2_Hist = momentumDF.Histo2D({"log10Q2_diff_vs_log10Q2", "log10(Q^{2}) difference vs MC log10(Q^{2}); log10(Q^{2}) MC; log10(Q^{2}) difference", log10Q2Bins, log10Q2Range[0], log10Q2Range[1], log10DiffBins, log10Q2DiffRange[0], log10Q2DiffRange[1]}, "log10Q2_rec", "log10Q2_diff"); + auto log10x_diff_vs_log10x_Hist = momentumDF.Histo2D({"log10x_diff_vs_log10x", "log10(x) difference vs MC log10(x); log10(x) MC; log10(x) difference", log10xBins, log10xRange[0], log10xRange[1], log10DiffBins, log10xDiffRange[0], log10xDiffRange[1]}, "log10x_rec", "log10x_diff"); + + // Q2 relative resolution histogram for threshold determination + auto Q2_rel_res_vs_log10Q2_Hist = momentumDF.Histo2D({"Q2_rel_res_vs_log10Q2", "Q^{2} Relative Resolution vs log_{10}(Q^{2}) MC; log_{10}(Q^{2}) MC [GeV^{2}]; Q^{2} Relative Resolution", log10Q2Bins, log10Q2Range[0], log10Q2Range[1], Q2RelResBins, Q2RelResRange[0], Q2RelResRange[1]}, "log10Q2_rec", "Q2_rel_resolution"); + + // New Q2 resolution analysis histograms + auto x_vs_Q2_all_Hist = momentumDF.Histo2D({"x_vs_Q2_all", "x vs Q^{2} (all events); Q^{2}_{MC} [GeV^{2}]; x_{MC}", Q2Bins, Q2Range[0], Q2Range[1], xBins, xRange[0], xRange[1]}, "Q2_mc", "x_bj"); + + // Acceptance histograms (denominator vs numerator) + auto E_all_Hist = denomDF.Histo1D({"E_all", "MC Electron Energy; E [GeV]; Entries", energyBins, energyRange[0], energyRange[1]}, "E_mc"); + auto theta_all_Hist = denomDF.Histo1D({"theta_all", "MC Electron Scattering Angle; scattering angle [mrad]; Entries", thetaBins, thetaRange[0], thetaRange[1]}, "theta_mc_mrad"); + auto phi_all_Hist = denomDF.Histo1D({"phi_all", "MC Electron Phi; phi [deg]; Entries", phiBins, phiRange[0], phiRange[1]}, "phi_mc"); + auto log10Q2_all_Hist = denomDF.Histo1D({"log10Q2_all", "MC log10(Q^{2}); log10(Q^{2}) [GeV^{2}]; Entries", log10Q2Bins, log10Q2Range[0], log10Q2Range[1]}, "log10Q2_mc"); + auto W_all_Hist = denomDF.Histo1D({"W_all", "MC W; W [GeV]; Entries", energyBins, 0.0, 5.0}, "W_mc"); + + auto E_acc_Hist = filterDF.Histo1D({"E_acc", "Accepted Electron Energy; E [GeV]; Entries", energyBins, energyRange[0], energyRange[1]}, "E_mc"); + auto theta_acc_Hist = filterDF.Histo1D({"theta_acc", "Accepted Electron Scattering Angle; scattering angle [mrad]; Entries", thetaBins, thetaRange[0], thetaRange[1]}, "theta_mc_mrad"); + auto phi_acc_Hist = filterDF.Histo1D({"phi_acc", "Accepted Electron Phi; phi [deg]; Entries", phiBins, phiRange[0], phiRange[1]}, "phi_mc"); + auto log10Q2_acc_Hist = filterDF.Histo1D({"log10Q2_acc", "Accepted log10(Q^{2}); log10(Q^{2}) [GeV^{2}]; Entries", log10Q2Bins, log10Q2Range[0], log10Q2Range[1]}, "log10Q2_mc"); + auto W_acc_Hist = filterDF.Histo1D({"W_acc", "Accepted W; W [GeV]; Entries", energyBins, 0.0, 5.0}, "W_mc"); + + // 2D acceptance: log10Q2 vs E + auto log10Q2_vs_E_all_Hist = denomDF.Histo2D({"log10Q2_vs_E_all", "MC log10(Q^{2}) vs E; E [GeV]; log10(Q^{2})", energyBins, energyRange[0], energyRange[1], log10Q2Bins, log10Q2Range[0], log10Q2Range[1]}, "E_mc", "log10Q2_mc"); + auto log10Q2_vs_E_acc_Hist = filterDF.Histo2D({"log10Q2_vs_E_acc", "Accepted log10(Q^{2}) vs E; E [GeV]; log10(Q^{2})", energyBins, energyRange[0], energyRange[1], log10Q2Bins, log10Q2Range[0], log10Q2Range[1]}, "E_mc", "log10Q2_mc"); + + // 2D acceptance: log10Q2 vs log10x + auto log10Q2_vs_log10x_all_Hist = denomDF.Histo2D({"log10Q2_vs_log10x_all", "MC log10(Q^{2}) vs log10(x); log10(x); log10(Q^{2})", log10xBins, log10xRange[0], log10xRange[1], log10Q2Bins, log10Q2Range[0], log10Q2Range[1]}, "log10x_mc", "log10Q2_mc"); + auto log10Q2_vs_log10x_acc_Hist = filterDF.Histo2D({"log10Q2_vs_log10x_acc", "Accepted log10(Q^{2}) vs log10(x); log10(x); log10(Q^{2})", log10xBins, log10xRange[0], log10xRange[1], log10Q2Bins, log10Q2Range[0], log10Q2Range[1]}, "log10x_mc", "log10Q2_mc"); + + TH2D* hLog10Q2_vs_log10x_acceptance = (TH2D*)log10Q2_vs_log10x_acc_Hist->Clone("hLog10Q2_vs_log10x_acceptance"); + hLog10Q2_vs_log10x_acceptance->Divide(log10Q2_vs_log10x_all_Hist.GetPtr()); + + TH1D* hE_acceptance = (TH1D*)E_acc_Hist->Clone("hE_acceptance"); hE_acceptance->Divide(E_all_Hist.GetPtr()); + TH1D* hTheta_acceptance = (TH1D*)theta_acc_Hist->Clone("hTheta_acceptance"); hTheta_acceptance->Divide(theta_all_Hist.GetPtr()); + TH1D* hPhi_acceptance = (TH1D*)phi_acc_Hist->Clone("hPhi_acceptance"); hPhi_acceptance->Divide(phi_all_Hist.GetPtr()); + TH1D* hLog10Q2_acceptance = (TH1D*)log10Q2_acc_Hist->Clone("hLog10Q2_acceptance"); hLog10Q2_acceptance->Divide(log10Q2_all_Hist.GetPtr()); + TH1D* hW_acceptance = (TH1D*)W_acc_Hist->Clone("hW_acceptance"); hW_acceptance->Divide(W_all_Hist.GetPtr()); + + // 2D acceptance + TH2D* hLog10Q2_vs_E_acceptance = (TH2D*)log10Q2_vs_E_acc_Hist->Clone("hLog10Q2_vs_E_acceptance"); + hLog10Q2_vs_E_acceptance->Divide(log10Q2_vs_E_all_Hist.GetPtr()); + + auto styleAcc = [](TH1* h,const char* title){ + h->SetTitle(title); + h->GetYaxis()->SetTitle("Acceptance"); + h->SetMinimum(0.0); + h->SetMaximum(1.05); + h->SetMarkerStyle(20); + h->SetMarkerColor(kBlue+1); + h->SetLineColor(kBlue+1); + }; + styleAcc(hE_acceptance, "Electron Acceptance vs E; E [GeV]"); + styleAcc(hTheta_acceptance, "Electron Acceptance vs scattering angle; scattering angle [mrad]"); + styleAcc(hPhi_acceptance, "Electron Acceptance vs phi; phi [deg]"); + styleAcc(hLog10Q2_acceptance, "Electron Acceptance vs log10(Q^{2}); log10(Q^{2})"); + styleAcc(hW_acceptance, "Electron Acceptance vs W; W [GeV]"); + + // Style 2D acceptance + hLog10Q2_vs_E_acceptance->SetTitle("Electron Acceptance: log10(Q^{2}) vs E; E [GeV]; log10(Q^{2})"); + hLog10Q2_vs_E_acceptance->GetZaxis()->SetTitle("Acceptance"); + hLog10Q2_vs_E_acceptance->SetMinimum(0.0); + hLog10Q2_vs_E_acceptance->SetMaximum(1.0); + + hLog10Q2_vs_log10x_acceptance->SetTitle("Electron Acceptance: log10(Q^{2}) vs log10(x); log10(x); log10(Q^{2})"); + hLog10Q2_vs_log10x_acceptance->GetZaxis()->SetTitle("Acceptance"); + hLog10Q2_vs_log10x_acceptance->SetMinimum(0.0); + hLog10Q2_vs_log10x_acceptance->SetMaximum(1.0); + + TCanvas *cAcceptance1D = new TCanvas("acceptance_1D_canvas", "Electron Acceptance (1D)", 2400, 1600); + cAcceptance1D->Divide(3,2); + cAcceptance1D->cd(1); hE_acceptance->Draw("E1"); + cAcceptance1D->cd(2); hTheta_acceptance->Draw("E1"); + cAcceptance1D->cd(3); hPhi_acceptance->Draw("E1"); + cAcceptance1D->cd(4); hLog10Q2_acceptance->Draw("E1"); + cAcceptance1D->cd(5); hW_acceptance->Draw("E1"); + cAcceptance1D->SetGrid(); + cAcceptance1D->Update(); + cAcceptance1D->SaveAs(acceptanceCanvasName); + + TCanvas *cAcceptanceXQ2 = new TCanvas("acceptance_xQ2_canvas", "Electron Acceptance: log10(Q^{2}) vs log10(x)", 1200, 1000); + cAcceptanceXQ2->cd(); + gStyle->SetOptStat(0); + hLog10Q2_vs_log10x_acceptance->Draw("COLZ"); + cAcceptanceXQ2->SetRightMargin(0.15); + cAcceptanceXQ2->Update(); + cAcceptanceXQ2->SaveAs(xQ2acceptanceCanvasName); + + // Create canvas for Q2 and x resolution + TCanvas *cQ2xResolution = new TCanvas("Q2_x_resolution_canvas", "Q^{2} and x Resolution", 3000, 1600); + cQ2xResolution->Divide(4, 2); + cQ2xResolution->cd(1); + log10Q2_comp_Hist->Draw("colz"); + gPad->SetLogz(); + cQ2xResolution->cd(2); + log10Q2_diff_Hist->Draw(); + cQ2xResolution->cd(3); + log10Q2_diff_vs_log10Q2_Hist->Draw("colz"); + gPad->SetLogz(); + cQ2xResolution->cd(4); + Q2_Hist->Draw("colz"); + gPad->SetLogz(); + cQ2xResolution->cd(5); + log10x_comp_Hist->Draw("colz"); + gPad->SetLogz(); + cQ2xResolution->cd(6); + log10x_diff_Hist->Draw(); + cQ2xResolution->cd(7); + log10x_diff_vs_log10x_Hist->Draw("colz"); + gPad->SetLogz(); + cQ2xResolution->cd(8); + x_Hist->Draw("colz"); + gPad->SetLogz(); + cQ2xResolution->SetGrid(); + cQ2xResolution->Update(); + cQ2xResolution->SaveAs(Q2xResolutionCanvasName); + // Create canvas for momentum component plots TCanvas *cMomentum = new TCanvas("momentum_canvas", "Momentum Resolution", 3000, 1600); cMomentum->Divide(3, 2); @@ -143,30 +421,7 @@ void reconstructionAnalysis(TString inFile = "/home/simong/EIC/detec cMomentum->Update(); // Save the canvas as a PNG file cMomentum->SaveAs(momentumCanvasName); - - // Create canvas for energy, theta and phi resolution plots - TCanvas *cEnergyThetaPhi = new TCanvas("energy_theta_phi_canvas", "Energy, Theta and Phi Resolution", 3000, 1600); - cEnergyThetaPhi->Divide(3, 2); - cEnergyThetaPhi->cd(1); - E_Hist->Draw("colz"); - gPad->SetLogz(); - cEnergyThetaPhi->cd(2); - theta_Hist->Draw("colz"); - gPad->SetLogz(); - cEnergyThetaPhi->cd(3); - phi_Hist->Draw("colz"); - gPad->SetLogz(); - cEnergyThetaPhi->cd(4); - E_res_Hist->Draw(); - cEnergyThetaPhi->cd(5); - theta_diff_Hist->Draw(); - cEnergyThetaPhi->cd(6); - phi_diff_Hist->Draw(); - cEnergyThetaPhi->SetGrid(); - cEnergyThetaPhi->Update(); - // Save the canvas as a PNG file - cEnergyThetaPhi->SaveAs(energyThetaPhiCanvasName); - + // Create canvas for resolution vs MC values TCanvas *cResolutionVsMC = new TCanvas("resolution_vs_mc_canvas", "Resolution vs MC Values", 3000, 1600); cResolutionVsMC->Divide(3, 3); @@ -202,13 +457,15 @@ void reconstructionAnalysis(TString inFile = "/home/simong/EIC/detec // Save the canvas as a PNG file cResolutionVsMC->SaveAs(relationCanvasName); - // Fit Gaussians to the E vs E histogram slices + // Fit Gaussians to histogram slices including new Q2 resolution analysis TObjArray* fitE_vs_E = new TObjArray(); TObjArray* fitTheta_vs_E = new TObjArray(); TObjArray* fitPhi_vs_theta = new TObjArray(); + TObjArray* fitQ2_resolution = new TObjArray(); E_res_vs_E_Hist->FitSlicesY(nullptr, 1, -1, 0, "Q", fitE_vs_E); theta_diff_vs_E_Hist->FitSlicesY(nullptr, 1, -1, 0, "Q", fitTheta_vs_E); phi_diff_vs_theta_Hist->FitSlicesY(nullptr, 1, -1, 0, "Q", fitPhi_vs_theta); + Q2_rel_res_vs_log10Q2_Hist->FitSlicesY(nullptr, 1, -1, 0, "Q", fitQ2_resolution); // Create graphs containing the fitted means and standard deviations TH1* hE_vs_E_mean = (TH1*)fitE_vs_E->At(1); // mean values (index 1) @@ -217,10 +474,12 @@ void reconstructionAnalysis(TString inFile = "/home/simong/EIC/detec TH1* hTheta_vs_E_stddev = (TH1*)fitTheta_vs_E->At(2); // stddev values (index 2) TH1* hPhi_vs_theta_mean = (TH1*)fitPhi_vs_theta->At(1); // mean values (index 1) TH1* hPhi_vs_theta_stddev = (TH1*)fitPhi_vs_theta->At(2); // stddev values (index 2) + TH1* hQ2_resolution_mean = (TH1*)fitQ2_resolution->At(1); // Q2 bias vs log10(Q2) + TH1* hQ2_resolution_sigma = (TH1*)fitQ2_resolution->At(2); // Q2 resolution vs log10(Q2) // Create a canvas for the resolution graphs - TCanvas *cResolutionGraphs = new TCanvas("resolution_graphs_canvas", "Resolution Graphs", 1200, 800); - cResolutionGraphs->Divide(3, 2); + TCanvas *cResolutionGraphs = new TCanvas("resolution_graphs_canvas", "Resolution Graphs", 1800, 1200); + cResolutionGraphs->Divide(4, 2); cResolutionGraphs->cd(1); hE_vs_E_mean->SetTitle("Mean Energy Offset vs E MC; Energy MC [GeV]; Mean Energy Offset [GeV]"); hE_vs_E_mean->SetMarkerStyle(20); @@ -236,33 +495,45 @@ void reconstructionAnalysis(TString inFile = "/home/simong/EIC/detec hE_vs_E_stddev->SetMinimum(0); // Adjust minimum for better visibility hE_vs_E_stddev->Draw(); cResolutionGraphs->cd(2); - hTheta_vs_E_mean->SetTitle("Mean Theta Offset vs E MC; Energy MC [GeV]; Mean Theta Offset [rad]"); + hTheta_vs_E_mean->SetTitle("Mean Theta Offset vs E MC; Energy MC [GeV]; Mean Theta Offset [mrad]"); hTheta_vs_E_mean->SetMarkerStyle(20); hTheta_vs_E_mean->SetMarkerColor(kBlue); - hTheta_vs_E_mean->SetMaximum(0.0003); // Adjust maximum for better visibility - hTheta_vs_E_mean->SetMinimum(-0.0003); // Adjust minimum for better visibility + hTheta_vs_E_mean->SetMaximum(0.3); // Adjust maximum for better visibility + hTheta_vs_E_mean->SetMinimum(-0.3); // Adjust minimum for better visibility hTheta_vs_E_mean->Draw(); cResolutionGraphs->cd(5); - hTheta_vs_E_stddev->SetTitle("Std Dev Theta Offset vs E MC; Energy MC [GeV]; Std Dev Theta Offset [rad]"); + hTheta_vs_E_stddev->SetTitle("Std Dev Theta Offset vs E MC; Energy MC [GeV]; Std Dev Theta Offset [mrad]"); hTheta_vs_E_stddev->SetMarkerStyle(20); hTheta_vs_E_stddev->SetMarkerColor(kRed); - hTheta_vs_E_stddev->SetMaximum(0.0005); // Adjust maximum for better visibility + hTheta_vs_E_stddev->SetMaximum(0.5); // Adjust maximum for better visibility hTheta_vs_E_stddev->SetMinimum(0); // Adjust minimum for better visibility hTheta_vs_E_stddev->Draw(); cResolutionGraphs->cd(3); - hPhi_vs_theta_mean->SetTitle("Mean Phi Offset vs theta MC; theta MC [rad]; Mean Phi Offset [rad]"); + hPhi_vs_theta_mean->SetTitle("Mean Phi Offset vs theta MC; theta MC [mrad]; Mean Phi Offset [deg]"); hPhi_vs_theta_mean->SetMarkerStyle(20); hPhi_vs_theta_mean->SetMarkerColor(kBlue); hPhi_vs_theta_mean->SetMaximum(5); // Adjust maximum for better visibility hPhi_vs_theta_mean->SetMinimum(-5); // Adjust minimum for better visibility hPhi_vs_theta_mean->Draw(); cResolutionGraphs->cd(6); - hPhi_vs_theta_stddev->SetTitle("Std Dev Phi Offset vs theta MC; theta MC [rad]; Std Dev Phi Offset [rad]"); + hPhi_vs_theta_stddev->SetTitle("Std Dev Phi Offset vs theta MC; theta MC [mrad]; Std Dev Phi Offset [deg]"); hPhi_vs_theta_stddev->SetMarkerStyle(20); hPhi_vs_theta_stddev->SetMarkerColor(kRed); hPhi_vs_theta_stddev->SetMaximum(60); // Adjust maximum for better visibility hPhi_vs_theta_stddev->SetMinimum(0); // Adjust minimum for better visibility hPhi_vs_theta_stddev->Draw(); + // Add Q2 resolution plots + cResolutionGraphs->cd(7); + hQ2_resolution_mean->SetTitle("Mean Q^{2} Bias vs log_{10}(Q^{2}) Reco; log_{10}(Q^{2}) Reco [GeV^{2}]; Mean Q^{2} Bias [GeV^{2}]"); + hQ2_resolution_mean->SetMarkerStyle(20); + hQ2_resolution_mean->SetMarkerColor(kBlue); + hQ2_resolution_mean->Draw(); + cResolutionGraphs->cd(8); + hQ2_resolution_sigma->SetTitle("Q^{2} Resolution vs log_{10}(Q^{2}) Reco; log_{10}(Q^{2}) Reco [GeV^{2}]; Q^{2} Resolution [GeV^{2}]"); + hQ2_resolution_sigma->SetMarkerStyle(20); + hQ2_resolution_sigma->SetMarkerColor(kRed); + hQ2_resolution_sigma->SetMinimum(0); + hQ2_resolution_sigma->Draw(); cResolutionGraphs->SetGrid(); cResolutionGraphs->Update(); // Save the canvas as a PNG file @@ -270,9 +541,14 @@ void reconstructionAnalysis(TString inFile = "/home/simong/EIC/detec TFile *f = new TFile(outFile,"RECREATE"); cMomentum->Write(); - cEnergyThetaPhi->Write(); + // cEnergyThetaPhi->Write(); + // cEnergyThetaPhiFiltered->Write(); cResolutionVsMC->Write(); cResolutionGraphs->Write(); + cAcceptance1D->Write(); + // cAcceptance2D->Write(); + cAcceptanceXQ2->Write(); + cQ2xResolution->Write(); px_Hist->Write(); py_Hist->Write(); pz_Hist->Write(); @@ -294,6 +570,29 @@ void reconstructionAnalysis(TString inFile = "/home/simong/EIC/detec phi_diff_vs_E_Hist->Write(); phi_diff_vs_theta_Hist->Write(); phi_diff_vs_phi_Hist->Write(); + phi_filtered_Hist->Write(); + phi_filtered_diff_Hist->Write(); + // Write Q2 and x resolution histograms + Q2_Hist->Write(); + x_Hist->Write(); + log10Q2_comp_Hist->Write(); + log10x_comp_Hist->Write(); + Q2_diff_Hist->Write(); + x_diff_Hist->Write(); + log10Q2_diff_Hist->Write(); + log10x_diff_Hist->Write(); + Q2_diff_vs_Q2_Hist->Write(); + x_diff_vs_x_Hist->Write(); + log10Q2_diff_vs_log10Q2_Hist->Write(); + log10x_diff_vs_log10x_Hist->Write(); + // Write acceptance related histograms + E_all_Hist->Write(); E_acc_Hist->Write(); hE_acceptance->Write(); + theta_all_Hist->Write(); theta_acc_Hist->Write(); hTheta_acceptance->Write(); + phi_all_Hist->Write(); phi_acc_Hist->Write(); hPhi_acceptance->Write(); + log10Q2_all_Hist->Write(); log10Q2_acc_Hist->Write(); hLog10Q2_acceptance->Write(); + W_all_Hist->Write(); W_acc_Hist->Write(); hW_acceptance->Write(); + log10Q2_vs_E_all_Hist->Write(); log10Q2_vs_E_acc_Hist->Write(); hLog10Q2_vs_E_acceptance->Write(); + log10Q2_vs_log10x_all_Hist->Write(); log10Q2_vs_log10x_acc_Hist->Write(); hLog10Q2_vs_log10x_acceptance->Write(); hE_vs_E_mean->Write(); hE_vs_E_stddev->Write(); @@ -301,6 +600,13 @@ void reconstructionAnalysis(TString inFile = "/home/simong/EIC/detec hTheta_vs_E_stddev->Write(); hPhi_vs_theta_mean->Write(); hPhi_vs_theta_stddev->Write(); + + // Write new Q2 resolution analysis histograms + x_vs_Q2_all_Hist->Write(); + hQ2_resolution_mean->Write(); + hQ2_resolution_sigma->Write(); + + Q2_rel_res_vs_log10Q2_Hist->Write(); f->Close();