add first version of using event weights for scaling etc

analyzers/dataframe/src/Analysis_FCChh.cc

@@ -924,6 +924,7 @@ AnalysisFCChh::get_tagged_jets(
   for (size_t pid_index = 0; pid_index < pids.size(); ++pid_index) {
 
     const auto tag =
+        // static_cast<unsigned>(tag_values[pid_index]);
         static_cast<unsigned>(tag_values[pids[pid_index].parameters_begin]);
 
     // std::cout << "Tag = " << tag << std::endl;
@@ -1145,7 +1146,8 @@ ROOT::VecOps::RVec<RecoParticlePair> AnalysisFCChh::getOSPairs(
   // FOR DEBUG?
   //  if (OS_pairs.size() > 1){
   //  	std::cout << "Number of possible OS pairs: " << OS_pairs.size() <<
-  //  std::endl; 	std::cout << "Build from: " << leptons_pos.size() << " pos, "
+  //  std::endl; 	std::cout << "Build from: " << leptons_pos.size() << "
+  //  pos, "
   //  << leptons_neg.size() << " neg." << std::endl;
   //  }
 
@@ -1570,8 +1572,8 @@ ROOT::VecOps::RVec<float> AnalysisFCChh::get_mT_pseudo(
 
 // 	//if one of the input particles is empty, just fill a default value of
 // -999 as mT 	if (particle_1.size() < 1 || particle_2.size() < 1||
-// MET_obj.size() < 1 ){ 		m_strans_vector.push_back(-999.); 		return
-// m_strans_vector;
+// MET_obj.size() < 1 ){ 		m_strans_vector.push_back(-999.);
+// return m_strans_vector;
 // 	}
 
 // 	TLorentzVector tlv_vis1 = getTLV_reco(particle_1.at(0));
@@ -1618,8 +1620,8 @@ ROOT::VecOps::RVec<float> AnalysisFCChh::get_mT_pseudo(
 
 // 	//if one of the input particles is empty, just fill a default value of
 // -999 as mT 	if (particle_1.size() < 1 || particle_2.size() < 1||
-// MET_obj.size() < 1 ){ 		m_strans_vector.push_back(-999.); 		return
-// m_strans_vector;
+// MET_obj.size() < 1 ){ 		m_strans_vector.push_back(-999.);
+// return m_strans_vector;
 // 	}
 
 // 	TLorentzVector tlv_vis1 = getTLV_reco(particle_1.at(0));
@@ -2496,7 +2498,7 @@ ROOT::VecOps::RVec<edm4hep::MCParticleData> AnalysisFCChh::getTruthTau(
           !isFromHiggsDirect(
               truth_part, parent_ids,
               truth_particles)) { //&& isFromHadron(truth_part, parent_ids,
-                                  //truth_particles) ) {
+                                  // truth_particles) ) {
         continue;
       }
       // select both from higgs or from hadrons
@@ -2541,7 +2543,7 @@ ROOT::VecOps::RVec<edm4hep::MCParticleData> AnalysisFCChh::getTruthTauLeps(
           !isFromHiggsDirect(
               truth_part, parent_ids,
               truth_particles)) { //&& isFromHadron(truth_part, parent_ids,
-                                  //truth_particles) ) {
+                                  // truth_particles) ) {
         continue;
       }
       // select both from higgs or from hadrons
@@ -2618,7 +2620,7 @@ ROOT::VecOps::RVec<edm4hep::MCParticleData> AnalysisFCChh::getTruthTauHads(
           !isFromHiggsDirect(
               truth_part, parent_ids,
               truth_particles)) { //&& isFromHadron(truth_part, parent_ids,
-                                  //truth_particles) ) {
+                                  // truth_particles) ) {
         continue;
       }
       // select both from higgs or from hadrons
@@ -3336,7 +3338,7 @@ AnalysisFCChh::find_reco_matches_exclusive(
     }
     // if (reco_match_vector.size() > 1){
     //	std::cout << "Warning in AnalysisFCChh::find_reco_matches() - Truth
-    //particle matched to more than one reco particle." << std::endl;
+    // particle matched to more than one reco particle." << std::endl;
     // }
 
     // check that the reco particle is not already in the out_vector

examples/FCChh/ggHH_bbyy/analysis_final.py

@@ -22,6 +22,9 @@
 
 saveTabular = True
 
+# Optional: Use weighted events
+do_weighted = True 
+
 ###Dictionnay of the list of cuts. The key is the name of the selection that will be added to the output file
 cutList = {
             "sel0_myy":"m_yy[0] > 100. && m_yy[0] < 180.",

examples/FCChh/ggHH_bbyy/analysis_stage1.py

@@ -48,6 +48,9 @@ def __init__(self, cmdline_args):
         # Optional: running on HTCondor, default is False
         # self.run_batch = False
 
+        # Optional: Use weighted events
+        self.do_weighted = True 
+
         # Optional: test file that is used if you run with the --test argument (fccanalysis run ./examples/FCChh/ggHH_bbyy/analysis_stage1.py --test)
         self.test_file = 'root://eospublic.cern.ch//eos/experiment/fcc/hh/' \
                          'tutorials/edm4hep_tutorial_data/' \
@@ -130,8 +133,8 @@ def analyzers(self, dframe):
             ########################################### APPLY PRE-SELECTION ########################################### 
 
             #require at least two b-jets and two photons, both with invariant masses compatible with the Higgs mass
-            .Filter("sel_bjets.size()>1")
-            .Filter("sel_gamma.size()>1") 
+            # .Filter("sel_bjets.size()>1")
+            # .Filter("sel_gamma.size()>1") 
             # .Filter("m_bb[0] < 200.") 
             # .Filter("m_bb[0] > 80.") 
             # .Filter("m_yy[0] < 180.") 
@@ -147,6 +150,7 @@ def output(self):
         Output variables which will be saved to output root file.
         '''
         branch_list = [
+            'weight',
             # Photons and H(yy) system:
             'ngamma', 'g1_pt', 'g2_pt', 'g1_eta', 'g2_eta', 'm_yy',
             # b-jets and H(bb) system:

python/process.py

@@ -7,6 +7,7 @@
 import json
 import glob
 import logging
+from typing import Optional
 import urllib.request
 import yaml  # type: ignore
 import ROOT  # type: ignore
@@ -30,6 +31,68 @@ def get_entries(inpath: str) -> int:
     return nevents
 
 
+def get_entries_sow(infilepath: str, nevents_max: Optional[int] = None, get_local: bool = True, weight_name: str = "EventHeader.weight") -> tuple[int, int, float, float]:
+    '''
+    Get number of original entries and number of actual entries in the file, as well as the sum of weights
+    '''
+
+    infile = ROOT.TFile.Open(infilepath)
+    infile.cd()
+
+    processEvents = 0
+    processSumOfWeights = 0.
+    try:
+        processEvents = infile.Get('eventsProcessed').GetVal()
+    except AttributeError:
+        print('----> Warning: Input file is missing information about ' # should these warnings be kept? 
+              'original number of events!')
+    try:
+        processSumOfWeights = infile.Get('SumOfWeights').GetVal()
+    except AttributeError:
+        print('----> Warning: Input file is missing information about '
+              'original sum of weights!')
+
+    if not get_local:
+        return processEvents, None, processSumOfWeights, None
+
+    eventsTTree = 0
+    sumOfWeightsTTree = 0.
+
+     # check for empty chunk (can this be improved? exception from RDF cannot be caught it seems?)
+    tree =infile.Get("events")
+    if not tree:
+        print("Tree not found in file", infilepath, " possibly empty chunk - continuing with next one.")
+        infile.Close()
+        return processEvents, eventsTTree, processSumOfWeights, sumOfWeightsTTree
+
+    try:
+
+         #use a RDF here too so the nevents restriction option can be imposed easily for the local events
+        rdf_tmp = ROOT.ROOT.RDataFrame("events", infilepath)
+
+        if nevents_max:
+            rdf_tmp = rdf_tmp.Range(0, nevents_max)
+
+        eventsTTree = rdf_tmp.Count().GetValue()
+
+        # eventsTTree = infile.Get("events").GetEntries()
+        ROOT.gROOT.SetBatch(True)
+        try:
+            # infile.Get("events").Draw('EventHeader.weight[0]>>histo')
+            # histo=ROOT.gDirectory.Get('histo')
+            histo = rdf_tmp.Histo1D(weight_name)
+            sumOfWeightsTTree=float(eventsTTree)*histo.GetMean()
+        except AttributeError:
+            print('----> Warning: Input file has no event weights.')
+    except AttributeError:
+        print('----> Error: Input file is missing events TTree! Probably empty chunk.')
+        infile.Close()
+
+    infile.Close()
+
+    return processEvents, eventsTTree, processSumOfWeights, sumOfWeightsTTree
+
+
 def get_process_info(process: str,
                      prod_tag: str,
                      input_dir: str) -> tuple[list[str], list[int]]: