packages/CaloReco/RawClusterBuilderTopo.h

0001 #ifndef CALORECO_RAWCLUSTERTOPO_H
0002 #define CALORECO_RAWCLUSTERTOPO_H
0003
0004 //===========================================================
0005 /// \file RawClusterBuilderTopo.h
0006 /// \brief 3-D topoClustering across calorimeter layers
0007 /// \author Dennis V. Perepelitsa
0008 //===========================================================
0009
0010 #include <fun4all/SubsysReco.h>
0011
0012 #include <map>
0013 #include <string>
0014 #include <utility>  // for pair
0015 #include <vector>
0016
0017 class PHCompositeNode;
0018 class RawClusterContainer;
0019 class RawTowerGeomContainer;
0020
0021 class RawClusterBuilderTopo : public SubsysReco
0022 {
0023  public:
0024   explicit RawClusterBuilderTopo(const std::string &name = "RawClusterBuilderTopo");
0025   ~RawClusterBuilderTopo() override {}
0026
0027   int InitRun(PHCompositeNode *topNode) override;
0028   int process_event(PHCompositeNode *topNode) override;
0029   int End(PHCompositeNode *topNode) override;
0030
0031   void set_nodename(const std::string &nodename)
0032   {
0033     ClusterNodeName = nodename;
0034   }
0035
0036   void set_noise(float noise_0 = 0.0025, float noise_1 = 0.006, float noise_2 = 0.03)
0037   {
0038     _noise_LAYER[0] = noise_0;
0039     _noise_LAYER[1] = noise_1;
0040     _noise_LAYER[2] = noise_2;
0041   }
0042
0043   void set_significance(float seed, float grow, float peri)
0044   {
0045     _sigma_seed = seed;
0046     _sigma_grow = grow;
0047     _sigma_peri = peri;
0048   }
0049
0050   void set_absE(bool allow)
0051   {
0052     _use_absE = allow;
0053   }
0054
0055   void allow_corner_neighbor(bool allow)
0056   {
0057     _allow_corner_neighbor = allow;
0058   }
0059
0060   void set_enable_HCal(bool enable_HCal)
0061   {
0062     _enable_HCal = enable_HCal;
0063   }
0064
0065   void set_enable_EMCal(bool enable_EMCal)
0066   {
0067     _enable_EMCal = enable_EMCal;
0068   }
0069
0070   void set_do_split(bool do_split)
0071   {
0072     _do_split = do_split;
0073   }
0074
0075   void set_minE_local_max(float minE_0 = 1, float minE_1 = 1, float minE_2 = 1)
0076   {
0077     _local_max_minE_LAYER[0] = minE_0;
0078     _local_max_minE_LAYER[1] = minE_1;
0079     _local_max_minE_LAYER[2] = minE_2;
0080   }
0081
0082   void set_R_shower(float R_shower)
0083   {
0084     _R_shower = R_shower;
0085   }
0086
0087   void set_use_only_good_towers(bool b)
0088   {
0089     _only_good_towers = b;
0090   }
0091   bool get_use_only_good_towers()
0092   {
0093     return _only_good_towers;
0094   }
0095
0096   void set_min_cluster_E_saved(float min_cluster_E)
0097   {
0098     _min_cluster_E = min_cluster_E;
0099   }
0100
0101   void setInputTowerNodePrefix(const std::string& inputPrefix)
0102   {
0103     _inputnodeprefix = inputPrefix;
0104   }
0105
0106  private:
0107   void CreateNodes(PHCompositeNode *topNode);
0108
0109   // geometric constants to express IHCal<->EMCal overlap in eta
0110   static int RawClusterBuilderTopo_constants_EMCal_eta_start_given_IHCal[];
0111
0112   static int RawClusterBuilderTopo_constants_EMCal_eta_end_given_IHCal[];
0113
0114   static int RawClusterBuilderTopo_constants_IHCal_eta_given_EMCal[];
0115
0116   // utility functions to express IHCal<->EMCal overlap in phi
0117   int get_first_matching_EMCal_phi_from_IHCal(int index_hcal_phi)
0118   {
0119     return (4 * index_hcal_phi + 5) % _EMCAL_NPHI;
0120   }
0121
0122   int get_matching_HCal_phi_from_EMCal(int index_emcal_phi)
0123   {
0124     return ((index_emcal_phi + 251) / 4) % _HCAL_NPHI;
0125   }
0126
0127   std::vector<int> get_adjacent_towers_by_ID(int ID);
0128
0129   static float calculate_dR(float, float, float, float);
0130
0131   void export_single_cluster(const std::vector<int> &);
0132
0133   void export_clusters(const std::vector<int> &, std::map<int, std::pair<int, int> >, unsigned int, const std::vector<float> &, const std::vector<float> &, const std::vector<float> &);
0134
0135   int get_ID(int ilayer, int ieta, int iphi)
0136   {
0137     if (ilayer < 2)
0138     {
0139       return ilayer * _HCAL_NETA * _HCAL_NPHI + ieta * _HCAL_NPHI + iphi;
0140     }
0141     return _EMCAL_NPHI * _EMCAL_NETA + ieta * _EMCAL_NPHI + iphi;
0142   }
0143
0144   int get_ilayer_from_ID(int ID)
0145   {
0146     if (ID < _EMCAL_NPHI * _EMCAL_NETA)
0147     {
0148       return ((int) (ID / (_HCAL_NETA * _HCAL_NPHI)));
0149     }
0150     else
0151     {
0152       return 2;
0153     }
0154   }
0155
0156   int get_ieta_from_ID(int ID)
0157   {
0158     if (ID < _EMCAL_NPHI * _EMCAL_NETA)
0159     {
0160       return ((int) ((ID % (_HCAL_NETA * _HCAL_NPHI)) / (_HCAL_NPHI)));
0161     }
0162     else
0163     {
0164       return ((int) ((ID - _EMCAL_NPHI * _EMCAL_NETA) / _EMCAL_NPHI));
0165     }
0166   }
0167
0168   int get_iphi_from_ID(int ID)
0169   {
0170     if (ID < _EMCAL_NPHI * _EMCAL_NETA)
0171     {
0172       return ((int) (ID % _HCAL_NPHI));
0173     }
0174     else
0175     {
0176       return ((int) ((ID - _EMCAL_NPHI * _EMCAL_NETA) % _EMCAL_NPHI));
0177     }
0178   }
0179
0180   int get_status_from_ID(int ID)
0181   {
0182     if (ID < _EMCAL_NPHI * _EMCAL_NETA)
0183     {
0184       return _TOWERMAP_STATUS_LAYER_ETA_PHI[get_ilayer_from_ID(ID)][get_ieta_from_ID(ID)][get_iphi_from_ID(ID)];
0185     }
0186     return _EMTOWERMAP_STATUS_ETA_PHI[get_ieta_from_ID(ID)][get_iphi_from_ID(ID)];
0187   }
0188
0189   float get_E_from_ID(int ID)
0190   {
0191     if (ID < _EMCAL_NPHI * _EMCAL_NETA)
0192     {
0193       return _TOWERMAP_E_LAYER_ETA_PHI[get_ilayer_from_ID(ID)][get_ieta_from_ID(ID)][get_iphi_from_ID(ID)];
0194     }
0195     return _EMTOWERMAP_E_ETA_PHI[get_ieta_from_ID(ID)][get_iphi_from_ID(ID)];
0196   }
0197
0198   void set_status_by_ID(int ID, int status)
0199   {
0200     if (ID < _EMCAL_NPHI * _EMCAL_NETA)
0201     {
0202       _TOWERMAP_STATUS_LAYER_ETA_PHI[get_ilayer_from_ID(ID)][get_ieta_from_ID(ID)][get_iphi_from_ID(ID)] = status;
0203     }
0204     else
0205     {
0206       _EMTOWERMAP_STATUS_ETA_PHI[get_ieta_from_ID(ID)][get_iphi_from_ID(ID)] = status;
0207     }
0208   }
0209
0210   RawClusterContainer *_clusters {nullptr};
0211
0212   RawTowerGeomContainer *_geom_containers[3]{};
0213
0214   // geometric parameters defined at runtime
0215   int _EMCAL_NETA {-1};
0216   int _EMCAL_NPHI{-1};
0217
0218   int _HCAL_NETA{-1};
0219   int _HCAL_NPHI{-1};
0220
0221   float _noise_LAYER[3]{};
0222
0223   float _sigma_seed {4.0};
0224   float _sigma_grow {2.0};
0225   float _sigma_peri {0.0};
0226   float _local_max_minE_LAYER[3]{};
0227   float _R_shower {0.025};
0228   float _min_cluster_E{0.0};
0229
0230   bool _allow_corner_neighbor {true};
0231   bool _use_absE {true};
0232
0233   bool _enable_HCal {true};
0234   bool _enable_EMCal {true};
0235
0236   bool _do_split {true};
0237   bool _only_good_towers {true};
0238
0239   std::vector<std::vector<std::vector<float> > > _TOWERMAP_E_LAYER_ETA_PHI;
0240   std::vector<std::vector<std::vector<int> > > _TOWERMAP_KEY_LAYER_ETA_PHI;
0241   std::vector<std::vector<std::vector<int> > > _TOWERMAP_STATUS_LAYER_ETA_PHI;
0242
0243   std::vector<std::vector<float> > _EMTOWERMAP_E_ETA_PHI;
0244   std::vector<std::vector<int> > _EMTOWERMAP_KEY_ETA_PHI;
0245   std::vector<std::vector<int> > _EMTOWERMAP_STATUS_ETA_PHI;
0246
0247   std::string _inputnodeprefix;
0248   std::string ClusterNodeName {"TOPOCLUSTER_HCAL"};
0249 };
0250
0251 #endif