packages/Prototype3/CaloCalibration.cc

0001 #include "CaloCalibration.h"
0002
0003 #include "PROTOTYPE3_FEM.h"
0004 #include "RawTower_Prototype3.h"
0005
0006 #include <calobase/RawTower.h>  // for RawTower
0007 #include <calobase/RawTowerContainer.h>
0008 #include <calobase/RawTowerDefs.h>  // for keytype
0009
0010 #include <phparameter/PHParameters.h>  // for PHParameters
0011
0012 #include <fun4all/Fun4AllReturnCodes.h>
0013 #include <fun4all/SubsysReco.h>  // for SubsysReco
0014
0015 #include <phool/PHCompositeNode.h>
0016 #include <phool/PHIODataNode.h>    // for PHIODataNode
0017 #include <phool/PHNode.h>          // for PHNode
0018 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
0019 #include <phool/PHObject.h>        // for PHObject
0020 #include <phool/getClass.h>
0021
0022 #include <boost/format.hpp>
0023
0024 #include <cassert>
0025 #include <cmath>  // for NAN, isnan
0026 #include <iostream>
0027 #include <map>        // for _Rb_tree_iterator
0028 #include <stdexcept>  // for runtime_error
0029 #include <string>
0030 #include <utility>  // for pair
0031 #include <vector>   // for vector
0032
0033 using namespace std;
0034
0035 //____________________________________
0036 CaloCalibration::CaloCalibration(const std::string &name)
0037   :  //
0038   SubsysReco(string("CaloCalibration_") + name)
0039   ,  //
0040   _calib_towers(nullptr)
0041   , _raw_towers(nullptr)
0042   , detector(name)
0043   ,  //
0044   _calib_tower_node_prefix("CALIB")
0045   ,  //
0046   _raw_tower_node_prefix("RAW")
0047   ,  //
0048   _calib_params(name)
0049 {
0050   SetDefaultParameters(_calib_params);
0051 }
0052
0053 //____________________________________
0054 int CaloCalibration::Init(PHCompositeNode *topNode)
0055 {
0056   return Fun4AllReturnCodes::EVENT_OK;
0057 }
0058
0059 //_____________________________________
0060 int CaloCalibration::InitRun(PHCompositeNode *topNode)
0061 {
0062   CreateNodeTree(topNode);
0063
0064   if (Verbosity())
0065   {
0066     std::cout << Name() << "::" << detector << "::" << __PRETTY_FUNCTION__
0067               << " - print calibration parameters: " << endl;
0068     _calib_params.Print();
0069   }
0070
0071   return Fun4AllReturnCodes::EVENT_OK;
0072 }
0073
0074 //____________________________________
0075 int CaloCalibration::process_event(PHCompositeNode *topNode)
0076 {
0077   if (Verbosity())
0078   {
0079     std::cout << Name() << "::" << detector << "::" << __PRETTY_FUNCTION__
0080               << "Process event entered" << std::endl;
0081   }
0082
0083   const double calib_const_scale =
0084       _calib_params.get_double_param("calib_const_scale");
0085   const bool use_chan_calibration =
0086       _calib_params.get_int_param("use_chan_calibration") > 0;
0087
0088   RawTowerContainer::Range begin_end = _raw_towers->getTowers();
0089   RawTowerContainer::Iterator rtiter;
0090   for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
0091   {
0092     RawTowerDefs::keytype key = rtiter->first;
0093     RawTower_Prototype3 *raw_tower =
0094         dynamic_cast<RawTower_Prototype3 *>(rtiter->second);
0095     assert(raw_tower);
0096
0097     double calibration_const = calib_const_scale;
0098
0099     if (use_chan_calibration)
0100     {
0101       // channel to channel calibration.
0102       const int column = raw_tower->get_column();
0103       const int row = raw_tower->get_row();
0104
0105       assert(column >= 0);
0106       assert(row >= 0);
0107
0108       string calib_const_name(
0109           (boost::format("calib_const_column%d_row%d") % column % row).str());
0110
0111       calibration_const *= _calib_params.get_double_param(calib_const_name);
0112     }
0113
0114     vector<double> vec_signal_samples;
0115     for (int i = 0; i < RawTower_Prototype3::NSAMPLES; i++)
0116     {
0117       vec_signal_samples.push_back(raw_tower->get_signal_samples(i));
0118     }
0119
0120     double peak = NAN;
0121     double peak_sample = NAN;
0122     double pedstal = NAN;
0123
0124     PROTOTYPE3_FEM::SampleFit_PowerLawExp(vec_signal_samples, peak, peak_sample,
0125                                           pedstal, Verbosity());
0126
0127     // store the result - raw_tower
0128     if (std::isnan(raw_tower->get_energy()))
0129     {
0130       // Raw tower was never fit, store the current fit
0131
0132       raw_tower->set_energy(peak);
0133       raw_tower->set_time(peak_sample);
0134     }
0135
0136     // store the result - calib_tower
0137     RawTower_Prototype3 *calib_tower = new RawTower_Prototype3(*raw_tower);
0138     calib_tower->set_energy(peak * calibration_const);
0139     calib_tower->set_time(peak_sample);
0140
0141     for (int i = 0; i < RawTower_Prototype3::NSAMPLES; i++)
0142     {
0143       calib_tower->set_signal_samples(i, (vec_signal_samples[i] - pedstal) *
0144                                              calibration_const);
0145     }
0146
0147     _calib_towers->AddTower(key, calib_tower);
0148
0149   }  //  for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
0150
0151   if (Verbosity())
0152   {
0153     std::cout << Name() << "::" << detector << "::" << __PRETTY_FUNCTION__
0154               << "input sum energy = " << _raw_towers->getTotalEdep()
0155               << ", output sum digitalized value = "
0156               << _calib_towers->getTotalEdep() << std::endl;
0157   }
0158
0159   return Fun4AllReturnCodes::EVENT_OK;
0160 }
0161
0162 //_______________________________________
0163 void CaloCalibration::CreateNodeTree(PHCompositeNode *topNode)
0164 {
0165   PHNodeIterator iter(topNode);
0166
0167   PHCompositeNode *dstNode =
0168       dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0169   if (!dstNode)
0170   {
0171     std::cerr << Name() << "::" << detector << "::" << __PRETTY_FUNCTION__
0172               << "DST Node missing, doing nothing." << std::endl;
0173     throw std::runtime_error(
0174         "Failed to find DST node in RawTowerCalibration::CreateNodes");
0175   }
0176
0177   RawTowerNodeName = "TOWER_" + _raw_tower_node_prefix + "_" + detector;
0178   _raw_towers =
0179       findNode::getClass<RawTowerContainer>(dstNode, RawTowerNodeName.c_str());
0180   if (!_raw_towers)
0181   {
0182     std::cerr << Name() << "::" << detector << "::" << __PRETTY_FUNCTION__
0183               << " " << RawTowerNodeName << " Node missing, doing bail out!"
0184               << std::endl;
0185     throw std::runtime_error("Failed to find " + RawTowerNodeName +
0186                              " node in RawTowerCalibration::CreateNodes");
0187   }
0188
0189   // Create the tower nodes on the tree
0190   PHNodeIterator dstiter(dstNode);
0191   PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(
0192       dstiter.findFirst("PHCompositeNode", detector));
0193   if (!DetNode)
0194   {
0195     DetNode = new PHCompositeNode(detector);
0196     dstNode->addNode(DetNode);
0197   }
0198
0199   // Be careful as a previous calibrator may have been registered for this
0200   // detector
0201   CaliTowerNodeName = "TOWER_" + _calib_tower_node_prefix + "_" + detector;
0202   _calib_towers =
0203       findNode::getClass<RawTowerContainer>(DetNode, CaliTowerNodeName.c_str());
0204   if (!_calib_towers)
0205   {
0206     _calib_towers = new RawTowerContainer(_raw_towers->getCalorimeterID());
0207     PHIODataNode<PHObject> *towerNode = new PHIODataNode<PHObject>(
0208         _calib_towers, CaliTowerNodeName.c_str(), "PHObject");
0209     DetNode->addNode(towerNode);
0210   }
0211
0212   // update the parameters on the node tree
0213   PHCompositeNode *parNode =
0214       dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
0215   assert(parNode);
0216   const string paramnodename = string("Calibration_") + detector;
0217
0218   //   this step is moved to after detector construction
0219   //   save updated persistant copy on node tree
0220   _calib_params.SaveToNodeTree(parNode, paramnodename);
0221 }
0222
0223 //___________________________________
0224 int CaloCalibration::End(PHCompositeNode *topNode)
0225 {
0226   return Fun4AllReturnCodes::EVENT_OK;
0227 }
0228
0229 void CaloCalibration::SetDefaultParameters(PHParameters &param)
0230 {
0231   param.set_int_param("use_chan_calibration", 0);
0232
0233   // additional scale for the calibration constant
0234   // negative pulse -> positive with -1
0235   param.set_double_param("calib_const_scale", -1);
0236 }