12 #include "TFitResult.h"
46 const char*
const muon_t =
"muon";
56 int main(
int argc,
char **argv)
60 using namespace KM3NETDAQ;
63 typedef JTriggeredFileScanner_t::multi_pointer_type multi_pointer_type;
67 JTriggeredFileScanner_t inputFile;
82 JParser<> zap(
"Utility program to determine energy correction.");
84 zap[
'f'] =
make_field(inputFile,
"event file(s) or histogram file");
86 zap[
'n'] =
make_field(numberOfEvents) = JLimit::max();
88 zap[
'F'] =
make_field(formula,
"fit formula, e.g: \"[0]+[1]*x\"") =
"";
92 zap[
'O'] =
make_field(option,
"Fit option") =
"";
99 catch(
const exception& error) {
100 FATAL(error.what() << endl);
107 if (inputFile.size() == 1) {
109 in = TFile::Open(inputFile[0].c_str(),
"READ");
111 if (
in != NULL &&
in->IsOpen()) {
112 h2 = (TH2D*)
in->Get(
"h2");
120 for (JMultipleFileScanner_t::const_iterator
i = inputFile.begin();
i != inputFile.end(); ++
i) {
128 if (p == enigma.end())
129 enigma.push_back(head);
133 catch(
const exception& error) {
138 double Emin = numeric_limits<double>::max();
143 if (!
i->is_valid(&JHead::cut_nu)) {
144 FATAL(
"Missing neutrino parameters." << endl);
147 if (
i->cut_nu.E.getLowerLimit() < Emin) { Emin =
i->cut_nu.E.getLowerLimit(); }
148 if (
i->cut_nu.E.getUpperLimit() > Emax) { Emax =
i->cut_nu.E.getUpperLimit(); }
150 if (
i->is_valid(&JHead::genvol)) {
151 cout <<
"Neutrino energy range: " <<
i->cut_nu.E << endl;
152 cout <<
"Number of generated events: " <<
i->genvol.numberOfEvents << endl;
163 Vec offset(0.0, 0.0, 0.0);
165 for (
string filename; inputFile.hasNext(); ) {
167 STATUS(
"event: " << setw(10) << inputFile.getCounter() <<
'\r');
DEBUG(endl);
169 if (filename != inputFile.getFilename()) {
171 filename = inputFile.getFilename();
181 if (p->is_valid(&JHead::genvol)) {
182 W = 1.0 / (double) p->genvol.numberOfEvents;
185 cout <<
"File: " << filename <<
" weight " << W << endl;
187 catch(
const exception& error) {
192 multi_pointer_type ps = inputFile.next();
226 if (evt->begin() == __end) {
230 JEvt::iterator p = min_element(evt->begin(), __end,
quality_sorter);
240 if (ta.getE() > 0.0 && tb.
getE() > 0.0) {
241 h2->Fill(
log10(tb.
getE()),
log10(ta.getE()), (enigma.size() > (size_t) 1 ? W * event->w[2] : 1.0));
249 const Int_t nx = h2->GetXaxis()->GetNbins();
250 const Double_t
xmin = h2->GetXaxis()->GetXmin();
251 const Double_t
xmax = h2->GetXaxis()->GetXmax();
253 const Int_t ny = h2->GetYaxis()->GetNbins();
255 TH1D h1(
"h1", NULL, nx,
xmin,
xmax);
263 for (Int_t ix = 1; ix <= nx; ++ix) {
269 for (
int iy = 1; iy <= ny; ++iy) {
271 const Double_t
x = h2->GetYaxis()->GetBinCenter(iy);
272 const Double_t
y = h2->GetBinContent(ix,iy);
273 const Double_t z = h2->GetBinError (ix,iy);
277 h0->SetBinContent(iy, y);
278 h0->SetBinError (iy, z);
290 TF1
f1(
"f1",
"[0]*TMath::Gaus(x,[1],[2]) + [3]");
295 f1.SetParameter(0, Q.
getWmax());
296 f1.SetParameter(1, x0);
297 f1.SetParameter(2, sigma);
298 f1.SetParameter(3, 0.0);
300 TFitResultPtr
result = h0->Fit(&f1,
"SQ",
"same");
302 if (result.Get() == NULL) {
303 FATAL(
"Invalid TFitResultPtr " << h0->GetName() << endl);
308 << setw(4) << ix <<
' '
309 <<
FIXED(7,3) << h2->GetXaxis()->GetBinCenter(ix) <<
' '
310 <<
FIXED(7,3) << result->Chi2() <<
'/'
311 << result->Ndf() <<
' '
312 << (result->IsValid() ?
"" :
"failed") << endl;
315 if (result->IsValid()) {
317 const Double_t
x = h2->GetXaxis()->GetBinCenter(ix);
319 h1.SetBinContent(ix, f1.GetParameter(1));
320 h1.SetBinError (ix, f1.GetParError (1));
322 g1.put(x,
pow(10.0, x) -
pow(10.0, f1.GetParameter(1)));
323 g2.put(x,
pow(10.0, abs(f1.GetParameter(2)) - f1.GetParameter(1)));
324 g3.put(x, f1.GetParameter(1), 0.0, f1.GetParameter(2));
333 f1 =
new TF1(
"f1", formula.c_str());
335 if (
f1->IsZombie()) {
336 FATAL(
"Function: " << formula <<
" is zombie." << endl);
350 FATAL(error << endl);
353 if (option.find(
'S') == string::npos) { option +=
"S"; }
355 TFitResultPtr
result = h1.Fit(
f1, option.c_str(),
"same",
X.getLowerLimit(),
X.getUpperLimit());
361 Double_t z[2] = { 0.00, 0.00 };
373 for (Int_t ix = nx; ix >= 1; --ix) {
375 x[0] = h1.GetXaxis()->GetBinCenter(ix);
376 y[0] = h1.GetBinContent(ix);
377 z[0] = h1.GetBinError (ix);
379 if (z[0] == 0.0 || !
X(x[0])) {
385 <<
"x >= " <<
FIXED(5,3) << x[0]
387 <<
"x < " <<
FIXED(5,3) << x[1]
390 <<
"(" <<
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]) <<
")";
403 <<
"x < " <<
FIXED(5,3) << x[1]
406 <<
"(" <<
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]) <<
")";
409 string buffer = os.str();
413 for (string::iterator
i = buffer.begin();
i != buffer.end(); ++
i) {
419 f1 =
new TF1(
"f1", buffer.c_str(),
xmin,
xmax);
421 h1.GetListOfFunctions()->Add(
f1);
427 out <<
JMeta(argc, argv);
429 out << *h2 << H0 << h1;
439 out << *(
f1->GetFormula());
Utility class to parse command line options.
Q(UTCMax_s-UTCMin_s)-livetime_s
double getValue(const JScale_t scale)
Get numerical value corresponding to scale.
int main(int argc, char *argv[])
int getParameter(const std::string &text)
Get parameter number from text string.
ROOT TTree parameter settings of various packages.
Vec getOffset(const JHead &header)
Get offset.
JTrack3E getTrack(const Trk &track)
Get track.
static const int JENERGY_ENERGY
uncorrected energy [GeV] from JEnergy.cc
Synchronously read DAQ events and Monte Carlo events (and optionally other events).
Data structure for graph data.
const char *const neutrino_t
bool has_neutrino(const Evt &evt)
Test whether given event has an incoming neutrino.
double getIntersection(const JVector3D &pos) const
Get longitudinal position along axis of position of closest approach with given position.
bool is_muon(const Trk &track)
Test whether given track is a (anti-)muon.
General purpose sorter of fit results.
Auxiliary class to synchronously read DAQ events and Monte Carlo events (and optionally other events)...
Data structure for graph data.
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
static const JGeaneWater gWater
Function object for energy loss of muon in sea water.
Dynamic ROOT object management.
Auxiliary data structure for floating point format specification.
Auxiliary data structure to build TGraph.
double E
Energy [GeV] (either MC truth or reconstructed)
void move(const double step, const double velocity, const JGeane &geane)
Move vertex along this track with given velocity.
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
bool is_neutrino(const Trk &track)
Test whether given track is a neutrino.
const JPolynome f1(1.0, 2.0, 3.0)
Function.
Auxiliary class for defining the range of iterations of objects.
The Vec class is a straightforward 3-d vector, which also works in pyroot.
Type definition of range.
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys...
double getE() const
Get energy.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
set_variable E_E log10(E_{fit}/E_{#mu})"
Auxiliary class to test history.
void setE(const double E)
Set energy.
JPosition3D getPosition(const Vec &pos)
Get position.
T pow(const T &x, const double y)
Power .
Reconstruction type dependent comparison of track quality.
General purpose messaging.
Auxiliary data structure to build TGraphErrors.
const double getSpeedOfLight()
Get speed of light.
then fatal The output file must have the wildcard in the e g root fi eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY JAcoustics sh $DETECTOR_ID source JAcousticsToolkit sh CHECK_EXIT_CODE typeset A EMITTERS get_tripods $WORKDIR tripod txt EMITTERS get_transmitters $WORKDIR transmitter txt EMITTERS for EMITTER in
Utility class to parse command line options.
Wrapper class around string.
Exception for parsing value.
no fit printf nominal n $STRING awk v X
const char * getName()
Get ROOT name of given data type.
const JLimit & getLimit() const
Get limit.
const Trk & get_neutrino(const Evt &evt)
Get incoming neutrino.
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
int numberOfBins
number of bins for average CDF integral of optical module
static const int JMUONENERGY
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
#define DEBUG(A)
Message macros.
Double_t g1(const Double_t x)
Function.