src/TowerRho/DetermineTowerRho.h

0001 #ifndef RHOBASE_DETERMINETOWERRHO_H
0002 #define RHOBASE_DETERMINETOWERRHO_H
0003
0004 //===========================================================
0005 /// \file DetermineTowerRho.h
0006 /// \brief UE background rho calculator
0007 /// \author Tanner Mengel
0008 //===========================================================
0009
0010 #include "TowerRho.h"
0011 #include "TowerRhov1.h"
0012
0013 // fun4all includes
0014 #include <fun4all/SubsysReco.h>
0015
0016 // system includes
0017 #include <string>
0018 #include <vector>
0019 #include <ostream>
0020
0021 #include <jetbase/Jet.h>
0022 #include <jetbase/JetAlgo.h>
0023 #include <jetbase/JetInput.h>
0024
0025 #include <fastjet/JetDefinition.hh>
0026 #include <fastjet/PseudoJet.hh>
0027
0028
0029 // forward declarations
0030 class PHCompositeNode;
0031
0032
0033 /// \class DetermineTowerRho
0034 ///
0035 /// \brief UE background calculator
0036 ///
0037 /// This module estimates rho for the area and multiplicty methods using kt jets
0038 ///
0039
0040
0041 class DetermineTowerRho : public SubsysReco
0042 {
0043   public:
0044     DetermineTowerRho(const std::string &name = "DetermineTowerRho");
0045     ~DetermineTowerRho() override;
0046
0047     // standard Fun4All methods
0048     int InitRun(PHCompositeNode *topNode) override;
0049     int process_event(PHCompositeNode *topNode) override;
0050
0051     // add rho method (Area or Multiplicity)
0052     void add_method(TowerRho::Method rho_method, std::string output= ""){
0053
0054       // get method name
0055       std::string method_name = TowerRhov1::get_method_string(rho_method);
0056
0057       // check if method already exists
0058       for (auto &method : _rho_methods)
0059       {
0060         if (method == rho_method)
0061         {
0062           std::cout << "WARNING: rho method " << method_name << " already exists" << std::endl;
0063           return ;
0064         }
0065       }
0066
0067       _rho_methods.push_back(rho_method);
0068
0069       // if no output name is specified, use default
0070       if(output == "")
0071       {
0072         output = "TowerRho_" + method_name;
0073       }
0074
0075       // add output name
0076       _outputs.push_back(output);
0077
0078   }
0079
0080     // inputs for background jets
0081     void add_tower_input(JetInput* input){ _inputs.push_back(input); }
0082
0083     void set_algo(Jet::ALGO algo) { m_bkgd_jet_algo = algo; }
0084     Jet::ALGO get_algo() { return m_bkgd_jet_algo; }
0085
0086     // set the jet algorithm parameter used to calculate the background jets
0087     void set_par(float par) { m_par = par; } // default is 0.4
0088     float get_par() { return m_par; }
0089
0090     // set the absolute eta range for the background jet calculation // default is 1.1
0091     void set_tower_abs_eta(float abseta){ m_abs_tower_eta_range = abseta; }
0092     float get_tower_abs_eta() const { return m_abs_tower_eta_range; }
0093
0094     // set the absolute eta range for the background jet calculation // default is 1.1 - m_par (set in the init_algo method)
0095     void set_jet_abs_eta(float abseta){ m_abs_jet_eta_range = abseta; }
0096     float get_jet_abs_eta() const { return m_abs_jet_eta_range; }
0097
0098     // set the number of hardest towers to omit from the background jet calculation // default is 2
0099     void set_omit_nhardest(int omit_nhardest) { m_omit_nhardest = omit_nhardest; }
0100     int get_omit_nhardest() const { return m_omit_nhardest; }
0101
0102     // set the minimum pT for towers and jets used in the background jet calculation
0103     void set_tower_min_pT(float min_pT) { m_tower_min_pT = min_pT; } // default is 0.0
0104     float get_tower_min_pT() const { return m_tower_min_pT; }
0105
0106     // set the minimum pT for towers and jets used in the background jet calculation
0107     void set_jet_min_pT(float min_pT) { m_jet_min_pT = min_pT; } // default is 1.0
0108     float get_jet_min_pT() const { return m_jet_min_pT; }
0109
0110     // set the ghost area for the background jet calculation
0111     void set_ghost_area(float ghost_area) { m_ghost_area = ghost_area; } // default is 0.01
0112     float get_ghost_area() const { return m_ghost_area; }
0113
0114     // print settings
0115     void print_settings(std::ostream& os = std::cout) const;
0116
0117
0118     // tower cut
0119     void do_tower_cut(bool b = true) { m_do_tower_cut = b; }
0120
0121     void set_tower_threshold(float threshold)
0122     {
0123       m_tower_threshold = threshold;
0124     }
0125
0126
0127
0128 private:
0129
0130   // internal methods
0131   int CreateNode(PHCompositeNode *topNode);
0132   void FillNode(PHCompositeNode *topNode,  unsigned int ipos, const float rho, const float sigma);
0133
0134   // variables
0135   std::vector<JetInput *> _inputs;
0136   std::vector<std::string> _outputs;
0137   std::vector<TowerRho::Method> _rho_methods;
0138   Jet::ALGO _algo{Jet::ALGO::KT};
0139
0140   int m_verbosity { 0 };
0141
0142   // settings for the background jet calculation
0143   Jet::ALGO m_bkgd_jet_algo { Jet::ALGO::KT }; // default is KT
0144
0145   float m_par { 0.4 }; // default is 0.4
0146
0147   // eta ranges for the background jet calculation
0148   float m_abs_tower_eta_range { 1.1 }; // default is 1.1
0149   float m_abs_jet_eta_range { 5.0 }; // default is 1.1 - m_par (set in the init_algo method)
0150
0151   float m_abs_ghost_eta {5.0}; // set in the init_algo method
0152
0153   int m_omit_nhardest { 2 };
0154
0155   float m_tower_min_pT { 0.0 };
0156   float m_jet_min_pT { 0.0 };
0157
0158   float m_ghost_area { 0.01 };
0159
0160
0161   // tower threshold
0162   bool m_do_tower_cut { false };
0163   float m_tower_threshold { 0.0 };
0164
0165   fastjet::JetAlgorithm get_fj_algo();
0166   std::string get_fj_algo_name();
0167   std::string m_fj_algo_name;
0168   std::vector<fastjet::PseudoJet> get_pseudojets(PHCompositeNode *topNode);
0169
0170   float _percentile(const std::vector<float> & sorted_vec,
0171         const float percentile,
0172         const float nempty) const ;
0173   void _median_stddev(const std::vector<float> & vec,
0174         float n_empty_jets,
0175         float & median,
0176         float & std_dev) const ;
0177
0178 };
0179
0180 #endif // RHOBASE_DETERMINETOWERRHO_H