packages/Prototype3/EventInfoSummary.cc

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