diff --git a/analysis/data/2016_mc_ana_corrections.json b/analysis/data/2016_mc_ana_corrections.json new file mode 100644 index 000000000..5660078d5 --- /dev/null +++ b/analysis/data/2016_mc_ana_corrections.json @@ -0,0 +1,7 @@ +{ + "corrections": { + "track_time": -5.5, + "track_cluster_dt": -0.7, + "track_z0": -0.06 + } +} diff --git a/analysis/data/track_bias_corrections_data_2016.json b/analysis/data/track_bias_corrections_data_2016.json new file mode 100644 index 000000000..3f52f0ca3 --- /dev/null +++ b/analysis/data/track_bias_corrections_data_2016.json @@ -0,0 +1,4 @@ +{ + "track_time": -1.5 +} + diff --git a/analysis/data/track_bias_corrections_tritrig_2016.json b/analysis/data/track_bias_corrections_tritrig_2016.json new file mode 100644 index 000000000..cda3a44b4 --- /dev/null +++ b/analysis/data/track_bias_corrections_tritrig_2016.json @@ -0,0 +1,5 @@ +{ + "track_time": -5.5, + "track_z0": -0.06 +} + diff --git a/analysis/data/v0_projection_2016_mc_7800_config.json b/analysis/data/v0_projection_2016_mc_7800_config.json deleted file mode 100644 index cd6f5875c..000000000 --- a/analysis/data/v0_projection_2016_mc_7800_config.json +++ /dev/null @@ -1,10 +0,0 @@ -{ - "7800": { - "target_position": -4.3, - "rotated_mean_x": -2.31164e-01, - "rotated_mean_y": -5.15632e-02, - "rotated_sigma_x": 2.01293e-01, - "rotated_sigma_y": 8.16385e-02, - "rotation_angle_mrad": -120.719 - } -} diff --git a/analysis/data/v0_projection_2016_mc_config.json b/analysis/data/v0_projection_2016_mc_config.json new file mode 100644 index 000000000..9dd1d5b87 --- /dev/null +++ b/analysis/data/v0_projection_2016_mc_config.json @@ -0,0 +1,14 @@ +{ + "7984": { + "target_position": -4.3, + "rotated_mean_x": -0.22984095484241324, + "rotated_mean_y": -0.027767688133970195, + "rotated_sigma_x": 0.22379266968589445, + "rotated_sigma_y": 0.08329193763283319, + "rotation_angle_mrad": -89.95092043750218, + "unrotated_mean_x": -0.22718604898785316, + "unrotated_mean_y": -0.048091565786742764, + "unrotated_sigma_x": 0.26342442474167516, + "unrotated_sigma_y": 0.10030767283985152 + } +} \ No newline at end of file diff --git a/analysis/data/v0_projection_2016_mc_signal_config.json b/analysis/data/v0_projection_2016_mc_signal_config.json deleted file mode 100644 index 32f1a431a..000000000 --- a/analysis/data/v0_projection_2016_mc_signal_config.json +++ /dev/null @@ -1,10 +0,0 @@ -{ - "7800": { - "target_position": -4.3, - "rotated_mean_x": 0.00032589410352742634, - "rotated_mean_y": 0.0005273697881972964, - "rotated_sigma_x": 0.19701452441142028, - "rotated_sigma_y": 0.09424252765235556, - "rotation_angle_mrad": -171.21036357655905 - } -} diff --git a/analysis/include/AnaHelpers.h b/analysis/include/AnaHelpers.h index 241df6f2b..430fd76c0 100644 --- a/analysis/include/AnaHelpers.h +++ b/analysis/include/AnaHelpers.h @@ -83,6 +83,17 @@ class AnaHelpers { * @return false */ bool GetParticlesFromVtx(Vertex* vtx, Particle*& ele, Particle*& pos); + + /** + * @brief Get the Particles From Vtx object + * + * @param vtx + * @param ele1 + * @param ele2 + * @return true + * @return false + */ + bool GetSameParticlesFromVtx(Vertex* vtx, Particle*& ele1, Particle*& ele2); /** * @brief brief description diff --git a/analysis/plotconfigs/tracking/simps/vtxAnalysis_2016_simp_reach_light.json b/analysis/plotconfigs/tracking/simps/vtxAnalysis_2016_simp_reach_light.json new file mode 100644 index 000000000..9b477afd9 --- /dev/null +++ b/analysis/plotconfigs/tracking/simps/vtxAnalysis_2016_simp_reach_light.json @@ -0,0 +1,324 @@ +{ + "ele_track_cluster_EoverP_h": { + "bins" : 1000, + "minX" : 0.0, + "maxX" : 10.0, + "xtitle" : "e^{-} Track Cluster E/P", + "ytitle" : "Events" + }, + "pos_track_cluster_EoverP_h": { + "bins" : 1000, + "minX" : 0.0, + "maxX" : 10.0, + "xtitle" : "e^{+} Track Cluster E/P", + "ytitle" : "Events" + }, + "corr_eleClus_t_h": { + "bins" : 600, + "minX" : -30.0, + "maxX" : 30.0, + "xtitle" : "e^{-}_{clus} corr time [ns]", + "ytitle" : "Events" + }, + "corr_posClus_t_h": { + "bins" : 600, + "minX" : -30.0, + "maxX" : 30.0, + "xtitle" : "e^{+}_{clus} corr time [ns]", + "ytitle" : "Events" + }, + "ele_hitlayers_h": { + "bins" : 12, + "minX" : -0.5, + "maxX": 11.5, + "xtitle" : "hit layer", + "ytitle" : "nHits on Track" + }, + "pos_hitlayers_h": { + "bins" : 12, + "minX" : -0.5, + "maxX": 11.5, + "xtitle" : "hit layer", + "ytitle" : "nHits on Track" + }, + "ele_pos_clusTimeDiff_h": { + "bins" : 320, + "minX" : -16, + "maxX" : 16, + "xtitle" : "#Delta_{t}(e^{-},e^{+})_{clus} [ns]", + "ytitle" : "Events" + }, + "ele_track_clus_dt_h": { + "bins" : 200, + "minX" : -20, + "maxX" : 20, + "xtitle" : "e^{-} Track Cluster #Delta_{t} [ns]", + "ytitle" : "Events" + }, + "pos_track_clus_dt_h": { + "bins" : 200, + "minX" : -20, + "maxX" : 20, + "xtitle" : "e^{+} Track Cluster #Delta_{t} [ns]", + "ytitle" : "Events" + }, + "ele_d0_h" : { + "bins" : 200, + "minX" : -10, + "maxX" : 10, + "xtitle" : "d_{0} [mm]", + "ytitle" : "Tracks" + }, + "ele_Phi_h" : { + "bins" : 100, + "minX" : -0.4, + "maxX" : 0.4, + "xtitle" : "#phi_{0}", + "ytitle" : "Tracks" + }, + "ele_Omega_h" : { + "bins" : 100, + "minX" : -0.001, + "maxX" : 0.001, + "xtitle" : "#omega", + "ytitle" : "Tracks" + }, + "ele_TanLambda_h" : { + "bins" : 200, + "minX" : -0.2, + "maxX" : 0.2, + "xtitle" : "tan(#lambda)", + "ytitle" : "Tracks" + }, + "ele_Z0_h" : { + "bins" : 200, + "minX" : -5, + "maxX" : 5, + "xtitle" : "z_{0} [mm]", + "ytitle" : "Tracks" + }, + "ele_time_h" : { + "bins" : 200, + "minX" : -10, + "maxX" : 10, + "xtitle" : "track time [ns]", + "ytitle" : "Tracks" + }, + "ele_chi2ndf_h" : { + "bins" : 200, + "minX" : 0, + "maxX" : 30, + "xtitle" : "ele track #chi^{2} / ndf", + "ytitle" : "Tracks" + }, + "ele_chi2_h" : { + "bins" : 200, + "minX" : 0, + "maxX" : 20, + "xtitle" : "track #chi^{2}", + "ytitle" : "Tracks" + }, + "ele_p_h" : { + "bins" : 250, + "minX" : 0, + "maxX" : 2.5, + "xtitle" : "p_{e^{-}} [GeV]", + "ytitle" : "Tracks" + }, + "pos_d0_h" : { + "bins" : 200, + "minX" : -10, + "maxX" : 10, + "xtitle" : "d_{0} [mm]", + "ytitle" : "Tracks" + }, + "pos_Phi_h" : { + "bins" : 100, + "minX" : -0.4, + "maxX" : 0.4, + "xtitle" : "#phi_{0}", + "ytitle" : "Tracks" + }, + "pos_Omega_h" : { + "bins" : 100, + "minX" : -0.001, + "maxX" : 0.001, + "xtitle" : "#omega", + "ytitle" : "Tracks" + }, + "pos_TanLambda_h" : { + "bins" : 200, + "minX" : -0.2, + "maxX" : 0.2, + "xtitle" : "tan(#lambda)", + "ytitle" : "Tracks" + }, + "pos_Z0_h" : { + "bins" : 200, + "minX" : -5, + "maxX" : 5, + "xtitle" : "z_{0} [mm]", + "ytitle" : "Tracks" + }, + "pos_time_h" : { + "bins" : 200, + "minX" : -10, + "maxX" : 10, + "xtitle" : "track time [ns]", + "ytitle" : "Tracks" + }, + "pos_chi2ndf_h" : { + "bins" : 200, + "minX" : 0, + "maxX" : 30, + "xtitle" : "pos track #chi^{2} / ndf", + "ytitle" : "Tracks" + }, + "pos_chi2_h" : { + "bins" : 200, + "minX" : 0, + "maxX" : 20, + "xtitle" : "track #chi^{2}", + "ytitle" : "Tracks" + }, + "pos_p_h" : { + "bins" : 250, + "minX" : 0, + "maxX" : 2.5, + "xtitle" : "p_{e^{+}} [GeV]", + "ytitle" : "Tracks" + }, + "vtx_chi2_h" : { + "bins" : 200, + "minX" : 0, + "maxX" : 30, + "xtitle" : "vtx #chi^{2}", + "ytitle" : "Vertices" + }, + "vtx_X_svt_h" : { + "bins" : 200, + "minX" : -5, + "maxX" : 5, + "xtitle" : "vtx X pos [mm]", + "ytitle" : "Vertices" + }, + "vtx_Y_svt_h" : { + "bins" : 200, + "minX" : -5, + "maxX" : 5, + "xtitle" : "vtx Y pos [mm]", + "ytitle" : "Vertices" + }, + "vtx_Z_svt_h" : { + "bins" : 200, + "minX" : -50, + "maxX" : 50, + "xtitle" : "vtx Z pos [mm]", + "ytitle" : "Vertices" + }, + "vtx_sigma_X_h" : { + "bins" : 100, + "minX" : 0, + "maxX" : 5, + "xtitle" : "vtx #sigma_{x} [mm]", + "ytitle" : "Vertices" + }, + "vtx_sigma_Y_h" : { + "bins" : 100, + "minX" : 0, + "maxX" : 5, + "xtitle" : "vtx #sigma_{y} [mm]", + "ytitle" : "Vertices" + }, + "vtx_sigma_Z_h" : { + "bins" : 100, + "minX" : 0, + "maxX" : 5, + "xtitle" : "vtx #sigma_{z} [mm]", + "ytitle" : "Vertices" + }, + "vtx_InvM_h" : { + "bins" : 200, + "minX" : 0, + "maxX" : 0.2, + "xtitle" : "vtx Mass [GeV]", + "ytitle" : "Vertices" + }, + "vtx_InvMErr_Z_h" : { + "bins" : 100, + "minX" : 0, + "maxX" : 0.05, + "xtitle" : "vtx #sigma_{z} [mm]", + "ytitle" : "Vertices" + }, + "vtx_px_h" : { + "bins" : 300, + "minX" : -1.5, + "maxX" : 1.5, + "xtitle" : "vtx p_{x} [GeV]", + "ytitle" : "Vertices" + }, + "vtx_py_h" : { + "bins" : 300, + "minX" : -1.5, + "maxX" : 1.5, + "xtitle" : "vtx p_{y} [GeV]", + "ytitle" : "Vertices" + }, + "vtx_pz_h" : { + "bins" : 350, + "minX" : 0, + "maxX" : 3.5, + "xtitle" : "vtx p_{z} [GeV]", + "ytitle" : "Vertices" + }, + "vtx_p_h" : { + "bins" : 350, + "minX" : 0, + "maxX" : 3.5, + "xtitle" : "vtx p [GeV]", + "ytitle" : "Vertices" + }, + "vtx_Psum_h" : { + "bins" : 350, + "minX" : 0, + "maxX" : 3.5, + "xtitle" : "vtx p_{sum} [GeV]", + "ytitle" : "Vertices" + }, + "vtx_Esum_h" : { + "bins" : 350, + "minX" : 0, + "maxX" : 3.5, + "xtitle" : "vtx E_{sum} [GeV]", + "ytitle" : "Vertices" + }, + "Esum_h" : { + "bins" : 350, + "minX" : 0, + "maxX" : 3.5, + "xtitle" : "E_{e^{-}} + E_{e^{+}} [GeV]", + "ytitle" : "Events" + }, + "Psum_h" : { + "bins" : 350, + "minX" : 0, + "maxX" : 3.5, + "xtitle" : "p_{e^{-}} + p_{e^{+}} [GeV]", + "ytitle" : "Events" + }, + "hitCode_h" : { + "bins" : 16, + "minX" : -0.5, + "maxX" : 15.5, + "xtitle" : "hit code", + "ytitle" : "Tracks" + }, + "vtx_proj_significance_h": { + "bins" : 200, + "minX" : 0.0, + "maxX" : 10.0, + "xtitle" : "vtx proj significance N#sigma", + "ytitle" : "Events" + } +} diff --git a/analysis/selections/simps/Tight_2016_simp_SR_analysis.json b/analysis/selections/simps/Tight_2016_simp_SR_analysis.json new file mode 100644 index 000000000..76eee9311 --- /dev/null +++ b/analysis/selections/simps/Tight_2016_simp_SR_analysis.json @@ -0,0 +1,22 @@ +{ + "L1Requirement_eq" : { + "cut" : 1, + "id" : 0, + "info" : "L1L1" + }, + "pSum_lt" : { + "cut" : 1.9, + "id" : 1, + "info" : "P_{e^{-}} + P_{e^{+}} < 1.9 [GeV]" + }, + "pSum_gt" : { + "cut" : 0.4, + "id" : 2, + "info" : "P_{e^{-}} + P_{e^{+}} > 0.4 GeV" + }, + "nVtxs_eq" : { + "cut" : 1, + "id" : 3, + "info" : "N_{vtx}=1" + } +} diff --git a/analysis/selections/simps/radMatchTight_2016_simp_SR_analysis.json b/analysis/selections/simps/radMatchTight_2016_simp_SR_analysis.json new file mode 100644 index 000000000..fca0a9275 --- /dev/null +++ b/analysis/selections/simps/radMatchTight_2016_simp_SR_analysis.json @@ -0,0 +1,27 @@ +{ + "L1Requirement_eq" : { + "cut" : 1, + "id" : 0, + "info" : "L1L1" + }, + "pSum_lt" : { + "cut" : 1.9, + "id" : 1, + "info" : "P_{e^{-}} + P_{e^{+}} < 1.9 [GeV]" + }, + "pSum_gt" : { + "cut" : 0.4, + "id" : 2, + "info" : "P_{e^{-}} + P_{e^{+}} > 0.4 GeV" + }, + "isRadEle_eq" : { + "cut" : 1, + "id" : 3, + "info" : "isRadEle" + }, + "nVtxs_eq" : { + "cut" : 1, + "id" : 4, + "info" : "N_{vtx}=1" + } +} diff --git a/analysis/selections/simps/radMatchTight_L1L1_nvtx1.json b/analysis/selections/simps/radMatchTight_L1L1_nvtx1.json new file mode 100644 index 000000000..add042862 --- /dev/null +++ b/analysis/selections/simps/radMatchTight_L1L1_nvtx1.json @@ -0,0 +1,17 @@ +{ + "L1Requirement_eq": { + "cut": 1, + "id": 0, + "info": "L1L1" + }, + "isRadEle_eq" : { + "cut" : 1, + "id" : 1, + "info" : "isRadEle" + }, + "nVtxs_eq": { + "cut": 1, + "id": 2, + "info": "N_{vtx}=1" + } +} diff --git a/analysis/selections/simps/radMatchTight_nocuts.json b/analysis/selections/simps/radMatchTight_nocuts.json new file mode 100644 index 000000000..2989e5ddd --- /dev/null +++ b/analysis/selections/simps/radMatchTight_nocuts.json @@ -0,0 +1,7 @@ +{ + "isRadEle_eq" : { + "cut" : 1, + "id" : 0, + "info" : "isRadEle" + } +} diff --git a/analysis/selections/trackHit/trackHitAna.json b/analysis/selections/trackHit/trackHitAna.json index 9f8ba19cb..83cd1b868 100644 --- a/analysis/selections/trackHit/trackHitAna.json +++ b/analysis/selections/trackHit/trackHitAna.json @@ -5,18 +5,38 @@ "info" : "N Hits >= 10" }, "chi2ndf_lt" : { - "cut" : 10.0, + "cut" : 6.0, "id" : 1, - "info" : "#chi^{2}/ndf <= 10.0" + "info" : "#chi^{2}/ndf <= 6.0" }, "p_gt" : { - "cut" : 0.4, + "cut" : 1.0, "id" : 2, - "info" : "p > 0.4" + "info" : "p > 1.0" }, "p_lt" : { - "cut" : 6.0, + "cut" : 4.0, "id" : 3, - "info" : "p < 6.0" + "info" : "p < 4.0" + }, + "trk_ecal_lt" : { + "cut" : 50.0, + "id" : 4, + "info" : "trk ecal x < 50 mm" + }, + "trk_ecal_gt" : { + "cut" : -100.0, + "id" : 5, + "info" : "trk ecal x > -100 mm" + }, + "trk_time_gt" : { + "cut" : -20.0, + "id" : 6, + "info" : "trk t > -20 ns" + }, + "trk_time_lt" : { + "cut" : 20.0, + "id" : 7, + "info" : "trk t < 20 ns" } } diff --git a/analysis/src/AnaHelpers.cxx b/analysis/src/AnaHelpers.cxx index 74f6dee13..07f954acc 100644 --- a/analysis/src/AnaHelpers.cxx +++ b/analysis/src/AnaHelpers.cxx @@ -87,6 +87,38 @@ bool AnaHelpers::MatchToGBLTracks(int ele_id, int pos_id, Track* & ele_trk, Trac return foundele * foundpos; } +//Use this to get two electrons +bool AnaHelpers::GetSameParticlesFromVtx(Vertex* vtx, Particle*& ele1, Particle*& ele2) { + + + bool foundele1 = false; + bool foundele2 = false; + + for (int ipart = 0; ipart < vtx->getParticles().GetEntries(); ++ipart) { + + + int pdg_id = ((Particle*)vtx->getParticles().At(ipart))->getPDG(); + if (debug_) std::cout<<"In Loop "<getParticles().At(ipart)); + foundele1=true; + } + else{ + ele2 = ((Particle*)vtx->getParticles().At(ipart)); + foundele2=true; + } + } + } + + if (!ele1 || !ele2) { + std::cout<<"Vertex formed without ele/ele. Skip."< data', color='darkblue') +plt.scatter(runs, ypositions, marker='o', s=150,label='Beamspot data', color='teal') +plt.axhline(mc_unrotated_mean_x, linestyle='--', linewidth=5.0, color='darkblue', label ='Beamspot MC') +plt.axhline(mc_unrotated_mean_y, linestyle='--', linewidth=5.0, color='teal', label ='Beamspot MC') +plt.xlabel('Run Number') +plt.ylabel('Beamspot Position [mm]') +plt.legend(fontsize=40) +plt.savefig(f'{outdir}/data_fitted_beamspot_positions_unrotated.pdf') + +# Plot data and MC beamspot x and y widths +fig, ax = plt.subplots(figsize=(30,20)) +plt.scatter(runs, xwidths, marker='^', s=150,label='Beamspot $\sigma_{x}$ data', color='darkblue') +plt.scatter(runs, ywidths, marker='^', s=150,label='Beamspot $\sigma_{y}$ data', color='teal') +plt.axhline(mc_unrotated_sigma_x, linestyle='--', linewidth=5.0, color='darkblue', label ='Beamspot $\sigma_{x}$ MC') +plt.axhline(mc_unrotated_sigma_y, linestyle='--', linewidth=5.0, color='teal', label ='Beamspot $\sigma_{x}$ MC') +plt.xlabel('Run Number') +plt.ylabel('Beamspot Width [um]') +plt.legend(fontsize=40) +plt.savefig(f'{outdir}/data_fitted_beamspot_widths_unrotated.pdf') + +# Plot data and MC beamspot rotation angles +fig, ax = plt.subplots(figsize=(30,20)) +plt.scatter(runs, angles, marker='*', s=250,label='Beamspot Rotation Angle Data', color='darkred') +plt.axhline(mc_rotation_angle_mrad, linestyle='--', linewidth=5.0, color='darkblue', label ='Beamspot Rotation Angle MC') +plt.xlabel('Run Number') +plt.ylabel('Beam Rotation Angle [mrad]') +plt.legend(fontsize=40) +plt.savefig(f'{outdir}/data_fitted_beamspot_rotation_angles.pdf') + diff --git a/plotUtils/simps/gen_oim_lookuptable.py b/plotUtils/simps/gen_oim_lookuptable.py new file mode 100644 index 000000000..ece546641 --- /dev/null +++ b/plotUtils/simps/gen_oim_lookuptable.py @@ -0,0 +1,342 @@ +#!/usr/bin/env python3 +# coding: utf-8 + +# # $\bar{C}_\text{max}(1 - \alpha, \mu)$ + +# Computing $\bar{C}_\text{max}(C, \mu)$ for optimum interval calculation, where $\mu$ is the number of expected events and $1 - \alpha$ is how frequently you reject the null hypothesis when it is true. + +# The single-event energy spectrum, that is, the probability density function which tells us which energy depositions are likely to occur, is independent of the chosen WIMP model -- we always expect a simple exponential recoil spectrum. +# +# The number of dark matter events detected does depend on the WIMP mass and cross-section. We know, however, that it must follow a Poisson distribution, which leaves the Poisson mean (which equals the expected number of events) as the only parameter left to estimate. From an upper limit on this mean, an upper limit curve in the dark matter mass – cross-section plane can be computed. +# +# * A list_of_energies list of reconstructed energy depositions of single events (from here on simply ‘energies’), either measured during some run of an actual detector, or generated using Monte Carlo.) +# * An interval is an interval in energy space. +# * The size of an interval is the fraction of energies expected in that interval. Clearly, this depends on which energy spectrum we assume, but is independent of the Poisson mean we are trying to constrain. By definition this is a number between 0 and 1. +# * The K-largest interval of a run is the largest interval containing K events in that run. Recall our definition of size: a ‘large’ interval is one which is unusually empty in that run. Clearly k-largest intervals will terminate at (or technically, just before) an observed energy, or at one of the boundaries of our energy space. Again, which interval in a run is the k–largest, depends on our energy spectrum, but not on our Poisson mean. +# * The extremeness of a K-largest interval is the probability of finding the K-largest interval in a run to be smaller. This clearly does depend on the Poisson mean: if we expect very few events, large gap sizes are more likely. Clearly extremeness is a number between 0 and 1; values close to 1 indicate unusually large intervals, that is, usually large (almost-)empty regions in the measured energies. +# For example, if the extremeness of a k-largest interval in a run is 0.8, that means that 80% of runs have k-largest intervals which are smaller than the k-largest interval in this run. +# * The optimum interval statistic of a run is extremity of the most extreme k-largest interval in a run. +# * The extremeness of the optimum interval statistic is the probability of finding a lower optimum interval statistic, that is, of finding the optimum interval in a run to be less extreme. +# +# The max gap method rejects a theory (places a mean outside the upper limit) based on a run if the 0-largest interval (the largest gap) is too extreme. +# +# The optimum interval method rejects a theory based on a run if the optimum interval statistic is too large. +# +# * The energy cumulant $\epsilon(E)$ is the fraction of energies expected below the energy $E$. Whatever the (1-normalized) energy distribution $dN/dE$, $dN/d\epsilon$ is uniform[0,1], where $0$ and $1$ correspond to the boundaries of our experimental range. +# +# + +import functools +from scipy.optimize import brenth +import matplotlib.pyplot as plt +import numpy as np +import pickle +import os +import sys +import time + + +# Disable +def blockPrint(): + sys.stdout = open(os.devnull, 'w') + +# Restore +def enablePrint(): + sys.stdout = sys.__stdout__ + +def kLargestIntervals(list_of_energies, spectrumCDF = lambda x: x): + """ + Returns a list of the sizes of the K-largest intervals in that run according to the energy spectrum (given as a CDF). + That is, kLargestIntervals(...)[i] is the size of the largest interval containing i events, where ‘largest’ is defined above. + + * Transform energies to energy cumulants + * Add events at 0 and 1 + * Foreach k, compute interval sizes, take max + """ + answer = {} + + list_of_energies.sort() + #print('list of energies', list_of_energies) + + energy_cumulants = spectrumCDF(list_of_energies) + #print('cumulants', energy_cumulants) + + for i,interval_size in enumerate(range(len(energy_cumulants))): + if (1 + interval_size) >= len(energy_cumulants): + continue + + temp_data = energy_cumulants.copy() + gap_sizes = temp_data[(1+interval_size):] - temp_data[0:-1*(1 + interval_size)] + #print(f'gap sizes in interval_size {i+1}:', gap_sizes) + + answer[interval_size] = np.max(gap_sizes) + return answer + +def extremenessOfInterval(x, k, mu): + """ + Returns the extremeness of a k-largest interval of size, if the poisson mean is mu. + + (Number of itvSizes[mu][k] smaller than size) / mcTrials[mu] + + x - also size in above comment + k - gap (rename k) + """ + #print('extremenessOfInterval') + #print(f'mu is {mu}') + # [0] is because where returns list, where [0] is answer + if k not in itvSizes[mu]: + return 0 + ##print('extremeness of x:', x) + #print(f'mu is {mu}') + #print(f'k is {k}') + #print('mcTrials[mu]:', mcTrials[mu]) + #print(itvSizes[mu][k]) + #print(np.where(itvSizes[mu][k] < x)[0]) + #print(f'number of gaps with size less than {x}:', np.where(itvSizes[mu][k] < x)[0].size) + #print(f'fraction of gaps less than {x}:', np.where(itvSizes[mu][k] < x)[0].size / mcTrials[mu]) + return np.where(itvSizes[mu][k] < x)[0].size / mcTrials[mu] + + +def optimumItvStatistic(list_of_energies, mu, spectrumCDF = lambda x: x): + """ + Returns the optimum interval statistic of the run. + + Max of extremenssOfInterval's + """ + #print('running optimumItvStatistic') + #print(f'mu is {mu}') + #print(f'list of energies: {list_of_energies}') + #print('check klargest intervals:', kLargestIntervals(list_of_energies, spectrumCDF).items()) + #print('max is', np.max([extremenessOfInterval(x, k, mu) for k, x in kLargestIntervals(list_of_energies, spectrumCDF).items()])) + return np.max([extremenessOfInterval(x, k, mu) for k, x in kLargestIntervals(list_of_energies, spectrumCDF).items()]) + + +def extremenessOfOptItvStat(stat, mu): + """ + Returns the extremeness of the optimum interval statistic stat, given mu + + (Number of optItvs[mu] smaller than stat) / mcTrials[mu] + """ + #print('extremenessOfOptItvStat') + #print(f'mu is {mu}') + return np.where(optItvs[mu] < stat)[0].size / mcTrials[mu] + +def optItvUpperLimit(list_of_energies, c, spectrumCDF = lambda x: x, + n = 1000): + """ + Returns the c- confidence upper limit on mu using optimum interval + + For which mu is extremenessOfOptItvStat( optimumItvStatistic(run), mu ) = c + + c - e.g., 0.9 + """ + def f(mu, list_of_energies, c, spectrumCDF, n): + generate_table(mu, n) + x = optimumItvStatistic(list_of_energies, mu, spectrumCDF) + prob = extremenessOfOptItvStat(x, mu) + return prob - c + + mu = 0 + + for mu in np.arange(10, 2 * list_of_energies.size): + if f(mu, list_of_energies, c, spectrumCDF, n) > 0: + #print('Found seed mu=%f' % mu) + break + + try: + xsec = brenth(f, mu - 5, mu + 5, + args=(list_of_energies, c, spectrumCDF, n), + xtol=1e-2) + #print('Improved xsec:', xsec) + except: + #print("ERROR: could not minimize", mu) + return mu + return xsec + +def generate_trial_experiment(mu, n): + #print(f'Generate Trial Experiment with mu:{mu} and n{n}') + trials = [] + + for index in range(n): + this_mu = np.random.poisson(mu) + #print(f'this_mu: {this_mu}') + rand_numbers = np.random.random(size=this_mu) + #print(f'rand numbers, {rand_numbers}') + rand_numbers = np.append(rand_numbers, [0.0, 1.0]) + rand_numbers.sort() + #print(f'rand numbers, {rand_numbers}') + trials.append(rand_numbers) + + return trials + + +# ## Monte Carlo for populating itvSizes[$\mu$][$k$] and optItvs[$\mu$] + + +def get_filename(): + #return '/fs/ddn/sdf/group/hps/users/alspellm/mc_storage/opt_int_lookuptable_max50_10ktoys.p' + #return '/fs/ddn/sdf/group/hps/users/alspellm/mc_storage/opt_int_lookuptable_max25_100ktoys.p' + #return '/fs/ddn/sdf/group/hps/users/alspellm/mc_storage/opt_int_lookuptable_max25_100ktoys_0.05steps.p' + return '/fs/ddn/sdf/group/hps/users/alspellm/mc_storage/opt_int_lookuptable_max25_10ktoys_0.05steps_v2.p' + +def load_table_from_disk(): + global itvSizes + global optItvs + global mcTrials + + if os.path.exists(get_filename()): + f = open(get_filename(), 'rb') + itvSizes = pickle.load(f) + optItvs = pickle.load(f) + mcTrials = pickle.load(f) + f.close() + + +def write_table_to_disk(): + f = open(get_filename(), 'wb') + pickle.dump(itvSizes, f) + pickle.dump(optItvs, f) + pickle.dump(mcTrials, f) + f.close() + +itvSizes = {} +optItvs = {} +mcTrials = {} +load_table_from_disk() + +def generate_table_new(mu,n): + """ #Generate trial runs""" + if mu in mcTrials and mcTrials[mu] >= n: + return + + mcTrials[mu] = n + itvSizes[mu] = {} + optItvs[mu] = [] + + mu_trials = np.random.poisson(mu, size=n) + trials = [np.sort(np.append(np.random.random(mu),[0.0,1.0])) for mu in mu_trials] + + for i,trial in enumerate(trials): + + intermediate_result = kLargestIntervals(trial) + for k, v in intermediate_result.items(): + if k not in itvSizes[mu]: + itvSizes[mu][k] = [] + + itvSizes[mu][k].append(v) + #print('interm results') + #print(itvSizes) + # Numpy-ize it + for k, array in itvSizes[mu].items(): + itvSizes[mu][k] = np.array(array) + + for i,trial in enumerate(trials): + #print('\n') + #print(f'trial {i}: {trial}') + #print(f'mu is {mu}') + optItvs[mu].append(optimumItvStatistic(trial, mu)) + #print('trial result: ', optItvs[mu]) + #print(f'trial result: {optimumItvStatistic(trial, mu)}') + #print('\n') + + #print('summarize trials') + + # Numpy-ize it + optItvs[mu] = np.array(optItvs[mu]) + #print(optItvs[mu]) + +def generate_table(mu, n): + """ #Generate trial runs""" + if mu in mcTrials and mcTrials[mu] >= n: + return + + #print("Generating mu=", mu) + + mcTrials[mu] = n + trials = generate_trial_experiment(mu, mcTrials[mu]) + #print(trials[0:10]) + itvSizes[mu] = {} + optItvs[mu] = [] + + for i,trial in enumerate(trials): + + intermediate_result = kLargestIntervals(trial) + for k, v in intermediate_result.items(): + if k not in itvSizes[mu]: + itvSizes[mu][k] = [] + + itvSizes[mu][k].append(v) + #print('interm results') + #print(itvSizes) + # Numpy-ize it + for k, array in itvSizes[mu].items(): + itvSizes[mu][k] = np.array(array) + + for i,trial in enumerate(trials): + #print('\n') + #print(f'trial {i}: {trial}') + #print(f'mu is {mu}') + optItvs[mu].append(optimumItvStatistic(trial, mu)) + #print('trial result: ', optItvs[mu]) + #print(f'trial result: {optimumItvStatistic(trial, mu)}') + #print('\n') + + #print('summarize trials') + + # Numpy-ize it + optItvs[mu] = np.array(optItvs[mu]) + #print(optItvs[mu]) + + + +def cache_values(my_max=200, n=100): + for i in range(3, my_max): + generate_table(i, n) + write_table_to_disk() + + +import time +def plot_something(): + x, y = [], [] + + for i,mu in enumerate(np.linspace(0.0, 200.0,1000)): + start_time = time.time() + generate_table(mu, 10000) + x.append(mu) + a = brenth(lambda x: extremenessOfOptItvStat(x, mu) - 0.9, + 0, + 1, + xtol=1e-2) + + y.append(a) + + end_time = time.time() + print('time: ', end_time - start_time) + + plt.plot(x,y) + plt.xscale('log') + plt.xlim(0.0, 200.0) + + +#plot_something() +total_time = 0.0 +for i,mu in enumerate(np.linspace(0.0, 25.0,500)): + print(f'Running mu = {mu} | {100.*i/500}% complete') + start_time = time.time() + generate_table(mu, 10000) + #generate_table_new(mu, 100000) + end_time = time.time() + elapsed_time = end_time - start_time + total_time = total_time + elapsed_time + avg_time = total_time/(i+1) + print(f'Time to generate table entry: {elapsed_time}') + print(f'Average time: {avg_time}') + print(f'Estimated completion: {avg_time*(500-i+1)}') + +write_table_to_disk() + + + + + + diff --git a/plotUtils/simps/mass_resolution/.moller_ana.py.swp b/plotUtils/simps/mass_resolution/.moller_ana.py.swp new file mode 100644 index 000000000..5b60fb8b7 Binary files /dev/null and b/plotUtils/simps/mass_resolution/.moller_ana.py.swp differ diff --git a/plotUtils/simps/mass_resolution/fee_smearing_nhits.py b/plotUtils/simps/mass_resolution/fee_smearing_nhits.py new file mode 100644 index 000000000..681b0c088 --- /dev/null +++ b/plotUtils/simps/mass_resolution/fee_smearing_nhits.py @@ -0,0 +1,189 @@ +#!/usr/bin/python3 +""" +This script calculates the FEE momentum smearing required to match MC and data. +FEE's are selected in MC and data, and their vertex Psum FEE peak is fit with a Gaussian function to measure +the FEE momentum resolution. +A smearing factor is calculated as a function of the number of hits on track to match MC to data. +""" +import os +import numpy as np +import math +import ROOT as r + +#======================================================================================================================================= +# INITIALIZE +#======================================================================================================================================= +# --outdir: Specify output directory. + +import argparse +parser = argparse.ArgumentParser(description='Process some inputs.') +parser.add_argument('--outdir', type=str, default='moller_mass_fits') +args = parser.parse_args() +outdir = args.outdir + +#======================================================================================================================================= +# FUNCTIONS +#======================================================================================================================================= +def gaus_fit(histo, xmin, xmax, smean, swidth, snorm, nsigma=2.0, isData=False): + + #initial fit with seeds + fitfunc = r.TF1("gaus","gaus") + fitfunc.SetParameter(0, snorm) + fitfunc.SetParameter(1, smean) + fitfunc.SetParameter(2, swidth) + fitRes = histo.Fit(fitfunc,"QLES","", xmin, xmax) + params = fitRes.Parameters() + chi2 = fitRes.Chi2() + ndf = fitRes.Ndf() + + #set first fist to be best fit + best_chi2 = chi2/ndf + best_params = params + + #iterate over randomly fluctuated fit parameters. Keep the best resulting fit + niters = 100 + for n in range(niters): + norm = params[0]*np.random.uniform(80,120)*0.01 + mu = params[1]*np.random.uniform(80,120)*0.01 + sigma = params[2]*np.random.uniform(80,120)*0.01 + + #Data has shoulders, so we can specify the xmin and xmax to do an asymmetric fit window + if isData: + xminx = mu - nsigma*sigma + xmaxx = mu + nsigma*sigma + if xminx < xmin: + xminx = xmin + if xmaxx > xmax: + xmaxx = xmax + fitfunc.SetParameter(0, norm) + fitfunc.SetParameter(1, mu) + fitfunc.SetParameter(2, sigma) + fitRes = histo.Fit(fitfunc,"QLES","", xminx, xmaxx) + #If fit fails, skip + try: + if fitRes.Parameters()[1] < xmin or fitRes.Parameters()[1] > xmax or fitRes.Ndf() < 1: + continue + except: + continue + + params = fitRes.Parameters() #these results seed the next iteration...maybe should only do if improved? + chi2 = fitRes.Chi2() + ndf = fitRes.Ndf() + #replace best fit + if chi2/ndf < best_chi2: + best_params = params + + #Do the final fit using the best parameters found + fitfunc.SetParameter(0, best_params[0]) + fitfunc.SetParameter(1, best_params[1]) + fitfunc.SetParameter(2, best_params[2]) + xminx = best_params[1] - nsigma*best_params[2] + xmaxx = best_params[1] + nsigma*best_params[2] + + #again, if data, use asymmetric fit window to avoid the left shoulder + if isData: + if xminx < xmin: + xminx = xmin + if xmaxx > xmax: + xmaxx = xmax + + fitRes = histo.Fit(fitfunc,"QLES","", xminx, xmaxx) + params = fitRes.Parameters() + errors = fitRes.Errors() + chi2 = fitRes.Chi2() + ndf = fitRes.Ndf() + return histo, params, errors, chi2/ndf + +#======================================================================================================================================= +# LOAD DATA AND FIT FEE PEAK +#======================================================================================================================================= + +# Load FEEs tracks in data +data_results = {} +infilename = '/sdf/group/hps/user-data/alspellm/2016/fee_smearing/run7800/hadd/hadd_fee_2pt3_recon_fee_histos.root' #FEE skimmed tracks from hpstr track analysis processor + +# Read track hit histograms +histoname = 'KalmanFullTracks/KalmanFullTracks_p_vs_nHits_top_hh' #top +infile = r.TFile(f'{infilename}',"READ") +top_h = copy.deepcopy(infile.Get(f'{histoname}')) +histoname = 'KalmanFullTracks/KalmanFullTracks_p_vs_nHits_bot_hh' #bot +bot_h = copy.deepcopy(infile.Get(f'{histoname}')) +infile.Close() + +# Change the names to use as keys +top_h.SetName('top') +bot_h.SetName('bot') + +# Fit the FEE peak for each category of nhits. Just have access to 10, 11, 12 for now +for h in [top_h, bot_h]: + histo = h + for nhits in [10, 11, 12]: + # Get the nhits momentum projection + proj = histo.ProjectionY(f'proj_{h.GetName()}_{nhits}hits', histo.GetXaxis().FindBin(nhits), histo.GetXaxis().FindBin(nhits),"") + # Fit the data + _, params, errors, chi2ndf = gaus_fit(proj, 2.0, 2.5, 2.4, 0.47, 12000, nsigma=1.5, isData=True) + + # Store the results [mu,sigma] for top/bot nhits= + data_results[f'{h.GetName()}_nhits_{nhits}'] = [params[1], params[2]] + +#======================================================================================================================================= +# LOAD MC AND FIT FEE PEAK +#======================================================================================================================================= + +# Load MC FEE's from hpstr track analysis processor +mc_results = {} +infilename= '/sdf/group/hps/user-data/alspellm/2016/fee_smearing/tritrig/hadd/hadd_fee_2pt3_recon_tritrig_histos.root' + +# Read track hit histograms +histoname = 'KalmanFullTracks/KalmanFullTracks_p_vs_nHits_top_hh' #top +infile = r.TFile(f'{infilename}',"READ") +top_h = copy.deepcopy(infile.Get(f'{histoname}')) +histoname = 'KalmanFullTracks/KalmanFullTracks_p_vs_nHits_bot_hh' #bot +bot_h = copy.deepcopy(infile.Get(f'{histoname}')) +infile.Close() + +#Change the names to use as keys +top_h.SetName('top') +bot_h.SetName('bot') + +for h in [top_h, bot_h]: + histo = h + for nhits in [10, 11, 12]: + # Get the nhits momentum projection + proj = histo.ProjectionY(f'proj_{h.GetName()}_{nhits}hits', histo.GetXaxis().FindBin(nhits), histo.GetXaxis().FindBin(nhits),"") + # Fit the data + _, params, errors, chi2ndf = gaus_fit(proj, 2.1, 2.5, 2.2, 0.1, proj.GetMaximum(), nsigma=1.5) + # Store the results [mu, sigma] for top/bot nhits= + mc_results[f'{h.GetName()}_nhits_{nhits}'] = [params[1], params[2]] + +#======================================================================================================================================= +# CALCULATE MOMENTUM SMEARING FACTORS +#======================================================================================================================================= +# Store momentum smearing factors in ROOT file +outfile = r.TFile(f'{outdir}/smearingFile_2016_nhits.root',"RECREATE") +outfile.cd() +# Calculate smearing for Top and Bot +smtop_h = r.TH1F('KalmanFullTracks_p_vs_nHits_hh_smearing_rel_top','p_vs_nHits_smearing_rel_top;nhits;smear factor', 3, 9.5, 12.5) +smbot_h = r.TH1F('KalmanFullTracks_p_vs_nHits_hh_smearing_rel_bot','p_vs_nHits_smearing_rel_bot;nhits;smear factor', 3, 9.5, 12.5) + +# Calculate smearing factor according to 2016 Bump Hunt +smear_fac = lambda mu_data, sig_data, mu_mc, sig_mc : np.sqrt(np.square(sig_data/mu_data) - np.square(sig_mc/mu_mc)) +for key, vals in data_results.items(): + istop = False + if 'top' in key: + istop = True + nhits = float(key.split('_')[2]) + mu_data = vals[0] + sig_data = vals[1] + mu_mc = mc_results[key][0] + sig_mc = mc_results[key][1] + sf = smear_fac(mu_data, sig_data, mu_mc, sig_mc) + print(f'{key} sf={sf}') + + #save results + if istop: + smtop_h.SetBinContent(smtop_h.GetXaxis().FindBin(nhits), sf) + else: + smbot_h.SetBinContent(smbot_h.GetXaxis().FindBin(nhits), sf) + +outfile.Write() diff --git a/plotUtils/simps/mass_resolution/fit_mass_resolution.py b/plotUtils/simps/mass_resolution/fit_mass_resolution.py new file mode 100644 index 000000000..44fc06f75 --- /dev/null +++ b/plotUtils/simps/mass_resolution/fit_mass_resolution.py @@ -0,0 +1,342 @@ +#!/usr/bin/python3 +""" +This script fits the MC signal mass resolution as a function of invariant mass. +The mass resolution uses SIMP signal reconstructed vertex invariant mass with radMatchTight selection (truth matched ele). +Fill a histogram with invariant mass and fit the distribution with a Gaussian fit function. +Fit the mass resolution as a function of invariant mass with a polynomial. Use p-test to find best fit order. +""" +import os +import awkward as ak +import numpy as np +import hist +from hist import Hist +import uproot +import ROOT as r +import matplotlib.pyplot as plt +import matplotlib as mpl + +# SIMP tools defined in hpstr +import sys +hpstr_base = os.getenv('HPSTR_BASE') +sys.path.append(f'{hpstr_base}/plotUtils/simps') +import simp_plot_utils as utils + +plt.rcParams.update({'font.size': 40, # Font size for text + 'axes.titlesize': 40, # Font size for titles + 'axes.labelsize': 40, # Font size for axis labels + 'xtick.labelsize': 40, # Font size for x-axis tick labels + 'ytick.labelsize': 40, # Font size for y-axis tick labels + 'lines.linewidth':4.0, + 'legend.fontsize': 40}) # Font size for legend + +#======================================================================================================================================= +# INITIALIZE +#======================================================================================================================================= +# --outdir: Specify output directory. + +import argparse +parser = argparse.ArgumentParser(description='Process some inputs.') +parser.add_argument('--outdir', type=str, default='moller_mass_fits') +parser.add_argument('--outfilename', type=str, default='mass_resolution') +args = parser.parse_args() +outdir = args.outdir +outfilename = args.outfilename + +#======================================================================================================================================= +# FUNCTIONS +#======================================================================================================================================= + +def load_signal(filepath, selection, cuts=None, expressions=None): + #load signal with the hpstr ana processor flat tuple struct + with uproot.open(filepath) as f: + events = f[f'{selection}/{selection}_tree'].arrays( + expressions = expressions, cut = cuts + ) + return events + +def gaus_fit(histo, xmin, xmax, smean, swidth, snorm, nsigma=2.0, isData=False): + + print('seeds:', xmin, xmax, smean, swidth, snorm) + #initial fit with seeds + fitfunc = r.TF1("gaus","gaus") + fitfunc.SetParameter(0, snorm) + fitfunc.SetParameter(1, smean) + fitfunc.SetParameter(2, swidth) + fitRes = histo.Fit(fitfunc,"QLES","", xmin, xmax) + try: + params = fitRes.Parameters() + chi2 = fitRes.Chi2() + ndf = fitRes.Ndf() + + best_chi2 = chi2/ndf + best_params = params + except: + best_chi2 = 999999.9 + params = [snorm, smean, swidth] + best_params=params + + niters = 100 + for n in range(niters): + norm = params[0]#*np.random.uniform(80,120)*0.01 + mu = params[1]#*np.random.uniform(80,120)*0.01 + sigma = params[2]#*np.random.uniform(80,120)*0.01 + + xminx = mu - nsigma*sigma + xmaxx = mu + nsigma*sigma + if isData: + if xminx < xmin: + xminx = xmin + if xmaxx > xmax: + xmaxx = xmax + fitfunc.SetParameter(0, norm) + fitfunc.SetParameter(1, mu) + fitfunc.SetParameter(2, sigma) + fitRes = histo.Fit(fitfunc,"QLES","", xminx, xmaxx) + try: + if fitRes.Parameters()[1] < xminx or fitRes.Parameters()[1] > xmaxx or fitRes.Ndf() < 1: + continue + except: + continue + + params = fitRes.Parameters() + #print(params) + chi2 = fitRes.Chi2() + ndf = fitRes.Ndf() + if chi2/ndf < best_chi2: + best_params = params + + fitfunc.SetParameter(0, best_params[0]) + fitfunc.SetParameter(1, best_params[1]) + fitfunc.SetParameter(2, best_params[2]) + xminx = best_params[1] - nsigma*best_params[2] + xmaxx = best_params[1] + nsigma*best_params[2] + + if isData: + if xminx < xmin: + xminx = xmin + if xmaxx > xmax: + xmaxx = xmax + + print('result: ', xminx, xmaxx, best_params[1], best_params[2], best_params[0]) + fitRes = histo.Fit(fitfunc,"QLES","", xminx, xmaxx) + params = fitRes.Parameters() + errors = fitRes.Errors() + chi2 = fitRes.Chi2() + ndf = fitRes.Ndf() + if ndf > 0: + chi2 = chi2/ndf + else: + chi2 = 99999.9 + return histo, params, errors, chi2 + +def fit_with_poly(tgrapherrs): + polys = [] + chi2s = [] + fstats = [] + fitResults = [] + npoints = len(tgrapherrs.GetX()) + + for n in range(10): + fitfunc = r.TF1(f'pol{n}',f'pol{n}', tgrapherrs.GetX()[0], tgrapherrs.GetX()[-1]) + fitfunc.SetRange(30.0, 124) + fitfunc.SetLineColor(r.kRed) + fitfunc.SetMarkerSize(0.0) + fitRes = tgrapherrs.Fit(fitfunc,"SRQ") + chi2s.append(fitRes.Chi2()) + polys.append(n) + fitResults.append(fitRes) + + #Perform fstat test to see how much fit improves with additional order (why does this work?) + if n > 0: + fstats.append( (chi2s[n-1]-chi2s[n])*(npoints-n-1)/(chi2s[n])) + else: + fstats.append(0.0) + + #Pick the order that shows greatest positive improvement in fstat + best_diff = 0.0 + best_n = None + for n,fstat in enumerate(fstats): + if n == 0: + continue + diff = fstats[n-1] - fstat + if diff > 0 and diff > best_diff: + best_diff = diff + best_n = n + + print(f'best n: {best_n}') + return fitResults[best_n] + +#======================================================================================================================================= +# LOAD MC SIGNAL +#======================================================================================================================================= + +# Load MC signal with momentum smearing +indir = '/sdf/group/hps/user-data/alspellm/2016/simp_mc/pass4b/beam/smeared' +infilename = lambda mass: f'mass_{mass}_hadd-simp-beam_ana_smeared_corr.root' +# Masses of each MC generated signal mass +masses = [x for x in range(30,126,2)] + +# Specify any cuts. Here we put a cut on the number of hits on track based on the 2016 SIMP L1L1 analysis +inv_masses_h={} +cuts = '( (unc_vtx_psum > 1.0) & (unc_vtx_ele_track_nhits >= 7) & (unc_vtx_pos_track_nhits >= 7) )' +# Load each MC signal mass +for mass in masses: + infile = os.path.join(indir, infilename(mass)) + selection = 'vtxana_radMatchTight_2016_simp_SR_analysis' # Use radMatchTight + branches = ['unc_vtx_mass', 'unc_vtx_ele_track_nhits', 'unc_vtx_pos_track_nhits'] + # Load events from hpstr vertex analysis processor + tree = load_signal(infile,selection, expressions=branches, cuts=cuts) + + # Initialize histogram + inv_masses_h[f'mass_{mass}'] = ( + hist.Hist.new + .Reg(800,17.5,217.5,label='Invariant Mass [MeV]') + .Double() + ) + + # Fill histogram + inv_masses_h[f'mass_{mass}'].fill(tree.unc_vtx_mass*1000.) + +# Load MC signal WITHOUT momentum smearing +indir = '/sdf/group/hps/user-data/alspellm/2016/simp_mc/pass4b/beam/nosmearing' +infilename = lambda mass: f'mass_{mass}_hadd-simp-beam-ana-nosmearing.root' +inv_masses_unsm_h={} +for mass in masses: + infile = os.path.join(indir, infilename(mass)) + selection = 'vtxana_radMatchTight_2016_simp_SR_analysis' + branches = ['unc_vtx_mass' ] + tree = load_signal(infile,selection, expressions=branches, cuts=cuts) + + inv_masses_unsm_h[f'mass_{mass}'] = ( + hist.Hist.new + .Reg(800,17.5,217.5,label='Invariant Mass [MeV]') + .Double() + ) + + inv_masses_unsm_h[f'mass_{mass}'].fill(tree.unc_vtx_mass*1000.) + + +# Convert Histograms to ROOT for convenient fitting +for key, histo in inv_masses_h.items(): + h = utils.cnvHistogramToROOT(histo) + h.SetName(f'{key}_smeared') + inv_masses_h[key] = h + +for key, histo in inv_masses_unsm_h.items(): + h = utils.cnvHistogramToROOT(histo) + h.SetName(f'{key}_unsmeared') + inv_masses_unsm_h[key] = h + + +# Initialize output ROOT file to save fit results +outfile = r.TFile(f'{outdir}/{outfilename}',"RECREATE") +outfile.cd() +fit_results = {} + +#======================================================================================================================================= +# CALCULATE MC SIGNAL MASS RESOLUTION WITH SMEARING +#======================================================================================================================================= + +# Run over smeared MC signal and fit +for mass, histo in inv_masses_h.items(): + print(f'Fitting {mass}') + fit_histo, params, errors, chi2ndf = gaus_fit(histo, histo.GetXaxis().GetBinLowEdge(histo.FindFirstBinAbove(0)), histo.GetXaxis().GetBinLowEdge(histo.FindLastBinAbove(0)), histo.GetMean(), histo.GetRMS(), histo.GetMaximum(), nsigma=2.0, isData=False) + + #Draw Fit result to canvas, and save to png and to root file + text = [f'\mu = {round(params[1],2)}',f'\sigma = {round(params[2],2)}',f'\chi^{2}/n.d.f = {round(chi2ndf,2)}'] + xmin = params[1] - 4.0*params[2] + xmax = params[1] + 4.0*params[2] + ymax = params[0]*1.1 + c = utils.drawTH1s([histo], histo.GetName(), drawOpts=['hist'],xrange=[xmin,xmax], yrange=[0.0,ymax],size=(2000,1500), text=text, text_pos=[0.2,0.78], line_spacing=0.05) + c.Write() + c.SaveAs(f'{outdir}/{histo.GetName()}.png') + c.Close() + #save fit results + fit_results[mass] = [fit_histo, params, errors] + +# Fit the fitted mass resolution results as a function of invariant mass +fit_masses = sorted(np.array([float(x.replace('mass_','')) for x in fit_results.keys()], dtype=float)) +fit_res = np.array([fit_results[f'mass_{int(mass)}'][1][2] for mass in fit_masses], dtype=float) +fit_errs = np.array([fit_results[f'mass_{int(mass)}'][2][2] for mass in fit_masses], dtype=float) +zeros = np.array([0.0 for mass in fit_masses], dtype=float) + +# Make TGraphErrors to fit with polynomial and get mass resolution as function of mass +massRes_smeared_ge = r.TGraphErrors(len(fit_masses), np.array(fit_masses), np.array(fit_res), np.array(zeros), np.array(fit_errs)) +massRes_smeared_ge.GetXaxis().SetTitle('Invariant Mass [MeV]') +massRes_smeared_ge.GetYaxis().SetTitle('Invariant Mass Resolution [MeV]') +massRes_smeared_ge.GetYaxis().SetTitleOffset(1.) +fitResult_smeared = fit_with_poly(massRes_smeared_ge) +text_smeared = [f'Smeared: {fitResult_smeared.Parameters()}'] +c_smeared = utils.drawTGraphs([massRes_smeared_ge], 'smeared_mc', drawOpts=['AP'], text=text_smeared, + text_pos=[0.2,0.8]) +c_smeared.SaveAs(f'{outdir}/invariant_mass_resolution_function_smeared.png') +c_smeared.Write() + +#======================================================================================================================================= +# CALCULATE MC SIGNAL MASS RESOLUTION WITH SMEARING +#======================================================================================================================================= + +# Run over unsmeared MC signal and fit +for mass, histo in inv_masses_unsm_h.items(): + print(f'Fitting {mass}') + fit_histo, params, errors, chi2ndf = gaus_fit(histo, histo.GetXaxis().GetBinLowEdge(histo.FindFirstBinAbove(0)), histo.GetXaxis().GetBinLowEdge(histo.FindLastBinAbove(0)), histo.GetMean(), histo.GetRMS(), histo.GetMaximum(), nsigma=2.0, isData=False) + + # Draw Fit result to canvas, and save to png and to root file + text = [f'\mu = {round(params[1],2)}',f'\sigma = {round(params[2],2)}',f'\chi^{2}/n.d.f = {round(chi2ndf,2)}'] + xmin = params[1] - 4.0*params[2] + xmax = params[1] + 4.0*params[2] + ymax = params[0]*1.1 + c = utils.drawTH1s([histo], f'{histo.GetName()}_unsmeared', drawOpts=['hist'],xrange=[xmin,xmax], yrange=[0.0,ymax],size=(2000,1500), text=text, text_pos=[0.2,0.78], line_spacing=0.05) + c.Write() + c.SaveAs(f'{outdir}/{histo.GetName()}.png') + c.Close() + + # Save fit results + fit_results[mass] = [fit_histo, params, errors] + +# Fit the fitted mass resolution results as a function of invariant mass +fit_masses = sorted(np.array([float(x.replace('mass_','')) for x in fit_results.keys()], dtype=float)) +fit_res = np.array([fit_results[f'mass_{int(mass)}'][1][2] for mass in fit_masses], dtype=float) +fit_errs = np.array([fit_results[f'mass_{int(mass)}'][2][2] for mass in fit_masses], dtype=float) +zeros = np.array([0.0 for mass in fit_masses], dtype=float) + +# Make TGraphErrors to fit with polynomial and get mass resolution as function of mass +massRes_unsmeared_ge = r.TGraphErrors(len(fit_masses), np.array(fit_masses), np.array(fit_res), np.array(zeros), np.array(fit_errs)) +massRes_unsmeared_ge.GetXaxis().SetTitle('Invariant Mass [MeV]') +massRes_unsmeared_ge.GetYaxis().SetTitle('Invariant Mass Resolution [MeV]') +massRes_unsmeared_ge.GetYaxis().SetTitleOffset(1.0) +fitResult_unsmeared = fit_with_poly(massRes_unsmeared_ge) +text_unsmeared = [f'unsmeared: {fitResult_unsmeared.Parameters()}'] +c_unsmeared = utils.drawTGraphs([massRes_unsmeared_ge], 'unsmeared_mc', drawOpts=['AP'], text=text_unsmeared, + text_pos=[0.2,0.8]) +print(fitResult_unsmeared.Parameters()) +c_unsmeared.SaveAs(f'{outdir}/invariant_mass_resolution_function_unsmeared.png') +c_unsmeared.Write() + + +#======================================================================================================================================= +# SUMMARY PLOT +#======================================================================================================================================= + +colors = utils.getColorsHPS() +#format smeared +utils.format_th1(massRes_smeared_ge, title='Smeared MC', linecolor=r.kBlack, markerstyle=20, markercolor=r.kBlack) +massRes_smeared_ge.SetMarkerSize(2) +massRes_smeared_ge.GetYaxis().SetTitleOffset(0.5) +massRes_smeared_ge.GetYaxis().SetRangeUser(0.0,8.0) +#format unsmeared +utils.format_th1(massRes_unsmeared_ge, title='Un-smeared MC', linecolor=r.kBlack, markerstyle=4, markercolor=r.kBlack) +massRes_unsmeared_ge.SetMarkerSize(2) +massRes_unsmeared_ge.GetListOfFunctions().At(0).SetLineColor(colors[1]) + + +c = r.TCanvas('c','c',2500,1500) +c.cd() +massRes_smeared_ge.Draw('AP') +massRes_unsmeared_ge.Draw('PSAME') +c.Draw() +legend = utils.buildLegend([massRes_smeared_ge, massRes_unsmeared_ge], position=(0.3,0.7, 0.4, 0.8), titles=['Smeared MC', 'Un-smeared MC']) +legend.SetTextSize(0.03) +legend.Draw() +c.SaveAs(f'{outdir}/invariant_mass_resolution_fits.png') + diff --git a/plotUtils/simps/mass_resolution/moller_ana.py b/plotUtils/simps/mass_resolution/moller_ana.py new file mode 100644 index 000000000..93bb8d248 --- /dev/null +++ b/plotUtils/simps/mass_resolution/moller_ana.py @@ -0,0 +1,364 @@ +#!/usr/bin/python3 +""" +This script selects Moller events in MC and data, plots the invariant mass distributions, and fits them with a Gaussian +fit function to calculate the Moller mass resolutions. +This script compares data, unsmeared MC, and smeared MC. +""" +import os +import awkward as ak +import numpy as np +import hist +from hist import Hist +import uproot +import math +import ROOT as r +import matplotlib as mpl +import matplotlib.pyplot as plt +import sys + +# SIMP tools defined in hpstr +hpstr_base = os.getenv('HPSTR_BASE') +sys.path.append(f'{hpstr_base}/plotUtils/simps') +import simp_plot_utils as utils + + +plt.rcParams.update({'font.size': 40, # Font size for text + 'axes.titlesize': 40, # Font size for titles + 'axes.labelsize': 40, # Font size for axis labels + 'xtick.labelsize': 40, # Font size for x-axis tick labels + 'ytick.labelsize': 40, # Font size for y-axis tick labels + 'lines.linewidth':4.0, + 'legend.fontsize': 40}) # Font size for legend + +#======================================================================================================================================= +# INITIALIZE +#======================================================================================================================================= +# --outdir: Specify output directory. + +import argparse +parser = argparse.ArgumentParser(description='Process some inputs.') +parser.add_argument('--outdir', type=str, default='moller_mass_fits') +args = parser.parse_args() +outdir = args.outdir + +#======================================================================================================================================= +# FUNCTIONS +#======================================================================================================================================= + +def load_data(filepath, selection, cut_expression = None, expressions=None): + with uproot.open(filepath) as f: + events = f[f'{selection}/{selection}_tree'].arrays( + cut=cut_expression, + expressions = expressions + ) + return events + +def TrackFiducial(array, isData=False, isMC=False): + """ + The Moller track fiducial regions are defined in the 2016 Bump Hunt analysis note. + The selection is slightly different between MC and data. + """ + ele_fid = ak.full_like(array.unc_vtx_ele_track_ecal_y, False, dtype=bool) + pos_fid = ak.full_like(array.unc_vtx_pos_track_ecal_y, False, dtype=bool) + + if isData: + ele_condition_1 = ( + (array.unc_vtx_ele_track_ecal_y > 0) & + (array.unc_vtx_ele_track_ecal_y < 42.0) & + (array.unc_vtx_ele_track_ecal_y < 13 - 0.26 * array.unc_vtx_ele_track_ecal_x) & + (array.unc_vtx_ele_track_ecal_y > 18 - 0.08 * array.unc_vtx_ele_track_ecal_x) & + (((array.unc_vtx_ele_track_ecal_x > -125) & (array.unc_vtx_ele_track_ecal_x < -95)) | + ((array.unc_vtx_ele_track_ecal_x > -85) & (array.unc_vtx_ele_track_ecal_x < -55))) + ) + ele_condition_2 = ( + (array.unc_vtx_ele_track_ecal_y < 0) & + (array.unc_vtx_ele_track_ecal_y < -23) & + (array.unc_vtx_ele_track_ecal_y < -15 + 0.08 * array.unc_vtx_ele_track_ecal_x) & + (array.unc_vtx_ele_track_ecal_y > -18 + 0.22 * array.unc_vtx_ele_track_ecal_x) & + (((array.unc_vtx_ele_track_ecal_x > -75) & (array.unc_vtx_ele_track_ecal_x < -45)) | + ((array.unc_vtx_ele_track_ecal_x > -110) & (array.unc_vtx_ele_track_ecal_x < -95))) + ) + ele_fid = ele_condition_1 | ele_condition_2 + + pos_condition_1 = ( + (array.unc_vtx_pos_track_ecal_y > 0) & + (array.unc_vtx_pos_track_ecal_y < 42.0) & + (array.unc_vtx_pos_track_ecal_y < 13 - 0.26 * array.unc_vtx_pos_track_ecal_x) & + (array.unc_vtx_pos_track_ecal_y > 18 - 0.08 * array.unc_vtx_pos_track_ecal_x) & + (((array.unc_vtx_pos_track_ecal_x > -125) & (array.unc_vtx_pos_track_ecal_x < -95)) | + ((array.unc_vtx_pos_track_ecal_x > -85) & (array.unc_vtx_pos_track_ecal_x < -55))) + ) + pos_condition_2 = ( + (array.unc_vtx_pos_track_ecal_y < 0) & + (array.unc_vtx_pos_track_ecal_y < -23) & + (array.unc_vtx_pos_track_ecal_y < -15 + 0.08 * array.unc_vtx_pos_track_ecal_x) & + (array.unc_vtx_pos_track_ecal_y > -18 + 0.22 * array.unc_vtx_pos_track_ecal_x) & + (((array.unc_vtx_pos_track_ecal_x > -75) & (array.unc_vtx_pos_track_ecal_x < -45)) | + ((array.unc_vtx_pos_track_ecal_x > -110) & (array.unc_vtx_pos_track_ecal_x < -95))) + ) + pos_fid = pos_condition_1 | pos_condition_2 + + elif isMC: + ele_condition_1 = ( + (array.unc_vtx_ele_track_ecal_y > 0) & + (array.unc_vtx_ele_track_ecal_y > 23) & + (array.unc_vtx_ele_track_ecal_y > 15 - 0.1 * array.unc_vtx_ele_track_ecal_x) & + (array.unc_vtx_ele_track_ecal_y < 12 - 0.3 * array.unc_vtx_ele_track_ecal_x) & + (((array.unc_vtx_ele_track_ecal_x > -75) & (array.unc_vtx_ele_track_ecal_x < -50)) | + ((array.unc_vtx_ele_track_ecal_x > -130) & (array.unc_vtx_ele_track_ecal_x < -95))) + ) + ele_condition_2 = ( + (array.unc_vtx_ele_track_ecal_y < 0) & + (array.unc_vtx_ele_track_ecal_y < -22) & + (array.unc_vtx_ele_track_ecal_y < -15 + 0.1 * array.unc_vtx_ele_track_ecal_x) & + (array.unc_vtx_ele_track_ecal_y > -15 + 0.25 * array.unc_vtx_ele_track_ecal_x) & + (((array.unc_vtx_ele_track_ecal_x > -120) & (array.unc_vtx_ele_track_ecal_x < -94)) | + ((array.unc_vtx_ele_track_ecal_x > -75) & (array.unc_vtx_ele_track_ecal_x < -50))) + ) + ele_fid = ele_condition_1 | ele_condition_2 + + pos_condition_1 = ( + (array.unc_vtx_pos_track_ecal_y > 0) & + (array.unc_vtx_pos_track_ecal_y > 23) & + (array.unc_vtx_pos_track_ecal_y > 15 - 0.1 * array.unc_vtx_pos_track_ecal_x) & + (array.unc_vtx_pos_track_ecal_y < 12 - 0.3 * array.unc_vtx_pos_track_ecal_x) & + (((array.unc_vtx_pos_track_ecal_x > -75) & (array.unc_vtx_pos_track_ecal_x < -50)) | + ((array.unc_vtx_pos_track_ecal_x > -130) & (array.unc_vtx_pos_track_ecal_x < -95))) + ) + pos_condition_2 = ( + (array.unc_vtx_pos_track_ecal_y < 0) & + (array.unc_vtx_pos_track_ecal_y < -22) & + (array.unc_vtx_pos_track_ecal_y < -15 + 0.1 * array.unc_vtx_pos_track_ecal_x) & + (array.unc_vtx_pos_track_ecal_y > -15 + 0.25 * array.unc_vtx_pos_track_ecal_x) & + (((array.unc_vtx_pos_track_ecal_x > -120) & (array.unc_vtx_pos_track_ecal_x < -94)) | + ((array.unc_vtx_pos_track_ecal_x > -75) & (array.unc_vtx_pos_track_ecal_x < -50))) + ) + pos_fid = pos_condition_1 | pos_condition_2 + + array['ele_fid'] = ele_fid + array['pos_fid'] = pos_fid + + filtered_array = array[ele_fid & pos_fid] + return filtered_array + +def fit_w_gaussian(histo, nsigma=2.0): + + #Seed initial fit with a simple gaussian fit over xmin and xmax + xmin = histo.GetXaxis().GetBinLowEdge(histo.FindFirstBinAbove(0)) + xmax = histo.GetXaxis().GetBinLowEdge(histo.FindLastBinAbove(0)) + fitfunc = r.TF1("gaus","gaus", xmin, xmax) + fitRes = histo.Fit(fitfunc,"QMES") + params = fitRes.Parameters() + errors = fitRes.Errors() + chi2 = fitRes.Chi2() + ndf = fitRes.Ndf() + + #Iterate fit by randomly seeding sigma, keep fit with best chi2 + best_chi2ndf = chi2/ndf + best_params = params + best_errors = errors + for i in range(80): + norm = best_params[0] + mu = best_params[1] + sigma = np.random.uniform(100,200)*0.01 + + #establish new fit range + xmin = mu - nsigma*sigma + xmax = mu + nsigma*sigma + fitfunc.SetParameter(0, norm) + fitfunc.SetParameter(1, mu) + fitfunc.SetParameter(2, sigma) + fitRes = histo.Fit(fitfunc,"QMES","", xmin, xmax) + + params = fitRes.Parameters() + errors = fitRes.Errors() + chi2ndf = fitRes.Chi2()/fitRes.Ndf() + + if chi2ndf < best_chi2ndf: + best_params = params + best_chi2ndf = chi2ndf + best_errors = errors + + #Refit with best parameters + fitfunc.FixParameter(0, best_params[0]) + fitfunc.FixParameter(1, best_params[1]) + fitfunc.FixParameter(2, best_params[2]) + sigma = best_params[2] + xmin = mu - nsigma*sigma + xmax = mu + nsigma*sigma + histo.Fit(fitfunc,"QMES","", xmin, xmax) + print(fitfunc.GetParError(2)) + return histo, best_params, best_errors, best_chi2ndf + +#======================================================================================================================================= +# LOAD DATA +#======================================================================================================================================= + +samples = {} +branches = ["unc_vtx_mass","unc_vtx_psum", "unc_vtx_ele_track_t", "unc_vtx_pos_track_t","unc_vtx_ele_track_ecal_x", + "unc_vtx_ele_track_ecal_y", "unc_vtx_pos_track_ecal_x", "unc_vtx_pos_track_ecal_y", + "unc_vtx_ele_clust_x", "unc_vtx_ele_clust_y", "unc_vtx_pos_clust_x", "unc_vtx_pos_clust_y", + "unc_vtx_ele_track_px","unc_vtx_ele_track_py","unc_vtx_pos_track_px","unc_vtx_pos_track_py", + "unc_vtx_ele_track_pz", "unc_vtx_pos_track_pz", + "unc_vtx_ele_track_nhits", "unc_vtx_pos_track_nhits", + "unc_vtx_px", "unc_vtx_py", "unc_vtx_pz"] + +# Read Data Mollers Unconstrained +infile = '/sdf/group/hps/user-data/alspellm/run/new_run_scripts/mollers/hadd/hadd-sample0-moller-ana.root' +cut_expression = ('( (unc_vtx_ele_track_t > -3.0) & (unc_vtx_pos_track_t > -3.0) & (unc_vtx_ele_track_t < 2.5) & (unc_vtx_pos_track_t < 2.5) & (unc_vtx_psum > 2.1) & (unc_vtx_psum < 2.45))') +selection = 'vtxana_Tight_nocuts' +samples['data'] = load_data(infile,selection, cut_expression=cut_expression, expressions=branches) +samples['data']['vertex_psum'] = np.sqrt(np.square(samples['data'].unc_vtx_px) + np.square(samples['data'].unc_vtx_py) + np.square(samples['data'].unc_vtx_pz)) + +# Read MC mollers unconstrained no smearing +infile = '/sdf/group/hps/user-data/alspellm/run/new_run_scripts/mollers/hadd/hadd-molv4-beamv6-PhysicsRun2016-Pass2_iss650_singles0_ana.root' +cut_expression = ('( abs(unc_vtx_ele_track_t - unc_vtx_pos_track_t) < 2.5) & (unc_vtx_psum > 2.1) & (unc_vtx_psum < 2.45) ') +samples['mc'] = load_data(infile,selection, cut_expression = cut_expression, expressions=branches) +samples['mc']['vertex_psum'] = np.sqrt(np.square(samples['mc'].unc_vtx_px) + np.square(samples['mc'].unc_vtx_py) + np.square(samples['mc'].unc_vtx_pz)) + +# Read MC mollers unconstrained with smearing +infile = '/sdf/group/hps/user-data/alspellm/run/new_run_scripts/mollers/hadd/hadd-molv4-beamv6_HPS-PhysicsRun2016-Pass2_iss650_singles0_ana_smeared_topbot_corr.root' +samples['mc_smear'] = load_data(infile,selection, cut_expression=cut_expression, expressions=branches) +samples['mc_smear']['vertex_psum'] = np.sqrt(np.square(samples['mc_smear'].unc_vtx_px) + np.square(samples['mc_smear'].unc_vtx_py) + np.square(samples['mc_smear'].unc_vtx_pz)) + + +# Apply fiducial cuts +samples['data'] = TrackFiducial(samples['data'], isData=True) +samples['data_cons'] = TrackFiducial(samples['data_cons'], isData=True) +samples['mc'] = TrackFiducial(samples['mc'], isData=False, isMC=True) +samples['mc_smear'] = TrackFiducial(samples['mc_smear'], isData=False, isMC=True) + +#Cut on angular relationship? +samples['data']['theta_1'] = np.arctan( np.sqrt(np.square(samples['data'].unc_vtx_ele_track_py) + np.square(samples['data'].unc_vtx_ele_track_px))/samples['data'].unc_vtx_ele_track_pz ) +samples['data']['theta_2'] = np.arctan( np.sqrt(np.square(samples['data'].unc_vtx_pos_track_py) + np.square(samples['data'].unc_vtx_pos_track_px))/samples['data'].unc_vtx_pos_track_pz ) +samples['data_cons']['theta_1'] = np.arctan( np.sqrt(np.square(samples['data_cons'].unc_vtx_ele_track_py) + np.square(samples['data_cons'].unc_vtx_ele_track_px))/samples['data_cons'].unc_vtx_ele_track_pz ) +samples['data_cons']['theta_2'] = np.arctan( np.sqrt(np.square(samples['data_cons'].unc_vtx_pos_track_py) + np.square(samples['data_cons'].unc_vtx_pos_track_px))/samples['data_cons'].unc_vtx_pos_track_pz ) +samples['mc']['theta_1'] = np.arctan( np.sqrt(np.square(samples['mc'].unc_vtx_ele_track_py) + np.square(samples['mc'].unc_vtx_ele_track_px))/samples['mc'].unc_vtx_ele_track_pz ) +samples['mc']['theta_2'] = np.arctan( np.sqrt(np.square(samples['mc'].unc_vtx_pos_track_py) + np.square(samples['mc'].unc_vtx_pos_track_px))/samples['mc'].unc_vtx_pos_track_pz ) +samples['mc_smear']['theta_1'] = np.arctan( np.sqrt(np.square(samples['mc_smear'].unc_vtx_ele_track_py) + np.square(samples['mc_smear'].unc_vtx_ele_track_px))/samples['mc_smear'].unc_vtx_ele_track_pz ) +samples['mc_smear']['theta_2'] = np.arctan( np.sqrt(np.square(samples['mc_smear'].unc_vtx_pos_track_py) + np.square(samples['mc_smear'].unc_vtx_pos_track_px))/samples['mc_smear'].unc_vtx_pos_track_pz ) + +#======================================================================================================================================= +# PLOT INVARIANT MASS +#======================================================================================================================================= + +# Initialize invariant mass histograms +invm_h = ( + hist.Hist.new + .StrCategory(list(samples.keys()), name='sel') + .Reg(100, 30.0, 70,label='Invariant Mass [MeV]') + .Double() +) + +# Fill histograms +invm_h.fill('data', samples['data'].unc_vtx_mass*1000.)#, weight = 1./len(samples['data'].unc_vtx_mass)) +invm_h.fill('data_cons', samples['data_cons'].unc_vtx_mass*1000.)#, weight = 1./len(samples['data_cons'].unc_vtx_mass)) +invm_h.fill('mc', samples['mc'].unc_vtx_mass*1000.)#, weight = 0.7/len(samples['mc'].unc_vtx_mass)) +invm_h.fill('mc_smear', samples['mc_smear'].unc_vtx_mass*1000.)#, weight = 1.7/len(samples['mc_smear'].unc_vtx_mass)) + +# Plot invariant mass histograms +fig, ax = plt.subplots(figsize=(25,15)) +invm_h.plot() +plt.legend() +plt.xlabel('Invariant Mass [MeV]') + + +# Plot showing fiducial cuts +fig, ax = plt.subplots(figsize=(25,15)) +plt.scatter(samples['mc'].unc_vtx_ele_track_ecal_x, samples['mc'].unc_vtx_ele_track_ecal_y) +plt.scatter(samples['mc'].unc_vtx_pos_track_ecal_x, samples['mc'].unc_vtx_pos_track_ecal_y) +plt.xlabel('Track at Ecal x [mm]') +plt.ylabel('Track at Ecal y [mm]') +# Generate x data points +x = np.linspace(-140, -40, 1000) + +# Data Cuts +# Positive Y Region +y1_data_pos = 13 - 0.26 * x +y2_data_pos = 18 - 0.08 * x + +# Negative Y Region +y1_data_neg = -15 + 0.08 * x +y2_data_neg = -18 + 0.22 * x +# Plot Data Cuts +plt.plot(x, y1_data_pos, label='Data Cut y < 13 - 0.26 * x',color='black') +plt.plot(x, y2_data_pos, label='Data Cut y > 18 - 0.08 * x',color='black') +plt.plot(x, y1_data_neg, label='Data Cut y < -15 + 0.08 * x',color='black') +plt.plot(x, y2_data_neg, label='Data Cut y > -18 + 0.22 * x',color='black') +plt.text(-100,0,'Data') +plt.xlim(-160,-35) +plt.ylim(-50,50) + + +# Plot angles theta1 and theta2 around beam axis +# Possibly useful to get higher purity Moller sample +fig, ax = plt.subplots(figsize=(25,15)) +plt.scatter(samples['mc'].theta_1, samples['mc'].theta_2) +plt.xlabel('arctan(sqrt(px^2+py^2)/pz) track 1') +plt.ylabel('arctan(sqrt(px^2+py^2)/pz) track 2') +plt.plot([0.02, 0.05],[0.05,0.02], color='red') +plt.plot([0.02, 0.058],[0.058,0.02], color='red') +coefficients = np.polyfit([0.02, 0.05],[0.05,0.02], 1) +slope, intercept = coefficients +coefficients = np.polyfit([0.02, 0.058],[0.058,0.02], 1) +slope, intercept = coefficients + + +fig, ax = plt.subplots(figsize=(25,15)) +plt.scatter(samples['data'].theta_1, samples['data'].theta_2) +plt.xlabel('arctan(sqrt(px^2+py^2)/pz) track 1') +plt.ylabel('arctan(sqrt(px^2+py^2)/pz) track 2') +plt.plot([0.02, 0.05],[0.05,0.02], color='red') +plt.plot([0.02, 0.058],[0.058,0.02], color='red') +coefficients = np.polyfit([0.02, 0.05],[0.05,0.02], 1) +slope, intercept = coefficients +coefficients = np.polyfit([0.02, 0.058],[0.058,0.02], 1) +slope, intercept = coefficients + +#======================================================================================================================================= +# FIT MOLLER MASS +#======================================================================================================================================= + +# Convert invariant mass histograms to ROOT histograms +moller_masses_h = {} +for sname, sample in samples.items(): + h = utils.cnvHistogramToROOT(invm_h[sname,:]) + moller_masses_h[sname] = h + + +# Fit the data Moller mass peak +histo = moller_masses_h['data'] +fit_histo, params, errors, chi2ndf = fit_w_gaussian(histo, nsigma=2.0) + +# Plot data fit +text = [f'\mu = {round(params[1],2)} \pm {round(errors[1],3)}',f'\sigma = {round(params[2],2)} \pm {round(errors[2],3)}',f'\chi^{2}/n.d.f = {round(chi2ndf,4)}'] +xmin = params[1] - 4.0*params[2] +xmax = params[1] + 4.0*params[2] +ymax = params[0]*1.1 +c = utils.drawTH1s([histo], histo.GetName(), drawOpts=['hist'],xrange=[xmin,xmax], yrange=[0.0,ymax],size=(2040,1080), text=text, text_pos=[0.2,0.78], line_spacing=0.05) +c.SaveAs('moller_fit_data.png') + +# Fit the MC Moller mass peak without FEE calibrated momentum smearing +histo = moller_masses_h['mc'] +fit_histo, params, errors, chi2ndf = fit_w_gaussian(histo, nsigma=2.0) + +# Plot MC Moller unsmeared fit +text = [f'\mu = {round(params[1],2)} \pm {round(errors[1],3)}',f'\sigma = {round(params[2],2)} \pm {round(errors[2],3)}',f'\chi^{2}/n.d.f = {round(chi2ndf,4)}'] +xmin = params[1] - 4.0*params[2] +xmax = params[1] + 4.0*params[2] +ymax = params[0]*1.1 +c = utils.drawTH1s([histo], histo.GetName(), drawOpts=['hist'],xrange=[xmin,xmax], yrange=[0.0,ymax],size=(2040,1080), text=text, text_pos=[0.2,0.78], line_spacing=0.05) +c.SaveAs('moller_fit_mc_unsmeared.png') + + +# Fit the MC Moller mass peak WITH FEE calibrated momentum smearing applied +histo = moller_masses_h['mc_smear'] +fit_histo, params, errors, chi2ndf = fit_w_gaussian(histo, nsigma=2.0) + +# Plot MC Moller SMEARED fit +text = [f'\mu = {round(params[1],2)} \pm {round(errors[1],3)}',f'\sigma = {round(params[2],2)} \pm {round(errors[2],3)}',f'\chi^{2}/n.d.f = {round(chi2ndf,4)}'] +xmin = params[1] - 4.0*params[2] +xmax = params[1] + 4.0*params[2] +ymax = params[0]*1.1 +c = utils.drawTH1s([histo], histo.GetName(), drawOpts=['hist'],xrange=[xmin,xmax], yrange=[0.0,ymax],size=(2040,1080), text=text, text_pos=[0.2,0.78], line_spacing=0.05) +c.SaveAs('moller_fit_mc_smeared.png') diff --git a/plotUtils/simps/plot_final_results.py b/plotUtils/simps/plot_final_results.py new file mode 100644 index 000000000..cd0772ba3 --- /dev/null +++ b/plotUtils/simps/plot_final_results.py @@ -0,0 +1,208 @@ +#!/usr/bin/python3 +import os +import awkward as ak +import numpy as np +import hist +from hist import Hist +import uproot +import ROOT as r +import copy +import matplotlib.pyplot as plt +import matplotlib as mpl +import matplotlib.gridspec as gridspec +import sys +import math + +#Import simp class where systematics info is stored +import simp_signal_2016 + +#format mpl plots. Should make a template... +plt.rcParams.update({'font.size': 60, # Font size for text + 'axes.titlesize': 60, # Font size for titles + 'axes.labelsize': 60, # Font size for axis labels + 'xtick.labelsize': 60, # Font size for x-axis tick labels + 'ytick.labelsize': 60, # Font size for y-axis tick labels + 'lines.linewidth':5.0, + 'legend.fontsize': 60}) # Font size for legend +plt.rcParams['font.family'] = 'DejaVu Sans' + +import argparse +parser = argparse.ArgumentParser(description='') +parser.add_argument('--infile_signal_search', type=str, default='.simp_2016_results/unblinded_result_100pct.root') +parser.add_argument('--infile_signal_opt_interval', type=str, default='.simp_2016_results/unblinded_result_100pct.root') +parser.add_argument('--outdir', type=str, default='./search_results') +parser.add_argument('--mpifpi', type=float, default=4.*np.pi) + +args = parser.parse_args() +outdir = args.outdir +######################################################################################################################################## + +#Grab search results (exp bkg, obs, local pvalues) +with uproot.open(args.infile_signal_search) as f: + expected_bkg = f['expected_background'].values() + expected_bkg_errlow = f['expected_background'].errors('low')[1] + expected_bkg_errhigh = f['expected_background'].errors('high')[1] + Nobs = f['Nobs'].values() + local_pvalue = f['local_pvalue'].values() + local_pvalue_errlow = f['local_pvalue'].errors('low')[1] + local_pvalue_errhigh = f['local_pvalue'].errors('high')[1] + +#Calculate the "Look-Elsewhere Effect" correction. Total search window mass range divided by avg mass resolution +search_window = 1.5 #used in final analysis +signalProcessor = simp_signal_2016.SignalProcessor(args.mpifpi, search_window) +masses = expected_bkg[0] +masses = masses[np.where(masses <= 124.0)[0]] #I restricted search for signal to 124 MeV +avg_resolution = np.average( np.array([signalProcessor.mass_resolution(x) for x in masses]) ) +look_elsewhere = (masses[-1] - masses[0])/(avg_resolution) +print(f'Average mass resolution: {avg_resolution}') +print(f'Look elsewhere effect: {look_elsewhere}') +thresholds = [] +thresholds_lew = [] #thresholds corrected for look elsewhere +from scipy.stats import norm +for nsigma in [1,2,3,4,5]: + gaus_cdf = norm.cdf(nsigma) + threshold = (1.0 - gaus_cdf)/look_elsewhere + thresholds_lew.append(threshold) + thresholds.append((1.0 - gaus_cdf)) + +#top plot: expected background and observed events for all search windows +fig, ax = plt.subplots(2, 1, figsize=(60,60), gridspec_kw={'height_ratios': [3, 2]}) +gs = gridspec.GridSpec(2, 1, height_ratios=[10, 1]) +plt.subplot(2,1,1) +plt.errorbar(expected_bkg[0], expected_bkg[1], yerr=(expected_bkg_errlow, expected_bkg_errhigh), marker='o', markersize=20, color='blue', label='Estimated Background') +plt.plot(Nobs[0], Nobs[1], marker='o', color='black', markersize=20, label='Observed Events') +plt.legend() +plt.xlabel('$V_{D}$ Invariant Mass Search Window [MeV]') +plt.ylabel('Number Of Events') +plt.xlim(29.5,124.0) + +#bottom plot: local p-values and significance thresholds +plt.subplot(2,1,2) +plt.errorbar(local_pvalue[0], local_pvalue[1], yerr=(local_pvalue_errlow, local_pvalue_errhigh), marker='o', markersize=20, color='black', label='Local p-value') +from scipy.stats import norm +for n,nsigma in enumerate(thresholds): + if n < len(thresholds)-1: + plt.axhline(y=nsigma, linestyle='--', linewidth=4.0, color='black') + plt.axhline(y=thresholds_lew[n], linestyle='--', linewidth=4.0, color='red') + else: + plt.axhline(y=nsigma, linestyle='--', linewidth=4.0, color='black', label='Local $N\sigma$') + plt.axhline(y=thresholds_lew[n], linestyle='--', linewidth=4.0, color='red', label='Global LEE $N\sigma$') +plt.xlabel('$V_{D}$ Invariant Mass Search Window [MeV]') +plt.ylabel('p-value') +plt.legend() + +#Identify search window with larges p-value fluctuation (smallest p-value) +pvalmin = np.min(local_pvalue[1]) +pvalmin_mass = local_pvalue[0][np.argmin(local_pvalue[1])] +nsigma_local = np.round(norm.ppf(1 - pvalmin / 2),1) +nsigma_global = np.round(norm.ppf(1 - (pvalmin*look_elsewhere) / 2),1) +t = plt.text(80.0, .9e-4, f'Smallest p-value: {pvalmin} at {pvalmin_mass} MeV\nLocal Significance: {nsigma_local}$\sigma$\nGlobal Significance: {nsigma_global}$\sigma$') +t.set_bbox(dict(facecolor='white', alpha=1.0, edgecolor='black')) +plt.yscale('log') +plt.ylim(1e-6,1.1) +plt.xlim(29.5,124.0) +plt.savefig(f'{outdir}/signal_search_results.png') + +########################################################################################################################## +#Calculate exclusion contour with systematic uncertainties included +infile_signal_opt_interval = args.infile_signal_opt_interval +#no systematics included in these results. They are applied in this script by reducing the expected signal rate. +#Published result may prefer to generate more MC signal masses, re-calculate OIM, use ratio (OIM w sys)/(OIM w/o sys)... +with uproot.open(infile_signal_opt_interval) as f: + exclusion_h = f['sensitivity_ap_h'].to_hist() + expected_signal_h = f['total_yield_ap_h'].to_hist() + excluded_signal_h = f['excluded_signal_ap_h'].to_hist() + +sysProc = simp_signal_2016.SignalProcessor(args.mpifpi, search_window) +sysProc.systematic_uncertainties() +masses = expected_signal_h.axes[0].centers +sys_uncertainties = [] +for mass in masses: + #systematics are accessed through SignalProcessor + rad_targ_nom, rad_targ_mpt5, rad_targ_ppt5, simp_targ, mass_unc, radfrac_unc = sysProc.evaluate_polynomials(mass) + #for rad targ, take ratio of off-nominal to nominal. Max value is used. + rad_targ = np.maximum(rad_targ_mpt5/rad_targ_nom, rad_targ_ppt5/rad_targ_nom) - 1.0 + if rad_targ < 0.0: + rad_targ = 0.0 + #for simp targ, if ratio is greater than 1, set to 1 + simp_targ = 1.0-simp_targ + if simp_targ < 0.0: + simp_targ = 0.0 + sys_unc = np.sqrt( rad_targ**2 + simp_targ**2 + radfrac_unc**2 + mass_unc**2 ) + sys_uncertainties.append(sys_unc) + +#Plot systematic uncertaint as function of mass +fig, ax = plt.subplots(figsize=(40,20)) +plt.plot(masses, sys_uncertainties, marker='o',markersize=10, color='black') +plt.xlabel('A\' Invariant Mass [MeV]') +plt.ylabel('Systematic Uncertainty') +plt.ylim(0.0, 0.2) +plt.axhline(y=0.05, linestyle='--', color='darkred') +plt.axhline(y=0.1, linestyle='--', color='darkred') +plt.text(160, 0.15, 'HPS Preliminary\n(incomplete)', horizontalalignment='center') +plt.savefig(f'{outdir}/systematic_uncertainty_summary.png') + +#Rescale the expected signal according to the systematic uncertainty +values = expected_signal_h.values() +rescaled_values = np.zeros_like(values) +for m, mass in enumerate(masses): + rescaled_values[m, :] = values[m, :] * (1.0-sys_uncertainties[m]) +rescaled_signal_h = expected_signal_h.copy() +rescaled_signal_h[...] = rescaled_values + +#Divide rescaled signal by upper limit to get exclusion contour +new_exclusion = rescaled_values/excluded_signal_h.values() +rescaled_exclusion_h = exclusion_h.copy() +rescaled_exclusion_h[...] = new_exclusion + +#Apply systematics to expected signal... +fig, ax = plt.subplots(figsize=(40,30)) +rescaled_signal_h.plot() +plt.xlim(50,180) +plt.ylim(-6.5, -4.0) +plt.xlabel('A\' Invariant Mass [MeV]', fontsize=80) +plt.ylabel('$\log{\epsilon^2}$', fontsize=80) +plt.text(150, -4.4,f'HPS PRELIMINARY\n(Partial Systematics)', color='white', weight='bold', fontsize=80, horizontalalignment='center') +plt.text(156, -3.97,'Expected Signal', fontsize=80) +plt.savefig(f'{outdir}/expected_signal_2d.png') + +#save upper limit +fig, ax = plt.subplots(figsize=(40,30)) +excluded_signal_h.plot() +plt.xlim(50,180) +plt.ylim(-6.5, -4.0) +plt.xlabel('A\' Invariant Mass [MeV]', fontsize=80) +plt.ylabel('$\log{\epsilon^2}$', fontsize=80) +plt.text(165, -3.97,'Upper Limit', fontsize=80) +plt.text(150, -4.4,f'HPS PRELIMINARY', color='white', weight='bold', fontsize=80, horizontalalignment='center') +plt.savefig(f'{outdir}/upper_limit_2d.png') + +#Plot exclusion contour with systematics +fig, ax = plt.subplots(figsize=(40,30)) +rescaled_exclusion_h.plot(cmin=0.0) +plt.xlim(50,180) +plt.ylim(-6.5, -4.0) +plt.xlabel('A\' Invariant Mass [MeV]', fontsize=80) +plt.ylabel('$\log{\epsilon^2}$', fontsize=80) +plt.text(150, -4.4,f'HPS PRELIMINARY\n(Partial Systematics)', color='white', weight='bold', fontsize=80, horizontalalignment='center') +recenters_x = rescaled_exclusion_h.axes[0].centers +recenters_y = rescaled_exclusion_h.axes[1].centers +revalues = rescaled_exclusion_h.values() +reX, reY = np.meshgrid(recenters_x, recenters_y) +contour_levels = [1.0] +colors = ['red'] +recontour = ax.contour(reX, reY, revalues.T, levels=contour_levels, colors=colors) # Transpose values for correct orientation +ax.clabel(recontour, inline=1, fontsize=60) +plt.savefig(f'{outdir}/exclusion_contour_2d.png') + +#Write exclusion contour to .dat files +import csv +for l,level in enumerate(contour_levels): + for line in range(len(recontour.allsegs[l])): + print(line) + with open(f'{outdir}/exclusion_mpifpi_{args.mpifpi}_contour_line{line+1}.dat', 'w', newline='') as csvfile: + contour_line = recontour.allsegs[l][line] + data = [[x[0], np.sqrt(10**x[1])] for x in contour_line] + data.append([contour_line[0][0], np.sqrt(10**contour_line[0][1])]) + writer = csv.writer(csvfile, delimiter=' ') + writer.writerows(data) diff --git a/plotUtils/simps/run_opt_interval.py b/plotUtils/simps/run_opt_interval.py new file mode 100644 index 000000000..86349c19c --- /dev/null +++ b/plotUtils/simps/run_opt_interval.py @@ -0,0 +1,478 @@ +#!/usr/bin/python3 +""" +This script runs the optimum interval method to calculate a 90% confidence upper limit on the SIMP signal rate as a function of +mass and epsilon. + +Load in the data (in the form of the flat tuple output by the hpstr vertex analysis processor) and apply all selection criteria. + *The vertex ana processor applies Preselection, and a few Tight cuts, but you need apply all remaining tight cuts. +Load in the MC signal (for each generated mass), and apply all selection criteria. + +The reconstructed vertex z distribution of the remaining data events is transformed into a normalized uniform distribution according +to the expected MC signal shape in reconstructed vertex z. Since this shape is a function of epsilon^2, upper limit depends on both +mass and epsilon^2. + +Calculating the upper limit here requires an external lookuptable that is generated using cmax.py. +""" +import os +import awkward as ak +import numpy as np +import hist +from hist import Hist +import uproot +import math +import ROOT as r +import simp_signal_2016 +from simp_theory_equations import SimpEquations as simpeqs +import copy +import pickle + +#======================================================================================================================================= +# FUNCTIONS +#======================================================================================================================================= +def kLargestIntervals(list_of_energies, spectrumCDF = lambda x: x): + """ + Returns a list of the sizes of the K-largest intervals in that run according to the energy spectrum (given as a CDF). + That is, kLargestIntervals(...)[i] is the size of the largest interval containing i events, where ‘largest’ is defined above. + + * Transform energies to energy cumulants + * Add events at 0 and 1 + * Foreach k, compute interval sizes, take max + """ + answer = {} + list_of_energies.sort() + energy_cumulants = spectrumCDF(list_of_energies) + for i,interval_size in enumerate(range(len(energy_cumulants))): + if (1 + interval_size) >= len(energy_cumulants): + continue + + temp_data = energy_cumulants.copy() + gap_sizes = temp_data[(1+interval_size):] - temp_data[0:-1*(1 + interval_size)] + + answer[interval_size] = np.max(gap_sizes) + return answer +#======================================================================================================================================= +# INITIALIZE +#======================================================================================================================================= +# --outfilename: Specify output file name. +# --tenpct: If True run OIM on 10% data (or single hpstr vertex ana output tuple). +# --highPsum: If True, run OIM in high Psum (CR). +# --mpifpi: Ratio of dark pion mass to dark pion decay constant (benchmarks are 3 and 4pi). +# --signal_sf: Scale the signal. +# --nsigma: Size of the signal invariant mass search window (+-nsigma) + +import argparse +parser = argparse.ArgumentParser(description='Process some inputs.') +parser.add_argument('--outfilename', type=str, default='oim_results') +parser.add_argument('--tenpct', type=int, default=1) +parser.add_argument('--highPsum', type=int, default=0) +parser.add_argument('--mpifpi', type=float, default=4*np.pi) +parser.add_argument('--signal_sf', type=float, default=1.0) +parser.add_argument('--nsigma', type=float, default=1.5) + +args = parser.parse_args() +outfilename = args.outfilename +mpifpi = args.mpifpi +nsigma = args.nsigma +signal_sf = args.signal_sf +tenpct = args.tenpct +print(f'Search Window Size: +-', nsigma) + +#======================================================================================================================================= +# LOAD DATA +#======================================================================================================================================= +#Initialize Signal Processor +signalProcessor = simp_signal_2016.SignalProcessor(mpifpi=mpifpi, nsigma=nsigma) + +data = ak.Array([]) +if args.tenpct: + # If tenpct True, run OIM on 10% data (or a single file) + outfilename = f'{outfilename}_10pct' + inv_mass_range = (30,124) + print('Loading 10% Data') + branches = ["unc_vtx_mass","unc_vtx_psum", "unc_vtx_ele_track_z0", "unc_vtx_pos_track_z0", "unc_vtx_z", "unc_vtx_proj_sig"] + infile = '/sdf/group/hps/user-data/alspellm/2016/data/full_hadd_blpass4c_ana.root' + selection = 'vtxana_Tight_2016_simp_reach_SR' + data = signalProcessor.load_data(infile,selection, expressions=branches) + data['weight'] = 1.0 + +else: + # Run OIM on 100% data (multiple files in a single directory) + outfilename = f'{outfilename}_100pct' + print('Loading 100% Data') + inv_mass_range = (30,200) + branches = ["unc_vtx_mass","unc_vtx_psum", "unc_vtx_ele_track_z0", "unc_vtx_pos_track_z0", "unc_vtx_z", "unc_vtx_proj_sig"] + indir = '/fs/ddn/sdf/group/hps/users/alspellm/data_storage/pass4kf/pass4kf_ana_20240513' + + # If highPsum is True, can look at all masses + if args.highPsum: + selection = 'vtxana_Tight_2016_simp_reach_CR' + mass_safety = 'unc_vtx_mass*1000. >= 0' + + # If highPsum is False, can only look at 10% data before unblinding + else: + selection = 'vtxana_Tight_2016_simp_reach_SR' + #mass_safety = 'unc_vtx_mass*1000. > 135' #CANT LOOK BELOW THIS MASS UNTIL UNBLINDING! + mass_safety = 'unc_vtx_mass*1000. > 0.0' #UNBLINDED! + inv_mass_range = (30, 124) + + # Loop over all input files and combine into single data array + for filename in sorted(os.listdir(indir)): + if not filename.endswith('.root'): + continue + run = filename.split('_')[4] + print('Loading Run ', run) + infile = os.path.join(indir,filename) + data = ak.concatenate([data, signalProcessor.load_data(infile, selection, cut_expression=mass_safety, expressions=branches)]) + data['weight'] = 1.0 + + +# Load the differential radiative trident rate lokup table. This scales the expected signal to the data +print('Load lookup table') +cr_data = '/sdf/group/hps/user-data/alspellm/2016/data/hadd_BLPass4c_1959files.root' # If using 10% data. +full_lumi_path = '/fs/ddn/sdf/group/hps/users/alspellm/data_storage/pass4kf/pass4kf_ana_20240513' # If using 100% data. +preselection = "vtxana_Tight_nocuts" +signal_mass_range = [x for x in range(20,130,1)] +signalProcessor.set_diff_prod_lut(cr_data, preselection, signal_mass_range, tenpct, full_lumi_path) + +#======================================================================================================================================= +# INITIALIZE HISTOGRAMS +#======================================================================================================================================= + +masses = [x for x in range(inv_mass_range[0], inv_mass_range[-1]+2,2)] +ap_masses = [round(x*signalProcessor.mass_ratio_ap_to_vd,1) for x in masses] +eps2_range = np.logspace(-4.0,-8.0,num=1000) +logeps2_range = np.log10(eps2_range) +min_eps = min(np.log10(eps2_range)) +max_eps = max(np.log10(eps2_range)) +num_bins = len(eps2_range) + +exclusion_conf_h = ( + hist.Hist.new + .Reg(len(masses)-1, np.min(masses),np.max(masses),label='v_{D} Invariant Mass [MeV]') + .Reg(len(eps2_range), min_eps,max_eps,label=r'$log10(\epsilon^2)$') + .Double() +) +exclusion_bestk_h = ( + hist.Hist.new + .Reg(len(masses)-1, np.min(masses),np.max(masses),label='v_{D} Invariant Mass [MeV]') + .Reg(len(eps2_range), min_eps,max_eps,label=r'$log10(\epsilon^2)$') + .Double() +) +total_yield_h = ( + hist.Hist.new + .Reg(len(masses)-1, np.min(masses),np.max(masses),label='v_{D} Invariant Mass [MeV]') + .Reg(len(eps2_range), min_eps,max_eps,label=r'$log10(\epsilon^2)$') + .Double() +) +excluded_signal_h = ( + hist.Hist.new + .Reg(len(masses)-1, np.min(masses),np.max(masses),label='v_{D} Invariant Mass [MeV]') + .Reg(len(eps2_range), min_eps,max_eps,label=r'$log10(\epsilon^2)$') + .Double() +) +sensitivity_h = ( + hist.Hist.new + .Reg(len(masses)-1, np.min(masses),np.max(masses),label='v_{D} Invariant Mass [MeV]') + .Reg(len(eps2_range), min_eps,max_eps,label=r'$log10(\epsilon^2)$') + .Double() +) + +ap_masses = [signalProcessor.mass_ratio_ap_to_vd*x for x in range(inv_mass_range[0], inv_mass_range[-1]+2,2)] +total_yield_ap_h = ( + hist.Hist.new + .Reg(len(ap_masses), np.min(ap_masses),np.max(ap_masses),label='A\' Invariant Mass [MeV]') + .Reg(len(eps2_range), min_eps,max_eps,label=r'$log10(\epsilon^2)$') + .Double() +) +excluded_signal_ap_h = ( + hist.Hist.new + .Reg(len(ap_masses), np.min(ap_masses),np.max(ap_masses),label='A\' Invariant Mass [MeV]') + .Reg(len(eps2_range), min_eps,max_eps,label=r'$log10(\epsilon^2)$') + .Double() +) +sensitivity_ap_h = ( + hist.Hist.new + .Reg(len(ap_masses), np.min(ap_masses),np.max(ap_masses),label='A\' Invariant Mass [MeV]') + .Reg(len(eps2_range), min_eps,max_eps,label=r'$log10(\epsilon^2)$') + .Double() +) + +#Plot the excluded signal value right before reaching 90% confidence. Debugging purposes +excluded_signal_minus1_h = ( + hist.Hist.new + .Reg(len(masses)-1, np.min(masses),np.max(masses),label='v_{D} Invariant Mass [MeV]') + .Reg(len(eps2_range), min_eps,max_eps,label=r'$log10(\epsilon^2)$') + .Double() +) +exclusion_conf_minus1_h = ( + hist.Hist.new + .Reg(len(masses)-1, np.min(masses),np.max(masses),label='v_{D} Invariant Mass [MeV]') + .Reg(len(eps2_range), min_eps,max_eps,label=r'$log10(\epsilon^2)$') + .Double() +) + +#======================================================================================================================================= +# RUN OPTIMUM INTERVAL METHOD +#======================================================================================================================================= + +# Load OIM lookup table generated using cmax.py +lookuptable_path = '/fs/ddn/sdf/group/hps/users/alspellm/mc_storage/opt_int_lookuptable_max25_10ktoys_0.05steps_v2.p' + +# Number of toy events thrown for each mu in the loaded lookup table +ntrials = 10000 # This value is defined in cmax.py when the lookup table is generated. It MUST MATCH the value used. +with open(lookuptable_path, 'rb') as f: + lookupTable = pickle.load(f) + +# Open an output file to store the results +outfile = uproot.recreate(f'{outfilename}.root') + +# Run OIM for each MC generated signal mass (30-124 @ 2 MeV intervals) +for signal_mass in masses: + + # Initialize histograms + confidence_level_mass_h = ( + hist.Hist.new + .Reg(300, 0, 30.0,label='mu') + .Reg(len(eps2_range), min_eps,max_eps,label=r'$log10(\epsilon^2)$') + .Double() + ) + best_kvalue_mass_h = ( + hist.Hist.new + .Reg(300, 0, 30.0,label='mu') + .Reg(len(eps2_range), min_eps,max_eps,label=r'$log10(\epsilon^2)$') + .Double() + ) + + print(f'Signal Mass {signal_mass}') + + # Define the invariant mass search window boundaries based on the search window size + mass_low = signal_mass - signalProcessor.mass_resolution(signal_mass)*nsigma + mass_high = signal_mass + signalProcessor.mass_resolution(signal_mass)*nsigma + + # Build the final selection used in the signal search and apply to the data + zcut_sel = signalProcessor.zcut_sel(data) # zcut on target position at -4.3 mm. + vprojsig_sel = signalProcessor.vprojsig_sel(data) # Require target projected vertex significance < 2.0 + minz0_sel = signalProcessor.minz0_sel(data) # Cut on minimum track vertical impact parameter z0 (aka y0) + sameside_sel = signalProcessor.sameside_z0_cut(data) # Cut events where both tracks have same side z0 + masswindow_sel = signalProcessor.mass_sel(data, signal_mass) # Define search window mass boundaries + + # Set the Psum selection + if not args.tenpct and args.highPsum: + psum_sel = signalProcessor.psum_sel(data, case='cr') + elif args.tenpct and not args.highPsum: + psum_sel = signalProcessor.psum_sel(data, case='sr') + elif args.tenpct and args.highPsum: + psum_sel = signalProcessor.psum_sel(data, case='cr') + else: + psum_sel = signalProcessor.psum_sel(data, case='sr') + + # Combine the selections and apply to data + tight_sel = np.logical_and.reduce([zcut_sel, vprojsig_sel, sameside_sel, psum_sel, minz0_sel, masswindow_sel]) + data_z = data[tight_sel].unc_vtx_z + + #================================================================================================================================== + # LOAD MC SIGNAL + #================================================================================================================================== + + indir = '/sdf/group/hps/user-data/alspellm/2016/simp_mc/pass4b/beam/smeared_fixbeamspot' + + # hpstr MC ana processor output file that stores the pre-readout MC signal truth information + signal_pre_readout_path = lambda mass: f'/sdf/group/hps/user-data/alspellm/2016/simp_mc/pass4b/nobeam/mass_{mass}_simp_2pt3_slic_hadd_ana.root' + # hpstr vertex ana processor output tuple + signal_path = lambda mass: f'{indir}/mass_{mass}_hadd-simp-beam_ana_smeared_corr.root' + signal_selection = 'vtxana_radMatchTight_2016_simp_SR_analysis' + + # Calculate the total A' production rate per epsilon^2 + total_yield_per_epsilon2 = signalProcessor.total_signal_production_per_epsilon2(signal_mass) + print('Total Yield Per eps2: ', total_yield_per_epsilon2) + + # Load signal before tight selection + print('Load Signal ', signal_path(signal_mass)) + signal = signalProcessor.load_signal(signal_path(signal_mass), signal_pre_readout_path(signal_mass), signal_mass, signal_selection) + + # Build final selection for signal + zcut_sel = signalProcessor.zcut_sel(signal) + vprojsig_sel = signalProcessor.vprojsig_sel(signal) + minz0_sel = signalProcessor.minz0_sel(signal) + sameside_sel = signalProcessor.sameside_z0_cut(signal) + masswindow_sel = signalProcessor.mass_sel(signal, signal_mass) + + # Set the Psum selection + if not args.tenpct and args.highPsum: + psum_sel = signalProcessor.psum_sel(signal, case='cr') + elif args.tenpct and not args.highPsum: + psum_sel = signalProcessor.psum_sel(signal, case='sr') + elif args.tenpct and args.highPsum: + psum_sel = signalProcessor.psum_sel(signal, case='cr') + else: + psum_sel = signalProcessor.psum_sel(signal, case='sr') + print('UNBLINDED!') + + # Combine the selections and apply to MC signal + tight_sel = np.logical_and.reduce([zcut_sel, vprojsig_sel, sameside_sel, psum_sel, minz0_sel, masswindow_sel]) + signal = signal[tight_sel] + + #================================================================================================================================== + # CALCULATE UPPER LIMIT ON SIGNAL: As function of epsilon^2 + #================================================================================================================================== + for i, eps2 in enumerate(eps2_range): + signal = signalProcessor.get_exp_sig_eps2(signal_mass, signal, eps2) + total_yield = ak.sum(signal['reweighted_accxEff'])*total_yield_per_epsilon2*eps2 + if i%20 == 0: + print(f'eps2 = {eps2}') + print(total_yield) + + # Signal acceptance*efficiency*dark_vector_probability in reconstructed vertex z. + # This represents the shape of the signal in 1D. + exp_sig_eff_z = ( + hist.Hist.new + .Reg(140, -40.0,100.0,label=r'Recon z [mm]') + .Double() + ) + exp_sig_eff_z.fill(signal.unc_vtx_z, weight=signal.reweighted_accxEff*total_yield_per_epsilon2*eps2) + + # Convert the remaining events in data to a normalized uniform distribution in reconstructed vertex z according to signal shape. + data_uniform_z = ( + hist.Hist.new + .Reg(101, -0.005,1.005,label=r'Recon z [mm]') + .Double() + ) + + # Initialize an array to store the new data events that are transformed according to the signal shape. + # Add endpoints to the new uniform data array, 0 in front, and 1.0 and the end. + dataArray = np.zeros(len(data_z)+2) + dataArray[0] = 0.0 + for k in range (0, len(data_z)): + thisX = data_z[k] + dataArray[k+1] = total_yield - exp_sig_eff_z[hist.loc(thisX)::sum] #transformation + + dataArray[len(data_z)+1] = total_yield + dataArray = dataArray/total_yield # normalize distribution based on total signal rate + dataArray = np.nan_to_num(dataArray, nan=1.0) + dataArray[0] = 0.0 + dataArray.sort() + data_uniform_z.fill(dataArray) + + # Calculate maximum gaps with k events allowed between events + kints = kLargestIntervals(dataArray) + + # Loop through the lookup table to find what upper limit on the signal rate results in 90% confidence + mu_90p = 99999.9 + k_90p = -1 + conf_90p = -1.0 + + previous_mu = 999999.9 + previous_conf = -9.9 + # Loop over values of mu (mean expected signal rate) + for i,mu in enumerate(sorted(lookupTable.keys())): + best_k = -1 + best_conf = -1.0 # Store best confidence level across all values of k + + # Loop over all values of k (k events allowed in gap between data events) + for k in sorted(lookupTable[mu].keys()): + if k > len(kints)-1: + break + x = np.max(kints[k]) + conf = np.where(lookupTable[mu][k] < x)[0].size / (ntrials) # Confidence level for this mu and k + if conf > best_conf: + best_k = k + best_conf = conf + + # Debug histos + confidence_level_mass_h.fill(mu, np.log10(eps2), weight=best_conf) + best_kvalue_mass_h.fill(mu, np.log10(eps2), weight=best_k) + + # If the condience level is >= 90%, this value of mu is the upper limit + if best_conf >= 0.9: + mu_90p = mu + k_90p = best_k + conf_90p = best_conf + + # Fill debug histos + excluded_signal_minus1_h.fill(signal_mass, np.log10(eps2), weight=previous_mu) + exclusion_conf_minus1_h.fill(signal_mass, np.log10(eps2), weight=previous_conf) + break + + # More debug + previous_mu = mu + previous_conf = best_conf + + + # Fill OIM results in histograms + exclusion_conf_h.fill(signal_mass, np.log10(eps2), weight=conf_90p) + exclusion_bestk_h.fill(signal_mass, np.log10(eps2), weight=k_90p) + total_yield_h.fill(signal_mass, np.log10(eps2), weight=total_yield) + excluded_signal_h.fill(signal_mass, np.log10(eps2), weight=mu_90p) + sensitivity_h.fill(signal_mass, np.log10(eps2), weight=(total_yield/mu_90p)) + + total_yield_ap_h.fill(signalProcessor.mass_ratio_ap_to_vd*signal_mass, np.log10(eps2), weight=total_yield) + excluded_signal_ap_h.fill(signalProcessor.mass_ratio_ap_to_vd*signal_mass, np.log10(eps2), weight=mu_90p) + sensitivity_ap_h.fill(signalProcessor.mass_ratio_ap_to_vd*signal_mass, np.log10(eps2), weight=(total_yield/mu_90p)) + + # Save mass dependent histograms + outfile[f'masses/confidence_levels_{signal_mass}_h'] = confidence_level_mass_h + outfile[f'masses/best_kvalues_{signal_mass}_h'] = best_kvalue_mass_h + +# Save results across all masses +outfile['total_yield_h'] = total_yield_h +outfile['excluded_signal_h'] = excluded_signal_h +outfile['sensitivity_h'] = sensitivity_h +outfile['confidence_level_h'] = exclusion_conf_h +outfile['best_exclusion_k_h'] = exclusion_bestk_h + +outfile['total_yield_ap_h'] = total_yield_ap_h +outfile['excluded_signal_ap_h'] = excluded_signal_ap_h +outfile['sensitivity_ap_h'] = sensitivity_ap_h + +# Save debug plots +outfile['excluded_signal_minus1_h'] = excluded_signal_minus1_h +outfile['exclusion_conf_minus1_h'] = exclusion_conf_minus1_h + +''' +#Get sensitivity contour +sens = sensitivity_h.values() +print(sens) +indices = np.where(sens >= 0.1) +print(indices) +if len(indices) > 0: + massx = + x = np.array(sensitivity_h.axes[0].centers[indices[0]]) + y = np.array(sensitivity_h.axes[1].centers[indices[1]]) + #get the contour + cupper = [] + clower = [] + for i in set(x): + cupper = y[np.where(x==i)[0][-1]] + clower = y[np.where(x==i)[0][0]] + cupper = cupper + list(reversed(clower)) + cupper[0] + + contour_g = r.TGraph(len(x), np.array(cupper,dtype=float), np.array(clower,dtype=float)) + contour_g.SetName('exclusion_contour_g') + contour_g.SetTitle('Exclusion Contour;Invariant Mass [MeV];#epsilon^{2}') + outfile['contour_g'] = contour_g +''' + +''' +fig, ax = plt.subplots(figsize=(30,20)) +total_yield_h.plot() +#plt.xlim(20.0, 126.0) +#plt.ylim(-10.8, -3.8) +plt.show() + +fig, ax = plt.subplots(figsize=(30,20)) +excluded_signal_h.plot() +#plt.xlim(20.0, 126.0) +#plt.ylim(-10.8, -3.8) +#plt.show() + +fig, ax = plt.subplots(figsize=(30,20)) +sensitivity_h.plot(cmin=1.0) +#plt.xlim(20.0, 126.0) +#plt.ylim(-10.8, -3.8) +plt.show() +''' + + + + + diff --git a/plotUtils/simps/run_signal_search.py b/plotUtils/simps/run_signal_search.py new file mode 100644 index 000000000..1e75aa4ee --- /dev/null +++ b/plotUtils/simps/run_signal_search.py @@ -0,0 +1,464 @@ +#!/usr/bin/python3 +import os +import awkward as ak +import numpy as np +import hist +from hist import Hist +import uproot +import math +import ROOT as r +import simp_signal_2016 +from simp_theory_equations import SimpEquations as simpeqs +import copy +#======================================================================================================================================= +# Initialize +#======================================================================================================================================= +# --outfilename: Specify output file name. +# --trials: Number of toy mc trials used to calculate data significance (typically millions+). +# --highPsum: If set, use high psum selection (control region). +# --tenpct: If True search for signal in 10% data. If False search for signal in 100% data. + +import argparse +parser = argparse.ArgumentParser(description='Process some inputs.') +parser.add_argument('--outfilename', type=str, default='expected_background_test') +parser.add_argument('--trials', type=int, default=100000) +parser.add_argument('--highPsum', type=int, default=0) +parser.add_argument('--tenpct', type=int, default=1) + +args = parser.parse_args() +outfilename = args.outfilename +t0_trials = args.trials #number of toy mc trials used to calculat p-value +nsigma=1.5 #+-1.5*mass_res (3sigma wide in total) +left_nsigma=nsigma+4.0 #mass sideband edge (sideband is 4sigma wide) +right_nsigma=nsigma+4.0 +z0_floor_threshold=1000 #Defines the sideband size in min_z0. 1000 is arbitrary choice that performs well (no signal contam) +print(f'Search Window Size: +-', nsigma) + +#======================================================================================================================================= +# Read data +#======================================================================================================================================= + +#Import all definitions from simp_signal_2016, so that these results are consistent with the expected signal calculations +signalProcessor = simp_signal_2016.SignalProcessor(mpifpi=4.*np.pi, nsigma=1.5) + +data = ak.Array([]) + +if args.tenpct: + # Search for signal in 10% data + print('Loading 10% Data') + outfilename = f'{outfilename}_10pct' + inv_mass_range = (26,200) + branches = ["unc_vtx_mass","unc_vtx_psum", "unc_vtx_ele_track_z0", "unc_vtx_pos_track_z0", "unc_vtx_z", "unc_vtx_proj_sig"] + infile = '/sdf/group/hps/user-data/alspellm/2016/data/full_hadd_blpass4c_ana.root' + selection = 'vtxana_Tight_2016_simp_reach_SR' + data = signalProcessor.load_data(infile,selection, expressions=branches) + data['weight'] = 1.0 + +else: + # Search for signal in 100% data + print('Loading 100% Data') + outfilename = f'{outfilename}_100pct' + inv_mass_range = (30,124) + branches = ["unc_vtx_mass","unc_vtx_psum", "unc_vtx_ele_track_z0", "unc_vtx_pos_track_z0", "unc_vtx_z", "unc_vtx_proj_sig"] + indir = '/fs/ddn/sdf/group/hps/users/alspellm/data_storage/pass4kf/pass4kf_ana_20240513' + mass_safety = 'unc_vtx_mass*1000. >= 0' #choose to blind certain masses if necessary + + #If high psum, can look at all masses + if args.highPsum: + selection = 'vtxana_Tight_2016_simp_reach_CR' + else: + selection = 'vtxana_Tight_2016_simp_reach_SR' + + #Loop over all 100% data vertex ana tuples + for filename in sorted(os.listdir(indir)): + if not filename.endswith('.root'): + continue + run = filename.split('_')[4] + print('Loading Run ', run) + infile = os.path.join(indir,filename) + + #build awkward array for 100% data + data = ak.concatenate([data, signalProcessor.load_data(infile, selection, cut_expression=mass_safety, expressions=branches)]) + data['weight'] = 1.0 + + +####################################################################################################################################### +#======================================================================================================================================= +# Useful Functions +#======================================================================================================================================= +poisson_low_err = lambda n : np.sqrt(n - 0.25) if n >= 0.25 else 0.0 +poisson_up_err = lambda n : np.sqrt(n+0.75) + 1 +def get_abcd_error(A, B, C, D, E): + """ + Calculate the estimated background error bars using the ABCD method. + + Args: + A,B,C,D,E (int): Number of events in each region. + + Returns: + sigma_F_up, sigma_F_low (float, float): Error bars on the estimated background in region F (signal region). + """ + + #Calculate poisson errors on regions A and E combined. + sigma_AE_up = poisson_up_err(A+E) + sigma_AE_low = poisson_low_err(A+E) + + #Calculate Gaussian errors on regions B+D and region C. + sigma_BD = math.sqrt(B+D) + sigma_C = math.sqrt(C) + + # Calculate the partial derivatives + partial_AE = C / (B + D) + partial_BD = -((A + E) * C) / ((B + D) ** 2) + partial_C = (A + E) / (B + D) + + # Calculate the propagated uncertainty + sigma_F_up = math.sqrt( + (partial_AE * sigma_AE_up) ** 2 + + (partial_BD * sigma_BD) ** 2 + + (partial_C * sigma_C) ** 2 + ) + + sigma_F_low = math.sqrt( + (partial_AE * sigma_AE_low) ** 2 + + (partial_BD * sigma_BD) ** 2 + + (partial_C * sigma_C) ** 2 + ) + return sigma_F_up, sigma_F_low + +def calc_exp_bkg(A, B, C, D, E): + """ + Calculates the expected background using the counts in regions A, B, C, D, and E. + + Args: + A, B, C, D, E (int): Number of events in each region. + + Returns: + exp_bkg (float): Expected background using ABCD estimation. + counts (list[int]): List of integers corresponding to counts in each region. + + Notes: + If the counts in A+E < 1, set A+E = 0.4. + 0.4 is the Poisson mean where 0 is observed ~68% of the time. + """ + # Set the minimum count rate in A+E + if A+E < 0.4: + A = 0.4 + E = 0.0 + + #Expected background ABCD calculation + exp_bkg = C*((A+E)/(B+D)) + counts = [A, B, C, D, E] + return exp_bkg, counts + +def run_abcd_method(data, signal_mass): + """ + Calculate the expected background, with errors, in the invariant mass window signal region by defining the regions ABCDEF. + Also return the number of events observed in the signal region. + + Args: + data (awkward array): Input data array. + signal_mass (float): Center of the invariant mass search window in dark vector meson mass. + + Returns: + expected_bkg (float): Expected background in the search window signal region estimated using ABCD method. + nobs (float): Number of events observed in the search window signal region. + bkg_error ([float,float]): [upper,lower] error bars on the expected background. + counts (list[float]): Number of events in each region except search window signal region [A, B, C, D, E]. + """ + + # Define the boundaries of the invariant mass search window + mass_low = signal_mass - signalProcessor.mass_resolution(signal_mass)*nsigma + mass_high = signal_mass + signalProcessor.mass_resolution(signal_mass)*nsigma + # Define the upper and lower boundaries of the mass sidebands. + left_xlow = signal_mass - signalProcessor.mass_resolution(signal_mass)*left_nsigma + right_xhigh = signal_mass + signalProcessor.mass_resolution(signal_mass)*right_nsigma + print(f'Left:{left_xlow}-{mass_low} | Search:{mass_low}-{mass_high} | Right:{mass_high}-{right_xhigh}') + + # Apply the invariant mass search window mass selection on the data, as well as the sidebands. + mass_sel = {} + mass_sel[f'search_window'] = (data.unc_vtx_mass * 1000. <= mass_high) & (data.unc_vtx_mass * 1000. >= mass_low) + mass_sel[f'left_sideband'] = (data.unc_vtx_mass * 1000. <= mass_low) & (data.unc_vtx_mass * 1000. >= left_xlow) + mass_sel[f'right_sideband'] = (data.unc_vtx_mass * 1000. <= right_xhigh) & (data.unc_vtx_mass * 1000. >= mass_high) + + # ABCD method extends the sidebands and search window into the minimum z0 (aka y0) parameter distribution. + # Define and fill the z0 (aka y0) distributions for each region. + min_z0_h = ( + hist.Hist.new + .StrCategory(list(mass_sel.keys()), name='mass selection') + .Reg(400,-0.005,2.995,label='Vtx Min Z0 [mm]') + .Double() + ) + min_z0_h.fill(f'search_window', data[mass_sel[f'search_window']].unc_vtx_min_z0,weight=data[mass_sel[f'search_window']].weight ) + min_z0_h.fill(f'left_sideband', data[mass_sel[f'left_sideband']].unc_vtx_min_z0, weight=data[mass_sel[f'left_sideband']].weight) + min_z0_h.fill(f'right_sideband', data[mass_sel[f'right_sideband']].unc_vtx_min_z0, weight=data[mass_sel[f'right_sideband']].weight) + # Regions A E and F are defined based on the value of the minimum z0 cut in region F + minz0_coeffs = signalProcessor.get_minz0_cut() + min_z0_cut = signalProcessor.polynomial(minz0_coeffs[0],minz0_coeffs[1],minz0_coeffs[2])(signal_mass) + + # Define the minimum z0 floor used to count events in regions B, C, and D. + # This floor is defined so that the background estimate is weighted towards the tails of the minimum z0 distributions rather + # than the core of the distribution. However, the tails cannot be so small that signal contamination is an issue in C. + xwidth = min_z0_h[f'search_window',:].axes[0].widths[0] + xmax = min_z0_cut - 2*xwidth + threshold = z0_floor_threshold + while xmax > 0.0: + sig_int = min_z0_h[f'search_window',:][hist.loc(xmax):hist.loc(min_z0_cut):sum] + if sig_int < threshold: + xmax = xmax - xwidth + else: + break + z0_floor = round(xmax,2) + + # Count the events in each region + B = min_z0_h[f'left_sideband',:][hist.loc(z0_floor):hist.loc(min_z0_cut):sum] + A = min_z0_h[f'left_sideband',:][hist.loc(min_z0_cut)+1::sum] + D = min_z0_h[f'right_sideband',:][hist.loc(z0_floor):hist.loc(min_z0_cut):sum] + E = min_z0_h[f'right_sideband',:][hist.loc(min_z0_cut)+1::sum] + C = min_z0_h[f'search_window',:][hist.loc(z0_floor):hist.loc(min_z0_cut):sum] + counts = [A, B, C, D, E] + + #Calculate the error on the background estimate according to the counts + bkg_error = get_abcd_error(A, B, C, D, E) + #Calculate the expected background + expected_bkg, counts = calc_exp_bkg(A, B, C, D, E) + + # Determine the observed number of events after applying the last selection criteria, mass and minimum z0. + minz0_sel = signalProcessor.minz0_sel(data) + masswindow_sel = signalProcessor.mass_sel(data, signal_mass) + final_sel = np.logical_and.reduce([masswindow_sel, minz0_sel]) + nobs = ak.sum(data[final_sel].weight) + + min_z0_h.reset() + return expected_bkg, nobs, bkg_error, counts + +def get_t0(A, B, C, D, E, ntrials=100000): + """ + Constructs the background-only test-statistic distribution used to calculate the significance of the data in each search window. + This test-statistic distribution has the statistical uncertainty on the background baked into it. + + Args: + A, B, C, D, E (int): Counts in each region. + ntrials (int, optional): Specify the number of toy MC trials used to construct the test-stat distribution. + This should be large enough to achieve a statistically significant right tail. + + Returns: + t0_distribution (hist): Hist histogram of the background-only test-statistic distribution. + Performing a one-sided tail integral of this histogram starting from the number of observed events + gives the local p-value of the data in the search window signal region. + """ + + # Define the background-only test-statistic distribution + t0_distribution = ( + hist.Hist.new + .Reg(500, 0.0, 500.0, label='Expected Background Toy MC Trials') + .Double() + ) + + # Sample the parent distributions of each region used to estimate the background rate + A_E_s = np.random.poisson(lam=(A+E), size=ntrials) + B_D_s = np.random.normal(loc=(B+D), scale=np.sqrt(B+D), size=ntrials) + C_s = np.random.normal(loc=C, scale=np.sqrt(C), size=ntrials) + + # Calculate the expected background (b) and t0 for all trials + b = (A_E_s / B_D_s) * C_s + + #Use the expected background b as the true mean of the background parent distribution (which is an approximation) + t0 = np.random.poisson(lam=b) + + # Fill background-only test-statistic distribution + t0_distribution.fill(t0) + + return t0_distribution + +def get_pvalue(test_stat_h, nobs): + """ + Calculate the local p-value of the data in the search window signal region. + + Args: + test_stat_h (hist histogram): Background-only test-statistic distribution returned by get_t0(). + nobs (float0: Number of events observed in the search window signal region. + + Returns: + mean (float): Local p-value mean. + low_err (float): Lower error bar on local p-value mean. + up_err (float): Upper error bar on local p-value mean. + """ + + # Count the number of events in the background-only test-statistic distribution greater than the number of events observed. + # This will be normalized to the total distribution. + try: + nover = test_stat_h[hist.loc(nobs)::sum] + except: + nover = 0.0 + + # Define hist histograms to store the total number of events in the background-only test-stat distribution, and nover. + numer_h = ( + hist.Hist.new + .Reg(1, 0.0, 1.1, label='Events past nobs') + .Double() + ) + numer_h.fill(np.ones(int(nover))) + + denom_h = ( + hist.Hist.new + .Reg(1, 0.0, 1.1, label='Total events thrown') + .Double() + ) + denom_h.fill(np.ones(int(test_stat_h[::sum]))) + test_stat_h.reset() + + # Convert the hist histograms to ROOT histograms so that they can be divided and have the correct error bars. + # Seems messy, but filling the hist histograms and then converting to ROOT is *way* faster when you have 100 million events. + # The reason we use histograms is that we get the Clopper Pearson errors when we divide two ROOT histograms. + histos = [numer_h, denom_h] + uproot_file = uproot.recreate(f'tmp_cnv_histos.root') + for i, histo in enumerate(histos): + uproot_file[f'histo_{i}'] = histo + uproot_file.close() + infile = r.TFile(f'tmp_cnv_histos.root',"READ") + for i, histo in enumerate(histos): + histos[i] = copy.deepcopy(infile.Get(f'histo_{i}')) + infile.Close() + + # The local p-value of the search window signal region is calculated by taking the ratio of nover to the total number of toys + histos[0].SetBinErrorOption(1) + histos[1].SetBinErrorOption(1) + result_g = r.TGraphAsymmErrors() + result_g.Divide(histos[0], histos[1], opt="cp") # Specify clopper pearson errors for correct error bars + + # Reset the histograms to clear mem + numer_h.reset() + denom_h.reset() + + mean = result_g.GetY()[0] + up_err = result_g.GetErrorYhigh(0) + low_err = result_g.GetErrorYlow(0) + + return mean, low_err, up_err + +#======================================================================================================================================= +# SEARCH FOR SIGNAL: Calculates the local p-value of each invariant mass search window signal region. +#======================================================================================================================================= + +#======================================================================================================================================= +# INITIALIZE +#======================================================================================================================================= + +# Define the invariant mass range to perform the search for signal +inv_masses = np.array([x for x in range(inv_mass_range[0], inv_mass_range[-1])]) + +# Save the results to these arrays +exp_bkg_mev=[] +nobs_mev=[] +bkg_uperror_mev = [] +bkg_lowerror_mev = [] +abcd_counts_mev = [] +pvalue_mev = [] +pvalue_uperr_mev = [] +pvalue_lowerr_mev = [] + +# Select the data set (10% or 100%) and the Psum region to search for signal +if not args.tenpct and args.highPsum: + psum_sel = signalProcessor.psum_sel(data, case='cr') +elif args.tenpct and not args.highPsum: + psum_sel = signalProcessor.psum_sel(data, case='sr') +elif args.tenpct and args.highPsum: + psum_sel = signalProcessor.psum_sel(data, case='cr') +else: + psum_sel = signalProcessor.psum_sel(data, case='sr') + +# Initialize subset of the pre-defined selections. These do not yet include mass selection or minimum z0 (aka y0) cut. +# Those cuts can't be applied before the background estimation is completed. +zcut_sel = signalProcessor.zcut_sel(data) # Cut on reconstructed vertex z < -4.3 mm +vprojsig_sel = signalProcessor.vprojsig_sel(data) # Cut events with target projected vertex significance > 2.0 +sameside_sel = signalProcessor.sameside_z0_cut(data) # Cut events where both tracks have same-side z0 (degenerate cut...) +initial_sel = np.logical_and.reduce([zcut_sel, vprojsig_sel, psum_sel, sameside_sel]) # Combine initial selections + +#======================================================================================================================================= +# RUN SEARCH FOR SIGNAL +#======================================================================================================================================= + +# Loop over invariant mass range +for m,mass in enumerate(inv_masses): + print(f'Running Signal Mass Window Center {mass}') + + # Estimate the background in the search window signal region using ABCD method + exp_bkg, nobs, bkg_error, counts = run_abcd_method(data[initial_sel], mass) + bkg_lowerror_mev.append(bkg_error[1]) + bkg_uperror_mev.append(bkg_error[0]) + exp_bkg_mev.append(exp_bkg) + nobs_mev.append(nobs) + abcd_counts_mev.append(counts) + print(f'background estimate: {exp_bkg} | nobs: {nobs} | counts: {counts} | bkg error: {bkg_error}') + + #Calculate the p-value by building the test statistic distribution t0 + # Calculate the significance of the data (local p-value) by constructing the background-only test-statistic distribution + # and then performing a one-tailed integral starting at the observed number of events in the search window signal region + t0_distr_h = get_t0(counts[0], counts[1], counts[2], counts[3], counts[4], ntrials=t0_trials) + pmean, plow_err, pup_err = get_pvalue(t0_distr_h, nobs) + pvalue_uperr_mev.append(pup_err) + pvalue_lowerr_mev.append(plow_err) + pvalue_mev.append(pmean) + print('local pvalue: ', pmean) + +# Convert the search results into numpy arrays +inv_masses = np.array(inv_masses, dtype=float) +exp_bkg_mev = np.array(exp_bkg_mev, dtype=float) +nobs_mev = np.array(nobs_mev, dtype=float) +bkg_uperror_mev = np.array(bkg_uperror_mev, dtype=float) +bkg_lowerror_mev = np.array(bkg_lowerror_mev, dtype=float) +pvalue_mev = np.array(pvalue_mev , dtype=float) +pvalue_uperr_mev = np.array(pvalue_uperr_mev, dtype=float) +pvalue_lowerr_mev = np.array(pvalue_lowerr_mev, dtype=float) + + +# Save the expected background with errors in a ROOT TGraphAsymmErrors +expected_bkg_g = r.TGraphAsymmErrors(len(inv_masses), inv_masses, exp_bkg_mev, np.zeros(len(inv_masses)), np.zeros(len(inv_masses)), bkg_lowerror_mev, bkg_uperror_mev) +expected_bkg_g.SetName('expected_background') +expected_bkg_g.SetTitle('Expected Background;Vd Invariant Mass [MeV]; Events') + +# Save the number of observed events in each search window signal region as a ROOT Tgraph +nobs_g = r.TGraph(len(inv_masses), inv_masses, nobs_mev) +nobs_g.SetName('Nobs') +nobs_g.SetTitle('Observed;Vd Invariant Mass [MeV]; Events') + +# Save the local p-value with errors for each search window signal region as a ROOT TGraphAsymmErrors +pvalue_g = r.TGraphAsymmErrors(len(inv_masses), inv_masses, pvalue_mev, np.zeros(len(inv_masses)), np.zeros(len(inv_masses)), pvalue_lowerr_mev, pvalue_uperr_mev) +pvalue_g.SetName('local_pvalue') +pvalue_g.SetTitle('Local P-Value;Vd Invariant Mass [MeV]; local p-value') + +# Calcualte the global p-value corrected for the "Look Elsewhere Effect" (LEE) +avg_resolution = np.average(np.array([signalProcessor.mass_resolution(x) for x in inv_masses])) +look_elsewhere = np.array((inv_masses[-1] - inv_masses[0])/avg_resolution) +print(f'Average mass resolution: {avg_resolution}') +print(f'Look elsewhere effect: {look_elsewhere}') + +# Global p-value with LEE correction +pvalue_global_g = r.TGraphAsymmErrors(len(inv_masses), inv_masses, pvalue_mev*look_elsewhere, np.zeros(len(inv_masses)), np.zeros(len(inv_masses)), pvalue_lowerr_mev*look_elsewhere, pvalue_uperr_mev*look_elsewhere) +pvalue_global_g.SetName('global_pvalue') +pvalue_global_g.SetTitle('Global P-Value;Vd Invariant Mass [MeV]; global p-value') + +# Save the search results +outfile = r.TFile(f'{outfilename}.root', "RECREATE") +outfile.cd() +expected_bkg_g.Write() +nobs_g.Write() +pvalue_g.Write() +pvalue_global_g.Write() +outfile.Close() + +# Print the local and global thresholds for reference/recording +thresholds = [] +thresholds_lew = [] +from scipy.stats import norm +for nsigma in [1,2,3,4,5]: + gaus_cdf = norm.cdf(nsigma) + threshold = (1.0 - gaus_cdf)/look_elsewhere + thresholds_lew.append(threshold) + thresholds.append((1.0 - gaus_cdf)) + +print('Local Nsigma thresholds: ', thresholds) +print('Global LE thresholds: ', thresholds_lew) + + diff --git a/plotUtils/simps/simp_plot_utils.py b/plotUtils/simps/simp_plot_utils.py new file mode 100644 index 000000000..208155460 --- /dev/null +++ b/plotUtils/simps/simp_plot_utils.py @@ -0,0 +1,597 @@ +import ROOT as r +import numpy as np +from array import array +import copy +import uproot + +def getColorsHPS(): + colors = [r.kBlue+2, r.kCyan+2, r.kRed+2, r.kOrange+10, r.kYellow+2, r.kGreen-1, r.kAzure-2, r.kGreen-8, r.kOrange+3, r.kYellow+2, r.kRed+2, r.kBlue+2, r.kGreen-8, r.kOrange+3, r.kYellow+2, r.kRed+2, r.kBlue+2, r.kGreen-8, r.kOrange+3, r.kYellow+2, r.kRed+2, r.kBlue+2, r.kGreen-8, r.kOrange+3, r.kYellow+2, r.kRed+2, r.kBlue+2, r.kGreen-8, r.kOrange+3] + return colors + +def plot_zbi_multigraph(graph_nsig, graph_nbkg, graph_zbi, canvas_name, leftyaxis_title, rightyaxis_title): + # Create TMultiGraph to hold all TGraphs + multi_graph = r.TMultiGraph() + + # Add nsig and nbkg graphs to the multi_graph (use left y-axis) + multi_graph.Add(graph_nsig) + multi_graph.Add(graph_nbkg) + + # Create a TCanvas to draw the multi_graph + canvas = r.TCanvas(canvas_name, canvas_name, 2040, 1080) + + # Draw the multi_graph with left y-axis (default) + multi_graph.Draw("AP") + + c.Draw() + # Set axis labels and titles for the left y-axis + multi_graph.GetYaxis().SetTitle(leftyaxis_title) + multi_graph.GetYaxis().SetTitleOffset(1.2) + + # Create a second y-axis for the ZBi graph + axis_zbi = multi_graph.GetHistogram().GetXaxis() # Copy the x-axis + axis_zbi.SetTitle(rightyaxis_title) # Set title for the right y-axis + axis_zbi.SetTitleOffset(1.2) # Adjust title offset for the right y-axis + axis_zbi.SetLabelOffset(999) # Hide labels for the right y-axis + + # Draw the ZBi graph with the right y-axis + graph_zbi.SetLineColor(r.kRed) # Set color for better visibility + graph_zbi.Draw("same") + graph_zbi.SetMarkerColor(r.kRed) # Set marker color + + # Set the second y-axis for the ZBi graph + multi_graph.GetHistogram().GetListOfFunctions().Add(axis_zbi) + + # Update the canvas + canvas.Update() + + # Return the canvas + return canvas + +#def SetMyStyle(tsize=0.025, tzsize=0.025, font=42, setOptTitle=0, setOptStat=0, setOptFit=0): +def SetMyStyle(tsize=0.05, tzsize=0.05, font=42, setOptTitle=0, setOptStat=0, setOptFit=0): + print("SETTING MY STYLE") + + colors = getColorsHPS() + r.gROOT.SetBatch(1) + + myStyle = r.TStyle("myStyle", "my style") + + # Set your custom attributes here + myStyle.SetOptTitle(setOptTitle) + myStyle.SetOptStat(setOptStat) + myStyle.SetOptFit(setOptFit) + myStyle.SetTitleFont(font) + myStyle.SetTitleSize(tsize) + #myStyle.SetTitleX(0.5) + #myStyle.SetTitleY(0.98) + + #Set legend text size + myStyle.SetLegendTextSize(0.02) + + # Set the title text color to black + myStyle.SetTitleTextColor(r.kBlack) + + # use plain black on white colors + icol = 0 + myStyle.SetFrameBorderMode(icol) + myStyle.SetCanvasBorderMode(icol) + myStyle.SetPadBorderMode(icol) + myStyle.SetPadColor(icol) + myStyle.SetCanvasColor(icol) + myStyle.SetStatColor(icol) + + # set the paper & margin sizes + myStyle.SetPaperSize(20, 26) + myStyle.SetPadTopMargin(0.10) + myStyle.SetPadRightMargin(0.05) + myStyle.SetPadRightMargin(0.10) + myStyle.SetPadBottomMargin(0.15) + myStyle.SetPadLeftMargin(0.10) + + myStyle.SetTextSize(tsize) + myStyle.SetLabelFont(font, "x") + myStyle.SetTitleFont(font, "x") + myStyle.SetLabelFont(font, "y") + myStyle.SetTitleFont(font, "y") + myStyle.SetLabelFont(font, "z") + myStyle.SetTitleFont(font, "z") + + myStyle.SetLabelSize(tsize, "x") + myStyle.SetTitleSize(tsize, "x") + myStyle.SetLabelSize(tsize, "y") + myStyle.SetTitleSize(tsize, "y") + myStyle.SetLabelSize(tzsize, "z") + myStyle.SetTitleSize(tzsize, "z") + + myStyle.SetTitleOffset(1.25, "y") + myStyle.SetTitleOffset(1.5, "x") + + #use bold lines and markers + myStyle.SetMarkerSize(1.0) + myStyle.SetMarkerStyle(8) + myStyle.SetMarkerColor(1) + myStyle.SetLineColor(1) + myStyle.SetHistLineWidth(3) + #myStyle.SetLineStyleString(2, "[12 12]") # postscript dashes + + # put tick marks on top and RHS of plots + #myStyle.SetPadTickX(1) + #myStyle.SetPadTickY(1) + + r.gROOT.SetStyle("myStyle") + r.gROOT.ForceStyle() + + NRGBs = 5 + NCont = 255 + + stops = array("d", [0.00, 0.34, 0.61, 0.84, 1.00]) + red = array("d", [0.00, 0.00, 0.87, 1.00, 0.51]) + green = array("d", [0.00, 0.81, 1.00, 0.20, 0.00]) + blue = array("d", [0.51, 1.00, 0.12, 0.00, 0.00]) + r.TColor.CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont) + + return myStyle + + + +def buildLegend(graphs,titles=[], position=(0.50,0.6,0.85,0.9),clear_legend=True, text_size=0.030, entry_format=None): + legend = r.TLegend(*position) + #legend = r.TLegend() + # Set the legend to transparent (clear) if the option is specified + if clear_legend: + legend.SetFillStyle(0) + legend.SetFillColor(0) + legend.SetLineColor(0) + legend.SetBorderSize(0) + legend.SetTextSize(text_size) + + for i,graph in enumerate(graphs): + if entry_format is None: + if len(titles) < 1: + legend.AddEntry(graph, graph.GetTitle()) + else: + legend.AddEntry(graph, titles[i]) + else: + if len(titles) < 1: + legend.AddEntry(graph, graph.GetTitle(),"%s"%(entry_format[i])) + else: + legend.AddEntry(graph, titles[i],"%s"%(entry_format[i])) + return legend + + +def setTitle(hist, title=None, titlex=None, titley=None, showTitle=False): + if title is not None: + hist.SetTitle(title) + if titlex is not None: + hist.GetXaxis().SetTitle(titlex) + if titley is not None: + hist.GetYaxis().SetTitle(titley) + if showTitle: + r.gStyle.SetOptTitle(1) + else: + r.gStyle.SetOptTitle(0) + +def drawStatsBox(hist, statsPos=[0.7,0.9,0.7,0.9]): + statsbox = hist.FindObject("stats") + statsbox.SetX1NDC(statsPos[0]) + statsbox.SetX2NDC(statsPos[1]) + statsbox.SetY1NDC(statsPos[2]) + statsbox.SetY2NDC(statsPos[3]) + +def drawTH2(hist, canvas_name, drawOpt='colz', xrange = (None,None), yrange=(None,None),rebinx=None, + rebiny=None, size=(2200,1600), logX=False, logY=False, logZ=False, save=False, outdir='.', + drawStats=False, statsPos=[0.7,0.9,0.7,0.9], + text=[], text_pos = [0.2, 0.8], line_spacing=0.03, text_size=0.025, Hps=True): + + c = r.TCanvas(f'{canvas_name}',f'{canvas_name}',size[0], size[1]) + c.cd() + r.gROOT.ForceStyle() + hist.Draw(f'{drawOpt}') + r.gPad.Update() + c.UseCurrentStyle() + if rebinx is not None: + hist.RebinX(rebinx) + if rebiny is not None: + hist.RebinY(rebiny) + + if drawStats: + drawStatsBox(hist, statsPos) + else: + hist.SetStats(0) + + #Set Xrange + xmin = hist.GetXaxis().GetXmin() + xmax = hist.GetXaxis().GetXmax() + if xrange[0] is not None: + xmin = xrange[0] + if xrange[1] is not None: + xmax = xrange[1] + hist.GetXaxis().SetRangeUser(0.9*xmin,1.1*xmax) + + #Set Yrange + ymin = hist.GetYaxis().GetXmin() + ymax = hist.GetYaxis().GetXmax() + if yrange[0] is not None: + ymin = yrange[0] + if yrange[1] is not None: + ymax = yrange[1] + hist.GetYaxis().SetRangeUser(0.9*ymin,1.1*ymax) + + #Set Log + if logX: + c.SetLogx(1) + else: + c.SetLogx(0) + if logY: + c.SetLogy(1) + else: + c.SetLogy(0) + if logZ: + c.SetLogz(1) + else: + c.SetLogz(0) + + #Draw latex text + if len(text) > 0 or Hps: + drawText = insertText(text, text_pos, line_spacing, text_size, Hps) + drawText.Draw() + r.gPad.Update() + + if save: + c.SaveAs(f'{outdir}/{canvas_name}.png') + c.Close() + else: + return c + +def insertText(text=[], text_pos = [0.2, 0.98], line_spacing=0.04, text_size=0.035, Hps=True): + latex = r.TLatex() + latex.SetTextSize(text_size) + text_x = text_pos[0] + text_y = text_pos[1] + if (Hps): + latex.DrawLatexNDC(text_x, text_y,'#bf{#it{HPS}} Internal') + text_y = text_y - line_spacing + + for line in text: + latex.DrawLatexNDC(text_x, text_y,line) + text_y = text_y - line_spacing + return latex + +def drawTH1s(histograms, canvas_name, drawOpts=[], xrange = (None,None), yrange=(None,None),rebinx=None, legend=None, + rebiny=None, size=(1280,720), logX=False, logY=False, logZ=False, save=False, outdir='.', + drawStats=False, statsPos=[0.7,0.9,0.7,0.9], text=[], text_pos = [0.15, 0.8], line_spacing=0.03, text_size=0.025, Hps=True, freezeXaxis=True): + + c = r.TCanvas(f'{canvas_name}',f'{canvas_name}',size[0], size[1]) + c.cd() + r.gROOT.ForceStyle() + c.UseCurrentStyle() + + # Find the maximum x and y values among all histograms + min_x = min(h.GetBinLowEdge(h.FindFirstBinAbove(0.0)) for h in histograms) + max_x, max_y = max(h.GetBinLowEdge(h.FindLastBinAbove(0.0)) + h.GetBinWidth(0) for h in histograms), max(h.GetMaximum() for h in histograms) + min_y = 1e10 + for h in histograms: + local_miny = 1e10 + min_y_l = h.GetBinContent(h.FindFirstBinAbove(0.0)) + min_y_u = h.GetBinContent(h.FindLastBinAbove(0.0)) + if min_y_l < min_y_u: + local_miny = min_y_l + else: + local_miny = min_y_u + if local_miny < min_y: + min_y = local_miny + if(freezeXaxis == False): + h.SetAxisRange(0.95*min_x, 1.05*max_x, "X") + h.GetXaxis().SetRangeUser(0.95*min_x,1.05*max_x) + h.SetAxisRange(min_y, 1.05 * max_y, "Y") + h.GetYaxis().SetRangeUser(min_y, 1.05 * max_y) + + for i, gr in enumerate(histograms): + if xrange[0] is not None and xrange[0] is not None: + gr.GetXaxis().SetRangeUser(0.95*xrange[0], 1.05*xrange[1]) + if yrange[0] is not None: + min_y = yrange[0] + if yrange[1] is not None: + max_y = yrange[1] + if min_y < 0: + gr.GetYaxis().SetRangeUser(0.0, 1.05*max_y) + else: + gr.GetYaxis().SetRangeUser(0.95*min_y, 1.05*max_y) + + #Draw histograms + r.gPad.Update() + if i < 1: + gr.Draw('%s'%(drawOpts[i])) + gr.SetLineColor(colors[i]) + r.gPad.Update() + if drawStats: + drawStatsBox(gr, statsPos) + else: + gr.SetStats(0) + else: + gr.Draw('%sSAME'%(drawOpts[i])) + gr.SetLineColor(colors[i]) + r.gPad.Update() + if drawStats: + drawStatsBox(gr, statsPos) + else: + gr.SetStats(0) + + #Draw fit functions if present + if len(gr.GetListOfFunctions()) > 0: + func_name = gr.GetListOfFunctions().At(0).GetName() + func = gr.GetFunction("%s"%(func_name)) + func.Draw("%sSAME"%(drawOpts[i])) + func.SetLineColor(r.kRed) + + #Draw legend + if legend is not None: + legend.Draw() + c.Update() + + #Set Log + if logX: + c.SetLogx(1) + else: + c.SetLogx(0) + if logY: + c.SetLogy(1) + else: + c.SetLogy(0) + if logZ: + c.SetLogz(1) + else: + c.SetLogz(0) + + #Draw latex text + if len(text) > 0 or Hps: + drawText = insertText(text, text_pos, line_spacing, text_size, Hps) + drawText.Draw() + c.Update() + + c.Draw() + if save: + c.SaveAs(f'{outdir}/{canvas_name}.png') + c.Close() + else: + return c + +def drawTGraphs(graphs, canvas_name, drawOpts=[], xrange = (None,None), yrange=(None,None),rebinx=None, legend=None, + rebiny=None, size=(2200,1600), logX=False, logY=False, logZ=False, save=False, outdir='.', + drawStats=False, statsPos=[0.7,0.9,0.7,0.9], text=[], text_pos = [0.15, 0.8], line_spacing=0.03, text_size=0.025, Hps=True): + + c = r.TCanvas(f'{canvas_name}',f'{canvas_name}',size[0], size[1]) + c.cd() + r.gROOT.ForceStyle() + c.UseCurrentStyle() + + #Set Range + ymin = 9999.9 + ymax = -9999.9 + for i, gr in enumerate(graphs): + num_points = gr.GetN() + y_values = np.array([gr.GetY()[i] for i in range(num_points)]) + local_ymax = np.max(y_values) + if local_ymax > ymax: + ymax = local_ymax + local_ymin = np.min(y_values) + if local_ymin < ymin: + ymin = local_ymin + + for i, gr in enumerate(graphs): + if xrange[0] is not None and xrange[0] is not None: + gr.GetHistogram().GetXaxis().SetRangeUser(0.9*xrange[0], 1.1*xrange[1]) + if yrange[0] is not None: + ymin = yrange[0] + if yrange[1] is not None: + ymax = yrange[1] + if ymin < 0: + gr.GetHistogram().GetYaxis().SetRangeUser(1.2*ymin, 1.2*ymax) + else: + gr.GetHistogram().GetYaxis().SetRangeUser(0.8*ymin, 1.2*ymax) + + #Draw Graphs + if i < 1: + gr.Draw('A%s'%(drawOpts[i])) + r.gPad.Update() + if drawStats: + drawStatsBox(gr, statsPos) + else: + gr.GetHistogram().SetStats(0) + else: + gr.Draw('%sSAME'%(drawOpts[i])) + r.gPad.Update() + if drawStats: + drawStatsBox(gr, statsPos) + else: + gr.GetHistogram().SetStats(0) + + #Draw fit functions if present + if len(gr.GetListOfFunctions()) > 0: + func_name = gr.GetListOfFunctions().At(0).GetName() + func = gr.GetFunction("%s"%(func_name)) + func.Draw("%sSAME"%(drawOpts[i])) + + #Draw legend + if legend is not None: + legend.Draw() + c.Update() + + #Set Log + if logX: + c.SetLogx(1) + else: + c.SetLogx(0) + if logY: + c.SetLogy(1) + else: + c.SetLogy(0) + if logZ: + c.SetLogz(1) + else: + c.SetLogz(0) + + #Draw latex text + if len(text) > 0 or Hps: + drawText = insertText(text, text_pos, line_spacing, text_size, Hps) + drawText.Draw() + c.Update() + + c.Draw() + if save: + c.SaveAs(f'{outdir}/{canvas_name}.png') + c.Close() + else: + return c + +def readROOTHisto(infilename, histoname): + infile = r.TFile(f'{infilename}',"READ") + histo = copy.deepcopy(infile.Get(f'{histoname}')) + infile.Close() + return histo + +def cnvHistosToROOT(histos=[], tempname='temporary_uproot'): + return_histos = [] + uproot_file = uproot.recreate(f'trash_{tempname}.root') + for i, histo in enumerate(histos): + uproot_file[f'histo_{i}'] = histo + uproot_file.close() + infile = r.TFile(f'trash_{tempname}.root',"READ") + for i, histo in enumerate(histos): + return_histos.append(copy.deepcopy(infile.Get(f'histo_{i}'))) + infile.Close() + return return_histos + +def cnvHistogramToROOT(histo, tempname='temporary_uproot'): + uproot_file = uproot.recreate(f'trash_{tempname}.root') + uproot_file['histogram'] = histo + root_hist = readROOTHisto(f'trash_{tempname}.root', 'histogram;1') + uproot_file.close() + return root_hist + +def cnvHistoToROOT(histoname, histo, tempname='temporary_uproot'): + uproot_file = uproot.recreate(f'trash_{tempname}.root') + uproot_file[histoname] = histo + infile = r.TFile(f'trash_{tempname}.root',"READ") + root_hist = copy.deepcopy(infile.Get(f'{histoname}')) + infile.Close() + uproot_file.close() + return root_hist + +def quickDraw(plot, name='c', drawOpts=""): + c = r.TCanvas(name,name,1400,700) + c.cd() + plot.Draw(drawOpts) + c.Draw() + return c + +def format_th1(histo, name=None, title=None, xlabel=None, ylabel=None, linecolor=None, linewidth=None, markerstyle=None, + markercolor=None, xrang=(None, None), yrang=(None,None), rebin=None, labelsize=None, + titlesize=None, xtitleoffset=None, ytitleoffset=None): + if name: + histo.SetName(name) + if title: + histo.SetTitle(title) + if xlabel: + histo.GetXaxis().SetTitle(xlabel) + if ylabel: + histo.GetYaxis().SetTitle(ylabel) + if linecolor: + histo.SetLineColor(linecolor) + if linewidth: + histo.SetLineWidth(linewidth) + if markerstyle: + histo.SetMarkerStyle(markerstyle) + if markercolor: + histo.SetMarkerColor(markercolor) + if xrang[0]: + histo.GetXaxis().SetRangeUser(xrang[0],xrang[1]) + if yrang[0]: + histo.GetYaxis().SetRangeUser(yrang[0], yrang[1]) + if rebin: + histo.Rebin(rebin) + if labelsize: + histo.GetXaxis().SetLabelSize(labelsize) + histo.GetYaxis().SetLabelSize(labelsize) + if titlesize: + histo.GetXaxis().SetTitleSize(titlesize) + histo.GetYaxis().SetTitleSize(titlesize) + if xtitleoffset: + histo.GetXaxis().SetTitleOffset(xtitleoffset) + if ytitleoffset: + histo.GetYaxis().SetTitleOffset(ytitleoffset) + +def plot_ratio_th1s(top_plots, bot_plots, cname, size=(1800,1400), top_drawOpts="pe", bot_drawOpts="pe", + topx1=0.0, topy1=0.4, topx2=1.0, topy2=1.0, top_logY=False, + botx1=0.0, boty1=0.0, botx2=1.0, boty2=0.4, bot_logY=False, + top_botmarg=0.1, top_topmarg=None, bot_topmarg=None, bot_botmarg=0.1, leftmarg=0.1): + + c = r.TCanvas(cname, cname, size[0], size[1]) + top = r.TPad("%s_top"%(cname), "%s_top"%(cname), topx1, topy1, topx2, topy2) + top.SetBottomMargin(top_botmarg) + top.SetLeftMargin(leftmarg) + if top_topmarg: + top.SetTopMargin(top_topmarg) + top.Draw() + + + + bot = r.TPad("%s_bot"%(cname), "%s_bot"%(cname), botx1, boty1, botx2, boty2) + bot.SetBottomMargin(bot_botmarg) + bot.SetLeftMargin(leftmarg) + if bot_topmarg: + bot.SetTopMargin(bot_topmarg) + bot.Draw() + + top.cd() + for i,plot in enumerate(top_plots): + if i < 1: + plot.Draw(top_drawOpts) + else: + plot.Draw("%ssame"%(top_drawOpts)) + + r.gPad.SetLogy(top_logY) + + + bot.cd() + for i,plot in enumerate(bot_plots): + if i < 1: + plot.Draw(bot_drawOpts) + + else: + plot.Draw("%ssame"%(bot_drawOpts)) + + #Draw fit functions if present + if len(plot.GetListOfFunctions()) > 0: + func_name = plot.GetListOfFunctions().At(0).GetName() + func = plot.GetFunction("%s"%(func_name)) + func.Draw("%ssame"%(bot_drawOpts)) + + return c, top, bot + + + +def draw_th1(histo, cname, size=(1800,1400), drawOpts="", logY=False, color=None): + + c = r.TCanvas(cname, cname, size[0], size[1]) + c.cd() + if color: + histo.SetLineColor(color) + histo.Draw("%s"%(drawOpts)) + + #Draw fit functions if present + if len(histo.GetListOfFunctions()) > 0: + func_name = histo.GetListOfFunctions().At(0).GetName() + func = histo.GetFunction("%s"%(func_name)) + if color: + func.SetLineColor(color) + func.Draw("same") + + if logY: + c.SetLogy(1) + + return c + +style = SetMyStyle(setOptStat=0) +colors = getColorsHPS() + diff --git a/plotUtils/simps/simp_signal_2016.py b/plotUtils/simps/simp_signal_2016.py new file mode 100644 index 000000000..955d10584 --- /dev/null +++ b/plotUtils/simps/simp_signal_2016.py @@ -0,0 +1,905 @@ +#!/usr/bin/python3 +#======================================================================================================================================= +""" +SignalProcessor Class +--------------------- +This script defines the 'SignalProcessor' class, which handles the expected signal calculations used in the 2016 SIMP L1L1 analysis. +This script also imports the necessary simp equations to perform the calculations. +This processor operates on the output of the hpstr vertex analysis processor, which is a flat tuple of events. + +Modules: + - os + - awkward as ak + - numpy as np + - hist + - uproot + - ROOT as r + - argparse + - simp_theory_equations (imported as SimpEquations) +""" + +import os +import awkward as ak +import numpy as np +import hist +from hist import Hist +import uproot +import math +import ROOT as r +import copy +import argparse +from simp_theory_equations import SimpEquations as simpeqs +#======================================================================================================================================= + +class SignalProcessor: + """ + A class for handling signal data processing and physics parameter calculations from SIMPs. + + Attributes: + nsigma (float): Size of the invariant mass search window in terms of the mass resolution (default is 1.5). + alpha_dark (float): The dark sector fine structure constant. + mass_ratio_ap_to_vd (float): Ratio of dark photon mass to dark vector meson mass. + mass_ratio_ap_to_pid (float): Ratio of dark photon mass to pion mass. + mass_lepton (float): ele mass (in MeV). + target_pos (float): Target position in mm (default is -4.3 mm) + mpifpi (float): ratio of dark pion mass and dark pion decay constant, default value i 4*pi + mass_resolution (function): Polynomial function representing mass resolution. + radiative_fraction (function): Polynomial function for radiative fraction. + radiative_acceptance (function): Polynomial function for radiative acceptance. + psum_reweighting (function): Polynomial function for reweighting based on MC bkg and data Psum. + minz0_cut_poly (function): Polynomial function for minimum z0 cut. + cr_psum_low (float): Lower bound for CR psum. + cr_psum_high (float): Upper bound for CR psum. + sr_psum_low (float): Lower bound for SR psum. + sr_psum_high (float): Upper bound for SR psum. + trident_differential_production: Holds the differential trident production rate + """ + + def __init__(self, mpifpi=4*np.pi, nsigma=1.5): + """ + Initializes the SignalProcessor with default SIMP model parameters and physics constants. + + Args: + mpifpi (float, optional): benchmark values are 3 and 4pi. + nsigma (float, optional): Size of the invariant mass search window in terms of the mass resolution (default is 1.5). + """ + + #search window size + self.nsigma = nsigma + + #SIMP parameters + self.alpha_dark = 0.01 + self.mass_ratio_ap_to_vd = 1.66 + self.mass_ratio_ap_to_pid = 3.0 + self.mass_lepton = 0.511 + self.target_pos = -4.3 + self.mpifpi = mpifpi + + #fit functions calculated externally + self.mass_resolution = self.polynomial(.75739851, 0.031621002, 5.2949672e-05) + self.radiative_fraction = self.polynomial(0.10541434, -0.0011737697, 7.4487930e-06, -1.6766332e-08) + self.radiative_acceptance = self.polynomial(-0.48922505, 0.073733061, -0.0043873158, 0.00013455495, -2.3630535e-06, 2.5402516e-08, -1.7090900e-10, 7.0355585e-13, -1.6215982e-15, 1.6032317e-18) + self.psum_reweighting = self.polynomial(0.094272950, 0.87334446, -0.19641796) #GeV argument + self.minz0_cut_poly = self.polynomial(1.07620094e+00 + 0.1, -7.44533811e-03, 1.58745903e-05) + + #signal and control region boundaries + self.cr_psum_low = 1.9 + self.cr_psum_high = 2.4 + self.sr_psum_low = 1.0 + self.sr_psum_high = 1.9 + + #scales the expected signal to data + self.trident_differential_production = None + + + def set_radiative_acceptance(self, *coefficients): + """ + Sets the polynomial coefficients for the radiative acceptance function. + Used for systematic studies where the radiative acceptance changes from nominal. + + Args: + *coefficients (float): Coefficients of the polynomial function. + """ + self.radiative_acceptance = self.polynomial(*coefficients) + + @staticmethod + def polynomial(*coefficients): + def _implementation(x): + return sum([ + coefficient * x**power + for power, coefficient in enumerate(coefficients) + ]) + return _implementation + + def load_data(self, filepath, selection, cut_expression=None, expressions=None): + """ + Loads data from hpstr vertex ana processor output file using Uproot, applying selection and cuts if provided. + + Args: + filepath (str): Path to the ROOT file. + selection (str): The dataset selection (subdir representing hpstr lvl selection) + cut_expression (str, optional): A cut expression to filter data. + expressions (list, optional): List of specific expressions (branches) to load. + + Returns: + awkward.Array: An array of selected data. + """ + with uproot.open(filepath) as f: + events = f[f'{selection}/{selection}_tree'].arrays( + cut=cut_expression, + expressions=expressions + ) + try: + events['unc_vtx_min_z0'] = np.minimum(abs(events.unc_vtx_ele_track_z0), abs(events.unc_vtx_pos_track_z0)) + except: + pass + return events + + @staticmethod + def safe_divide(numerator, denominator, default=0.0): + result = np.full(numerator.shape, default) + result[denominator > 0] = numerator[denominator > 0] / denominator[denominator > 0] + return result + + def load_pre_readout_signal_z_distribution(self, filepath): + """ + Loads MC signal vertex z distribution output from the hpstr mcana processor + + Args: + filepath (str): Path to ROOT file. + """ + with uproot.open(filepath) as sig_f: + return sig_f['mcAna/mcAna_mc625Z_h'].to_hist() + + def load_signal(self, filepath, pre_readout_filepath, mass_vd, selection, cut_expression=None, branches=None): + """ + Loads MC signal from hpstr vertex ana processor output. + + Args: + filepath (str): Path to ROOT file. + pre_readout_filepath (str): mc ana file corresponding to signal mass. + mass_vd (float): dark vector mass. + selection (str): subdir where ttree is located in ROOT file. + cut_expression (str, optional): A cut expression to filter data (Ex: "unc_vtx_psum < 1.9") + branches (str, optional): List of branches to load (default ones necessary for analysis) + + """ + with uproot.open(filepath) as sig_f: + if branches: + events = sig_f[f'{selection}/{selection}_tree'].arrays( + expressions=branches, + cut=cut_expression + ) + else: + events = sig_f[f'{selection}/{selection}_tree'].arrays( + expressions=[ + 'vd_true_vtx_z', 'unc_vtx_z', 'unc_vtx_ele_track_z0', + 'unc_vtx_pos_track_z0', 'unc_vtx_mass', 'unc_vtx_deltaZ', + 'unc_vtx_proj_sig', 'vd_true_vtx_energy', 'unc_vtx_psum' + ], + cut=cut_expression + ) + events['vd_true_gamma'] = events.vd_true_vtx_energy * 1000 / mass_vd + events['unc_vtx_min_z0'] = np.minimum(abs(events.unc_vtx_ele_track_z0), abs(events.unc_vtx_pos_track_z0)) + events['psum_reweight'] = self.psum_reweighting(events.unc_vtx_psum) + events['psum_reweight'] = ak.where(events['psum_reweight'] > 1., 1., events['psum_reweight']) + + #load the mc truth signal vertex distribution. Used to calculate signal acceptance*eff as F(z) + not_rebinned_pre_readout_z_h = self.load_pre_readout_signal_z_distribution(pre_readout_filepath) + + def sample_pre_readout_probability(z): + """ + Calculates the signal acceptance*efficiency as a function of F(z) + + Args: + z (float): truth vertex z bin. + """ + if z < not_rebinned_pre_readout_z_h.axes[0].edges[0]: + return 0. + if z > not_rebinned_pre_readout_z_h.axes[0].edges[-1]: + return 0. + index = not_rebinned_pre_readout_z_h.axes.index(z) + return not_rebinned_pre_readout_z_h.axes.widths[0][index] / not_rebinned_pre_readout_z_h[index].value + + events['event_weight_by_uniform_z'] = [ + sample_pre_readout_probability(z) + for z in events.vd_true_vtx_z + ] + return events + + def _load_trident_differential_production_lut(self, background_file, selection, signal_mass_range, mass_window_width, tenpct=True, full_lumi_path=None): + """ + Loads lookup table that stores the differential radiative trident rate as a function of A' mass. + This rate scales the expected signal to the data set. + + Args: + background_file (str): hpstr vertex ana ROOT file containing reconstructed+selected background events in CR. + selection (str): subdir where ttree is located in ROOT file. + signal_mass_range (tuple): range of dark vector masses that expected signal will be calculated for. + mass_window_width (float): width of the mass window used to calculate the reconstructed bkg rate (it gets averaged) + tenpct (bool, optional): If True, use 10% data to normalize signal. If False, use 100% data (must provide full lumi path). + full_lumi_path (str, optional): If tenpct=False, provide path to 100% data. + + Returns: + The look-up table for radiative trident differential production rates as function of A' mass. + + """ + dNdm_by_mass_vd = {} + bkgd_CR = ak.Array([]) + + if tenpct: + with uproot.open(background_file) as bkgd_f: + bkgd_CR = bkgd_f[f'{selection}/{selection}_tree'].arrays( + cut=f'( (unc_vtx_psum > {self.cr_psum_low}) & (unc_vtx_psum < {self.cr_psum_high}) )', + expressions=['unc_vtx_mass', 'unc_vtx_z'], + ) + else: + for filename in sorted(os.listdir(full_lumi_path)): + if not filename.endswith('.root'): + continue + run = filename.split('_')[4] + print('Loading Run ', run) + + + background_file = os.path.join(full_lumi_path,filename) + with uproot.open(background_file) as bkgd_f: + bkgd_CR_per_run = bkgd_f[f'{selection}/{selection}_tree'].arrays( + cut=f'( (unc_vtx_psum > {self.cr_psum_low}) & (unc_vtx_psum < {self.cr_psum_high}) )', + expressions=['unc_vtx_mass', 'unc_vtx_z'], + ) + bkgd_CR = ak.concatenate([bkgd_CR, bkgd_CR_per_run]) + + for mass_vd in signal_mass_range: + window_half_width = mass_window_width * self.mass_resolution(mass_vd) / 2. + dNdm_by_mass_vd[mass_vd] = ak.sum( + (bkgd_CR.unc_vtx_mass * 1000 > self.mass_ratio_ap_to_vd * (mass_vd - window_half_width)) & + (bkgd_CR.unc_vtx_mass * 1000 < self.mass_ratio_ap_to_vd * (mass_vd + window_half_width)) + ) / (2 * window_half_width * self.mass_ratio_ap_to_vd) + return dNdm_by_mass_vd + + def trident_differential_production(self, mass_vd): + """ + Returns the radiative trident differential production rate. + + Args: + mass_vd (float): dark vector mass + + Notes: + If the dark vector mass isn't found, could be the result of the mvd -> map conversion. + """ + if int(mass_vd) in self.trident_differential_production.keys(): + return self.trident_differential_production[mass_vd] + raise ValueError(f'The dark vector mass {mass_vd} is not found in the trident differential production look-up table.') + + def set_diff_prod_lut(self,infile, preselection, signal_mass_range, tenpct=True, full_lumi_path=None): + """ + Initializes the trident differential production lookup table. + + Args: + infile (str): hpstr vertex ana ROOT file containing reconstructed+selected background events in CR. + preselection (str): ROOT file subdir containing preselection ttree. + signal_mass_range (tuple): dark vector meson mass range + tenpct (bool, optional): If True use 10% data sample. If False, use 100% data (requires full_lumi_path). + full_lumi_path (str, optional): If tenpct == False, provide path to 100% data. + + Notes: + This was all developed using 10% data stored in a single ROOT file, hence the tenpct option. But you can pass any + single ROOT file with tenpct set to True. + The final analysis uses 100% data, which consists of multiple separate root files located at full_lumi_path. + """ + + #Initialize the lookup table to calculate the expected signal scale factor + self.trident_differential_production = self._load_trident_differential_production_lut(infile, preselection, signal_mass_range, 2.0*self.nsigma, tenpct=tenpct, full_lumi_path=full_lumi_path) + + def total_signal_production_per_epsilon2(self, signal_mass): + """ + Calculates the total dark vector meson production rate as a function of epsilon^2 + + Args: + signal_mass (float): dark vector meson mass. + + Notes: + Notice that you pass the dark vector meson mass, and it gets converted to an A' mass. + The radiative fraction, acceptance, and trident prod are all calculated using the A' mass, NOT the vector mass! + The dark vector to A' mass ratio is set upon initializing the instance of this class. + """ + mass_ap = self.mass_ratio_ap_to_vd*signal_mass + return ( + (3. * (137. / 2.) * np.pi) + * mass_ap * self.radiative_fraction(mass_ap) + * self.trident_differential_production[(int((mass_ap / self.mass_ratio_ap_to_vd)))] + / self.radiative_acceptance(mass_ap) + ) + + def get_exp_sig_eps2(self, signal_mass, signal_array, eps2): + """ + Calculates the neutral dark vector meson signal acceptance*efficiency*probability for each MC signal event for value eps^2. + This function reweights the MC signal events by considering the dark vector (rho and phi) decay probabilities. + + Args: + signal_mass (float): The generated dark vector meson mass. + signal_array (awkward array): MC signal array loaded by load_signal. + eps2 (float): The square of the kineitc mixing strength parameter (affects decay probabilities). + + Returns: + signal_array (awkward array): The input signal array with added columns for the reweighted signal + acceptance*efficiency*probability ('reweighted_accxEff'). + + Notes: + The vector decay probabilities depend on the A' mass, dark pion mass, and mpi/fpi, all of which are initialized + with the instance of this class. + + """ + + #Define SIMP masses based on mass ratios and constants + mass_ap = self.mass_ratio_ap_to_vd*signal_mass # A' mass from Vd mass + mass_pid = mass_ap / self.mass_ratio_ap_to_pid # Dark pion mass from A' mass + fpid = mass_pid / self.mpifpi # Dark pion decay constant from ratio of dark pion mass to dark pion decay constant + + #Calculate the decay length in the lab frame for rho and phi + rho_gctau = signal_array.vd_true_gamma * simpeqs.getCtau(mass_ap, mass_pid, signal_mass, np.sqrt(eps2), self.alpha_dark, fpid, self.mass_lepton, True) + phi_gctau = signal_array.vd_true_gamma * simpeqs.getCtau(mass_ap, mass_pid, signal_mass, np.sqrt(eps2), self.alpha_dark, fpid, self.mass_lepton, False) + + #Calculate the decay weights for rho and phi based on the vertex z position and gammactau + rho_decay_weight = (np.exp((self.target_pos - signal_array.vd_true_vtx_z) / rho_gctau) / rho_gctau) + phi_decay_weight = (np.exp((self.target_pos - signal_array.vd_true_vtx_z) / phi_gctau) / phi_gctau) + #signal_array['reweighted_accxEff_rho'] = simpeqs.br_Vrho_pi(mass_ap, mass_pid, signal_mass, self.alpha_dark, fpid)*rho_decay_weight*signal_array.event_weight_by_uniform_z*signal_array.psum_reweight + #signal_array['reweighted_accxEff_phi'] = simpeqs.br_Vphi_pi(mass_ap, mass_pid, signal_mass, self.alpha_dark, fpid)*phi_decay_weight*signal_array.event_weight_by_uniform_z*signal_array.psum_reweight + + #Calculate the combined decay weight for both rho and phi mesons. + #This result represenets the overall expected signal decay. + combined_decay_weight = ( + (rho_decay_weight * simpeqs.br_Vrho_pi(mass_ap, mass_pid, signal_mass, self.alpha_dark, fpid)) + + (phi_decay_weight * simpeqs.br_Vphi_pi(mass_ap, mass_pid, signal_mass, self.alpha_dark, fpid)) + ) + #The final reweighting includes both rho and phi, and includes the signal acceptance*efficiency as a function of z via + #signal_array.event_weight_by_uniform_z. Psum re-weighting (calculated externally, saved in this classs) is also included. + signal_array['reweighted_accxEff'] = combined_decay_weight*signal_array.event_weight_by_uniform_z*signal_array.psum_reweight + + return signal_array + + @staticmethod + def get_minz0_cut(): + """ + Defines minimum z0 cut polynomial coefficients. + These coefficients were determined by optimizing this cut using ZBi as a function of MC signal mass. + The optimized cut values as a function of mass were then fit with a polynomial. + The cut shape was then tightened by 0.1 mm. + + Returns: + coeffs (list): minimum z0 cut function coefficients. + + + """ + coeffs = [1.07620094e+00 + 0.1, -7.44533811e-03, 1.58745903e-05] + return coeffs + + def minz0_sel(self,array): + """ + Applies the minimum z0 cut. + + Args: + array (awkward array): data array. + + Returns: + sel (awkward array): Boolean array representing events passing the minimum z0 cut. + """ + # Retrieve coefficients for the minimum z0 cut + coeffs = self.get_minz0_cut() + p0 = coeffs[0] + p1 = coeffs[1] + p2 = coeffs[2] + + # Get boolean mask of events that pass + sel = ( + ( array.unc_vtx_min_z0 > (p0 + p1*array.unc_vtx_mass*1000 + (p2*np.square(array.unc_vtx_mass*1000.))) ) + ) + return sel + + def mass_sel(self,array, signal_mass): + """ + Applies a mass window cut around the signal mass. + The width of the window is determined by the mass resolution of the number of sigma. + + Args: + array (awkward array): Data array. + signal_mass (float): The dark vector meson search window mass center. + + Returns: + sel (awkward array): Boolean array representing events inside the mass window. + """ + # Calcualte the lower and upper bounds of the mass search window based on the mass resolution + mass_low = signal_mass - self.nsigma*self.mass_resolution(signal_mass) + mass_high = signal_mass + self.nsigma*self.mass_resolution(signal_mass) + + # Get the boolean mask of events that pass + sel = ( + ( array.unc_vtx_mass*1000. >= {mass_low}) & (array.unc_vtx_mass*1000. <= {mass_high}) + ) + return sel + + @staticmethod + def psum_sel(array, case='sr'): + """ + Applies the Psum selection (signal region or control region). + + Args: + array (awkward array): Data array. + case (str, optional): Specify signal region ('sr') or control region ('cr'). + + Returns: + sel (awkward array): Boolean array representing events in Psum region. + """ + + if case == 'sr': + # Select momentum sum for signal region (between 1.0 and 1.9 GeV) + sel = ( + (array.unc_vtx_psum > 1.0) & (array.unc_vtx_psum < 1.9) + ) + elif case == 'cr': + # Select momentum sum for control region (between 1.9 and 2.4 GeV) + sel = ( + (array.unc_vtx_psum > 1.9) & (array.unc_vtx_psum < 2.4) + ) + else: + # No selection if invalid case + sel = () + return sel + + @staticmethod + def vprojsig_sel(array): + """ + Applies cut on the target projected vertex significance. + + Args: + array (awkward array): Data array. + + Returns: + sel (awkward array): Boolean array representing events that pass this cut. + """ + # Select events with v0projsig < 2.0 + sel = ( + (array.unc_vtx_proj_sig < 2) + ) + return sel + + @staticmethod + def sameside_z0_cut(array): + """ + Applies a cut to remove events where both tracks have z0 (aka y0) with the same sign. + This cut doesn't seem to do much after the final cut, but could be useful earlier on in the analysis. + + Args: + array (awkward array): Data array. + + Returns: + sel (awkward array): Boolean array representing events that pass this cut. + """ + sel = ( + (-1.*(array.unc_vtx_ele_track_z0*array.unc_vtx_pos_track_z0) > 0) + ) + return sel + + @staticmethod + def zcut_sel(array): + """ + Applies a cut on the reconstructed z vertex position. + + Args: + array (awkward array): Data array. + + Returns: + sel (awkward array): Boolean array representing events passing this cut. + """ + + # Select events where the reconstructed vertex z position is greater than -4.3 mm (target location) + sel = ( + (array.unc_vtx_z > -4.8) + ) + return sel + + + #Combine all of the selections into one function if you like. + def tight_selection(self, array, signal_mass, case=1): + coeffs = self.get_minz0_cut() + p0 = coeffs[0] + p1 = coeffs[1] + p2 = coeffs[2] + mass_low = signal_mass - self.nsigma*self.mass_resolution(signal_mass) + mass_high = signal_mass + self.nsigma*self.mass_resolution(signal_mass) + + if case == 1: #full tight analysis selection + sel = ( + ( array.unc_vtx_min_z0 > (p0 + p1*array.unc_vtx_mass*1000 + (p2*np.square(array.unc_vtx_mass*1000.))) ) & + ( array.unc_vtx_mass*1000. >= {mass_low}) & (array.unc_vtx_mass*1000. <= {mass_high}) & + (array.unc_vtx_proj_sig < 2) & (array.unc_vtx_z > -4.3) & (array.unc_vtx_psum > 1.0) & (array.unc_vtx_psum < 1.9) + ) + if case == 2: #tight selection without minz0 cut + sel = ( + ( array.unc_vtx_mass*1000. >= {mass_low}) & (array.unc_vtx_mass*1000. <= {mass_high}) & + (array.unc_vtx_proj_sig < 2) & (array.unc_vtx_z > -4.3) & (array.unc_vtx_psum > 1.0) & (array.unc_vtx_psum < 1.9) + ) + + if case == 3: #tight selection without mass and minz0 cut + sel = ( + (array.unc_vtx_proj_sig < 2) & (array.unc_vtx_z > -4.3) & (array.unc_vtx_psum > 1.0) & (array.unc_vtx_psum < 1.9) + ) + + return sel + + @staticmethod + def readROOTHisto(infilename, histoname): + """ + Quickly read a ROOT histogram from a root file, and make a deep copy. + + Args: + infilename (str): Input ROOT file. + histoname (str): Name of histogram being loaded. + + Returns: + histo (ROOT obj): Graph or histogram. + """ + #Open ROOT file + infile = r.TFile(f'{infilename}',"READ") + #Make copy of histogram + histo = copy.deepcopy(infile.Get(f'{histoname}')) + infile.Close() + return histo + + + @staticmethod + def cnvHistosToROOT(histos=[], tempname='temporary_uproot'): + """ + Converts hist histograms to ROOT histograms. Convenient because ROOT automatically calculates errors. + Also nice if you want to present result using ROOT histogram formatting. + + Args: + histos (list[Hist]): List of hist histograms. + tempname (str): Name of temporary ROOT file. Necessary to save Hist using uproot, and then read ROOT out. + + Returns: + return_histos (list[ROOT histograms]): ROOT histograms. + """ + return_histos = [] + uproot_file = uproot.recreate(f'trash_{tempname}.root') + for i, histo in enumerate(histos): + uproot_file[f'histo_{i}'] = histo + uproot_file.close() + infile = r.TFile(f'trash_{tempname}.root',"READ") + for i, histo in enumerate(histos): + return_histos.append(copy.deepcopy(infile.Get(f'histo_{i}'))) + infile.Close() + return return_histos + + @staticmethod + def cnvHistoToROOT(histo, tempname='temporary_uproot'): + uproot_file = uproot.recreate(f'trash_{tempname}.root') + uproot_file['histogram'] = histo + uproot_file.close() + infile = r.TFile(f'trash_{tempname}.root',"READ") + root_hist = (copy.deepcopy(infile.Get(f'histogram'))) + infile.Close() + return root_hist + + @staticmethod + def cnv_root_to_np(histo): + """ + Extract values from ROOT histogram to enable analysis or plotting with other tools. + + Args: + hist (TH1): Input ROOT histogram. + + """ + nbins = histo.GetNbinsX() + xvals = np.array([histo.GetBinCenter(x+1) for x in range(nbins+1)]) + yvals = np.array([histo.GetBinContent(x+1) for x in range(nbins+1)]) + errors = np.array([histo.GetBinError(x+1) for x in range(nbins+1)]) + underflow = histo.GetBinContent(0) + overflow = histo.GetBinContent(nbins+1) + + #add over/underflow + xvals = np.insert(xvals, 0, xvals[0]-histo.GetBinWidth(1)) + yvals = np.insert(yvals, 0, underflow) + xvals = np.append(xvals, xvals[-1]+histo.GetBinWidth(1)) + yvals = np.append(yvals, overflow) + errors = np.insert(errors, 0, 0.0) + errors = np.append(errors, 0.0) + + #get fit function if it exist + x_fit = None + y_fit = None + if len(histo.GetListOfFunctions()) > 0: + fitfunc = histo.GetListOfFunctions()[0] + x_fit = np.linspace(fitfunc.GetXmin(), fitfunc.GetXmax(), int((fitfunc.GetXmax()-fitfunc.GetXmin())/histo.GetBinWidth(1))) + y_fit = np.array([fitfunc.Eval(x) for x in x_fit]) + + return (xvals, yvals, errors), (x_fit, y_fit) + + @staticmethod + def cnv_tgraph_to_np(tgraph): + """ + Extract values from ROOT Tgraph. + + Args: + tgraph (TGraph): Input ROOT Tgraph + """ + # Number of points in the TGraph + npoints = tgraph.GetN() + + # Retrieve X and Y values + xvals = np.array([tgraph.GetX()[i] for i in range(npoints)]) + yvals = np.array([tgraph.GetY()[i] for i in range(npoints)]) + + #Errors not available in standard TGraph, set to zero. + errors = np.zeros(npoints) + + # Handle fit function if it exists + x_fit = None + y_fit = None + if len(tgraph.GetListOfFunctions()) > 0: + fitfunc = tgraph.GetListOfFunctions()[0] + x_fit = np.linspace(fitfunc.GetXmin(), fitfunc.GetXmax(), 100) # 100 points for the fit + y_fit = np.array([fitfunc.Eval(x) for x in x_fit]) + + return (xvals, yvals, errors), (x_fit, y_fit) + + + @staticmethod + def fit_plot_with_poly(plot, tgraph=False, specify_n=None, set_xrange=False, xrange=(0.0, 1.0)): + """ + Fit TH1 or TGraph with polynomial fit function. Uses fstat test to guide choice of polynomial degree. + + Args: + plot (TH1 or TGraph): Input plot to be fitted. + tgraph (boolean, optional): If True, access npoints through TGraph method. + specify_n (int, optional): If None, fit plot using all n degrees, calculate and print fstat for each n. + If int, fit plot using n degree polynomial. + set_xrange (boolean, optional): If True, set fit range according to xrange. + xrange (float, float, optional): Define fit range. + + Returns: + params (list): List of fit function parameters. + errors (list): List of fit parameter errors. + + Notes: + If trying to determine what order polynomial to fit plot, specify_n=None. Function will print out fstat results. + Once best order determined, set specifcy_n= to get fit resulst. + """ + polys = [] + chi2s = [] + fstats = [] + fit_resultults = [] + + if tgraph: + npoints = plot.GetN() + else: + npoints = 0 + nBins = plot.GetNbinsX() + for ibin in range(nBins): + if plot.GetBinContent(ibin) > 0: + npoints += 1 + pass + + + if not specify_n: + for n in range(11): + fitfunc = r.TF1(f'pol{n}',f'pol{n}') + fitfunc.SetLineColor(r.kRed) + if set_xrange: + fitfunc.SetRange(xrange[0], xrange[1]) + fit_result = plot.Fit(fitfunc,"RSQ") + else: + fit_result = plot.Fit(fitfunc,"SQ") + fitfunc.SetLineColor(r.kRed) + fitfunc.SetMarkerSize(0.0) + chi2s.append(fit_result.Chi2()) + polys.append(n) + fit_resultults.append(fit_result) + + #Perform fstat test to see how much fit improves with additional order + if n > 0: + fstats.append( (chi2s[n-1]-chi2s[n])*(npoints-n-1)/(chi2s[n])) + else: + fstats.append(0.0) + + print(fstats) + return None, None + else: + fitfunc = r.TF1(f'pol{specify_n}',f'pol{specify_n}') + fitfunc.SetLineColor(r.kRed) + fitfunc.SetLineWidth(5) + if set_xrange: + fitfunc.SetRange(xrange[0], xrange[1]) + fit_result = plot.Fit(fitfunc,"RSQ") + else: + fit_result = plot.Fit(fitfunc,"SQ") + params = fit_result.Parameters() + errors = fit_result.Errors() + #return fit_result + return params, errors + + + + def systematic_uncertainties(self): + """ + This method contains all of the systematic uncertainties applied in Alic's 2016 SIMP L1L1 dissertation. + The detector misalignment systematic is not included as of 09/06/2024. + These functions and values were calculated externally using a +- 1.5 sigma search window, with cuts frozen per dissertation. + All polynomials are a function of A' mass. + """ + self.radacc_targetz_nominal = self.polynomial(0.24083419, -0.017612076, 0.00037553660, -1.0223921e-06, -3.8793240e-08, + 4.2199609e-10, -1.6641414e-12, 2.3433278e-15) + self.radacc_targetz_Mpt5 = self.polynomial(0.22477846, -0.015984559, 0.00030943435, 3.6182165e-07, -5.4820194e-08, + 5.2531952e-10, -2.0102027e-12, 2.8109430e-15) + self.radacc_targetz_Ppt5 = self.polynomial( 0.22779999, -0.016020742, 0.00029960205, 7.6823260e-07, -6.0956281e-08, + 5.6914810e-10, -2.1602258e-12, 3.0094146e-15) + self.simp_targetz = self.polynomial(-1.38077250e+00, 8.00749424e-02, -9.78327706e-04, 5.13294008e-06, -9.77393492e-09) + self.mass_unc = 0.043 + self.radfrac = 0.07 + + def evaluate_polynomials(self, mass): + """ + Returns the systematic uncertainties defined in systematic_uncertainties() as a function of mass. + + Args: + mass (float): A' mass. + + Note: + You need to understand each systematic in order to correctly combine these numbers. + """ + nominal_values = self.radacc_targetz_nominal(mass) + Mpt5_values = self.radacc_targetz_Mpt5(mass) + Ppt5_values = self.radacc_targetz_Ppt5(mass) + simp_values = self.simp_targetz(mass) + + return nominal_values, Mpt5_values, Ppt5_values, simp_values, self.mass_unc, self.radfrac + + @staticmethod + def inject_signal_mc(signal, data, nevents=100): + """ + Randomly selects MC signal events and injects them into data array. + + Args: + signal (awkward array): MC signal array loaded using load_signal(). + data (awwkard array): Data array loaded using load_data(). + nevents (int): Specify number of MC signal events to inject into data. + + Returns: + injected_data (awkward array): Copy of data array with injected signal. + thrown_events (int): Number of MC signal events that were thrown. Useful sanity check. + """ + # Identify the maximum signal event weight in the MC signal array + max_weight = np.max(signal.expected_signal_weight) + + # Randomly sample the MC signal array until the specified number of events is thrown + events_thrown = 0 + thrown_mask = [] + while events_thrown < nevents: + + # Randomly select a signal event from the array + rint = np.random.randint(0,len(signal.expected_signal_weight)-1) + random_event = signal[rint] + + # Randomly sample the uniform distribution between 0-maximum signal event weight + # If the uniform sample weight is less than the randomly selected event weight, thrown the event + rweight = np.random.uniform(0, max_weight) + if rweight < random_event.expected_signal_weight: + events_thrown += 1 + thrown_mask.append(rint) + thrown_events = signal[thrown_mask] + thrown_events['weight'] = 1.0 + + # Inject the randomly selected signal events into the data by combining the awkward arrays + injected_data = ak.concatenate([data, thrown_events]) + return injected_data, thrown_events + +#======================================================================================================================================= +# MAIN: Calculate the expected signal +#======================================================================================================================================= + +if __name__ == '__main__': + + #parse input arguments + parser = argparse.ArgumentParser(description='Process some inputs.') + parser.add_argument('--outfilename', type=str, default='expected_signal_output.root') + parser.add_argument('--mpifpi', type=float, default=4*np.pi) + parser.add_argument('--signal_sf', type=float, default=1.0) + parser.add_argument('--nsigma', type=float, default=1.5) + parser.add_argument('--tenpct', type=int, default=0) + args = parser.parse_args() + + mpifpi = args.mpifpi + nsigma = args.nsigma + signal_sf = args.signal_sf + outfilename = args.outfilename + tenpct = args.tenpct + + # Initialize signal processor + processor = SignalProcessor(mpifpi=mpifpi, nsigma=nsigma) + + # Load either 10% data or 100% data, use the CR reconstructed bkg rate to scale the expected signal rate + cr_data = '/sdf/group/hps/user-data/alspellm/2016/data/hadd_BLPass4c_1959files.root' + full_lumi_path = '/fs/ddn/sdf/group/hps/users/alspellm/data_storage/pass4kf/pass4kf_ana_20240513' + preselection = "vtxana_Tight_nocuts" + signal_mass_range = [x for x in range(30,130,1)] + processor.set_diff_prod_lut(cr_data, preselection, signal_mass_range, tenpct, full_lumi_path) + + # Initialize the mass and epsilon^2 range for the expected signal calculation + mass_max = 124 + mass_min = 30 + mass_step = 2 # MC signal files were generated at 2 MeV increments + ap_step = round(mass_step*processor.mass_ratio_ap_to_vd,1) + masses = np.array([x for x in range(mass_min, mass_max+mass_step, mass_step)]) + ap_masses = np.array([round(x*processor.mass_ratio_ap_to_vd,1) for x in masses]) + eps2_range = np.logspace(-4.0,-8.0,num=40) + logeps2_range = np.log10(eps2_range) + min_eps = min(np.log10(eps2_range)) + max_eps = max(np.log10(eps2_range)) + num_bins = len(eps2_range) + + # Initialize the histograms used to store the expected signal. One for Vd mass, one for Ap mass. + expected_signal_vd_h = ( + hist.Hist.new + .Reg(len(masses), np.min(masses), np.max(masses)+mass_step, label='Vd Invariant Mass [MeV]') + .Reg(num_bins, min_eps, max_eps,label=f'$log_{10}(\epsilon^2)$') + .Double() + ) + expected_signal_ap_h = ( + hist.Hist.new + .Reg(len(ap_masses), np.min(ap_masses), np.max(ap_masses)+ap_step, label='A\' Invariant Mass [MeV]') + .Reg(num_bins, min_eps, max_eps,label=f'$log_{10}(\epsilon^2)$') + .Double() + ) + + # Calculate expected signal for each MC generated mass + for signal_mass in masses: + + #Load MC Signal + indir = '/sdf/group/hps/user-data/alspellm/2016/simp_mc/pass4b/beam/smeared_fixbeamspot' + signal_pre_readout_path = lambda mass: f'/sdf/group/hps/user-data/alspellm/2016/simp_mc/pass4b/nobeam/mass_{mass}_simp_2pt3_slic_hadd_ana.root' + signal_path = lambda mass: f'{indir}/mass_{mass}_hadd-simp-beam_ana_smeared_corr.root' + signal_selection = 'vtxana_radMatchTight_2016_simp_SR_analysis' + + #Get the total signal yield as a function of eps2 + total_yield_per_epsilon2 = processor.total_signal_production_per_epsilon2(signal_mass) + print('Total Yield Per eps2: ', total_yield_per_epsilon2) + + print('Load Signal ', signal_path(signal_mass)) + signal = processor.load_signal(signal_path(signal_mass), signal_pre_readout_path(signal_mass), signal_mass, signal_selection) + + #Get Tight selection + psum_sel = processor.psum_sel(signal, case='sr') + zcut_sel = processor.zcut_sel(signal) + vprojsig_sel = processor.vprojsig_sel(signal) + minz0_sel = processor.minz0_sel(signal) + masswindow_sel = processor.mass_sel(signal, signal_mass) + #sameside_sel = processor.sameside_z0_cut(signal) + + # Combine selections + #tight_sel = np.logical_and.reduce([psum_sel,zcut_sel, vprojsig_sel, psum_sel, minz0_sel, masswindow_sel, sameside_sel]) + tight_sel = np.logical_and.reduce([psum_sel,zcut_sel, vprojsig_sel, psum_sel, minz0_sel, masswindow_sel]) + + for l, eps2 in enumerate(eps2_range): + signal = processor.get_exp_sig_eps2(signal_mass, signal, eps2) + total_yield = signal_sf*ak.sum(signal['reweighted_accxEff'][tight_sel])*total_yield_per_epsilon2*eps2 + expected_signal_vd_h.fill(signal_mass, logeps2_range[l], weight=total_yield) + expected_signal_ap_h.fill(signal_mass*processor.mass_ratio_ap_to_vd, logeps2_range[l], weight=total_yield) + + outfile = uproot.recreate(outfilename) + outfile['expected_signal_vd_h'] = expected_signal_vd_h + outfile['expected_signal_ap_h'] = expected_signal_ap_h + + + + + + + + + + + + + + + diff --git a/plotUtils/simps/simp_theory_equations.py b/plotUtils/simps/simp_theory_equations.py new file mode 100644 index 000000000..d6cf16a56 --- /dev/null +++ b/plotUtils/simps/simp_theory_equations.py @@ -0,0 +1,172 @@ +import math +import ROOT as r +import numpy as np +class SimpEquations: + + def __init__(self, year = 2016, alpha_dark = 0.01, mass_ratio_Ap_to_Vd = 1.66, mass_ratio_Ap_to_Pid = 3.0, + ratio_mPi_to_fPi = 12.566, lepton_mass = 0.511): + self.year = year + self.alpha_dark = alpha_dark + self.mass_ratio_Ap_to_Vd = mass_ratio_Ap_to_Vd + self.mass_ratio_Ap_to_Pid = mass_ratio_Ap_to_Pid + self.ratio_mPi_to_fPi = ratio_mPi_to_fPi + self.lepton_mass = lepton_mass + + @staticmethod + def rate_Ap_ee(m_Ap, eps): + ml = 0.511 + r = ml/m_Ap + coeff1 = ((1.0/137.0)*eps**2)/3.0 + coeff2 = (1.0 - 4.0*(r**2))**(0.5) + coeff3 = (1.0 + 2.0*(r**2))*m_Ap + return coeff1*coeff2*coeff3 + + @staticmethod + def rate_2pi(m_Ap, m_pi, m_V, alpha_dark): + coeff = (2.0 * alpha_dark / 3.0) * m_Ap + pow1 = np.power((1 - (4 * m_pi * m_pi / (m_Ap * m_Ap))), 3 / 2.0) + pow2 = np.power(((m_V * m_V) / ((m_Ap * m_Ap) - (m_V * m_V))), 2) + return coeff * pow1 * pow2 + + @staticmethod + def rate_Vrho_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi): + x = m_pi / m_Ap + y = m_V / m_Ap + Tv = 3.0/4.0 + coeff = alpha_dark * Tv / (192.0 * np.power(math.pi, 4)) + return coeff * np.power((m_Ap / m_pi), 2) * np.power(m_V / m_pi, 2) * np.power((m_pi / f_pi), 4) * m_Ap * np.power(SimpEquations.Beta(x, y), 3 / 2.0) + + def rate_Vphi_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi): + x = m_pi / m_Ap + y = m_V / m_Ap + Tv = 3.0/2.0 + coeff = alpha_dark * Tv / (192.0 * np.power(math.pi, 4)) + return coeff * np.power((m_Ap / m_pi), 2) * np.power(m_V / m_pi, 2) * np.power((m_pi / f_pi), 4) * m_Ap * np.power(SimpEquations.Beta(x, y), 3 / 2.0) + + def rate_Vcharged_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi): + x = m_pi / m_Ap + y = m_V / m_Ap + Tv = 18.0 - ((3.0/2.0)+(3.0/4.0)) + coeff = alpha_dark * Tv / (192.0 * np.power(math.pi, 4)) + return coeff * np.power((m_Ap / m_pi), 2) * np.power(m_V / m_pi, 2) * np.power((m_pi / f_pi), 4) * m_Ap * np.power(SimpEquations.Beta(x, y), 3 / 2.0) + + @staticmethod + def br_2pi(m_Ap, m_pi, m_V, alpha_dark, f_pi): + total_rate = (SimpEquations.rate_Vrho_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) + + SimpEquations.rate_Vphi_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) + + SimpEquations.rate_Vcharged_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) + + SimpEquations.rate_2pi(m_Ap, m_pi, m_V, alpha_dark)) + + # Add the rate_2V contribution element-wise where m_Ap > 2.0 * m_V + total_rate += np.where(m_Ap > 2.0 * m_V, SimpEquations.rate_2V(m_Ap, m_V, alpha_dark), 0.0) + + return SimpEquations.rate_2pi(m_Ap, m_pi, m_V, alpha_dark) / total_rate + + @staticmethod + def br_Vrho_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi): + total_rate = (SimpEquations.rate_Vrho_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) + + SimpEquations.rate_Vphi_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) + + SimpEquations.rate_Vcharged_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) + + SimpEquations.rate_2pi(m_Ap, m_pi, m_V, alpha_dark)) + + # Add the rate_2V contribution element-wise where m_Ap > 2.0 * m_V + total_rate += np.where(m_Ap > 2.0 * m_V, SimpEquations.rate_2V(m_Ap, m_V, alpha_dark), 0.0) + + return SimpEquations.rate_Vrho_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) / total_rate + + @staticmethod + def br_Vphi_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi): + total_rate = (SimpEquations.rate_Vrho_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) + + SimpEquations.rate_Vphi_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) + + SimpEquations.rate_Vcharged_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) + + SimpEquations.rate_2pi(m_Ap, m_pi, m_V, alpha_dark)) + + # Add the rate_2V contribution element-wise where m_Ap > 2.0 * m_V + total_rate += np.where(m_Ap > 2.0 * m_V, SimpEquations.rate_2V(m_Ap, m_V, alpha_dark), 0.0) + + return SimpEquations.rate_Vphi_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) / total_rate + + @staticmethod + def br_Vcharged_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi): + total_rate = (SimpEquations.rate_Vrho_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) + + SimpEquations.rate_Vphi_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) + + SimpEquations.rate_Vcharged_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) + + SimpEquations.rate_2pi(m_Ap, m_pi, m_V, alpha_dark)) + + # Add the rate_2V contribution element-wise where m_Ap > 2.0 * m_V + total_rate += np.where(m_Ap > 2.0 * m_V, SimpEquations.rate_2V(m_Ap, m_V, alpha_dark), 0.0) + + return SimpEquations.rate_Vcharged_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) / total_rate + + @staticmethod + def br_2V(m_Ap, m_pi, m_V, alpha_dark, f_pi): + total_rate = (SimpEquations.rate_Vrho_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) + + SimpEquations.rate_Vphi_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) + + SimpEquations.rate_Vcharged_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) + + SimpEquations.rate_2pi(m_Ap, m_pi, m_V, alpha_dark)) + + # Calculate the rate_2V contribution only where m_Ap > 2.0 * m_V and add it to total_rate + total_rate += np.where(m_Ap > 2.0 * m_V, SimpEquations.rate_2V(m_Ap, m_V, alpha_dark), 0.0) + + # Return 0.0 where 2 * m_V >= m_Ap + rate_2V_contrib = np.where((2.0 * m_V >= m_Ap) | (m_Ap <= 2.9 * m_V), 0.0, SimpEquations.rate_2V(m_Ap, m_V, alpha_dark)) + + return rate_2V_contrib / total_rate + + @staticmethod + def Tv(rho, phi): + if rho: + return 3.0 / 4.0 + elif phi: + return 3.0 / 2.0 + else: + return 18.0 + + @staticmethod + def Beta(x, y): + return (1 + np.power(y, 2) - np.power(x, 2) - 2 * y) * (1 + np.power(y, 2) - np.power(x, 2) + 2 * y) + + @staticmethod + def rate_2V(m_Ap, m_V, alpha_dark): + r = m_V / m_Ap + return alpha_dark / 6.0 * m_Ap * SimpEquations.f(r) + + @staticmethod + def f(r): + # Define your function f(r) here + # Example: return some_expression + return -1. + + @staticmethod + def rate_2l(m_Ap, m_pi, m_V, eps, alpha_dark, f_pi, m_l, rho): + alpha = 1.0 / 137.0 + coeff = (16 * math.pi * alpha_dark * alpha * eps**2 * f_pi**2) / (3 * m_V**2) + term1 = (m_V**2 / (m_Ap**2 - m_V**2))**2 + term2 = (1 - (4 * m_l**2 / m_V**2))**0.5 + term3 = 1 + (2 * m_l**2 / m_V**2) + constant = 1 if not rho else 2 + return coeff * term1 * term2 * term3 * m_V * constant + + @staticmethod + def getCtau(m_Ap, m_pi, m_V, eps, alpha_dark, f_pi, m_l, rho): + #c = 3.00e10 # cm/s 0: + fitfunc = histo.GetListOfFunctions()[0] + x_fit = np.linspace(fitfunc.GetXmin(), fitfunc.GetXmax(), int((fitfunc.GetXmax()-fitfunc.GetXmin())/histo.GetBinWidth(1))) + y_fit = np.array([fitfunc.Eval(x) for x in x_fit]) + + return (xvals, yvals, errors), (x_fit, y_fit) +#======================================================================================================================================= +# LOAD DATA +#======================================================================================================================================= + +# Initialize signal processor +search_window = 1.5 #used in final search +signalProcessor = simp_signal_2016.SignalProcessor(args.mpifpi, search_window) + +# Load data, MC background, and signal +samples = {} +branches = ["unc_vtx_ele_track_z0","unc_vtx_pos_track_z0"] + +# Load 10% Data +infile = '/sdf/group/hps/user-data/alspellm/2016/data/hadd_BLPass4c_1959files.root' +selection = 'vtxana_Tight_L1L1_nvtx1' +samples['data'] = signalProcessor.load_data(infile,selection, expressions=branches, cut_expression='((unc_vtx_psum > 1.0) & (unc_vtx_psum < 1.9) )') +samples['data']['weight'] = 1.0 #Assign weight of 10 to scale up to full lumi + +# Load MC background +lumi = 10.7*.1 #pb-1 +mc_scale = {'data' : 1.0, + 'tritrig' : 1.416e9*lumi/(50000*10000), + 'wab' : 0.1985e12*lumi/(100000*10000)} + +# Load MC tritrig +infile = '/sdf/group/hps/user-data/alspellm/2016/tritrig_mc/pass4b/hadded_tritrig-beam-10kfiles-ana-smeared-corr_beamspotfix.root' +samples['tritrig'] = signalProcessor.load_data(infile, selection, cut_expression='((unc_vtx_psum > 1.0) & (unc_vtx_psum < 1.9) )', expressions=branches) +samples['tritrig']['weight'] = mc_scale['tritrig'] + +# Load MC wab +infile = '/sdf/group/hps/user-data/alspellm/2016/wab_mc/pass4b/hadded_wab-beam-10kfiles-ana-smeared-corr_beamspotfix.root' +samples['wab'] = signalProcessor.load_data(infile, selection, cut_expression='((unc_vtx_psum > 1.0) & (unc_vtx_psum < 1.9) )', expressions=branches) +samples['wab']['weight'] = mc_scale['wab'] + +#======================================================================================================================================= +# Check smearing after the fact +#======================================================================================================================================= + +# After the smearing factor has been calculated, set to True to compare data and smeared MC background +smear = True +if smear: + mc_sigma = 0.1251 # Calculated in this script + data_sigma = 0.1348 # Calculated in this script + smearF = np.sqrt(data_sigma**2 - mc_sigma**2) + + # Smear MC tritrig + rel_smear = np.random.normal(0.0, 1.0, len(samples['tritrig'].unc_vtx_min_z0)) + smearfactors = rel_smear*smearF + samples['tritrig']['unc_vtx_ele_track_z0'] = smearfactors + samples['tritrig']['unc_vtx_ele_track_z0'] + rel_smear = np.random.normal(0.0, 1.0, len(samples['tritrig'].unc_vtx_min_z0)) + smearfactors = rel_smear*smearF + samples['tritrig']['unc_vtx_pos_track_z0'] = smearfactors + samples['tritrig']['unc_vtx_pos_track_z0'] + + # Smear MC wab + rel_smear = np.random.normal(0.0, 1.0, len(samples['wab'].unc_vtx_min_z0)) + smearfactors = rel_smear*smearF + samples['wab']['unc_vtx_ele_track_z0'] = smearfactors + samples['wab']['unc_vtx_ele_track_z0'] + rel_smear = np.random.normal(0.0, 1.0, len(samples['wab'].unc_vtx_min_z0)) + smearfactors = rel_smear*smearF + samples['wab']['unc_vtx_pos_track_z0'] = smearfactors + samples['wab']['unc_vtx_pos_track_z0'] + +#======================================================================================================================================= +# Calculate z0 width in data and MC background +#======================================================================================================================================= + +# Initialize z0 histograms +z0_h = ( + hist.Hist.new + .StrCategory(list(samples.keys()), name='samples') + .Reg(200, -1.,1.,label=r'Vertical Track Impact Parameter [mm]') + .Double() +) + +# Fill without weights, so that histos can be converted to ROOT and retain statistical uncertainty +z0_histos = {} +for sname, sample in samples.items(): + z0_h.fill(sname, sample.unc_vtx_ele_track_z0) + z0_h.fill(sname, sample.unc_vtx_pos_track_z0) + z0_histos[sname] = signalProcessor.cnvHistoToROOT(z0_h[sname,:]) # Convert hist histogram to ROOT + z0_histos[sname].Scale(mc_scale[sname]) + +# Scale Tritrig and WAB and combine with proper errors +z0_histos['tritrig_wab'] = z0_histos['tritrig'].Clone() +z0_histos['tritrig_wab'].Add(z0_histos['wab']) + +# Normalize histograms +for sname, sample in z0_histos.items(): + print(z0_histos[sname].Integral(0,-1)) + z0_histos[sname].Scale(1./z0_histos[sname].Integral(0,-1)) + +# Plot data and MC background z0 distributions +# Fit each with Gaussian to determine z0 width +# MC width narrower than data. Smear MC background to match data +fig, ax = plt.subplots(2,1, figsize=(25,30)) +#Data +plt.subplot(2,1,1) +plt.xlabel('Vertical Track Impact Parameter [mm]') +plt.ylabel('Normalized Events') +_, params, fiterrors, _ = fit_root_gaussian(z0_histos['data'], nsigma=1.5) +(xvals, yvals, errors), (x_fit, y_fit) = cnv_root_to_np(z0_histos['data']) +plt.errorbar(xvals, yvals, yerr=errors, marker='o', color='black', linestyle='', label='10% Data') +plt.plot(x_fit, y_fit, color='red', label='Fit 10% Data', linewidth=2.0) +plt.text(0.25, 0.015, f'Norm={params[0]}$\pm${fiterrors[0]} \n$\mu$={params[1]}$\pm${fiterrors[1]} \n$\sigma$={params[2]}$\pm${fiterrors[2]}') +plt.legend() +#MC +plt.subplot(2,1,2) +plt.xlabel('Vertical Track Impact Parameter [mm]') +plt.ylabel('Normalized Events') +_, params, fiterrors, _ = fit_root_gaussian(z0_histos['tritrig_wab'], nsigma=1.5) +(xvals, yvals, errors), (x_fit, y_fit) = cnv_root_to_np(z0_histos['tritrig_wab']) +plt.errorbar(xvals, yvals, yerr=errors, marker='o', color='darkblue', linestyle='', label='MC Background') +plt.plot(x_fit, y_fit, color='red', label='Fit MC Background', linewidth=2.0) +plt.text(0.25, 0.015, f'Norm={params[0]}$\pm${fiterrors[0]} \n$\mu$={params[1]}$\pm${fiterrors[1]} \n$\sigma$={params[2]}$\pm${fiterrors[2]}') +plt.legend() +plt.ylim(0.0, 0.03) +plt.savefig(f'{outdir}/impact_parameter_data_v_mc_smeared_{smear}.png') + +#smearing factors calculated from comparing data and MC bkg z0 widths +mc_sigma = 0.1251 +data_sigma = 0.1348 +smearF = np.sqrt(data_sigma**2 - mc_sigma**2) + +#======================================================================================================================================= +# Smear MC signal and calculate change in signal efficiency +#======================================================================================================================================= + +sysvals = [] +masses = [] +# Directory containing MC signal hpstr vertex ana processor tuples +indir = '/sdf/group/hps/user-data/alspellm/2016/simp_mc/pass4b/beam/smeared_fixbeamspot' +for mass in range(30,120,4): + masses.append(mass) + # MC signal hpstr MC ana processor (truth vertex z information) + signal_pre_readout_path = lambda mass: f'/sdf/group/hps/user-data/alspellm/2016/simp_mc/pass4b/nobeam/mass_{mass}_simp_2pt3_slic_hadd_ana.root' + signal_path = lambda mass: f'{indir}/mass_{mass}_hadd-simp-beam_ana_smeared_corr.root' + signal_selection = 'vtxana_radMatchTight_2016_simp_SR_analysis' + signal = signalProcessor.load_signal(signal_path(signal_mass), signal_pre_readout_path(signal_mass), signal_mass, signal_selection) + signal['weight']=1.0 + psum_sel = signalProcessor.psum_sel(signal, case='sr') #Psum signal region selection + + # Smear electron and positron track z0 + rel_smear = np.random.normal(0.0, 1.0, len(signal.unc_vtx_min_z0)) + smearfactors = rel_smear*smearF + signal['unc_vtx_ele_track_z0_smeared'] = smearfactors + signal['unc_vtx_ele_track_z0'] + + rel_smear = np.random.normal(0.0, 1.0, len(signal.unc_vtx_min_z0)) + smearfactors = rel_smear*smearF + signal['unc_vtx_pos_track_z0_smeared'] = smearfactors + signal['unc_vtx_pos_track_z0'] + + # Calculate smeared minz0 + signal['unc_vtx_min_z0_smeared'] = np.minimum(abs(signal['unc_vtx_ele_track_z0_smeared']), abs(signal['unc_vtx_pos_track_z0_smeared'])) + + #Calculate change in efficiency + unsmeared_sel = signal.unc_vtx_min_z0 > signalProcessor.minz0_cut_poly(signal.unc_vtx_mass*1000.) + smeared_sel = signal.unc_vtx_min_z0_smeared > signalProcessor.minz0_cut_poly(signal.unc_vtx_mass*1000.) + unsmeared_eff = ak.sum(signal[unsmeared_sel].weight)/ak.sum(signal.weight) + smeared_eff = ak.sum(signal[smeared_sel].weight)/ak.sum(signal.weight) + systematic = smeared_eff/unsmeared_eff + sysvals.append(systematic) + + #final systematic stored in sysvals + print(sysvals) + diff --git a/plotUtils/simps/systematics/radiative_acceptance_systematic-misalignment.py b/plotUtils/simps/systematics/radiative_acceptance_systematic-misalignment.py new file mode 100644 index 000000000..11e1ac94f --- /dev/null +++ b/plotUtils/simps/systematics/radiative_acceptance_systematic-misalignment.py @@ -0,0 +1,283 @@ +#!/usr/bin/python3 +#================================================================================================================================== +# Description: Compares the radiative acceptance between nominal and misaligned HPS detector versions +# Author: Alic Spellman +# Date: 09/05/2024 +# Script to load MC samples, plot and compute misalignment systematics for 2016 simp L1L1 analysis +# Calculates both radiative trident acceptance and signal acceptance +import os +import awkward as ak +import numpy as np +import hist +from hist import Hist +import uproot +import ROOT as r +import copy +import matplotlib.pyplot as plt +import matplotlib as mpl +import matplotlib.gridspec as gridspec +import sys +import math + +#simp tools defined in hpstr +hpstr_base = os.getenv('HPSTR_BASE') +sys.path.append(f'{hpstr_base}/plotUtils/simps') +import simp_signal_2016 + +#====================================================================================================================================== +#INITIALIZATION +#======================================================================================================================================= +# Set plotting parameters for matplotlib +plt.rcParams.update({'font.size': 40, 'axes.titlesize': 40, 'axes.labelsize': 40, 'xtick.labelsize': 40, 'ytick.labelsize': 40, 'lines.linewidth': 3.0, 'legend.fontsize': 40}) +plt.rcParams['font.family'] = 'DejaVu Sans' + +#parse input arguments +import argparse +parser = argparse.ArgumentParser(description='') +parser.add_argument('--outdir', type=str, default='./search_results') +parser.add_argument('--mpifpi', type=float, default=4.*np.pi) + +args = parser.parse_args() +outdir = args.outdir + +#======================================================================================================================================= +# LOAD DATA: Initialize signal processor and load radiative trident MC samples +#======================================================================================================================================= +search_window = 1.5 +signalProcessor = simp_signal_2016.SignalProcessor(args.mpifpi, search_window) + +##Load MC samples for nominal and misaligned detectors +samples = {} +mcsamples = {} +branches = ["unc_vtx_mass", "unc_vtx_psum"] + +#Load reconstructed and selected radiative events for nominal detector +infile = '/sdf/group/hps/user-data/alspellm/2016/systematics/radacc/misalignments/hadd_2kfiles_rad_nobeam_nominal_recon_ana.root' +selection = 'vtxana_radMatchTight_nocuts' #USE RADMATCHTIGHT! +samples['nominal'] = signalProcessor.load_data(infile, selection, cut_expression='((unc_vtx_psum > 1.9) & (unc_vtx_psum < 2.4) )', expressions=branches) + +#Load generated events (mc ana) for nominal detector +infile = '/sdf/group/hps/user-data/alspellm/2016/systematics/radacc/misalignments/hadd_2kfiles_rad_nobeam_nominal_mc_ana.root' +slicfile = r.TFile(infile, "READ") +mcsamples['nominal'] = copy.deepcopy(slicfile.Get('mcAna/mcAna_mc622Mass_h')) +slicfile.Close() + +#Load reconstructed and selected radiative events for misaligned detector +infile = '/sdf/group/hps/user-data/alspellm/2016/systematics/radacc/misalignments/hadd_2kfiles_rad_nobeam_misalignments_1_recon_ana.root' +selection = 'vtxana_radMatchTight_nocuts' #USE RADMATCHTIGHT! +samples['misaligned_v1'] = signalProcessor.load_data(infile, selection, cut_expression='((unc_vtx_psum > 1.9) & (unc_vtx_psum < 2.4) )', expressions=branches) + +#Load generated events (mc ana) for misaligned detector +infile = '/sdf/group/hps/user-data/alspellm/2016/systematics/radacc/misalignments/hadd_2kfiles_rad_nobeam_misalignments_1_mc_ana.root' +slicfile = r.TFile(infile, "READ") +mcsamples['misaligned_v1'] = copy.deepcopy(slicfile.Get('mcAna/mcAna_mc622Mass_h')) +slicfile.Close() + +#output file to store systematic results +outfile = uproot.recreate(f'{outdir}/radacc_misalignment_systematic_results.root') + +#======================================================================================================================================= +# CHECK RADIATIVE PEAK: Ensure that the radiative peaks for misaligned and nominal are commensurate. +# Peak will be off if overly-misaligned +#======================================================================================================================================= +psum_h = ( + hist.Hist.new + .StrCategory(list(samples.keys()), name='samples') + .Reg(50,1.9, 2.4,label='Psum [GeV]') + .Double() +) +colors=['black', 'darkred'] +fig, ax = plt.subplots(figsize=(25,15)) +for i,(sname, sample) in enumerate(samples.items()): + psum_h.fill(sname, sample.unc_vtx_psum) +psum_h.plot(color=['black','darkred'], linewidth=3.0) +plt.legend() +plt.ylabel('Events') +plt.savefig(f'{outdir}/radiative_peak_misaligned.png') + +#======================================================================================================================================= +# RADIATIVE ACCEPTANCE HISTOGRAMS: Initialize invariant mass histograms, rebin, and convert them to ROOT histograms +#======================================================================================================================================= +#invariant mass histogram binning must match mc ana invariant mass to take ratio +nbinsx = mcsamples['nominal'].GetNbinsX() +first_bin = mcsamples['nominal'].GetBinLowEdge(1) +last_bin = nbinsx*mcsamples['nominal'].GetBinWidth(1) +invmass_h = ( + hist.Hist.new + .StrCategory(list(samples.keys()), name='samples') + .Reg(nbinsx,first_bin,last_bin,label='Invariant Mass [MeV]') + .Double() +) + +#Fill mc components without weights and convert to ROOT, rebin +invmass_histos = {} +for sname, sample in samples.items(): + invmass_h.fill(sname, sample.unc_vtx_mass*1000.) + invmass_histos[sname] = signalProcessor.cnvHistoToROOT(invmass_h[sname,:]) + invmass_histos[sname].Rebin(2) + mcsamples[sname].Rebin(2) + +def nonUniBinning(histo, start, size): + edges_a = np.arange(histo.GetBinLowEdge(1),start+histo.GetBinWidth(1),histo.GetBinWidth(1)) + edges_b = np.arange(start,histo.GetBinLowEdge(histo.GetNbinsX()), size) + bin_edges = np.concatenate([edges_a, edges_b[1:]]) + histo_rebinned = r.TH1F(f'{histo.GetName()}_rebinned', f'{histo.GetTitle()}', len(bin_edges)-1, bin_edges) + for bin in range(1, histo.GetNbinsX() + 1): + content = histo.GetBinContent(bin) + center = histo.GetBinCenter(bin) + error = histo.GetBinError(bin) + new_bin = histo_rebinned.FindBin(center) + histo_rebinned.SetBinContent(new_bin, histo_rebinned.GetBinContent(new_bin)+content) + histo_rebinned.SetBinError(new_bin, np.sqrt(histo_rebinned.GetBinError(new_bin)**2 + error**2)) + return histo_rebinned + +#enable non-uniform binning +for sname, sample in samples.items(): + invmass_histos[sname] = nonUniBinning(invmass_histos[sname], 150, 4) + mcsamples[sname] = nonUniBinning(mcsamples[sname], 150, 4) + outfile[f'recon_{sname}'] = invmass_histos[sname] + outfile[f'mc_{sname}'] = mcsamples[sname] + +#======================================================================================================================================= +# RADIATIVE ACCEPTANCE: Compute and plot the radiative acceptance for nominal and misaligned +#======================================================================================================================================= + +fits = {} +colors = ['#d62728', '#bcbd22', '#2ca02c', '#17becf', '#1f77b4', '#9467bd', '#7f7f7f'] +colors = ['black', 'darkred', 'darkblue', 'darkgreen', 'darkorange'] + +#Figure to plot radiative acceptance and systematic uncertainty +fig, ax = plt.subplots(2,1,figsize=(20,20)) +plt.subplot(2,1,1) +plt.xlabel('Invariant Mass [MeV]') +plt.ylabel('Radiative Acceptance') +plt.ylim(0.0, .15) +plt.xlim(20.0,206.0) + +#Calculate radiative acceptance as ratio of recon+sel/generated for each detector +#these are root histograms. +for i,(sname, histo) in enumerate(invmass_histos.items()): + + #divide recon+sel by generated + ratio = invmass_histos[sname].Clone() + ratio.Divide(mcsamples[sname]) + + #fit radiative acceptance + fitparams, _ = signalProcessor.fit_plot_with_poly(ratio, specify_n=7, set_xrange=True, xrange=(30.0, 220.0)) + outfile[f'rad_acc_{sname}'] = ratio + print(sname,fitparams) + + #convert root histograms to numpy data for convenient plotting using mpl + (xvals, yvals, errors), (x_fit, y_fit) = cnv_root_to_np(ratio) + fits[sname] = (x_fit, y_fit) + + #plot + plt.errorbar(xvals, yvals, yerr=errors, linestyle='', marker='o', color=colors[i], label=sname) + plt.plot(x_fit, y_fit, linewidth=3.0, color=colors[i]) + +plt.legend(fontsize=20) +#plot the real radiative acceptance (includes beam) +#radacc_off = polynomial(-0.48922505, 0.073733061, -0.0043873158, 0.00013455495, -2.3630535e-06, 2.5402516e-08, -1.7090900e-10, 7.0355585e-13, -1.6215982e-15, 1.6032317e-18) +#plt.plot(xvals, radacc_off(xvals), label='rad+beam', marker='o', color='blue') + +#======================================================================================================================================= +# CALCULATE SYSTEMATIC UNCERTAINTY: Using the ratio of the nominal and misaligned radiative acceptance functions +# If the radiative acceptance increases with misalignment, that represents a decrease in expected signal (no systematic) +# If the radiative acceptance decreases with misalignment, that will boost expected signal and must be accounted for. +#======================================================================================================================================= + +#this is a subfigure of the figure above +plt.subplot(2,1,2) + +#calculate ratio of the two radiative acceptance fits +#if ratio < 1.0, apply systematic to expected signal +fit_ratio = fits['misaligned_v1'][1]/fits['nominal'][1] +xvalues = fits['nominal'][0] + +#plot the ratio +plt.plot(xvalues, fit_ratio, color='black', marker = '+', mew=5) +plt.axhline(y=1.0, linestyle='--', color='black') +plt.axhline(y=0.8, linestyle='--', color='black') +plt.xlim(20.0,206.) +plt.ylim(0.6,1.1) +plt.xlabel('A\' Invariant Mass [MeV]') +plt.ylabel('Systematic Uncertainty') +plt.savefig(f'{outdir}/radiative_acceptance_misalignment.png') + + +#fit the systematic uncertainty results and save the fit +sys_gr = r.TGraph(len(xvalues), xvalues, fit_ratio) +params_sys, errors_sys = signalProcessor.fit_plot_with_poly(sys_gr, tgraph=True, specify_n = 9, set_xrange=True, xrange=(50.0, 220.0)) +(xvals, yvals, errors), (x_fit, y_fit) = signalProcessor.cnv_tgraph_to_np(sys_gr) +fig, ax = plt.subplots(figsize=(20,10)) +plt.plot(xvals, yvals, marker='+', mew=3, markersize=10, color='darkblue') +plt.plot(x_fit, y_fit, linewidth=3.0, color='red') +outfile['misalignment_systematic'] = sys_gr + +#======================================================================================================================================= +# SIGNAL MISALIGNMENT: calculate using the radiative acceptance fits from above. +# 1. Use nominal signal (NO BEAM) and the nominal radiative acceptance (NO BEAM) to calculate expected signal. +# 2. Use misaligned signal (NO BEAM) and the misaligned radiative acceptance (NO BEAM) to calculate expected signal. +# 3. Calculate ratio between nominal and misaligned expected signal rates. +# *I've already done steps 1 and 2 externally, and am loading the results below. +#======================================================================================================================================= +#Load expected signal using nominal detector and nominal radiative acceptance (NO BEAM) +infile = '/sdf/group/hps/user-data/alspellm/2016/systematics/simp_radacc_misalignment_v1/simp_systematic_nominal.root' +with uproot.open(infile) as f: + nominal_h = f['expected_signal_ap_h'].to_hist() + +#Load expected signal using misaligned detector and misaligned radiative acceptance (NO BEAM) +infile = '/sdf/group/hps/user-data/alspellm/2016/systematics/simp_radacc_misalignment_v1/simp_systematic_misaligned.root' +with uproot.open(infile) as f: + misaligned_h = f['expected_signal_ap_h'].to_hist() + exp_sig_ratio_h = f['expected_signal_ap_h'].to_hist().reset() #make copy to use as ratio plot +outfile['expected_signal_nominal'] = nominal_h +outfile['expected_signal_misaligned'] = misaligned_h + +# Calculate expected signal ratio between nominal and misaligned +ratio = nominal_h.values()/misaligned_h.values() +xbins = exp_sig_ratio_h.axes[0].centers +ybins = exp_sig_ratio_h.axes[1].centers +xgrid, ygrid = np.meshgrid(xbins, ybins, indexing='ij') +exp_sig_ratio_h.fill(xgrid.flatten(), ygrid.flatten(), weight=ratio.flatten()) +outfile['signal_systematic_ratio'] = exp_sig_ratio_h + +#Plot ratio in 2d +fig, ax = plt.subplots(figsize=(25,15)) +exp_sig_ratio_h.plot(cmin=np.min(ratio.flatten()[ratio.flatten()>0.0]), cmax=np.max(ratio.flatten())) +plt.savefig(f'{outdir}/simp_radacc_misaligned_v1.png') + +#======================================================================================================================================= +# CALCULATE SYSTEMATIC UNCERTAINTY: Using the 2d expected signal plot of nominal/misaligned, decide how to get systematic. +#======================================================================================================================================= + +#Foddr each A' invariant mass, take the largest value, since this is in the numerator and reduces the expected signal rate +# The systematic uncertainty on the signal acceptance resulting from detector misaligned is calculated as a function of A' mass. +# For each mass, take the maximum ratio across the relevant range of epsilon^2 to be the uncertainty. +sigsys_y = [] +sigsys_x = [] +for xbin,mass in enumerate(exp_sig_ratio_h.axes[0].centers): + sigsys = np.max(exp_sig_ratio_h.values()[xbin]) + if sigsys == 0.0: + continue + sigsys_x.append(mass) + sigsys_y.append(sigsys) + +sigsys_gr = r.TGraph(len(sigsys_x), np.array(sigsys_x), np.array(sigsys_y)) +params_sigsys, errors_sigsys = fit_plot_with_poly(sigsys_gr, tgraph=True, specify_n = 5, set_xrange=True, xrange=(sigsys_x[0], sigsys_x[-1])) +print(params_sigsys) +(sigsys_x, sigsys_y, sigsys_errors), (sigsys_xfit, sigsys_yfit) = cnv_tgraph_to_np(sigsys_gr) +plt.plot(sigsys_x, sigsys_y) +plt.plot(sigsys_xfit, sigsys_yfit) +sigsys_final = np.max(sigsys_yfit) +print(f'Signal misalignment acceptance systematic: {sigsys_final}') + + +#Combine the signal and radiative acceptance systematics +radsys_fitpoly = polynomial(-10.307720, 0.97578691, -0.036585723, 0.00077903787, -1.0393704e-05, 9.0187487e-08, -5.0948313e-10, 1.8078746e-12, -3.6566050e-15, 3.2111742e-18) +masses = np.array([float(x) for x in range(60,230,1)]) +#Divide signal systematic by radiative acceptance +misalignmentsys = sigsys_final/radsys_fitpoly(masses) +fig, ax = plt.subplots(figsize=(25,15)) +plt.plot(masses, misalignmentsys, marker='+', markersize=10, mew=3, color='black') + diff --git a/plotUtils/simps/systematics/radiative_acceptance_target_position_systematic.py b/plotUtils/simps/systematics/radiative_acceptance_target_position_systematic.py new file mode 100644 index 000000000..e47ce397d --- /dev/null +++ b/plotUtils/simps/systematics/radiative_acceptance_target_position_systematic.py @@ -0,0 +1,167 @@ +#!/usr/bin/python3 +#======================================================================================================================================= +# Description: Calculates systematic uncertainty associated with target position uncertainty (0.5 mm) +# related to radiative trident acceptance. +# MC radiative tridents+beam at nominal, -0.5 mm, and +0.5 mm + +import os +import awkward as ak +import numpy as np +import hist +from hist import Hist +import uproot +import ROOT as r +import copy +import matplotlib.pyplot as plt +import matplotlib as mpl +import mplhep +import matplotlib.gridspec as gridspec +import sys +import math + +# SIMP tools in hpstr +hpstr_base = os.getenv('HPSTR_BASE') +sys.path.append(f'{hpstr_base}/plotUtils/simps') +import simp_signal_2016 + +#====================================================================================================================================== +#INITIALIZATION +#======================================================================================================================================= +# Set plotting parameters for matplotlib +plt.rcParams.update({'font.size': 40, 'axes.titlesize': 40, 'axes.labelsize': 40, 'xtick.labelsize': 40, 'ytick.labelsize': 40, 'lines.linewidth': 3.0, 'legend.fontsize': 40}) +plt.rcParams['font.family'] = 'DejaVu Sans' + +#parse input arguments +import argparse +parser = argparse.ArgumentParser(description='') +parser.add_argument('--outdir', type=str, default='./search_results') +parser.add_argument('--mpifpi', type=float, default=4.*np.pi) + +args = parser.parse_args() +outdir = args.outdir + +#======================================================================================================================================= +# LOAD MC RADIATIVE TRIDENTS FOR EACH TARGET POSITION +#======================================================================================================================================= + +# Invariant mass search window size +search_window = 1.5 +# Initialize signal processor that contains everything needed to calculate expected signal +signalProcessor = simp_signal_2016.SignalProcessor(args.mpifpi, search_window) + +samples = {} +mcsamples = {} +branches = ["unc_vtx_mass"] + +# Load nominal radiative tridents + beam MC +infile = '/sdf/group/hps/user-data/alspellm/2016/rad_mc/pass4b/rad_beam/rad-beam-hadd-10kfiles-ana-smeared-corr.root' +selection = 'vtxana_radMatchTight_nocuts' +samples['nominal_beam'] = signalProcessor.load_data(infile, selection, cut_expression='((unc_vtx_psum > 1.9) & (unc_vtx_psum < 2.4) )', expressions=branches) +#mc ana +infile = '/sdf/group/hps/user-data/alspellm/2016/rad_mc/pass4b/rad_nobeam/rad_nobeam_slic_hadd10ktuples_ana.root' +slicfile = r.TFile(infile, "READ") +mcsamples['nominal_beam'] = copy.deepcopy(slicfile.Get('mcAna/mcAna_mc622Mass_h')) +slicfile.Close() + +# Load nominal-0.5 mm radiative tridents + beam MC +infile = '/sdf/group/hps/user-data/alspellm/2016/systematics/radacc/hadd_1937files_rad_beam_targetz_Mpt5_recon_ana.root' +selection = 'vtxana_radMatchTight_nocuts' #USE RADMATCHTIGHT! +samples['targetz_Mpt5_beam'] = signalProcessor.load_data(infile, selection, cut_expression='((unc_vtx_psum > 1.9) & (unc_vtx_psum < 2.4) )', expressions=branches) +# Load radiative trident mc ana (generated rate) +infile = '/sdf/group/hps/user-data/alspellm/2016/systematics/radacc/hadd_1937files_rad_beam_targetz_Mpt5_mc_ana.root' +slicfile = r.TFile(infile, "READ") +mcsamples['targetz_Mpt5_beam'] = copy.deepcopy(slicfile.Get('mcAna/mcAna_mc622Mass_h')) +slicfile.Close() + +# Load nominal-0.5 mm radiative tridents + beam MC +infile = '/sdf/group/hps/user-data/alspellm/2016/systematics/radacc/hadd_1937files_rad_beam_targetz_Ppt5_recon_ana.root' +selection = 'vtxana_radMatchTight_nocuts' #USE RADMATCHTIGHT! +samples['targetz_Ppt5_beam'] = signalProcessor.load_data(infile, selection, cut_expression='((unc_vtx_psum > 1.9) & (unc_vtx_psum < 2.4) )', expressions=branches) +# Load radiative trident mc ana (generated rate) +infile = '/sdf/group/hps/user-data/alspellm/2016/systematics/radacc/hadd_1937files_rad_beam_targetz_Ppt5_mc_ana.root' +slicfile = r.TFile(infile, "READ") +mcsamples['targetz_Ppt5_beam'] = copy.deepcopy(slicfile.Get('mcAna/mcAna_mc622Mass_h')) +slicfile.Close() + +#======================================================================================================================================= +# CALCULATE RADIATIVE TRIDENT ACCEPTANCE +#======================================================================================================================================= + +# Initialize histogram used to calculate radiative trident acceptance +nbinsx = mcsamples['nominal_beam'].GetNbinsX() +first_bin = mcsamples['nominal_beam'].GetBinLowEdge(1) +last_bin = nbinsx*mcsamples['nominal_beam'].GetBinWidth(1) +invmass_h = ( + hist.Hist.new + .StrCategory(list(samples.keys()), name='samples') + .Reg(nbinsx,first_bin,last_bin,label='Invariant Mass [MeV]') + .Double() +) + +# Fill without weights, so that histos can be converted to ROOT and retain statistical uncertainty +invmass_histos = {} +for sname, sample in samples.items(): + invmass_h.fill(sname, sample.unc_vtx_mass*1000.) + invmass_histos[sname] = signalProcessor.cnvHistoToROOT(invmass_h[sname,:]) + invmass_histos[sname].Rebin(2) + mcsamples[sname].Rebin(2) + +def nonUniBinning(histo, start, size): + edges_a = np.arange(histo.GetBinLowEdge(1),start+histo.GetBinWidth(1),histo.GetBinWidth(1)) + edges_b = np.arange(start,histo.GetBinLowEdge(histo.GetNbinsX()), size) + bin_edges = np.concatenate([edges_a, edges_b[1:]]) + histo_rebinned = r.TH1F(f'{histo.GetName()}_rebinned', f'{histo.GetTitle()}', len(bin_edges)-1, bin_edges) + for bin in range(1, histo.GetNbinsX() + 1): + content = histo.GetBinContent(bin) + center = histo.GetBinCenter(bin) + error = histo.GetBinError(bin) + new_bin = histo_rebinned.FindBin(center) + histo_rebinned.SetBinContent(new_bin, histo_rebinned.GetBinContent(new_bin)+content) + histo_rebinned.SetBinError(new_bin, np.sqrt(histo_rebinned.GetBinError(new_bin)**2 + error**2)) + return histo_rebinned +#nonUniBinning(invmass_histos['nominal'], 150, 5) +#for sname, sample in samples.items(): +# invmass_histos[sname] = nonUniBinning(invmass_histos[sname], 150, 4) +# mcsamples[sname] = nonUniBinning(mcsamples[sname], 150, 4) + + +# Calculate radiative trident acceptance for each case (nominal, -0.5 mm, +0.5 mm) +fits = {} +colors = ['#d62728', '#bcbd22', '#2ca02c', '#17becf', '#1f77b4', '#9467bd', '#7f7f7f'] +colors = ['black', 'darkred', 'darkblue', 'darkgreen', 'darkorange'] +fig, ax = plt.subplots(2,1,figsize=(35,25),gridspec_kw={'height_ratios': [3, 2]}) +plt.subplot(2,1,1) +plt.xlabel('A\' Invariant Mass [MeV]') +plt.ylabel('Radiative Acceptance') +labels = ['Nominal (-4.3 mm)', '-4.8 mm', '-3.8 mm'] +for i,(sname, histo) in enumerate(invmass_histos.items()): + ratio = invmass_histos[sname].Clone() + ratio.Divide(mcsamples[sname]) # radiative trident acceptance + + # Fit the radiative trident acceptance ratio with a polynomial fit function + fit_params,_ = signalProcessor.fit_plot_with_poly(ratio, specify_n=7, set_xrange=True, xrange=(30.0, 220.0)) + print(sname, fit_params) + (xvals, yvals, errors), (x_fit, y_fit) = signalProcessor.cnv_root_to_np(ratio) + plt.errorbar(xvals, yvals, yerr=errors, linestyle='', marker='o', color=colors[i], label=labels[i]) + plt.plot(x_fit, y_fit, linewidth=3.0, color=colors[i]) + fits[sname] = (x_fit, y_fit) +plt.ylim(0.0, .15) +plt.xlim(30.0,220.0) +plt.legend(fontsize=50) +plt.text(53,0.125,'MC Radiative Tridents\n + Beam', horizontalalignment='center', fontsize=40) + +plt.subplot(2,1,2) +fit_ratio_Mpt5 = fits['targetz_Mpt5_beam'][1]/fits['nominal_beam'][1] +fit_ratio_Ppt5 = fits['targetz_Ppt5_beam'][1]/fits['nominal_beam'][1] +xvalues = fits['nominal_beam'][0] +plt.plot(xvalues, fit_ratio_Mpt5, color='darkred', marker='o', label='-4.8 mm : Nominal') +plt.plot(xvalues, fit_ratio_Ppt5, color='darkblue', marker='o', label='-3.8 mm : Nominal') +plt.axhline(y=1.0, linestyle='--', color='black') +plt.axhline(y=1.05, linestyle='--', color='black') +plt.xlim(30.0,220.) +plt.ylim(0.85,1.075) +plt.xlabel('A\' Invariant Mass [MeV]') +plt.ylabel('Ratio') +plt.legend() + +plt.savefig(f'{outdir}/radiative_acceptance_target_deltaz.png') diff --git a/plotUtils/simps/systematics/signal_misalignment_systematic.py b/plotUtils/simps/systematics/signal_misalignment_systematic.py new file mode 100644 index 000000000..af9515918 --- /dev/null +++ b/plotUtils/simps/systematics/signal_misalignment_systematic.py @@ -0,0 +1,84 @@ +#!/usr/bin/env python +# coding: utf-8 + +# In[1]: + + +import os +import awkward as ak +import numpy as np +import hist +from hist import Hist +import uproot +import ROOT as r +import copy + +import matplotlib.pyplot as plt +import matplotlib as mpl +import mplhep +import matplotlib.gridspec as gridspec + +import sys +sys.path.append('/sdf/group/hps/user-data/alspellm/2016/plotting') +import hps_plot_utils as utils + +get_ipython().run_line_magic('matplotlib', 'inline') +mpl.style.use(mplhep.style.ROOT) +import math +import pickle + +sys.path.append('/sdf/home/a/alspellm/src/hpstr_v62208/plotUtils/simps') +import simp_signal_2016 +from simp_theory_equations import SimpEquations as simpeqs +import copy +# Set global font sizes +plt.rcParams.update({'font.size': 40, # Font size for text + 'axes.titlesize': 40, # Font size for titles + 'axes.labelsize': 40, # Font size for axis labels + 'xtick.labelsize': 40, # Font size for x-axis tick labels + 'ytick.labelsize': 40, # Font size for y-axis tick labels + 'lines.linewidth':3.0, + 'legend.fontsize': 40}) # Font size for legend +plt.rcParams['font.family'] = 'DejaVu Sans' + + +# In[ ]: + + +infile = '/sdf/group/hps/user-data/alspellm/2016/systematics/simps/expected_signals/simp_systematic_nominal.root' +with uproot.open(infile) as f: + nominal_h = f['expected_signal_ap_h'].to_hist() + +infile = '/sdf/group/hps/user-data/alspellm/2016/systematics/simps/expected_signals/simp_systematic_misaligned.root' +with uproot.open(infile) as f: + misaligned_h = f['expected_signal_ap_h'].to_hist() + ratio_h = f['expected_signal_ap_h'].to_hist().reset() + +nominal_h.plot() +plt.show() +misaligned_h.plot() +plt.show() + +#take ratio of densities, misaligned to nominal +ratio = misaligned_h.values()/nominal_h.values() +# Find where nominal_h values are less than 1.0 +mask = nominal_h.values() < 1.0 +# Set corresponding ratio values to 0 +ratio[mask] = 0 + + +xbins = ratio_h.axes[0].centers +ybins = ratio_h.axes[1].centers +xgrid, ygrid = np.meshgrid(xbins, ybins, indexing='ij') +ratio_h.reset() +ratio_h.fill(xgrid.flatten(), ygrid.flatten(), weight=ratio.flatten()) +fig, ax = plt.subplots(figsize=(25,15)) +ratio_h.plot(cmin=0.0, cmax=np.max(ratio.flatten())) + + + +# In[ ]: + + +# + diff --git a/plotUtils/simps/systematics/signal_target_position_systematic.py b/plotUtils/simps/systematics/signal_target_position_systematic.py new file mode 100644 index 000000000..116b870f7 --- /dev/null +++ b/plotUtils/simps/systematics/signal_target_position_systematic.py @@ -0,0 +1,102 @@ +#!/usr/bin/python3 +#======================================================================================================================================= +# Description: Calculates systematic uncertainty associated with target position uncertainty (0.5 mm) +# related to MC signal acceptance +# MC signal (NO beam) at nominal, -0.5 mm, and +0.5 mm + +import os +import awkward as ak +import numpy as np +import hist +from hist import Hist +import uproot +import ROOT as r +import copy +import matplotlib.pyplot as plt +import matplotlib as mpl +import mplhep +import matplotlib.gridspec as gridspec +import sys +import math + +# SIMP tools in hpstr +hpstr_base = os.getenv('HPSTR_BASE') +sys.path.append(f'{hpstr_base}/plotUtils/simps') +import simp_signal_2016 + +# Load nominal +infile = '/sdf/group/hps/user-data/alspellm/2016/systematics/simp_target_unc/nov0proj/nominal.root' +with uproot.open(infile) as f: + nominal_h = f['expected_signal_ap_h'].to_hist() + test_h = f['expected_signal_vd_h'].to_hist() + +# Load -0.5 mm +infile = '/sdf/group/hps/user-data/alspellm/2016/systematics/simp_target_unc/nov0proj/target_Mpt5.root' +with uproot.open(infile) as f: + Mpt5_h = f['expected_signal_ap_h'].to_hist() + ratio_Mpt5_h = f['expected_signal_ap_h'].to_hist().reset() + +# Load +0.5 mm +infile = '/sdf/group/hps/user-data/alspellm/2016/systematics/simp_target_unc/nov0proj/target_Ppt5.root' +with uproot.open(infile) as f: + Ppt5_h = f['expected_signal_ap_h'].to_hist() + ratio_Ppt5_h = f['expected_signal_ap_h'].to_hist().reset() + + + + +#take ratio of densities, misaligned to nominal + +# Estimate systematic in a simple way by just taking the ratio of the expected signal rate between the two shifted target positions. +# This is probably too conservative, but what I did to finish my dissertation in time (Alic) + +ratio_PM = Ppt5_h.values()/Mpt5_h.values() + +xbins = ratio_Ppt5_h.axes[0].centers +ybins = ratio_Ppt5_h.axes[1].centers +xgrid, ygrid = np.meshgrid(xbins, ybins, indexing='ij') +ratio_Ppt5_h.reset() +ratio_Ppt5_h.fill(xgrid.flatten(), ygrid.flatten(), weight=ratio_PM.flatten()) +fig, ax = plt.subplots(figsize=(25,15)) +ratio_Ppt5_h.plot(cmin=0.90, cmax=1.1, cmap='seismic') +plt.text(124, -5.4,'Expected Signal Ratio\n Between Off-Nominal Targets' , horizontalalignment='center') +plt.ylim(-6.25, -4.7) +plt.xlim(50.0,210.0) +plt.ylabel(r'$\log{\epsilon^2}$', fontsize=50) +plt.xlabel('A\' Invariant Mass [MeV]') +plt.savefig('signal_target_uncertainty_offnominal_ratio_2d.png') + +# The systematic uncertainty is a function of both mass and epsilon. +# I decided to just take the worst case scenario for each MC mass across all relevent values of epsilon (where by relevent, I mean +# values of epsilon where we were able to put a 90% upper limit on the signal rate). +masses = [] +minvalues = [] +for m, mass in enumerate(ratio_Ppt5_h.axes[0].centers): + values = ratio_Ppt5_h.values()[m] + ybins = ratio_Ppt5_h.axes[1].centers + mask = np.where((ybins > -6.25) & (ybins < -4.7) ) + values_masked = values[mask] + minv = np.min(values_masked) + if not np.isfinite(minv): + continue + masses.append(mass) + minvalues.append(minv) + +# Fit the systematic uncertainty as a function of mass +coefficients = np.polyfit(masses, minvalues, 4) +print(coefficients) +fitfunc = np.poly1d(coefficients) +xfit = np.linspace(min(masses), max(masses),100) +yfit = fitfunc(xfit) + +fig, ax = plt.subplots(figsize=(25,15)) +plt.scatter(masses, minvalues, s=800, marker='+', color='black', label='Minimum Ratio') +plt.plot(xfit, yfit, color='darkred', label='4th-Order Fit') +plt.axhline(y=1.0, linestyle='--', color='black') +plt.ylim(0.5,1.2) +plt.xlim(50.0,210.) +plt.xlabel('A\' Invariant Mass [MeV]') +plt.ylabel('Expected Signal Ratio') +plt.legend() +plt.savefig('signal_target_uncertainty_offnominal_v2.png') + diff --git a/plotUtils/simps/systematics/v0projsig_systematic.py b/plotUtils/simps/systematics/v0projsig_systematic.py new file mode 100644 index 000000000..0a6cc596e --- /dev/null +++ b/plotUtils/simps/systematics/v0projsig_systematic.py @@ -0,0 +1,102 @@ +#!/usr/bin/python3 +#======================================================================================================================================= +# Description: Calculates systematic uncertainty associated with the target projected vertex significance cut. +# The systematic is calculated by comparing the resulting cut efficiency in data and MC background. + +import os +import awkward as ak +import numpy as np +import hist +from hist import Hist +import uproot +import ROOT as r +import copy +import matplotlib.pyplot as plt +import matplotlib as mpl +import mplhep +import matplotlib.gridspec as gridspec +import sys +import math + +# SIMP tools in hpstr +hpstr_base = os.getenv('HPSTR_BASE') +sys.path.append(f'{hpstr_base}/plotUtils/simps') +import simp_signal_2016 + +#====================================================================================================================================== +#INITIALIZATION +#======================================================================================================================================= +# Set plotting parameters for matplotlib +plt.rcParams.update({'font.size': 40, 'axes.titlesize': 40, 'axes.labelsize': 40, 'xtick.labelsize': 40, 'ytick.labelsize': 40, 'lines.linewidth': 3.0, 'legend.fontsize': 40}) +plt.rcParams['font.family'] = 'DejaVu Sans' + +#parse input arguments +import argparse +parser = argparse.ArgumentParser(description='') +parser.add_argument('--outdir', type=str, default='./search_results') +parser.add_argument('--mpifpi', type=float, default=4.*np.pi) + +args = parser.parse_args() +outdir = args.outdir + +#======================================================================================================================================= +# LOAD DATA AND MC BACKGROUND +#======================================================================================================================================= + +# Initialize the signal processor +search_window = 1.5 #used in final search +signalProcessor = simp_signal_2016.SignalProcessor(args.mpifpi, search_window) + +samples = {} +branches = ["unc_vtx_proj_sig","unc_vtx_ele_track_z0","unc_vtx_pos_track_z0"] + +# Load 10% data +infile = '/sdf/group/hps/user-data/alspellm/2016/data/hadd_BLPass4c_1959files.root' +selection = 'vtxana_Tight_L1L1_nvtx1' +# Select Psum signal region +samples['data'] = signalProcessor.load_data(infile,selection, expressions=branches, cut_expression='((unc_vtx_psum > 1.0) & (unc_vtx_psum < 1.9) )') +samples['data']['weight'] = 1.0 #Assign weight of 10 to scale up to full lumi + +# Load MC background +lumi = 10.7*.1 #pb-1 +mc_scale = {'tritrig' : 1.416e9*lumi/(50000*10000), + 'wab' : 0.1985e12*lumi/(100000*10000)} +#Load tritrig +infile = '/sdf/group/hps/user-data/alspellm/2016/tritrig_mc/pass4b/hadded_tritrig-beam-10kfiles-ana-smeared-corr_beamspotfix.root' +samples['tritrig'] = signalProcessor.load_data(infile, selection, cut_expression='((unc_vtx_psum > 1.2) & (unc_vtx_psum < 1.9) )', expressions=branches) +samples['tritrig']['weight'] = mc_scale['tritrig'] + +#Load wab +infile = '/sdf/group/hps/user-data/alspellm/2016/wab_mc/pass4b/hadded_wab-beam-10kfiles-ana-smeared-corr_beamspotfix.root' +samples['wab'] = signalProcessor.load_data(infile, selection, cut_expression='((unc_vtx_psum > 1.2) & (unc_vtx_psum < 1.9) )', expressions=branches) +samples['wab']['weight'] = mc_scale['wab'] + +#Combine tritrig and wab +samples['tritrig+wab+beam'] = ak.concatenate([samples['tritrig'], samples['wab']]) + +# Initialize histograms to calculate cut efficiencies +v0projsig_h = ( + hist.Hist.new + .StrCategory(list(samples.keys()), name='samples') + .Reg(300, 0.0,30,label=r'Target Projected Vertex Significance $[N\sigma_{\text{V0proj}}]$') + .Double() +) + +# Fill without weights, so that histos can be converted to ROOT and retain statistical uncertainty +for sname, sample in samples.items(): + v0projsig_h.fill(sname, sample.unc_vtx_proj_sig, weight=sample.weight/ak.sum(sample.weight)) + +# Calculate efficiencies +eff_mc = round(v0projsig_h['tritrig+wab+beam',:][:hist.loc(2.0):sum]/v0projsig_h['tritrig+wab+beam',:][::sum],2) +eff_data = round(v0projsig_h['data',:][:hist.loc(2.0):sum]/v0projsig_h['data',:][::sum],2) + +fig, ax = plt.subplots(figsize=(60,40)) +v0projsig_h['data',:].plot(linewidth=3.0, label='10% Data', color='black') +v0projsig_h['tritrig+wab+beam',:].plot(linewidth=3.0, label='Tritrig+WAB+Beam', color='blue') +plt.axvline(x=2.0, linestyle='--', color='red', label='Cut > 2.0') +plt.text(15.0, 3e-3, f'Data Eff: {eff_data}\nMC Bkg Eff: {eff_mc}') +plt.legend() +plt.ylabel('Normalized Events') +plt.yscale('log') +plt.savefig(f'{outdir}/v0projsig_systematic_lowpsum.png') + diff --git a/processors/config/anaFeeSmearing_cfg.py b/processors/config/anaFeeSmearing_cfg.py index e045112f4..2151d52c8 100644 --- a/processors/config/anaFeeSmearing_cfg.py +++ b/processors/config/anaFeeSmearing_cfg.py @@ -34,7 +34,7 @@ anaTrks.parameters["isData"] = options.isData #SmearingClosureTest -anaTrks.parameters["pSmearingFile"] = os.environ['HPSTR_BASE']+"/utils/data/smearingFile_2016_all_12112023.root" +anaTrks.parameters["pSmearingFile"] = os.environ['HPSTR_BASE']+"/utils/data/smearingFile_2016_all_20240620.root" RegionPath = os.environ['HPSTR_BASE']+"/analysis/selections/feeSmearing/" anaTrks.parameters["regionDefinitions"] = [] diff --git a/processors/config/anaSimps_2016_cfg.py b/processors/config/anaSimps_2016_cfg.py index 911b66881..956371794 100644 --- a/processors/config/anaSimps_2016_cfg.py +++ b/processors/config/anaSimps_2016_cfg.py @@ -7,6 +7,9 @@ help="Make True to make vertex ana flat tuple", metavar="makeFlatTuple", default=1) base.parser.add_argument("-r", "--isRadPDG", type=int, dest="isRadPDG", help="Set radiative trident PDG ID", metavar="isRadPDG", default=622) +base.parser.add_argument("--pSmearingSeed", type=int, dest="pSmearingSeed", + help="Set job dependent momentum smearing seed", metavar="pSmearingSeed", default=42) + options = base.parser.parse_args() # Use the input file to set the output file name @@ -30,7 +33,7 @@ # Processors # ############################### vtxana = HpstrConf.Processor('vtxana', 'NewVertexAnaProcessor') -mcana = HpstrConf.Processor('mcpartana', 'MCAnaProcessor') +#mcana = HpstrConf.Processor('mcpartana', 'MCAnaProcessor') ############################### # Processor Configuration # ############################### @@ -43,7 +46,9 @@ vtxana.parameters["hitColl"] = "SiClustersOnTrackOnPartOnUVtx" vtxana.parameters["analysis"] = "vertex" vtxana.parameters["vtxSelectionjson"] = os.environ['HPSTR_BASE']+"/analysis/selections/simps/vertexSelection_2016_simp_preselection.json" -vtxana.parameters["histoCfg"] = os.environ['HPSTR_BASE']+"/analysis/plotconfigs/tracking/simps/vtxAnalysis_2016_simp_reach.json" +#vtxana.parameters["vtxSelectionjson"] = os.environ['HPSTR_BASE']+"/analysis/selections/simps/vertexSelection_2016_simp_nocuts.json" +vtxana.parameters["histoCfg"] = os.environ['HPSTR_BASE']+"/analysis/plotconfigs/tracking/simps/vtxAnalysis_2016_simp_reach_light.json" +#vtxana.parameters["histoCfg"] = os.environ['HPSTR_BASE']+"/analysis/plotconfigs/tracking/simps/vtxAnalysis_2016_simp_reach.json" vtxana.parameters["mcHistoCfg"] = os.environ['HPSTR_BASE']+'/analysis/plotconfigs/mc/basicMC.json' ##### vtxana.parameters["beamE"] = base.beamE[str(options.year)] @@ -51,22 +56,15 @@ vtxana.parameters["isRadPDG"] = options.isRadPDG vtxana.parameters["makeFlatTuple"] = options.makeFlatTuple vtxana.parameters["beamPosCfg"] = "" -vtxana.parameters["pSmearingFile"] = os.environ['HPSTR_BASE']+"/utils/data/smearingFile_2016_all_12112023.root" -if options.isData: +if options.isData and options.year == 2016: vtxana.parameters["v0ProjectionFitsCfg"] = os.environ['HPSTR_BASE']+'/analysis/data/v0_projection_2016_config.json' -else: - vtxana.parameters["v0ProjectionFitsCfg"] = os.environ['HPSTR_BASE']+'/analysis/data/v0_projection_2016_mc_7800_config.json' #For tritrig and wab mc - #vtxana.parameters["v0ProjectionFitsCfg"] = os.environ['HPSTR_BASE']+'/analysis/data/v0_projection_2016_mc_signal_config.json' #For signal (accidentally gen with bspt=(0,0) - -if options.isData: - vtxana.parameters["eleTrackTimeBias"] = -1.5 - vtxana.parameters["posTrackTimeBias"] = -1.5 -else: - vtxana.parameters["eleTrackTimeBias"] = -2.2 #MC - vtxana.parameters["posTrackTimeBias"] = -2.2 #MC - #vtxana.parameters["eleTrackTimeBias"] = -5.5 #MC For TTs new smearing samples...due to readout bug - #vtxana.parameters["posTrackTimeBias"] = -5.5 #MC For TTs new smearing samples...due to readout bug - + vtxana.parameters["trackBiasCfg"] = os.environ['HPSTR_BASE']+'/analysis/data/track_bias_corrections_data_2016.json' +elif not options.isData and options.year == 2016: + print('Running MC') + vtxana.parameters["trackBiasCfg"] = os.environ['HPSTR_BASE']+'/analysis/data/track_bias_corrections_tritrig_2016.json' + vtxana.parameters["pSmearingFile"] = os.environ['HPSTR_BASE']+'/utils/data/fee_smearing/smearingFile_2016_all_20240620.root' + vtxana.parameters["pSmearingSeed"] = options.pSmearingSeed + vtxana.parameters["v0ProjectionFitsCfg"] = os.environ['HPSTR_BASE']+'/analysis/data/v0_projection_2016_mc_config.json' CalTimeOffset = -999 @@ -75,7 +73,7 @@ print("Running on data file: Setting CalTimeOffset %d" % CalTimeOffset) elif (options.isData == 0): - CalTimeOffset = 43. + CalTimeOffset = 42.4 print("Running on MC file: Setting CalTimeOffset %d" % CalTimeOffset) else: print("Specify which type of ntuple you are running on: -t 1 [for Data] / -t 0 [for MC]") @@ -85,22 +83,25 @@ RegionPath = os.environ['HPSTR_BASE']+"/analysis/selections/simps/" if (options.year == 2016): RegionDefinitions = [RegionPath+'Tight_2016_simp_reach_CR.json', - RegionPath+'Tight_2016_simp_reach_SR.json', + RegionPath+'Tight_2016_simp_SR_analysis.json', RegionPath+'Tight_nocuts.json', RegionPath+'Tight_L1L1_nvtx1.json'] if(options.isData != 1): RegionDefinitions.extend([RegionPath+'radMatchTight_2016_simp_reach_CR.json', - RegionPath+'radMatchTight_2016_simp_reach_SR.json']) + RegionPath+'radMatchTight_2016_simp_SR_analysis.json', + RegionPath+'radMatchTight_nocuts.json', + RegionPath+'radMatchTight_L1L1_nvtx1.json'] + ) vtxana.parameters["regionDefinitions"] = RegionDefinitions #MCParticleAna -mcana.parameters["debug"] = 0 -mcana.parameters["anaName"] = "mcAna" -mcana.parameters["partColl"] = "MCParticle" -mcana.parameters["trkrHitColl"] = "TrackerHits" -mcana.parameters["ecalHitColl"] = "EcalHits" -mcana.parameters["histCfg"] = os.environ['HPSTR_BASE']+'/analysis/plotconfigs/mc/basicMC.json' +#mcana.parameters["debug"] = 0 +#mcana.parameters["anaName"] = "mcAna" +#mcana.parameters["partColl"] = "MCParticle" +#mcana.parameters["trkrHitColl"] = "TrackerHits" +#mcana.parameters["ecalHitColl"] = "EcalHits" +#mcana.parameters["histCfg"] = os.environ['HPSTR_BASE']+'/analysis/plotconfigs/mc/basicMC.json' # Sequence which the processors will run. p.sequence = [vtxana] # ,mcana] diff --git a/processors/include/NewVertexAnaProcessor.h b/processors/include/NewVertexAnaProcessor.h index 0c2a5df45..bc08b679e 100644 --- a/processors/include/NewVertexAnaProcessor.h +++ b/processors/include/NewVertexAnaProcessor.h @@ -117,6 +117,7 @@ class NewVertexAnaProcessor : public Processor { TTree* tree_{nullptr}; //!< description std::string pSmearingFile_{""}; + int pSmearingSeed_{42}; std::shared_ptr smearingTool_; std::string pBiasingFile_{""}; std::shared_ptr biasingTool_; @@ -150,6 +151,9 @@ class NewVertexAnaProcessor : public Processor { std::vector beamPosCorrections_ = {0.0,0.0,0.0}; //!< holds beam position corrections std::string v0ProjectionFitsCfg_{""};//!< json file w run dependent v0 projection fits json v0proj_fits_;//!< json object v0proj + std::string trackBiasCfg_{""}; //!< json containing track bias corrections + json tbc_configs_; // trackBiasCorrections_; double eleTrackTimeBias_ = 0.0; double posTrackTimeBias_ = 0.0; diff --git a/processors/src/NewVertexAnaProcessor.cxx b/processors/src/NewVertexAnaProcessor.cxx index b042b248a..70121bb23 100644 --- a/processors/src/NewVertexAnaProcessor.cxx +++ b/processors/src/NewVertexAnaProcessor.cxx @@ -39,7 +39,8 @@ void NewVertexAnaProcessor::configure(const ParameterSet& parameters) { analysis_ = parameters.getString("analysis"); pSmearingFile_ = parameters.getString("pSmearingFile",pSmearingFile_); - pBiasingFile_ = parameters.getString("pBiasingFile",pBiasingFile_); + pSmearingSeed_ = parameters.getInteger("pSmearingSeed", pSmearingSeed_); + pBiasingFile_ = parameters.getString("pBiasingFile",pBiasingFile_); //region definitions regionSelections_ = parameters.getVString("regionDefinitions",regionSelections_); @@ -47,15 +48,18 @@ void NewVertexAnaProcessor::configure(const ParameterSet& parameters) { //v0 projection fits v0ProjectionFitsCfg_ = parameters.getString("v0ProjectionFitsCfg", v0ProjectionFitsCfg_); - //beamspot positions + //beamspot positions beamPosCfg_ = parameters.getString("beamPosCfg", beamPosCfg_); + + trackBiasCfg_ = parameters.getString("trackBiasCfg", trackBiasCfg_); + //track time bias corrections eleTrackTimeBias_ = parameters.getDouble("eleTrackTimeBias",eleTrackTimeBias_); posTrackTimeBias_ = parameters.getDouble("posTrackTimeBias",posTrackTimeBias_); - //bField scale factor (ONLY to correct for mistakes in ntuplizing). If < 0 is not used - bFieldScaleFactor_ = parameters.getDouble("bFieldScaleFactor",-1); - + //bField scale factor (ONLY to correct for mistakes in ntuplizing). If < 0 is not used + bFieldScaleFactor_ = parameters.getDouble("bFieldScaleFactor",-1); + } catch (std::runtime_error& error) { @@ -82,6 +86,14 @@ void NewVertexAnaProcessor::initialize(TTree* tree) { _mc_vtx_histos->Define2DHistos(); } + //Ana corrections to misc parameters + /* + if(!anaCorrectionsCfg_.empty()){ + std::ifstream anac_file(anaCorrectionsCfg_); + anac_file >> anac_configs_; + anac_file.close(); + }*/ + //Load Run Dependent V0 target projection fits from json if(!v0ProjectionFitsCfg_.empty()){ std::ifstream v0proj_file(v0ProjectionFitsCfg_); @@ -96,6 +108,16 @@ void NewVertexAnaProcessor::initialize(TTree* tree) { bpc_file >> bpc_configs_; bpc_file.close(); } + + //Track parameter bias corrections + if(!trackBiasCfg_.empty()){ + std::ifstream tbc_file(trackBiasCfg_); + tbc_file >> tbc_configs_; + tbc_file.close(); + } + + + // histos = new MCAnaHistos(anaName_); //histos->loadHistoConfig(histCfgFilename_) //histos->DefineHistos(); @@ -122,7 +144,7 @@ void NewVertexAnaProcessor::initialize(TTree* tree) { _reg_mc_vtx_histos[regname]->DefineHistos(); } - //Build a flat tuple for vertex and track params + //Build a flat tuple for vertex and track params if (makeFlatTuple_){ _reg_tuples[regname] = std::make_shared(anaName_+"_"+regname+"_tree"); @@ -208,10 +230,12 @@ void NewVertexAnaProcessor::initialize(TTree* tree) { //clust vars _reg_tuples[regname]->addVariable("unc_vtx_ele_clust_E"); _reg_tuples[regname]->addVariable("unc_vtx_ele_clust_x"); + _reg_tuples[regname]->addVariable("unc_vtx_ele_clust_y"); _reg_tuples[regname]->addVariable("unc_vtx_ele_clust_corr_t"); _reg_tuples[regname]->addVariable("unc_vtx_pos_clust_E"); _reg_tuples[regname]->addVariable("unc_vtx_pos_clust_x"); + _reg_tuples[regname]->addVariable("unc_vtx_pos_clust_y"); _reg_tuples[regname]->addVariable("unc_vtx_pos_clust_corr_t"); if(!isData_) @@ -252,20 +276,20 @@ void NewVertexAnaProcessor::initialize(TTree* tree) { if(!isData_ && !mcColl_.empty()) tree_->SetBranchAddress(mcColl_.c_str() , &mcParts_, &bmcParts_); if (not pSmearingFile_.empty()) { - // just using the same seed=42 for now - smearingTool_ = std::make_shared(pSmearingFile_,true); + smearingTool_ = std::make_shared(pSmearingFile_,true, pSmearingSeed_); } if (not pBiasingFile_.empty()) { - biasingTool_ = std::make_shared(pBiasingFile_); + biasingTool_ = std::make_shared(pBiasingFile_); } - + } bool NewVertexAnaProcessor::process(IEvent* ievent) { if(debug_) { std:: cout << "----------------- Event " << evth_->getEventNumber() << " -----------------" << std::endl; } + HpsEvent* hps_evt = (HpsEvent*) ievent; double weight = 1.; int run_number = evth_->getRunNumber(); @@ -280,9 +304,15 @@ bool NewVertexAnaProcessor::process(IEvent* ievent) { closest_run = check_run; } } - beamPosCorrections_ = {bpc_configs_[std::to_string(closest_run)]["beamspot_x"], - bpc_configs_[std::to_string(closest_run)]["beamspot_y"], - bpc_configs_[std::to_string(closest_run)]["beamspot_z"]}; + beamPosCorrections_ = {bpc_configs_[std::to_string(closest_run)]["unrotated_mean_x"], + bpc_configs_[std::to_string(closest_run)]["unrotated_mean_y"]}; + } + + //Load track parameter bias corrections if specified + if(!tbc_configs_.empty()){ + for(auto entry : tbc_configs_.items()){ + trackBiasCorrections_[entry.key()] = entry.value(); + } } @@ -337,7 +367,7 @@ bool NewVertexAnaProcessor::process(IEvent* ievent) { // Loop over vertices in event and make selections for ( int i_vtx = 0; i_vtx < vtxs_->size(); i_vtx++ ) { vtxSelector->getCutFlowHisto()->Fill(0.,weight); - + Vertex* vtx = vtxs_->at(i_vtx); Particle* ele = nullptr; Particle* pos = nullptr; @@ -347,9 +377,29 @@ bool NewVertexAnaProcessor::process(IEvent* ievent) { if (isData_) { if (!vtxSelector->passCutEq("Pair1_eq",(int)evth_->isPair1Trigger(),weight)) break; + if (!vtxSelector->passCutEq("Single0_eq",(int)evth_->isSingle0Trigger(),weight)) + break; } - bool foundParts = _ah->GetParticlesFromVtx(vtx,ele,pos); + //05072023 Test Cam's hypothesis. My understanding is that signal is injected every 500ns (250 events), so + //Triggered event time mod 500 should be ~0. + //MC beam background is inserted uniformly in all the events. + //Cam and Sarah find that, looking at triggered pulser data, the triggered events fall off in time away from event t mod 500 + //However, the MC beam is uniform across this distribution, insinuating that there are way more triggered backgrounds in MC + //than we actually expect in data...I think... + + if(!isData_){ + if (!vtxSelector->passCutLt("evtTimeMod500_lt", fabs(evth_->getEventTime()%500),weight)) + break; + } + + + bool foundParts = false; + if(vtxColl_ == "UnconstrainedMollerVertices" || vtxColl_ == "BeamspotConstrainedMollerVertices"){ + foundParts = _ah->GetSameParticlesFromVtx(vtx, ele, pos); + } + else + foundParts = _ah->GetParticlesFromVtx(vtx,ele,pos); if (!foundParts) { if(debug_) std::cout<<"NewVertexAnaProcessor::WARNING::Found vtx without ele/pos. Skip."<getTrack(); Track pos_trk = pos->getTrack(); - - if (debug_) { - std::cout<<"Check Ele/Pos Track momenta"< 0) { - biasingTool_->updateWithBiasP(ele_trk,bFieldScaleFactor_); - biasingTool_->updateWithBiasP(pos_trk,bFieldScaleFactor_); - } - - biasingTool_->updateWithBiasP(ele_trk); - biasingTool_->updateWithBiasP(pos_trk); - } - - if (debug_) { - std::cout<<"Corrected Ele/Pos Track momenta"<updateWithSmearP(ele_trk); - smearingTool_->updateWithSmearP(pos_trk); - double smeared_prod = ele_trk.getP()*pos_trk.getP(); - invm_smear = sqrt(smeared_prod/unsmeared_prod); + + if (debug_) { + std::cout<<"Check Ele/Pos Track momenta"< 0) { + biasingTool_->updateWithBiasP(ele_trk,bFieldScaleFactor_); + biasingTool_->updateWithBiasP(pos_trk,bFieldScaleFactor_); + } + + biasingTool_->updateWithBiasP(ele_trk); + biasingTool_->updateWithBiasP(pos_trk); + } + + if (debug_) { + std::cout<<"Corrected Ele/Pos Track momenta"<updateWithSmearP(ele_trk); + double pos_smf = smearingTool_->updateWithSmearP(pos_trk); + smearingTool_->updateVertexWithSmearP(vtx, ele_smf, pos_smf); } - + + //After all modifications to the track are made, update the memory address that the particle gets the track from + //with the modified track, that way changes here persist throughout this processor + ele->setTrack(&ele_trk); + pos->setTrack(&pos_trk); //Add the momenta to the tracks - do not do that //ele_trk.setMomentum(ele->getMomentum()[0],ele->getMomentum()[1],ele->getMomentum()[2]); @@ -584,7 +635,6 @@ bool NewVertexAnaProcessor::process(IEvent* ievent) { _vtx_histos->Fill1DHisto("pos_track_n2dhits_h", pos2dHits, weight); _vtx_histos->Fill1DHisto("vtx_Psum_h", p_ele.P()+p_pos.P(), weight); _vtx_histos->Fill1DHisto("vtx_Esum_h", ele_E + pos_E, weight); - _vtx_histos->Fill1DHisto("vtx_smear_InvM_h", invm_smear*(vtx->getInvMass()), weight); _vtx_histos->Fill1DHisto("ele_pos_clusTimeDiff_h", (corr_eleClusterTime - corr_posClusterTime), weight); _vtx_histos->Fill2DHisto("ele_vtxZ_iso_hh", TMath::Min(ele_trk.getIsolation(0), ele_trk.getIsolation(1)), vtx->getZ(), weight); _vtx_histos->Fill2DHisto("pos_vtxZ_iso_hh", TMath::Min(pos_trk.getIsolation(0), pos_trk.getIsolation(1)), vtx->getZ(), weight); @@ -644,7 +694,7 @@ bool NewVertexAnaProcessor::process(IEvent* ievent) { //TODO Clean this up => Cuts should be implemented in each region? //TODO Bring the preselection out of this stupid loop - + if (debug_) std::cout << "start regions" << std::endl; //TODO add yields. => Quite terrible way to loop. @@ -666,7 +716,11 @@ bool NewVertexAnaProcessor::process(IEvent* ievent) { Particle* ele = nullptr; Particle* pos = nullptr; - _ah->GetParticlesFromVtx(vtx,ele,pos); + if(vtxColl_ == "UnconstrainedMollerVertices" || vtxColl_ == "BeamspotConstrainedMollerVertices"){ + _ah->GetSameParticlesFromVtx(vtx, ele, pos); + } + else + _ah->GetParticlesFromVtx(vtx,ele,pos); CalCluster eleClus = ele->getCluster(); CalCluster posClus = pos->getCluster(); @@ -690,32 +744,16 @@ bool NewVertexAnaProcessor::process(IEvent* ievent) { Track ele_trk = ele->getTrack(); Track pos_trk = pos->getTrack(); - //Beam Position Corrections - ele_trk.applyCorrection("z0", beamPosCorrections_.at(1)); - pos_trk.applyCorrection("z0", beamPosCorrections_.at(1)); - //Track Time Corrections - ele_trk.applyCorrection("track_time",eleTrackTimeBias_); - pos_trk.applyCorrection("track_time", posTrackTimeBias_); - - if (biasingTool_) { - - //Correct the wrong Bfield first - if (bFieldScaleFactor_ > 0) { - biasingTool_->updateWithBiasP(ele_trk,bFieldScaleFactor_); - biasingTool_->updateWithBiasP(pos_trk,bFieldScaleFactor_); - } - - biasingTool_->updateWithBiasP(ele_trk); - biasingTool_->updateWithBiasP(pos_trk); - } - - double invm_smear = 1.; - if (smearingTool_) { - double unsmeared_prod = ele_trk.getP()*pos_trk.getP(); - smearingTool_->updateWithSmearP(ele_trk); - smearingTool_->updateWithSmearP(pos_trk); - double smeared_prod = ele_trk.getP()*pos_trk.getP(); - invm_smear = sqrt(smeared_prod/unsmeared_prod); + if (biasingTool_) { + + //Correct the wrong Bfield first + if (bFieldScaleFactor_ > 0) { + biasingTool_->updateWithBiasP(ele_trk,bFieldScaleFactor_); + biasingTool_->updateWithBiasP(pos_trk,bFieldScaleFactor_); + } + + biasingTool_->updateWithBiasP(ele_trk); + biasingTool_->updateWithBiasP(pos_trk); } //Add the momenta to the tracks @@ -1083,7 +1121,14 @@ bool NewVertexAnaProcessor::process(IEvent* ievent) { Particle* ele = nullptr; Particle* pos = nullptr; - if (!vtx || !_ah->GetParticlesFromVtx(vtx,ele,pos)) + bool foundParts = false; + if(vtxColl_ == "UnconstrainedMollerVertices" || vtxColl_ == "BeamspotConstrainedMollerVertices"){ + foundParts = _ah->GetSameParticlesFromVtx(vtx, ele, pos); + } + else + foundParts = _ah->GetParticlesFromVtx(vtx,ele,pos); + + if (!vtx || !foundParts) continue; CalCluster eleClus = ele->getCluster(); @@ -1099,40 +1144,27 @@ bool NewVertexAnaProcessor::process(IEvent* ievent) { Track ele_trk = ele->getTrack(); Track pos_trk = pos->getTrack(); //Get the shared info - TODO change and improve - - //Track Time Corrections - ele_trk.applyCorrection("track_time",eleTrackTimeBias_); - pos_trk.applyCorrection("track_time", posTrackTimeBias_); - - - // Track Momentum bias - - if (biasingTool_) { - - // Correct for wrong track momentum - Bug Fix - // In case there was mis-configuration during reco/hpstr-ntuple step, correct - // the momentum magnitude here using the right bField for the data taking year - - if (bFieldScaleFactor_ > 0) { - biasingTool_->updateWithBiasP(ele_trk,bFieldScaleFactor_); - biasingTool_->updateWithBiasP(pos_trk,bFieldScaleFactor_); - } - - - biasingTool_->updateWithBiasP(ele_trk); - biasingTool_->updateWithBiasP(pos_trk); - } - - - double invm_smear = 1.; - if (smearingTool_) { - double unsmeared_prod = ele_trk.getP()*pos_trk.getP(); - smearingTool_->updateWithSmearP(ele_trk); - smearingTool_->updateWithSmearP(pos_trk); - double smeared_prod = ele_trk.getP()*pos_trk.getP(); - invm_smear = sqrt(smeared_prod/unsmeared_prod); + + // Track Momentum bias + + if (biasingTool_) { + + // Correct for wrong track momentum - Bug Fix + // In case there was mis-configuration during reco/hpstr-ntuple step, correct + // the momentum magnitude here using the right bField for the data taking year + + if (bFieldScaleFactor_ > 0) { + biasingTool_->updateWithBiasP(ele_trk,bFieldScaleFactor_); + biasingTool_->updateWithBiasP(pos_trk,bFieldScaleFactor_); + } + + + biasingTool_->updateWithBiasP(ele_trk); + biasingTool_->updateWithBiasP(pos_trk); } + + double invm_smear = 1.; //Get the layers hit on each track std::vector ele_hit_layers = ele_trk.getHitLayers(); int ele_Si0 = 0; @@ -1280,7 +1312,7 @@ bool NewVertexAnaProcessor::process(IEvent* ievent) { //DeltaZ double deltaZ = std::abs( (ele_trk_z0/ele_trk.getTanLambda()) - (pos_trk_z0/pos_trk.getTanLambda()) ); - + //Project vertex to target double vtx_proj_x = -999.9; double vtx_proj_y = -999.9; @@ -1348,7 +1380,7 @@ bool NewVertexAnaProcessor::process(IEvent* ievent) { //1d histos _reg_vtx_histos[region]->Fill1DHisto("ele_track_clus_dt_h", ele_trk.getTrackTime() - corr_eleClusterTime, weight); _reg_vtx_histos[region]->Fill1DHisto("pos_track_clus_dt_h", pos_trk.getTrackTime() - corr_posClusterTime, weight); - + //TODO put this in the Vertex! TVector3 vtxPosSvt; @@ -1437,10 +1469,12 @@ bool NewVertexAnaProcessor::process(IEvent* ievent) { //clust vars _reg_tuples[region]->setVariableValue("unc_vtx_ele_clust_E", eleClus.getEnergy()); _reg_tuples[region]->setVariableValue("unc_vtx_ele_clust_x", eleClus.getPosition().at(0)); + _reg_tuples[region]->setVariableValue("unc_vtx_ele_clust_y", eleClus.getPosition().at(1)); _reg_tuples[region]->setVariableValue("unc_vtx_ele_clust_corr_t",corr_eleClusterTime); _reg_tuples[region]->setVariableValue("unc_vtx_pos_clust_E", posClus.getEnergy()); _reg_tuples[region]->setVariableValue("unc_vtx_pos_clust_x", posClus.getPosition().at(0)); + _reg_tuples[region]->setVariableValue("unc_vtx_pos_clust_y", posClus.getPosition().at(1)); _reg_tuples[region]->setVariableValue("unc_vtx_pos_clust_corr_t",corr_posClusterTime); _reg_tuples[region]->setVariableValue("run_number", evth_->getRunNumber()); diff --git a/processors/src/utilities.cxx b/processors/src/utilities.cxx index 4ba738858..c69e75975 100644 --- a/processors/src/utilities.cxx +++ b/processors/src/utilities.cxx @@ -296,7 +296,7 @@ Track* utils::buildTrack(EVENT::Track* lc_track, // Check that the TrackData data structure is correct. If it's // not, throw a runtime exception. - if (track_datum->getNDouble() > 14 || track_datum->getNFloat() > 7 || track_datum->getNInt() != 1) { + if (track_datum->getNDouble() > 14 || track_datum->getNFloat() > 8 || track_datum->getNInt() != 1) { throw std::runtime_error("[ TrackingProcessor ]: The collection " + std::string(Collections::TRACK_DATA) + " has the wrong structure."); diff --git a/utils/data/fee_smearing/fee_smearing_nhits_2016.py b/utils/data/fee_smearing/fee_smearing_nhits_2016.py new file mode 100644 index 000000000..846368986 --- /dev/null +++ b/utils/data/fee_smearing/fee_smearing_nhits_2016.py @@ -0,0 +1,157 @@ +#!/usr/bin/env python +import os +import numpy as np +import math +import ROOT as r + +def gaus_fit(histo, xmin, xmax, smean, swidth, snorm, nsigma=2.0, isData=False): + + #initial fit with seeds + fitfunc = r.TF1("gaus","gaus") + fitfunc.SetParameter(0, snorm) + fitfunc.SetParameter(1, smean) + fitfunc.SetParameter(2, swidth) + fitRes = histo.Fit(fitfunc,"QLES","", xmin, xmax) + params = fitRes.Parameters() + chi2 = fitRes.Chi2() + ndf = fitRes.Ndf() + + #set first fist to be best fit + best_chi2 = chi2/ndf + best_params = params + + #iterate over randomly fluctuated fit parameters. Keep the best resulting fit + niters = 100 + for n in range(niters): + norm = params[0]#*np.random.uniform(80,120)*0.01 + mu = params[1]#*np.random.uniform(80,120)*0.01 + sigma = params[2]#*np.random.uniform(80,120)*0.01 + + #Data has shoulders, so we can specify the xmin and xmax to do an asymmetric fit window + xminx = mu - nsigma*sigma + xmaxx = mu + nsigma*sigma + if isData: + if xminx < xmin: + xminx = xmin + if xmaxx > xmax: + xmaxx = xmax + fitfunc.SetParameter(0, norm) + fitfunc.SetParameter(1, mu) + fitfunc.SetParameter(2, sigma) + fitRes = histo.Fit(fitfunc,"QLES","", xminx, xmaxx) + #If fit fails, skip + try: + if fitRes.Parameters()[1] < xmin or fitRes.Parameters()[1] > xmax or fitRes.Ndf() < 1: + continue + except: + continue + + params = fitRes.Parameters() #these results seed the next iteration...maybe should only do if improved? + chi2 = fitRes.Chi2() + ndf = fitRes.Ndf() + #replace best fit + if chi2/ndf < best_chi2: + best_params = params + + #Do the final fit using the best parameters found + fitfunc.SetParameter(0, best_params[0]) + fitfunc.SetParameter(1, best_params[1]) + fitfunc.SetParameter(2, best_params[2]) + xminx = best_params[1] - nsigma*best_params[2] + xmaxx = best_params[1] + nsigma*best_params[2] + + #again, if data, use asymmetric fit window to avoid the left shoulder + if isData: + if xminx < xmin: + xminx = xmin + if xmaxx > xmax: + xmaxx = xmax + + fitRes = histo.Fit(fitfunc,"QLES","", xminx, xmaxx) + params = fitRes.Parameters() + errors = fitRes.Errors() + chi2 = fitRes.Chi2() + ndf = fitRes.Ndf() + return histo, params, errors, chi2/ndf + +#Load Data Run 7800 FEEs +data_results = {} +infilename = '/sdf/group/hps/user-data/alspellm/2016/fee_smearing/run7800/hadd/hadd_fee_2pt3_recon_fee_histos.root' #FEE skimmed track ana + +#Read track hit histograms +histoname = 'KalmanFullTracks/KalmanFullTracks_p_vs_nHits_top_hh' #top +infile = r.TFile(f'{infilename}',"READ") +top_h = copy.deepcopy(infile.Get(f'{histoname}')) +histoname = 'KalmanFullTracks/KalmanFullTracks_p_vs_nHits_bot_hh' #bot +bot_h = copy.deepcopy(infile.Get(f'{histoname}')) +infile.Close() + +#Change the names to use as keys +top_h.SetName('top') +bot_h.SetName('bot') + +#Fit the FEE peak for each category of nhits. Just have access to 10, 11, 12 for now +for h in [top_h, bot_h]: + histo = h + for nhits in [10, 11, 12]: + #Get the nhits momentum projection + proj = histo.ProjectionY(f'proj_{h.GetName()}_{nhits}hits', histo.GetXaxis().FindBin(nhits), histo.GetXaxis().FindBin(nhits),"") + #Fit the data + _, params, errors, chi2ndf = gaus_fit(proj, 2.0, 2.5, 2.4, 0.47, 12000, nsigma=1.5, isData=True) + + #store the results [mu,sigma] for top/bot nhits= + data_results[f'{h.GetName()}_nhits_{nhits}'] = [params[1], params[2]] + + +################################################################################## +#Read in MC FEE samples...Yes I could have combined these, but I separated them in development. Feel free to improve :) +mc_results = {} +infilename= '/sdf/group/hps/user-data/alspellm/2016/fee_smearing/tritrig/hadd/hadd_fee_2pt3_recon_tritrig_histos.root' + +#Read track hit histograms +histoname = 'KalmanFullTracks/KalmanFullTracks_p_vs_nHits_top_hh' #top +infile = r.TFile(f'{infilename}',"READ") +top_h = copy.deepcopy(infile.Get(f'{histoname}')) +histoname = 'KalmanFullTracks/KalmanFullTracks_p_vs_nHits_bot_hh' #bot +bot_h = copy.deepcopy(infile.Get(f'{histoname}')) +infile.Close() + +#Change the names to use as keys +top_h.SetName('top') +bot_h.SetName('bot') + +for h in [top_h, bot_h]: + histo = h + for nhits in [10, 11, 12]: + print(nhits) + proj = histo.ProjectionY(f'proj_{h.GetName()}_{nhits}hits', histo.GetXaxis().FindBin(nhits), histo.GetXaxis().FindBin(nhits),"") + _, params, errors, chi2ndf = gaus_fit(proj, 2.1, 2.5, 2.2, 0.1, proj.GetMaximum(), nsigma=1.5) + mc_results[f'{h.GetName()}_nhits_{nhits}'] = [params[1], params[2]] + +#Create output file to save results +outfile = r.TFile('smearingFile_2016_nhits.root',"RECREATE") +outfile.cd() +smtop_h = r.TH1F('KalmanFullTracks_p_vs_nHits_hh_smearing_rel_top','p_vs_nHits_smearing_rel_top;nhits;smear factor', 3, 9.5, 12.5) +smbot_h = r.TH1F('KalmanFullTracks_p_vs_nHits_hh_smearing_rel_bot','p_vs_nHits_smearing_rel_bot;nhits;smear factor', 3, 9.5, 12.5) + +#Calculate smearing factor according to 2016 Bump Hunt +smear_fac = lambda mu_data, sig_data, mu_mc, sig_mc : np.sqrt(np.square(sig_data/mu_data) - np.square(sig_mc/mu_mc)) +for key, vals in data_results.items(): + istop = False + if 'top' in key: + istop = True + nhits = float(key.split('_')[2]) + mu_data = vals[0] + sig_data = vals[1] + mu_mc = mc_results[key][0] + sig_mc = mc_results[key][1] + sf = smear_fac(mu_data, sig_data, mu_mc, sig_mc) + print(f'{key} sf={sf}') + + #save results + if istop: + smtop_h.SetBinContent(smtop_h.GetXaxis().FindBin(nhits), sf) + else: + smbot_h.SetBinContent(smbot_h.GetXaxis().FindBin(nhits), sf) + +outfile.Write() diff --git a/utils/data/fee_smearing/smearingFile_2016_all_20240620.root b/utils/data/fee_smearing/smearingFile_2016_all_20240620.root new file mode 100644 index 000000000..f0ec08b67 Binary files /dev/null and b/utils/data/fee_smearing/smearingFile_2016_all_20240620.root differ diff --git a/utils/data/smearingFile_2016_all_12112023.root b/utils/data/smearingFile_2016_all_12112023.root deleted file mode 100644 index 63031e483..000000000 Binary files a/utils/data/smearingFile_2016_all_12112023.root and /dev/null differ diff --git a/utils/include/TrackSmearingTool.h b/utils/include/TrackSmearingTool.h index 819be035e..c0a223af2 100644 --- a/utils/include/TrackSmearingTool.h +++ b/utils/include/TrackSmearingTool.h @@ -12,6 +12,8 @@ //------------------// #include "Track.h" +#include "Vertex.h" +#include "Particle.h" class TFile; class TH1D; @@ -27,7 +29,8 @@ class TrackSmearingTool { const std::string& tracks = "KalmanFullTracks"); double smearTrackP(const Track& trk); - void updateWithSmearP(Track& trk); + double updateWithSmearP(Track& trk); + void updateVertexWithSmearP(Vertex* vtx, double ele_smear_factor, double pos_smear_factor); private: diff --git a/utils/src/TrackSmearingTool.cxx b/utils/src/TrackSmearingTool.cxx index 769e1dedb..844f39dbc 100644 --- a/utils/src/TrackSmearingTool.cxx +++ b/utils/src/TrackSmearingTool.cxx @@ -5,85 +5,105 @@ #include TrackSmearingTool::TrackSmearingTool(const std::string& smearingfile, - const bool relSmearing, - const int seed, - const std::string& tracks){ - - - relSmearing_ = relSmearing; - std::string hsuffix = relSmearing_ ? "_rel" : ""; - smearingfile_ = std::make_shared(smearingfile.c_str()); - - if (!smearingfile_) - throw std::invalid_argument( "Provided input smearing file does not exists"); - - //cache the smearing histograms - smearing_histo_top_ = (TH1D*) smearingfile_->Get((tracks+"_p_vs_nHits_hh_smearing"+hsuffix).c_str()); - smearing_histo_bot_ = (TH1D*) smearingfile_->Get((tracks+"_p_vs_nHits_hh_smearing"+hsuffix).c_str()); - - if (!smearing_histo_top_ || !smearing_histo_bot_) - throw std::invalid_argument("Top and Bottom smearing histograms not found in smearing file"); - - //setup random engine - if (debug_) - std::cout<<"Setting up random engine with seed "<(seed); - - normal_ = std::make_shared>(0.,1.); - + const bool relSmearing, + const int seed, + const std::string& tracks){ + + + relSmearing_ = relSmearing; + std::string hsuffix = relSmearing_ ? "_rel" : ""; + smearingfile_ = std::make_shared(smearingfile.c_str()); + + if (!smearingfile_) + throw std::invalid_argument( "Provided input smearing file does not exists"); + + //cache the smearing histograms + smearing_histo_top_ = (TH1D*) smearingfile_->Get((tracks+"_p_vs_nHits_hh_smearing"+hsuffix+"_top").c_str()); + smearing_histo_bot_ = (TH1D*) smearingfile_->Get((tracks+"_p_vs_nHits_hh_smearing"+hsuffix+"_bot").c_str()); + + if (!smearing_histo_top_ || !smearing_histo_bot_) + throw std::invalid_argument("Top and Bottom smearing histograms not found in smearing file"); + + //setup random engine + if (debug_) + std::cout<<"Setting up random engine with seed "<(seed); + + normal_ = std::make_shared>(0.,1.); + } double TrackSmearingTool::smearTrackP(const Track& track) { - - double p = track.getP(); - double nhits = track.getTrackerHitCount(); - bool isTop = track.getTanLambda() > 0. ? true : false; - int binN = smearing_histo_top_->FindBin(nhits); - - if (debug_) - std::cout<<"Track nhits="< smearing_histo_top_->GetXaxis()->GetNbins()) { - throw std::invalid_argument("Bin not found in smearing histogram"); - } - - double rel_smear = (*normal_)(*generator_); - double sp = 0.; - - if (isTop) - sp = rel_smear * smearing_histo_top_->GetBinContent(binN); - else - sp = rel_smear * smearing_histo_bot_->GetBinContent(binN); - - double psmear = 0.; - - if (relSmearing_) - psmear = p + sp*p; - else - psmear = p + sp; - - - if (debug_) { - std::cout<<"Track isTop: "< smearing_histo_top_->GetXaxis()->GetNbins()) { + throw std::invalid_argument("Bin not found in smearing histogram"); + } + + double rel_smear = (*normal_)(*generator_); + double sp = 0.; + + if (isTop) + sp = rel_smear * smearing_histo_top_->GetBinContent(binN); + else + sp = rel_smear * smearing_histo_bot_->GetBinContent(binN); + + double psmear = 0.; + + if (relSmearing_) + psmear = p + sp*p; + else + psmear = p + sp; + + + if (debug_) { + std::cout<<"Track isTop: "<