Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CaloTemplateFit.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file CaloTemplateFit.cc
1 #include "CaloTemplateFit.h"
2 
3 #include "PROTOTYPE4_FEM.h"
4 #include "RawTower_Prototype4.h"
5 
6 #include <calobase/RawTower.h> // for RawTower
7 #include <calobase/RawTowerContainer.h>
8 #include <calobase/RawTowerDefs.h> // for keytype
9 
10 #include <phparameter/PHParameters.h> // for PHParameters
11 
12 #include <fun4all/Fun4AllReturnCodes.h>
13 #include <fun4all/SubsysReco.h> // for SubsysReco
14 
15 #include <phool/PHCompositeNode.h>
16 #include <phool/PHIODataNode.h> // for PHIODataNode
17 #include <phool/PHNode.h> // for PHNode
18 #include <phool/PHNodeIterator.h> // for PHNodeIterator
19 #include <phool/PHObject.h> // for PHObject
20 #include <phool/getClass.h>
21 
22 #include <boost/format.hpp>
23 
24 #include <cassert>
25 #include <cmath> // for NAN, isnan
26 #include <cstdlib> // for exit
27 #include <iostream>
28 #include <limits> // for numeric_limits
29 #include <map> // for map, _Rb_tree_iter...
30 #include <stdexcept> // for runtime_error
31 #include <string>
32 #include <utility> // for pair
33 #include <vector> // for vector
34 #include <TProfile.h>
35 
36 #include <pthread.h>
37 #include <TFile.h>
38 #include <ROOT/TThreadExecutor.hxx>
39 #include <ROOT/TThreadedObject.hxx>
40 #include <iostream>
41 #include "TThread.h"
42 #include "TMutex.h"
43 #include "TMath.h"
44 
45 #include "Math/WrappedTF1.h"
46 #include "Math/WrappedMultiTF1.h"
47 #include "Fit/BinData.h"
48 #include "Fit/UnBinData.h"
49 #include "HFitInterface.h"
50 #include "Fit/Fitter.h"
51 #include "Fit/Chi2FCN.h"
52 #include <pthread.h>
53 #include "Math/Functor.h"
54 #include <TH1F.h>
55 using namespace std;
56 
57 
58 
59 
60 
61 double CaloTemplateFit::template_function(double *x, double *par)
62 {
63  Double_t v1 = par[0]*h_template->Interpolate(x[0]-par[1])+par[2];
64  // Double_t v1 = par[0]*TMath::Power(x[0],par[1])+par[2];
65  return v1;
66 }
67 
68 
69 
70 
71 std::vector<std::vector<float>> CaloTemplateFit::calo_processing_perchnl(std::vector<std::vector<float>> chnlvector)
72 {
73  ROOT::TThreadExecutor t(_nthreads);
74  auto func = [&](std::vector<float> &v) {
75  int size1 = v.size()-1;
76  auto h = new TH1F(Form("h_%d",(int)round(v.at(size1))),"",size1,0,size1);
77  float maxheight = 0;
78  int maxbin = 0;
79  for (int i = 0; i < size1;i++)
80  {
81  h->SetBinContent(i+1,v.at(i));
82  if (v.at(i)>maxheight)
83  {
84  maxheight = v.at(i);
85  maxbin = i;
86  }
87  }
88  float pedestal = 1500;
89  if (maxbin > 4)
90  {
91  pedestal= 0.5* ( v.at(maxbin - 4) + v.at(maxbin-5));
92  }
93  else if (maxbin > 3)
94  {
95  pedestal=( v.at(maxbin - 4));
96  }
97  else
98  {
99  pedestal = 0.5* ( v.at(size1-3) + v.at(size1-2));
100  }
101  auto f = new TF1(Form("f_%d",(int)round(v.at(size1))),template_function,0,31,3);
102  ROOT::Math::WrappedMultiTF1 * fitFunction = new ROOT::Math::WrappedMultiTF1(*f, 3 );
103  ROOT::Fit::BinData data(v.size()-1,1);
104  ROOT::Fit::FillData(data,h);
105  ROOT::Fit::Chi2Function *EPChi2 = new ROOT::Fit::Chi2Function(data, *fitFunction);
106  ROOT::Fit::Fitter *fitter = new ROOT::Fit::Fitter();
107  fitter->Config().MinimizerOptions().SetMinimizerType("GSLMultiFit");
108  double params[] = {static_cast<double>(maxheight),static_cast<double>(maxbin-5),static_cast<double>(pedestal)};
109  fitter->Config().SetParamsSettings(3,params);
110  fitter->FitFCN(*EPChi2,0,data.Size(),true);
111  for (int i =0;i<3;i++)
112  {
113  v.push_back(f->GetParameter(i));
114  }
115  h->Delete();
116  f->Delete();
117  delete fitFunction;
118  delete fitter;
119  delete EPChi2;
120  // v.clear();
121  };
122 
123  t.Foreach(func, chnlvector);
124 
125  int size3 = chnlvector.size();
126  std::vector<std::vector<float>> fit_params;
127  std::vector<float> fit_params_tmp;
128  for (int i = 0; i < size3;i++)
129  {
130  std::vector<float> tv = chnlvector.at(i);
131  int size2 = tv.size();
132  for (int q = 3; q > 0;q--)
133  {
134  fit_params_tmp.push_back(tv.at(size2-q));
135  }
136  fit_params.push_back(fit_params_tmp);
137  fit_params_tmp.clear();
138  }
139 
140 
141  chnlvector.clear();
142  return fit_params;
143  fit_params.clear();
144 }
145 
146 
147 
148 
149 std::vector<float> CaloTemplateFit::calo_processing_singlethread(std::vector<float> chnlvector)
150 {
151  int size1 = chnlvector.size();
152  float maxheight = 0;
153  int maxbin = 0;
154  for (int i = 0; i < size1;i++)
155  {
156  h_data->SetBinContent(i+1,chnlvector.at(i));
157  if (chnlvector.at(i)>maxheight)
158  {
159  maxheight = chnlvector.at(i);
160  maxbin = i;
161  }
162  }
163  float pedestal = 1500;
164  if (maxbin > 4)
165  {
166  pedestal= 0.5* ( chnlvector.at(maxbin - 4) + chnlvector.at(maxbin-5));
167  }
168  else if (maxbin > 3)
169  {
170  pedestal=( chnlvector.at(maxbin - 4));
171  }
172  else
173  {
174  pedestal = 0.5* ( chnlvector.at(size1-3) + chnlvector.at(size1-2));
175  }
176  ROOT::Math::WrappedMultiTF1 fitFunction(*f_fit, 3 );
177  ROOT::Fit::BinData data(chnlvector.size()-1,1);
178  ROOT::Fit::FillData(data,h_data);
179  ROOT::Fit::Chi2Function EPChi2(data, fitFunction);
180  ROOT::Fit::Fitter fitter;
181  fitter.Config().MinimizerOptions().SetMinimizerType("GSLMultiFit");
182  double params[] = {static_cast<double>(maxheight),static_cast<double>(maxbin-5),static_cast<double>(pedestal)};
183  fitter.Config().SetParamsSettings(3,params);
184  fitter.FitFCN(EPChi2,0,data.Size(),true);
185 
186  std::vector<float> fit_params;
187  for (int q = 0; q < 3;q++)
188  {
189  fit_params.push_back(f_fit->GetParameter(q));
190  }
191  return fit_params;
192 
193 
194 }
195 
196 //TProfile for the template
197 TProfile* CaloTemplateFit::h_template = nullptr;
198 
199 
200 
201 //____________________________________
202 CaloTemplateFit::CaloTemplateFit(const std::string &name)
203  : SubsysReco(string("CaloCalibration_") + name)
204  , _calib_towers(nullptr)
205  , _raw_towers(nullptr)
206  , detector(name)
207  , _calib_tower_node_prefix("CALIB")
208  , _raw_tower_node_prefix("RAW")
209  , _nthreads(1)
210  , _nsamples(0)
211  , template_input_file("/gpfs/mnt/gpfs02/sphenix/user/trinn/fitting_algorithm_playing/prdfcode/prototype/offline/packages/Prototype4/templates.root")
212  , _calib_params(name)
213  , _fit_type(kPowerLawDoubleExpWithGlobalFitConstraint)
214 {
215  SetDefaultParameters(_calib_params);
216 }
217 
218 //_____________________________________
219 int CaloTemplateFit::InitRun(PHCompositeNode *topNode)
220 {
221  CreateNodeTree(topNode);
222 
223  if (Verbosity())
224  {
225  std::cout << Name() << "::" << detector << "::" << __PRETTY_FUNCTION__
226  << " - print calibration parameters: " << endl;
227  _calib_params.Print();
228  }
229 
230 
231  TFile* fin = TFile::Open(template_input_file.c_str());
232  assert(fin);
233  assert(fin->IsOpen());
234  h_template = static_cast<TProfile*>(fin->Get("hp_electrons_fine_emcal_36_8GeV"));
235  h_template->SetDirectory(0);
236  fin->Close();
237  delete fin;
238 
239 
240  h_data = new TH1F("h_data","",32,0,32);
241  f_fit = new TF1("f_fit",template_function,0,31,3);
242 
243 
244  ROOT::EnableThreadSafety();
245 
246  return Fun4AllReturnCodes::EVENT_OK;
247 }
248 
249 //____________________________________
250 int CaloTemplateFit::process_event(PHCompositeNode *topNode)
251 {
252  if (Verbosity())
253  {
254  std::cout << Name() << "::" << detector << "::" << __PRETTY_FUNCTION__
255  << "Process event entered" << std::endl;
256  }
257 
258  map<int, double> parameters_constraints;
259 
260  if ( _raw_towers->size() > 1)
261  {
262  if (Verbosity())
263  {
264  std::cout
265  << Name() << "::" << detector << "::" << __PRETTY_FUNCTION__
266  << "Extract global fit parameter for constraining individual fits"
267  << std::endl;
268  }
269  // signal template
270  vector<float> waveforms;
271  vector<vector<float>> waveforms2;
272  int count = 0;
273  RawTowerContainer::Range begin_end = _raw_towers->getTowers();
274  RawTowerContainer::Iterator rtiter;
275  for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
276  {
277  RawTower_Prototype4 *raw_tower =
278  dynamic_cast<RawTower_Prototype4 *>(rtiter->second);
279  assert(raw_tower);
280  ++count;
281 
282  for (int i = 0; i < RawTower_Prototype4::NSAMPLES; i++)
283  {
284  if (i >= _nsamples && _nsamples != 0){continue;}
285  waveforms.push_back(raw_tower->get_signal_samples(i));
286  }
287 
288  // std::cout << waveforms.size() << std::endl;
289  waveforms.push_back(count);
290  waveforms2.push_back(waveforms);
291  waveforms.clear();
292  }
293  std::vector< std::vector<float>> pulse_quantification = calo_processing_perchnl(waveforms2);
294 
295  int towernumber = 0;
296 
297  for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
298  {
299  // RawTowerDefs::keytype key = rtiter->first;
300  RawTower_Prototype4 *raw_tower =
301  dynamic_cast<RawTower_Prototype4 *>(rtiter->second);
302  assert(raw_tower);
303  // store the result - raw_tower
304  if (std::isnan(raw_tower->get_energy()))
305  {
306  // Raw tower was never fit, store the current fit
307  std::vector<float> values = pulse_quantification.at(towernumber);
308  raw_tower->set_energy(values.at(0));
309  raw_tower->set_time(values.at(1));
310  values.clear();
311  }
312  towernumber++;
313  }
314  pulse_quantification.clear();
315  waveforms2.clear();
316  // std::cout <<waveforms.size() << " , " << waveforms2.size() << " , "<< pulse_quantification.size() << std::endl;
317  }
318 
319 
320  /*
321  if (count > 0)
322  {
323  for (int i = 0; i < RawTower_Prototype4::NSAMPLES; i++)
324  {
325  vec_signal_samples[i] /= count;
326  }
327 
328  double peak = NAN;
329  double peak_sample = NAN;
330  double pedstal = NAN;
331  map<int, double> parameters_io;
332 
333  PROTOTYPE4_FEM::SampleFit_PowerLawDoubleExp(vec_signal_samples, peak,
334  peak_sample, pedstal,
335  parameters_io, Verbosity());
336  // std::map<int, double> &parameters_io, //! IO for fullset of
337  // parameters. If a parameter exist and not an NAN, the fit parameter
338  // will be fixed to that value. The order of the parameters are
339  // ("Amplitude", "Sample Start", "Power", "Peak Time 1", "Pedestal",
340  // "Amplitude ratio", "Peak Time 2")
341 
342  parameters_constraints[1] = parameters_io[1];
343  parameters_constraints[2] = parameters_io[2];
344  parameters_constraints[3] = parameters_io[3];
345  parameters_constraints[5] = parameters_io[5];
346  parameters_constraints[6] = parameters_io[6];
347 
348  // //special constraint if Peak Time 1 == Peak Time 2
349  // if (abs(parameters_constraints[6] - parameters_constraints[3]) <
350  // 0.1)
351  // {
352  // const double average_peak_time = (parameters_constraints[6] +
353  // parameters_constraints[3]) / 2.;
354  //
355  // std::cout << Name() << "::" << detector << "::" <<
356  // __PRETTY_FUNCTION__
357  // << ": two shaping time are too close "
358  // << parameters_constraints[3] << " VS " <<
359  // parameters_constraints[6]
360  // << ". Use average peak time instead: " <<
361  // average_peak_time
362  // << std::endl;
363  //
364  // parameters_constraints[6] = average_peak_time;
365  // parameters_constraints[3] = average_peak_time;
366  // parameters_constraints[5] = 0;
367  // }
368  }
369  else
370  {
371  std::cout << Name() << "::" << detector << "::" << __PRETTY_FUNCTION__
372  << ": Failed to build signal template! Fit each channel "
373  "individually instead"
374  << std::endl;
375  }
376  }
377  */
378  /*
379  const double calib_const_scale =
380  _calib_params.get_double_param("calib_const_scale");
381  const bool use_chan_calibration =
382  _calib_params.get_int_param("use_chan_calibration") > 0;
383 
384  RawTowerContainer::Range begin_end = _raw_towers->getTowers();
385  RawTowerContainer::Iterator rtiter;
386  for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
387  {
388  RawTowerDefs::keytype key = rtiter->first;
389  RawTower_Prototype4 *raw_tower =
390  dynamic_cast<RawTower_Prototype4 *>(rtiter->second);
391  assert(raw_tower);
392 
393  double calibration_const = calib_const_scale;
394 
395  if (use_chan_calibration)
396  {
397  // channel to channel calibration.
398  const int column = raw_tower->get_column();
399  const int row = raw_tower->get_row();
400 
401  assert(column >= 0);
402  assert(row >= 0);
403 
404  string calib_const_name(
405  (boost::format("calib_const_column%d_row%d") % column % row).str());
406 
407  calibration_const *= _calib_params.get_double_param(calib_const_name);
408  }
409 
410  vector<double> vec_signal_samples;
411  for (int i = 0; i < RawTower_Prototype4::NSAMPLES; i++)
412  {
413  vec_signal_samples.push_back(raw_tower->get_signal_samples(i));
414  }
415 
416  double peak = NAN;
417  double peak_sample = NAN;
418  double pedstal = NAN;
419 
420  switch (_fit_type)
421  {
422  case kPowerLawExp:
423  PROTOTYPE4_FEM::SampleFit_PowerLawExp(vec_signal_samples, peak,
424  peak_sample, pedstal, Verbosity());
425  break;
426 
427  case kPeakSample:
428  PROTOTYPE4_FEM::SampleFit_PeakSample(vec_signal_samples, peak,
429  peak_sample, pedstal, Verbosity());
430  break;
431 
432  case kPowerLawDoubleExp:
433  {
434  map<int, double> parameters_io;
435 
436  PROTOTYPE4_FEM::SampleFit_PowerLawDoubleExp(vec_signal_samples, peak,
437  peak_sample, pedstal,
438  parameters_io, Verbosity());
439  }
440  break;
441 
442  case kPowerLawDoubleExpWithGlobalFitConstraint:
443  {
444  map<int, double> parameters_io(parameters_constraints);
445 
446  PROTOTYPE4_FEM::SampleFit_PowerLawDoubleExp(vec_signal_samples, peak,
447  peak_sample, pedstal,
448  parameters_io, Verbosity());
449  }
450  break;
451  default:
452  cout << __PRETTY_FUNCTION__ << " - FATAL error - unkown fit type "
453  << _fit_type << endl;
454  exit(3);
455  break;
456  }
457 
458  // store the result - raw_tower
459  if (std::isnan(raw_tower->get_energy()))
460  {
461  // Raw tower was never fit, store the current fit
462 
463  raw_tower->set_energy(peak);
464  raw_tower->set_time(peak_sample);
465 
466  if (Verbosity())
467  {
468  raw_tower->identify();
469  }
470  }
471 
472  // store the result - calib_tower
473  RawTower_Prototype4 *calib_tower = new RawTower_Prototype4(*raw_tower);
474  calib_tower->set_energy(peak * calibration_const);
475  calib_tower->set_time(peak_sample);
476 
477  for (int i = 0; i < RawTower_Prototype4::NSAMPLES; i++)
478  {
479  calib_tower->set_signal_samples(i, (vec_signal_samples[i] - pedstal) *
480  calibration_const);
481  }
482 
483  _calib_towers->AddTower(key, calib_tower);
484 
485  } // for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
486 
487  if (Verbosity())
488  {
489  std::cout << Name() << "::" << detector << "::" << __PRETTY_FUNCTION__
490  << "input sum energy = " << _raw_towers->getTotalEdep()
491  << ", output sum digitalized value = "
492  << _calib_towers->getTotalEdep() << std::endl;
493  }
494 */
495 
496 
497  return Fun4AllReturnCodes::EVENT_OK;
498 }
499 
500 //_______________________________________
501 void CaloTemplateFit::CreateNodeTree(PHCompositeNode *topNode)
502 {
503  PHNodeIterator iter(topNode);
504 
505  PHCompositeNode *dstNode =
506  dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
507  if (!dstNode)
508  {
509  std::cerr << Name() << "::" << detector << "::" << __PRETTY_FUNCTION__
510  << "DST Node missing, doing nothing." << std::endl;
511  throw std::runtime_error(
512  "Failed to find DST node in RawTowerCalibration::CreateNodes");
513  }
514 
515  RawTowerNodeName = "TOWER_" + _raw_tower_node_prefix + "_" + detector;
516  _raw_towers =
517  findNode::getClass<RawTowerContainer>(dstNode, RawTowerNodeName.c_str());
518  if (!_raw_towers)
519  {
520  std::cerr << Name() << "::" << detector << "::" << __PRETTY_FUNCTION__
521  << " " << RawTowerNodeName << " Node missing, doing bail out!"
522  << std::endl;
523  throw std::runtime_error("Failed to find " + RawTowerNodeName +
524  " node in RawTowerCalibration::CreateNodes");
525  }
526 
527  // Create the tower nodes on the tree
528  PHNodeIterator dstiter(dstNode);
529  PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(
530  dstiter.findFirst("PHCompositeNode", detector));
531  if (!DetNode)
532  {
533  DetNode = new PHCompositeNode(detector);
534  dstNode->addNode(DetNode);
535  }
536 
537  // Be careful as a previous calibrator may have been registered for this
538  // detector
539  CaliTowerNodeName = "TOWER_" + _calib_tower_node_prefix + "_" + detector;
540  _calib_towers =
541  findNode::getClass<RawTowerContainer>(DetNode, CaliTowerNodeName.c_str());
542  if (!_calib_towers)
543  {
544  _calib_towers = new RawTowerContainer(_raw_towers->getCalorimeterID());
545  PHIODataNode<PHObject> *towerNode = new PHIODataNode<PHObject>(
546  _calib_towers, CaliTowerNodeName.c_str(), "PHObject");
547  DetNode->addNode(towerNode);
548  }
549 
550  // update the parameters on the node tree
551  PHCompositeNode *parNode =
552  dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
553  assert(parNode);
554  const string paramnodename = string("Calibration_") + detector;
555 
556  // this step is moved to after detector construction
557  // save updated persistant copy on node tree
558  _calib_params.SaveToNodeTree(parNode, paramnodename);
559 }
560 
561 void CaloTemplateFit::SetDefaultParameters(PHParameters &param)
562 {
563  param.set_int_param("use_chan_calibration", 0);
564 
565  // additional scale for the calibration constant
566  // negative pulse -> positive with -1
567  param.set_double_param("calib_const_scale", 1);
568 }