EMCTemplateFitting/src/CaloTemplateFit.cc

0001 #include "CaloTemplateFit.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 #include <TProfile.h>
0035
0036 #include <pthread.h>
0037 #include <TFile.h>
0038 #include <ROOT/TThreadExecutor.hxx>
0039 #include <ROOT/TThreadedObject.hxx>
0040 #include <iostream>
0041 #include "TThread.h"
0042 #include "TMutex.h"
0043 #include "TMath.h"
0044
0045 #include "Math/WrappedTF1.h"
0046 #include "Math/WrappedMultiTF1.h"
0047 #include "Fit/BinData.h"
0048 #include "Fit/UnBinData.h"
0049 #include "HFitInterface.h"
0050 #include "Fit/Fitter.h"
0051 #include "Fit/Chi2FCN.h"
0052 #include <pthread.h>
0053 #include "Math/Functor.h"
0054 #include <TH1F.h>
0055 using namespace std;
0056
0057
0058
0059
0060
0061 double CaloTemplateFit::template_function(double *x, double *par)
0062 {
0063   Double_t v1 = par[0]*h_template->Interpolate(x[0]-par[1])+par[2];
0064   // Double_t v1 = par[0]*TMath::Power(x[0],par[1])+par[2];
0065   return v1;
0066 }
0067
0068
0069
0070
0071 std::vector<std::vector<float>> CaloTemplateFit::calo_processing_perchnl(std::vector<std::vector<float>> chnlvector)
0072 {
0073   ROOT::TThreadExecutor t(_nthreads);
0074   auto func = [&](std::vector<float> &v) {
0075     int size1 = v.size()-1;
0076     auto h = new TH1F(Form("h_%d",(int)round(v.at(size1))),"",size1,0,size1);
0077     float maxheight = 0;
0078     int maxbin = 0;
0079     for (int i = 0; i < size1;i++)
0080       {
0081     h->SetBinContent(i+1,v.at(i));
0082     if (v.at(i)>maxheight)
0083       {
0084         maxheight = v.at(i);
0085         maxbin = i;
0086       }
0087       }
0088     float pedestal = 1500;
0089     if (maxbin > 4)
0090       {
0091     pedestal=  0.5* ( v.at(maxbin - 4) + v.at(maxbin-5));
0092       }
0093     else if (maxbin > 3)
0094       {
0095     pedestal=( v.at(maxbin - 4));
0096       }
0097     else
0098       {
0099     pedestal = 0.5* ( v.at(size1-3) + v.at(size1-2));
0100       }
0101     auto f = new TF1(Form("f_%d",(int)round(v.at(size1))),template_function,0,31,3);
0102     ROOT::Math::WrappedMultiTF1 * fitFunction = new ROOT::Math::WrappedMultiTF1(*f, 3 );
0103     ROOT::Fit::BinData data(v.size()-1,1);
0104     ROOT::Fit::FillData(data,h);
0105     ROOT::Fit::Chi2Function *EPChi2 = new ROOT::Fit::Chi2Function(data, *fitFunction);
0106     ROOT::Fit::Fitter *fitter = new ROOT::Fit::Fitter();
0107     fitter->Config().MinimizerOptions().SetMinimizerType("GSLMultiFit");
0108     double params[] = {static_cast<double>(maxheight),static_cast<double>(maxbin-5),static_cast<double>(pedestal)};
0109     fitter->Config().SetParamsSettings(3,params);
0110     fitter->FitFCN(*EPChi2,0,data.Size(),true);
0111     for (int i =0;i<3;i++)
0112       {
0113     v.push_back(f->GetParameter(i));
0114       }
0115     h->Delete();
0116     f->Delete();
0117     delete fitFunction;
0118     delete fitter;
0119     delete EPChi2;
0120     // v.clear();
0121   };
0122
0123   t.Foreach(func, chnlvector);
0124
0125   int size3 = chnlvector.size();
0126   std::vector<std::vector<float>> fit_params;
0127   std::vector<float> fit_params_tmp;
0128   for (int i = 0; i < size3;i++)
0129     {
0130       std::vector<float> tv = chnlvector.at(i);
0131       int size2 = tv.size();
0132       for (int q = 3; q > 0;q--)
0133     {
0134       fit_params_tmp.push_back(tv.at(size2-q));
0135     }
0136       fit_params.push_back(fit_params_tmp);
0137       fit_params_tmp.clear();
0138     }
0139
0140
0141   chnlvector.clear();
0142   return fit_params;
0143   fit_params.clear();
0144 }
0145
0146
0147
0148
0149 std::vector<float> CaloTemplateFit::calo_processing_singlethread(std::vector<float> chnlvector)
0150 {
0151   int size1 = chnlvector.size();
0152   float maxheight = 0;
0153   int maxbin = 0;
0154   for (int i = 0; i < size1;i++)
0155     {
0156       h_data->SetBinContent(i+1,chnlvector.at(i));
0157       if (chnlvector.at(i)>maxheight)
0158     {
0159       maxheight = chnlvector.at(i);
0160       maxbin = i;
0161     }
0162     }
0163   float pedestal = 1500;
0164   if (maxbin > 4)
0165     {
0166       pedestal=  0.5* ( chnlvector.at(maxbin - 4) + chnlvector.at(maxbin-5));
0167     }
0168   else if (maxbin > 3)
0169     {
0170       pedestal=( chnlvector.at(maxbin - 4));
0171     }
0172   else
0173     {
0174       pedestal = 0.5* ( chnlvector.at(size1-3) + chnlvector.at(size1-2));
0175     }
0176   ROOT::Math::WrappedMultiTF1 fitFunction(*f_fit, 3 );
0177   ROOT::Fit::BinData data(chnlvector.size()-1,1);
0178   ROOT::Fit::FillData(data,h_data);
0179   ROOT::Fit::Chi2Function EPChi2(data, fitFunction);
0180   ROOT::Fit::Fitter fitter;
0181   fitter.Config().MinimizerOptions().SetMinimizerType("GSLMultiFit");
0182   double params[] = {static_cast<double>(maxheight),static_cast<double>(maxbin-5),static_cast<double>(pedestal)};
0183   fitter.Config().SetParamsSettings(3,params);
0184   fitter.FitFCN(EPChi2,0,data.Size(),true);
0185
0186   std::vector<float> fit_params;
0187   for (int q = 0; q < 3;q++)
0188     {
0189       fit_params.push_back(f_fit->GetParameter(q));
0190     }
0191   return fit_params;
0192
0193
0194 }
0195
0196 //TProfile for the template
0197 TProfile* CaloTemplateFit::h_template = nullptr;
0198
0199
0200
0201 //____________________________________
0202 CaloTemplateFit::CaloTemplateFit(const std::string &name)
0203   : SubsysReco(string("CaloCalibration_") + name)
0204   , _calib_towers(nullptr)
0205   , _raw_towers(nullptr)
0206   , detector(name)
0207   , _calib_tower_node_prefix("CALIB")
0208   , _raw_tower_node_prefix("RAW")
0209   , _nthreads(1)
0210   , _nsamples(0)
0211   , template_input_file("/gpfs/mnt/gpfs02/sphenix/user/trinn/fitting_algorithm_playing/prdfcode/prototype/offline/packages/Prototype4/templates.root")
0212  , _calib_params(name)
0213   , _fit_type(kPowerLawDoubleExpWithGlobalFitConstraint)
0214 {
0215   SetDefaultParameters(_calib_params);
0216 }
0217
0218 //_____________________________________
0219 int CaloTemplateFit::InitRun(PHCompositeNode *topNode)
0220 {
0221   CreateNodeTree(topNode);
0222
0223   if (Verbosity())
0224   {
0225     std::cout << Name() << "::" << detector << "::" << __PRETTY_FUNCTION__
0226               << " - print calibration parameters: " << endl;
0227     _calib_params.Print();
0228   }
0229
0230
0231   TFile* fin = TFile::Open(template_input_file.c_str());
0232   assert(fin);
0233   assert(fin->IsOpen());
0234   h_template = static_cast<TProfile*>(fin->Get("hp_electrons_fine_emcal_36_8GeV"));
0235   h_template->SetDirectory(0);
0236   fin->Close();
0237   delete fin;
0238
0239
0240   h_data = new TH1F("h_data","",32,0,32);
0241   f_fit = new TF1("f_fit",template_function,0,31,3);
0242
0243
0244   ROOT::EnableThreadSafety();
0245
0246   return Fun4AllReturnCodes::EVENT_OK;
0247 }
0248
0249 //____________________________________
0250 int CaloTemplateFit::process_event(PHCompositeNode *topNode)
0251 {
0252   if (Verbosity())
0253   {
0254     std::cout << Name() << "::" << detector << "::" << __PRETTY_FUNCTION__
0255               << "Process event entered" << std::endl;
0256   }
0257
0258   map<int, double> parameters_constraints;
0259
0260   if ( _raw_towers->size() > 1)
0261     {
0262       if (Verbosity())
0263     {
0264       std::cout
0265         << Name() << "::" << detector << "::" << __PRETTY_FUNCTION__
0266         << "Extract global fit parameter for constraining individual fits"
0267         << std::endl;
0268     }
0269       // signal template
0270       vector<float> waveforms;
0271       vector<vector<float>> waveforms2;
0272       int count = 0;
0273       RawTowerContainer::Range begin_end = _raw_towers->getTowers();
0274       RawTowerContainer::Iterator rtiter;
0275       for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
0276     {
0277       RawTower_Prototype4 *raw_tower =
0278         dynamic_cast<RawTower_Prototype4 *>(rtiter->second);
0279       assert(raw_tower);
0280       ++count;
0281
0282       for (int i = 0; i < RawTower_Prototype4::NSAMPLES; i++)
0283         {
0284           if (i >= _nsamples && _nsamples != 0){continue;}
0285           waveforms.push_back(raw_tower->get_signal_samples(i));
0286         }
0287
0288       // std::cout << waveforms.size() << std::endl;
0289       waveforms.push_back(count);
0290       waveforms2.push_back(waveforms);
0291       waveforms.clear();
0292     }
0293        std::vector< std::vector<float>> pulse_quantification = calo_processing_perchnl(waveforms2);
0294
0295        int towernumber = 0;
0296
0297        for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
0298      {
0299        // RawTowerDefs::keytype key = rtiter->first;
0300        RawTower_Prototype4 *raw_tower =
0301          dynamic_cast<RawTower_Prototype4 *>(rtiter->second);
0302        assert(raw_tower);
0303        // store the result - raw_tower
0304        if (std::isnan(raw_tower->get_energy()))
0305          {
0306            // Raw tower was never fit, store the current fit
0307            std::vector<float> values = pulse_quantification.at(towernumber);
0308            raw_tower->set_energy(values.at(0));
0309            raw_tower->set_time(values.at(1));
0310            values.clear();
0311          }
0312        towernumber++;
0313      }
0314        pulse_quantification.clear();
0315        waveforms2.clear();
0316        // std::cout <<waveforms.size() << " , " <<  waveforms2.size() << " , "<< pulse_quantification.size() << std::endl;
0317    }
0318
0319
0320     /*
0321     if (count > 0)
0322     {
0323       for (int i = 0; i < RawTower_Prototype4::NSAMPLES; i++)
0324       {
0325         vec_signal_samples[i] /= count;
0326       }
0327
0328       double peak = NAN;
0329       double peak_sample = NAN;
0330       double pedstal = NAN;
0331       map<int, double> parameters_io;
0332
0333       PROTOTYPE4_FEM::SampleFit_PowerLawDoubleExp(vec_signal_samples, peak,
0334                                                   peak_sample, pedstal,
0335                                                   parameters_io, Verbosity());
0336       //    std::map<int, double> &parameters_io,  //! IO for fullset of
0337       //    parameters. If a parameter exist and not an NAN, the fit parameter
0338       //    will be fixed to that value. The order of the parameters are
0339       //    ("Amplitude", "Sample Start", "Power", "Peak Time 1", "Pedestal",
0340       //    "Amplitude ratio", "Peak Time 2")
0341
0342       parameters_constraints[1] = parameters_io[1];
0343       parameters_constraints[2] = parameters_io[2];
0344       parameters_constraints[3] = parameters_io[3];
0345       parameters_constraints[5] = parameters_io[5];
0346       parameters_constraints[6] = parameters_io[6];
0347
0348       //      //special constraint if Peak Time 1 == Peak Time 2
0349       //      if (abs(parameters_constraints[6] - parameters_constraints[3]) <
0350       //      0.1)
0351       //      {
0352       //        const double average_peak_time = (parameters_constraints[6] +
0353       //        parameters_constraints[3]) / 2.;
0354       //
0355       //        std::cout << Name() << "::" << detector << "::" <<
0356       //        __PRETTY_FUNCTION__
0357       //                  << ": two shaping time are too close "
0358       //                  << parameters_constraints[3] << " VS " <<
0359       //                  parameters_constraints[6]
0360       //                  << ". Use average peak time instead: " <<
0361       //                  average_peak_time
0362       //                  << std::endl;
0363       //
0364       //        parameters_constraints[6] = average_peak_time;
0365       //        parameters_constraints[3] = average_peak_time;
0366       //        parameters_constraints[5] = 0;
0367       //      }
0368     }
0369     else
0370     {
0371       std::cout << Name() << "::" << detector << "::" << __PRETTY_FUNCTION__
0372                 << ": Failed to build signal template! Fit each channel "
0373                    "individually instead"
0374                 << std::endl;
0375     }
0376   }
0377     */
0378   /*
0379   const double calib_const_scale =
0380       _calib_params.get_double_param("calib_const_scale");
0381   const bool use_chan_calibration =
0382       _calib_params.get_int_param("use_chan_calibration") > 0;
0383
0384   RawTowerContainer::Range begin_end = _raw_towers->getTowers();
0385   RawTowerContainer::Iterator rtiter;
0386   for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
0387   {
0388     RawTowerDefs::keytype key = rtiter->first;
0389     RawTower_Prototype4 *raw_tower =
0390         dynamic_cast<RawTower_Prototype4 *>(rtiter->second);
0391     assert(raw_tower);
0392
0393     double calibration_const = calib_const_scale;
0394
0395     if (use_chan_calibration)
0396     {
0397       // channel to channel calibration.
0398       const int column = raw_tower->get_column();
0399       const int row = raw_tower->get_row();
0400
0401       assert(column >= 0);
0402       assert(row >= 0);
0403
0404       string calib_const_name(
0405           (boost::format("calib_const_column%d_row%d") % column % row).str());
0406
0407       calibration_const *= _calib_params.get_double_param(calib_const_name);
0408     }
0409
0410     vector<double> vec_signal_samples;
0411     for (int i = 0; i < RawTower_Prototype4::NSAMPLES; i++)
0412     {
0413       vec_signal_samples.push_back(raw_tower->get_signal_samples(i));
0414     }
0415
0416     double peak = NAN;
0417     double peak_sample = NAN;
0418     double pedstal = NAN;
0419
0420     switch (_fit_type)
0421     {
0422     case kPowerLawExp:
0423       PROTOTYPE4_FEM::SampleFit_PowerLawExp(vec_signal_samples, peak,
0424                                             peak_sample, pedstal, Verbosity());
0425       break;
0426
0427     case kPeakSample:
0428       PROTOTYPE4_FEM::SampleFit_PeakSample(vec_signal_samples, peak,
0429                                            peak_sample, pedstal, Verbosity());
0430       break;
0431
0432     case kPowerLawDoubleExp:
0433     {
0434       map<int, double> parameters_io;
0435
0436       PROTOTYPE4_FEM::SampleFit_PowerLawDoubleExp(vec_signal_samples, peak,
0437                                                   peak_sample, pedstal,
0438                                                   parameters_io, Verbosity());
0439     }
0440     break;
0441
0442     case kPowerLawDoubleExpWithGlobalFitConstraint:
0443     {
0444       map<int, double> parameters_io(parameters_constraints);
0445
0446       PROTOTYPE4_FEM::SampleFit_PowerLawDoubleExp(vec_signal_samples, peak,
0447                                                   peak_sample, pedstal,
0448                                                   parameters_io, Verbosity());
0449     }
0450     break;
0451     default:
0452       cout << __PRETTY_FUNCTION__ << " - FATAL error - unkown fit type "
0453            << _fit_type << endl;
0454       exit(3);
0455       break;
0456     }
0457
0458     // store the result - raw_tower
0459     if (std::isnan(raw_tower->get_energy()))
0460     {
0461       // Raw tower was never fit, store the current fit
0462
0463       raw_tower->set_energy(peak);
0464       raw_tower->set_time(peak_sample);
0465
0466       if (Verbosity())
0467       {
0468         raw_tower->identify();
0469       }
0470     }
0471
0472     // store the result - calib_tower
0473     RawTower_Prototype4 *calib_tower = new RawTower_Prototype4(*raw_tower);
0474     calib_tower->set_energy(peak * calibration_const);
0475     calib_tower->set_time(peak_sample);
0476
0477     for (int i = 0; i < RawTower_Prototype4::NSAMPLES; i++)
0478     {
0479       calib_tower->set_signal_samples(i, (vec_signal_samples[i] - pedstal) *
0480                                              calibration_const);
0481     }
0482
0483     _calib_towers->AddTower(key, calib_tower);
0484
0485   }  //  for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
0486
0487   if (Verbosity())
0488   {
0489     std::cout << Name() << "::" << detector << "::" << __PRETTY_FUNCTION__
0490               << "input sum energy = " << _raw_towers->getTotalEdep()
0491               << ", output sum digitalized value = "
0492               << _calib_towers->getTotalEdep() << std::endl;
0493   }
0494 */
0495
0496
0497   return Fun4AllReturnCodes::EVENT_OK;
0498 }
0499
0500 //_______________________________________
0501 void CaloTemplateFit::CreateNodeTree(PHCompositeNode *topNode)
0502 {
0503   PHNodeIterator iter(topNode);
0504
0505   PHCompositeNode *dstNode =
0506       dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0507   if (!dstNode)
0508   {
0509     std::cerr << Name() << "::" << detector << "::" << __PRETTY_FUNCTION__
0510               << "DST Node missing, doing nothing." << std::endl;
0511     throw std::runtime_error(
0512         "Failed to find DST node in RawTowerCalibration::CreateNodes");
0513   }
0514
0515   RawTowerNodeName = "TOWER_" + _raw_tower_node_prefix + "_" + detector;
0516   _raw_towers =
0517       findNode::getClass<RawTowerContainer>(dstNode, RawTowerNodeName.c_str());
0518   if (!_raw_towers)
0519   {
0520     std::cerr << Name() << "::" << detector << "::" << __PRETTY_FUNCTION__
0521               << " " << RawTowerNodeName << " Node missing, doing bail out!"
0522               << std::endl;
0523     throw std::runtime_error("Failed to find " + RawTowerNodeName +
0524                              " node in RawTowerCalibration::CreateNodes");
0525   }
0526
0527   // Create the tower nodes on the tree
0528   PHNodeIterator dstiter(dstNode);
0529   PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(
0530       dstiter.findFirst("PHCompositeNode", detector));
0531   if (!DetNode)
0532   {
0533     DetNode = new PHCompositeNode(detector);
0534     dstNode->addNode(DetNode);
0535   }
0536
0537   // Be careful as a previous calibrator may have been registered for this
0538   // detector
0539   CaliTowerNodeName = "TOWER_" + _calib_tower_node_prefix + "_" + detector;
0540   _calib_towers =
0541       findNode::getClass<RawTowerContainer>(DetNode, CaliTowerNodeName.c_str());
0542   if (!_calib_towers)
0543   {
0544     _calib_towers = new RawTowerContainer(_raw_towers->getCalorimeterID());
0545     PHIODataNode<PHObject> *towerNode = new PHIODataNode<PHObject>(
0546         _calib_towers, CaliTowerNodeName.c_str(), "PHObject");
0547     DetNode->addNode(towerNode);
0548   }
0549
0550   // update the parameters on the node tree
0551   PHCompositeNode *parNode =
0552       dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
0553   assert(parNode);
0554   const string paramnodename = string("Calibration_") + detector;
0555
0556   //   this step is moved to after detector construction
0557   //   save updated persistant copy on node tree
0558   _calib_params.SaveToNodeTree(parNode, paramnodename);
0559 }
0560
0561 void CaloTemplateFit::SetDefaultParameters(PHParameters &param)
0562 {
0563   param.set_int_param("use_chan_calibration", 0);
0564
0565   // additional scale for the calibration constant
0566   // negative pulse -> positive with -1
0567   param.set_double_param("calib_const_scale", 1);
0568 }