From 157bed5a43cd6ef8ce18239de8f4a06c78b147c9 Mon Sep 17 00:00:00 2001 From: mchiu-bnl Date: Thu, 7 May 2026 13:53:21 -0400 Subject: [PATCH 1/6] mods for mbd calib passes to work with waveform fit dsts --- offline/packages/mbd/Makefile.am | 2 + offline/packages/mbd/MbdCalibReco.cc | 800 +++++++++++++++++++++++++++ offline/packages/mbd/MbdCalibReco.h | 77 +++ offline/packages/mbd/MbdEvent.cc | 73 ++- 4 files changed, 942 insertions(+), 10 deletions(-) create mode 100644 offline/packages/mbd/MbdCalibReco.cc create mode 100644 offline/packages/mbd/MbdCalibReco.h diff --git a/offline/packages/mbd/Makefile.am b/offline/packages/mbd/Makefile.am index be588c0adb..af9a9d2598 100644 --- a/offline/packages/mbd/Makefile.am +++ b/offline/packages/mbd/Makefile.am @@ -70,6 +70,7 @@ pkginclude_HEADERS = \ MbdPmtHit.h \ MbdPmtHitV1.h \ MbdPmtSimHitV1.h \ + MbdCalibReco.h \ MbdRawContainer.h \ MbdRawContainerV1.h \ MbdRawContainerV2.h \ @@ -178,6 +179,7 @@ libmbd_io_la_SOURCES = \ MbdSig.cc libmbd_la_SOURCES = \ + MbdCalibReco.cc \ MbdEvent.cc \ MbdCalib.cc \ MbdReco.cc \ diff --git a/offline/packages/mbd/MbdCalibReco.cc b/offline/packages/mbd/MbdCalibReco.cc new file mode 100644 index 0000000000..0d97a0a9d9 --- /dev/null +++ b/offline/packages/mbd/MbdCalibReco.cc @@ -0,0 +1,800 @@ +#include "MbdCalibReco.h" +#include "MbdCalib.h" +#include "MbdDefs.h" +//#include "MbdRawContainer.h" +//#include "MbdRawHit.h" +#include "MbdPmtContainer.h" +#include "MbdPmtHit.h" +#include "MbdOut.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +MbdCalibReco::MbdCalibReco(const std::string &name) + : SubsysReco(name) +{ +} + +int MbdCalibReco::Init(PHCompositeNode * /*topNode*/) +{ + _mbdcal = new MbdCalib(); + _mbdcal->Verbosity( Verbosity() ); + return Fun4AllReturnCodes::EVENT_OK; +} + +int MbdCalibReco::InitRun(PHCompositeNode *topNode) +{ + _runheader = findNode::getClass(topNode, "RunHeader"); + if (!_runheader) + { + std::cout << PHWHERE << " RunHeader node not found, will use run number 0" << std::endl; + } + + _runnumber = _runheader ? _runheader->get_RunNumber() : 0; + + // Build run directory path and create it + std::ostringstream oss; + oss << _caldir << "/" << _runnumber; + _rundir = oss.str(); + gSystem->Exec(("mkdir -p " + _rundir).c_str()); + + if (!_cdbtag.empty()) + { + // Download baseline calibrations from CDB + recoConsts::instance()->set_StringFlag("CDB_GLOBALTAG", _cdbtag); + CDBInterface* cdb = CDBInterface::instance(); + std::string url; + + url = cdb->getUrl("MBD_SAMPMAX"); + if (!url.empty()) { _mbdcal->Download_SampMax(url); } + + url = cdb->getUrl("MBD_PED"); + if (!url.empty()) { _mbdcal->Download_Ped(url); } + + url = cdb->getUrl("MBD_TIMECORR"); + if (!url.empty()) { _mbdcal->Download_TimeCorr(url); } + + url = cdb->getUrl("MBD_SLEWCORR"); + if (!url.empty()) { _mbdcal->Download_SlewCorr(url); } + + std::cout << Name() << ": loaded calibrations from CDB tag " << _cdbtag << std::endl; + } + else + { + // Load baseline calibrations from local files if they exist + std::string calfile = _rundir + "/mbd_sampmax.calib"; + if (gSystem->AccessPathName(calfile.c_str()) == 0) + { + _mbdcal->Download_SampMax(calfile); + } + calfile = _rundir + "/mbd_ped.calib"; + if (gSystem->AccessPathName(calfile.c_str()) == 0) + { + _mbdcal->Download_Ped(calfile); + } + + // Load slew correction for subpass >= 2 + if (_subpass >= 2) + { + calfile = _rundir + "/mbd_slewcorr.calib"; + if (gSystem->AccessPathName(calfile.c_str()) == 0) + { + _mbdcal->Download_SlewCorr(calfile); + std::cout << Name() << ": loaded " << calfile << std::endl; + } + else + { + std::cout << Name() << ": WARNING: " << calfile << " not found" << std::endl; + } + } + } + + // Load t0 offsets for subpass >= 1 (always from local files — outputs of previous subpass) + if (_subpass >= 1) + { + std::string prevpass = "pass" + std::to_string(_subpass - 1) + "_"; + + std::string calfile = _rundir + "/" + prevpass + "mbd_tq_t0.calib"; + if (gSystem->AccessPathName(calfile.c_str()) == 0) + { + _mbdcal->Download_TQT0(calfile); + std::cout << Name() << ": loaded " << calfile << std::endl; + } + else + { + std::cout << Name() << ": WARNING: " << calfile << " not found" << std::endl; + } + + calfile = _rundir + "/" + prevpass + "mbd_tt_t0.calib"; + if (gSystem->AccessPathName(calfile.c_str()) == 0) + { + _mbdcal->Download_TTT0(calfile); + std::cout << Name() << ": loaded " << calfile << std::endl; + } + else + { + std::cout << Name() << ": WARNING: " << calfile << " not found" << std::endl; + } + } + + // Build bitmask of scaled triggers whose names begin with "MBD N&S" + _mbias_trigger_mask = 0xfc00; + + InitHistos(); + + // Open output ROOT file + std::string outfname = _rundir + "/calmbdpass2." + std::to_string(_subpass); + if (_subpass == 0) + { + outfname += "_time-" + std::to_string(_runnumber) + ".root"; + } + else if (_subpass == 1 || _subpass == 2) + { + outfname += "_slew-" + std::to_string(_runnumber) + ".root"; + } + else + { + outfname += "_q-" + std::to_string(_runnumber) + ".root"; + } + _outfile = std::make_unique(outfname.c_str(), "RECREATE"); + std::cout << Name() << ": output file " << outfname << std::endl; + + return Fun4AllReturnCodes::EVENT_OK; +} + +int MbdCalibReco::getNodes(PHCompositeNode *topNode) +{ + _evtheader = findNode::getClass(topNode, "EventHeader"); + if (!_evtheader) + { + std::cout << PHWHERE << " EvtHeader not found, will use run number 0" << std::endl; + } + + _gl1packet = findNode::getClass(topNode,14001); + if (!_gl1packet) + { + _gl1packet = findNode::getClass(topNode, "GL1Packet"); + static int counter = 0; + if ( counter<4 ) + { + std::cout << PHWHERE << " GL1Packet not found" << std::endl; + } + } + + /* + _mbdraws = findNode::getClass(topNode, "MbdRawContainer"); + if (!_mbdraws) + { + static int counter = 0; + if ( counter<4 ) + { + std::cout << PHWHERE << " MbdRawContainer not found" << std::endl; + } + } + */ + + _mbdpmts = findNode::getClass(topNode, "MbdPmtContainer"); + if (!_mbdpmts) + { + static int counter = 0; + if ( counter<4 ) + { + std::cout << PHWHERE << " MbdPmtContainer not found" << std::endl; + } + } + + _mbdout = findNode::getClass(topNode, "MbdOut"); + if (!_mbdout) + { + static int counter = 0; + if ( counter<4 ) + { + std::cout << PHWHERE << " MbdOut not found" << std::endl; + } + } + + _mbdgeom = findNode::getClass(topNode, "MbdGeom"); + if (!_mbdgeom) + { + static int counter = 0; + if ( counter<4 ) + { + std::cout << PHWHERE << " MbdGeom not found" << std::endl; + } + } + + return Fun4AllReturnCodes::EVENT_OK; +} + +void MbdCalibReco::InitHistos() +{ + // Histograms must not be associated with the output TFile at creation + // time (InitRun happens before _outfile is opened above, but we call + // InitHistos before opening the file, so ROOT's current directory is + // gROOT or whichever file is current from the framework). + gROOT->cd(); + + for (int ipmt = 0; ipmt < MbdDefs::MBD_N_PMT; ipmt++) + { + std::string sn = std::to_string(ipmt); + + h_tt[ipmt] = new TH1F(("h_tt" + sn).c_str(), ("tt" + sn).c_str(), 7000, -30., 30.); + h_tt[ipmt]->SetXTitle("ns"); + + h_tq[ipmt] = new TH1F(("h_tq" + sn).c_str(), ("tq" + sn).c_str(), 7000, -150., 31. * 17.7623); + h_tq[ipmt]->SetXTitle("ns"); + + h_qp[ipmt] = new TH1F(("h_q" + sn).c_str(), ("q" + sn).c_str(), 3000, -100., 14900.); + h_qp[ipmt]->SetXTitle("ADC"); + + if (_subpass >= 1) + { + h2_slew[ipmt] = new TH2F(("h2_slew" + sn).c_str(), ("slew curve, ch " + sn).c_str(), 4000, -0.5, 16000. - 0.5, 1100, -5., 6.); + h2_slew[ipmt]->SetXTitle("ADC"); + h2_slew[ipmt]->SetYTitle("#Delta T (ns)"); + } + else + { + h2_slew[ipmt] = nullptr; + } + } + + h2_tt = new TH2F("h2_tt", "ch vs tt", 900, -150., 150., MbdDefs::MBD_N_PMT, -0.5, MbdDefs::MBD_N_PMT - 0.5); + h2_tt->SetXTitle("tt [ns]"); + h2_tt->SetYTitle("pmt ch"); + + h2_tq = new TH2F("h2_tq", "ch vs tq", 900, -150., 150., MbdDefs::MBD_N_PMT, -0.5, MbdDefs::MBD_N_PMT - 0.5); + h2_tq->SetXTitle("tq [ns]"); + h2_tq->SetYTitle("pmt ch"); +} + +int MbdCalibReco::process_event(PHCompositeNode *topNode) +{ + getNodes(topNode); + + // Require a scaled "MBD N&S" trigger + if (_mbias_trigger_mask != 0) + { + uint64_t strig = _gl1packet->getScaledVector(); // scaled trigger only + if ( (strig&_mbias_trigger_mask)==0 ) + { + return Fun4AllReturnCodes::ABORTEVENT; + } + } + + // Per-event arrays for corrected times + /* + std::array ttcorr{}; + std::array tqcorr{}; + std::array adc_arr{}; + ttcorr.fill(std::numeric_limits::quiet_NaN()); + tqcorr.fill(std::numeric_limits::quiet_NaN()); + adc_arr.fill(0); + */ + + std::array armtime{}; + armtime.fill(0); + std::array nhit{}; + nhit.fill(0); + + Float_t zvtx = _mbdout->get_zvtx(); + // Vertex cut for subpass >= 1 + if ( _subpass >= 1 ) + { + if (std::abs(zvtx) > 60.) + { + return Fun4AllReturnCodes::EVENT_OK; + } + } + + for (int iarm=0; iarm<2; iarm++) + { + armtime[iarm] = _mbdout->get_time(iarm); + nhit[iarm] = _mbdout->get_npmt(iarm); + } + + for (int ipmt=0; ipmt < _mbdpmts->get_npmt(); ipmt++) + { + MbdPmtHit *pmt = _mbdpmts->get_pmt(ipmt); + if ( !pmt ) + { + continue; + } + + Short_t pmtno = pmt->get_pmt(); + Float_t q = pmt->get_q(); + Float_t tt = pmt->get_tt(); + Float_t tq = pmt->get_tq(); + + h_tt[pmtno]->Fill( tt ); + h2_tt->Fill( tt, pmtno ); + h_tq[pmtno]->Fill( tq ); + h2_tq->Fill( tq, pmtno ); + + // Fill charge histogram for in-time hits + if ( std::abs(tt)<26.0 && q > 0.) + { + h_qp[pmtno]->Fill( q ); + } + + int arm = _mbdgeom->get_arm( pmtno ); + + // Fill slew histogram for subpass >= 1 + if (_subpass >= 1 && h2_slew[pmtno]) + { + if (nhit[arm] >= 2. && q > 0.) + { + float dt = tt - armtime[arm]; + h2_slew[pmtno]->Fill(q, dt); + } + } + } + + return Fun4AllReturnCodes::EVENT_OK; +} + +int MbdCalibReco::End(PHCompositeNode * /*topNode*/) +{ + if (!_outfile) + { + return Fun4AllReturnCodes::EVENT_OK; + } + + // Write histograms to output file + _outfile->cd(); + h2_tt->Write(); + h2_tq->Write(); + for (int ipmt = 0; ipmt < MbdDefs::MBD_N_PMT; ipmt++) + { + h_tt[ipmt]->Write(); + h_tq[ipmt]->Write(); + h_qp[ipmt]->Write(); + if (h2_slew[ipmt]) + { + h2_slew[ipmt]->Write(); + } + } + + // Always fit and write t0 (done at every subpass from the accumulated histograms) + //FitAndWriteT0(); + + /* + if (_subpass == 1 || _subpass == 2) + { + FitAndWriteSlew(); + } + */ + + _outfile->Close(); + + return Fun4AllReturnCodes::EVENT_OK; +} + +int MbdCalibReco::getRunType() const +{ + // Run number → collision system (mirrors get_runtype() from get_runstr.h) + if (_runnumber <= 30000) + { + return 3; // SIMAUAU200 + } + if (_runnumber <= 53880) + { + return 1; // PP200 (Run24) + } + if (_runnumber <= 54962) + { + return 0; // AUAU200 (Run24) + } + if (_runnumber <= 78954) + { + return 0; // AUAU200 (Run25) + } + if (_runnumber <= 81667) + { + return 1; // PP200 (Run25) + } + if (_runnumber <= 82703) + { + return 2; // OO200 (Run25) + } + return -1; +} + +// --------------------------------------------------------------------------- +// FitAndWriteT0 — Gaussian fit to h_tt and h_tq, write *_t0.calib files +// --------------------------------------------------------------------------- +void MbdCalibReco::FitAndWriteT0() +{ + std::string passprefix = "pass" + std::to_string(_subpass) + "_"; + std::string tt_fname = _rundir + "/" + passprefix + "mbd_tt_t0.calib"; + std::string tq_fname = _rundir + "/" + passprefix + "mbd_tq_t0.calib"; + + std::ofstream tt_file(tt_fname); + std::ofstream tq_file(tq_fname); + if (!tt_file.is_open() || !tq_file.is_open()) + { + std::cout << Name() << "::FitAndWriteT0 ERROR cannot open calib files" << std::endl; + return; + } + + TF1 gaussian("mbdcal_gaus", "gaus", -25., 25.); + gaussian.SetLineColor(2); + + double min_twindow = -25.; + double max_twindow = 25.; + + // --- tt_t0 --- + for (int ipmt = 0; ipmt < MbdDefs::MBD_N_PMT; ipmt++) + { + if (ipmt == 0 || ipmt == 64) + { + h_tt[ipmt]->SetAxisRange(-25., 25.); + } + else + { + h_tt[ipmt]->SetAxisRange(min_twindow, max_twindow); + } + + int peakbin = h_tt[ipmt]->GetMaximumBin(); + double mean = h_tt[ipmt]->GetBinCenter(peakbin); + double peak = h_tt[ipmt]->GetMaximum(); + + gaussian.SetParameters(peak, mean, 5.); + gaussian.SetRange(mean - 3., mean + 3.); + h_tt[ipmt]->Fit(&gaussian, "RQ"); + + mean = gaussian.GetParameter(1); + double meanerr = gaussian.GetParError(1); + double sigma = gaussian.GetParameter(2); + double sigmaerr = gaussian.GetParError(2); + + if (ipmt == 0 || ipmt == 64) + { + min_twindow = mean - 3. * sigma; + max_twindow = mean + 3. * sigma; + } + + tt_file << ipmt << "\t" << mean << "\t" << meanerr << "\t" + << sigma << "\t" << sigmaerr << "\n"; + + // Normalise h2_tt row by fit peak amplitude + double fitpeak = gaussian.GetParameter(0); + if (fitpeak != 0.) + { + int nbinsx = h2_tt->GetNbinsX(); + for (int ibinx = 1; ibinx <= nbinsx; ibinx++) + { + float bc = h2_tt->GetBinContent(ibinx, ipmt + 1); + h2_tt->SetBinContent(ibinx, ipmt + 1, bc / fitpeak); + } + } + } + tt_file.close(); + + // Write canonical CDB ROOT file + { + MbdCalib tmpcal; + tmpcal.Download_TTT0(tt_fname); + std::string cdb_fname = tt_fname; + cdb_fname.replace(cdb_fname.rfind(".calib"), 6, ".root"); + tmpcal.Write_CDB_TTT0(cdb_fname); + } + + // --- tq_t0 --- + min_twindow = -25.; + max_twindow = 25.; + + for (int ipmt = 0; ipmt < MbdDefs::MBD_N_PMT; ipmt++) + { + if (ipmt == 0 || ipmt == 64) + { + h_tq[ipmt]->SetAxisRange(-25., 25.); + } + else + { + h_tq[ipmt]->SetAxisRange(min_twindow, max_twindow); + } + + int peakbin = h_tq[ipmt]->GetMaximumBin(); + double mean = h_tq[ipmt]->GetBinCenter(peakbin); + double peak = h_tq[ipmt]->GetMaximum(); + + gaussian.SetParameters(peak, mean, 5.); + gaussian.SetRange(mean - 3., mean + 3.); + h_tq[ipmt]->Fit(&gaussian, "RQ"); + + mean = gaussian.GetParameter(1); + double meanerr = gaussian.GetParError(1); + double sigma = gaussian.GetParameter(2); + double sigmaerr = gaussian.GetParError(2); + + if (ipmt == 0 || ipmt == 64) + { + min_twindow = mean - 3. * sigma; + max_twindow = mean + 3. * sigma; + } + + tq_file << ipmt << "\t" << mean << "\t" << meanerr << "\t" + << sigma << "\t" << sigmaerr << "\n"; + + // Normalise h2_tq row + double fitpeak = gaussian.GetParameter(0); + if (fitpeak != 0.) + { + int nbinsx = h2_tq->GetNbinsX(); + for (int ibinx = 1; ibinx <= nbinsx; ibinx++) + { + float bc = h2_tq->GetBinContent(ibinx, ipmt + 1); + h2_tq->SetBinContent(ibinx, ipmt + 1, bc / fitpeak); + } + } + } + tq_file.close(); + + { + MbdCalib tmpcal; + tmpcal.Download_TQT0(tq_fname); + std::string cdb_fname = tq_fname; + cdb_fname.replace(cdb_fname.rfind(".calib"), 6, ".root"); + tmpcal.Write_CDB_TQT0(cdb_fname); + } + + std::cout << Name() << ": wrote " << tt_fname << " and " << tq_fname << std::endl; +} + +// --------------------------------------------------------------------------- +// FindTH2Ridge — column-by-column Gaussian fits to find slew-correction ridge +// --------------------------------------------------------------------------- +void MbdCalibReco::FindTH2Ridge(const TH2 *h2, TGraphErrors *&gridge, + TGraphErrors *&grms) const +{ + int nbinsx = h2->GetNbinsX(); + double min_yrange = h2->GetYaxis()->GetBinLowEdge(1); + double max_yrange = h2->GetYaxis()->GetBinLowEdge(h2->GetNbinsY() + 1); + + gridge = new TGraphErrors(); + gridge->SetName("gridge"); + gridge->SetTitle("ridge"); + grms = new TGraphErrors(); + grms->SetName("grms"); + grms->SetTitle("rms of ridge"); + + TH1 *h_projx = h2->ProjectionX("_projx_tmp"); + TF1 gaussian("_slew_gaus", "gaus", min_yrange, max_yrange); + gaussian.SetLineColor(4); + + TH1 *h_projy = nullptr; + double adcmean = 0.; + double adcnum = 0.; + + for (int ibin = 1; ibin <= nbinsx; ibin++) + { + std::string projname = "_hproj_" + std::to_string(ibin); + if (!h_projy) + { + h_projy = h2->ProjectionY(projname.c_str(), ibin, ibin); + adcmean = h_projx->GetBinCenter(ibin); + adcnum = 1.; + } + else + { + TH1 *hadd = h2->ProjectionY(projname.c_str(), ibin, ibin); + h_projy->Add(hadd); + delete hadd; + adcmean += h_projx->GetBinCenter(ibin); + adcnum += 1.; + } + + if (h_projy->Integral() > 2000. || ibin == nbinsx) + { + adcmean /= adcnum; + + int maxbin = h_projy->GetMaximumBin(); + double xmax_g = h_projy->GetBinCenter(maxbin); + double ymax_g = h_projy->GetBinContent(maxbin); + gaussian.SetParameter(0, ymax_g); + gaussian.SetParameter(1, xmax_g); + gaussian.SetRange(xmax_g - 0.6, xmax_g + 0.6); + h_projy->Fit(&gaussian, "RWWQ"); + + double mean = gaussian.GetParameter(1); + double meanerr = gaussian.GetParError(1); + double rms = gaussian.GetParameter(2); + double rmserr = gaussian.GetParError(2); + + if (meanerr < 1.0) + { + int n = gridge->GetN(); + gridge->SetPoint(n, adcmean, mean); + gridge->SetPointError(n, 0., meanerr); + } + if (rmserr < 0.01) + { + int n = grms->GetN(); + grms->SetPoint(n, adcmean, rms); + grms->SetPointError(n, 0., rmserr); + } + + delete h_projy; + h_projy = nullptr; + adcmean = 0.; + adcnum = 0.; + } + } + + gridge->SetBit(TGraph::kIsSortedX); + grms->SetBit(TGraph::kIsSortedX); + delete h_projx; +} + +// --------------------------------------------------------------------------- +// FitAndWriteSlew — build slew-correction LUT from h2_slew ridge +// --------------------------------------------------------------------------- +void MbdCalibReco::FitAndWriteSlew() +{ + const int NPOINTS = 16000; + const int MINADC = 0; + const int MAXADC = 15999; + + std::string scorr_fname = _rundir + "/mbd_slewcorr.calib"; + std::ofstream scorr_file(scorr_fname); + if (!scorr_file.is_open()) + { + std::cout << Name() << "::FitAndWriteSlew ERROR cannot open " << scorr_fname << std::endl; + return; + } + + std::string trms_fname = _rundir + "/mbd_timerms.calib"; + std::ofstream trms_file(trms_fname); + + // Arrays of slew/trms graph pointers, indexed by feech (only T-channels used) + std::array g_slew{}; + std::array g_trms{}; + g_slew.fill(nullptr); + g_trms.fill(nullptr); + + for (int ipmt = 0; ipmt < MbdDefs::MBD_N_PMT; ipmt++) + { + if (!h2_slew[ipmt]) + { + continue; + } + + int feech_t = (ipmt / 8) * 16 + ipmt % 8; + + TGraphErrors *gr = nullptr; + TGraphErrors *grms_tmp = nullptr; + FindTH2Ridge(h2_slew[ipmt], gr, grms_tmp); + + g_slew[feech_t] = gr; + g_trms[feech_t] = grms_tmp; + + if (gr) + { + gr->SetName(("g_slew" + std::to_string(ipmt)).c_str()); + gr->SetMarkerStyle(20); + gr->SetMarkerSize(0.25); + } + if (grms_tmp) + { + grms_tmp->SetName(("g_trms" + std::to_string(ipmt)).c_str()); + grms_tmp->SetMarkerStyle(20); + grms_tmp->SetMarkerSize(0.25); + } + } + + // Write slew correction LUT (one T-channel feech at a time) + for (int ifeech = 0; ifeech < MbdDefs::MBD_N_FEECH; ifeech++) + { + // Only T-channels: type = (feech/8) % 2 == 0 + if ((ifeech / 8) % 2 == 1) + { + continue; + } + + if (!g_slew[ifeech]) + { + continue; + } + + scorr_file << ifeech << "\t" << NPOINTS << "\t" << MINADC << "\t" << MAXADC << "\n"; + int step = (MAXADC - MINADC) / (NPOINTS - 1); + for (int iadc = MINADC; iadc <= MAXADC; iadc += step) + { + scorr_file << g_slew[ifeech]->Eval(iadc) << " "; + if (iadc % 10 == 9) + { + scorr_file << "\n"; + } + } + } + scorr_file.close(); + + // Write time-RMS LUT + if (trms_file.is_open()) + { + for (int ifeech = 0; ifeech < MbdDefs::MBD_N_FEECH; ifeech++) + { + if ((ifeech / 8) % 2 == 1) + { + continue; + } + if (!g_trms[ifeech]) + { + continue; + } + + trms_file << ifeech << "\t" << NPOINTS << "\t" << MINADC << "\t" << MAXADC << "\n"; + int step = (MAXADC - MINADC) / (NPOINTS - 1); + for (int iadc = MINADC; iadc <= MAXADC; iadc += step) + { + trms_file << g_trms[ifeech]->Eval(iadc) << " "; + if (iadc % 10 == 9) + { + trms_file << "\n"; + } + } + } + trms_file.close(); + } + + // Write graphs to ROOT file and create CDB ROOT file + _outfile->cd(); + for (auto *g : g_slew) + { + if (g) + { + g->Write(); + } + } + for (auto *g : g_trms) + { + if (g) + { + g->Write(); + } + } + + { + MbdCalib tmpcal; + tmpcal.Download_SlewCorr(scorr_fname); + std::string cdb_fname = scorr_fname; + cdb_fname.replace(cdb_fname.rfind(".calib"), 6, ".root"); + tmpcal.Write_CDB_SlewCorr(cdb_fname); + } + + // Clean up + for (auto *g : g_slew) + { + delete g; + } + for (auto *g : g_trms) + { + delete g; + } + + std::cout << Name() << ": wrote " << scorr_fname << std::endl; +} + diff --git a/offline/packages/mbd/MbdCalibReco.h b/offline/packages/mbd/MbdCalibReco.h new file mode 100644 index 0000000000..f6d6448932 --- /dev/null +++ b/offline/packages/mbd/MbdCalibReco.h @@ -0,0 +1,77 @@ +#ifndef MBD_MBDCALIBRECO_H +#define MBD_MBDCALIBRECO_H + +#include "MbdDefs.h" + +#include + +#include +#include +#include +#include + +class PHCompositeNode; +class MbdCalib; +class MbdPmtContainer; +class MbdOut; +class MbdGeom; +class Gl1Packet; +class EventHeader; +class RunHeader; +class TH1; +class TH2; +class TFile; +class TGraphErrors; + +class MbdCalibReco : public SubsysReco +{ + public: + MbdCalibReco(const std::string& name = "MbdCalibReco"); + ~MbdCalibReco() override = default; + + int Init(PHCompositeNode* topNode) override; + int InitRun(PHCompositeNode* topNode) override; + int process_event(PHCompositeNode* topNode) override; + int End(PHCompositeNode* topNode) override; + + void SetSubPass(const int s) { _subpass = s; } + void SetCalDir(const std::string& d) { _caldir = d; } + void SetCDBTag(const std::string& t) { _cdbtag = t; } + + private: + int getNodes(PHCompositeNode* topNode); + void InitHistos(); + int getRunType() const; + + void FitAndWriteT0(); + void FitAndWriteSlew(); + + void FindTH2Ridge(const TH2* h2, TGraphErrors*& gridge, TGraphErrors*& grms) const; + + uint64_t _mbias_trigger_mask{0}; + + int _subpass{0}; + int _runnumber{0}; + std::string _caldir{"results"}; + std::string _rundir; // _caldir// + std::string _cdbtag{}; // non-empty → download from CDB instead of local files + + MbdCalib* _mbdcal{nullptr}; + MbdPmtContainer* _mbdpmts{nullptr}; + MbdOut* _mbdout{nullptr}; + MbdGeom* _mbdgeom{nullptr}; + EventHeader* _evtheader{nullptr}; + RunHeader* _runheader{nullptr}; + Gl1Packet* _gl1packet{nullptr}; + + std::array h_tt{}; + std::array h_tq{}; + std::array h_qp{}; + std::array h2_slew{}; + TH2* h2_tt{nullptr}; + TH2* h2_tq{nullptr}; + + std::unique_ptr _outfile{nullptr}; +}; + +#endif // MBD_MBDCALIBRECO_H diff --git a/offline/packages/mbd/MbdEvent.cc b/offline/packages/mbd/MbdEvent.cc index b97d750f29..5c1e9ed087 100644 --- a/offline/packages/mbd/MbdEvent.cc +++ b/offline/packages/mbd/MbdEvent.cc @@ -33,6 +33,7 @@ #include #include #include +#include #include #include #include @@ -165,18 +166,32 @@ int MbdEvent::InitRun() if ( _calpass>1 ) { std::string calfname = "results/"; calfname += std::to_string(_runnum); calfname += "/mbd_sampmax.calib"; - std::cout << "Loading local sampmax, " << calfname << std::endl; - _mbdcal->Download_SampMax( calfname ); + if ( std::filesystem::exists(calfname) ) + { + std::cout << "Loading local sampmax, " << calfname << std::endl; + _mbdcal->Download_SampMax( calfname ); + } + else + { + std::cout << PHWHERE << "local sampmax not found, skipping: " << calfname << std::endl; + } calfname = "results/"; calfname += std::to_string(_runnum); calfname += "/mbd_ped.calib"; - std::cout << "Loading local ped, " << calfname << std::endl; - _mbdcal->Download_Ped( calfname ); + if ( std::filesystem::exists(calfname) ) + { + std::cout << "Loading local ped, " << calfname << std::endl; + _mbdcal->Download_Ped( calfname ); + } + else + { + std::cout << PHWHERE << "local ped not found, skipping: " << calfname << std::endl; + } } // check if sampmax and ped calibs exist int scheck = _mbdcal->get_sampmax(0); - if ( (scheck<0 || _is_online) && _calpass!=1 ) + if ( (scheck<0 || _is_online) && _calpass==0 ) { _no_sampmax = 1000; // num events for on the fly calculation _calib_done = 0; @@ -281,11 +296,49 @@ int MbdEvent::InitRun() if ( _calpass == 2 ) { - // zero out the tt_t0, tq_t0, and gains to produce uncalibrated time and charge std::cout << "MBD Cal Pass 2" << std::endl; - _mbdcal->Reset_TTT0(); - _mbdcal->Reset_TQT0(); - _mbdcal->Reset_Gains(); + + // zero out the tt_t0, tq_t0, and gains to produce uncalibrated time and charge + // or load pass2 calibs from local file for calpass2+, if local files exist + if ( _calpass>1 ) + { + std::string calfname = "results/"; calfname += std::to_string(_runnum); calfname += "/mbd_tt_t0.calib"; + if ( std::filesystem::exists(calfname) ) + { + std::cout << "Loading local mbd_tt_t0, " << calfname << std::endl; + _mbdcal->Download_TTT0( calfname ); + } + else + { + _mbdcal->Reset_TTT0(); + std::cout << PHWHERE << "local mbd_tt_t0 not found, reset to 0: " << calfname << std::endl; + } + + calfname = "results/"; calfname += std::to_string(_runnum); calfname += "/mbd_tq_t0.calib"; + if ( std::filesystem::exists(calfname) ) + { + std::cout << "Loading local mbd_tq_t0, " << calfname << std::endl; + _mbdcal->Download_TQT0( calfname ); + } + else + { + _mbdcal->Reset_TQT0(); + std::cout << PHWHERE << "local mbd_tq_t0 not found, reset to 0: " << calfname << std::endl; + } + + calfname = "results/"; calfname += std::to_string(_runnum); calfname += "/mbd_qfit.calib"; + if ( std::filesystem::exists(calfname) ) + { + std::cout << "Loading local mbd_qfit, " << calfname << std::endl; + _mbdcal->Download_Gains( calfname ); + } + else + { + _mbdcal->Reset_Gains(); + std::cout << PHWHERE << "local mbd_gains not found, reset to 1: " << calfname << std::endl; + } + + } TDirectory *orig_dir = gDirectory; @@ -1381,7 +1434,7 @@ int MbdEvent::FillSampMaxCalib() // _no_sampmax keeps track of how many events to use for on-the-fly calibration _no_sampmax--; - if ( _no_sampmax==0 && _calpass != 1 ) + if ( _no_sampmax==0 && _calpass==0 ) { CalcSampMaxCalib(); _calib_done = 1; From 3a6becdf82fb2b9aa2bdcacf780d8479081c082b Mon Sep 17 00:00:00 2001 From: mchiu-bnl Date: Thu, 7 May 2026 15:43:33 -0400 Subject: [PATCH 2/6] rabbit fixes --- offline/packages/mbd/MbdCalibReco.cc | 488 ++------------------------- offline/packages/mbd/MbdCalibReco.h | 13 +- 2 files changed, 35 insertions(+), 466 deletions(-) diff --git a/offline/packages/mbd/MbdCalibReco.cc b/offline/packages/mbd/MbdCalibReco.cc index 0d97a0a9d9..10ed14b559 100644 --- a/offline/packages/mbd/MbdCalibReco.cc +++ b/offline/packages/mbd/MbdCalibReco.cc @@ -1,8 +1,6 @@ #include "MbdCalibReco.h" #include "MbdCalib.h" #include "MbdDefs.h" -//#include "MbdRawContainer.h" -//#include "MbdRawHit.h" #include "MbdPmtContainer.h" #include "MbdPmtHit.h" #include "MbdOut.h" @@ -24,6 +22,7 @@ #include #include #include +#include #include #include @@ -39,7 +38,7 @@ MbdCalibReco::MbdCalibReco(const std::string &name) int MbdCalibReco::Init(PHCompositeNode * /*topNode*/) { - _mbdcal = new MbdCalib(); + _mbdcal = std::make_unique(); _mbdcal->Verbosity( Verbosity() ); return Fun4AllReturnCodes::EVENT_OK; } @@ -54,6 +53,8 @@ int MbdCalibReco::InitRun(PHCompositeNode *topNode) _runnumber = _runheader ? _runheader->get_RunNumber() : 0; + getNodes(topNode); + // Build run directory path and create it std::ostringstream oss; oss << _caldir << "/" << _runnumber; @@ -142,9 +143,9 @@ int MbdCalibReco::InitRun(PHCompositeNode *topNode) // Build bitmask of scaled triggers whose names begin with "MBD N&S" _mbias_trigger_mask = 0xfc00; - InitHistos(); - // Open output ROOT file + TDirectory *origdir = gDirectory; + std::string outfname = _rundir + "/calmbdpass2." + std::to_string(_subpass); if (_subpass == 0) { @@ -159,8 +160,18 @@ int MbdCalibReco::InitRun(PHCompositeNode *topNode) outfname += "_q-" + std::to_string(_runnumber) + ".root"; } _outfile = std::make_unique(outfname.c_str(), "RECREATE"); + if (!_outfile || _outfile->IsZombie()) + { + std::cerr << PHWHERE << " ERROR: cannot open output file " << outfname << std::endl; + _outfile.reset(); + return Fun4AllReturnCodes::ABORTRUN; + } std::cout << Name() << ": output file " << outfname << std::endl; + InitHistos(); + + origdir->cd(); + return Fun4AllReturnCodes::EVENT_OK; } @@ -183,18 +194,6 @@ int MbdCalibReco::getNodes(PHCompositeNode *topNode) } } - /* - _mbdraws = findNode::getClass(topNode, "MbdRawContainer"); - if (!_mbdraws) - { - static int counter = 0; - if ( counter<4 ) - { - std::cout << PHWHERE << " MbdRawContainer not found" << std::endl; - } - } - */ - _mbdpmts = findNode::getClass(topNode, "MbdPmtContainer"); if (!_mbdpmts) { @@ -225,17 +224,16 @@ int MbdCalibReco::getNodes(PHCompositeNode *topNode) } } + if ( !_mbdgeom || !_mbdout || !_mbdpmts || !_gl1packet || !_evtheader ) + { + return Fun4AllReturnCodes::ABORTRUN; + } + return Fun4AllReturnCodes::EVENT_OK; } void MbdCalibReco::InitHistos() { - // Histograms must not be associated with the output TFile at creation - // time (InitRun happens before _outfile is opened above, but we call - // InitHistos before opening the file, so ROOT's current directory is - // gROOT or whichever file is current from the framework). - gROOT->cd(); - for (int ipmt = 0; ipmt < MbdDefs::MBD_N_PMT; ipmt++) { std::string sn = std::to_string(ipmt); @@ -251,9 +249,12 @@ void MbdCalibReco::InitHistos() if (_subpass >= 1) { - h2_slew[ipmt] = new TH2F(("h2_slew" + sn).c_str(), ("slew curve, ch " + sn).c_str(), 4000, -0.5, 16000. - 0.5, 1100, -5., 6.); - h2_slew[ipmt]->SetXTitle("ADC"); - h2_slew[ipmt]->SetYTitle("#Delta T (ns)"); + const int nbins[2] = {4000, 1100}; + const double xmin[2] = {-0.5, -5.}; + const double xmax[2] = {16000. - 0.5, 6.}; + h2_slew[ipmt] = new THnSparseF(("h2_slew" + sn).c_str(), ("slew curve, ch " + sn).c_str(), 2, nbins, xmin, xmax); + h2_slew[ipmt]->GetAxis(0)->SetTitle("ADC"); + h2_slew[ipmt]->GetAxis(1)->SetTitle("#Delta T (ns)"); } else { @@ -270,10 +271,8 @@ void MbdCalibReco::InitHistos() h2_tq->SetYTitle("pmt ch"); } -int MbdCalibReco::process_event(PHCompositeNode *topNode) +int MbdCalibReco::process_event(PHCompositeNode * /*topNode*/) { - getNodes(topNode); - // Require a scaled "MBD N&S" trigger if (_mbias_trigger_mask != 0) { @@ -284,16 +283,6 @@ int MbdCalibReco::process_event(PHCompositeNode *topNode) } } - // Per-event arrays for corrected times - /* - std::array ttcorr{}; - std::array tqcorr{}; - std::array adc_arr{}; - ttcorr.fill(std::numeric_limits::quiet_NaN()); - tqcorr.fill(std::numeric_limits::quiet_NaN()); - adc_arr.fill(0); - */ - std::array armtime{}; armtime.fill(0); std::array nhit{}; @@ -347,7 +336,8 @@ int MbdCalibReco::process_event(PHCompositeNode *topNode) if (nhit[arm] >= 2. && q > 0.) { float dt = tt - armtime[arm]; - h2_slew[pmtno]->Fill(q, dt); + const double coords[2] = {q, dt}; + h2_slew[pmtno]->Fill(coords); } } } @@ -355,7 +345,7 @@ int MbdCalibReco::process_event(PHCompositeNode *topNode) return Fun4AllReturnCodes::EVENT_OK; } -int MbdCalibReco::End(PHCompositeNode * /*topNode*/) +int MbdCalibReco::EndRun(const int /*runnumber*/) { if (!_outfile) { @@ -377,424 +367,8 @@ int MbdCalibReco::End(PHCompositeNode * /*topNode*/) } } - // Always fit and write t0 (done at every subpass from the accumulated histograms) - //FitAndWriteT0(); - - /* - if (_subpass == 1 || _subpass == 2) - { - FitAndWriteSlew(); - } - */ - _outfile->Close(); return Fun4AllReturnCodes::EVENT_OK; } -int MbdCalibReco::getRunType() const -{ - // Run number → collision system (mirrors get_runtype() from get_runstr.h) - if (_runnumber <= 30000) - { - return 3; // SIMAUAU200 - } - if (_runnumber <= 53880) - { - return 1; // PP200 (Run24) - } - if (_runnumber <= 54962) - { - return 0; // AUAU200 (Run24) - } - if (_runnumber <= 78954) - { - return 0; // AUAU200 (Run25) - } - if (_runnumber <= 81667) - { - return 1; // PP200 (Run25) - } - if (_runnumber <= 82703) - { - return 2; // OO200 (Run25) - } - return -1; -} - -// --------------------------------------------------------------------------- -// FitAndWriteT0 — Gaussian fit to h_tt and h_tq, write *_t0.calib files -// --------------------------------------------------------------------------- -void MbdCalibReco::FitAndWriteT0() -{ - std::string passprefix = "pass" + std::to_string(_subpass) + "_"; - std::string tt_fname = _rundir + "/" + passprefix + "mbd_tt_t0.calib"; - std::string tq_fname = _rundir + "/" + passprefix + "mbd_tq_t0.calib"; - - std::ofstream tt_file(tt_fname); - std::ofstream tq_file(tq_fname); - if (!tt_file.is_open() || !tq_file.is_open()) - { - std::cout << Name() << "::FitAndWriteT0 ERROR cannot open calib files" << std::endl; - return; - } - - TF1 gaussian("mbdcal_gaus", "gaus", -25., 25.); - gaussian.SetLineColor(2); - - double min_twindow = -25.; - double max_twindow = 25.; - - // --- tt_t0 --- - for (int ipmt = 0; ipmt < MbdDefs::MBD_N_PMT; ipmt++) - { - if (ipmt == 0 || ipmt == 64) - { - h_tt[ipmt]->SetAxisRange(-25., 25.); - } - else - { - h_tt[ipmt]->SetAxisRange(min_twindow, max_twindow); - } - - int peakbin = h_tt[ipmt]->GetMaximumBin(); - double mean = h_tt[ipmt]->GetBinCenter(peakbin); - double peak = h_tt[ipmt]->GetMaximum(); - - gaussian.SetParameters(peak, mean, 5.); - gaussian.SetRange(mean - 3., mean + 3.); - h_tt[ipmt]->Fit(&gaussian, "RQ"); - - mean = gaussian.GetParameter(1); - double meanerr = gaussian.GetParError(1); - double sigma = gaussian.GetParameter(2); - double sigmaerr = gaussian.GetParError(2); - - if (ipmt == 0 || ipmt == 64) - { - min_twindow = mean - 3. * sigma; - max_twindow = mean + 3. * sigma; - } - - tt_file << ipmt << "\t" << mean << "\t" << meanerr << "\t" - << sigma << "\t" << sigmaerr << "\n"; - - // Normalise h2_tt row by fit peak amplitude - double fitpeak = gaussian.GetParameter(0); - if (fitpeak != 0.) - { - int nbinsx = h2_tt->GetNbinsX(); - for (int ibinx = 1; ibinx <= nbinsx; ibinx++) - { - float bc = h2_tt->GetBinContent(ibinx, ipmt + 1); - h2_tt->SetBinContent(ibinx, ipmt + 1, bc / fitpeak); - } - } - } - tt_file.close(); - - // Write canonical CDB ROOT file - { - MbdCalib tmpcal; - tmpcal.Download_TTT0(tt_fname); - std::string cdb_fname = tt_fname; - cdb_fname.replace(cdb_fname.rfind(".calib"), 6, ".root"); - tmpcal.Write_CDB_TTT0(cdb_fname); - } - - // --- tq_t0 --- - min_twindow = -25.; - max_twindow = 25.; - - for (int ipmt = 0; ipmt < MbdDefs::MBD_N_PMT; ipmt++) - { - if (ipmt == 0 || ipmt == 64) - { - h_tq[ipmt]->SetAxisRange(-25., 25.); - } - else - { - h_tq[ipmt]->SetAxisRange(min_twindow, max_twindow); - } - - int peakbin = h_tq[ipmt]->GetMaximumBin(); - double mean = h_tq[ipmt]->GetBinCenter(peakbin); - double peak = h_tq[ipmt]->GetMaximum(); - - gaussian.SetParameters(peak, mean, 5.); - gaussian.SetRange(mean - 3., mean + 3.); - h_tq[ipmt]->Fit(&gaussian, "RQ"); - - mean = gaussian.GetParameter(1); - double meanerr = gaussian.GetParError(1); - double sigma = gaussian.GetParameter(2); - double sigmaerr = gaussian.GetParError(2); - - if (ipmt == 0 || ipmt == 64) - { - min_twindow = mean - 3. * sigma; - max_twindow = mean + 3. * sigma; - } - - tq_file << ipmt << "\t" << mean << "\t" << meanerr << "\t" - << sigma << "\t" << sigmaerr << "\n"; - - // Normalise h2_tq row - double fitpeak = gaussian.GetParameter(0); - if (fitpeak != 0.) - { - int nbinsx = h2_tq->GetNbinsX(); - for (int ibinx = 1; ibinx <= nbinsx; ibinx++) - { - float bc = h2_tq->GetBinContent(ibinx, ipmt + 1); - h2_tq->SetBinContent(ibinx, ipmt + 1, bc / fitpeak); - } - } - } - tq_file.close(); - - { - MbdCalib tmpcal; - tmpcal.Download_TQT0(tq_fname); - std::string cdb_fname = tq_fname; - cdb_fname.replace(cdb_fname.rfind(".calib"), 6, ".root"); - tmpcal.Write_CDB_TQT0(cdb_fname); - } - - std::cout << Name() << ": wrote " << tt_fname << " and " << tq_fname << std::endl; -} - -// --------------------------------------------------------------------------- -// FindTH2Ridge — column-by-column Gaussian fits to find slew-correction ridge -// --------------------------------------------------------------------------- -void MbdCalibReco::FindTH2Ridge(const TH2 *h2, TGraphErrors *&gridge, - TGraphErrors *&grms) const -{ - int nbinsx = h2->GetNbinsX(); - double min_yrange = h2->GetYaxis()->GetBinLowEdge(1); - double max_yrange = h2->GetYaxis()->GetBinLowEdge(h2->GetNbinsY() + 1); - - gridge = new TGraphErrors(); - gridge->SetName("gridge"); - gridge->SetTitle("ridge"); - grms = new TGraphErrors(); - grms->SetName("grms"); - grms->SetTitle("rms of ridge"); - - TH1 *h_projx = h2->ProjectionX("_projx_tmp"); - TF1 gaussian("_slew_gaus", "gaus", min_yrange, max_yrange); - gaussian.SetLineColor(4); - - TH1 *h_projy = nullptr; - double adcmean = 0.; - double adcnum = 0.; - - for (int ibin = 1; ibin <= nbinsx; ibin++) - { - std::string projname = "_hproj_" + std::to_string(ibin); - if (!h_projy) - { - h_projy = h2->ProjectionY(projname.c_str(), ibin, ibin); - adcmean = h_projx->GetBinCenter(ibin); - adcnum = 1.; - } - else - { - TH1 *hadd = h2->ProjectionY(projname.c_str(), ibin, ibin); - h_projy->Add(hadd); - delete hadd; - adcmean += h_projx->GetBinCenter(ibin); - adcnum += 1.; - } - - if (h_projy->Integral() > 2000. || ibin == nbinsx) - { - adcmean /= adcnum; - - int maxbin = h_projy->GetMaximumBin(); - double xmax_g = h_projy->GetBinCenter(maxbin); - double ymax_g = h_projy->GetBinContent(maxbin); - gaussian.SetParameter(0, ymax_g); - gaussian.SetParameter(1, xmax_g); - gaussian.SetRange(xmax_g - 0.6, xmax_g + 0.6); - h_projy->Fit(&gaussian, "RWWQ"); - - double mean = gaussian.GetParameter(1); - double meanerr = gaussian.GetParError(1); - double rms = gaussian.GetParameter(2); - double rmserr = gaussian.GetParError(2); - - if (meanerr < 1.0) - { - int n = gridge->GetN(); - gridge->SetPoint(n, adcmean, mean); - gridge->SetPointError(n, 0., meanerr); - } - if (rmserr < 0.01) - { - int n = grms->GetN(); - grms->SetPoint(n, adcmean, rms); - grms->SetPointError(n, 0., rmserr); - } - - delete h_projy; - h_projy = nullptr; - adcmean = 0.; - adcnum = 0.; - } - } - - gridge->SetBit(TGraph::kIsSortedX); - grms->SetBit(TGraph::kIsSortedX); - delete h_projx; -} - -// --------------------------------------------------------------------------- -// FitAndWriteSlew — build slew-correction LUT from h2_slew ridge -// --------------------------------------------------------------------------- -void MbdCalibReco::FitAndWriteSlew() -{ - const int NPOINTS = 16000; - const int MINADC = 0; - const int MAXADC = 15999; - - std::string scorr_fname = _rundir + "/mbd_slewcorr.calib"; - std::ofstream scorr_file(scorr_fname); - if (!scorr_file.is_open()) - { - std::cout << Name() << "::FitAndWriteSlew ERROR cannot open " << scorr_fname << std::endl; - return; - } - - std::string trms_fname = _rundir + "/mbd_timerms.calib"; - std::ofstream trms_file(trms_fname); - - // Arrays of slew/trms graph pointers, indexed by feech (only T-channels used) - std::array g_slew{}; - std::array g_trms{}; - g_slew.fill(nullptr); - g_trms.fill(nullptr); - - for (int ipmt = 0; ipmt < MbdDefs::MBD_N_PMT; ipmt++) - { - if (!h2_slew[ipmt]) - { - continue; - } - - int feech_t = (ipmt / 8) * 16 + ipmt % 8; - - TGraphErrors *gr = nullptr; - TGraphErrors *grms_tmp = nullptr; - FindTH2Ridge(h2_slew[ipmt], gr, grms_tmp); - - g_slew[feech_t] = gr; - g_trms[feech_t] = grms_tmp; - - if (gr) - { - gr->SetName(("g_slew" + std::to_string(ipmt)).c_str()); - gr->SetMarkerStyle(20); - gr->SetMarkerSize(0.25); - } - if (grms_tmp) - { - grms_tmp->SetName(("g_trms" + std::to_string(ipmt)).c_str()); - grms_tmp->SetMarkerStyle(20); - grms_tmp->SetMarkerSize(0.25); - } - } - - // Write slew correction LUT (one T-channel feech at a time) - for (int ifeech = 0; ifeech < MbdDefs::MBD_N_FEECH; ifeech++) - { - // Only T-channels: type = (feech/8) % 2 == 0 - if ((ifeech / 8) % 2 == 1) - { - continue; - } - - if (!g_slew[ifeech]) - { - continue; - } - - scorr_file << ifeech << "\t" << NPOINTS << "\t" << MINADC << "\t" << MAXADC << "\n"; - int step = (MAXADC - MINADC) / (NPOINTS - 1); - for (int iadc = MINADC; iadc <= MAXADC; iadc += step) - { - scorr_file << g_slew[ifeech]->Eval(iadc) << " "; - if (iadc % 10 == 9) - { - scorr_file << "\n"; - } - } - } - scorr_file.close(); - - // Write time-RMS LUT - if (trms_file.is_open()) - { - for (int ifeech = 0; ifeech < MbdDefs::MBD_N_FEECH; ifeech++) - { - if ((ifeech / 8) % 2 == 1) - { - continue; - } - if (!g_trms[ifeech]) - { - continue; - } - - trms_file << ifeech << "\t" << NPOINTS << "\t" << MINADC << "\t" << MAXADC << "\n"; - int step = (MAXADC - MINADC) / (NPOINTS - 1); - for (int iadc = MINADC; iadc <= MAXADC; iadc += step) - { - trms_file << g_trms[ifeech]->Eval(iadc) << " "; - if (iadc % 10 == 9) - { - trms_file << "\n"; - } - } - } - trms_file.close(); - } - - // Write graphs to ROOT file and create CDB ROOT file - _outfile->cd(); - for (auto *g : g_slew) - { - if (g) - { - g->Write(); - } - } - for (auto *g : g_trms) - { - if (g) - { - g->Write(); - } - } - - { - MbdCalib tmpcal; - tmpcal.Download_SlewCorr(scorr_fname); - std::string cdb_fname = scorr_fname; - cdb_fname.replace(cdb_fname.rfind(".calib"), 6, ".root"); - tmpcal.Write_CDB_SlewCorr(cdb_fname); - } - - // Clean up - for (auto *g : g_slew) - { - delete g; - } - for (auto *g : g_trms) - { - delete g; - } - - std::cout << Name() << ": wrote " << scorr_fname << std::endl; -} - diff --git a/offline/packages/mbd/MbdCalibReco.h b/offline/packages/mbd/MbdCalibReco.h index f6d6448932..6a2c49c027 100644 --- a/offline/packages/mbd/MbdCalibReco.h +++ b/offline/packages/mbd/MbdCalibReco.h @@ -9,6 +9,7 @@ #include #include #include +#include class PHCompositeNode; class MbdCalib; @@ -32,7 +33,7 @@ class MbdCalibReco : public SubsysReco int Init(PHCompositeNode* topNode) override; int InitRun(PHCompositeNode* topNode) override; int process_event(PHCompositeNode* topNode) override; - int End(PHCompositeNode* topNode) override; + int EndRun(const int runnumber) override; void SetSubPass(const int s) { _subpass = s; } void SetCalDir(const std::string& d) { _caldir = d; } @@ -41,12 +42,6 @@ class MbdCalibReco : public SubsysReco private: int getNodes(PHCompositeNode* topNode); void InitHistos(); - int getRunType() const; - - void FitAndWriteT0(); - void FitAndWriteSlew(); - - void FindTH2Ridge(const TH2* h2, TGraphErrors*& gridge, TGraphErrors*& grms) const; uint64_t _mbias_trigger_mask{0}; @@ -56,7 +51,7 @@ class MbdCalibReco : public SubsysReco std::string _rundir; // _caldir// std::string _cdbtag{}; // non-empty → download from CDB instead of local files - MbdCalib* _mbdcal{nullptr}; + std::unique_ptr _mbdcal; MbdPmtContainer* _mbdpmts{nullptr}; MbdOut* _mbdout{nullptr}; MbdGeom* _mbdgeom{nullptr}; @@ -67,7 +62,7 @@ class MbdCalibReco : public SubsysReco std::array h_tt{}; std::array h_tq{}; std::array h_qp{}; - std::array h2_slew{}; + std::array h2_slew{}; TH2* h2_tt{nullptr}; TH2* h2_tq{nullptr}; From 94ccd3502c6acd8184d8bdd5b7f095becd5e5867 Mon Sep 17 00:00:00 2001 From: mchiu-bnl Date: Thu, 7 May 2026 16:20:37 -0400 Subject: [PATCH 3/6] rabbit fixes 2 --- offline/packages/mbd/MbdCalibReco.cc | 29 +++++++++++- offline/packages/mbd/MbdCalibReco.h | 3 +- offline/packages/mbd/MbdEvent.cc | 67 +++++++++++++--------------- 3 files changed, 61 insertions(+), 38 deletions(-) diff --git a/offline/packages/mbd/MbdCalibReco.cc b/offline/packages/mbd/MbdCalibReco.cc index 10ed14b559..e102649332 100644 --- a/offline/packages/mbd/MbdCalibReco.cc +++ b/offline/packages/mbd/MbdCalibReco.cc @@ -168,7 +168,7 @@ int MbdCalibReco::InitRun(PHCompositeNode *topNode) } std::cout << Name() << ": output file " << outfname << std::endl; - InitHistos(); + BookHistograms(); origdir->cd(); @@ -232,8 +232,15 @@ int MbdCalibReco::getNodes(PHCompositeNode *topNode) return Fun4AllReturnCodes::EVENT_OK; } -void MbdCalibReco::InitHistos() +void MbdCalibReco::BookHistograms() { + // Delete histograms if they have already have been booked. + if ( h2_tt ) + { + DeleteHistograms(); + return; + } + for (int ipmt = 0; ipmt < MbdDefs::MBD_N_PMT; ipmt++) { std::string sn = std::to_string(ipmt); @@ -271,6 +278,24 @@ void MbdCalibReco::InitHistos() h2_tq->SetYTitle("pmt ch"); } +void MbdCalibReco::DeleteHistograms() +{ + for (int ipmt = 0; ipmt < MbdDefs::MBD_N_PMT; ipmt++) + { + if ( h_tt[ipmt] ) delete h_tt[ipmt]; + if ( h_tq[ipmt] ) delete h_tq[ipmt]; + if ( h_qp[ipmt] ) delete h_qp[ipmt]; + + if ( h2_slew[ipmt] ) + { + delete h2_slew[ipmt]; + } + } + + if ( h2_tt ) delete h2_tt; + if ( h2_tq ) delete h2_tq; +} + int MbdCalibReco::process_event(PHCompositeNode * /*topNode*/) { // Require a scaled "MBD N&S" trigger diff --git a/offline/packages/mbd/MbdCalibReco.h b/offline/packages/mbd/MbdCalibReco.h index 6a2c49c027..eed96688cd 100644 --- a/offline/packages/mbd/MbdCalibReco.h +++ b/offline/packages/mbd/MbdCalibReco.h @@ -41,7 +41,8 @@ class MbdCalibReco : public SubsysReco private: int getNodes(PHCompositeNode* topNode); - void InitHistos(); + void BookHistograms(); + void DeleteHistograms(); uint64_t _mbias_trigger_mask{0}; diff --git a/offline/packages/mbd/MbdEvent.cc b/offline/packages/mbd/MbdEvent.cc index 5c1e9ed087..f768dfa16d 100644 --- a/offline/packages/mbd/MbdEvent.cc +++ b/offline/packages/mbd/MbdEvent.cc @@ -300,45 +300,42 @@ int MbdEvent::InitRun() // zero out the tt_t0, tq_t0, and gains to produce uncalibrated time and charge // or load pass2 calibs from local file for calpass2+, if local files exist - if ( _calpass>1 ) + std::string calfname = "results/"; calfname += std::to_string(_runnum); calfname += "/mbd_tt_t0.calib"; + if ( std::filesystem::exists(calfname) ) { - std::string calfname = "results/"; calfname += std::to_string(_runnum); calfname += "/mbd_tt_t0.calib"; - if ( std::filesystem::exists(calfname) ) - { - std::cout << "Loading local mbd_tt_t0, " << calfname << std::endl; - _mbdcal->Download_TTT0( calfname ); - } - else - { - _mbdcal->Reset_TTT0(); - std::cout << PHWHERE << "local mbd_tt_t0 not found, reset to 0: " << calfname << std::endl; - } - - calfname = "results/"; calfname += std::to_string(_runnum); calfname += "/mbd_tq_t0.calib"; - if ( std::filesystem::exists(calfname) ) - { - std::cout << "Loading local mbd_tq_t0, " << calfname << std::endl; - _mbdcal->Download_TQT0( calfname ); - } - else - { - _mbdcal->Reset_TQT0(); - std::cout << PHWHERE << "local mbd_tq_t0 not found, reset to 0: " << calfname << std::endl; - } + std::cout << "Loading local mbd_tt_t0, " << calfname << std::endl; + _mbdcal->Download_TTT0( calfname ); + } + else + { + _mbdcal->Reset_TTT0(); + std::cout << PHWHERE << "local mbd_tt_t0 not found, reset to 0: " << calfname << std::endl; + } - calfname = "results/"; calfname += std::to_string(_runnum); calfname += "/mbd_qfit.calib"; - if ( std::filesystem::exists(calfname) ) - { - std::cout << "Loading local mbd_qfit, " << calfname << std::endl; - _mbdcal->Download_Gains( calfname ); - } - else - { - _mbdcal->Reset_Gains(); - std::cout << PHWHERE << "local mbd_gains not found, reset to 1: " << calfname << std::endl; - } + calfname = "results/"; calfname += std::to_string(_runnum); calfname += "/mbd_tq_t0.calib"; + if ( std::filesystem::exists(calfname) ) + { + std::cout << "Loading local mbd_tq_t0, " << calfname << std::endl; + _mbdcal->Download_TQT0( calfname ); + } + else + { + _mbdcal->Reset_TQT0(); + std::cout << PHWHERE << "local mbd_tq_t0 not found, reset to 0: " << calfname << std::endl; + } + calfname = "results/"; calfname += std::to_string(_runnum); calfname += "/mbd_qfit.calib"; + if ( std::filesystem::exists(calfname) ) + { + std::cout << "Loading local mbd_qfit, " << calfname << std::endl; + _mbdcal->Download_Gains( calfname ); } + else + { + _mbdcal->Reset_Gains(); + std::cout << PHWHERE << "local mbd_gains not found, reset to 1: " << calfname << std::endl; + } + TDirectory *orig_dir = gDirectory; From 5bbb0c2aeab8a1aeffc98ef9699109adcbccc1bd Mon Sep 17 00:00:00 2001 From: mchiu-bnl Date: Thu, 7 May 2026 16:47:50 -0400 Subject: [PATCH 4/6] clang-tidy fixes --- offline/packages/mbd/MbdCalibReco.cc | 20 ++++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) diff --git a/offline/packages/mbd/MbdCalibReco.cc b/offline/packages/mbd/MbdCalibReco.cc index e102649332..d2acec905d 100644 --- a/offline/packages/mbd/MbdCalibReco.cc +++ b/offline/packages/mbd/MbdCalibReco.cc @@ -282,18 +282,26 @@ void MbdCalibReco::DeleteHistograms() { for (int ipmt = 0; ipmt < MbdDefs::MBD_N_PMT; ipmt++) { - if ( h_tt[ipmt] ) delete h_tt[ipmt]; - if ( h_tq[ipmt] ) delete h_tq[ipmt]; - if ( h_qp[ipmt] ) delete h_qp[ipmt]; - + if ( h_tt[ipmt] ) + { + delete h_tt[ipmt]; + } + if ( h_tq[ipmt] ) + { + delete h_tq[ipmt]; + } + if ( h_qp[ipmt] ) + { + delete h_qp[ipmt]; + } if ( h2_slew[ipmt] ) { delete h2_slew[ipmt]; } } - if ( h2_tt ) delete h2_tt; - if ( h2_tq ) delete h2_tq; + delete h2_tt; + delete h2_tq; } int MbdCalibReco::process_event(PHCompositeNode * /*topNode*/) From 4134ac0bf46aaf1f7c7ce15fb4a455c3e865d121 Mon Sep 17 00:00:00 2001 From: mchiu-bnl Date: Thu, 7 May 2026 17:06:28 -0400 Subject: [PATCH 5/6] code rabbit fix 3 --- offline/packages/mbd/MbdCalibReco.cc | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/offline/packages/mbd/MbdCalibReco.cc b/offline/packages/mbd/MbdCalibReco.cc index d2acec905d..f131a0d6c2 100644 --- a/offline/packages/mbd/MbdCalibReco.cc +++ b/offline/packages/mbd/MbdCalibReco.cc @@ -346,6 +346,17 @@ int MbdCalibReco::process_event(PHCompositeNode * /*topNode*/) } Short_t pmtno = pmt->get_pmt(); + if ( pmtno<0 || pmtno>128 ) + { + static int counter = 0; + if ( counter<10 ) + { + std::cerr << PHWHERE << " invalide pmt no " << pmtno << std::endl; + counter++; + } + continue; + } + Float_t q = pmt->get_q(); Float_t tt = pmt->get_tt(); Float_t tq = pmt->get_tq(); From c7e5387d425f592730b11c385a61624a9ad46d2d Mon Sep 17 00:00:00 2001 From: mchiu-bnl Date: Thu, 7 May 2026 17:09:08 -0400 Subject: [PATCH 6/6] code rabbit fix 3.1 --- offline/packages/mbd/MbdCalibReco.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/offline/packages/mbd/MbdCalibReco.cc b/offline/packages/mbd/MbdCalibReco.cc index f131a0d6c2..3b64e0f63e 100644 --- a/offline/packages/mbd/MbdCalibReco.cc +++ b/offline/packages/mbd/MbdCalibReco.cc @@ -346,7 +346,7 @@ int MbdCalibReco::process_event(PHCompositeNode * /*topNode*/) } Short_t pmtno = pmt->get_pmt(); - if ( pmtno<0 || pmtno>128 ) + if ( pmtno<0 || pmtno>=MbdDefs::MBD_N_PMT ) { static int counter = 0; if ( counter<10 )