Skip to content

Commit

Permalink
Merge pull request #1291 from edcallaghan/update_TrkStrawHitMC_proven…
Browse files Browse the repository at this point in the history 
…ance

Add TrkStrawHitMC provenance and guard against empty StrawDigiMCs
  • Loading branch information
@kutschke
kutschke authored Jul 1, 2024
2 parents 465477c + 4ead963 commit 34a99ea
Show file tree
Hide file tree
Showing 14 changed files with 247 additions and 162 deletions.
192 changes: 103 additions & 89 deletions CommonMC/src/SelectRecoMC_module.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,6 +53,7 @@
#include <iostream>
#include <set>
#include <string>
#include <algorithm>

namespace mu2e {
class SelectRecoMC : public art::EDProducer {
Expand Down Expand Up @@ -233,14 +234,17 @@ namespace mu2e {
// find the referenced sim particle
int spref(-1);
auto const& sdmc = sdmcc.at(hit.index()); // bounds-check for security;
for(size_t isp=0;isp < spcc.size(); isp++){
auto const& spc = spcc[isp];
if(sdmc.earlyStrawGasStep()->simParticle() == spc._spp){
spref = isp;
break;
// if mc info is not meaningful, do not try to inspect any SimParticles
if (sdmc.containsSimulation()){
for(size_t isp=0;isp < spcc.size(); isp++){
auto const& spc = spcc[isp];
if(sdmc.earlyStrawGasStep()->simParticle() == spc._spp){
spref = isp;
break;
}
}
if(spref < 0)throw cet::exception("Reco")<<"mu2e::SelectRecoMC: missing index"<< std::endl;
}
if(spref < 0)throw cet::exception("Reco")<<"mu2e::SelectRecoMC: missing index"<< std::endl;
TrkStrawHitMC tshmc;
fillTSHMC(tshmc,hit.index(),spref,sdmc,tracker,srep);
mcseed._tshmcs.push_back(tshmc);
Expand All @@ -250,26 +254,30 @@ namespace mu2e {
void SelectRecoMC::fillUnusedTSHMC( SPCC const& spcc, StrawDigiMCCollection const& sdmcc,
Tracker const& tracker, std::shared_ptr<const StrawResponse>const& srep,
KalSeedMC& mcseed) {
// either keep hits only from the primary or from all contributing
size_t ispmax = _saveallunused? spcc.size() : 1;
// either keep hits only from the primary or all contributing particles,
// if any exist. if they do not, then there is no primary to default to.
size_t ispmax = _saveallunused? spcc.size() : std::min<int>(spcc.size(), 1);
for(size_t isp=0; isp < ispmax; ++isp){
auto const& spc = spcc[isp];
for (size_t isdmc=0; isdmc < sdmcc.size(); isdmc++){
auto const& sdmc = sdmcc[isdmc];
auto const& sgs = *(sdmc.earlyStrawGasStep());
if(sgs.simParticle() == spc._spp){
// search to see if the associated digi is already on the track
bool used(false);
for(auto const& tshmc : mcseed._tshmcs ) {
if(isdmc == tshmc.strawDigiMCIndex()){
used = true;
break;
// if this contains no MC information, then we cannot inspect it further
if (sdmc.containsSimulation()){
auto const& sgs = *(sdmc.earlyStrawGasStep());
if(sgs.simParticle() == spc._spp){
// search to see if the associated digi is already on the track
bool used(false);
for(auto const& tshmc : mcseed._tshmcs ) {
if(isdmc == tshmc.strawDigiMCIndex()){
used = true;
break;
}
}
if(!used){
TrkStrawHitMC tshmc;
fillTSHMC(tshmc,isdmc,isp,sdmc,tracker,srep);
mcseed._tshmcs.push_back(tshmc);
}
}
if(!used){
TrkStrawHitMC tshmc;
fillTSHMC(tshmc,isdmc,isp,sdmc,tracker,srep);
mcseed._tshmcs.push_back(tshmc);
}
}
}
Expand All @@ -278,78 +286,84 @@ namespace mu2e {

void SelectRecoMC::fillTSHMC(TrkStrawHitMC& tshmc, size_t isdmc, size_t isp,StrawDigiMC const& sdmc,
Tracker const& tracker, std::shared_ptr<const StrawResponse>const& srep ) {
// always record the index, to avoid downstream corruptions
tshmc._sdmcindex = isdmc;
tshmc._spindex = isp;
tshmc._energySum = sdmc.triggerEnergySum(sdmc.earlyEnd());
const auto& sgs = *(sdmc.earlyStrawGasStep());
tshmc._cpos = XYZVectorF(sdmc.clusterPosition(sdmc.earlyEnd()));
tshmc._mom = sgs.momentum();
tshmc._time = sgs.time();
if (_onSpill){
tshmc._time = fmod(tshmc._time,_mbtime);
// fix for DAQ wrapping
if(tshmc._time < -_pbtimemc)tshmc._time += _mbtime;
}
tshmc._strawId = sdmc.strawId();
tshmc._earlyend = sdmc.earlyEnd();
// compute the signal propagation time and drift time
double vprop = 2.0*srep->halfPropV(sdmc.strawId(),1000.0*tshmc._energySum);
auto const& straw = tracker.straw(sdmc.strawId());
auto wdir = XYZVectorF(straw.wireDirection());
auto tdir = sgs.momentum().Unit();
double pdist = (straw.wireEnd(sdmc.earlyEnd())-sdmc.clusterPosition(sdmc.earlyEnd())).dot(straw.wireDirection());
tshmc._tprop = fabs(pdist)/vprop;
tshmc._tdrift = sdmc.wireEndTime(sdmc.earlyEnd()) -tshmc._time - tshmc._tprop - _pbtimemc - 2.4; // temporary kludge offset FIXME!
if (_onSpill)
tshmc._tdrift = fmod(tshmc._tdrift,_mbtime);
const static XYZVectorF bdir(0.0,0.0,1.0);
// propagate interpretability of StrawDigiMC to TrkStrawHitMC
tshmc._provenance = sdmc.provenance();
// if there is no MC information at all, leave TrkStrawHitMC empty
if (tshmc.containsSimulation()){
tshmc._spindex = isp;
tshmc._energySum = sdmc.triggerEnergySum(sdmc.earlyEnd());
const auto& sgs = *(sdmc.earlyStrawGasStep());
tshmc._cpos = XYZVectorF(sdmc.clusterPosition(sdmc.earlyEnd()));
tshmc._mom = sgs.momentum();
tshmc._time = sgs.time();
if (_onSpill){
tshmc._time = fmod(tshmc._time,_mbtime);
// fix for DAQ wrapping
if(tshmc._time < -_pbtimemc)tshmc._time += _mbtime;
}
tshmc._strawId = sdmc.strawId();
tshmc._earlyend = sdmc.earlyEnd();
// compute the signal propagation time and drift time
double vprop = 2.0*srep->halfPropV(sdmc.strawId(),1000.0*tshmc._energySum);
auto const& straw = tracker.straw(sdmc.strawId());
auto wdir = XYZVectorF(straw.wireDirection());
auto tdir = sgs.momentum().Unit();
double pdist = (straw.wireEnd(sdmc.earlyEnd())-sdmc.clusterPosition(sdmc.earlyEnd())).dot(straw.wireDirection());
tshmc._tprop = fabs(pdist)/vprop;
tshmc._tdrift = sdmc.wireEndTime(sdmc.earlyEnd()) -tshmc._time - tshmc._tprop - _pbtimemc - 2.4; // temporary kludge offset FIXME!
if (_onSpill)
tshmc._tdrift = fmod(tshmc._tdrift,_mbtime);
const static XYZVectorF bdir(0.0,0.0,1.0);

TwoLinePCA_XYZ wirepca( XYZVectorF(straw.wirePosition()), XYZVectorF(straw.wireDirection()),
tshmc._cpos, tdir );
TwoLinePCA_XYZ strawpca( XYZVectorF(straw.strawPosition()), XYZVectorF(straw.strawDirection()),
tshmc._cpos, tdir );
TwoLinePCA_XYZ wirepca( XYZVectorF(straw.wirePosition()), XYZVectorF(straw.wireDirection()),
tshmc._cpos, tdir );
TwoLinePCA_XYZ strawpca( XYZVectorF(straw.strawPosition()), XYZVectorF(straw.strawDirection()),
tshmc._cpos, tdir );

// vector from middle of wire to POCA
auto wireMid_to_POCA = wirepca.point2() - XYZVectorF(straw.wirePosition());
// longitudinal distance along wire to POCA
tshmc._wireLen = wireMid_to_POCA.Dot(wdir);
// perpedicular vector from wire to POCA
auto POCA_delta = wirepca.point2() - wirepca.point1();
// define DOCA to be negative in direction of track cross wire
auto doca_sign_dir = tdir.Cross(wdir).Unit();
tshmc._wireDOCA = -1*doca_sign_dir.Dot(POCA_delta);
// phi is defined by vector to POCA, with 0 in 'Bdir'
// and +pi/2 in V dir (radially out)
auto Vdir = wdir.Cross(bdir);
// ensure Vdir is pointing radially out
if (Vdir.Dot(XYZVectorF(straw.wirePosition())) < 0.0) Vdir *= -1.0;
tshmc._wirePhi = atan2(POCA_delta.Dot(Vdir),POCA_delta.Dot(bdir));
// rdrift is the expected T2D given tdrift and phi
tshmc._rdrift = srep->strawDrift().T2D(tshmc._tdrift,tshmc._wirePhi);
// wireDot is cos angle between track and wire
tshmc._wireDot = tdir.Dot(wdir);
// wireTau is delta distance perpedicular to the wire from POCA to trigger cluster
// first define unit vector perpendicular to particle track and wire direction
auto track_cross_wire = tdir.Cross(wdir).Unit();
// then this cross wdir is the vector perpendicular to the wire and the vector to POCA
auto wire_cross_delta = wdir.Cross(track_cross_wire);
// relative position of the trigger cluster
auto wireMid_to_cluster = tshmc._cpos-XYZVectorF(straw.wirePosition());
// wireTau is then the relative position of the trigger cluster in this direction
tshmc._wireTau = wireMid_to_cluster.Dot(wire_cross_delta);
// vector from middle of wire to POCA
auto wireMid_to_POCA = wirepca.point2() - XYZVectorF(straw.wirePosition());
// longitudinal distance along wire to POCA
tshmc._wireLen = wireMid_to_POCA.Dot(wdir);
// perpedicular vector from wire to POCA
auto POCA_delta = wirepca.point2() - wirepca.point1();
// define DOCA to be negative in direction of track cross wire
auto doca_sign_dir = tdir.Cross(wdir).Unit();
tshmc._wireDOCA = -1*doca_sign_dir.Dot(POCA_delta);
// phi is defined by vector to POCA, with 0 in 'Bdir'
// and +pi/2 in V dir (radially out)
auto Vdir = wdir.Cross(bdir);
// ensure Vdir is pointing radially out
if (Vdir.Dot(XYZVectorF(straw.wirePosition())) < 0.0) Vdir *= -1.0;
tshmc._wirePhi = atan2(POCA_delta.Dot(Vdir),POCA_delta.Dot(bdir));
// rdrift is the expected T2D given tdrift and phi
tshmc._rdrift = srep->strawDrift().T2D(tshmc._tdrift,tshmc._wirePhi);
// wireDot is cos angle between track and wire
tshmc._wireDot = tdir.Dot(wdir);
// wireTau is delta distance perpedicular to the wire from POCA to trigger cluster
// first define unit vector perpendicular to particle track and wire direction
auto track_cross_wire = tdir.Cross(wdir).Unit();
// then this cross wdir is the vector perpendicular to the wire and the vector to POCA
auto wire_cross_delta = wdir.Cross(track_cross_wire);
// relative position of the trigger cluster
auto wireMid_to_cluster = tshmc._cpos-XYZVectorF(straw.wirePosition());
// wireTau is then the relative position of the trigger cluster in this direction
tshmc._wireTau = wireMid_to_cluster.Dot(wire_cross_delta);

auto sdir = XYZVectorF(straw.strawDirection());
// perpedicular vector from straw center to POCA
auto sPOCA_delta = strawpca.point2() - strawpca.point1();
// define DOCA to be negative in direction of track cross straw
auto sdoca_sign_dir = tdir.Cross(sdir).Unit();
tshmc._strawDOCA = -1*sdoca_sign_dir.Dot(sPOCA_delta);
// phi is defined by vector to POCA, with 0 in 'Bdir'
// and +pi/2 in V dir (radially out)
auto sVdir = sdir.Cross(bdir);
// ensure Vdir is pointing radially out
if (sVdir.Dot(XYZVectorF(straw.strawPosition())) < 0.0) sVdir *= -1.0;
tshmc._strawPhi = atan2(sPOCA_delta.Dot(sVdir),sPOCA_delta.Dot(bdir));
auto sdir = XYZVectorF(straw.strawDirection());
// perpedicular vector from straw center to POCA
auto sPOCA_delta = strawpca.point2() - strawpca.point1();
// define DOCA to be negative in direction of track cross straw
auto sdoca_sign_dir = tdir.Cross(sdir).Unit();
tshmc._strawDOCA = -1*sdoca_sign_dir.Dot(sPOCA_delta);
// phi is defined by vector to POCA, with 0 in 'Bdir'
// and +pi/2 in V dir (radially out)
auto sVdir = sdir.Cross(bdir);
// ensure Vdir is pointing radially out
if (sVdir.Dot(XYZVectorF(straw.strawPosition())) < 0.0) sVdir *= -1.0;
tshmc._strawPhi = atan2(sPOCA_delta.Dot(sVdir),sPOCA_delta.Dot(bdir));
}
}

void SelectRecoMC::fillSDMCI(KalSeedMC const& mcseed, SDIS& sdindices) {
Expand Down
33 changes: 18 additions & 15 deletions Filters/src/StrawDigiMCFilter_module.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -66,21 +66,24 @@ namespace mu2e {
// keep count of digis produced by specific particle
std::map<art::Ptr<SimParticle>,unsigned> pmap;
for(auto const& mcdigi : *mcdigis) {
// look at the early end
StrawEnd fend = mcdigi.earlyEnd();
auto const& step = mcdigi.strawGasStep(fend);
art::Ptr<SimParticle> const& sp = step->simParticle();
auto const& mom = step->momentum(); // cast to 3-vector
if(debug_ > 0)std::cout <<"SimParticle PDG = " << sp->pdgId() << " Mom = " << sqrt(mom.mag2()) << std::endl;
bool goodpdg(true);
if(pdgs_.size() > 0)
goodpdg = std::find(pdgs_.begin(),pdgs_.end(),sp->pdgId()) != pdgs_.end();
if(goodpdg && sqrt(mom.mag2()) > minpmom_ && sqrt(mom.mag2()) < maxpmom_ ){
auto mapfnd = pmap.find(sp);
if(mapfnd == pmap.end())
pmap[sp] = 1;
else
++mapfnd->second;
// do not inspect truth info for digis not produced in simulation
if (mcdigi.containsSimulation()){
// look at the early end
StrawEnd fend = mcdigi.earlyEnd();
auto const& step = mcdigi.strawGasStep(fend);
art::Ptr<SimParticle> const& sp = step->simParticle();
auto const& mom = step->momentum(); // cast to 3-vector
if(debug_ > 0)std::cout <<"SimParticle PDG = " << sp->pdgId() << " Mom = " << sqrt(mom.mag2()) << std::endl;
bool goodpdg(true);
if(pdgs_.size() > 0)
goodpdg = std::find(pdgs_.begin(),pdgs_.end(),sp->pdgId()) != pdgs_.end();
if(goodpdg && sqrt(mom.mag2()) > minpmom_ && sqrt(mom.mag2()) < maxpmom_ ){
auto mapfnd = pmap.find(sp);
if(mapfnd == pmap.end())
pmap[sp] = 1;
else
++mapfnd->second;
}
}
}
// check if any single particle generated enough digis. Save All the particles
Expand Down
28 changes: 28 additions & 0 deletions MCDataProducts/inc/DigiProvenance.hh
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,28 @@
// Flag to designate a digi as having been produced from simulation,
// read from preexisting data, or a hybridization fo the two
// Ed Callaghan, 2024

#ifndef MCDataProducts_DigiProvenance_hh
#define MCDataProducts_DigiProvenance_hh

// stl
#include <map>
#include <string>

// mu2e
#include "Offline/GeneralUtilities/inc/EnumToStringSparse.hh"

namespace mu2e{
// enum equipped with std::string descriptions
class DigiProvenanceDetail{
public:
enum enum_type {unknown=0, Simulation, Mixed, External};
static std::string const& typeName();
static std::map<enum_type, std::string> const& names();
};
using DigiProvenance = EnumToStringSparse<DigiProvenanceDetail>;

bool containsSimulation(const DigiProvenance&);
}

#endif
7 changes: 7 additions & 0 deletions MCDataProducts/inc/KalSeedMC.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,12 +10,14 @@
#include "Offline/DataProducts/inc/GenVector.hh"
#include "Offline/DataProducts/inc/StrawEnd.hh"
#include "Offline/DataProducts/inc/VirtualDetectorId.hh"
#include <Offline/GeneralUtilities/inc/EnumToStringSparse.hh>
#include "Offline/RecoDataProducts/inc/KalSeed.hh"
#include "Offline/MCDataProducts/inc/SimParticle.hh"
#include "Offline/MCDataProducts/inc/ProcessCode.hh"
#include "Offline/MCDataProducts/inc/GenId.hh"
#include "Offline/MCDataProducts/inc/MCRelationship.hh"
#include "Offline/MCDataProducts/inc/CaloClusterMC.hh"
#include "Offline/MCDataProducts/inc/DigiProvenance.hh"
#include "art/Framework/Principal/Handle.h"
#include "cetlib/map_vector.h"
#include <Rtypes.h>
Expand Down Expand Up @@ -70,6 +72,8 @@ namespace mu2e {
//
// MC information for TrackStrawHits on this fit
struct TrkStrawHitMC {
TrkStrawHitMC(): _provenance(DigiProvenance::Simulation) {}
bool containsSimulation() const;
StrawHitIndex strawDigiMCIndex() const { return _sdmcindex; }
StrawHitIndex simPartStubIndex() const { return _spindex; }
StrawId const& strawid() const { return _strawId; }
Expand All @@ -95,15 +99,18 @@ namespace mu2e {
float _wireTau; // threshold cluster distance to the wire along the perpedicular particle path
float _strawDOCA; // signed doca to straw
float _strawPhi; // cylindrical phi from -pi to pi with 0 in Z direction
DigiProvenance _provenance; // origin/validity of MC info object
};

struct KalSeedMC {
bool containsSimulation() const;
SimPartStub const& simParticle(size_t index=0) const { return _simps.at(index); }
std::vector<SimPartStub> const& simParticles() const { return _simps; }
std::vector<TrkStrawHitMC> const& trkStrawHitMCs() const { return _tshmcs; }
std::vector<VDStep> const& vdSteps() const { return _vdsteps; }
TrkStrawHitMC const& trkStrawHitMC(size_t index) const { return _tshmcs.at(index); }
SimPartStub const& simParticle(TrkStrawHitMC const& tshmc) const { return simParticle(tshmc.simPartStubIndex()); }
bool ContainsSimulation() const;
// data products
std::vector<SimPartStub> _simps; // associated sim particles, and their relationship
std::vector<TrkStrawHitMC> _tshmcs; // MC info for each TrkStrawHitSeed
Expand Down
Loading

0 comments on commit 34a99ea

Please sign in to comment.