Small fixes to the cross section code (#963)

* First completed code for cross section. To be tested

* Fixes for working hf_pt_spectrum

* Database for Nb method

* Add histogram for fractions

* Back to dd config

* Use configurable nevents for crosssec calculations

* Fix nevents for Nb method. Fix input for dd method

* Remove redundant hfraction. Add propagated statistical uncertainties to dd result.

* Narrow range for low pt reproduction

* Add dd_N method to calculate cross section from final corrected yields from cut variation

* Configure proper merged input files for dd in database and analyser

* Working cross section estimators. Only multiclass optimization remains

* Defer opening of nevents files

* Remove spurious lines

* Exclude changes in database and submission files

* Fixed next lines

Author: saganatt

machine_learning_hep/analysis/analyzerdhadrons.py

@@ -418,7 +418,7 @@ def efficiency(self):
         print(self.n_fileff)
         lfileeff = TFile.Open(self.n_fileff)
         lfileeff.ls()
-        fileouteff = TFile.Open(f"{self.d_resultsallpmc}/efficiencies{self.case}{self.typean}.root", "recreate")
+        fileouteff = TFile.Open(f"{self.d_resultsallpmc}/{self.efficiency_filename}{self.case}{self.typean}.root", "recreate")
         cEff = TCanvas("cEff", "The Fit Canvas")
         cEff.SetCanvasSize(1900, 1500)
         cEff.SetWindowSize(500, 500)
@@ -494,7 +494,7 @@ def makenormyields(self):  # pylint: disable=import-outside-toplevel, too-many-b
         if not os.path.exists(yield_filename):
             self.logger.fatal("Yield file %s could not be found", yield_filename)
 
-        fileouteff = f"{self.d_resultsallpmc}/efficiencies{self.case}{self.typean}.root"
+        fileouteff = f"{self.d_resultsallpmc}/{self.efficiency_filename}{self.case}{self.typean}.root"
         if not os.path.exists(fileouteff):
             self.logger.fatal("Efficiency file %s could not be found", fileouteff)
 
@@ -506,13 +506,12 @@ def makenormyields(self):  # pylint: disable=import-outside-toplevel, too-many-b
 
         histonorm = TH1F("histonorm", "histonorm", 1, 0, 1)
 
-        filemass = TFile.Open(self.n_filemass)
-        hevents = filemass.Get("all_events")
-        hselevents = filemass.Get("sel_events")
-
         if self.p_nevents is not None:
             selnorm = self.p_nevents
         else:
+            filemass = TFile.Open(self.n_filemass)
+            hevents = filemass.Get("all_events")
+            hselevents = filemass.Get("sel_events")
             norm, selnorm = self.calculate_norm(self.logger, hevents, hselevents)
             histonorm.SetBinContent(1, selnorm)
             self.logger.warning("Number of events %d", norm)
@@ -555,7 +554,9 @@ def makenormyields(self):  # pylint: disable=import-outside-toplevel, too-many-b
         f_fileoutcross = TFile.Open(fileoutcross)
         if f_fileoutcross:
             hcross = f_fileoutcross.Get("hptspectrum")
+            hcrossbr = f_fileoutcross.Get("hptspectrum_wo_br")
             fileoutcrosstot.cd()
             hcross.Write()
+            hcrossbr.Write()
         histonorm.Write()
         fileoutcrosstot.Close()

machine_learning_hep/hf_analysis_utils.py

@@ -268,8 +268,6 @@ def compute_fraction_dd(acc_eff_same, acc_eff_other, corryields_same, corryields
     yield_times_acceff_same = corryields_same * acc_eff_same
     yield_times_acceff_other = corryields_other * acc_eff_other
     frac_v = yield_times_acceff_same / (yield_times_acceff_same + yield_times_acceff_other)
-    # print(f"same yield times acceff: {yield_times_acceff_same} " \
-    #      f"other {yield_times_acceff_other} final frac: {frac_v}")
 
     denom = (yield_times_acceff_same + yield_times_acceff_other) ** 2
     der_same_same = (
@@ -279,8 +277,6 @@ def compute_fraction_dd(acc_eff_same, acc_eff_other, corryields_same, corryields
     unc = np.sqrt(
         der_same_same**2 * cov_same + der_same_other**2 * cov_other + 2 * der_same_same * der_same_other * cov_comb
     )
-    # print(f"denom {denom} der_same_same {der_same_same} der_same_other {der_same_other} " \
-    #      f"cov same {cov_same} cov other {cov_other} cov comb {cov_comb} final unc {unc}")
 
     return [frac_v, frac_v - unc, frac_v + unc]
 

machine_learning_hep/hf_pt_spectrum.py

@@ -106,24 +106,22 @@ def hf_pt_spectrum(
 
     histos = {}
 
-    infile_pred = TFile.Open(input_fonll_or_fdd_pred)
-    if frac_method in ("dd", "dd_N"):
-        histos["corryields_fdd"] = [infile_pred.Get("hCorrYieldsPrompt"), infile_pred.Get("hCorrYieldsNonPrompt")]
-        histos["covariances"] = [
-            infile_pred.Get("hCovPromptPrompt"),
-            infile_pred.Get("hCovNonPromptNonPrompt"),
-            infile_pred.Get("hCovPromptNonPrompt"),
-        ]
-    else:
-        histos["FONLL"] = {"prompt": {}, "nonprompt": {}}
-        for pred in ("central", "min", "max"):
-            histos["FONLL"]["nonprompt"][pred] = infile_pred.Get(f"{fonll_hist_name[channel]}fromBpred_{pred}_corr")
-            histos["FONLL"]["nonprompt"][pred].SetDirectory(0)
-            if frac_method == "fc":
-                histos["FONLL"]["prompt"][pred] = infile_pred.Get(f"{fonll_hist_name[channel]}pred_{pred}")
-                histos["FONLL"]["prompt"][pred].SetDirectory(0)
-
-    infile_pred.Close()
+    with TFile.Open(input_fonll_or_fdd_pred) as infile_pred:
+        if frac_method in ("dd", "dd_N"):
+            histos["corryields_fdd"] = [infile_pred.Get("hCorrYieldsPrompt"), infile_pred.Get("hCorrYieldsNonPrompt")]
+            histos["covariances"] = [
+                infile_pred.Get("hCovPromptPrompt"),
+                infile_pred.Get("hCovNonPromptNonPrompt"),
+                infile_pred.Get("hCovPromptNonPrompt"),
+            ]
+        else:
+            histos["FONLL"] = {"prompt": {}, "nonprompt": {}}
+            for pred in ("central", "min", "max"):
+                histos["FONLL"]["nonprompt"][pred] = infile_pred.Get(f"{fonll_hist_name[channel]}fromBpred_{pred}_corr")
+                histos["FONLL"]["nonprompt"][pred].SetDirectory(0)
+                if frac_method == "fc":
+                    histos["FONLL"]["prompt"][pred] = infile_pred.Get(f"{fonll_hist_name[channel]}pred_{pred}")
+                    histos["FONLL"]["prompt"][pred].SetDirectory(0)
 
     infile_rawy = TFile.Open(yield_filename)
     histos["rawyields"] = infile_rawy.Get(yield_histoname)
@@ -235,10 +233,6 @@ def hf_pt_spectrum(
                 histos["covariances"][1 - pnp_ind].GetBinContent(i_pt + 1),
                 histos["covariances"][2].GetBinContent(i_pt + 1),
             )
-            print("FRACTION HEREEEEEEEE: ")
-            print(histos["corryields_fdd"][pnp_ind].GetBinContent(i_pt + 1))
-            print(histos["corryields_fdd"][1 - pnp_ind].GetBinContent(i_pt + 1))
-            print("FRACTION HEREEEEEEEE: ")
         elif frac_method == "dd_N":
             pnp_ind = 0 if crosssec_prompt else 1
             frac = [histos["corryields_fdd"][pnp_ind].GetBinContent(i_pt + 1)] * 3