12 #include "TFitResult.h"
46 int main(
int argc,
char **argv)
50 using namespace KM3NETDAQ;
53 typedef JTriggeredFileScanner_t::multi_pointer_type multi_pointer_type;
56 JTriggeredFileScanner_t inputFile;
68 JParser<> zap(
"Utility program to determine energy correction.");
70 zap[
'f'] =
make_field(inputFile,
"event file(s) or histogram file");
72 zap[
'n'] =
make_field(numberOfEvents) = JLimit::max();
74 zap[
'F'] =
make_field(formula,
"fit formula, e.g: \"[0]+[1]*x\"") =
"";
77 zap[
'O'] =
make_field(option,
"Fit option") =
"";
82 catch(
const exception& error) {
83 FATAL(error.what() << endl);
90 if (inputFile.size() == 1) {
92 in = TFile::Open(inputFile[0].c_str(),
"READ");
94 if (
in != NULL &&
in->IsOpen()) {
95 h2 = (TH2D*)
in->Get(
"h2");
110 if (!head.
is_valid(&JHead::cut_nu)) {
111 FATAL(
"Missing neutrino parameters." << endl);
120 const Double_t xmin = log10(Emin);
121 const Double_t xmax = log10(Emax);
123 h2 =
new TH2D(
"h2", NULL, nx, xmin, xmax, nx, xmin, xmax);
125 while (inputFile.hasNext()) {
127 STATUS(
"event: " << setw(10) << inputFile.getCounter() <<
'\r');
DEBUG(endl);
129 multi_pointer_type ps = inputFile.next();
140 if (track->E > Emax) {
146 if (muon == event->mc_trks.end()) {
152 if (evt->begin() == __end) {
156 JEvt::iterator p = min_element(evt->begin(), __end,
quality_sorter);
166 if (ta.getE() > 0.0 && tb.
getE() > 0.0) {
167 h2->Fill(log10(tb.
getE()), log10(ta.getE()));
175 const Int_t nx = h2->GetXaxis()->GetNbins();
176 const Double_t xmin = h2->GetXaxis()->GetXmin();
177 const Double_t xmax = h2->GetXaxis()->GetXmax();
179 const Int_t ny = h2->GetYaxis()->GetNbins();
180 const Double_t ymin = h2->GetYaxis()->GetXmin();
181 const Double_t ymax = h2->GetYaxis()->GetXmax();
185 TH1D
h1(
"h1", NULL, nx, xmin, xmax);
187 for (Int_t ix = 1; ix <= nx; ++ix) {
189 TH1D h0(
"h0", NULL, ny, ymin, ymax);
193 for (
int iy = 1; iy <= ny; ++iy) {
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);
201 h0.SetBinContent(iy, y);
202 h0.SetBinError (iy, z);
214 TF1 f1(
"f1",
"[0]*TMath::Gaus(x,[1],[2])");
216 f1.SetParameter(0, Q.
getWmax());
217 f1.SetParameter(1, Q.
getMean());
220 TFitResultPtr
result = h0.Fit(&f1,
"SLLQ",
"same",
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;
233 if (result->IsValid()) {
234 h1.SetBinContent(ix, f1.GetParameter(1));
235 h1.SetBinError (ix, f1.GetParError (1));
244 f1 =
new TF1(
"f1", formula.c_str());
246 if (f1->IsZombie()) {
247 FATAL(
"Function: " << formula <<
" is zombie." << endl);
261 FATAL(error << endl);
264 if (option.find(
'S') == string::npos) { option +=
"S"; }
266 TFitResultPtr
result = h1.Fit(f1, option.c_str(),
"same");
270 Double_t x[2] = { xmin, xmax };
271 Double_t y[2] = { xmin, xmax };
272 Double_t z[2] = { 0.00, 0.00 };
278 <<
"x >= " <<
FIXED(5,3) << xmax
284 for (Int_t ix = nx; ix >= 1; --ix) {
286 x[0] = h1.GetXaxis()->GetBinCenter(ix);
287 y[0] = h1.GetBinContent(ix);
288 z[0] = h1.GetBinError (ix);
290 if (z[0] == 0.0 || y[1] < y[0]) {
303 <<
"x >= " <<
FIXED(5,3) << x[0]
305 <<
"x < " <<
FIXED(5,3) << x[1]
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]) <<
")";
318 <<
"x < " <<
FIXED(5,3) << xmin
324 string buffer = os.str();
328 for (string::iterator i = buffer.begin(); i != buffer.end(); ++i) {
334 f1 =
new TF1(
"f1", buffer.c_str(), xmin, xmax);
336 h1.GetListOfFunctions()->Add(f1);
342 out <<
JMeta(argc, argv);
349 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.
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).
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)...
then for HISTOGRAM in h0 h1
double Emax
Maximal energy [GeV].
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.
Auxiliary data structure for floating point format specification.
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.
Auxiliary class for defining the range of iterations of objects.
double getE() const
Get energy.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Auxiliary class to test history.
void setE(const double E)
Set energy.
Reconstruction type dependent comparison of track quality.
General purpose messaging.
const double getSpeedOfLight()
Get speed of light.
Utility class to parse command line options.
Wrapper class around string.
then usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
Exception for parsing value.
const char * getName()
Get ROOT name of given data type.
Data structure for position in three dimensions.
const JLimit & getLimit() const
Get limit.
bool is_valid(T JHead::*pd) const
Check validity of given data member in JHead.
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
double Emin
Minimal energy [GeV].
static const int JMUONENERGY
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.