packages/CaloReco/RawTowerCombiner.cc

0001 #include "RawTowerCombiner.h"
0002
0003 #include <calobase/RawTower.h>
0004 #include <calobase/RawTowerContainer.h>
0005 #include <calobase/RawTowerDefs.h>
0006 #include <calobase/RawTowerGeomContainer.h>
0007 #include <calobase/RawTowerGeomv1.h>
0008 #include <calobase/RawTowerv1.h>
0009
0010 #include <fun4all/Fun4AllBase.h>
0011 #include <fun4all/Fun4AllReturnCodes.h>
0012 #include <fun4all/SubsysReco.h>
0013
0014 #include <phool/PHCompositeNode.h>
0015 #include <phool/PHNode.h>
0016 #include <phool/PHNodeIterator.h>
0017 #include <phool/getClass.h>
0018
0019 #include <CLHEP/Vector/ThreeVector.h>
0020
0021 #include <algorithm>
0022 #include <cassert>
0023 #include <cmath>
0024 #include <cstdlib>
0025 #include <exception>
0026 #include <iostream>
0027 #include <map>
0028 #include <stdexcept>
0029 #include <utility>
0030 #include <vector>
0031
0032 RawTowerCombiner::RawTowerCombiner(const std::string &name)
0033   : SubsysReco(name)
0034   , _tower_node_prefix("SIM")
0035   , _n_combine_eta(2)
0036   , _n_combine_phi(2)
0037   , _towers(nullptr)
0038   , detector("NONE")
0039 {
0040 }
0041
0042 int RawTowerCombiner::InitRun(PHCompositeNode *topNode)
0043 {
0044   if (_n_combine_eta == 0)
0045   {
0046     std::cout << __PRETTY_FUNCTION__ << " Fatal error _n_combine_eta==0" << std::endl;
0047
0048     exit(1243);
0049   }
0050
0051   if (_n_combine_phi == 0)
0052   {
0053     std::cout << __PRETTY_FUNCTION__ << " Fatal error _n_combine_phi==0" << std::endl;
0054
0055     exit(1243);
0056   }
0057
0058   PHNodeIterator iter(topNode);
0059
0060   // Looking for the DST node
0061   PHCompositeNode *dstNode;
0062   dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode",
0063                                                            "DST"));
0064   if (!dstNode)
0065   {
0066     std::cout << __PRETTY_FUNCTION__ << "DST Node missing, doing nothing."
0067               << std::endl;
0068     exit(1);
0069   }
0070
0071   try
0072   {
0073     CreateNodes(topNode);
0074   }
0075   catch (std::exception &e)
0076   {
0077     std::cout << e.what() << std::endl;
0078     // exit(1);
0079   }
0080
0081   return Fun4AllReturnCodes::EVENT_OK;
0082 }
0083
0084 int RawTowerCombiner::process_event(PHCompositeNode * /*topNode*/)
0085 {
0086   assert(_towers);
0087
0088   const double input_e_sum = _towers->getTotalEdep();
0089   const double input_n_tower = _towers->size();
0090
0091   if (Verbosity())
0092   {
0093     std::cout << __PRETTY_FUNCTION__ << "Process event entered" << std::endl;
0094   }
0095
0096   using tower_id_t = std::pair<int, int>;
0097   using new_tower_map_t = std::map<tower_id_t, RawTower *>;
0098   new_tower_map_t new_tower_map;
0099
0100   RawTowerContainer::ConstRange all_towers = _towers->getTowers();
0101   for (RawTowerContainer::ConstIterator it = all_towers.first;
0102        it != all_towers.second; ++it)
0103   {
0104     const RawTower *input_tower = it->second;
0105     assert(input_tower);
0106
0107     const int intput_eta = input_tower->get_bineta();
0108     const int intput_phi = input_tower->get_binphi();
0109     const int output_eta = get_output_bin_eta(intput_eta);
0110     const int output_phi = get_output_bin_phi(intput_phi);
0111
0112     RawTower *output_tower = nullptr;
0113     new_tower_map_t::iterator it_new = new_tower_map.find(
0114         std::make_pair(output_eta, output_phi));
0115
0116     if (it_new == new_tower_map.end())
0117     {
0118       output_tower = new RawTowerv1(*input_tower);
0119       assert(output_tower);
0120       new_tower_map[std::make_pair(output_eta, output_phi)] = output_tower;
0121
0122       if (Verbosity() >= VERBOSITY_MORE)
0123       {
0124         std::cout << __PRETTY_FUNCTION__ << "::" << detector << "::"
0125                   << " new output tower (prior to tower ID assignment): ";
0126         output_tower->identify();
0127       }
0128     }
0129     else
0130     {
0131       output_tower = it_new->second;
0132
0133       output_tower->set_energy(
0134           output_tower->get_energy() + input_tower->get_energy());
0135
0136       output_tower->set_time(
0137           (std::abs(output_tower->get_energy()) * output_tower->get_time() + std::abs(input_tower->get_energy()) * input_tower->get_time())  //
0138           / (std::abs(output_tower->get_energy()) + std::abs(input_tower->get_energy()) + 1e-9)                                              // avoid devide 0
0139       );
0140
0141       RawTower::CellConstRange cell_range = input_tower->get_g4cells();
0142
0143       for (RawTower::CellConstIterator cell_iter = cell_range.first;
0144            cell_iter != cell_range.second; ++cell_iter)
0145       {
0146         output_tower->add_ecell(cell_iter->first, cell_iter->second);
0147       }
0148
0149       RawTower::ShowerConstRange shower_range =
0150           input_tower->get_g4showers();
0151
0152       for (RawTower::ShowerConstIterator shower_iter = shower_range.first;
0153            shower_iter != shower_range.second; ++shower_iter)
0154       {
0155         output_tower->add_eshower(shower_iter->first,
0156                                   shower_iter->second);
0157       }
0158
0159       if (Verbosity() >= VERBOSITY_MORE)
0160       {
0161         std::cout << __PRETTY_FUNCTION__ << "::" << detector << "::"
0162                   << " merget into output tower (prior to tower ID assignment) : ";
0163         output_tower->identify();
0164       }
0165     }
0166   }
0167
0168   // replace content in tower container
0169   _towers->Reset();
0170
0171   for (auto &it : new_tower_map)
0172   {
0173     const int eta = it.first.first;
0174     const int phi = it.first.second;
0175     RawTower *tower = it.second;
0176
0177     _towers->AddTower(eta, phi, tower);
0178   }
0179
0180   if (Verbosity())
0181   {
0182     std::cout << Name() << "::" << detector << "::" << __PRETTY_FUNCTION__
0183               << "input sum energy = " << input_e_sum << " from " << input_n_tower << " towers, merged sum energy = "
0184               << _towers->getTotalEdep() << " from " << _towers->size() << " towers" << std::endl;
0185     assert(input_e_sum == _towers->getTotalEdep());
0186   }
0187
0188   return Fun4AllReturnCodes::EVENT_OK;
0189 }
0190
0191 int RawTowerCombiner::End(PHCompositeNode * /*topNode*/)
0192 {
0193   return Fun4AllReturnCodes::EVENT_OK;
0194 }
0195
0196 void RawTowerCombiner::CreateNodes(PHCompositeNode *topNode)
0197 {
0198   PHNodeIterator iter(topNode);
0199   PHCompositeNode *runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst(
0200       "PHCompositeNode", "RUN"));
0201   if (!runNode)
0202   {
0203     std::cerr << __PRETTY_FUNCTION__ << "Run Node missing, doing nothing."
0204               << std::endl;
0205     throw std::runtime_error(
0206         "Failed to find Run node in RawTowerCombiner::CreateNodes");
0207   }
0208
0209   const RawTowerDefs::CalorimeterId caloid =
0210       RawTowerDefs::convert_name_to_caloid(detector);
0211
0212   const std::string iTowerGeomNodeName = "TOWERGEOM_" + detector;
0213   RawTowerGeomContainer *towergeom = findNode::getClass<
0214       RawTowerGeomContainer>(topNode, iTowerGeomNodeName);
0215   if (!towergeom)
0216   {
0217     std::cerr << __PRETTY_FUNCTION__ << " - " << iTowerGeomNodeName
0218               << " Node missing, doing nothing." << std::endl;
0219     throw std::runtime_error(
0220         "Failed to find input tower geometry node in RawTowerCombiner::CreateNodes");
0221   }
0222   assert(towergeom);
0223
0224   const double r = towergeom->get_radius();
0225
0226   const int new_phibins = ceil(
0227       double(towergeom->get_phibins()) / double(_n_combine_phi));
0228   const int new_etabins = ceil(
0229       double(towergeom->get_etabins()) / double(_n_combine_eta));
0230
0231   using bound_t = std::pair<double, double>;
0232   using bound_map_t = std::vector<bound_t>;
0233
0234   bound_map_t eta_bound_map;
0235   bound_map_t phi_bound_map;
0236
0237   for (int ibin = 0; ibin < new_phibins; ibin++)
0238   {
0239     const int first_bin = ibin * _n_combine_phi;
0240     assert(first_bin >= 0 && first_bin < towergeom->get_phibins());
0241
0242     int last_bin = ((ibin + 1) * _n_combine_phi) - 1;
0243     last_bin = std::min(last_bin, towergeom->get_phibins());
0244
0245     const std::pair<double, double> range1 = towergeom->get_phibounds(
0246         first_bin);
0247     const std::pair<double, double> range2 = towergeom->get_phibounds(
0248         last_bin);
0249
0250     phi_bound_map.emplace_back(range1.first, range2.second);
0251   }
0252
0253   for (int ibin = 0; ibin < new_etabins; ibin++)
0254   {
0255     const int first_bin = ibin * _n_combine_eta;
0256     assert(first_bin >= 0 && first_bin < towergeom->get_etabins());
0257
0258     int last_bin = ((ibin + 1) * _n_combine_eta) - 1;
0259     last_bin = std::min(last_bin, towergeom->get_etabins());
0260
0261     const std::pair<double, double> range1 = towergeom->get_etabounds(
0262         first_bin);
0263     const std::pair<double, double> range2 = towergeom->get_etabounds(
0264         last_bin);
0265
0266     eta_bound_map.emplace_back(range1.first, range2.second);
0267   }
0268
0269   // now update the tower geometry object with the new tower structure.
0270   towergeom->Reset();
0271
0272   towergeom->set_phibins(new_phibins);
0273   towergeom->set_etabins(new_etabins);
0274
0275   for (int ibin = 0; ibin < new_phibins; ibin++)
0276   {
0277     towergeom->set_phibounds(ibin, phi_bound_map[ibin]);
0278   }
0279   for (int ibin = 0; ibin < new_etabins; ibin++)
0280   {
0281     towergeom->set_etabounds(ibin, eta_bound_map[ibin]);
0282   }
0283
0284   // setup location of all towers
0285   for (int iphi = 0; iphi < towergeom->get_phibins(); iphi++)
0286   {
0287     for (int ieta = 0; ieta < towergeom->get_etabins(); ieta++)
0288     {
0289       RawTowerGeomv1 *tg = new RawTowerGeomv1(
0290           RawTowerDefs::encode_towerid(caloid, ieta, iphi));
0291
0292       CLHEP::Hep3Vector tower_pos;
0293       tower_pos.setRhoPhiEta(r, towergeom->get_phicenter(iphi), towergeom->get_etacenter(ieta));
0294
0295       tg->set_center_x(tower_pos.x());
0296       tg->set_center_y(tower_pos.y());
0297       tg->set_center_z(tower_pos.z());
0298
0299       towergeom->add_tower_geometry(tg);
0300     }
0301   }
0302   if (Verbosity() >= VERBOSITY_SOME)
0303   {
0304     towergeom->identify();
0305   }
0306
0307   const std::string input_TowerNodeName = "TOWER_" + _tower_node_prefix + "_" + detector;
0308   _towers = findNode::getClass<RawTowerContainer>(topNode,
0309                                                   input_TowerNodeName);
0310   if (!_towers)
0311   {
0312     std::cerr << Name() << "::" << detector << "::" << __PRETTY_FUNCTION__
0313               << " " << input_TowerNodeName << " Node missing, doing bail out!"
0314               << std::endl;
0315     throw std::runtime_error(
0316         "Failed to find " + input_TowerNodeName + " node in RawTowerCombiner::CreateNodes");
0317   }
0318
0319   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst(
0320       "PHCompositeNode", "DST"));
0321   if (!dstNode)
0322   {
0323     std::cerr << __PRETTY_FUNCTION__ << "DST Node missing, doing nothing."
0324               << std::endl;
0325     throw std::runtime_error(
0326         "Failed to find DST node in RawTowerCombiner::CreateNodes");
0327   }
0328
0329   return;
0330 }