49int main(
int argc,
char **argv)
56 typedef JTriggeredFileScanner_t::multi_pointer_type multi_pointer_type;
57 typedef JToken<
';'> JToken_t;
60 JTriggeredFileScanner_t inputFile;
74 JParser<> zap(
"Utility program to determine energy correction.");
76 zap[
'f'] =
make_field(inputFile,
"event file(s) or histogram file");
80 zap[
'F'] =
make_field(formula,
"fit formula, e.g: \"[0]+[1]*x\"") =
"";
84 zap[
'O'] =
make_field(option,
"Fit option") =
"";
90 catch(
const exception& error) {
91 FATAL(error.what() << endl);
95 const double NUMBER_OF_ENTRIES = 100.0;
96 const double THRESHOLD = 0.05;
101 if (inputFile.size() == 1) {
103 in = TFile::Open(inputFile[0].c_str(),
"READ");
105 if (in != NULL && in->IsOpen()) {
106 h2 = (TH2D*) in->Get(
"h2");
118 while (inputFile.hasNext()) {
120 STATUS(
"event: " << setw(10) << inputFile.getCounter() <<
'\r');
DEBUG(endl);
122 multi_pointer_type ps = inputFile.next();
129 for (vector<Trk>::const_iterator i = event->mc_trks.begin(); i != event->mc_trks.end(); ++i) {
143 if (evt->begin() == __end) {
147 JEvt::iterator p = min_element(evt->begin(), __end,
quality_sorter);
155 if (ta.
getE() > 0.0 && tb.
getE() > 0.0) {
156 h2->Fill(log10(tb.
getE()), log10(ta.
getE()));
164 const Int_t nx = h2->GetXaxis()->GetNbins();
165 const Double_t xmin = h2->GetXaxis()->GetXmin();
166 const Double_t xmax = h2->GetXaxis()->GetXmax();
167 const Int_t ny = h2->GetYaxis()->GetNbins();
168 const Double_t ymin = h2->GetYaxis()->GetXmin();
169 const Double_t ymax = h2->GetYaxis()->GetXmax();
171 TH1D h1(
"h1", NULL, nx, xmin, xmax);
175 for (Int_t ix = 1; ix <= nx; ++ix) {
184 for (
int iy = 1; iy <= ny; ++iy) {
186 const Double_t x = h2->GetYaxis()->GetBinCenter(iy);
187 const Double_t y = h2->GetBinContent(ix,iy);
188 const Double_t z = h2->GetBinError (ix,iy);
191 h0->SetBinContent(iy, y);
192 h0->SetBinError (iy, z);
203 for (
int iy = 1 ; iy <= ny; ++iy) {
205 const Double_t x = h0->GetXaxis()->GetBinCenter(iy);
206 const Double_t y = h0->GetBinContent(iy);
208 if (y > THRESHOLD * y0) {
219 TF1 f1(
"f1",
"(x < [1]) * [0]*TMath::Gaus(x,[1],[2]) + (x >= [1]) * ([0] * (1.0 - [3])*TMath::Gaus(x,[1],[2]) + [3]*[0])");
221 const double sigma = 0.7 * Q.
getSTDev();
223 f1.SetParameter(0, y0);
224 f1.SetParameter(1, x0);
225 f1.SetParameter(2, sigma);
226 f1.SetParameter(3, 0.05);
227 f1.SetParLimits(3, 0.0, 1.0);
229 TFitResultPtr
result = h0->Fit(&f1,
"SQL",
"same", x0 - 2.5*sigma, x0 + 2.5*sigma);
231 if (
result.Get() == NULL) {
232 FATAL(
"Invalid TFitResultPtr " << h0->GetName() << endl);
237 << setw(3) << ix <<
' '
238 <<
FIXED(6,3) << h2->GetXaxis()->GetBinCenter(ix) <<
' '
239 <<
FIXED(6,3) << x0 <<
" -> "
240 <<
FIXED(6,3) << f1.GetParameter(1) <<
' '
241 <<
FIXED(6,3) << sigma <<
" -> "
242 <<
FIXED(6,3) << f1.GetParameter(2) <<
' '
245 << (
result->IsValid() ?
"" :
"failed") << endl;
249 h1.SetBinContent(ix, f1.GetParameter(1));
250 h1.SetBinError (ix, f1.GetParError (1));
259 f1 =
new TF1(
"f1", formula.c_str());
261 if (f1->IsZombie()) {
262 FATAL(
"Function: " << formula <<
" is zombie." << endl);
267 for (vector<JToken_t>::const_iterator i = startValues.begin(); i != startValues.end(); ++i) {
271 for (vector<JToken_t>::const_iterator i = fixedValues.begin(); i != fixedValues.end(); ++i) {
276 FATAL(error << endl);
279 if (option.find(
'S') == string::npos) { option +=
"S"; }
285 Double_t x[2] = { xmin, xmax };
286 Double_t y[2] = { xmin, xmax };
287 Double_t z[2] = { 0.00, 0.00 };
293 <<
"x >= " <<
FIXED(5,3) << xmax
299 for (Int_t ix = nx; ix >= 1; --ix) {
301 x[0] = h1.GetXaxis()->GetBinCenter(ix);
302 y[0] = h1.GetBinContent(ix);
303 z[0] = h1.GetBinError (ix);
305 if (z[0] == 0.0 || !X(x[0])) {
311 <<
"x >= " <<
FIXED(5,3) << x[0]
313 <<
"x < " <<
FIXED(5,3) << x[1]
316 <<
"(" <<
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]) <<
")";
329 <<
"x < " <<
FIXED(5,3) << x[1]
332 <<
"(" <<
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]) <<
")";
335 string buffer = os.str();
339 for (string::iterator i = buffer.begin(); i != buffer.end(); ++i) {
345 f1 =
new TF1(
"f1", buffer.c_str(), xmin, xmax);
347 h1.GetListOfFunctions()->Add(f1);
353 out <<
JMeta(argc, argv);
355 out << *h2 << H0 << h1;
361 out << *(f1->GetFormula());