StreakEventsIdentifier/src/StreakEventsIdentifier.h

0001 #ifndef STREAK_EVENTS_IDENTIFIER_H
0002 #define STREAK_EVENTS_IDENTIFIER_H
0003
0004
0005 // sPHENIX core
0006 #include <fun4all/SubsysReco.h>
0007 #include <phool/PHCompositeNode.h>
0008
0009 // Cluster and tower types
0010 #include <calobase/RawCluster.h>
0011 #include <calobase/RawClusterContainer.h>
0012 #include <calobase/RawTowerGeomContainer.h>
0013 #include <calobase/TowerInfoContainer.h>
0014
0015 // Jet types
0016 #include <jetbase/JetContainer.h>
0017
0018 // ROOT
0019 #include <TH1F.h>
0020 #include <TH2F.h>
0021 #include <TFile.h>
0022
0023 #include <string>
0024 #include <map>
0025 #include <vector>
0026 #include <utility>
0027 #include <cmath>
0028 #include <memory>
0029
0030 class StreakEventsIdentifier : public SubsysReco
0031 {
0032 public:
0033   explicit StreakEventsIdentifier(const std::string& name = "StreakID",
0034                                           const std::string& outprefix = "streak_fromDST");
0035   virtual ~StreakEventsIdentifier() override;
0036
0037   // SubsysReco
0038   int Init(PHCompositeNode* topNode) override;
0039   int InitRun(PHCompositeNode* topNode) override;
0040   int process_event(PHCompositeNode* topNode) override;
0041   int End(PHCompositeNode* topNode) override;
0042
0043   // Config
0044   void SetOutputPrefix(const std::string& p) { m_outprefix = p; }
0045   void SetPngDir(const std::string& d) { m_png_dir = d; }
0046   void SetVertexCut(float cut) { m_vertex_cut = cut; }
0047   void SetEnergyThreshold(float etMin) { m_et_min = etMin; }
0048   void SetWetaCut(float wetaMin) { m_weta_min = wetaMin; }
0049   void SetTimingCuts(float absTimeNs, float rmsMinNs) { m_abs_time_cut_ns = absTimeNs; m_time_spread_min = rmsMinNs; }
0050   void SetTimeWeightEthresh(float e) { m_time_weight_Ethresh = e; }
0051
0052
0053   void SetJetContainer(const std::string& name) { m_jet_container_hint = name; }
0054
0055     // to control timing requirement
0056   void SetRequireValidTiming(bool require) { m_require_valid_timing = require; }
0057
0058
0059   //Adding MBD timing cuts
0060   void SetExcludeStreaks(bool exclude) { m_exclude_streaks = exclude; }
0061   void SetMbdTimingCuts(float minNs, float maxNs) {
0062     m_mbd_time_min = minNs;
0063     m_mbd_time_max = maxNs;
0064   }
0065
0066   //trigger bits to use for event selection
0067   void set_using_trigger_bits(const std::vector<int>& trigger_bits) {
0068     m_using_trigger_bits = trigger_bits;
0069   }
0070
0071   // Alternative: add single trigger bit
0072   void add_trigger_bit(int bit) {
0073     m_using_trigger_bits.push_back(bit);
0074   }
0075
0076 private:
0077
0078 std::pair<float, float> extractCaloTiming(TowerInfoContainer* towers,
0079                                             float energy_threshold = 0.1) const;
0080
0081   // Total calorimeter energies
0082   float m_totalEMCal_energy = 0.f;
0083   float m_totalIHCal_energy = 0.f;
0084   float m_totalOHCal_energy = 0.f;
0085
0086   static inline float safe_div(float n, float d) { return (std::fabs(d) > 1e-9f) ? n/d : 0.f; }
0087   static inline float wrap_dphi(float a, float b) {
0088     float d = a - b; while (d > M_PI) d -= 2.f*M_PI; while (d <= -M_PI) d += 2.f*M_PI; return std::fabs(d);
0089   }
0090
0091
0092 static void shift_tower_index(int& ieta, int& iphi, int neta, int nphi)
0093 {
0094   if (ieta < 0) { ieta += neta; }
0095   if (ieta >= neta) { ieta -= neta; }
0096
0097   if (iphi < 0) { iphi += nphi; }
0098   if (iphi >= nphi) { iphi -= nphi; }
0099 }
0100
0101 // Trigger selection
0102   std::vector<int> m_using_trigger_bits;
0103
0104   // Trigger data storage
0105   bool m_scaledtrigger[64] = {false};
0106   bool m_livetrigger[64] = {false};
0107   int m_nscaledtrigger[64] = {0};
0108   int m_nlivetrigger[64] = {0};
0109   float m_trigger_prescale[64] = {-1.0};
0110
0111   // Scaler tracking
0112   bool m_initilized = false;
0113   long long m_initscaler[64][3] = {{0}};
0114   long long m_currentscaler[64][3] = {{0}};
0115   long long m_currentscaler_raw[64] = {0};
0116   long long m_currentscaler_live[64] = {0};
0117   long long m_currentscaler_scaled[64] = {0};
0118
0119   // Event tracking
0120   int m_eventnumber = 0;
0121
0122   // Streaks + MBD timing exclusion
0123   bool m_exclude_streaks = false;
0124   float m_mbd_time_min = -10.0f;  // ns
0125   float m_mbd_time_max = 10.0f;   // ns
0126   int m_excluded_count = 0;
0127
0128   // Node getters
0129   RawClusterContainer* getCemcClusterContainer(PHCompositeNode* topNode) const;
0130   TowerInfoContainer*  getCemcCalibTowers(PHCompositeNode* topNode) const;
0131   TowerInfoContainer*  getCemcRawTowers(PHCompositeNode* topNode) const;
0132   RawTowerGeomContainer* getCemcGeom(PHCompositeNode* topNode) const;
0133   JetContainer*        getJetContainer(PHCompositeNode* topNode) const;
0134
0135
0136   // IHCal
0137   TowerInfoContainer*    getIhcalCalibTowers(PHCompositeNode* topNode) const;
0138   TowerInfoContainer*    getIhcalRawTowers  (PHCompositeNode* topNode) const;
0139   RawTowerGeomContainer* getIhcalGeom       (PHCompositeNode* topNode) const;
0140
0141   // OHCal
0142   TowerInfoContainer*    getOhcalCalibTowers(PHCompositeNode* topNode) const;
0143   TowerInfoContainer*    getOhcalRawTowers  (PHCompositeNode* topNode) const;
0144   RawTowerGeomContainer* getOhcalGeom       (PHCompositeNode* topNode) const;
0145
0146
0147   float getVertexZ(PHCompositeNode* topNode) const;
0148
0149   // Per-cluster calculations (7x7 around centroid)
0150   struct ClusterWindows {
0151   float E77[7][7]{};   // calibrated energy (thresholded)
0152   float T77[7][7]{};   // time (ns)
0153   int   Own77[7][7]{}; // ownership mask
0154
0155   inline float E(int i, int j) const { return E77[i][j]; }
0156   inline float T(int i, int j) const { return T77[i][j]; }
0157   inline int   Own(int i, int j) const { return Own77[i][j]; }
0158   };
0159
0160
0161   bool fillClusterWindows(const RawCluster* cl,
0162                           TowerInfoContainer* emcCalib,
0163                           TowerInfoContainer* emcRaw,
0164                           ClusterWindows& win,
0165                           int& maxieta, int& maxiphi,
0166                           float& cog_eta, float& cog_phi) const;
0167
0168   void computeWidths(const ClusterWindows& w, float cog_eta, float cog_phi,
0169                    float& weta, float& weta_cogx,
0170                    float& wphi, float& wphi_cogx) const;
0171
0172   void computeTimeMoments(const ClusterWindows& win,
0173                           float& t_avg, float& t_rms) const;
0174
0175   // Book / save histos
0176   void bookHistos();
0177   void savePNGs() const;
0178
0179 private:
0180   std::string m_outprefix;
0181   std::string m_png_dir;
0182
0183     //Timing requirement
0184   bool m_require_valid_timing = true;  //  require timing
0185   bool has_valid_timing = true;
0186
0187   // thresholds
0188   float m_vertex_cut = 200.0f; // z-verte cut in cm
0189   float m_et_min = 10.f;       // GeV
0190   float m_weta_min = 0.5f;      // Weta cut
0191   float m_abs_time_cut_ns = 10.0f; //10
0192   float m_time_spread_min = 5.0f; //5
0193   float m_min_tower_E_for_shapes = 0.07f; // GeV
0194   float m_time_weight_Ethresh = 0.50f;    // GeV
0195
0196   std::string m_jet_container_hint;
0197
0198   // bookkeeping
0199   int   m_runnumber = 0;
0200   long  m_evt_processed = 0;
0201   long  m_streak_count = 0;
0202   std::vector<std::pair<int,int>> m_streakEvents; // (run, evt)
0203   std::map<int,int> m_run_total, m_run_streak;
0204
0205   std::unique_ptr<TFile> m_out;
0206
0207   // Histograms
0208   TH2F *h_eta_phi_all = nullptr, *h_eta_phi_streak = nullptr;
0209   TH2F *h_phi_Et_all  = nullptr, *h_phi_Et_streak  = nullptr;
0210   TH1F *h_weta_all    = nullptr, *h_weta_streak    = nullptr;
0211   TH1F *h_weta_all_x    = nullptr, *h_weta_streak_x    = nullptr;
0212   TH1F *h_cluster_Et_all = nullptr, *h_cluster_Et_streak = nullptr;
0213   TH1F *h_asymmetry = nullptr, *h_jet_dphi = nullptr, *h_jet_asym = nullptr;
0214   TH2F *h_dphi_vs_pt_ratio = nullptr, *h_jet_pt_vs_phi = nullptr;
0215   TH1F *h_cluster_time_all = nullptr, *h_cluster_time_streak = nullptr;
0216   TH1F *h_time_spread_all = nullptr, *h_time_spread_streak = nullptr;
0217
0218   TH1F* h_wphi_all{nullptr};
0219   TH1F* h_wphi_streak{nullptr};
0220
0221
0222 };
0223
0224 #endif
0225 // -----------------------------------------------------------