Add cross section calculation. Update databases. Plotting. Fixes

machine_learning_hep/analysis/analyzer_jets.py

@@ -95,7 +95,9 @@ def __init__(self, datap, case, typean, period):
         self.hfeeddown_det = {'mc': {}, 'data': {}}
         self.h_reflcorr = create_hist('h_reflcorr', ';p_{T}^{HF} (GeV/#it{c})', self.bins_candpt)
         self.n_events = {}
-        self.n_colls = {}
+        self.n_colls_read = {}
+        self.n_colls_tvx = {}
+        self.n_bcs_tvx = {}
 
         self.path_fig = Path(f'{os.path.expandvars(self.d_resultsallpdata)}/fig')
         for folder in ['qa', 'fit', 'roofit', 'sideband', 'signalextr', 'sidesub', 'sigextr', 'fd', 'uf', 'eff']:
@@ -145,9 +147,13 @@ def init(self):
                 if not histonorm:
                     self.logger.critical('histonorm not found')
                 self.n_events[mcordata] = histonorm.GetBinContent(1)
-                self.n_colls[mcordata] = histonorm.GetBinContent(2)
-                self.logger.info('Number of sampled collisions for %s: %g', mcordata, self.n_colls[mcordata])
+                self.n_colls_read[mcordata] = histonorm.GetBinContent(2)
+                self.n_colls_tvx[mcordata] = histonorm.GetBinContent(3)
+                self.n_bcs_tvx[mcordata] = histonorm.GetBinContent(4)
                 self.logger.debug('Number of selected events for %s: %d', mcordata, self.n_events[mcordata])
+                self.logger.info('Number of sampled collisions for %s: %g', mcordata, self.n_colls_read[mcordata])
+                self.logger.info('Number of TVX collisions for %s: %g', mcordata, self.n_colls_tvx[mcordata])
+                self.logger.info('Number of TVX BCs for %s: %g', mcordata, self.n_bcs_tvx[mcordata])
 
     def qa(self): # pylint: disable=invalid-name
         self.logger.info("Producing basic QA histograms")
@@ -589,7 +595,7 @@ def _subtract_sideband(self, hist, var, mcordata, ipt):
         subtract_sidebands = False
         if mcordata == 'data' and self.cfg('sidesub_per_ptjet'):
             self.logger.info('Subtracting sidebands in pt jet bins')
-            for iptjet in range(1, get_nbins(fh_subtracted, 0)):
+            for iptjet in range(get_nbins(fh_subtracted, 0)):
                 if rws := self.roo_ws_ptjet[mcordata][iptjet][ipt]:
                     f = rws.pdf("bkg").asTF(self.roo_ws[mcordata][ipt].var("m"))
                 else:
@@ -617,9 +623,9 @@ def _subtract_sideband(self, hist, var, mcordata, ipt):
                 areaNormFactor = area['signal'] / (area['sideband_left'] + area['sideband_right'])
                 fh_sideband.Scale(areaNormFactor)
 
+        self._save_hist(fh_sideband,
+                        f'sideband/h_ptjet{label}_sideband_{string_range_pthf(range_pthf)}_{mcordata}.png')
         if subtract_sidebands:
-            self._save_hist(fh_sideband,
-                            f'sideband/h_ptjet{label}_sideband_{string_range_pthf(range_pthf)}_{mcordata}.png')
             fh_subtracted.Add(fh_sideband, -1.)
 
         self._clip_neg(fh_subtracted)
@@ -823,7 +829,7 @@ def _analyze(self, method = 'sidesub'):
                                 hproj = project_hist(h, [1], {0: (j+1, j+1)})
                                 empty = hproj.Integral() < 1.e-7
                                 if empty and i == 0:
-                                    self.logger.error("Projection %s %s is empty.", mcordata,
+                                    self.logger.error("Projection %s %s %s is empty.", var, mcordata,
                                                       string_range_ptjet(range_ptjet))
                                 self._save_hist(
                                     hproj,
@@ -916,9 +922,7 @@ def estimate_feeddown(self):
         with TFile(self.cfg('fd_root')) as rfile:
             powheg_xsection = rfile.Get('fHistXsection')
             powheg_xsection_scale_factor = powheg_xsection.GetBinContent(1) / powheg_xsection.GetEntries()
-        self.logger.info('powheg scale factor %g', powheg_xsection_scale_factor)
-        self.logger.info('number of collisions in data: %g', self.n_colls['data'])
-        self.logger.info('number of collisions in MC: %g', self.n_colls['mc'])
+        self.logger.info('POWHEG luminosity (mb^{-1}): %g', 1. / powheg_xsection_scale_factor)
 
         df = pd.read_parquet(self.cfg('fd_parquet'))
         col_mapping = {'dr': 'delta_r_jet', 'zpar': 'z'} # TODO: check mapping
@@ -1047,11 +1051,17 @@ def estimate_feeddown(self):
             hfeeddown_det.Scale(powheg_xsection_scale_factor * self.cfg('branching_ratio'))
             hfeeddown_det_mc = hfeeddown_det.Clone()
             hfeeddown_det_mc.SetName(hfeeddown_det_mc.GetName() + '_mc')
-            hfeeddown_det.Scale(self.n_colls['data'] / self.cfg('xsection_inel'))
-            hfeeddown_det_mc.Scale(self.n_colls['mc'] / self.cfg('xsection_inel_mc'))
-
-            self._save_hist(hfeeddown_det, f'fd/h_ptjet-{var}_feeddown_det_final.png')
-            self._save_hist(hfeeddown_det, f'fd/h_ptjet-{var}_feeddown_det_final_mc.png')
+            luminosity_data = (self.n_colls_read['data'] / self.n_colls_tvx['data'] *
+                               self.n_bcs_tvx['data'] / self.cfg('xsection_inel'))
+            self.logger.info("Scaling feed-down with data luminosity (mb^{-1}): %g", luminosity_data)
+            hfeeddown_det.Scale(luminosity_data)
+            luminosity_mc = (self.n_colls_read['mc'] / self.n_colls_tvx['mc'] *
+                             self.n_bcs_tvx['mc'] / self.cfg('xsection_inel') * self.cfg('lumi_scale_mc'))
+            self.logger.info("Scaling feed-down with MC luminosity (mb^{-1}): %g", luminosity_mc)
+            hfeeddown_det_mc.Scale(luminosity_mc)
+
+            self._save_hist(hfeeddown_det, f'fd/h_ptjet-{var}_feeddown_det_final_data.png')
+            self._save_hist(hfeeddown_det_mc, f'fd/h_ptjet-{var}_feeddown_det_final_mc.png')
             self.hfeeddown_det['data'][var] = hfeeddown_det
             self.hfeeddown_det['mc'][var] = hfeeddown_det_mc
 

machine_learning_hep/data/data_run3/database_ml_parameters_D0Jet_pp.yml

@@ -13,7 +13,7 @@
 
 D0Jet_pp:
   doml: true
-  mass: 1.864
+  mass: 1.86484
   sel_reco_unp: "fPt > 1."
   sel_reco_singletrac_unp: null
   sel_gen_unp: "fPt > 1. and fPt < 10."
@@ -76,7 +76,7 @@ D0Jet_pp:
           O2hfd0pbase: [fIndexHfD0McCollBases, fPt, fEta, fPhi, fFlagMcMatchGen, fOriginMcGen]
           O2d0cmcpjeto: [fIndexD0CMCPJETCOS, fJetPt, fJetPhi, fJetEta, fJetNConstituents, fJetR]
           O2d0cmcpjetmo: [fIndexArrayD0CMCDJETOS_hf, fIndexArrayD0CMCDJETOS_geo, fIndexArrayD0CMCDJETOS_pt]
-          O2d0cmcpjetsso: [fEnergyMother, fPtLeading, fPtSubLeading, fTheta, fNSub2DR, fNSub1, fNSub2]
+          O2d0cmcpjetsso: [fEnergyMother, fPtLeading, fPtSubLeading, fTheta, fNSub2DR, fNSub1, fNSub2, fAngularity, fPairPt, fPairEnergy, fPairTheta]
         extra:
           fY: log((sqrt(1.864**2 + (fPt * cosh(fEta))**2) + fPt * sinh(fEta)) / sqrt(1.864**2 + fPt**2)) #TODO : change mass or make sure Lc mass is updated
         tags:
@@ -110,7 +110,7 @@ D0Jet_pp:
           O2hfd0ml: [fMlScores]
           O2d0cmcdjeto: [fIndexD0CMCDJETCOS, fJetPt, fJetPhi, fJetEta, fJetNConstituents, fJetR]
           O2d0cmcdjetmo: [fIndexArrayD0CMCPJETOS_hf, fIndexArrayD0CMCPJETOS_geo, fIndexArrayD0CMCPJETOS_pt]
-          O2d0cmcdjetsso: [fEnergyMother, fPtLeading, fPtSubLeading, fTheta, fNSub2DR, fNSub1, fNSub2]
+          O2d0cmcdjetsso: [fEnergyMother, fPtLeading, fPtSubLeading, fTheta, fNSub2DR, fNSub1, fNSub2, fAngularity, fPairPt, fPairEnergy, fPairTheta]
         extra:
           fY: log((sqrt(1.864**2 + (fPt * cosh(fEta))**2) + fPt * sinh(fEta)) / sqrt(1.864**2 + fPt**2)) #TODO : change mass or make sure Lc mass is updated
           isd0: fFlagMcMatchRec == 1
@@ -150,7 +150,7 @@ D0Jet_pp:
           O2hfd0ml: [fMlScores]
           O2hfd0sel: [fCandidateSelFlag]
           O2d0cjeto: [fIndexD0CJETCOS, fJetPt, fJetPhi, fJetEta, fJetNConstituents, fJetR]
-          O2d0cjetsso: [fEnergyMother, fPtLeading, fPtSubLeading, fTheta, fNSub2DR, fNSub1, fNSub2]
+          O2d0cjetsso: [fEnergyMother, fPtLeading, fPtSubLeading, fTheta, fNSub2DR, fNSub1, fNSub2, fAngularity, fPairPt, fPairEnergy, fPairTheta]
         extra:
           fY: log((sqrt(1.864**2 + (fPt * cosh(fEta))**2) + fPt * sinh(fEta)) / sqrt(1.864**2 + fPt**2)) #TODO : change mass or make sure Lc mass is updated
         extract_component:
@@ -182,6 +182,9 @@ D0Jet_pp:
       collcnt:
         level: all
         file: AnalysisResultsCollCnt.parquet
+      bccnt:
+        level: all
+        file: AnalysisResultsBcCnt.parquet
 
   variables:
     var_all: [fDecayLength, fDecayLengthXY, fDecayLengthNormalised, fDecayLengthXYNormalised, fCpa, fCpaXY, fImpactParameter0, fImpactParameter1, fErrorImpactParameter0, fErrorImpactParameter1, fNSigTpcPiExpPi, fNSigTpcKaExpPi, fNSigTpcPiExpKa, fNSigTpcKaExpKa]
@@ -232,6 +235,7 @@ D0Jet_pp:
     namefile_reco: AnalysisResultsReco.parquet
     namefile_evt: AnalysisResultsEvt.parquet
     namefile_collcnt: AnalysisResultsCollCnt.parquet
+    namefile_bccnt: AnalysisResultsBcCnt.parquet
     namefile_evtvalroot: AnalysisResultsROOTEvtVal.root
     namefile_evtorig: AnalysisResultsEvtOrig.parquet
     namefile_gen: AnalysisResultsGen.parquet
@@ -375,13 +379,15 @@ D0Jet_pp:
     jet_obs: &jet_default
       sel_an_binmin:    [1,2,3,4,5,6,7,8,10,12,16,24] # hadron pt bins (sel_an_binmin bins)
       sel_an_binmax:    [2,3,4,5,6,7,8,10,12,16,24,36] # hadron pt bins (sel_an_binmin bins) # FIXME: move the last edge in sel_an_binmin
-      bins_ptjet: [5, 7, 15, 30, 50] # systematics
-      # TODO: split rec and gen binning
-      bins_ptjet_eff: [2, 5, 7, 15, 30, 50, 70]
-
+      bins_ptjet: [5, 7, 15, 30, 50] # systematics, TODO: split rec and gen binning
+      bins_ptjet_eff: [2, 5, 7, 15, 30, 50, 70] # systematics, TODO: split rec and gen binning
       cand_collidx: fIndexHfD0CollBases
-      # sigmamb: 57.8e-3 #NB: multiplied by 1e12 before giving to HFPtSpectrum!
-      cnt_events_read: fReadCountsWithTVXAndZVertexAndSel8
+      counter_read_data: fReadCountsWithTVXAndZVertexAndSel8
+      counter_read_mc: fReadCountsWithTVXAndZVertexAndSelMC
+      counter_tvx: fReadCountsWithTVX
+      xsection_inel: 59.4 # (mb) cross-section of minimum-bias events # used # systematics
+      lumi_scale_mc: 408  # charm enhancement factor in MC to scale the MC luminosity
+      branching_ratio: 3.947e-2 # used
 
       observables:
         zg:
@@ -696,10 +702,6 @@ D0Jet_pp:
       fd_root: '/data2/vkucera/powheg/trees_powheg_fd_central.root' # systematics
       fd_parquet: '/data2/jklein/powheg/trees_powheg_fd_central.parquet' # systematics
 
-      branching_ratio: 3.95e-2 # used
-      xsection_inel: 57.8 # (mb) cross-section of minimum-bias events # used # systematics
-      xsection_inel_mc: .65 # 1.68 # (mb) cross-section of minimum-bias events # used
-
       # obsolete?
       proc_type: Jets # used
       useperiod: [1] #list of periods # used
@@ -787,7 +789,6 @@ D0Jet_pp:
       # domodeldep: false
       # path_modeldep: /home/nzardosh/PYTHIA_Sim/PYTHIA8_Simulations/Plots/D0_Substructure_Simulations_Output.root
 
-
       # replace with fd_root...
       # powheg_path_nonprompt: /data/POWHEG/trees_powheg_fd_central.root
       # powheg_path_prompt: /data/POWHEG/trees_powheg_pr_central.root

machine_learning_hep/data/data_run3/database_ml_parameters_LcJet_pp.yml

@@ -13,7 +13,7 @@
 
 LcJet_pp:
   doml: true
-  mass: 2.286
+  mass: 2.28646
   sel_reco_unp: "fPt>0"
   sel_gen_unp: "fPt>0"
   sel_cen_unp: null
@@ -169,6 +169,9 @@ LcJet_pp:
       collcnt:
         level: all
         file: AnalysisResultsCollCnt.parquet
+      bccnt:
+        level: all
+        file: AnalysisResultsBcCnt.parquet
 
   variables:
     var_all: [fIndexCollisions, fFlagMcMatchRec, fCandidateSelFlag, fOriginMcRec, fIsCandidateSwapped, fNProngsContributorsPV,
@@ -222,6 +225,7 @@ LcJet_pp:
     namefile_reco: AnalysisResultsReco.parquet
     namefile_evt: AnalysisResultsEvt.parquet
     namefile_collcnt: AnalysisResultsCollCnt.parquet
+    namefile_bccnt: AnalysisResultsBcCnt.parquet
     namefile_evtvalroot: AnalysisResultsROOTEvtVal.root
     namefile_evtorig: AnalysisResultsEvtOrig.parquet
     namefile_gen: AnalysisResultsGen.parquet
@@ -310,7 +314,7 @@ LcJet_pp:
       FF: 0.204 # fragmentation fraction
       sigma_MB: 57.8e-3  # Minimum Bias cross section (pp) 50.87e-3 [b], 1 for Pb-Pb
       Taa: 1 # 23260 [b^-1] in 0-10% Pb-Pb, 3917 [b^-1] in 30-50% Pb-Pb, 1 for pp
-      BR: 6.23e-2 # branching ratio of the decay Lc -> p K- pi+
+      BR: 6.24e-2 # branching ratio of the decay Lc -> p K- pi+
       f_prompt: 0.9 # estimated fraction of prompt candidates
       bkg_data_fraction: 0.05 # fraction of real data used in the estimation
       num_steps: 111 # number of steps used in efficiency and signif. estimation
@@ -353,11 +357,15 @@ LcJet_pp:
     jet_obs: &jet_default
       sel_an_binmin:    [3,4,5,6,7,8,10,12,16] # hadron pt bins (sel_an_binmin bins)
       sel_an_binmax:      [4,5,6,7,8,10,12,16,24] # hadron pt bins (sel_an_binmin bins)
-      # binning_matching:   [0, 1, 1, 2, 2, 3, 3, 4, 4, 6, 6] # mapping to skimming pt bins (sel_an_binmin bins)
-      bins_ptjet: [5, 7, 15, 30, 50]
-      bins_ptjet_eff: [2, 5, 7, 15, 30, 50, 70]
+      bins_ptjet: [5, 7, 15, 30, 50] # systematics, TODO: split rec and gen binning
+      bins_ptjet_eff: [2, 5, 7, 15, 30, 50, 70] # systematics, TODO: split rec and gen binning
       cand_collidx: fIndexHf3PCollBases
-      cnt_events_read: fReadCountsWithTVXAndZVertexAndSel8
+      counter_read_data: fReadCountsWithTVXAndZVertexAndSel8
+      counter_read_mc: fReadCountsWithTVXAndZVertexAndSelMC
+      counter_tvx: fReadCountsWithTVX
+      xsection_inel: 59.4 # (mb) cross-section of minimum-bias events # used # systematics
+      lumi_scale_mc: 408  # charm enhancement factor in MC to scale the MC luminosity
+      branching_ratio: 6.24e-2 # used
 
       observables:
         zg:
@@ -553,10 +561,6 @@ LcJet_pp:
       niterunfolding: 15
       niterunfoldingchosen: 4
 
-      branching_ratio: 3.95e-2
-      xsection_inel: 57.8 # (mb) cross-section of minimum-bias events
-      xsection_inel_mc: 57.8 # (mb) cross-section of minimum-bias events
-
       doprior: false
       domodeldep: false
       path_modeldep: /home/nzardosh/PYTHIA_Sim/PYTHIA8_Simulations/Plots/D0_Substructure_Simulations_Output.root

machine_learning_hep/data/data_run3/database_variations_D0Jet_pp_jet_obs.yml

@@ -817,7 +817,7 @@ categories:
                 diffs:
                     analysis:
                         jet_obs:
-                            xsection_inel: [54.91, 60.69]  # TODO: update values for Run 3
+                            xsection_inel: [53.46, 65.34]
     tracking:
         activate: yes
         processor: true

machine_learning_hep/data/data_run3/database_variations_LcJet_pp_jet_obs.yml

@@ -460,7 +460,7 @@ categories:
                 diffs:
                     analysis:
                         jet_obs:
-                            xsection_inel: [54.91, 60.69]  # TODO: update values for Run 3
+                            xsection_inel: [53.46, 65.34]
     tracking:
         activate: yes
         processor: true