Skip to content

Add time calibration to reco digis #2030

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
michaelmackenzie merged 9 commits into from
Oct 13, 2026
Merged

Add time calibration to reco digis #2030

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
bc5f1cc
Add time calibration to reco digis
giro94 Oct 8, 2026
bec92d2
Merge branch 'Mu2e:main' into dev
giro94 Oct 8, 2026
0728f73
HitsFromDTCEvents now uses prodition, and applies both energy and tim…
giro94 Oct 13, 2026
0f1d931
CaloHitMakerFast now uses prodition, and applies both energy and time…
giro94 Oct 13, 2026
4a2217c
remove obsolete variables
giro94 Oct 13, 2026
fbf726a
bugfix
giro94 Oct 13, 2026
e4d763e
Forgot to forward calCalib inside function
giro94 Oct 13, 2026
5c0bc8e
removed unused beginRun
giro94 Oct 13, 2026
80044dc
Removed adc2mev fcl parameter
giro94 Oct 13, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 0 additions & 1 deletion CaloReco/fcl/prolog.fcl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -67,7 +67,6 @@ CaloReco : {
nPEperMeV : 30
nSigmaNoise : 4
noiseLevelMeV : 5.5e-1
ADCToMeV : @local::ADCToMeV
diagLevel : 0
}
}
Expand Down
24 changes: 15 additions & 9 deletions CaloReco/src/CaloHitMakerFast_module.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,8 +4,10 @@
#include "art/Framework/Principal/Handle.h"
#include "fhiclcpp/types/Sequence.h"

#include "Offline/CaloConditions/inc/CalCalib.hh"
#include "Offline/DataProducts/inc/CaloConst.hh"
#include "Offline/DataProducts/inc/CaloSiPMId.hh"
#include "Offline/ProditionsService/inc/ProditionsHandle.hh"
#include "Offline/RecoDataProducts/inc/CaloDigi.hh"
#include "Offline/RecoDataProducts/inc/CaloHit.hh"
#include "Offline/RecoDataProducts/inc/ProtonBunchTime.hh"
Expand Down Expand Up @@ -48,7 +50,6 @@ namespace mu2e {
fhicl::Atom<float> nSigmaNoise { Name("nSigmaNoise"), Comment("Maxnumber of sigma Noise to combine digi") };
fhicl::Atom<float> caphriEDepMax { Name("caphriEDepMax"), Comment("Maximum CAPHRI hit energy in MeV")};
fhicl::Atom<float> caphriEDepMin { Name("caphriEDepMin"), Comment("Minimum CAPHRI hit energy in MeV")};
fhicl::Atom<double> ADCToMeV { Name("ADCToMeV"), Comment("ADC to MeV conversion factor") };
fhicl::Atom<int> diagLevel { Name("diagLevel"), Comment("Diag Level"),0 };
};

Expand All @@ -62,7 +63,6 @@ namespace mu2e {
nSigmaNoise_ (config().nSigmaNoise()),
caphriEDepMax_ (config().caphriEDepMax()),
caphriEDepMin_ (config().caphriEDepMin()),
ADCToMeV_ (config().ADCToMeV()),
diagLevel_ (config().diagLevel())
{
produces<CaloHitCollection>("calo");
Expand All @@ -74,9 +74,11 @@ namespace mu2e {


private:
ProditionsHandle<mu2e::CalCalib> calCalibHandle_;

using pulseMapType = std::map<unsigned, std::vector<HitInfo>>;

void extractHits(const CaloDigiCollection& caloDigis, CaloHitCollection& caloHitsColl, CaloHitCollection& caphriHitsColl, IntensityInfoCalo& intCalo, double pbtOffset);
void extractHits(const CaloDigiCollection& caloDigis, CaloHitCollection& caloHitsColl, CaloHitCollection& caphriHitsColl, IntensityInfoCalo& intCalo, double pbtOffset, CalCalib const& calCalib);
void addPulse(pulseMapType& pulseMap, unsigned crystalID, float time, float eDep);

art::ProductToken<CaloDigiCollection> caloDigisToken_;
Expand All @@ -87,7 +89,6 @@ namespace mu2e {
float nSigmaNoise_;
float caphriEDepMax_;
float caphriEDepMin_;
float ADCToMeV_;
int diagLevel_;
};

Expand All @@ -101,14 +102,15 @@ namespace mu2e {
auto pbtH = event.getValidHandle(pbttoken_);
const ProtonBunchTime& pbt(*pbtH);
double pbtOffset = pbt.pbtime_;
CalCalib const& calCalib = calCalibHandle_.get(event.id());

@bechenard bechenard Oct 13, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

AArgh, the better style is obviously const CalCalib&....

@giro94 giro94 Oct 13, 2026

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Totally agree, I copied the style somewhere...
Will fix once for all

@giro94 giro94 Oct 13, 2026

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more. 

grep "const&" */*/*.cc | wc -l
2825

No I will not

const auto& caloDigisHandle = event.getValidHandle(caloDigisToken_);
const auto& caloDigis(*caloDigisHandle);

auto caloHitsColl = std::make_unique<CaloHitCollection>();
auto caphriHitsColl = std::make_unique<CaloHitCollection>();
auto intInfo = std::make_unique<IntensityInfoCalo>();
extractHits(caloDigis,*caloHitsColl,*caphriHitsColl,*intInfo,pbtOffset);
extractHits(caloDigis,*caloHitsColl,*caphriHitsColl,*intInfo,pbtOffset,calCalib);

event.put(std::move(caloHitsColl), "calo");
event.put(std::move(caphriHitsColl),"caphri");
Expand All @@ -120,7 +122,7 @@ namespace mu2e {


//--------------------------------------------------------------------------------------------------------------
void CaloHitMakerFast::extractHits(const CaloDigiCollection& caloDigis, CaloHitCollection& caloHitsColl, CaloHitCollection& caphriHitsColl, IntensityInfoCalo& intInfo, double pbtOffset)
void CaloHitMakerFast::extractHits(const CaloDigiCollection& caloDigis, CaloHitCollection& caloHitsColl, CaloHitCollection& caphriHitsColl, IntensityInfoCalo& intInfo, double pbtOffset, CalCalib const& calCalib)
{
std::vector<size_t> digis(caloDigis.size());
std::iota(digis.begin(),digis.end(),0);
Expand All @@ -143,9 +145,13 @@ namespace mu2e {
for (size_t i=0; i<nSamPed; ++i){ baseline += caloDigi.waveform().at(i);}
baseline /= nSamPed;

double eDep = (caloDigi.waveform().at(caloDigi.peakpos())-baseline)*ADCToMeV_;//FIXME! we should use the function ::Peak2MeV, I also think that we should: (i) discard the hit if eDep is <0 (noise/stange pulse), (ii) require a minimum pulse length. gianipez
double time = (caloDigi.t0() + caloDigi.peakpos())*CaloConst::_digitizationPeriod - pbtOffset; //Giani's definition; t0 in ticks
//double time = (caloDigi.t0() + caloDigi.peakpos()+0.5)*CaloConst::_digitizationPeriod - pbtOffset; //Bertrand's definition
double adc2MeV = calCalib.ADC2MeV(caloDigi.SiPMID());
double timeOffset = calCalib.timeoffset(caloDigi.SiPMID());
double eDep = (caloDigi.waveform().at(caloDigi.peakpos())-baseline)*adc2MeV;
double time = (caloDigi.t0() + caloDigi.peakpos())*CaloConst::_digitizationPeriod - pbtOffset - timeOffset; //Giani's definition; t0 in ticks
//double time = (caloDigi.t0() + caloDigi.peakpos()+0.5)*CaloConst::_digitizationPeriod - pbtOffset - timeOffset; //Bertrand's definition

//FIXME! I think that we should: (i) discard the hit if eDep is <0 (noise/stange pulse), (ii) require a minimum pulse length. gianipez

addPulse(pulseMap, crystalID, time, eDep);
if (diagLevel_ > 2) std::cout<<"[CaloHitMakerFast] extracted Digi with crystalID="<<crystalID<<" eDep="<<eDep<<"\t time=" <<time<<std::endl;
Expand Down
3 changes: 2 additions & 1 deletion CaloReco/src/CaloRecoDigiMaker_module.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -142,6 +142,7 @@ void CaloRecoDigiMaker::extractRecoDigi(const art::ValidHandle<CaloDigiCollectio
int SiPMID = caloDigi.SiPMID();
double t0 = caloDigi.t0() * CaloConst::_digitizationPeriod; // CaloDigi t0 is in digitizer clock ticks
double adc2MeV = calCalib.ADC2MeV(SiPMID);
double timeOffset = calCalib.timeoffset(SiPMID);

const std::vector<int>& waveform = caloDigi.waveform();

Expand All @@ -164,7 +165,7 @@ void CaloRecoDigiMaker::extractRecoDigi(const art::ValidHandle<CaloDigiCollectio
for (int i = 0; i < waveformProcessor_->nPeaks(); ++i) {
double eDep = waveformProcessor_->amplitude(i) * adc2MeV;
double eDepErr = waveformProcessor_->amplitudeErr(i) * adc2MeV;
double time = waveformProcessor_->time(i) - pbtOffset; // correct to time since protons
double time = waveformProcessor_->time(i) - pbtOffset - timeOffset;
double timeErr = waveformProcessor_->timeErr(i);
bool isPileUp = waveformProcessor_->isPileUp(i);
double chi2 = waveformProcessor_->chi2();
Expand Down
36 changes: 14 additions & 22 deletions DAQ/src/CaloHitsFromDTCEvents_module.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@
#include "artdaq-core-mu2e/Overlays/FragmentType.hh"
#include <artdaq-core/Data/Fragment.hh>

#include "Offline/CaloConditions/inc/CalCalib.hh"
#include "Offline/CaloConditions/inc/CalDAQMap.hh"
#include "Offline/DAQ/inc/CaloDAQUtilities.hh"
#include "Offline/DataProducts/inc/CaloConst.hh"
Expand Down Expand Up @@ -77,16 +78,16 @@ class art::CaloHitsFromDataDTCEvents : public EDProducer {
explicit CaloHitsFromDataDTCEvents(const art::EDProducer::Table<Config>& config);
virtual ~CaloHitsFromDataDTCEvents() {}

virtual void beginRun(art::Run&) override;

// --- Production:
virtual void produce(Event&);
virtual void endJob();

private:
mu2e::ProditionsHandle<mu2e::CalCalib> _caloCalib_h;
mu2e::ProditionsHandle<mu2e::CalDAQMap> _calodaqconds_h;

void analyze_calorimeter_(mu2e::CalDAQMap const& calodaqconds,
mu2e::CalCalib const& calCalib,
const mu2e::CalorimeterDataDecoder& cc,
std::unique_ptr<mu2e::CaloHitCollection> const& calo_hits,
std::unique_ptr<mu2e::CaloHitCollection> const& caphri_hits,
Expand All @@ -109,9 +110,6 @@ class art::CaloHitsFromDataDTCEvents : public EDProducer {
pulseMap_; // Temporary hack until the Calorimeter channel map is finalized
mu2e::CaloDAQUtilities caloDAQUtil_;

std::array<float, mu2e::CaloConst::_nCrystalChannel> peakADC2MeV_;
std::array<float, mu2e::CaloConst::_nCrystalChannel> timeCalib_;

long int total_events;
long int total_hits;
long int total_hits_good;
Expand All @@ -121,17 +119,6 @@ class art::CaloHitsFromDataDTCEvents : public EDProducer {

// ======================================================================

void art::CaloHitsFromDataDTCEvents::beginRun(art::Run& Run) {

// FIX ME!
// here we need to load the prodition-service with the calibrations
// for converting teh peakADC into MeV and sync times
for (size_t i = 0; i < peakADC2MeV_.size(); ++i) {
peakADC2MeV_[i] = 0.0461333;
timeCalib_[i] = 0.;
}
}

void art::CaloHitsFromDataDTCEvents::addPulse(
uint16_t& crystalID, float& time, float& eDep,
std::unique_ptr<mu2e::CaloHitCollection> const& hits_calo,
Expand Down Expand Up @@ -201,6 +188,7 @@ void art::CaloHitsFromDataDTCEvents::produce(Event& event) {
art::EventNumber_t eventNumber = event.event();

mu2e::CalDAQMap const& calodaqconds = _calodaqconds_h.get(event.id());
mu2e::CalCalib const& calCalib = _caloCalib_h.get(event.id());

// Collection of CaloHits for the event
std::unique_ptr<mu2e::CaloHitCollection> calo_hits(new mu2e::CaloHitCollection);
Expand All @@ -221,7 +209,7 @@ void art::CaloHitsFromDataDTCEvents::produce(Event& event) {
for (auto& subevent : caloSEvents) {
decoderColl->emplace_back(subevent);
auto& decoder = decoderColl->back();
analyze_calorimeter_(calodaqconds, decoder, calo_hits, caphri_hits, int_info);
analyze_calorimeter_(calodaqconds, calCalib, decoder, calo_hits, caphri_hits, int_info);
for (size_t i = 0; i < decoder.block_count(); ++i) {
totalSize += decoder.blockSizeBytes(i);
}
Expand Down Expand Up @@ -257,7 +245,7 @@ void art::CaloHitsFromDataDTCEvents::produce(Event& event) {
} // produce()

void art::CaloHitsFromDataDTCEvents::analyze_calorimeter_(
mu2e::CalDAQMap const& calodaqconds, const mu2e::CalorimeterDataDecoder& cc,
mu2e::CalDAQMap const& calodaqconds, mu2e::CalCalib const& calCalib, const mu2e::CalorimeterDataDecoder& cc,
std::unique_ptr<mu2e::CaloHitCollection> const& calo_hits,
std::unique_ptr<mu2e::CaloHitCollection> const& caphri_hits,
std::unique_ptr<mu2e::IntensityInfoCalo> const& int_info) {
Expand Down Expand Up @@ -310,6 +298,8 @@ void art::CaloHitsFromDataDTCEvents::analyze_calorimeter_(
mu2e::CaloRawSiPMId rawId(thisHitPacket.BoardID, thisHitPacket.ChannelID);
mu2e::CaloSiPMId offlineId = calodaqconds.offlineId(rawId);
uint16_t SiPMID = offlineId.id();
double adc2MeV = calCalib.ADC2MeV(SiPMID);
double timeOffset = calCalib.timeoffset(SiPMID);

//Check that the channel is valid in the offline world
if (SiPMID >= mu2e::CaloConst::_nChannel){
Expand All @@ -325,8 +315,8 @@ void art::CaloHitsFromDataDTCEvents::analyze_calorimeter_(

uint16_t crystalID = offlineId.crystal().id();
size_t peakIndex = thisHitPacket.IndexOfMaxDigitizerSample;
float eDep = thisHitPeak * peakADC2MeV_[SiPMID];
float time = (thisHitPacket.Time + peakIndex) * mu2e::CaloConst::_digitizationPeriod + timeCalib_[SiPMID]; // packet Time is in clock ticks
float eDep = thisHitPeak * adc2MeV;
float time = (thisHitPacket.Time + peakIndex) * mu2e::CaloConst::_digitizationPeriod - timeOffset; // packet Time is in clock ticks

bool isCaphri = offlineId.crystal().isCaphri();
if(!isCaphri) ++nhits[offlineId.crystal().disk()];
Expand Down Expand Up @@ -385,6 +375,8 @@ void art::CaloHitsFromDataDTCEvents::analyze_calorimeter_(
mu2e::CaloRawSiPMId rawId(thisHitPacket.BoardID, thisHitPacket.ChannelID);
mu2e::CaloSiPMId offlineId = calodaqconds.offlineId(rawId);
uint16_t SiPMID = offlineId.id();
double adc2MeV = calCalib.ADC2MeV(SiPMID);
double timeOffset = calCalib.timeoffset(SiPMID);

//Check that the channel is valid in the offline world
if (SiPMID >= mu2e::CaloConst::_nChannel){
Expand All @@ -400,8 +392,8 @@ void art::CaloHitsFromDataDTCEvents::analyze_calorimeter_(

uint16_t crystalID = offlineId.crystal().id();
size_t peakIndex = thisHitPacket.IndexOfMaxDigitizerSample;
float eDep = thisHitPeak * peakADC2MeV_[SiPMID];
float time = (thisHitPacket.Time + peakIndex) * mu2e::CaloConst::_digitizationPeriod + timeCalib_[SiPMID]; // packet Time is in clock ticks
float eDep = thisHitPeak * adc2MeV;
float time = (thisHitPacket.Time + peakIndex) * mu2e::CaloConst::_digitizationPeriod - timeOffset; // packet Time is in clock ticks

bool isCaphri = offlineId.crystal().isCaphri();
if(!isCaphri) ++nhits[offlineId.crystal().disk()];
Expand Down