packages/mbd/MbdCalibReco.cc

0001 #include "MbdCalibReco.h"
0002 #include "MbdCalib.h"
0003 #include "MbdDefs.h"
0004 #include "MbdPmtContainer.h"
0005 #include "MbdPmtHit.h"
0006 #include "MbdOut.h"
0007
0008 #include <ffamodules/CDBInterface.h>
0009 #include <ffaobjects/EventHeader.h>
0010 #include <ffaobjects/RunHeader.h>
0011 #include <ffarawobjects/Gl1Packet.h>
0012 #include <fun4all/Fun4AllReturnCodes.h>
0013 #include <phool/PHCompositeNode.h>
0014 #include <phool/getClass.h>
0015 #include <phool/phool.h>
0016 #include <phool/recoConsts.h>
0017
0018 #include <TF1.h>
0019 #include <TFile.h>
0020 #include <TGraphErrors.h>
0021 #include <TH1.h>
0022 #include <TH2.h>
0023 #include <TROOT.h>
0024 #include <TSystem.h>
0025 #include <TDirectory.h>
0026
0027 #include <cmath>
0028 #include <filesystem>
0029 #include <fstream>
0030 #include <iostream>
0031 #include <sstream>
0032 #include <string>
0033
0034 MbdCalibReco::MbdCalibReco(const std::string &name)
0035   : SubsysReco(name)
0036 {
0037 }
0038
0039 int MbdCalibReco::Init(PHCompositeNode * /*topNode*/)
0040 {
0041   _mbdcal = std::make_unique<MbdCalib>();
0042   _mbdcal->Verbosity( Verbosity() );
0043   return Fun4AllReturnCodes::EVENT_OK;
0044 }
0045
0046 int MbdCalibReco::InitRun(PHCompositeNode *topNode)
0047 {
0048   _runheader = findNode::getClass<RunHeader>(topNode, "RunHeader");
0049   if (!_runheader)
0050   {
0051     std::cout << PHWHERE << " RunHeader node not found, will use run number 0" << std::endl;
0052   }
0053
0054   _runnumber = _runheader ? _runheader->get_RunNumber() : 0;
0055
0056   getNodes(topNode);
0057
0058   // Build run directory path and create it
0059   std::ostringstream oss;
0060   oss << _caldir << "/" << _runnumber;
0061   _rundir = oss.str();
0062   gSystem->Exec(("mkdir -p " + _rundir).c_str());
0063
0064   if (!_cdbtag.empty())
0065   {
0066     // Download baseline calibrations from CDB
0067     recoConsts::instance()->set_StringFlag("CDB_GLOBALTAG", _cdbtag);
0068     CDBInterface* cdb = CDBInterface::instance();
0069     std::string url;
0070
0071     url = cdb->getUrl("MBD_SAMPMAX");
0072     if (!url.empty()) { _mbdcal->Download_SampMax(url); }
0073
0074     url = cdb->getUrl("MBD_PED");
0075     if (!url.empty()) { _mbdcal->Download_Ped(url); }
0076
0077     url = cdb->getUrl("MBD_TIMECORR");
0078     if (!url.empty()) { _mbdcal->Download_TimeCorr(url); }
0079
0080     url = cdb->getUrl("MBD_SLEWCORR");
0081     if (!url.empty()) { _mbdcal->Download_SlewCorr(url); }
0082
0083     std::cout << Name() << ": loaded calibrations from CDB tag " << _cdbtag << std::endl;
0084   }
0085   else
0086   {
0087     // Load baseline calibrations from local files if they exist
0088     std::string calfile = _rundir + "/mbd_sampmax.calib";
0089     if (gSystem->AccessPathName(calfile.c_str()) == 0)
0090     {
0091       _mbdcal->Download_SampMax(calfile);
0092     }
0093     calfile = _rundir + "/mbd_ped.calib";
0094     if (gSystem->AccessPathName(calfile.c_str()) == 0)
0095     {
0096       _mbdcal->Download_Ped(calfile);
0097     }
0098
0099     // Load slew correction for subpass >= 2
0100     if (_subpass >= 2)
0101     {
0102       calfile = _rundir + "/mbd_slewcorr.calib";
0103       if (gSystem->AccessPathName(calfile.c_str()) == 0)
0104       {
0105         _mbdcal->Download_SlewCorr(calfile);
0106         std::cout << Name() << ": loaded " << calfile << std::endl;
0107       }
0108       else
0109       {
0110         std::cout << Name() << ": WARNING: " << calfile << " not found" << std::endl;
0111       }
0112     }
0113   }
0114
0115   // Load t0 offsets for subpass >= 1 (always from local files — outputs of previous subpass)
0116   if (_subpass >= 1)
0117   {
0118     std::string prevpass = "pass" + std::to_string(_subpass - 1) + "_";
0119
0120     std::string calfile = _rundir + "/" + prevpass + "mbd_tq_t0.calib";
0121     if (gSystem->AccessPathName(calfile.c_str()) == 0)
0122     {
0123       _mbdcal->Download_TQT0(calfile);
0124       std::cout << Name() << ": loaded " << calfile << std::endl;
0125     }
0126     else
0127     {
0128       std::cout << Name() << ": WARNING: " << calfile << " not found" << std::endl;
0129     }
0130
0131     calfile = _rundir + "/" + prevpass + "mbd_tt_t0.calib";
0132     if (gSystem->AccessPathName(calfile.c_str()) == 0)
0133     {
0134       _mbdcal->Download_TTT0(calfile);
0135       std::cout << Name() << ": loaded " << calfile << std::endl;
0136     }
0137     else
0138     {
0139       std::cout << Name() << ": WARNING: " << calfile << " not found" << std::endl;
0140     }
0141   }
0142
0143   // Build bitmask of scaled triggers whose names begin with "MBD N&S"
0144   _mbias_trigger_mask = 0xfc00;
0145
0146   // Open output ROOT file
0147   TDirectory *origdir = gDirectory;
0148
0149   std::string outfname = _rundir + "/calmbdpass2." + std::to_string(_subpass);
0150   if (_subpass == 0)
0151   {
0152     outfname += "_time-" + std::to_string(_runnumber) + ".root";
0153   }
0154   else if (_subpass == 1 || _subpass == 2)
0155   {
0156     outfname += "_slew-" + std::to_string(_runnumber) + ".root";
0157   }
0158   else
0159   {
0160     outfname += "_q-" + std::to_string(_runnumber) + ".root";
0161   }
0162   _outfile = std::make_unique<TFile>(outfname.c_str(), "RECREATE");
0163   if (!_outfile || _outfile->IsZombie())
0164   {
0165     std::cerr << PHWHERE << " ERROR: cannot open output file " << outfname << std::endl;
0166     _outfile.reset();
0167     return Fun4AllReturnCodes::ABORTRUN;
0168   }
0169   std::cout << Name() << ": output file " << outfname << std::endl;
0170
0171   BookHistograms();
0172
0173   origdir->cd();
0174
0175   return Fun4AllReturnCodes::EVENT_OK;
0176 }
0177
0178 int MbdCalibReco::getNodes(PHCompositeNode *topNode)
0179 {
0180   _evtheader = findNode::getClass<EventHeader>(topNode, "EventHeader");
0181   if (!_evtheader)
0182   {
0183     std::cout << PHWHERE << " EvtHeader not found, will use run number 0" << std::endl;
0184   }
0185
0186   _gl1packet = findNode::getClass<Gl1Packet>(topNode,14001);
0187   if (!_gl1packet)
0188   {
0189     _gl1packet = findNode::getClass<Gl1Packet>(topNode, "GL1Packet");
0190     static int counter = 0;
0191     if ( counter<4 )
0192     {
0193       std::cout << PHWHERE << " GL1Packet not found" << std::endl;
0194     }
0195   }
0196
0197   _mbdpmts = findNode::getClass<MbdPmtContainer>(topNode, "MbdPmtContainer");
0198   if (!_mbdpmts)
0199   {
0200     static int counter = 0;
0201     if ( counter<4 )
0202     {
0203       std::cout << PHWHERE << " MbdPmtContainer not found" << std::endl;
0204     }
0205   }
0206
0207   _mbdout = findNode::getClass<MbdOut>(topNode, "MbdOut");
0208   if (!_mbdout)
0209   {
0210     static int counter = 0;
0211     if ( counter<4 )
0212     {
0213       std::cout << PHWHERE << " MbdOut not found" << std::endl;
0214     }
0215   }
0216
0217   _mbdgeom = findNode::getClass<MbdGeom>(topNode, "MbdGeom");
0218   if (!_mbdgeom)
0219   {
0220     static int counter = 0;
0221     if ( counter<4 )
0222     {
0223       std::cout << PHWHERE << " MbdGeom not found" << std::endl;
0224     }
0225   }
0226
0227   if ( !_mbdgeom || !_mbdout || !_mbdpmts || !_gl1packet || !_evtheader )
0228   {
0229     return Fun4AllReturnCodes::ABORTRUN;
0230   }
0231
0232   return Fun4AllReturnCodes::EVENT_OK;
0233 }
0234
0235 void MbdCalibReco::BookHistograms()
0236 {
0237   // Delete histograms if they have already have been booked.
0238   if ( h2_tt )
0239   {
0240     DeleteHistograms();
0241     return;
0242   }
0243
0244   for (int ipmt = 0; ipmt < MbdDefs::MBD_N_PMT; ipmt++)
0245   {
0246     std::string sn = std::to_string(ipmt);
0247
0248     h_tt[ipmt] = new TH1F(("h_tt" + sn).c_str(), ("tt" + sn).c_str(), 7000, -30., 30.);
0249     h_tt[ipmt]->SetXTitle("ns");
0250
0251     h_tq[ipmt] = new TH1F(("h_tq" + sn).c_str(), ("tq" + sn).c_str(), 7000, -150., 31. * 17.7623);
0252     h_tq[ipmt]->SetXTitle("ns");
0253
0254     h_qp[ipmt] = new TH1F(("h_q" + sn).c_str(), ("q" + sn).c_str(), 3000, -100., 14900.);
0255     h_qp[ipmt]->SetXTitle("ADC");
0256
0257     if (_subpass >= 1)
0258     {
0259       const int    nbins[2] = {4000, 1100};
0260       const double xmin[2]  = {-0.5, -5.};
0261       const double xmax[2]  = {16000. - 0.5, 6.};
0262       h2_slew[ipmt] = new THnSparseF(("h2_slew" + sn).c_str(), ("slew curve, ch " + sn).c_str(), 2, nbins, xmin, xmax);
0263       h2_slew[ipmt]->GetAxis(0)->SetTitle("ADC");
0264       h2_slew[ipmt]->GetAxis(1)->SetTitle("#Delta T (ns)");
0265     }
0266     else
0267     {
0268       h2_slew[ipmt] = nullptr;
0269     }
0270   }
0271
0272   h2_tt = new TH2F("h2_tt", "ch vs tt", 900, -150., 150., MbdDefs::MBD_N_PMT, -0.5, MbdDefs::MBD_N_PMT - 0.5);
0273   h2_tt->SetXTitle("tt [ns]");
0274   h2_tt->SetYTitle("pmt ch");
0275
0276   h2_tq = new TH2F("h2_tq", "ch vs tq", 900, -150., 150., MbdDefs::MBD_N_PMT, -0.5, MbdDefs::MBD_N_PMT - 0.5);
0277   h2_tq->SetXTitle("tq [ns]");
0278   h2_tq->SetYTitle("pmt ch");
0279 }
0280
0281 void MbdCalibReco::DeleteHistograms()
0282 {
0283   for (int ipmt = 0; ipmt < MbdDefs::MBD_N_PMT; ipmt++)
0284   {
0285     if ( h_tt[ipmt] )
0286     {
0287       delete h_tt[ipmt];
0288     }
0289     if ( h_tq[ipmt] )
0290     {
0291       delete h_tq[ipmt];
0292     }
0293     if ( h_qp[ipmt] )
0294     {
0295       delete h_qp[ipmt];
0296     }
0297     if ( h2_slew[ipmt] )
0298     {
0299       delete h2_slew[ipmt];
0300     }
0301   }
0302
0303   delete h2_tt;
0304   delete h2_tq;
0305 }
0306
0307 int MbdCalibReco::process_event(PHCompositeNode * /*topNode*/)
0308 {
0309   // Require a scaled "MBD N&S" trigger
0310   if (_mbias_trigger_mask != 0)
0311   {
0312     uint64_t strig = _gl1packet->getScaledVector();  // scaled trigger only
0313     if ( (strig&_mbias_trigger_mask)==0 )
0314     {
0315       return Fun4AllReturnCodes::ABORTEVENT;
0316     }
0317   }
0318
0319   std::array<Float_t, MbdDefs::MBD_N_ARMS> armtime{};
0320   armtime.fill(0);
0321   std::array<Float_t, MbdDefs::MBD_N_ARMS> nhit{};
0322   nhit.fill(0);
0323
0324   Float_t zvtx = _mbdout->get_zvtx();
0325   // Vertex cut for subpass >= 1
0326   if ( _subpass >= 1 )
0327   {
0328     if (std::abs(zvtx) > 60.)
0329     {
0330       return Fun4AllReturnCodes::EVENT_OK;
0331     }
0332   }
0333
0334   for (int iarm=0; iarm<2; iarm++)
0335   {
0336     armtime[iarm] = _mbdout->get_time(iarm);
0337     nhit[iarm] = _mbdout->get_npmt(iarm);
0338   }
0339
0340   for (int ipmt=0; ipmt < _mbdpmts->get_npmt(); ipmt++)
0341   {
0342     MbdPmtHit *pmt = _mbdpmts->get_pmt(ipmt);
0343     if ( !pmt )
0344     {
0345       continue;
0346     }
0347
0348     Short_t pmtno = pmt->get_pmt();
0349     if ( pmtno<0 || pmtno>=MbdDefs::MBD_N_PMT )
0350     {
0351       static int counter = 0;
0352       if ( counter<10 )
0353       {
0354         std::cerr << PHWHERE << " invalide pmt no " << pmtno << std::endl;
0355         counter++;
0356       }
0357       continue;
0358     }
0359
0360     Float_t q = pmt->get_q();
0361     Float_t tt  = pmt->get_tt();
0362     Float_t tq  = pmt->get_tq();
0363
0364     h_tt[pmtno]->Fill( tt );
0365     h2_tt->Fill( tt, pmtno );
0366     h_tq[pmtno]->Fill( tq );
0367     h2_tq->Fill( tq, pmtno );
0368
0369     // Fill charge histogram for in-time hits
0370     if ( std::abs(tt)<26.0 && q > 0.)
0371     {
0372       h_qp[pmtno]->Fill( q );
0373     }
0374
0375     int arm = _mbdgeom->get_arm( pmtno );
0376
0377     // Fill slew histogram for subpass >= 1
0378     if (_subpass >= 1 && h2_slew[pmtno])
0379     {
0380       if (nhit[arm] >= 2. && q > 0.)
0381       {
0382         float dt = tt - armtime[arm];
0383         const double coords[2] = {q, dt};
0384         h2_slew[pmtno]->Fill(coords);
0385       }
0386     }
0387   }
0388
0389   return Fun4AllReturnCodes::EVENT_OK;
0390 }
0391
0392 int MbdCalibReco::EndRun(const int /*runnumber*/)
0393 {
0394   if (!_outfile)
0395   {
0396     return Fun4AllReturnCodes::EVENT_OK;
0397   }
0398
0399   // Write histograms to output file
0400   _outfile->cd();
0401   h2_tt->Write();
0402   h2_tq->Write();
0403   for (int ipmt = 0; ipmt < MbdDefs::MBD_N_PMT; ipmt++)
0404   {
0405     h_tt[ipmt]->Write();
0406     h_tq[ipmt]->Write();
0407     h_qp[ipmt]->Write();
0408     if (h2_slew[ipmt])
0409     {
0410       h2_slew[ipmt]->Write();
0411     }
0412   }
0413
0414   _outfile->Close();
0415
0416   return Fun4AllReturnCodes::EVENT_OK;
0417 }
0418