Skip to content
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.

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

Implement digi mixing for tracker #1245

Merged
merged 13 commits into from
May 24, 2024
Merged
Show file tree
Hide file tree
Changes from 2 commits
Commits
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
52 changes: 52 additions & 0 deletions CommonMC/src/ProtonBunchTimeMCFromProtonBunchTime_module.cc
Original file line number Diff line number Diff line change
@@ -0,0 +1,52 @@
// Convert proton bunch time marker from data to MC data product
kutschke marked this conversation as resolved.
Show resolved Hide resolved
// To be used for mixing simulated signal into real data from detector
//
// Ed Callaghan, 2024

// art
#include <art/Framework/Core/EDProducer.h>
#include <art/Framework/Principal/Event.h>

// canvas
#include <canvas/Utilities/InputTag.h>

// fhiclcpp
#include <fhiclcpp/types/Atom.h>

// mu2e
#include <Offline/MCDataProducts/inc/ProtonBunchTimeMC.hh>
#include <Offline/RecoDataProducts/inc/ProtonBunchTime.hh>

namespace mu2e{
class ProtonBunchTimeMCFromProtonBunchTime: public art::EDProducer{
public:
struct Config{
fhicl::Atom<art::InputTag> pbttag{
fhicl::Name("PBTTag"),
fhicl::Comment("ProtonBunchTime provenance"),
"" // default, eventually to specify "from data"
};
};

using Parameters = art::EDProducer::Table<Config>;
explicit ProtonBunchTimeMCFromProtonBunchTime(Parameters const& config);
protected:
/**/
private:
void produce(art::Event& event) override;
art::InputTag pbttag_;
};

ProtonBunchTimeMCFromProtonBunchTime::ProtonBunchTimeMCFromProtonBunchTime(Parameters const& config): EDProducer(config), pbttag_(config().pbttag()){
this->produces<ProtonBunchTimeMC>();
}

void ProtonBunchTimeMCFromProtonBunchTime::produce(art::Event& event){
auto pbt = event.getHandle<ProtonBunchTime>(pbttag_);
std::unique_ptr<ProtonBunchTimeMC> pbtmc(new ProtonBunchTimeMC);
pbtmc->pbtime_ = pbt->pbtime_;
event.put(std::move(pbtmc));
}
}

DEFINE_ART_MODULE(mu2e::ProtonBunchTimeMCFromProtonBunchTime);
25 changes: 24 additions & 1 deletion EventMixing/inc/Mu2eProductMixer.hh
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,9 @@
#include "Offline/MCDataProducts/inc/CosmicLivetime.hh"
#include "Offline/MCDataProducts/inc/SimTimeOffset.hh"
#include "Offline/MCDataProducts/inc/PhysicalVolumeInfoMultiCollection.hh"

#include "Offline/RecoDataProducts/inc/StrawDigi.hh"
#include "Offline/MCDataProducts/inc/StrawDigiMC.hh"
#include "Offline/DataProducts/inc/EventWindowMarker.hh"


//================================================================
Expand Down Expand Up @@ -102,6 +104,10 @@ namespace mu2e {
fhicl::Table<CollectionMixerConfig> crvStepMixer { fhicl::Name("crvStepMixer") };
fhicl::Table<CollectionMixerConfig> extMonSimHitMixer { fhicl::Name("extMonSimHitMixer") };
fhicl::Table<CollectionMixerConfig> eventIDMixer { fhicl::Name("eventIDMixer") };
fhicl::Table<CollectionMixerConfig> strawDigiMixer { fhicl::Name("strawDigiMixer") };
fhicl::Table<CollectionMixerConfig> strawDigiADCWaveformMixer { fhicl::Name("strawDigiADCWaveformMixer") };
fhicl::Table<CollectionMixerConfig> strawDigiMCMixer { fhicl::Name("strawDigiMCMixer") };
fhicl::Table<CollectionMixerConfig> eventWindowMarkerMixer { fhicl::Name("eventWindowMarkerMixer") };
fhicl::OptionalTable<CosmicLivetimeMixerConfig> cosmicLivetimeMixer { fhicl::Name("cosmicLivetimeMixer") };
fhicl::OptionalTable<VolumeInfoMixerConfig> volumeInfoMixer { fhicl::Name("volumeInfoMixer") };
fhicl::OptionalAtom<art::InputTag> simTimeOffset { fhicl::Name("simTimeOffset"), fhicl::Comment("Simulation time offset to apply (optional)") };
Expand Down Expand Up @@ -148,6 +154,23 @@ namespace mu2e {
ExtMonFNALSimHitCollection& out,
art::PtrRemapper const& remap);

bool mixStrawDigis(std::vector<StrawDigiCollection const*> const& in,
StrawDigiCollection& out,
art::PtrRemapper const& remap);

bool mixStrawDigiADCWaveforms(std::vector<StrawDigiADCWaveformCollection const*> const& in,
StrawDigiADCWaveformCollection& out,
art::PtrRemapper const& remap);

bool mixStrawDigiMCs(std::vector<StrawDigiMCCollection const*> const& in,
StrawDigiMCCollection& out,
art::PtrRemapper const& remap);

bool mixEventWindowMarkers(std::vector<EventWindowMarker const*> const& in,
EventWindowMarker& out,
art::PtrRemapper const& remap);


bool mixEventIDs(std::vector<art::EventIDSequence const*> const &in,
art::EventIDSequence& out,
art::PtrRemapper const& remap);
Expand Down
113 changes: 113 additions & 0 deletions EventMixing/src/MixDigis_module.cc
Original file line number Diff line number Diff line change
@@ -0,0 +1,113 @@
// Mix digi-level data products, as opposed to simulation-level data products,
// into art events
//
// Ed Callaghan, 2023

// stl

// art
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/SubRun.h"
#include "art/Framework/IO/ProductMix/MixHelper.h"
#include "art/Framework/Modules/MixFilter.h"
#include "art_root_io/RootIOPolicy.h"
#include "fhiclcpp/types/Atom.h"
#include "fhiclcpp/types/Comment.h"
#include "fhiclcpp/types/Name.h"
#include "fhiclcpp/types/Sequence.h"
#include "fhiclcpp/types/Table.h"
#include "fhiclcpp/types/TupleAs.h"

// cetlib_except
#include "cetlib_except/exception.h"

// mu2e
#include "Offline/EventMixing/inc/Mu2eProductMixer.hh"

namespace mu2e{
// detail class for art::MixFilter<>
class MixDigisDetail{
public:
struct Mu2eConfig{
// configuration is simple --- which data products to mix in
fhicl::Table<Mu2eProductMixer::Config> products{
fhicl::Name("products"),
fhicl::Comment("Products to be mixed, as form of a mixingMap, i.e. tuples of InputTags to output instance names.")

};
};

struct Config{
fhicl::Table<Mu2eConfig> mu2e{ fhicl::Name("mu2e") };
};

using Parameters = art::MixFilterTable<Config>;
explicit MixDigisDetail(const Parameters& parameters,
art::MixHelper& helper);
size_t nSecondaries();

void processEventIDs(const art::EventIDSequence& seq);
void beginSubRun(const art::SubRun& subrun);
void endSubRun(art::SubRun& subrun);
void startEvent(const art::Event& event);
void finalizeEvent(art::Event& event);

protected:
/**/

private:
Mu2eProductMixer mixer;
rlcee marked this conversation as resolved.
Show resolved Hide resolved
art::EventIDSequence evids;
};

// implementation
MixDigisDetail::MixDigisDetail(const Parameters& parameters,
art::MixHelper& helper)
: mixer{parameters().mu2e().products(), helper}{
helper.produces<art::EventIDSequence>();
}

// always overlay simulated event onto one event window
size_t MixDigisDetail::nSecondaries(){
size_t rv = 1;
return rv;
}

// always keep track of secondary EventIDs
void MixDigisDetail::processEventIDs(const art::EventIDSequence& seq){
// intent is to overlay signals onto a nominally-complete event,
// so we should only mix in one event
if (seq.size() != 1){
throw cet::exception("MIX") << "mu2e::DixDigisDetail: mixing more than one digi frame" << std::endl;
}
mixer.processEventIDs(seq);
evids = seq;
}

// forward beginSubRun
void MixDigisDetail::beginSubRun(const art::SubRun& subrun){
mixer.beginSubRun(subrun);
}

// forward endSubRun
void MixDigisDetail::endSubRun(art::SubRun& subrun){
mixer.endSubRun(subrun);
}

// forward startEvent
void MixDigisDetail::startEvent(const art::Event& event){
mixer.startEvent(event);
}

void MixDigisDetail::finalizeEvent(art::Event& event){
// "manually" mix in EventIDs, as these do not preexist as data products
auto seq = std::make_unique<art::EventIDSequence>();
seq->swap(evids);
event.put(std::move(seq));
}

// define the module
typedef art::MixFilter<MixDigisDetail,art::RootIOPolicy> MixDigis;
}

DEFINE_ART_MODULE(mu2e::MixDigis);
58 changes: 58 additions & 0 deletions EventMixing/src/Mu2eProductMixer.cc
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@
#include <algorithm>
#include <iterator>
#include <iostream>
#include <string>

#include "cetlib_except/exception.h"

Expand Down Expand Up @@ -96,6 +97,26 @@ namespace mu2e {
(e.inTag, e.resolvedInstanceName(), &Mu2eProductMixer::mixEventIDs, *this);
}

for(const auto& e: conf.strawDigiMixer().mixingMap()) {
helper.declareMixOp
(e.inTag, e.resolvedInstanceName(), &Mu2eProductMixer::mixStrawDigis, *this);
}

for(const auto& e: conf.strawDigiADCWaveformMixer().mixingMap()) {
helper.declareMixOp
(e.inTag, e.resolvedInstanceName(), &Mu2eProductMixer::mixStrawDigiADCWaveforms, *this);
}

for(const auto& e: conf.strawDigiMCMixer().mixingMap()) {
helper.declareMixOp
(e.inTag, e.resolvedInstanceName(), &Mu2eProductMixer::mixStrawDigiMCs, *this);
}

for(const auto& e: conf.eventWindowMarkerMixer().mixingMap()) {
helper.declareMixOp
(e.inTag, e.resolvedInstanceName(), &Mu2eProductMixer::mixEventWindowMarkers, *this);
}

//----------------------------------------------------------------
// VolumeInfo handling

Expand Down Expand Up @@ -397,6 +418,43 @@ namespace mu2e {
return true;
}

bool Mu2eProductMixer::mixStrawDigis(std::vector<StrawDigiCollection const*> const& in,
StrawDigiCollection& out,
art::PtrRemapper const& remap)
{
art::flattenCollections(in, out);
return true;
}

bool Mu2eProductMixer::mixStrawDigiADCWaveforms(std::vector<StrawDigiADCWaveformCollection const*> const& in,
StrawDigiADCWaveformCollection& out,
art::PtrRemapper const& remap)
{
art::flattenCollections(in, out);
return true;
}

bool Mu2eProductMixer::mixStrawDigiMCs(std::vector<StrawDigiMCCollection const*> const& in,
StrawDigiMCCollection& out,
art::PtrRemapper const& remap)
{
art::flattenCollections(in, out);
return true;
}

bool Mu2eProductMixer::mixEventWindowMarkers(std::vector<EventWindowMarker const*> const& in,
EventWindowMarker& out,
art::PtrRemapper const& remap){
// assert that only one EventWindowMarker be present
if (in.size() != 1){
std::string msg = "attempting to mix more than 1 EventWindowMarker: ";
msg += std::to_string(in.size()) + " present";
throw cet::exception("BADINPUT") << msg << std::endl;
}
out = *in[0];
return true;
}

//----------------------------------------------------------------
bool Mu2eProductMixer::mixEventIDs(std::vector<art::EventIDSequence const*> const &in,
art::EventIDSequence& out,
Expand Down
8 changes: 7 additions & 1 deletion MCDataProducts/inc/StrawDigiMC.hh
Original file line number Diff line number Diff line change
Expand Up @@ -32,16 +32,21 @@ namespace mu2e {
typedef std::array<SGSP,StrawEnd::nends> SGSPA;
typedef std::array<XYZVectorF,StrawEnd::nends> PA;
typedef std::array<float,StrawEnd::nends> FA;
enum Validity {Valid, PartiallyValid, Invalid};
rlcee marked this conversation as resolved.
Show resolved Hide resolved
Copy link
Collaborator

Choose a reason for hiding this comment

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

I suggest renaming this enum 'Provenance', with states "sim", "mixed", "external"

Copy link
Collaborator

Choose a reason for hiding this comment

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

Should we use an EnumToString so that we have access to printed values?

Copy link
Collaborator

Choose a reason for hiding this comment

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

Copy link
Collaborator

@kutschke kutschke May 22, 2024

Choose a reason for hiding this comment

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

@brownd1978 why "external" and not "data"?

@edcallaghan There will be two ways to make pure sim digis. One is our present way and the other is to make background frames of digis from sim and then to overlay signal sims. Will both be of Provenance "sim" (in Dave's proposed language?). If so, should we distinguish them?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

IMO the "most correct" answer depends on what the conflict resolution looks like for the MC truth data products, which I haven't through through in at a detailed level. If there are no caveats to "merging" the MC truth info for the two constituent sim digis, then labeling the single output digi as "Sim" is fine. But if there is some kind of catch, then we probably want to label is as "Sim/Sim" as opposed to e.g. "Sim" or "Sim/Data".

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Oh, but also, this is not limited to just two digis input. It won't occur at scale, but in principle arbitrarily many digis can participate. So I'd err on the side of Dave's suggestion of renaming the three situations (purely new MC, mixed new MC with something preexisting, and purely something preexisting). In the case of the second situation, the onus is then on the user to know what it is that they have mixed (preexisting simulation, or real data), which I think is a reasonable assumption if they are looking at MC truth info.

Copy link
Collaborator

Choose a reason for hiding this comment

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

I think we should distinguish 3 different kinds of mixed: (sim+sim), (sim+data) and (data+data). Maybe this is not the right place to record that information - but we do need a plan for a robust audit trail so that an end user can determine what it is programatically and not by looking at a wiki page or a log book.

An example of data+data would be if we took at track from Offspill and embedded it in a random on-spill trigger as part of a study to understand reconstruction efficiency.

Copy link
Contributor Author

@edcallaghan edcallaghan May 22, 2024

Choose a reason for hiding this comment

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

I agree that that information (i.e. the mixing constituents/identities) should, for official production, be stored somewhere. I do not have a strong opinion on whether it should be part of the data structure or not. But I do not think that this class is the correct place for that book-keeping. For one thing, if someone is looking at a StrawDigiMC, then they are either going down the wrong path to begin with, or they are aware that they have run a simulation, and just need to know if the MC truth available gives the full story of the associated digi.

Also, as a breadcrumb for any future readers, because the actual "merging" of data products is done within the StrawDigisFromStrawGasSteps simulation module, this PR cannot create a single collection of digis from a "data + data" situation. Though, as Rob sas, we may want to accommodate that for future use cases --- in which case, the implementation will need to be developed.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

@brownd1978 @kutschke The latest push exchanges Validity::{Valid, PartiallyValid, Invalid} for Provenance::{unknown, Simulation, Mixed, External}, which is implemented as an EnumToStringSparse specialization. The existence of unknown is forced by that interface. In the latter case, the onus is on the user to understand that External may either signify that real data or preexisting MC (with or without MC info available yet undetermined, as that will depend on how the collision resolution ends up being performed). But the use case is really just to flag that one should not expect the MC truth object to actually represent MC truth information, as the digi came from an alternate source.

As to resolving the ambiguity of whether it was data or simulation which was provided externally, my current thinking is as follows: the EventID of the mixed event is required to be propagated as a data product so that the conditions can be matched. As long as we don't drop that data product from the output, this would satisfy the book-keeping (i.e., the run number of the mixed event tells us an analyst if it came from data or not).

Copy link
Collaborator

Choose a reason for hiding this comment

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

Thanks Ed. We have converged as well as we can for now. Your suggestion for documenting that data was the external might well work - we can work through the details when the time comes.


StrawDigiMC();
// construct from hitlets
StrawDigiMC(StrawId sid, PA cpos, FA ctime, FA wetime, SGSPA sgs);
StrawDigiMC(StrawId sid, PA cpos, FA ctime, FA wetime, SGSPA sgs, Validity=Valid);

// use compuater copy construcors
StrawDigiMC(const StrawDigiMC& rhs, SGSPA sgsp ); // update the Ptrs
StrawDigiMC(const StrawDigiMC& rhs, Validity validity ); // update validity

// Accessors
StrawId strawId() const { return _strawid; }

Validity validity() const { return _validity; }

SGSP const& strawGasStep(StrawEnd strawend) const { return _sgspa[strawend]; }
SGSPA const& strawGasSteps() const { return _sgspa; }
StrawEnd earlyEnd() const { return (_wtime[StrawEnd::cal] < _wtime[StrawEnd::hv]) ? StrawEnd::cal : StrawEnd::hv; }
Expand All @@ -68,6 +73,7 @@ namespace mu2e {
FA _ctime; // times of the trigger clusters
FA _wtime; // times at the wire ends of the signals which fired the TDC.
SGSPA _sgspa; // StrawGasStep that triggered each end
Validity _validity; // level of association with any true MC events
};

inline std::ostream& operator<<( std::ostream& ost,
Expand Down
8 changes: 7 additions & 1 deletion MCDataProducts/src/StrawDigiMC.cc
Original file line number Diff line number Diff line change
Expand Up @@ -18,16 +18,22 @@ namespace mu2e {
// Default constructor is required for persistable classes
StrawDigiMC::StrawDigiMC()
: _strawid(StrawId::_invalid)
, _validity(Invalid)
{}

StrawDigiMC::StrawDigiMC(StrawId sid, PA cpos, FA ctime, FA wetime, SGSPA sgs):
StrawDigiMC::StrawDigiMC(StrawId sid, PA cpos, FA ctime, FA wetime, SGSPA sgs, Validity validity):
_strawid(sid), _cpos(cpos), _ctime(ctime), _wtime(wetime), _sgspa(sgs)
, _validity(validity)
{}

StrawDigiMC::StrawDigiMC(const StrawDigiMC& rhs, SGSPA sgspa ) : StrawDigiMC(rhs) {
_sgspa = sgspa;
}

StrawDigiMC::StrawDigiMC(const StrawDigiMC& rhs, Validity validity): StrawDigiMC(rhs){
_validity = validity;
}

bool StrawDigiMC::isCrossTalk(StrawEnd strawend) const {
bool retval(false);
if(!_sgspa[strawend].isNull()){
Expand Down
35 changes: 35 additions & 0 deletions TrackerMC/inc/StrawDigiBundle.hh
Original file line number Diff line number Diff line change
@@ -0,0 +1,35 @@
// Bundle together references to various StrawDigi products, to circumvent
// triplicated-loops which pop up when iterating over digis
//
// Ed Callaghan, 2023

#ifndef TrackerMC_StrawDigiBundle_hh
#define TrackerMC_StrawDigiBundle_hh

// mu2e
#include <Offline/RecoDataProducts/inc/StrawDigi.hh>
#include <Offline/MCDataProducts/inc/StrawDigiMC.hh>

namespace mu2e{
class StrawDigiBundle{
public:
StrawDigiBundle(const StrawDigi digi, const StrawDigiADCWaveform adcs, const StrawDigiMC mc);
StrawDigiBundle(const StrawDigi digi, const StrawDigiADCWaveform adcs);
StrawDigiBundle(const StrawDigiBundle& bundle);

const StrawDigi GetStrawDigi() const;
rlcee marked this conversation as resolved.
Show resolved Hide resolved
const StrawDigiADCWaveform GetStrawDigiADCWaveform() const;
const StrawDigiMC GetStrawDigiMC() const;

// interface for sorting into buckets of overlapping digitization windows
const double time() const;
protected:
const StrawDigi digi;
Copy link
Collaborator

Choose a reason for hiding this comment

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

this payload is deep-copy. Does the content come from the event? If so, (const) reference payload would avoid copies.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

There are two reasons that the members are values instead of references. The first is that, to use references, the underlying types (specifically the waveforms, IIRC) would need non-trivial copy constructors implemented and some of the loops need to be restructured to deal with the limitations of constructing objects with reference members. Those are hurdles which are both completely surmountable and, IMO, annoying.

The second reason is that, while right now all of the digis do come "from an event," this will not be true in the future: when we implement "real" collision resolution, e.g. by redigitizing a reconstructed analog sum, then this collection will be responsible for transferring entirely new data products into the scope of the simulation module, and hence should be deep-copied. But as @rlcee pointed out above, the getters can at least return references.

Copy link
Collaborator

Choose a reason for hiding this comment

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

Looking ahead, will some of these data members be purely input and others output? If so, I suggest holding inputs by const& and then outputs must be by value.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

They will all be used bidirectionally. That could change if, in the future, the collision resolution is factored out of StrawDigisFromStrawGasSteps and into some standalone module.

const StrawDigiADCWaveform adcs;
const StrawDigiMC mc;
private:
/**/
};
}

#endif
Loading