#include #include #include #include "TROOT.h" #include "TFile.h" #include "TH1D.h" #include "TH2D.h" #include "TF1.h" #include "TFitResult.h" #include "TMath.h" #include "km3net-dataformat/online/JDAQ.hh" #include "km3net-dataformat/online/JDAQHeader.hh" #include "JLang/JLangToolkit.hh" #include "JPhysics/JConstants.hh" #include "JDetector/JDetector.hh" #include "JDetector/JDetectorToolkit.hh" #include "JDetector/JPMTParametersMap.hh" #include "JDetector/JPMTAnalogueSignalProcessor.hh" #include "JROOT/JRootToolkit.hh" #include "JTools/JRange.hh" #include "JSupport/JMeta.hh" #include "JSupport/JSupport.hh" #include "JSupport/JFileRecorder.hh" #include "JSupport/JSingleFileScanner.hh" #include "JCalibrate/JCalibrateToT.hh" #include "JCalibrate/JFitToT.hh" #include "Jeep/JPrint.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" namespace { /** * Wild card for histogram naming.\n * The wild card will replaced by the module identifier. */ static const char* const WILDCARD = "%"; } /** * \file * * Auxiliary program to fit time-over-threshold distributions. * The fitted parameters are the PMT gain and gain spread. * \author mkarel */ int main(int argc, char **argv) { using namespace std; using namespace JPP; using namespace KM3NETDAQ; typedef JRange JRange_t; typedef JTYPELIST::typelist typelist; string inputFile; JFileRecorder outputFile; string detectorFile; string pmtFile; bool writeFits; double Wmin = 5e3; // Minimal histogram weight in fit range double LeftMargin = 10.0; // Fit range left of time-over-threshold peak double RightMargin = 10.0; // Fit range right of time-over-threshold peak double gradientThreshold = -0.005; // Minimal right-to-left gradient value in time-over-threshold peak-search,\n // normalized with respect to the total number of histogram entries JRange_t gainLimits (FITTOT_GAIN_MIN, FITTOT_GAIN_MAX); // gain limits JRange_t gainSpreadLimits(FITTOT_GAINSPREAD_MIN, FITTOT_GAINSPREAD_MAX); // gainspread limits JRange_t ToTfitRange; string option; string regexp; int debug; try { JProperties properties; properties.insert(gmake_property(Wmin)); properties.insert(gmake_property(LeftMargin)); properties.insert(gmake_property(RightMargin)); properties.insert(gmake_property(gradientThreshold)); properties.insert(gmake_property(gainLimits)); properties.insert(gmake_property(gainSpreadLimits)); JParser<> zap("Auxiliary program to fit time-over-threshold distributions."); zap['f'] = make_field(inputFile, "input file (output from JCalibrateToT)."); zap['o'] = make_field(outputFile, "output file.") = "fit.root"; zap['a'] = make_field(detectorFile, "detector file."); zap['P'] = make_field(pmtFile, "specify PMT file name that can be given to JTriggerEfficiency.") = ""; zap['w'] = make_field(writeFits, "write fit results."); zap['O'] = make_field(option, "ROOT fit options, see TH1::Fit") = ""; zap['@'] = make_field(properties, "fit properties") = JPARSER::initialised(); zap['x'] = make_field(ToTfitRange, "ROOT fit range (time-over-threshold).") = JRange_t(0.0, 35.0); zap['R'] = make_field(regexp, "regular expression for histogram name.") = MAKE_CSTRING(WILDCARD << _2SToT); zap['d'] = make_field(debug, "debug.") = 1; zap(argc, argv); } catch(const exception &error) { FATAL(error.what() << endl); } //---------------------------------------------------------- // load detector file and PMT parameters //---------------------------------------------------------- JDetector detector; try { load(detectorFile, detector); } catch(const JException& error) { FATAL(error); } JPMTParametersMap parametersMap; if (!pmtFile.empty()) { try { parametersMap.load(pmtFile.c_str()); } catch(const JException& error) {} } parametersMap.comment.add(JMeta(argc, argv)); // Set ROOT fit options if (option.find('R') == string::npos) { option += 'R'; } if (option.find('S') == string::npos) { option += 'S'; } if (option.find('Q') == string::npos && debug < JEEP::debug_t) { option += 'Q'; } TFile* in = TFile::Open(inputFile.c_str(), "exist"); if (in == NULL || !in->IsOpen()) { FATAL("File: " << inputFile << " not opened." << endl); } TH2D gain("gain", NULL, 100, 0.0, 2.0, 100, 0.0, 1.0); TH1D chi2("chi2", NULL, 100, 0.0, 5.0); outputFile.open(); if (!outputFile.is_open()) { FATAL("Error opening file " << outputFile << endl); } outputFile.put(JMeta(argc, argv)); //---------------------------------------------------------- // loop over modules in detector file to fit histogram //---------------------------------------------------------- for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) { NOTICE("Module " << setw(10) << module->getID() << ' ' << module->getLocation() << endl); if (module->empty()) { continue; } //---------------------------------------------------------- // get histogram for this module //---------------------------------------------------------- TH2D* h2s = (TH2D*) in->Get(replace(regexp, WILDCARD, MAKE_STRING(module->getID())).c_str()); if (h2s == NULL) { WARNING("No histogram " << module->getID() << "; skip fit." << endl); continue; } const double ny = h2s->GetYaxis()->GetNbins(); const double ymin = h2s->GetYaxis()->GetXmin(); const double ymax = h2s->GetYaxis()->GetXmax(); TH1D chi2_1d (MAKE_CSTRING(module->getID() << ".1chi2"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS-0.5); TH1D gain_1d (MAKE_CSTRING(module->getID() << ".1gain"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS-0.5); TH1D gainspread_1d(MAKE_CSTRING(module->getID() << ".1gainspread"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS-0.5); for (int ix = 1; ix <= h2s->GetXaxis()->GetNbins(); ++ix) { const JPMTIdentifier id(module->getID(), ix-1); const string name = MAKE_STRING(module->getID() << '.' << FILL(2,'0') << id.getPMTAddress() << FITTOT_SUFFIX); DEBUG("Processing histogram " << name << endl); TH1D h1(name.c_str(), NULL, ny, ymin, ymax); h1.Sumw2(); for (int iy = h2s->GetNbinsY(); iy >= 1; --iy) { const double w = h1.GetBinWidth (iy); const double y = h2s->GetBinContent(ix, iy); h1.SetBinContent(iy, y/w); h1.SetBinError (iy, sqrt(y)/w); } JRange_t range(ToTfitRange); const double weight = h1.Integral(h1.FindBin(range.getLowerLimit()), h1.FindBin(range.getUpperLimit())); if (weight <= Wmin) { WARNING("Insufficient histogram entries for PMT " << id << "(" << weight << " < " << Wmin << "); skip fit." << endl); if (writeFits) { outputFile.put(h1); } continue; } //---------------------------------------------------------- // find the mode of the time-over-threshold distribution //---------------------------------------------------------- double max = 0.0; double ToTmode = JDETECTOR::TIME_OVER_THRESHOLD_NS; JPMTParameters& parameters = parametersMap[id]; DEBUG(RIGHT(10) << "ToT" << RIGHT(10) << "Counts" << RIGHT(10) << "gradient" << RIGHT(10) << "threshold" << RIGHT(10) << "Mode" << RIGHT(10) << "Max" << endl); for (int i = h1.GetBinCenter(ny-1); i > h1.GetBinCenter(h1.FindBin(parameters.mean_ns + parameters.sigma_ns)); --i) { const double x = h1.GetBinCenter (i); const double y = h1.GetBinContent(i); const double gradient = ( (h1.GetBinContent(i-1) - h1.GetBinContent(i+1)) / (h1.GetBinCenter (i+1) - h1.GetBinCenter (i-1)) ); DEBUG(FIXED(10,1) << x << FIXED(10,1) << y << FIXED(10,1) << gradient << FIXED(10,1) << -fabs(gradientThreshold) * h1.Integral() << FIXED(10,1) << ToTmode << FIXED(10,1) << max << endl); if (y > max) { ToTmode = x; max = y; } else if (gradient < -fabs(gradientThreshold) * h1.Integral()) { break; } } //---------------------------------------------------------- // set fit range //---------------------------------------------------------- if (!range.is_valid()) { range.setRange(ToTmode - LeftMargin, ToTmode + RightMargin); } range.setRange(std::max(h1.GetXaxis()->GetXmin(), range.getLowerLimit()), std::min(h1.GetXaxis()->GetXmax(), range.getUpperLimit())); DEBUG(LEFT(10) << "Fit-range: " << FIXED(20,1) << range << endl); //---------------------------------------------------------- // perform minimum chi-square fit //---------------------------------------------------------- JFitToT fit(parameters, range); setParLimits(fit, fit.getModelParameter(&JFitToT::JFitToTParameters::gain), gainLimits.getLowerLimit(), gainLimits.getUpperLimit()); setParLimits(fit, fit.getModelParameter(&JFitToT::JFitToTParameters::gainSpread), gainSpreadLimits.getLowerLimit(), gainSpreadLimits.getUpperLimit()); const double initGain = fit.getCPU().getNPE(ToTmode); const double initGainSpread = initGain / 3.0; if (!gainLimits(initGain)) { // Invalid initial gain value WARNING("Invalid initial gain for PMT " << id << "; set default values." << endl); parameters.gain = gainLimits .constrain(initGain); parameters.gainSpread = gainSpreadLimits.constrain(initGainSpread); if (writeFits) { outputFile.put(h1); } continue; } fit.gain = initGain; fit.gainSpread = initGainSpread; const TFitResultPtr result = fit(h1, option); if (result.Get() == NULL) { FATAL("Invalid TFitResultPtr " << h1.GetName() << endl); } if (int(result) < 0 || !result->IsValid()) { // Fit failed WARNING("Fit failure for PMT " << id << "; set initialization values." << endl); h1.GetListOfFunctions()->Clear(); parameters.gain = initGain; parameters.gainSpread = initGainSpread; if (writeFits) { outputFile.put(h1); } continue; } if ((ymin < range.getLowerLimit() || ymax > range.getUpperLimit()) && (option.find("0") == string::npos && option.find("N") == string::npos)) { // add function with full range TF1* f1 = new TF1(MAKE_CSTRING(FITTOT_FNAME << "_fullRange"), &fit, &JFitToT::getValue, ymin, ymax, JFitToT::getNumberOfModelParameters()); f1->SetParameters(fit.GetParameters()); f1->SetLineStyle(kDotted); f1->SetNpx(1000); h1.GetListOfFunctions()->Add(f1); } // store fit results const double reducedChi2 = result->Chi2() / result->Ndf(); const JFitToTParameters values(fit.GetParameters()); const JFitToTParameters errors(fit.GetParErrors()); chi2_1d. SetBinContent(ix, reducedChi2); gain_1d. SetBinContent(ix, values.gain); gain_1d. SetBinError (ix, errors.gain); gainspread_1d.SetBinContent(ix, values.gainSpread); gainspread_1d.SetBinError (ix, errors.gainSpread); gain.Fill(fit.gain, fit.gainSpread); chi2.Fill(TMath::Log10(reducedChi2)); NOTICE("PMT " << id << " chi2 / NDF " << result->Chi2() << ' ' << result->Ndf() << endl); parameters.gain = fit.gain; parameters.gainSpread = fit.gainSpread; if (writeFits) { outputFile.put(h1); } } if (writeFits) { outputFile.put(chi2_1d); outputFile.put(gain_1d); outputFile.put(gainspread_1d); } } for (JSingleFileScanner in(inputFile); in.hasNext(); ) { const JMeta& meta = *in.next(); parametersMap.comment.add(meta); outputFile.put(meta); } for (JRootFileReader in(inputFile.c_str()); in.hasNext(); ) { outputFile.put(*in.next()); } if (pmtFile != "") { parametersMap.store(pmtFile.c_str()); } outputFile.put(gain); outputFile.put(chi2); outputFile.close(); }