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