EMCTemplateFitting/src/EventInfoSummary.cc

0001 #include "EventInfoSummary.h"
0002
0003 #include "RawTower_Prototype4.h"
0004
0005 #include <calobase/RawTower.h>  // for RawTower
0006 #include <calobase/RawTowerContainer.h>
0007
0008 #include <phparameter/PHParameters.h>
0009
0010 #include <pdbcalbase/PdbParameterMap.h>
0011
0012 #include <fun4all/Fun4AllBase.h>  // for Fun4AllB...
0013 #include <fun4all/Fun4AllReturnCodes.h>
0014 #include <fun4all/SubsysReco.h>  // for SubsysReco
0015
0016 #include <phool/PHCompositeNode.h>
0017 #include <phool/PHIODataNode.h>    // for PHIOData...
0018 #include <phool/PHNodeIterator.h>  // for PHNodeIt...
0019 #include <phool/getClass.h>
0020
0021 #include <Event/Event.h>
0022 #include <Event/EventTypes.h>
0023 #include <Event/packet.h>
0024
0025 #include <boost/accumulators/accumulators.hpp>
0026 #include <boost/accumulators/statistics.hpp>
0027
0028 #include <cassert>
0029 #include <cmath>  // for NAN
0030 #include <iostream>
0031 #include <string>
0032 #include <utility>  // for pair
0033
0034 // with the exception of std,  please no using namespaces
0035 // it severely muddies the ability to figure out
0036 // where things come from
0037 using namespace std;
0038
0039 //____________________________________
0040 EventInfoSummary::EventInfoSummary()
0041   : SubsysReco("EventInfoSummary")
0042   , eventinfo_node_name("EVENT_INFO")
0043 {
0044 }
0045
0046 //_____________________________________
0047 int EventInfoSummary::InitRun(PHCompositeNode *topNode)
0048 {
0049   CreateNodeTree(topNode);
0050   return Fun4AllReturnCodes::EVENT_OK;
0051 }
0052
0053 //____________________________________
0054 int EventInfoSummary::process_event(PHCompositeNode *topNode)
0055 {
0056   Event *event = findNode::getClass<Event>(topNode, "PRDF");
0057   if (!event)
0058   {
0059     if (Verbosity() >= VERBOSITY_SOME)
0060       cout << "EventInfoSummary::Process_Event - Event not found" << endl;
0061     return Fun4AllReturnCodes::DISCARDEVENT;
0062   }
0063
0064   // search for run info
0065   if (event->getEvtType() != DATAEVENT)
0066     return Fun4AllReturnCodes::EVENT_OK;
0067   else  // DATAEVENT
0068   {
0069     if (Verbosity() >= VERBOSITY_SOME)
0070     {
0071       cout << "EventInfoSummary::process_event - with DATAEVENT events ";
0072       event->identify();
0073     }
0074
0075     map<int, Packet *> packet_list;
0076
0077     PHParameters Params("EventInfo");
0078
0079     // spill indicator
0080     {
0081       RawTowerContainer *TOWER_RAW_SPILL_WARBLER =
0082           findNode::getClass<RawTowerContainer>(topNode,
0083                                                 "TOWER_RAW_SPILL_WARBLER");
0084       assert(TOWER_RAW_SPILL_WARBLER);
0085
0086       RawTower_Prototype4 *raw_tower = dynamic_cast<RawTower_Prototype4 *>(
0087           TOWER_RAW_SPILL_WARBLER->getTower(0));
0088       assert(raw_tower);
0089
0090       boost::accumulators::accumulator_set<double, boost::accumulators::features<boost::accumulators::tag::variance>> acc;
0091
0092       for (int i = 0; i < RawTower_Prototype4::NSAMPLES; i++)
0093       {
0094         acc(raw_tower->get_signal_samples(i));
0095       }
0096
0097       const double warbler_rms = boost::accumulators::variance(acc);
0098       const bool is_in_spill = warbler_rms > (1000 * 1000);
0099       Params.set_double_param("beam_SPILL_WARBLER_RMS", warbler_rms);
0100       Params.set_double_param("beam_Is_In_Spill", is_in_spill);
0101       Params.set_int_param("beam_Is_In_Spill", is_in_spill);
0102     }
0103
0104     // energy sums
0105     {
0106       RawTowerContainer *TOWER_CALIB_CEMC =
0107           findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC");
0108       assert(TOWER_CALIB_CEMC);
0109
0110       RawTowerContainer *TOWER_CALIB_LG_HCALIN =
0111           findNode::getClass<RawTowerContainer>(topNode,
0112                                                 "TOWER_CALIB_LG_HCALIN");
0113
0114       RawTowerContainer *TOWER_CALIB_LG_HCALOUT =
0115           findNode::getClass<RawTowerContainer>(topNode,
0116                                                 "TOWER_CALIB_LG_HCALOUT");
0117
0118       // process inner HCAL
0119       if (TOWER_CALIB_CEMC)
0120       {
0121         double sum_energy_calib = 0;
0122
0123         auto range = TOWER_CALIB_CEMC->getTowers();
0124         for (auto it = range.first; it != range.second; ++it)
0125         {
0126           RawTower *tower = it->second;
0127           assert(tower);
0128
0129           const int col = tower->get_bineta();
0130           const int row = tower->get_binphi();
0131
0132           if (col < 0 or col >= 8)
0133             continue;
0134           if (row < 0 or row >= 8)
0135             continue;
0136
0137           const double energy_calib = tower->get_energy();
0138           sum_energy_calib += energy_calib;
0139
0140         }  //       for (auto it = range.first; it != range.second; ++it)
0141         Params.set_double_param("CALIB_CEMC_Sum", sum_energy_calib);
0142       }  // process inner HCAL
0143
0144       // process inner HCAL
0145       if (TOWER_CALIB_LG_HCALIN)
0146       {
0147         double sum_energy_calib = 0;
0148
0149         auto range = TOWER_CALIB_LG_HCALIN->getTowers();
0150         for (auto it = range.first; it != range.second; ++it)
0151         {
0152           RawTower *tower = it->second;
0153           assert(tower);
0154
0155           const int col = tower->get_bineta();
0156           const int row = tower->get_binphi();
0157
0158           if (col < 0 or col >= 4)
0159             continue;
0160           if (row < 0 or row >= 4)
0161             continue;
0162
0163           const double energy_calib = tower->get_energy();
0164           sum_energy_calib += energy_calib;
0165
0166         }  //       for (auto it = range.first; it != range.second; ++it)
0167         Params.set_double_param("CALIB_LG_HCALIN_Sum", sum_energy_calib);
0168       }  // process inner HCAL
0169
0170       // process outer HCAL
0171       if (TOWER_CALIB_LG_HCALOUT)
0172       {
0173         double sum_energy_calib = 0;
0174
0175         auto range = TOWER_CALIB_LG_HCALOUT->getTowers();
0176         for (auto it = range.first; it != range.second; ++it)
0177         {
0178           RawTower *tower = it->second;
0179           assert(tower);
0180
0181           const int col = tower->get_bineta();
0182           const int row = tower->get_binphi();
0183
0184           if (col < 0 or col >= 4)
0185             continue;
0186           if (row < 0 or row >= 4)
0187             continue;
0188
0189           const double energy_calib = tower->get_energy();
0190           sum_energy_calib += energy_calib;
0191
0192         }  //       for (auto it = range.first; it != range.second; ++it)
0193
0194         Params.set_double_param("CALIB_LG_HCALOUT_Sum", sum_energy_calib);
0195       }  // process outer HCAL
0196     }
0197
0198     // generic packets
0199     for (typ_channel_map::const_iterator it = channel_map.begin();
0200          it != channel_map.end(); ++it)
0201     {
0202       const string &name = it->first;
0203       const channel_info &info = it->second;
0204
0205       if (packet_list.find(info.packet_id) == packet_list.end())
0206       {
0207         packet_list[info.packet_id] = event->getPacket(info.packet_id);
0208       }
0209
0210       Packet *packet = packet_list[info.packet_id];
0211
0212       if (!packet)
0213       {
0214         //          if (Verbosity() >= VERBOSITY_SOME)
0215         cout << "EventInfoSummary::process_event - failed to locate packet "
0216              << info.packet_id << " from ";
0217         event->identify();
0218
0219         Params.set_double_param(name, NAN);
0220         continue;
0221       }
0222
0223       const int ivalue = packet->iValue(info.offset);
0224
0225       const double dvalue = ivalue * info.calibration_const;
0226
0227       if (Verbosity() >= VERBOSITY_SOME)
0228       {
0229         cout << "EventInfoSummary::process_event - " << name << " = " << dvalue
0230              << ", raw = " << ivalue << " @ packet " << info.packet_id
0231              << ", offset " << info.offset << endl;
0232       }
0233
0234       Params.set_double_param(name, dvalue);
0235     }
0236
0237     for (map<int, Packet *>::iterator it = packet_list.begin();
0238          it != packet_list.end(); ++it)
0239     {
0240       if (it->second)
0241         delete it->second;
0242     }
0243
0244     Params.SaveToNodeTree(topNode, eventinfo_node_name);
0245
0246     if (Verbosity() >= VERBOSITY_SOME)
0247       Params.Print();
0248   }
0249   return Fun4AllReturnCodes::EVENT_OK;
0250 }
0251
0252 //_______________________________________
0253 void EventInfoSummary::CreateNodeTree(PHCompositeNode *topNode)
0254 {
0255   PHNodeIterator nodeItr(topNode);
0256
0257   PHNodeIterator iter(topNode);
0258
0259   // DST node
0260   PHCompositeNode *dst_node =
0261       static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0262   if (!dst_node)
0263   {
0264     cout << "PHComposite node created: DST" << endl;
0265     dst_node = new PHCompositeNode("DST");
0266     topNode->addNode(dst_node);
0267   }
0268
0269   PdbParameterMap *nodeparams =
0270       findNode::getClass<PdbParameterMap>(dst_node, eventinfo_node_name);
0271   if (not nodeparams)
0272   {
0273     dst_node->addNode(new PHIODataNode<PdbParameterMap>(new PdbParameterMap(),
0274                                                         eventinfo_node_name));
0275   }
0276 }
0277
0278 void EventInfoSummary::add_channel(
0279     const std::string &name,        //! name of the channel
0280     const int packet_id,            //! packet id
0281     const unsigned int offset,      //! offset in packet data
0282     const double calibration_const  //! conversion constant from integer to
0283                                     //! meaningful value
0284 )
0285 {
0286   channel_map.insert(
0287       make_pair(name, channel_info(packet_id, offset, calibration_const)));
0288 }