packages/epd/EpdReco.cc

0001 #include "EpdReco.h"
0002
0003 #include "EpdGeomV2.h"
0004
0005 #include <calobase/TowerInfo.h>
0006 #include <calobase/TowerInfoContainer.h>
0007 #include <calobase/TowerInfoContainerv4.h>
0008 #include <calobase/TowerInfoDefs.h>
0009
0010 #include <cdbobjects/CDBTTree.h>  // for CDBTTree
0011
0012 #include <ffamodules/CDBInterface.h>
0013
0014 #include <fun4all/Fun4AllReturnCodes.h>
0015 #include <fun4all/SubsysReco.h>  // for SubsysReco
0016
0017 #include <phool/PHCompositeNode.h>
0018 #include <phool/PHIODataNode.h>
0019 #include <phool/PHNode.h>  // for PHNode
0020 #include <phool/PHNodeIterator.h>
0021 #include <phool/PHObject.h>  // for PHObject
0022 #include <phool/getClass.h>
0023 #include <phool/phool.h>  // for PHWHERE
0024
0025 #include <TSystem.h>
0026
0027 #include <array>    // for array
0028 #include <cstdlib>  // for exit
0029 #include <iostream>
0030
0031 EpdReco::EpdReco(const std::string &name)
0032   : SubsysReco(name)
0033 {
0034   FillTilePhi0Array();
0035   FillTilePhiArray();
0036 }
0037
0038 EpdReco::~EpdReco()
0039 {
0040   delete cdbttree;
0041 }
0042
0043 int EpdReco::InitRun(PHCompositeNode *topNode)
0044 {
0045   if (!m_overrideCalibName)
0046   {
0047     m_calibName = "SEPD_NMIP_CALIB";
0048   }
0049   if (!m_overrideFieldName)
0050   {
0051     m_fieldname = "sepd_calib";
0052   }
0053   std::string calibdir = CDBInterface::instance()->getUrl(m_calibName);
0054   if (!calibdir.empty())
0055   {
0056     cdbttree = new CDBTTree(calibdir);
0057   }
0058   else
0059   {
0060     std::cout << "EpdReco::::InitRun No calibration file for domain "
0061               << m_calibName << " found" << std::endl;
0062     exit(1);
0063   }
0064
0065   CreateNodes(topNode);
0066
0067   PHNodeIterator node_itr(topNode);
0068   PHCompositeNode *runNode = dynamic_cast<PHCompositeNode *>(node_itr.findFirst("PHCompositeNode", "RUN"));
0069   if (!runNode)
0070   {
0071     std::cout << PHWHERE << "RUN Node not found - that is fatal" << std::endl;
0072     gSystem->Exit(1);
0073     exit(1);
0074   }
0075
0076   PHNodeIterator runiter(runNode);
0077   PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(runiter.findFirst("PHCompositeNode", m_Detector));
0078   if (!DetNode)
0079   {
0080     DetNode = new PHCompositeNode(m_Detector);
0081     runNode->addNode(DetNode);
0082   }
0083
0084   EpdGeom *epdGeom = findNode::getClass<EpdGeom>(topNode, "TOWERGEOM_EPD");
0085   if (!epdGeom)
0086   {
0087     epdGeom = new EpdGeomV2();
0088     PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(epdGeom, "TOWERGEOM_EPD", "PHObject");
0089     DetNode->addNode(newNode);
0090   }
0091
0092   // fill epd geometry
0093   unsigned int epdchannels = 744;
0094   for (unsigned int ch = 0; ch < epdchannels; ch++)
0095   {
0096     unsigned int thiskey = TowerInfoDefs::encode_epd(ch);
0097     epdGeom->set_z(thiskey, GetTileZ(TowerInfoDefs::get_epd_arm(thiskey)));
0098     epdGeom->set_r(thiskey, GetTileR(TowerInfoDefs::get_epd_rbin(thiskey)));
0099     if (TowerInfoDefs::get_epd_rbin(thiskey) == 0)
0100     {
0101       epdGeom->set_phi0(thiskey, GetTilePhi0(TowerInfoDefs::get_epd_phibin(thiskey)));
0102     }
0103     else
0104     {
0105       epdGeom->set_phi(thiskey, GetTilePhi(TowerInfoDefs::get_epd_phibin(thiskey)));
0106     }
0107   }
0108
0109   return Fun4AllReturnCodes::EVENT_OK;
0110 }
0111
0112 int EpdReco::process_event(PHCompositeNode *topNode)
0113 {
0114   if (Verbosity() > 1)
0115   {
0116     std::cout << "EpdReco::process_event -- entered" << std::endl;
0117   }
0118
0119   TowerInfoContainer *_sepd_towerinfo = findNode::getClass<TowerInfoContainer>(topNode, "TOWERS_SEPD");
0120   unsigned int ntowers = 0;
0121   if (_sepd_towerinfo)
0122   {
0123     ntowers = _sepd_towerinfo->size();
0124   }
0125   if (ntowers != 744)
0126   {
0127     std::cout << "sEPD container has unexpected size - exiting now!"
0128               << std::endl;
0129     exit(1);
0130   }
0131
0132   TowerInfoContainer *m_TowerInfoContainer_calib = findNode::getClass<TowerInfoContainer>(topNode, m_TowerInfoNodeName_calib);
0133   if (!m_TowerInfoContainer_calib)
0134   {
0135     std::cout << PHWHERE << "Could not locate TowerInfoContainer node "
0136               << m_TowerInfoNodeName_calib << std::endl;
0137     exit(1);
0138   }
0139
0140   for (unsigned int ch = 0; ch < ntowers; ch++)
0141   {
0142     float ch_time = _sepd_towerinfo->get_tower_at_channel(ch)->get_time_float();
0143     float ch_adc = _sepd_towerinfo->get_tower_at_channel(ch)->get_energy();
0144     float ch_mpv = cdbttree->GetFloatValue(ch, m_fieldname);
0145     double ch_nmip = ch_adc / ch_mpv;
0146     m_TowerInfoContainer_calib->get_tower_at_channel(ch)->set_energy(ch_nmip);
0147     m_TowerInfoContainer_calib->get_tower_at_channel(ch)->set_time_float(ch_time);
0148   }
0149
0150   return Fun4AllReturnCodes::EVENT_OK;
0151 }
0152
0153 float EpdReco::GetTilePhi(int thisphi)
0154 {
0155   return tilephi.at(thisphi);
0156 }
0157
0158 float EpdReco::GetTilePhi0(int thisphi0)
0159 {
0160   return tilephi0.at(thisphi0);
0161 }
0162
0163 float EpdReco::GetTileR(int thisr)
0164 {
0165   static const std::array<float, 16> tileR{
0166       6.8, 11.2, 15.6, 20.565, 26.095, 31.625, 37.155, 42.685,
0167       48.215, 53.745, 59.275, 64.805, 70.335, 75.865, 81.395, 86.925};
0168   return tileR.at(thisr);
0169 }
0170
0171 float EpdReco::GetTileZ(int thisz)
0172 {
0173   static const std::array<float, 2> tileZ{-316.0, 316.0};
0174   return tileZ.at(thisz);
0175 }
0176
0177 void EpdReco::CreateNodes(PHCompositeNode *topNode)
0178 {
0179   PHNodeIterator iter(topNode);
0180   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0181   if (!dstNode)
0182   {
0183     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0184     gSystem->Exit(1);
0185     exit(1);
0186   }
0187
0188   PHNodeIterator dstiter(dstNode);
0189   PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", m_Detector));
0190   if (!DetNode)
0191   {
0192     DetNode = new PHCompositeNode(m_Detector);
0193     dstNode->addNode(DetNode);
0194   }
0195
0196   TowerInfoContainer *m_TowerInfoContainer_calib = findNode::getClass<TowerInfoContainer>(DetNode, m_TowerInfoNodeName_calib);
0197   if (m_TowerInfoContainer_calib == nullptr)
0198   {
0199     m_TowerInfoContainer_calib = new TowerInfoContainerv4(TowerInfoContainer::DETECTOR::SEPD);
0200     PHIODataNode<PHObject> *TowerInfoNodecalib = new PHIODataNode<PHObject>(m_TowerInfoContainer_calib, m_TowerInfoNodeName_calib, "PHObject");
0201     DetNode->addNode(TowerInfoNodecalib);
0202   }
0203
0204   return;
0205 }
0206
0207 void EpdReco::FillTilePhiArray()
0208 {
0209   size_t i = 0;
0210   for (float &iter : tilephi)
0211   {
0212     iter = ((2 * i + 1) * M_PI) / 24;
0213     i++;
0214   }
0215   return;
0216 }
0217
0218 void EpdReco::FillTilePhi0Array()
0219 {
0220   size_t i = 0;
0221   for (float &iter : tilephi0)
0222   {
0223     iter = ((2 * i + 1) * M_PI) / 12;
0224     i++;
0225   }
0226   return;
0227 }