packages/Prototype4/CaloCalibration.cc

0001 #include "CaloCalibration.h"
0002
0003 #include "PROTOTYPE4_FEM.h"
0004 #include "RawTower_Prototype4.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 <cstdlib>  // for exit
0027 #include <iostream>
0028 #include <limits>     // for numeric_limits
0029 #include <map>        // for map, _Rb_tree_iter...
0030 #include <stdexcept>  // for runtime_error
0031 #include <string>
0032 #include <utility>  // for pair
0033 #include <vector>   // for vector
0034
0035 using namespace std;
0036
0037 //____________________________________
0038 CaloCalibration::CaloCalibration(const std::string &name)
0039   : SubsysReco(string("CaloCalibration_") + name)
0040   , _calib_towers(nullptr)
0041   , _raw_towers(nullptr)
0042   , detector(name)
0043   , _calib_tower_node_prefix("CALIB")
0044   , _raw_tower_node_prefix("RAW")
0045   , _calib_params(name)
0046   , _fit_type(kPowerLawDoubleExpWithGlobalFitConstraint)
0047 {
0048   SetDefaultParameters(_calib_params);
0049 }
0050
0051 //_____________________________________
0052 int CaloCalibration::InitRun(PHCompositeNode *topNode)
0053 {
0054   CreateNodeTree(topNode);
0055
0056   if (Verbosity())
0057   {
0058     std::cout << Name() << "::" << detector << "::" << __PRETTY_FUNCTION__
0059               << " - print calibration parameters: " << endl;
0060     _calib_params.Print();
0061   }
0062
0063   return Fun4AllReturnCodes::EVENT_OK;
0064 }
0065
0066 //____________________________________
0067 int CaloCalibration::process_event(PHCompositeNode *topNode)
0068 {
0069   if (Verbosity())
0070   {
0071     std::cout << Name() << "::" << detector << "::" << __PRETTY_FUNCTION__
0072               << "Process event entered" << std::endl;
0073   }
0074
0075   map<int, double> parameters_constraints;
0076   if (_fit_type == kPowerLawDoubleExpWithGlobalFitConstraint and
0077       _raw_towers->size() > 1)
0078   {
0079     if (Verbosity())
0080     {
0081       std::cout
0082           << Name() << "::" << detector << "::" << __PRETTY_FUNCTION__
0083           << "Extract global fit parameter for constraining individual fits"
0084           << std::endl;
0085     }
0086
0087     // signal template
0088
0089     vector<double> vec_signal_samples(PROTOTYPE4_FEM::NSAMPLES, 0);
0090
0091     int count = 0;
0092
0093     RawTowerContainer::Range begin_end = _raw_towers->getTowers();
0094     RawTowerContainer::Iterator rtiter;
0095     for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
0096     {
0097       RawTower_Prototype4 *raw_tower =
0098           dynamic_cast<RawTower_Prototype4 *>(rtiter->second);
0099       assert(raw_tower);
0100
0101       //      bool signal_check_pass = true;
0102       //      for (int i = 0; i < RawTower_Prototype4::NSAMPLES; i++)
0103       //      {
0104       //        if (raw_tower->get_signal_samples(i) <= 10 or
0105       //        raw_tower->get_signal_samples(i) >= ((1 << 14) - 10))
0106       //        {
0107       //          signal_check_pass = false;
0108       //          break;
0109       //        }
0110       //      }
0111
0112       //      if (signal_check_pass)
0113       //      {
0114       ++count;
0115
0116       for (int i = 0; i < RawTower_Prototype4::NSAMPLES; i++)
0117       {
0118         if (raw_tower->get_signal_samples(i) <= 10 or
0119             raw_tower->get_signal_samples(i) >= ((1 << 14) - 10))
0120           vec_signal_samples[i] =
0121               numeric_limits<double>::quiet_NaN();  // invalidate this sample
0122         else
0123           vec_signal_samples[i] += raw_tower->get_signal_samples(i);
0124       }
0125       //      }
0126     }
0127
0128     if (count > 0)
0129     {
0130       for (int i = 0; i < RawTower_Prototype4::NSAMPLES; i++)
0131       {
0132         vec_signal_samples[i] /= count;
0133       }
0134
0135       double peak = NAN;
0136       double peak_sample = NAN;
0137       double pedstal = NAN;
0138       map<int, double> parameters_io;
0139
0140       PROTOTYPE4_FEM::SampleFit_PowerLawDoubleExp(vec_signal_samples, peak,
0141                                                   peak_sample, pedstal,
0142                                                   parameters_io, Verbosity());
0143       //    std::map<int, double> &parameters_io,  //! IO for fullset of
0144       //    parameters. If a parameter exist and not an NAN, the fit parameter
0145       //    will be fixed to that value. The order of the parameters are
0146       //    ("Amplitude", "Sample Start", "Power", "Peak Time 1", "Pedestal",
0147       //    "Amplitude ratio", "Peak Time 2")
0148
0149       parameters_constraints[1] = parameters_io[1];
0150       parameters_constraints[2] = parameters_io[2];
0151       parameters_constraints[3] = parameters_io[3];
0152       parameters_constraints[5] = parameters_io[5];
0153       parameters_constraints[6] = parameters_io[6];
0154
0155       //      //special constraint if Peak Time 1 == Peak Time 2
0156       //      if (abs(parameters_constraints[6] - parameters_constraints[3]) <
0157       //      0.1)
0158       //      {
0159       //        const double average_peak_time = (parameters_constraints[6] +
0160       //        parameters_constraints[3]) / 2.;
0161       //
0162       //        std::cout << Name() << "::" << detector << "::" <<
0163       //        __PRETTY_FUNCTION__
0164       //                  << ": two shaping time are too close "
0165       //                  << parameters_constraints[3] << " VS " <<
0166       //                  parameters_constraints[6]
0167       //                  << ". Use average peak time instead: " <<
0168       //                  average_peak_time
0169       //                  << std::endl;
0170       //
0171       //        parameters_constraints[6] = average_peak_time;
0172       //        parameters_constraints[3] = average_peak_time;
0173       //        parameters_constraints[5] = 0;
0174       //      }
0175     }
0176     else
0177     {
0178       std::cout << Name() << "::" << detector << "::" << __PRETTY_FUNCTION__
0179                 << ": Failed to build signal template! Fit each channel "
0180                    "individually instead"
0181                 << std::endl;
0182     }
0183   }
0184
0185   const double calib_const_scale =
0186       _calib_params.get_double_param("calib_const_scale");
0187   const bool use_chan_calibration =
0188       _calib_params.get_int_param("use_chan_calibration") > 0;
0189
0190   RawTowerContainer::Range begin_end = _raw_towers->getTowers();
0191   RawTowerContainer::Iterator rtiter;
0192   for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
0193   {
0194     RawTowerDefs::keytype key = rtiter->first;
0195     RawTower_Prototype4 *raw_tower =
0196         dynamic_cast<RawTower_Prototype4 *>(rtiter->second);
0197     assert(raw_tower);
0198
0199     double calibration_const = calib_const_scale;
0200
0201     if (use_chan_calibration)
0202     {
0203       // channel to channel calibration.
0204       const int column = raw_tower->get_column();
0205       const int row = raw_tower->get_row();
0206
0207       assert(column >= 0);
0208       assert(row >= 0);
0209
0210       string calib_const_name(
0211           (boost::format("calib_const_column%d_row%d") % column % row).str());
0212
0213       calibration_const *= _calib_params.get_double_param(calib_const_name);
0214     }
0215
0216     vector<double> vec_signal_samples;
0217     for (int i = 0; i < RawTower_Prototype4::NSAMPLES; i++)
0218     {
0219       vec_signal_samples.push_back(raw_tower->get_signal_samples(i));
0220     }
0221
0222     double peak = NAN;
0223     double peak_sample = NAN;
0224     double pedstal = NAN;
0225
0226     switch (_fit_type)
0227     {
0228     case kPowerLawExp:
0229       PROTOTYPE4_FEM::SampleFit_PowerLawExp(vec_signal_samples, peak,
0230                                             peak_sample, pedstal, Verbosity());
0231       break;
0232
0233     case kPeakSample:
0234       PROTOTYPE4_FEM::SampleFit_PeakSample(vec_signal_samples, peak,
0235                                            peak_sample, pedstal, Verbosity());
0236       break;
0237
0238     case kPowerLawDoubleExp:
0239     {
0240       map<int, double> parameters_io;
0241
0242       PROTOTYPE4_FEM::SampleFit_PowerLawDoubleExp(vec_signal_samples, peak,
0243                                                   peak_sample, pedstal,
0244                                                   parameters_io, Verbosity());
0245     }
0246     break;
0247
0248     case kPowerLawDoubleExpWithGlobalFitConstraint:
0249     {
0250       map<int, double> parameters_io(parameters_constraints);
0251
0252       PROTOTYPE4_FEM::SampleFit_PowerLawDoubleExp(vec_signal_samples, peak,
0253                                                   peak_sample, pedstal,
0254                                                   parameters_io, Verbosity());
0255     }
0256     break;
0257     default:
0258       cout << __PRETTY_FUNCTION__ << " - FATAL error - unkown fit type "
0259            << _fit_type << endl;
0260       exit(3);
0261       break;
0262     }
0263
0264     // store the result - raw_tower
0265     if (std::isnan(raw_tower->get_energy()))
0266     {
0267       // Raw tower was never fit, store the current fit
0268
0269       raw_tower->set_energy(peak);
0270       raw_tower->set_time(peak_sample);
0271
0272       if (Verbosity())
0273       {
0274         raw_tower->identify();
0275       }
0276     }
0277
0278     // store the result - calib_tower
0279     RawTower_Prototype4 *calib_tower = new RawTower_Prototype4(*raw_tower);
0280     calib_tower->set_energy(peak * calibration_const);
0281     calib_tower->set_time(peak_sample);
0282
0283     for (int i = 0; i < RawTower_Prototype4::NSAMPLES; i++)
0284     {
0285       calib_tower->set_signal_samples(i, (vec_signal_samples[i] - pedstal) *
0286                                              calibration_const);
0287     }
0288
0289     _calib_towers->AddTower(key, calib_tower);
0290
0291   }  //  for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
0292
0293   if (Verbosity())
0294   {
0295     std::cout << Name() << "::" << detector << "::" << __PRETTY_FUNCTION__
0296               << "input sum energy = " << _raw_towers->getTotalEdep()
0297               << ", output sum digitalized value = "
0298               << _calib_towers->getTotalEdep() << std::endl;
0299   }
0300
0301   return Fun4AllReturnCodes::EVENT_OK;
0302 }
0303
0304 //_______________________________________
0305 void CaloCalibration::CreateNodeTree(PHCompositeNode *topNode)
0306 {
0307   PHNodeIterator iter(topNode);
0308
0309   PHCompositeNode *dstNode =
0310       dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0311   if (!dstNode)
0312   {
0313     std::cerr << Name() << "::" << detector << "::" << __PRETTY_FUNCTION__
0314               << "DST Node missing, doing nothing." << std::endl;
0315     throw std::runtime_error(
0316         "Failed to find DST node in RawTowerCalibration::CreateNodes");
0317   }
0318
0319   RawTowerNodeName = "TOWER_" + _raw_tower_node_prefix + "_" + detector;
0320   _raw_towers =
0321       findNode::getClass<RawTowerContainer>(dstNode, RawTowerNodeName.c_str());
0322   if (!_raw_towers)
0323   {
0324     std::cerr << Name() << "::" << detector << "::" << __PRETTY_FUNCTION__
0325               << " " << RawTowerNodeName << " Node missing, doing bail out!"
0326               << std::endl;
0327     throw std::runtime_error("Failed to find " + RawTowerNodeName +
0328                              " node in RawTowerCalibration::CreateNodes");
0329   }
0330
0331   // Create the tower nodes on the tree
0332   PHNodeIterator dstiter(dstNode);
0333   PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(
0334       dstiter.findFirst("PHCompositeNode", detector));
0335   if (!DetNode)
0336   {
0337     DetNode = new PHCompositeNode(detector);
0338     dstNode->addNode(DetNode);
0339   }
0340
0341   // Be careful as a previous calibrator may have been registered for this
0342   // detector
0343   CaliTowerNodeName = "TOWER_" + _calib_tower_node_prefix + "_" + detector;
0344   _calib_towers =
0345       findNode::getClass<RawTowerContainer>(DetNode, CaliTowerNodeName.c_str());
0346   if (!_calib_towers)
0347   {
0348     _calib_towers = new RawTowerContainer(_raw_towers->getCalorimeterID());
0349     PHIODataNode<PHObject> *towerNode = new PHIODataNode<PHObject>(
0350         _calib_towers, CaliTowerNodeName.c_str(), "PHObject");
0351     DetNode->addNode(towerNode);
0352   }
0353
0354   // update the parameters on the node tree
0355   PHCompositeNode *parNode =
0356       dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
0357   assert(parNode);
0358   const string paramnodename = string("Calibration_") + detector;
0359
0360   //   this step is moved to after detector construction
0361   //   save updated persistant copy on node tree
0362   _calib_params.SaveToNodeTree(parNode, paramnodename);
0363 }
0364
0365 void CaloCalibration::SetDefaultParameters(PHParameters &param)
0366 {
0367   param.set_int_param("use_chan_calibration", 0);
0368
0369   // additional scale for the calibration constant
0370   // negative pulse -> positive with -1
0371   param.set_double_param("calib_const_scale", 1);
0372 }