Jpp  15.0.0
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
software/JReconstruction/JEnergyCorrection.cc File Reference

Utility program to determine energy correction. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TF1.h"
#include "TMath.h"
#include "TFitResult.h"
#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/definitions/reconstruction.hh"
#include "km3net-dataformat/definitions/fitparameters.hh"
#include "JAAnet/JHead.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JDAQ/JDAQEventIO.hh"
#include "JSupport/JTriggeredFileScanner.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JMeta.hh"
#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JEvtToolkit.hh"
#include "JReconstruction/JEnergyCorrection.hh"
#include "JGeometry3D/JTrack3E.hh"
#include "JPhysics/JGeane.hh"
#include "JTools/JQuantile.hh"
#include "JROOT/JRootToolkit.hh"
#include "JGizmo/JGizmoToolkit.hh"
#include "JLang/JToken.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Detailed Description

Utility program to determine energy correction.

Author
mdejong

Definition in file software/JReconstruction/JEnergyCorrection.cc.

Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 46 of file software/JReconstruction/JEnergyCorrection.cc.

47 {
48  using namespace std;
49  using namespace JPP;
50  using namespace KM3NETDAQ;
51 
52  typedef JTriggeredFileScanner<JEvt> JTriggeredFileScanner_t;
53  typedef JTriggeredFileScanner_t::multi_pointer_type multi_pointer_type;
54  typedef JToken<';'> JToken_t;
55 
56  JTriggeredFileScanner_t inputFile;
57  JLimit_t& numberOfEvents = inputFile.getLimit();
58  string outputFile;
59  JQualitySorter quality_sorter;
60  string formula;
61  vector<JToken_t> startValues;
62  vector<JToken_t> fixedValues;
63  string option;
64  int debug;
65 
66  try {
67 
68  JParser<> zap("Utility program to determine energy correction.");
69 
70  zap['f'] = make_field(inputFile, "event file(s) or histogram file");
71  zap['o'] = make_field(outputFile) = "energy_correction.root";
72  zap['n'] = make_field(numberOfEvents) = JLimit::max();
73  zap['L'] = make_field(quality_sorter) = JPARSER::initialised();
74  zap['F'] = make_field(formula, "fit formula, e.g: \"[0]+[1]*x\"") = "";
75  zap['@'] = make_field(startValues, "start values, e.g: \"p0 = GetMaximum;\"") = JPARSER::initialised();
76  zap['='] = make_field(fixedValues, "fixed values, e.g: \"p0 = GetMaximum;\"") = JPARSER::initialised();
77  zap['O'] = make_field(option, "Fit option") = "";
78  zap['d'] = make_field(debug) = 2;
79 
80  zap(argc, argv);
81  }
82  catch(const exception& error) {
83  FATAL(error.what() << endl);
84  }
85 
86 
87  TFile* in = NULL;
88  TH2D* h2 = NULL;
89 
90  if (inputFile.size() == 1) {
91 
92  in = TFile::Open(inputFile[0].c_str(), "READ");
93 
94  if (in != NULL && in->IsOpen()) {
95  h2 = (TH2D*) in->Get("h2");
96  }
97  }
98 
99  if (h2 == NULL) {
100 
101  JHead head;
102 
103  try {
104  head = getHeader(inputFile);
105  }
106  catch(const JException& error) {
107  FATAL(error);
108  }
109 
110  if (!head.is_valid(&JHead::cut_nu)) {
111  FATAL("Missing neutrino parameters." << endl);
112  }
113 
114  const JPosition3D center = get<JPosition3D>(head);
115 
116  const double Emin = head.cut_nu.Emin;
117  const double Emax = head.cut_nu.Emax;
118 
119  const Int_t nx = 40;
120  const Double_t xmin = log10(Emin);
121  const Double_t xmax = log10(Emax);
122 
123  h2 = new TH2D("h2", NULL, nx, xmin, xmax, nx, xmin, xmax);
124 
125  while (inputFile.hasNext()) {
126 
127  STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
128 
129  multi_pointer_type ps = inputFile.next();
130 
131  JEvt* evt = ps;
132  Evt* event = ps;
133 
134  double Emax = 0.0;
135 
136  vector<Trk>::const_iterator muon = event->mc_trks.end();
137 
138  for (vector<Trk>::const_iterator track = event->mc_trks.begin(); track != event->mc_trks.end(); ++track) {
139  if (is_muon(*track)) {
140  if (track->E > Emax) {
141  muon = track;
142  }
143  }
144  }
145 
146  if (muon == event->mc_trks.end()) {
147  continue;
148  }
149 
150  JEvt::iterator __end = partition(evt->begin(), evt->end(), JHistory::is_event(JMUONENERGY));
151 
152  if (evt->begin() == __end) {
153  continue;
154  }
155 
156  JEvt::iterator p = min_element(evt->begin(), __end, quality_sorter);
157 
158  JTrack3E ta = getTrack(*muon);
159  JTrack3E tb = getTrack(*p);
160 
161  ta.add(center);
162 
163  tb.setE(p->getW(JENERGY_ENERGY));
164  tb.move(tb.getIntersection(ta.getPosition()), getSpeedOfLight(), gWater);
165 
166  if (ta.getE() > 0.0 && tb.getE() > 0.0) {
167  h2->Fill(log10(tb.getE()), log10(ta.getE()));
168  }
169  }
170  STATUS(endl);
171  }
172 
173  // short hands
174 
175  const Int_t nx = h2->GetXaxis()->GetNbins();
176  const Double_t xmin = h2->GetXaxis()->GetXmin();
177  const Double_t xmax = h2->GetXaxis()->GetXmax();
178 
179  const Int_t ny = h2->GetYaxis()->GetNbins();
180  const Double_t ymin = h2->GetYaxis()->GetXmin();
181  const Double_t ymax = h2->GetYaxis()->GetXmax();
182 
183  // 1D profile
184 
185  TH1D h1("h1", NULL, nx, xmin, xmax);
186 
187  for (Int_t ix = 1; ix <= nx; ++ix) {
188 
189  TH1D h0("h0", NULL, ny, ymin, ymax);
190 
191  JQuantile Q; // start values
192 
193  for (int iy = 1; iy <= ny; ++iy) {
194 
195  const Double_t x = h2->GetYaxis()->GetBinCenter(iy);
196  const Double_t y = h2->GetBinContent(ix,iy);
197  const Double_t z = h2->GetBinError (ix,iy);
198 
199  if (y > 0.0) {
200 
201  h0.SetBinContent(iy, y);
202  h0.SetBinError (iy, z);
203 
204  Q.put(x,y);
205  }
206  }
207 
208  if (Q.getCount() >= 3) {
209 
210  if (debug >= debug_t) {
211  Q.print(cout);
212  }
213 
214  TF1 f1("f1", "[0]*TMath::Gaus(x,[1],[2])");
215 
216  f1.SetParameter(0, Q.getWmax());
217  f1.SetParameter(1, Q.getMean());
218  f1.SetParameter(2, Q.getSTDev());
219 
220  TFitResultPtr result = h0.Fit(&f1, "SLLQ", "same",
221  Q.getMean() - 3.0 * Q.getSTDev(),
222  Q.getMean() + 3.0 * Q.getSTDev());
223 
224  if (result.Get() == NULL) {
225  FATAL("Invalid TFitResultPtr " << h0.GetName() << endl);
226  }
227 
228  if (debug >= debug_t || !result->IsValid()) {
229  cout << "Bin: "
230  << setw(4) << ix << ' '
231  << FIXED(7,3) << h2->GetXaxis()->GetBinCenter(ix) << ' '
232  << FIXED(7,3) << result->Chi2() << '/'
233  << result->Ndf() << ' '
234  << (result->IsValid() ? "" : "failed") << endl;
235  }
236 
237  if (result->IsValid()) {
238  h1.SetBinContent(ix, f1.GetParameter(1));
239  h1.SetBinError (ix, f1.GetParError (1));
240  }
241  }
242  }
243 
244  TF1* f1 = NULL;
245 
246  if (formula != "") { // fit
247 
248  f1 = new TF1("f1", formula.c_str());
249 
250  if (f1->IsZombie()) {
251  FATAL("Function: " << formula << " is zombie." << endl);
252  }
253 
254  try {
255 
256  for (vector<JToken_t>::const_iterator i = startValues.begin(); i != startValues.end(); ++i) {
257  f1->SetParameter(getParameter(*i), getValue(*i,&h1));
258  }
259 
260  for (vector<JToken_t>::const_iterator i = fixedValues.begin(); i != fixedValues.end(); ++i) {
261  f1->FixParameter(getParameter(*i), getValue(*i,&h1));
262  }
263  }
264  catch(JLANG::JParseError& error) {
265  FATAL(error << endl);
266  }
267 
268  if (option.find('S') == string::npos) { option += "S"; }
269 
270  TFitResultPtr result = h1.Fit(f1, option.c_str(), "same");
271 
272  } else { // interpolate
273 
274  Double_t x[2] = { xmin, xmax };
275  Double_t y[2] = { xmin, xmax };
276  Double_t z[2] = { 0.00, 0.00 };
277 
278  ostringstream os;
279 
280  os << " " // extrapolate
281  << "("
282  << "x >= " << FIXED(5,3) << xmax
283  << ")"
284  << " * "
285  << "(x)";
286  os << "\n";
287 
288  for (Int_t ix = nx; ix >= 1; --ix) {
289 
290  x[0] = h1.GetXaxis()->GetBinCenter(ix);
291  y[0] = h1.GetBinContent(ix);
292  z[0] = h1.GetBinError (ix);
293 
294  if (z[0] == 0.0 || y[1] < y[0]) {
295 
296  if (x[1] == xmax)
297  continue; // continue
298  else
299  ix = 0; // terminate
300 
301  x[0] = xmin;
302  y[0] = xmin;
303  }
304 
305  os << " + "
306  << "("
307  << "x >= " << FIXED(5,3) << x[0]
308  << " && "
309  << "x < " << FIXED(5,3) << x[1]
310  << ")"
311  << " * "
312  << "(" << FIXED(5,3) << y[0] << " + " << FIXED(5,3) << (y[1] - y[0]) << " * (x - " << FIXED(5,3) << x[0] << ") / " << FIXED(5,3) << (x[1] - x[0]) << ")";
313  os << "\n";
314 
315  x[1] = x[0];
316  y[1] = y[0];
317  z[1] = z[0];
318  }
319 
320  os << " + " // extrapolate
321  << "("
322  << "x < " << FIXED(5,3) << xmin
323  << ")"
324  << " * "
325  << "(x)";
326  os << "\n";
327 
328  string buffer = os.str();
329 
330  NOTICE(buffer << endl);
331 
332  for (string::iterator i = buffer.begin(); i != buffer.end(); ++i) {
333  if (*i == '\n') {
334  *i = ' ';
335  }
336  }
337 
338  f1 = new TF1("f1", buffer.c_str(), xmin, xmax);
339 
340  h1.GetListOfFunctions()->Add(f1);
341  }
342 
343 
344  TFile out(outputFile.c_str(), "recreate");
345 
346  out << JMeta(argc, argv);
347  out << h1 << *h2;
348 
349  if (f1 != NULL) { // separately store energy correction
350 
351  f1->GetFormula()->SetName(JEnergyCorrection::getName());
352 
353  out << *(f1->GetFormula());
354  }
355 
356  out.Write();
357  out.Close();
358 
359  if (in != NULL) {
360  in->Close();
361  }
362 }
JSUPPORT::JMeta
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:70
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1500
JLANG::JException
General exception.
Definition: JException.hh:23
Q
Q(UTCMax_s-UTCMin_s)-livetime_s
JEEP::getValue
double getValue(const JScale_t scale)
Get numerical value corresponding to scale.
Definition: JScale.hh:47
JEEP::debug_t
debug
Definition: JMessage.hh:29
JGIZMO::getParameter
int getParameter(const std::string &text)
Get parameter number from text string.
Definition: JGizmoToolkit.hh:327
JAANET::getTrack
JTrack3E getTrack(const Trk &track)
Get track.
Definition: JAAnetToolkit.hh:246
JENERGY_ENERGY
static const int JENERGY_ENERGY
uncorrected energy [GeV] from JEnergy.cc
Definition: fitparameters.hh:16
JTOOLS::JQuantile::getSTDev
double getSTDev() const
Get standard deviation.
Definition: JQuantile.hh:266
JTOOLS::JQuantile
Auxiliary data structure for running average, standard deviation and quantiles.
Definition: JQuantile.hh:43
STATUS
#define STATUS(A)
Definition: JMessage.hh:63
JGEOMETRY3D::JAxis3D::getIntersection
double getIntersection(const JVector3D &pos) const
Get longitudinal position along axis of position of closest approach with given position.
Definition: JAxis3D.hh:146
JAANET::is_muon
bool is_muon(const Trk &track)
Test whether given track is a (anti-)muon.
Definition: JAAnetToolkit.hh:326
JRECONSTRUCTION::JQualitySorter
General purpose sorter of fit results.
Definition: JReconstruction/JEvtToolkit.hh:259
JSUPPORT::JTriggeredFileScanner
Auxiliary class to synchronously read DAQ events and Monte Carlo events (and optionally other events)...
Definition: JTriggeredFileScanner.hh:39
h1
then for HISTOGRAM in h0 h1
Definition: JMatrixNZ.sh:71
JAANET::cut::Emax
double Emax
Maximal energy [GeV].
Definition: JHead.hh:352
JPARSER::initialised
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:66
JPHYSICS::gWater
static const JGeaneWater gWater
Function object for energy loss of muon in sea water.
Definition: JGeane.hh:323
FIXED
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
JGEOMETRY3D::JTrack3E
3D track with energy.
Definition: JTrack3E.hh:30
JTOOLS::JQuantile::getMean
double getMean() const
Get mean value.
Definition: JQuantile.hh:252
JGEOMETRY3D::JTrack3E::move
void move(const double step, const double velocity, const JGeane &geane)
Move vertex along this track with given velocity.
Definition: JTrack3E.hh:117
JTOOLS::JQuantile::getWmax
double getWmax() const
Get maximum weight.
Definition: JQuantile.hh:241
JSUPPORT::getHeader
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
Definition: JMonteCarloFileSupportkit.hh:471
JSUPPORT::JLimit
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
JGEOMETRY3D::JTrack3E::getE
double getE() const
Get energy.
Definition: JTrack3E.hh:93
JACOUSTICS::JEvt
Acoustic event fit.
Definition: JAcoustics/JEvt.hh:154
JTOOLS::JQuantile::print
std::ostream & print(std::ostream &out, bool lpr=true) const
Print quantile.
Definition: JQuantile.hh:335
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
JFIT::JHistory::is_event
Auxiliary class to test history.
Definition: JHistory.hh:104
JTOOLS::result
return result
Definition: JPolint.hh:727
JGEOMETRY3D::JTrack3E::setE
void setE(const double E)
Set energy.
Definition: JTrack3E.hh:104
NOTICE
#define NOTICE(A)
Definition: JMessage.hh:64
JTOOLS::JQuantile::put
void put(const double x, const double w=1.0)
Put value.
Definition: JQuantile.hh:133
debug
int debug
debug level
Definition: JSirene.cc:63
JTOOLS::JQuantile::getCount
long long int getCount() const
Get total count.
Definition: JQuantile.hh:186
quality_sorter
Reconstruction type dependent comparison of track quality.
Definition: tools/reconstruction.hh:85
JAANET::JHead
Monte Carlo run header.
Definition: JHead.hh:1113
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
JPHYSICS::getSpeedOfLight
const double getSpeedOfLight()
Get speed of light.
Definition: JPhysics/JConstants.hh:105
JLANG::JToken
Wrapper class around string.
Definition: JToken.hh:23
DEBUG
then usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
Definition: JCalibratePMT.sh:21
JLANG::JParseError
Exception for parsing value.
Definition: JException.hh:180
JROOT::getName
const char * getName()
Get ROOT name of given data type.
Definition: JRootToolkit.hh:57
JGEOMETRY3D::JPosition3D
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
JSUPPORT::JLimit::getLimit
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:73
JAANET::JHead::cut_nu
JAANET::cut_nu cut_nu
Definition: JHead.hh:1404
JAANET::JHead::is_valid
bool is_valid(T JHead::*pd) const
Check validity of given data member in JHead.
Definition: JHead.hh:1183
in
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable INPUT_FILE $argv[2] eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY source JAcoustics sh $DETECTOR_ID CHECK_EXIT_CODE typeset A TRIPODS get_tripods $WORKDIR tripod txt TRIPODS for EMITTER in
Definition: JCanberra.sh:41
JAANET::cut::Emin
double Emin
Minimal energy [GeV].
Definition: JHead.hh:351
JMUONENERGY
static const int JMUONENERGY
Definition: definitions/reconstruction.hh:17
Evt
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition: Evt.hh:19
Generated by   doxygen 1.8.5