Pulling from newly added structures in ReconstructedParticle.h

hshaddix · April 2, 2023, 2:22am

I noticed a while back that the lorentz vector information is forgotten for the particles that are used to meet the resonanceBuilder requirement in ReconstructedParticle.cc. In an attempt to add those back in, I have added a structure in the .h file:

struct ResoWithLegs;
struct resonanceBuilder {
float m_resonance_mass;
resonanceBuilder(float arg_resonance_mass);
ROOT::VecOps::RVecFCCAnalyses::ReconstructedParticle::ResoWithLegs operator()(ROOT::VecOps::RVecedm4hep::ReconstructedParticleData legs);
};

struct ResoWithLegs {
edm4hep::ReconstructedParticleData reso;
std::pair<TLorentzVector, TLorentzVector> legs_lv;
};
with some further defining in the .cc file. My question in this case is how I’d go about calling the two legs of my ResoWithLegs.

For further explanation, I am trying to find the deltaR for the two muons specifically that meet the requirement for the resonanceBuilder. In order to do this I’d need to write something similar to:

.Define(“RC_Muminus_q”,“ReconstructedParticle::sel_charge(-1.0,false)(selected_muons)”)
.Define(“RC_Muplus_q”, "ReconstructedParticle::sel_charge(1.0,false)(selected_muons) ")
.Define(“RC_Muminus_tlv”, "ReconstructedParticle::get_tlv( RC_Muminus_q ) ")
.Define(“RC_Muplus_tlv”, "ReconstructedParticle::get_tlv( RC_Muplus_q ) ")
.Define(“RCdeltaR”, " if ( RC_Muminus_tlv.size() > 0 && RC_Muplus_tlv.size() > 0) return RC_Muminus_tlv[0].DeltaR( RC_Muplus_tlv[0] ) ; else return -9999. ; ")

but with the two Lorentz vectors being those from the ResoWithLegs. Thank you for any help or insight you are able to give. This has been causing me some headache.

Cheers,
Hayden Shaddix

jaeyserm · April 3, 2023, 12:28pm

Hi Hayden,

I had to implement something similar in the past for the Higgs mass analysis: i.e. to get the 2 muons that result in the highest quality resonance: the mass should be as close as possible to the Higgs mass and the recoil as close as possible to the Z mass.

I have an implementation here:

github.com

jeyserma/FCCAnalyzer/blob/main/analysis/include/zh_mass_xsec.h#L82


      
                  p.momentum.x = p.momentum.x*(1. + scaleunc);

                  p.momentum.y = p.momentum.y*(1. + scaleunc);

                  p.momentum.z = p.momentum.z*(1. + scaleunc);

                  result.emplace_back(p);

              }

              return result;

          }

          

          

          

          // build the Z resonance based on the available leptons. Returns the best lepton pair compatible with the Z mass and recoil at 125 GeV

          // technically, it returns a ReconstructedParticleData object with index 0 the di-lepton system, index and 2 the leptons of the pair

          struct resonanceBuilder_mass_recoil {

              float m_resonance_mass;

              float m_recoil_mass;

              float chi2_recoil_frac;

              float ecm;

              bool m_use_MC_Kinematics;

              resonanceBuilder_mass_recoil(float arg_resonance_mass, float arg_recoil_mass, float arg_chi2_recoil_frac, float arg_ecm, bool arg_use_MC_Kinematics);

              Vec_rp operator()(Vec_rp legs, Vec_i recind, Vec_i mcind, Vec_rp reco, Vec_mc mc, Vec_i parents, Vec_i daugthers) ;

          };

So it takes the usually reconstructed particles as input (+ the MC information for some truth matching), and it returns a vector of 2 reconstructed particles that satisfy the best resonance criteria.

Could you use this code and adapt it to your needs? Let me know.

Best,
Jan

hshaddix · April 4, 2023, 5:26pm

Jan,

I have been trying to use it in conjunction with some code similar to yours in FCCAnalyzers:

github.com

jeyserma/FCCAnalyzer/blob/main/analysis/histmakers/ZH_mass_xsec.py


import analysis, functions, helpers
import ROOT
import argparse

parser = argparse.ArgumentParser()
parser.add_argument("--nThreads", type=int, help="number of threads", default=None)
parser.add_argument("--maxFiles", type=int, help="Max number of files (per dataset)", default=-1)
parser.add_argument("--flavor", type=str, help="Flavor (mumu or ee)", choices=["mumu", "ee"], default="mumu")
parser.add_argument("--type", type=str, help="Run type (mass or xsec)", choices=["mass", "xsec"], default="mass")
args = parser.parse_args()

functions.set_threads(args)


# define histograms
bins_p_mu = (20000, 0, 200) # 10 MeV bins
bins_m_ll = (3000, 0, 300) # 100 MeV bins
bins_p_ll = (200, 0, 200) # 1 GeV bins
bins_recoil = (200000, 0, 200) # 1 MeV bins

This file has been truncated. show original

But when I clone the repository and move down into FCCanalyses and add the branches between around 233-244, I get errors about the namespace of resonanceBuilder_mass_recoil. It seems as if it is defined properly in the file, unless I am supposed to move the files into analyzers/dataframe/FCCAnalyses or something.

Thanks,
Hayden Shaddix

jaeyserm · April 5, 2023, 7:41am

Hi Hayden,

You should not clone and run the entire repo, it’s a slightly different framework. Just copy the C++ function to your analysis code, and adapt it to make it work. Let me know if you need help with that.

Best,
Jan

hshaddix · April 5, 2023, 10:20am

Jan,

I would like some help or pointers with that if possible. I’d tried taking the structure for resonanceBuilder_mass_recoil and putting it in my analyzers/dataframe/FCCAnalyses/ReconstructedParticle.h while I put the resonanceBuilder_mass_recoil::resonanceBuilder_mass_recoil section in my ReconstructedParticle.cc along with the defining for deltaR, but it seems there are issues with the namespace. I’m unsure if that requires me to make a structure for deltaR or there are other issues.

Thanks,
Hayden Shaddix

jaeyserm · April 5, 2023, 12:11pm

Yes, the namespace is different. What you need to do is to replace the following aliases:

Vec_rp → ROOT::VecOps::RVecedm4hep::ReconstructedParticleData
Vec_mc → ROOT::VecOps::RVecedm4hep::MCParticleData
Vec_i → ROOT::VecOps::RVec

To use it in the dataframe, you can look at the code here:

github.com

jeyserma/FCCAnalyzer/blob/main/analysis/histmakers/ZH_mass_xsec.py#L233-L240


      
          df = df.Define("zbuilder_result", "FCCAnalyses::resonanceBuilder_mass_recoil(91.2, 125, 0.4, 240, false)(leps, MCRecoAssociations0, MCRecoAssociations1, ReconstructedParticles, Particle, Particle0, Particle1)")

          df = df.Define("zll", "ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData>{zbuilder_result[0]}") # the Z

          df = df.Define("zll_leps", "ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData>{zbuilder_result[1],zbuilder_result[2]}") # the leptons

          df = df.Define("zll_m", "FCCAnalyses::ReconstructedParticle::get_mass(zll)[0]")

          df = df.Define("zll_p", "FCCAnalyses::ReconstructedParticle::get_p(zll)[0]")

          df = df.Define("zll_recoil", "FCCAnalyses::ReconstructedParticle::recoilBuilder(240)(zll)")

          df = df.Define("zll_recoil_m", "FCCAnalyses::ReconstructedParticle::get_mass(zll_recoil)[0]")

          df = df.Define("zll_category", "FCCAnalyses::polarAngleCategorization(0.8, 2.34)(zll_leps)")

Hope it helps.

Best,
Jan

hshaddix · April 5, 2023, 2:07pm

Jan,

So I would basically replace those aliases in all the places they are defined in the resonanceBuilder_mass_recoil and deltaR and then get rid of the defines.h header file.?

Thanks,
Hayden Shaddix

hshaddix · April 5, 2023, 2:27pm

Jan,

Something seems to have improved. Now I am getting an error message that at least two leptons are required. I am interested in the two muons that reconstruct to 5 GeV for the process I am looking at. So I have my resonanceBuilder_mass_recoil set up as follows:

resonanceBuilder_mass_recoil(5,125,0.4,240,false)

The only thing I can think is that the mass of the particle I want to reconstruct is so low that the code assumes it is not the product of two muons, unless I have one of the other parameters set up incorrectly. I appreciate how helpful you’ve been.

Cheers,
Hayden Shaddix

jaeyserm · April 5, 2023, 2:38pm

Hi Hayden,

Great, you’re making progress! Yes, the resonance builder requires the input lepton collection to be at least 2, otherwise, it can’t make a resonance. Before calling the resonance code, you should add a filter like df = df.Filter(“leps_no >= 2”).

Best,
Jan

hshaddix · April 5, 2023, 2:48pm

Jan,

It works!!! Sorry for making you walk me through things so thoroughly. Sometimes the errors given for simple things can be very off-putting and make the troubleshooting not as obvious since I am not entirely familiar with the code structures. Thank you very much.

Cheers,
Hayden Shaddix

jaeyserm · April 5, 2023, 6:44pm

Great that it is working. Eventually, I’ll add the code to the repository so other people also can use it.

Best,
Jan

hshaddix · April 5, 2023, 7:45pm

Jan,

I actually wanted to ask if there are any examples of your resonanceBuilder_mass_recoil that selects legs based off of just the mass itself of the resonanceBuilder mass. From what I see of the code for mass_recoil it’d be a bit simpler, just removing the recoil information and defining the d and d_min to allow looping over the mass values.

Thanks,
Hayden Shaddix

hshaddix · April 5, 2023, 10:08pm

To continue, I mangled the resonanceBuilder to become something like this:

massBuilder::massBuilder(float arg_resonance_mass, float arg_recoil_mass, float arg_chi2_recoil_frac, float arg_ecm, bool arg_use_MC_Kinematics){resonance_mass = arg_resonance_mass, recoil_mass = arg_recoil_mass, chi2_recoil_frac = arg_c
hi2_recoil_frac, ecm = arg_ecm, arg_use_MC_Kinematics = arg_use_MC_Kinematics;}

ROOT::VecOps::RVecedm4hep::ReconstructedParticleData massBuilder::massBuilder::operator()(ROOT::VecOps::RVecedm4hep::ReconstructedParticleData legs, ROOT::VecOps::RVec recind, ROOT::VecOps::RVec mcind, ROOT::VecOps::RVec<ed
m4hep::ReconstructedParticleData> reco, ROOT::VecOps::RVecedm4hep::MCParticleData mc, ROOT::VecOps::RVec parents, ROOT::VecOps::RVec daugthers) {

ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData> result;
result.reserve(3);
std::vector<std::vector<int>> pairs; // for each permutation, add the indices of the muons                                                                                                                                               \
                                                                                                                                                                                                                                          
int n = legs.size();

if(n > 1) {
  ROOT::VecOps::RVec<bool> v(n);
  std::fill(v.end() - 2, v.end(), true); // helper variable for permutations                                                                                                                                                             \
                                                                                                                                                                                                                                          
  do {
    std::vector<int> pair;
    edm4hep::ReconstructedParticleData reso;
    reso.charge = 0;
    TLorentzVector reso_lv;
    for(int i = 0; i < n; ++i) {
      if(v[i]) {
        pair.push_back(i);
        reso.charge += legs[i].charge;
        TLorentzVector leg_lv;

if(use_mc_kinematics) { // MC kinematics
int track_index = legs[i].tracks_begin; // index in the Track array \

          int mc_index = ReconstructedParticle2MC::getTrack2MC_index(track_index, recind, mcind, reco);
          if (mc_index >= 0 && mc_index < mc.size()) {
            leg_lv.SetXYZM(mc.at(mc_index).momentum.x, mc.at(mc_index).momentum.y, mc.at(mc_index).momentum.z, mc.at(mc_index).mass);
          }
        }
        else { // reco kinematics                                                                                                                                                                                                        \
                                                                                                                                                                                                                                          
          leg_lv.SetXYZM(legs[i].momentum.x, legs[i].momentum.y, legs[i].momentum.z, legs[i].mass);
        }

        reso_lv += leg_lv;
      }
    }
    if(reso.charge != 0) continue; // neglect non-zero charge pairs                                                                                                                                                                      \
                                                                                                                                                                                                                                          
    reso.momentum.x = reso_lv.Px();
    reso.momentum.y = reso_lv.Py();
    reso.momentum.z = reso_lv.Pz();
    reso.mass = reso_lv.M();
    result.emplace_back(reso);
    pairs.push_back(pair);

  } while(std::next_permutation(v.begin(), v.end()));
}

else {
std::cout << “ERROR: resonanceBuilder_mass_recoil, at least two leptons required.” << std::endl;
exit(1);
}

if(result.size() > 1) {

  ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData> bestReso;

  int idx_min = -1;
  float mass_min = 9e9;
  for (int i = 0; i < result.size(); ++i) {                                                                                                                                                                                              \

    TLorentzVector tg;
    tg.SetXYZM(result.at(i).momentum.x, result.at(i).momentum.y, result.at(i).momentum.z, result.at(i).mass);

    float boost = tg.P();
    float mass = (result.at(i).mass - resonance_mass); // mass                                                                                                                                                                            

    if(mass < mass_min) {
      mass_min = mass;
      idx_min = i;
    }

}
if(idx_min > -1) {
bestReso.push_back(result.at(idx_min));
auto & l1 = legs[pairs[idx_min][0]];
auto & l2 = legs[pairs[idx_min][1]];
bestReso.emplace_back(l1);
bestReso.emplace_back(l2);
}
else {
std::cout << “ERROR: resonanceBuilder_mass_recoil, no mininum found.” << std::endl;
exit(1);
}
return bestReso;
}
else {
auto & l1 = legs[0];
auto & l2 = legs[1];
result.emplace_back(l1);
result.emplace_back(l2);
return result;
}
}

I really was thinking it would be something like taking the recoil calculations and changing the definitions of d and d_min where i replaced with mass but it seems it would be more involved than that.

Thanks again,
Hayden Shaddix

jaeyserm · April 5, 2023, 10:35pm

Hi Hayden,

In case you only want to compare the mass values and not the recoil value, you can put the recoil fraction “chi2_recoil_frac” to zero: resonanceBuilder_mass_recoil(91.2, 125, 0.0, 240, false). So there is in principle no need to re-write the function then.

Best,
Jan

hshaddix · April 6, 2023, 1:18am

Jan,

I see now. Thank you for the clarification. That’s very convenient.

Best,
Hayden