count constituents for jet charge calculation

mverwe · Apr 15, 2021 · 2fa226a · 2fa226a
1 parent f0c88e2
commit 2fa226a
Show file tree

Hide file tree

Showing 2 changed files with 52 additions and 4 deletions.
diff --git a/include/JetCharge.hh b/include/JetCharge.hh
@@ -18,20 +18,30 @@ public:
   virtual double result(const fastjet::PseudoJet &jet) const {
     // check the jet is appropriate for computation
     if (!jet.has_constituents()) {
-      Printf("Angularities can only be applied on jets for which the constituents are known.");
+      Printf("Jet charge calculation can only be applied on jets for which the constituents are known.");
       return -999.;
     }
     vector<fastjet::PseudoJet> constits = jet.constituents();
     double sumcharge = 0.;
     double sumpt = 0.;
+    int nconst = 0;
     for(fastjet::PseudoJet p : constits) {
       if(p.perp()<_ptmin) continue;
-      const int & ch = p.user_info<PU14>().charge(); //three_charge()
+      const double & ch = p.user_info<PU14>().charge(); //three_charge()
       //std::cout << "charge: " << ch << " three_charge: " << p.user_info<PU14>().three_charge() << "  pdg: " << p.user_info<PU14>().pdg_id() << std::endl;
       sumcharge += ch*std::pow(p.pt(),_kappa);
       sumpt += std::pow(p.pt(),_kappa);
+      ++nconst;
     }
-    double charge = sumcharge/sumpt;
+    double charge = -999.;
+    if(sumpt>0.) charge = sumcharge/sumpt;
+    // if(abs(charge)<0.01) {
+    //   std::cout << "ptmin: " << _ptmin << "  charge: " << charge << "  jetpt: " << jet.perp() << "  sumcharge: " << sumcharge << " sumpt: " << sumpt << "  nconst: " << nconst<< std::endl;
+    //   for(fastjet::PseudoJet p : constits) {
+    //     if(p.perp()<_ptmin) continue;
+    //     std::cout << " const pt: " << p.perp() << " const charge: " << p.user_info<PU14>().charge() << "  pdg: " << p.user_info<PU14>().pdg_id() << std::endl;
+    //   }
+    // }
     return charge;
   }
 

diff --git a/runJetChargeQCDAware.cc b/runJetChargeQCDAware.cc
@@ -55,7 +55,7 @@ int main (int argc, char ** argv) {
   treeWriter trwSig("jetTreeSig");
 
   bool runFull = true;
-  bool runCharged = false;//true;
+  bool runCharged = true;
 
   //Jet definition
   double R                   = 0.4;
@@ -213,6 +213,25 @@ int main (int argc, char ** argv) {
       jetCollectionSig.addVector("sigJetChargePtmin1", jetChargePtmin1Sig);
       jetCollectionSig.addVector("sigJetChargePtmin2", jetChargePtmin2Sig);
 
+      //count constituents entering jet charge calculation
+      vector<int> jetNConstSig;       jetNConstSig.reserve(jetCollectionSig.getJet().size());
+      vector<int> jetNConstPtmin1Sig; jetNConstPtmin1Sig.reserve(jetCollectionSig.getJet().size());
+      vector<int> jetNConstPtmin2Sig; jetNConstPtmin2Sig.reserve(jetCollectionSig.getJet().size());
+      for(PseudoJet jet : jetCollectionSig.getJet()) {
+        int nconst[3] = {0,0,0};
+        for(fastjet::PseudoJet p : jet.constituents()) {
+          ++nconst[0];
+          if(p.perp()>1.) ++nconst[1];
+          if(p.perp()>2.) ++nconst[2];
+        }
+        jetNConstSig.push_back(nconst[0]);
+        jetNConstPtmin1Sig.push_back(nconst[1]);
+        jetNConstPtmin2Sig.push_back(nconst[2]);
+      }
+      jetCollectionSig.addVector("sigJetNConst", jetNConstSig);
+      jetCollectionSig.addVector("sigJetNConstPtmin1", jetNConstPtmin1Sig);
+      jetCollectionSig.addVector("sigJetNConstPtmin2", jetNConstPtmin2Sig);
+
       trwSig.addCollection("sigJet",        jetCollectionSig);
     }
 
@@ -259,6 +278,25 @@ int main (int argc, char ** argv) {
       jetCollectionSigCh.addVector("sigJetChChargePtmin1", jetChargePtmin1SigCh);
       jetCollectionSigCh.addVector("sigJetChChargePtmin2", jetChargePtmin2SigCh);
 
+      //count constituents entering jet charge calculation
+      vector<int> jetNConstSigCh;       jetNConstSigCh.reserve(jetCollectionSigCh.getJet().size());
+      vector<int> jetNConstPtmin1SigCh; jetNConstPtmin1SigCh.reserve(jetCollectionSigCh.getJet().size());
+      vector<int> jetNConstPtmin2SigCh; jetNConstPtmin2SigCh.reserve(jetCollectionSigCh.getJet().size());
+
+      for(PseudoJet jet : jetCollectionSigCh.getJet()) {
+        int nconst[3] = {0,0,0};
+        for(fastjet::PseudoJet p : jet.constituents()) {
+          ++nconst[0];
+          if(p.perp()>1.) ++nconst[1];
+          if(p.perp()>2.) ++nconst[2];
+        }
+        jetNConstSigCh.push_back(nconst[0]);
+        jetNConstPtmin1SigCh.push_back(nconst[1]);
+        jetNConstPtmin2SigCh.push_back(nconst[2]);
+      }
+      jetCollectionSigCh.addVector("sigJetChNConst", jetNConstSigCh);
+      jetCollectionSigCh.addVector("sigJetChNConstPtmin1", jetNConstPtmin1SigCh);
+      jetCollectionSigCh.addVector("sigJetChNConstPtmin2", jetNConstPtmin2SigCh);
 
       //---------------------------------------------------------------------------
       //   write tree