Adding ghost matching jet flavour defn

analyzers/dataframe/JetTaggingUtils.cc

@@ -61,6 +61,210 @@ JetTaggingUtils::get_flavour(ROOT::VecOps::RVec<fastjet::PseudoJet> in,
   return result;
 }
 
+
+get_ghostFlavour::get_ghostFlavour(int algo, float arg_radius, int arg_exclusive, float arg_cut, int arg_sorted, int arg_recombination, float arg_add1, float arg_add2)
+{m_algo = algo; m_radius = arg_radius; m_exclusive = arg_exclusive; m_cut = arg_cut; m_sorted = arg_sorted; m_recombination = arg_recombination; m_add1 = arg_add1; m_add2 = arg_add2;}
+
+ghostFlavour JetTaggingUtils::get_ghostFlavour::operator() (ROOT::VecOps::RVec<edm4hep::MCParticleData> Particle, ROOT::VecOps::RVec<int> ind, std::vector<fastjet::PseudoJet> pseudoJets, int partonFlag) {
+
+
+
+  ghostFlavour result;
+
+  unsigned int index = pseudoJets.size();
+  std::vector<float> pdg(pseudoJets.size(),0);
+  std::vector<float> ghostStatus(pseudoJets.size(),0);
+  std::vector<int> MCindex(pseudoJets.size(),-1);
+
+
+  std::vector<int> partonStatus = {20, 30};
+
+  // In loop below ghosts are selected from MCParticle collection
+  for (size_t i = 0; i < Particle.size(); ++i) {
+    bool isGhost = false;
+
+    if(!partonFlag){
+      // Ghost partons as any partons that do not have partons as daughters
+      if (std::abs(Particle[i].PDG)<=5||Particle[i].PDG==21) {
+        isGhost = true;
+        auto daughters = MCParticle::get_list_of_particles_from_decay(i, Particle, ind);
+        for(auto& daughter_index : daughters){
+          if (std::abs(Particle[daughter_index].PDG)<=5||Particle[daughter_index].PDG==21) {
+              isGhost = false;
+              break;
+            }
+        }
+        if (isGhost) ghostStatus.push_back(1);
+      }
+    }
+    else{
+      // Ghost partons are selected via chosen Pythia8 status codes
+      if ((Particle[i].generatorStatus<partonStatus[1]) && (Particle[i].generatorStatus>partonStatus[0])) {
+        isGhost = true;
+        ghostStatus.push_back(1);
+      }
+    }
+
+
+    // Ghost hadrons are selected as b/c hadrons that do not have b/c hadrons as daughters
+    if ((std::abs(int((Particle[i].PDG/100))%10)==5)||(std::abs(int((Particle[i].PDG/1000))%10)==5)){
+      isGhost = true;
+      auto daughters = MCParticle::get_list_of_particles_from_decay(i, Particle, ind);
+      for(auto& daughter_index : daughters){
+          if((std::abs(int((Particle[daughter_index].PDG/100))%10)==5)||(std::abs(int((Particle[daughter_index].PDG/1000))%10)==5)){
+            isGhost = false;
+            break;
+          }
+      }
+      if (isGhost) ghostStatus.push_back(2);
+
+    }
+    else if ((std::abs(int((Particle[i].PDG/100))%10)==4)||(std::abs(int((Particle[i].PDG/1000))%10)==4)){
+      isGhost = true;
+      auto daughters = MCParticle::get_list_of_particles_from_decay(i, Particle, ind);
+      for(auto& daughter_index : daughters){
+          if((std::abs(int((Particle[daughter_index].PDG/100))%10)==4)||(std::abs(int((Particle[daughter_index].PDG/1000))%10)==4)){
+            isGhost = false;
+            break;
+          }
+      }
+      if (isGhost) ghostStatus.push_back(2);
+    }
+
+    // Ghosts 4-mom is scaled by 10^-18  
+    if (isGhost){
+      pdg.push_back(Particle[i].PDG);
+      MCindex.push_back(i);
+      // the double conversion here is verbose but for precision if I recall...
+      double px = Particle[i].momentum.x;//*pow(10, -1);
+      double py = Particle[i].momentum.y;
+      double pz = Particle[i].momentum.z;
+      double m = Particle[i].mass;
+      double  E = sqrt(px*px + py*py + pz*pz + m*m);
+      pseudoJets.emplace_back(px*pow(10, -18), py*pow(10, -18), pz*pow(10, -18), E*pow(10, -18));
+      pseudoJets.back().set_user_index(index);
+      ++index;
+    }
+  }
+
+  result.ghostStatus = ghostStatus;
+  result.MCindex = MCindex;
+
+
+  // Jet clustering algorithm is run according to user choice m_algo
+  JetClusteringUtils::FCCAnalysesJet FCCAnalysesGhostJets;
+
+  if (m_algo == 0) FCCAnalysesGhostJets = JetClustering::clustering_kt(m_radius, m_exclusive, m_cut, m_sorted, m_recombination)(pseudoJets);
+  else if (m_algo == 1) FCCAnalysesGhostJets = JetClustering::clustering_antikt(m_radius, m_exclusive, m_cut, m_sorted, m_recombination)(pseudoJets);
+  else if (m_algo == 2) FCCAnalysesGhostJets = JetClustering::clustering_cambridge(m_radius, m_exclusive, m_cut, m_sorted, m_recombination)(pseudoJets);
+  else if (m_algo == 3) FCCAnalysesGhostJets = JetClustering::clustering_ee_kt(m_exclusive, m_cut, m_sorted, m_recombination)(pseudoJets);
+  else if (m_algo == 4) FCCAnalysesGhostJets = JetClustering::clustering_ee_genkt(m_radius, m_exclusive, m_cut, m_sorted, m_recombination, m_add1)(pseudoJets);
+  else if (m_algo == 5) FCCAnalysesGhostJets = JetClustering::clustering_genkt(m_radius, m_exclusive, m_cut, m_sorted, m_recombination, m_add1)(pseudoJets);
+  else if (m_algo == 6) FCCAnalysesGhostJets = JetClustering::clustering_valencia(m_radius, m_exclusive, m_cut, m_sorted, m_recombination, m_add1, m_add2)(pseudoJets);
+  else if (m_algo == 7) FCCAnalysesGhostJets = JetClustering::clustering_jade(m_radius, m_exclusive, m_cut, m_sorted, m_recombination)(pseudoJets);
+  else return result;
+
+
+  result.jets = FCCAnalysesGhostJets;
+
+  // Jet constituents and pseudojets are read from resulting jet struct
+  auto ghostJets_ee_kt = JetClusteringUtils::get_pseudoJets(FCCAnalysesGhostJets);
+
+  auto jetconstituents = JetClusteringUtils::get_constituents(FCCAnalysesGhostJets);
+
+
+
+
+
+
+  // Flav vector is defined before jets are checked for clustered ghosts
+  std::vector<std::vector<float>> flav_vector;
+
+
+  std::vector<float> partonFlavs;
+  std::vector<float> hadronFlavs;
+  for (auto& consti_index : jetconstituents) {
+
+  float partonFlav = 0;
+  float partonMom2 = 0;
+  float partonMom2_b = 0;
+  float partonMom2_c = 0;
+  float hadronFlav = 0;
+  float hadronMom2_b = 0;
+  float hadronMom2_c = 0;
+
+    for (auto& i : consti_index) {
+
+      //Parton-flav loop 
+      if (ghostStatus[i]==1) {
+          if (std::abs(pdg[i])==5){
+              if (pseudoJets[i].modp2()>partonMom2_b){
+                  partonFlav = pdg[i];
+                  partonMom2_b = pseudoJets[i].modp2();
+              }
+          }
+          else if ((std::abs(pdg[i])==4) && (std::abs(partonFlav)<5)) {
+              if (pseudoJets[i].modp2()>partonMom2_c){
+                  partonFlav = pdg[i];
+                  partonMom2_c = pseudoJets[i].modp2();
+              }
+          }
+          else if (((std::abs(pdg[i])==3) || (std::abs(pdg[i])==2) || (std::abs(pdg[i])==1) || (std::abs(pdg[i])==21)) && ((std::abs(partonFlav)<4) || (std::abs(partonFlav)==21))) {
+              if (pseudoJets[i].modp2()>partonMom2){
+                  partonFlav = pdg[i];
+                  partonMom2 = pseudoJets[i].modp2();
+              }
+          }
+      }
+
+      // Hadron-flav loop    
+      if (ghostStatus[i]==2) {
+          if ((std::abs(int((pdg[i]/100))%10)==5)||(std::abs(int((pdg[i]/1000))%10)==5)){
+              if (pseudoJets[i].modp2()>hadronMom2_b){
+                  hadronFlav = ((pdg[i]<0)-(pdg[i]>0))*5;
+                  hadronMom2_b = pseudoJets[i].modp2();
+              }
+          }
+          else if (((std::abs(int((pdg[i]/100))%10)==4)||(std::abs(int((pdg[i]/1000))%10)==4)) && (std::abs(hadronFlav)<5)) {
+              if (pseudoJets[i].modp2()>hadronMom2_c){
+                  hadronFlav = ((pdg[i]>0)-(pdg[i]<0))*4;
+                  hadronMom2_c = pseudoJets[i].modp2();
+              }
+          }
+      }
+    }
+    partonFlavs.push_back(partonFlav);
+    hadronFlavs.push_back(hadronFlav);
+  }
+
+
+  flav_vector.push_back(partonFlavs);
+  flav_vector.push_back(hadronFlavs);
+
+  result.flavour = flav_vector;
+  return result;
+
+}
+
+std::vector<std::vector<float>> JetTaggingUtils::get_flavour(ghostFlavour ghostStruct){
+  return ghostStruct.flavour;
+}
+
+JetClusteringUtils::FCCAnalysesJet JetTaggingUtils::get_jets(ghostFlavour ghostStruct){
+  return ghostStruct.jets;
+}
+
+std::vector<float> JetTaggingUtils::get_ghostStatus(ghostFlavour ghostStruct){
+  return ghostStruct.ghostStatus;
+}
+
+std::vector<int> JetTaggingUtils::get_MCindex(ghostFlavour ghostStruct){
+  return ghostStruct.MCindex;
+}
+
+
+
+
 ROOT::VecOps::RVec<int>
 JetTaggingUtils::get_btag(ROOT::VecOps::RVec<int> in,
                           float efficiency, float mistag_c,

analyzers/dataframe/JetTaggingUtils.h

@@ -5,8 +5,12 @@
 #include "Math/Vector4D.h"
 #include "ROOT/RVec.hxx"
 #include "edm4hep/MCParticleData.h"
+#include "edm4hep/ReconstructedParticleData.h"
 #include "fastjet/JetDefinition.hh"
 #include "TRandom3.h"
+#include "JetClustering.h"
+#include "JetClusteringUtils.h"
+#include "MCParticle.h"
 
 /** Jet tagging utilities interface.
 This represents a set functions and utilities to perfom jet tagging from a list of jets.
@@ -20,6 +24,44 @@ namespace JetTaggingUtils{
 
   //Get flavour association of jet 
   ROOT::VecOps::RVec<int> get_flavour(ROOT::VecOps::RVec<fastjet::PseudoJet> in, ROOT::VecOps::RVec<edm4hep::MCParticleData> MCin);
+
+  struct ghostFlavour {
+    std::vector<std::vector<float>> flavour;
+    JetClusteringUtils::FCCAnalysesJet jets;
+    std::vector<float> ghostStatus;
+    std::vector<int> MCindex;
+  };
+
+
+  //Get ghost flavour (MC flavour) of jet described here: ..
+  //a struct is returned containing
+  // - vector of flavours: <<partonFlav Jet 1, partonFlav Jet 2, ...>, <hadronFlav Jet 1, hadron Flav Jet 2, ...>>
+  // - FCCAnalysesJet struct containing resulting jets (on which functions defined in JetClusteringUtils.cc can be called)
+  // - vector of jet constituent indices <<constis jet 1>, <constis jet 2>, ...>
+  // - ghostStatus (0 = not ghost, 1 = ghost parton, 2 = ghost hadron)
+  // - MCindex vector of indices of ghosts to MC collection (and -1 if not a ghost) <part 1, ..., part n>
+  struct get_ghostFlavour {
+    get_ghostFlavour(int algo, float arg_radius, int arg_exclusive, float arg_cut, int arg_sorted, int arg_recombination, float arg_add1=0, float arg_add2=0);
+
+    int m_algo = 0; ///< flag to select jet clustering algorithm defined in JetClustering.cc (0 = kt, 1 = antikt, 2 = cambridge, 3 = eekt, 4 = ee genkt, 5 = genkt, 6 = valencia, 7 = jade) 
+    float m_radius = 0.5; ///< jet cone radius (note that this variable should be passed even when not used e.g. for eekt)
+    int   m_exclusive = 0; ///< flag for exclusive jet clustering. Possible choices are 0=inclusive clustering, 1=exclusive clustering that would be obtained when running the algorithm with the given dcut, 2=exclusive clustering when the event is clustered (in the exclusive sense) to exactly njets, 3=exclusive clustering when the event is clustered (in the exclusive sense) up to exactly njets, 4=exclusive jets obtained at the given ycut 
+    float m_cut = 5.; ///< pT cut for m_exclusive=0, dcut for m_exclusive=1, N jets for m_exlusive=2, N jets for m_exclusive=3, ycut for m_exclusive=4
+    int m_sorted = 0; ///< pT ordering=0, E ordering=1
+    int m_recombination = 0; ///< E_scheme=0, pt_scheme=1, pt2_scheme=2, Et_scheme=3, Et2_scheme=4, BIpt_scheme=5, BIpt2_scheme=6, E0_scheme=10, p_scheme=11
+    float m_add1 = 0.; /// first additional parameter (to be used in selected clustering scheme)
+    float m_add2 = 0.; /// second additional parameter
+    ghostFlavour operator() (ROOT::VecOps::RVec<edm4hep::MCParticleData> Particle, ROOT::VecOps::RVec<int> ind, std::vector<fastjet::PseudoJet> pseudoJets, int partonFlag);
+  };
+
+  std::vector<std::vector<float>> get_flavour(ghostFlavour ghostStruct);
+
+  JetClusteringUtils::FCCAnalysesJet get_jets(ghostFlavour ghostStruct);
+
+  std::vector<float> get_ghostStatus(ghostFlavour ghostStruct);
+
+  std::vector<int> get_MCindex(ghostFlavour ghostStruct);
+
   //Get b-tags with an efficiency applied
   ROOT::VecOps::RVec<int> get_btag(ROOT::VecOps::RVec<int> in, float efficiency, float mistag_c=0., float mistag_l=0., float mistag_g=0.);
   //Get c-tags with an efficiency applied