StreakAnalyzer/src/StreakAnalyzer.cc

0001 #include <calobase/PhotonClusterContainer.h>
0002 #include <caloreco/PhotonClusterBuilder.h>
0003 #include <Calobase/PhotonClusterv1.h>
0004
0005 #include <jetbase/JetContainer.h>
0006 #include <jetbase/Jet.h>
0007
0008 #include "StreakAnalyzer.h"
0009
0010 using namespace streak;
0011
0012 // -----------------------------------------------------------
0013 //Getting the Photon Node: PHOTONCLUSTER_CEMC and calcualating physics variables
0014 //-----------------------------------------------------------------------
0015 bool StreakAnalyzer::pass_cluster_selection(PHCompositeNode* topNode, std::string* why) const
0016 {
0017   auto* phCont = findNode::getClass<RawClusterContainer>(topNode, m_cfg.photon_node);
0018   if (!phCont) {
0019     if (m_cfg.verbose && why) (*why) += "PhotonClusterContainer '" + m_cfg.photon_node + "' not found.\n";
0020     return false;
0021   }
0022
0023   const auto range = phCont->getClusters();
0024   int n_total = 0, n_pass = 0;
0025
0026   for (auto it = range.first; it != range.second; ++it) {
0027     const RawCluster* rc = it->second;
0028     if (!ph) continue;
0029
0030     PhotonClusterv1* ph = dynamic_cast<PhotonClusterv1*>(rc);
0031     if (!ph) continue;
0032
0033     ++n_total;
0034
0035     //Getting the photon cluster eta, E and Et from the PhotonClusterBuilder
0036     const float eta = ph->get_shower_shape_parameter("cluster_eta");
0037     const float E   = ph->get_shower_shape_parameter("e77");
0038     const float Et  = (std::isfinite(eta) && std::cosh(eta) != 0.0) ? (E / std::cosh(eta)) : 0.f;
0039
0040     //Getting weta_Center-of-gravity in x direction
0041     float weta = ph->get_shower_shape_parameter("weta_cogx");
0042     if (!std::isfinite(weta) || weta <= 0.f){
0043
0044     // backup if 'weta_cogx' is not found, so we can take 'weta_cog' instead
0045       weta = ph->get_shower_shape_parameter("weta_cog");
0046     }
0047
0048       //cuts for Et and weta to identify streak clusters
0049     const bool pass = (Et >= m_cfg.et_min) && (weta >= m_cfg.weta_min);
0050     if (pass) ++n_pass;
0051
0052     if (m_cfg.verbose && why) {
0053       (*why) += "Cluster: Et=" + std::to_string(Et)
0054              + " GeV, weta=" + std::to_string(weta)
0055              + (pass ? " [PASS]\n" : " [FAIL]\n");
0056     }
0057   }
0058
0059   //If the event has been rejected due to no clusters passing, provide info
0060   if (m_cfg.verbose && why) {
0061     (*why) += "Clusters passing = " + std::to_string(n_pass) + " / " + std::to_string(n_total) + "\n";
0062   }
0063   return (n_pass > 0);
0064 }
0065
0066 // -----------------------------------------------------------
0067 //Dijet veto: looking for back-to-back jets in specified jet containers
0068 // -----------------------------------------------------------
0069 bool StreakAnalyzer::pass_dijet_veto(PHCompositeNode* topNode, std::string* why) const
0070 {
0071   if (!m_cfg.apply_dijet_veto) return true;
0072
0073   // find first available jet container from preference list
0074   JetContainer* jets = nullptr;
0075   std::string picked;
0076
0077   for (const auto& name : m_cfg.jet_nodes) {
0078
0079     if (auto* found = findNode::getClass<JetContainer>(topNode, name)) {
0080       jets = const_cast<JetContainer*>(found);
0081       picked = name;
0082       break;
0083     }
0084   }
0085
0086   // If no jets found, cannot veto
0087   if (!jets) {
0088     if (m_cfg.verbose && why) (*why) += "Dijet veto: no JetContainer found \n";
0089     return true;
0090   }
0091
0092   // Collect jet phi of jets above threshold
0093   std::vector<float> phi;
0094   phi.reserve(jets->size());
0095
0096
0097   for (Jet::IterJetTCA it = jets->begin(); it != jets->end(); ++it) {
0098     Jet* j = *it;
0099     if (!j) continue;
0100     if (j->get_pt() >= m_cfg.jet_pt_min) phi.push_back(j->get_phi());
0101   }
0102
0103   // If not enough jets above threshold, cannot veto
0104   if (m_cfg.verbose && why) {
0105     (*why) += "Dijet veto: picked '" + picked + "', jets above "
0106            + std::to_string(m_cfg.jet_pt_min) + " GeV = "
0107            + std::to_string(phi.size()) + "\n";
0108   }
0109
0110   // veto of back-to-back jets
0111   const int n = static_cast<int>(phi.size());
0112   for (int i = 0; i < n; ++i) {
0113     for (int k = i + 1; k < n; ++k) {
0114       if (dphi_abs(phi[i], phi[k]) >= m_cfg.back_to_back_min_dphi) {
0115         if (m_cfg.verbose && why)
0116           (*why) += "Dijet veto: back-to-back jets found -> FAIL.\n";
0117         return false;
0118       }
0119     }
0120   }
0121   return true; // when there is  no veto
0122 }
0123
0124 // -----------------------------------------------------------
0125 //Main function to determine if event is streak event
0126 // -----------------------------------------------------------
0127 bool StreakAnalyzer::isStreakEvent(PHCompositeNode* topNode, std::string* why) const
0128 {
0129
0130
0131   if (m_cfg.verbose && why) {
0132   (*why) += "Streak config: photon_node='" + m_cfg.photon_node
0133          + "', et_min="   + std::to_string(m_cfg.et_min)
0134          + " GeV, weta_min=" + std::to_string(m_cfg.weta_min)
0135          + ", dijet_veto=" + std::string(m_cfg.apply_dijet_veto ? "on" : "off")
0136          + ", jet_pt_min=" + std::to_string(m_cfg.jet_pt_min)
0137          + ", dphi_min="   + std::to_string(m_cfg.back_to_back_min_dphi) + "\n";
0138 }
0139
0140
0141   if (!topNode) return false;
0142   if (!pass_cluster_selection(topNode, why)) return false;
0143   if (!pass_dijet_veto(topNode, why))        return false;
0144   return true;
0145 }