Jpp  pmt_effective_area_update
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 (debug >= debug_t || !result->IsValid()) {
225  cout << "Bin: "
226  << setw(4) << ix << ' '
227  << FIXED(7,3) << h2->GetXaxis()->GetBinCenter(ix) << ' '
228  << FIXED(7,3) << result->Chi2() << '/'
229  << result->Ndf() << ' '
230  << (result->IsValid() ? "" : "failed") << endl;
231  }
232 
233  if (result->IsValid()) {
234  h1.SetBinContent(ix, f1.GetParameter(1));
235  h1.SetBinError (ix, f1.GetParError (1));
236  }
237  }
238  }
239 
240  TF1* f1 = NULL;
241 
242  if (formula != "") { // fit
243 
244  f1 = new TF1("f1", formula.c_str());
245 
246  if (f1->IsZombie()) {
247  FATAL("Function: " << formula << " is zombie." << endl);
248  }
249 
250  try {
251 
252  for (vector<JToken_t>::const_iterator i = startValues.begin(); i != startValues.end(); ++i) {
253  f1->SetParameter(getParameter(*i), getValue(*i,&h1));
254  }
255 
256  for (vector<JToken_t>::const_iterator i = fixedValues.begin(); i != fixedValues.end(); ++i) {
257  f1->FixParameter(getParameter(*i), getValue(*i,&h1));
258  }
259  }
260  catch(JLANG::JParseError& error) {
261  FATAL(error << endl);
262  }
263 
264  if (option.find('S') == string::npos) { option += "S"; }
265 
266  TFitResultPtr result = h1.Fit(f1, option.c_str(), "same");
267 
268  } else { // interpolate
269 
270  Double_t x[2] = { xmin, xmax };
271  Double_t y[2] = { xmin, xmax };
272  Double_t z[2] = { 0.00, 0.00 };
273 
274  ostringstream os;
275 
276  os << " " // extrapolate
277  << "("
278  << "x >= " << FIXED(5,3) << xmax
279  << ")"
280  << " * "
281  << "(x)";
282  os << "\n";
283 
284  for (Int_t ix = nx; ix >= 1; --ix) {
285 
286  x[0] = h1.GetXaxis()->GetBinCenter(ix);
287  y[0] = h1.GetBinContent(ix);
288  z[0] = h1.GetBinError (ix);
289 
290  if (z[0] == 0.0 || y[1] < y[0]) {
291 
292  if (x[1] == xmax)
293  continue; // continue
294  else
295  ix = 0; // terminate
296 
297  x[0] = xmin;
298  y[0] = xmin;
299  }
300 
301  os << " + "
302  << "("
303  << "x >= " << FIXED(5,3) << x[0]
304  << " && "
305  << "x < " << FIXED(5,3) << x[1]
306  << ")"
307  << " * "
308  << "(" << 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]) << ")";
309  os << "\n";
310 
311  x[1] = x[0];
312  y[1] = y[0];
313  z[1] = z[0];
314  }
315 
316  os << " + " // extrapolate
317  << "("
318  << "x < " << FIXED(5,3) << xmin
319  << ")"
320  << " * "
321  << "(x)";
322  os << "\n";
323 
324  string buffer = os.str();
325 
326  NOTICE(buffer << endl);
327 
328  for (string::iterator i = buffer.begin(); i != buffer.end(); ++i) {
329  if (*i == '\n') {
330  *i = ' ';
331  }
332  }
333 
334  f1 = new TF1("f1", buffer.c_str(), xmin, xmax);
335 
336  h1.GetListOfFunctions()->Add(f1);
337  }
338 
339 
340  TFile out(outputFile.c_str(), "recreate");
341 
342  out << JMeta(argc, argv);
343  out << h1 << *h2;
344 
345  if (f1 != NULL) { // separately store energy correction
346 
347  f1->GetFormula()->SetName(JEnergyCorrection::getName());
348 
349  out << *(f1->GetFormula());
350  }
351 
352  out.Write();
353  out.Close();
354 
355  if (in != NULL) {
356  in->Close();
357  }
358 }
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:241
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:321
JRECONSTRUCTION::JQualitySorter
General purpose sorter of fit results.
Definition: JReconstruction/JEvtToolkit.hh:265
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:54
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:40
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