packages/CaloEmbedding/caloTowerEmbed.cc

0001 #include "caloTowerEmbed.h"
0002
0003 #include <calobase/RawTowerGeom.h>
0004 #include <calobase/RawTowerGeomContainer.h>
0005 #include <calobase/TowerInfo.h>  // for TowerInfo
0006 #include <calobase/TowerInfoContainer.h>
0007 #include <calobase/TowerInfoContainerv2.h>
0008 #include <calobase/TowerInfoDefs.h>
0009
0010 #include <ffamodules/CDBInterface.h>
0011
0012 #include <ffaobjects/EventHeader.h>
0013
0014 #include <fun4all/Fun4AllReturnCodes.h>
0015 #include <fun4all/Fun4AllServer.h>
0016 #include <fun4all/SubsysReco.h>  // for SubsysReco
0017
0018 #include <phool/PHCompositeNode.h>
0019 #include <phool/PHIODataNode.h>    // for PHIODataNode
0020 #include <phool/PHNode.h>          // for PHNode
0021 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
0022 #include <phool/PHObject.h>        // for PHObject
0023 #include <phool/getClass.h>
0024 #include <phool/phool.h>
0025 #include <phool/recoConsts.h>
0026
0027 #include <cstdlib>    // for exit
0028 #include <exception>  // for exception
0029 #include <iostream>   // for operator<<, basic_ostream
0030 #include <stdexcept>  // for runtime_error
0031
0032 //____________________________________________________________________________..
0033 caloTowerEmbed::caloTowerEmbed(const std::string &name)
0034   : SubsysReco(name)
0035 {
0036   std::cout << "caloTowerEmbed::caloTowerEmbed(const std::string &name) Calling ctor" << std::endl;
0037 }
0038
0039 //____________________________________________________________________________..
0040 caloTowerEmbed::~caloTowerEmbed()
0041 {
0042   std::cout << "caloTowerEmbed::~caloTowerEmbed() Calling dtor" << std::endl;
0043 }
0044
0045 //____________________________________________________________________________..
0046 int caloTowerEmbed::InitRun(PHCompositeNode *topNode)
0047 {
0048   if (m_dettype == CaloTowerDefs::CEMC)
0049   {
0050     m_detector = "CEMC";
0051   }
0052   else if (m_dettype == CaloTowerDefs::HCALIN)
0053   {
0054     m_detector = "HCALIN";
0055   }
0056   else if (m_dettype == CaloTowerDefs::HCALOUT)
0057   {
0058     m_detector = "HCALOUT";
0059   }
0060   else if (m_dettype == CaloTowerDefs::ZDC)
0061   {
0062     m_detector = "ZDC";
0063   }
0064   else if (m_dettype == CaloTowerDefs::SEPD)
0065   {
0066     m_detector = "SEPD";
0067   }
0068
0069   try
0070   {
0071     CreateNodeTree(topNode);
0072   }
0073   catch (std::exception &e)
0074   {
0075     std::cout << e.what() << std::endl;
0076     return Fun4AllReturnCodes::ABORTRUN;
0077   }
0078
0079   return Fun4AllReturnCodes::EVENT_OK;
0080 }
0081
0082 //____________________________________________________________________________..
0083 int caloTowerEmbed::process_event(PHCompositeNode * topNode)
0084 {
0085   ++m_eventNumber;
0086
0087   if (Verbosity())
0088   {
0089     std::cout << "event " << m_eventNumber << " working on " << m_detector << std::endl;
0090   }
0091
0092   if (m_embedwaveform)
0093   {
0094     std::string ped_nodename = "PEDESTAL_" + m_detector;
0095     m_PedestalContainer = findNode::getClass<TowerInfoContainer>(topNode, ped_nodename);
0096
0097     if (!m_PedestalContainer)
0098     {
0099       std::cout << PHWHERE << " " << ped_nodename << " Node missing, doing nothing." << std::endl;
0100       return Fun4AllReturnCodes::ABORTEVENT;
0101     }
0102   }
0103
0104
0105   unsigned int ntowers = _data_towers->size();
0106   for (unsigned int channel = 0; channel < ntowers; channel++)
0107   {
0108
0109
0110     _sim_towers->get_tower_at_channel(channel)->set_status(_data_towers->get_tower_at_channel(channel)->get_status());
0111
0112     TowerInfo *caloinfo_data = _data_towers->get_tower_at_channel(channel);
0113     TowerInfo *caloinfo_sim = _sim_towers->get_tower_at_channel(channel);
0114
0115     if (m_embedwaveform)
0116     {
0117       // when the data is not ZS-ed
0118       if (!(caloinfo_data->get_isZS() || caloinfo_data->get_isNotInstr()))
0119       {
0120         // here we really don't want the m_samples being greater than what data has!
0121         for (int j = 0; j < m_nsamples; j++)
0122         {
0123           // superpose the waveforms
0124           caloinfo_sim->set_waveform_value(j, caloinfo_data->get_waveform_value(j) + caloinfo_sim->get_waveform_value(j));
0125         }
0126       }
0127       else
0128       {
0129         // if this is ZS-ed or empty(like the packet is gone)
0130         // we add the noise pedestal and add the post - pre to sample 6
0131         TowerInfo *pedestal_tower = m_PedestalContainer->get_tower_at_channel(channel);
0132         float pedestal_mean = 0;
0133         std::vector<float> m_waveform_pedestal;
0134         m_waveform_pedestal.resize(m_nsamples);
0135         // this is for pedestal scaling setting
0136         for (int j = 0; j < m_nsamples; j++)
0137         {
0138           m_waveform_pedestal.at(j) = (j < m_datasamples) ? pedestal_tower->get_waveform_value(j) : pedestal_tower->get_waveform_value(m_datasamples - 1);
0139           pedestal_mean += m_waveform_pedestal.at(j);
0140         }
0141         pedestal_mean /= m_nsamples;
0142         for (int j = 0; j < m_nsamples; j++)
0143         {
0144           m_waveform_pedestal.at(j) = (m_waveform_pedestal.at(j) - pedestal_mean) * m_pedestal_scale + pedestal_mean;
0145         }
0146         // add the pedestal
0147         for (int j = 0; j < m_nsamples; j++)
0148         {
0149           // superpose the waveforms
0150
0151           caloinfo_sim->set_waveform_value(j, caloinfo_sim->get_waveform_value(j) + m_waveform_pedestal.at(j));
0152         }
0153         // add the post - pre to sample 6
0154         int post_pre = caloinfo_data->get_waveform_value(1) - caloinfo_data->get_waveform_value(0);
0155
0156         caloinfo_sim->set_waveform_value(6, caloinfo_sim->get_waveform_value(6) + post_pre);
0157       }
0158     }
0159     else
0160     {
0161       float data_E = caloinfo_data->get_energy();
0162       float sim_E = caloinfo_sim->get_energy();
0163       float embed_E = data_E + sim_E;
0164
0165       _sim_towers->get_tower_at_channel(channel)->set_energy(embed_E);
0166       _sim_towers->get_tower_at_channel(channel)->set_time(caloinfo_data->get_time());
0167     }
0168
0169   }  // end loop over channels
0170
0171   return Fun4AllReturnCodes::EVENT_OK;
0172 }
0173
0174 int caloTowerEmbed::End(PHCompositeNode * /*topNode*/)
0175 {
0176   return Fun4AllReturnCodes::EVENT_OK;
0177 }
0178
0179 void caloTowerEmbed::CreateNodeTree(PHCompositeNode *topNode)
0180 {
0181   std::cout << "creating node" << std::endl;
0182
0183   std::string TowerNodeName = m_inputNodePrefix + m_detector;
0184   std::string SimTowerNodeName = TowerNodeName;
0185   if (m_embedwaveform)
0186   {
0187     SimTowerNodeName = m_waveformNodePrefix + m_detector;
0188   }
0189   std::string GeomNodeName = "TOWERGEOM_" + m_detector;
0190   if (m_useRetower && m_detector == "CEMC")
0191   {
0192     TowerNodeName = m_inputNodePrefix + m_detector + "_RETOWER";
0193     GeomNodeName = "TOWERGEOM_HCALIN";
0194   }
0195
0196   Fun4AllServer *se = Fun4AllServer::instance();
0197   PHCompositeNode *dataTopNode = se->topNode("TOPData");
0198
0199   tower_geom = findNode::getClass<RawTowerGeomContainer>(topNode, GeomNodeName);
0200   if (!tower_geom)
0201   {
0202     std::cerr << Name() << "::" << m_detector << "::" << __PRETTY_FUNCTION__
0203               << "tower geom " << GeomNodeName << " missing, doing nothing." << std::endl;
0204     throw std::runtime_error("Failed to find " + GeomNodeName + " node");
0205   }
0206
0207   PHNodeIterator simIter(topNode);
0208   PHNodeIterator dataIter(dataTopNode);
0209
0210   // data top node first
0211
0212   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(dataIter.findFirst(
0213       "PHCompositeNode", "DST"));
0214   if (!dstNode)
0215   {
0216     std::cerr << Name() << "::" << m_detector << "::" << __PRETTY_FUNCTION__
0217               << "DST Node missing, doing nothing." << std::endl;
0218     throw std::runtime_error(
0219         "Failed to find DST node in RawTowerCalibration::CreateNodes");
0220   }
0221
0222   PHCompositeNode *dstNodeSim = dynamic_cast<PHCompositeNode *>(simIter.findFirst(
0223       "PHCompositeNode", "DST"));
0224   if (!dstNodeSim)
0225   {
0226     std::cerr << Name() << "::" << m_detector << "::" << __PRETTY_FUNCTION__
0227               << "DSTsim Node missing, doing nothing." << std::endl;
0228     throw std::runtime_error(
0229         "Failed to find DSTsim node in RawTowerCalibration::CreateNodes");
0230   }
0231
0232   // data
0233   _data_towers = findNode::getClass<TowerInfoContainer>(dstNode, TowerNodeName);
0234   if (!_data_towers)
0235   {
0236     std::cerr << Name() << "::" << m_detector << "::" << __PRETTY_FUNCTION__
0237               << TowerNodeName << " Node missing, doing nothing." << std::endl;
0238     throw std::runtime_error(
0239         "Failed to find " + TowerNodeName + " node in caloTowerEmbed::CreateNodes");
0240   }
0241
0242   // sim
0243   _sim_towers = findNode::getClass<TowerInfoContainer>(dstNodeSim, SimTowerNodeName);
0244   if (!_sim_towers)
0245   {
0246     std::cerr << Name() << "::" << m_detector << "::" << __PRETTY_FUNCTION__
0247               << TowerNodeName << " Sim Node missing, doing nothing." << std::endl;
0248     throw std::runtime_error(
0249         "Failed to find " + TowerNodeName + " Sim node in caloTowerEmbed::CreateNodes");
0250   }
0251
0252   PHNodeIterator dstIterSim(dstNodeSim);
0253   PHCompositeNode *caloNode = dynamic_cast<PHCompositeNode *>(dstIterSim.findFirst("PHCompositeNode", m_detector));
0254
0255   if (!caloNode)
0256   {
0257     caloNode = new PHCompositeNode(m_detector);
0258     dstNodeSim->addNode(caloNode);
0259   }
0260
0261   return;
0262 }