Skip to content

Added crystalID field to caloHitMC and changed LY for LYSO #1756

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
oksuzian merged 4 commits into from
Mar 12, 2026
Merged

Added crystalID field to caloHitMC and changed LY for LYSO #1756

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
faa7503
Resolved conflict in CaloHitTruthMatch
bechenard Mar 10, 2026
e9b9527
Fixed bug and perf improvements
bechenard Mar 10, 2026
bbce729
Included Paolo wishes
bechenard Mar 10, 2026
64bab70
Fixed bug
bechenard Mar 10, 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
8 changes: 5 additions & 3 deletions CaloMC/fcl/prolog.fcl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ CaloMC : {
digiSampling : @local::HitMakerDigiSampling
nMaxFragment : 1000
minPeakADC : @local::HitMakerMinPeakADC
readoutPEPerMeV : @local::readoutPEPerMeV
readoutPEPerMeV : @local::readoutPEPerMeVCsI
MeVToADC : @local::MeVToADC
diagLevel : 0
}
Expand Down Expand Up @@ -55,7 +55,8 @@ CaloMC : { @table::CaloMC
digitizationEnd : @local::HitMakerDigitizationEnd
digitizationBuffer : 10
crystalNonUniformity : 1.05
readoutPEPerMeV : @local::readoutPEPerMeV
readoutPEPerMeVCsI : @local::readoutPEPerMeVCsI
readoutPEPerMeVLyso : @local::readoutPEPerMeVLyso
LRUCorrection : true
BirksCorrection : true
PEStatCorrection : true
Expand All @@ -81,7 +82,8 @@ CaloMC : { @table::CaloMC
minPeakADC : @local::HitMakerMinPeakADC
nBinsPeak : 2
bufferDigi : 16
readoutPEPerMeV : @local::readoutPEPerMeV
readoutPEPerMeVCsI : @local::readoutPEPerMeVCsI
readoutPEPerMeVLyso : @local::readoutPEPerMeVLyso
MeVToADC : @local::MeVToADC
diagLevel : 0
}
Expand Down
36 changes: 21 additions & 15 deletions CaloMC/src/CaloDigiMaker_module.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@
#include "Offline/RecoDataProducts/inc/CaloDigi.hh"
#include "Offline/SeedService/inc/SeedService.hh"
#include "Offline/MCDataProducts/inc/ProtonBunchTimeMC.hh"
#include "Offline/DataProducts/inc/CaloSiPMId.hh"

#include "CLHEP/Vector/ThreeVector.h"
#include "CLHEP/Random/RandPoissonQ.h"
Expand Down Expand Up @@ -60,13 +61,14 @@ namespace mu2e {
fhicl::Atom<art::InputTag> pbtmcTag { Name("protonBunchTimeMC"), Comment("ProtonBunchTimeMC producer") };
fhicl::Atom<std::string> pulseFileName { Name("pulseFileName"), Comment("Calo pulse file name") };
fhicl::Atom<std::string> pulseHistName { Name("pulseHistName"), Comment("Calo pulse hist name") };
fhicl::Atom<double> digitizationStart { Name("digitizationStart"), Comment("Start of digitization window relative to nominal pb time") };
fhicl::Atom<double> digitizationEnd { Name("digitizationEnd"), Comment("End of digitization window relative to nominal pb time")};
fhicl::Atom<float> digitizationStart { Name("digitizationStart"), Comment("Start of digitization window relative to nominal pb time") };
fhicl::Atom<float> digitizationEnd { Name("digitizationEnd"), Comment("End of digitization window relative to nominal pb time")};
fhicl::Atom<bool> addNoise { Name("addNoise"), Comment("Add noise to waveform") };
fhicl::Atom<bool> addRandomNoise { Name("addRandomNoise"), Comment("Add random salt and pepper noise") };
fhicl::Atom<double> digiSampling { Name("digiSampling"), Comment("Digitization time sampling") };
fhicl::Atom<double> pePerMeV { Name("readoutPEPerMeV"), Comment("Number of pe / MeV for Readout") };
fhicl::Atom<double> MeVToADC { Name("MeVToADC"), Comment("MeV to ADC conversion factor") };
fhicl::Atom<float> digiSampling { Name("digiSampling"), Comment("Digitization time sampling") };
fhicl::Atom<float> pePerMeVCsI { Name("readoutPEPerMeVCsI"), Comment("Number of pe / MeV for Readout for CsI") };
fhicl::Atom<float> pePerMeVLyso { Name("readoutPEPerMeVLyso"), Comment("Number of pe / MeV for Readout for LYSO") };
fhicl::Atom<float> MeVToADC { Name("MeVToADC"), Comment("MeV to ADC conversion factor") };
fhicl::Atom<int> nBits { Name("nBits"), Comment("ADC Number of bits") };
fhicl::Atom<unsigned> nBinsPeak { Name("nBinsPeak"), Comment("Window size for finding local maximum to digitize wf") };
fhicl::Atom<int> minPeakADC { Name("minPeakADC"), Comment("Minimum ADC hits of local peak to digitize") };
Expand All @@ -85,7 +87,8 @@ namespace mu2e {
bufferDigi_ (config().bufferDigi()),
startTimeBuffer_ (config().digiSampling()*config().bufferDigi()),
maxADCCounts_ ((1 << config().nBits()) -1),
pePerMeV_ (config().pePerMeV()),
pePerMeVCsI_ (config().pePerMeVCsI()),
pePerMeVLyso_ (config().pePerMeVLyso()),
MeVToADC_ (config().MeVToADC()),
pulseShape_ (CaloPulseShape(config().pulseFileName(),config().pulseHistName(),config().digiSampling())),
wfExtractor_ (config().bufferDigi(),config().nBinsPeak(),config().minPeakADC(),config().bufferDigi()),
Expand Down Expand Up @@ -122,15 +125,16 @@ namespace mu2e {
const art::ProductToken<CaloShowerROCollection> caloShowerToken_;
art::InputTag ewMarkerTag_;
art::InputTag pbtmcTag_;
double digitizationStart_;
double digitizationEnd_;
double timeFromProtonsToDRMarker_;
double digiSampling_;
float digitizationStart_;
float digitizationEnd_;
float timeFromProtonsToDRMarker_;
float digiSampling_;
unsigned bufferDigi_;
double startTimeBuffer_;
float startTimeBuffer_;
int maxADCCounts_;
double pePerMeV_;
double MeVToADC_;
float pePerMeVCsI_;
float pePerMeVLyso_;
float MeVToADC_;
CaloPulseShape pulseShape_;
CaloWFExtractor wfExtractor_;
CLHEP::HepRandomEngine& engine_;
Expand Down Expand Up @@ -227,8 +231,10 @@ namespace mu2e {
bool CaloDigiMaker::fillROHits(unsigned iRO, std::vector<double>& waveform, const CaloShowerROCollection& CaloShowerROs,
const ProtonBunchTimeMC& pbtmc)
{
bool isEmpty(true);
float scaleFactor(MeVToADC_/pePerMeV_);
bool isEmpty = true;
bool isCaphri = CaloSiPMId(iRO).crystal().isCaphri();
auto pePerMeV = isCaphri ? pePerMeVLyso_ : pePerMeVCsI_;
auto scaleFactor = MeVToADC_/pePerMeV;

for (const auto& CaloShowerRO : CaloShowerROs)
{
Expand Down
4 changes: 2 additions & 2 deletions CaloMC/src/CaloHitTruthMatch_module.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -150,8 +150,8 @@ namespace mu2e {

makeTruthMatch(event, *caloHitMCs, *caloHitMCTruth, *caloShowerMCTruth, primaryParticles);

event.put(std::move(caloHitMCTruth));
event.put(std::move(caloHitMCs));
event.put(std::move(caloHitMCTruth));
if (fillDetailedMC_) event.put(std::move(caloShowerMCTruth));
}

Expand Down Expand Up @@ -242,7 +242,7 @@ namespace mu2e {

//sort CaloEDepMC by decreasing energy, create caloHit and keep track of hit -> MChit association
std::sort(edeps.begin(),edeps.end(),[](const auto& a, const auto& b){return a.energyDep() > b.energyDep();});
caloHitMCs.emplace_back(CaloHitMC(std::move(edeps)));
caloHitMCs.emplace_back(CaloHitMC(std::move(edeps),hit.crystalID()));

art::Ptr<CaloHitMC> hitMCPtr = art::Ptr<CaloHitMC>(hitMCProductID, caloHitMCs.size()-1, hitMCProductGetter);
CaloHitTruthMatch.addSingle(hitPtr,hitMCPtr);
Expand Down
16 changes: 11 additions & 5 deletions CaloMC/src/CaloShowerROMaker_module.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -96,7 +96,8 @@ namespace mu2e {
fhicl::Atom<float> digitizationEnd { Name("digitizationEnd"), Comment("Maximum time of hit to be digitized") };
fhicl::Atom<float> digitizationBuffer { Name("digitizationBuffer"), Comment("Digi time buffer for photon propagation inside crystal") };
fhicl::Atom<float> crystalNonUniformity { Name("crystalNonUniformity"), Comment("LRU parameter 0") };
fhicl::Atom<float> pePerMeV { Name("readoutPEPerMeV"), Comment("Number of pe / MeV for Readout") };
fhicl::Atom<float> pePerMeVCsI { Name("readoutPEPerMeVCsI"), Comment("Number of pe / MeV for Readout for CsI") };
fhicl::Atom<float> pePerMeVLyso { Name("readoutPEPerMeVLyso"), Comment("Number of pe / MeV for Readout for LYSO") };
fhicl::Atom<bool> LRUCorrection { Name("LRUCorrection"), Comment("Include LRU corrections") };
fhicl::Atom<bool> BirksCorrection { Name("BirksCorrection"), Comment("Include Birks corrections") };
fhicl::Atom<bool> PEStatCorrection { Name("PEStatCorrection"), Comment("Include PE Poisson fluctuations") };
Expand All @@ -112,7 +113,8 @@ namespace mu2e {
digitizationEnd_ (config().digitizationEnd()),
digitizationBuffer_ (config().digitizationBuffer()),
crystalNonUniformity_(config().crystalNonUniformity()),
pePerMeV_ (config().pePerMeV()),
pePerMeVCsI_ (config().pePerMeVCsI()),
pePerMeVLyso_ (config().pePerMeVLyso()),
LRUCorrection_ (config().LRUCorrection()),
BirksCorrection_ (config().BirksCorrection()),
PEStatCorrection_ (config().PEStatCorrection()),
Expand Down Expand Up @@ -151,7 +153,8 @@ namespace mu2e {
float digitizationEnd_;
float digitizationBuffer_;
float crystalNonUniformity_;
float pePerMeV_;
float pePerMeVCsI_;
float pePerMeVLyso_;
bool LRUCorrection_;
bool BirksCorrection_;
bool PEStatCorrection_;
Expand Down Expand Up @@ -274,11 +277,14 @@ namespace mu2e {
if (BirksCorrection_) edep_corr = step.energyDepBirks();
if (LRUCorrection_) edep_corr = LRUCorrection(crystalID, posZ/cryhalflength, edep_corr);

bool isCaphri = CrystalId(crystalID).isCaphri();
float pePerMeV = isCaphri ? pePerMeVLyso_ : pePerMeVCsI_;

// Generate individual PEs and their arrival times
for (int i=0; i<nROs; ++i)
{
int SiPMID = SiPMIDBase + i;
int NPE = randPoisson_.fire(edep_corr*pePerMeV_);
int NPE = randPoisson_.fire(edep_corr*pePerMeV);
if (NPE==0) continue;

std::vector<float> PETime(NPE,hitTime);
Expand All @@ -293,7 +299,7 @@ namespace mu2e {
if (diagLevel_ > 1) for (const auto& time : PETime) hTime_->Fill(2.0*cryhalflength-posZ,time-hitTime);

diagSum.totNPE += NPE;
diagSum.totEdepNPE += double(NPE)/pePerMeV_/2.0; //average between the two RO
diagSum.totEdepNPE += double(NPE)/pePerMeV/2.0; //average between the two RO
}
diagSum.totEdep += step.energyDepG4();
diagSum.totEdepCorr += edep_corr;
Expand Down
3 changes: 2 additions & 1 deletion CaloReco/fcl/common.fcl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,8 @@ BEGIN_PROLOG
HitMakerMinPeakADC : 16
pulseFileName : "OfflineData/ConditionsService/CsI-waveform-2020-02-12.root"
pulseHistName : "h_waveform"
readoutPEPerMeV : 30.0
readoutPEPerMeVCsI : 30.0
readoutPEPerMeVLyso : 500.0
MeVToADC : 16.0
ADCToMeV : 0.0625
END_PROLOG
Expand Down
13 changes: 11 additions & 2 deletions MCDataProducts/inc/CaloHitMC.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,8 +18,15 @@ namespace mu2e
class CaloHitMC
{
public:
CaloHitMC() : edeps_() {};
CaloHitMC(const std::vector<CaloEDepMC>& edeps) : edeps_(edeps) {};
CaloHitMC() : edeps_(), crystalID_(-1) {};
CaloHitMC(const std::vector<CaloEDepMC>& edeps, int crystalId) :
edeps_(edeps),
crystalID_(crystalId)
{};
CaloHitMC(const std::vector<CaloEDepMC>&& edeps, int crystalId) :
edeps_(std::move(edeps)),
crystalID_(crystalId)
{};

void resetSim(SimParticleRemapping const& remap);

Expand All @@ -30,9 +37,11 @@ namespace mu2e
float time () const {return edeps_.empty() ? 0.0 : edeps_.at(0).time();}
float totalEnergyDep () const;
float totalEnergyDepG4() const;
int crystalID () const { return crystalID_;}

private:
std::vector<CaloEDepMC> edeps_;
int crystalID_;
};

using CaloHitMCCollection = std::vector<mu2e::CaloHitMC>;
Expand Down
2 changes: 1 addition & 1 deletion MCDataProducts/inc/CaloMCTruthAssns.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@
namespace mu2e
{
using CaloShowerMCTruthAssn = art::Assns<CaloHit, SimParticle, art::Ptr<CaloShowerSim>>;
using CaloHitMCTruthAssn = art::Assns<CaloHit, CaloHitMC> ;
using CaloHitMCTruthAssn = art::Assns<CaloHit, CaloHitMC> ;
using CaloClusterMCTruthAssn = art::Assns<CaloCluster,CaloClusterMC> ;
}

Expand Down