packages/Prototype2/CaloCalibration.cc

0001 #include "CaloCalibration.h"
0002
0003 #include "PROTOTYPE2_FEM.h"
0004 #include "RawTower_Prototype2.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>
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 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::InitRun(PHCompositeNode *topNode)
0055 {
0056   CreateNodeTree(topNode);
0057
0058   if (Verbosity())
0059   {
0060     cout << Name() << "::" << detector << "::" << __PRETTY_FUNCTION__
0061          << " - print calibration parameters: " << endl;
0062     _calib_params.Print();
0063   }
0064
0065   return Fun4AllReturnCodes::EVENT_OK;
0066 }
0067
0068 //____________________________________
0069 int CaloCalibration::process_event(PHCompositeNode */*topNode*/)
0070 {
0071   if (Verbosity())
0072   {
0073     cout << Name() << "::" << detector << "::" << __PRETTY_FUNCTION__
0074          << "Process event entered" << endl;
0075   }
0076
0077   const double calib_const_scale =
0078       _calib_params.get_double_param("calib_const_scale");
0079   const bool use_chan_calibration =
0080       _calib_params.get_int_param("use_chan_calibration") > 0;
0081
0082   RawTowerContainer::Range begin_end = _raw_towers->getTowers();
0083   RawTowerContainer::Iterator rtiter;
0084   for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
0085   {
0086     RawTowerDefs::keytype key = rtiter->first;
0087     RawTower_Prototype2 *raw_tower =
0088         dynamic_cast<RawTower_Prototype2 *>(rtiter->second);
0089     assert(raw_tower);
0090
0091     double calibration_const = calib_const_scale;
0092
0093     if (use_chan_calibration)
0094     {
0095       // channel to channel calibration.
0096       const int column = raw_tower->get_column();
0097       const int row = raw_tower->get_row();
0098
0099       assert(column >= 0);
0100       assert(row >= 0);
0101       string calib_const_name(
0102           (boost::format("calib_const_column%d_row%d") % column % row).str());
0103
0104       calibration_const *= _calib_params.get_double_param(calib_const_name);
0105     }
0106
0107     vector<double> vec_signal_samples;
0108     for (int i = 0; i < RawTower_Prototype2::NSAMPLES; i++)
0109     {
0110       vec_signal_samples.push_back(raw_tower->get_signal_samples(i));
0111     }
0112
0113     double peak = NAN;
0114     double peak_sample = NAN;
0115     double pedstal = NAN;
0116
0117     PROTOTYPE2_FEM::SampleFit_PowerLawExp(vec_signal_samples, peak, peak_sample,
0118                                           pedstal, Verbosity());
0119
0120     // store the result - raw_tower
0121     if (std::isnan(raw_tower->get_energy()))
0122     {
0123       // Raw tower was never fit, store the current fit
0124
0125       raw_tower->set_energy(peak);
0126       raw_tower->set_time(peak_sample);
0127     }
0128
0129     // store the result - calib_tower
0130     RawTower_Prototype2 *calib_tower = new RawTower_Prototype2(*raw_tower);
0131     calib_tower->set_energy(peak * calibration_const);
0132     calib_tower->set_time(peak_sample);
0133
0134     for (int i = 0; i < RawTower_Prototype2::NSAMPLES; i++)
0135     {
0136       calib_tower->set_signal_samples(i, (vec_signal_samples[i] - pedstal) *
0137                                              calibration_const);
0138     }
0139
0140     _calib_towers->AddTower(key, calib_tower);
0141
0142   }  //  for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
0143
0144   if (Verbosity())
0145   {
0146     cout << Name() << "::" << detector << "::" << __PRETTY_FUNCTION__
0147          << "input sum energy = " << _raw_towers->getTotalEdep()
0148          << ", output sum digitalized value = " << _calib_towers->getTotalEdep()
0149          << endl;
0150   }
0151
0152   return Fun4AllReturnCodes::EVENT_OK;
0153 }
0154
0155 //_______________________________________
0156 void CaloCalibration::CreateNodeTree(PHCompositeNode *topNode)
0157 {
0158   PHNodeIterator iter(topNode);
0159
0160   PHCompositeNode *dstNode =
0161       dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0162   if (!dstNode)
0163   {
0164     cout << Name() << "::" << detector << "::" << __PRETTY_FUNCTION__
0165          << "DST Node missing, doing nothing." << endl;
0166     throw runtime_error(
0167         "Failed to find DST node in RawTowerCalibration::CreateNodes");
0168   }
0169
0170   RawTowerNodeName = "TOWER_" + _raw_tower_node_prefix + "_" + detector;
0171   _raw_towers =
0172       findNode::getClass<RawTowerContainer>(dstNode, RawTowerNodeName.c_str());
0173   if (!_raw_towers)
0174   {
0175     cout << Name() << "::" << detector << "::" << __PRETTY_FUNCTION__ << " "
0176          << RawTowerNodeName << " Node missing, doing bail out!" << endl;
0177     throw runtime_error("Failed to find " + RawTowerNodeName +
0178                         " node in RawTowerCalibration::CreateNodes");
0179   }
0180
0181   // Create the tower nodes on the tree
0182   PHNodeIterator dstiter(dstNode);
0183   PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(
0184       dstiter.findFirst("PHCompositeNode", detector));
0185   if (!DetNode)
0186   {
0187     DetNode = new PHCompositeNode(detector);
0188     dstNode->addNode(DetNode);
0189   }
0190
0191   // Be careful as a previous calibrator may have been registered for this
0192   // detector
0193   CaliTowerNodeName = "TOWER_" + _calib_tower_node_prefix + "_" + detector;
0194   _calib_towers =
0195       findNode::getClass<RawTowerContainer>(DetNode, CaliTowerNodeName.c_str());
0196   if (!_calib_towers)
0197   {
0198     _calib_towers = new RawTowerContainer(_raw_towers->getCalorimeterID());
0199     PHIODataNode<PHObject> *towerNode = new PHIODataNode<PHObject>(
0200         _calib_towers, CaliTowerNodeName.c_str(), "PHObject");
0201     DetNode->addNode(towerNode);
0202   }
0203
0204   // update the parameters on the node tree
0205   PHCompositeNode *parNode =
0206       dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
0207   assert(parNode);
0208   const string paramnodename = string("Calibration_") + detector;
0209
0210   //   this step is moved to after detector construction
0211   //   save updated persistant copy on node tree
0212   _calib_params.SaveToNodeTree(parNode, paramnodename);
0213 }
0214
0215 void CaloCalibration::SetDefaultParameters(PHParameters &param)
0216 {
0217   param.set_int_param("use_chan_calibration", 0);
0218
0219   // additional scale for the calibration constant
0220   // negative pulse -> positive with -1
0221   param.set_double_param("calib_const_scale", -1);
0222 }