Public Member Functions | Private Member Functions | Private Attributes | List of all members
filter::MuonFilter Class Reference
Inheritance diagram for filter::MuonFilter:
art::EDFilter art::detail::Filter art::detail::LegacyModule art::Modifier art::ModuleBase art::ProductRegistryHelper

Public Member Functions

 MuonFilter (fhicl::ParameterSet const &)
 
- Public Member Functions inherited from art::EDFilter
 EDFilter (fhicl::ParameterSet const &pset)
 
template<typename Config >
 EDFilter (Table< Config > const &config)
 
std::string workerType () const
 
- Public Member Functions inherited from art::detail::Filter
virtual ~Filter () noexcept
 
 Filter (fhicl::ParameterSet const &)
 
 Filter (Filter const &)=delete
 
 Filter (Filter &&)=delete
 
Filter & operator= (Filter const &)=delete
 
Filter & operator= (Filter &&)=delete
 
void doBeginJob (SharedResources const &resources)
 
void doEndJob ()
 
void doRespondToOpenInputFile (FileBlock const &fb)
 
void doRespondToCloseInputFile (FileBlock const &fb)
 
void doRespondToOpenOutputFiles (FileBlock const &fb)
 
void doRespondToCloseOutputFiles (FileBlock const &fb)
 
bool doBeginRun (RunPrincipal &rp, ModuleContext const &)
 
bool doEndRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)
 
- Public Member Functions inherited from art::Modifier
 ~Modifier () noexcept
 
 Modifier ()
 
 Modifier (Modifier const &)=delete
 
 Modifier (Modifier &&)=delete
 
Modifier & operator= (Modifier const &)=delete
 
Modifier & operator= (Modifier &&)=delete
 
- Public Member Functions inherited from art::ModuleBase
virtual ~ModuleBase () noexcept
 
 ModuleBase ()
 
ModuleDescription const & moduleDescription () const
 
void setModuleDescription (ModuleDescription const &)
 
std::array< std::vector< ProductInfo >, NumBranchTypes > const & getConsumables () const
 
void sortConsumables (std::string const &current_process_name)
 
template<typename T , BranchType BT>
ViewToken< T > consumesView (InputTag const &tag)
 
template<typename T , BranchType BT>
ViewToken< T > mayConsumeView (InputTag const &tag)
 

Private Member Functions

bool filter (art::Event &evt) override
 

Private Attributes

std::string const fClusterModuleLabel
 
std::string const fLineModuleLabel
 
std::vector< double > const fCuts
 
double const fDCenter
 
double const fDelay
 
double const fTolerance
 
double const fMaxIon
 
double const fIonFactor
 
int const fDeltaWire
 allowed differences in wire number between 2 planes More...
 

Additional Inherited Members

- Public Types inherited from art::EDFilter
using ModuleType = EDFilter
 
using WorkerType = WorkerT< EDFilter >
 
- Public Types inherited from art::detail::Filter
template<typename UserConfig >
using Table = Modifier::Table< UserConfig >
 
- Public Types inherited from art::Modifier
template<typename UserConfig , typename UserKeysToIgnore = void>
using Table = ProducerTable< UserConfig, detail::ModuleConfig, UserKeysToIgnore >
 
- Static Public Attributes inherited from art::detail::Filter
static constexpr bool Pass {true}
 
static constexpr bool Fail {false}
 
- Protected Member Functions inherited from art::ModuleBase
ConsumesCollector & consumesCollector ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > consumes (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > consumesView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void consumesMany ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > mayConsume (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > mayConsumeView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void mayConsumeMany ()
 

Detailed Description

Definition at line 41 of file MuonFilter_module.cc.

Constructor & Destructor Documentation

filter::MuonFilter::MuonFilter ( fhicl::ParameterSet const &  pset)
explicit

Definition at line 61 of file MuonFilter_module.cc.

62  : EDFilter{pset}
63  , fClusterModuleLabel{pset.get<std::string>("ClusterModuleLabel")}
64  , fLineModuleLabel{pset.get<std::string>("LineModuleLabel")}
65  , fCuts{pset.get<std::vector<double>>("Cuts")}
66  , fDCenter{pset.get<double>("DCenter")}
67  , fDelay{pset.get<double>("Delay")}
68  , fTolerance{pset.get<double>("Tolerance")}
69  , fMaxIon{pset.get<double>("MaxIon")}
70  , fIonFactor{pset.get<double>("IonFactor")}
71  , fDeltaWire{pset.get<int>("DeltaWire")}
72  {}
string
std::string string
Definition: nybbler.cc:12
filter::MuonFilter::fCuts
std::vector< double > const fCuts
Definition: MuonFilter_module.cc:50
filter::MuonFilter::fDeltaWire
int const fDeltaWire
allowed differences in wire number between 2 planes
Definition: MuonFilter_module.cc:56
filter::MuonFilter::fLineModuleLabel
std::string const fLineModuleLabel
Definition: MuonFilter_module.cc:49
filter::MuonFilter::fClusterModuleLabel
std::string const fClusterModuleLabel
Definition: MuonFilter_module.cc:48
art::EDFilter::EDFilter
EDFilter(fhicl::ParameterSet const &pset)
Definition: EDFilter.h:21

Member Function Documentation

bool filter::MuonFilter::filter ( art::Event &  evt)
overrideprivatevirtual
Todo:
: the delta wire numbers seem a bit magic,
Todo:
: should also change to using Geometry::ChannelsIntersect method

Implements art::EDFilter.

Definition at line 76 of file MuonFilter_module.cc.

77  {
78  art::ServiceHandle<geo::Geometry const> geom;
79  auto const clockData =
80  art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
81  auto const detProp =
82  art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt, clockData);
83 
84  // Drift Velocity in cm/us Sampling rate in ns
85  double drift = detProp.DriftVelocity(detProp.Efield(), detProp.Temperature()) *
86  sampling_rate(clockData) / 1000.0;
87 
88  // This code only works comparing 2 planes so for now these are the
89  // last induction plane and collection plane
90  int vPlane = geom->Nplanes() - 1;
91  geo::View_t vView = geom->Plane(vPlane).View();
92  int uPlane = vPlane - 1;
93  geo::View_t uView = geom->Plane(uPlane).View();
94  art::Handle<std::vector<recob::Cluster>> clustHandle;
95  evt.getByLabel(fClusterModuleLabel, clustHandle);
96 
97  art::FindManyP<recob::Hit> fmh(clustHandle, evt, fClusterModuleLabel);
98 
99  art::PtrVector<recob::Cluster> clusters;
100  clusters.reserve(clustHandle->size());
101  for (unsigned int i = 0; i < clustHandle->size(); ++i) {
102  clusters.push_back(art::Ptr<recob::Cluster>(clustHandle, i));
103  }
104  double indIon(0), colIon(0);
105  std::map<int, int> indMap;
106  std::map<int, int> colMap;
107  std::vector<std::pair<int, int>> rLook;
108  int matchNum = 1;
109  std::vector<std::vector<double>> tGoing;
110  std::vector<std::vector<double>> matched;
111  std::vector<double> pointTemp(6);
112  std::pair<int, int> pairTemp;
113  double ionSum(0.0);
114  for (size_t cluster = 0; cluster < clusters.size(); ++cluster) {
115  ionSum = 0.0;
116  std::vector<art::Ptr<recob::Hit>> hits = fmh.at(cluster);
117  for (unsigned int hit = 0; hit < hits.size(); hit++) {
118  ionSum += hits[hit]->PeakAmplitude();
119  }
120  if (clusters[cluster]->View() == uView)
121  indIon += ionSum;
122  else if (clusters[cluster]->View() == vView)
123  colIon += ionSum;
124  }
125  mf::LogInfo("MuonFilter") << "Ionizations: " << indIon << " " << colIon;
126  art::Handle<std::vector<recob::Cluster>> lines;
127  art::PtrVector<recob::Cluster> inductionSegments, collectionSegments;
128  evt.getByLabel(fLineModuleLabel, lines);
129  art::PtrVector<recob::Cluster> lineVec;
130  lineVec.reserve(lines->size());
131  for (unsigned int i = 0; i < lines->size(); ++i) {
132  lineVec.push_back(art::Ptr<recob::Cluster>(lines, i));
133  }
134 
135  for (size_t cl = 0; cl < clusters.size(); cl++) {
136  std::vector<art::Ptr<recob::Hit>> hits = fmh.at(cl);
137  if (hits.size() > 0 && clusters[cl]->View() == uView)
138  inductionSegments.push_back(clusters[cl]);
139  else if (hits.size() > 0 && clusters[cl]->View() == vView)
140  collectionSegments.push_back(clusters[cl]);
141  }
142 
143  art::FindManyP<recob::Hit> fmhi(inductionSegments, evt, fClusterModuleLabel);
144  art::FindManyP<recob::Hit> fmhc(collectionSegments, evt, fClusterModuleLabel);
145 
146  if (inductionSegments.size() == 0 || collectionSegments.size() == 0) {
147  mf::LogInfo("MuonFilter") << "At least one plane with no track";
148  }
149  else {
150 
151  for (unsigned int i = 0; i < inductionSegments.size(); i++) {
152  if (indMap[i]) continue;
153  for (unsigned int j = 0; j < collectionSegments.size(); j++) {
154  if (colMap[j]) continue;
155 
156  art::Ptr<recob::Cluster> indSeg = inductionSegments[i];
157  art::Ptr<recob::Cluster> colSeg = collectionSegments[j];
158 
159  std::vector<art::Ptr<recob::Hit>> indHits = fmhi.at(i);
160 
161  std::vector<art::Ptr<recob::Hit>> colHits = fmhc.at(j);
162 
163  double trk1Start = indSeg->StartTick() + fDelay;
164  double trk1End = indSeg->EndTick() + fDelay;
165  double trk2Start = colSeg->StartTick();
166  double trk2End = colSeg->EndTick();
167 
168  int uPos1 = indSeg->StartWire();
169  int uPos2 = indSeg->EndWire();
170  int const vPos1 = colSeg->StartWire();
171  int const vPos2 = colSeg->EndWire();
172  mf::LogInfo("MuonFilter") << "I J " << i << " " << j;
173  mf::LogInfo("MuonFilter")
174  << "Start/end " << indSeg->StartWire() << " " << colSeg->StartWire() << " "
175  << indSeg->EndWire() << " " << colSeg->EndWire();
176  mf::LogInfo("MuonFilter")
177  << "U's " << uPos1 << " " << uPos2 << "V's " << vPos1 << " " << vPos2 << " times "
178  << trk1End << " " << trk2End << " " << trk1Start << " " << trk2Start;
179  // need to have the corresponding endpoints matched
180  // check if they match in this order else switch
181  // really should use the crossing function and then have limits
182  // on distance outide tpc, or some other way of dealing with
183  // imperfect matches
184  if ((TMath::Abs(uPos1 - vPos1) > fDeltaWire || TMath::Abs(uPos2 - vPos2) > fDeltaWire) &&
185  (TMath::Abs(uPos1 - vPos2) <= fDeltaWire &&
186  TMath::Abs(uPos2 - vPos1) <= fDeltaWire)) {
187  mf::LogInfo("MuonFilter") << "Swapped1";
188  std::swap(uPos1, uPos2);
189  }
190  // check for time tolerance
191  if ((TMath::Abs(trk1Start - trk2Start) > fTolerance &&
192  TMath::Abs(trk1End - trk2End) > fTolerance) &&
193  (TMath::Abs(trk1Start - trk2End) < fTolerance &&
194  TMath::Abs(trk1End - trk2Start) < fTolerance)) {
195  std::swap(trk1Start, trk1End);
196  std::swap(uPos1, uPos2);
197  mf::LogInfo("MuonFilter") << "Swapped2";
198  }
199  mf::LogInfo("MuonFilter")
200  << "Times: " << trk1Start << " " << trk2Start << " " << trk1End << " " << trk2End;
201  // again needs to be fixed
202  ///\todo: the delta wire numbers seem a bit magic,
203  ///\todo: should also change to using Geometry::ChannelsIntersect
204  ///method
205  if ((TMath::Abs(trk1Start - trk2Start) < fTolerance &&
206  TMath::Abs(trk1End - trk2End) < fTolerance) &&
207  (TMath::Abs(uPos1 - vPos1) <= fDeltaWire + 2 &&
208  TMath::Abs(uPos2 - vPos2) <= fDeltaWire + 2)) {
209  geo::WireID u_wID1(indSeg->Plane(), uPos1);
210  geo::WireID u_wID2(indSeg->Plane(), uPos2);
211  geo::WireID v_wID1(colSeg->Plane(), vPos1);
212  geo::WireID v_wID2(colSeg->Plane(), vPos2);
213 
214  double y1, y2, z1, z2;
215  geom->IntersectionPoint(u_wID1, v_wID1, y1, z1);
216  geom->IntersectionPoint(u_wID2, v_wID2, y2, z2);
217 
218  double const x1 = (trk1Start + trk2Start) / 2.0 * drift - fDCenter;
219  double const x2 = (trk1End + trk2End) / 2.0 * drift - fDCenter;
220  mf::LogInfo("MuonFilter") << "Match " << matchNum << " " << x1 << " " << y1 << " " << z1
221  << " " << x2 << " " << y2 << " " << z2;
222  bool x1edge, x2edge, y1edge, y2edge, z1edge, z2edge;
223  indMap[i] = matchNum;
224  colMap[j] = matchNum;
225  matchNum++;
226  pointTemp[0] = x1;
227  pointTemp[1] = y1;
228  pointTemp[2] = z1;
229  pointTemp[3] = x2;
230  pointTemp[4] = y2;
231  pointTemp[5] = z2;
232  x1edge = (TMath::Abs(x1) - fCuts[0] > 0);
233  x2edge = (TMath::Abs(x2) - fCuts[0] > 0);
234  y1edge = (TMath::Abs(y1) - fCuts[1] > 0);
235  y2edge = (TMath::Abs(y2) - fCuts[1] > 0);
236  z1edge = (TMath::Abs(z1) - fCuts[2] > 0);
237  z2edge = (TMath::Abs(z2) - fCuts[2] > 0);
238  if ((x1edge || y1edge || z1edge) && (x2edge || y2edge || z2edge)) {
239  tGoing.push_back(pointTemp);
240  mf::LogInfo("MuonFilter") << "outside Removed induction ion: ";
241 
242  for (size_t h = 0; h < indHits.size(); h++) {
243  mf::LogInfo("MuonFilter") << indHits[h]->PeakAmplitude() << " ";
244  indIon -= indHits[h]->PeakAmplitude();
245  }
246  mf::LogInfo("MuonFilter") << "Removed collection ion: ";
247 
248  for (size_t h = 0; h < colHits.size(); h++) {
249  mf::LogInfo("MuonFilter") << colHits[h]->PeakAmplitude() << " ";
250  colIon -= colHits[h]->PeakAmplitude();
251  }
252  mf::LogInfo("MuonFilter")
253  << "Ionization outside track I/C: " << indIon << " " << colIon;
254  }
255  else if ((x1edge || y1edge || z1edge) && !(x2edge || y2edge || z2edge) &&
256  (z2 - z1) > 1.2) {
257  tGoing.push_back(pointTemp);
258  mf::LogInfo("MuonFilter") << "stopping Removed induction ion: ";
259  for (size_t h = 0; h < indHits.size(); h++) {
260  mf::LogInfo("MuonFilter") << indHits[h]->PeakAmplitude() << " ";
261  indIon -= indHits[h]->PeakAmplitude();
262  }
263  mf::LogInfo("MuonFilter") << "Removed collection ion: ";
264  for (size_t h = 0; h < colHits.size(); h++) {
265  mf::LogInfo("MuonFilter") << colHits[h]->PeakAmplitude() << " ";
266  colIon -= colHits[h]->PeakAmplitude();
267  }
268  mf::LogInfo("MuonFilter")
269  << "Ionization outside track I/C: " << indIon << " " << colIon;
270  }
271  else {
272  pairTemp = std::make_pair(i, j);
273  mf::LogInfo("MuonFilter") << "rLook matchnum " << matchNum << " " << i << " " << j;
274  rLook.push_back(pairTemp);
275  matched.push_back(pointTemp);
276  }
277  break; // advances i, makes j=0;
278  }
279  }
280  }
281  }
282  // after all matches are made, remove deltas
283  double distance = 0;
284  for (unsigned int i = 0; i < tGoing.size(); i++)
285  for (unsigned int j = 0; j < matched.size(); j++) {
286  mf::LogInfo("MuonFilter") << tGoing.size() << " " << matched.size() << " " << i << " " << j;
287  // test if one is contained within the other in the z-direction
288  if ((tGoing[i][2] <= matched[j][2]) && (tGoing[i][5] >= matched[j][5])) {
289  TVector3 a1(&tGoing[i][0]);
290  TVector3 a2(&tGoing[i][3]);
291  TVector3 b1(&matched[j][0]);
292  distance = TMath::Abs((((a1 - a2).Cross((a1 - a2).Cross(a1 - b1))).Unit()).Dot(a1 - b1));
293  mf::LogInfo("MuonFilter") << "distance " << distance;
294  if (distance < 6) {
295  mf::LogInfo("MuonFilter")
296  << "Removing delta ion " << rLook.size() << " " << rLook[j].first << " " << matchNum;
297  std::vector<art::Ptr<recob::Hit>> temp = fmhi.at(rLook[j].first);
298  for (unsigned int h = 0; h < temp.size(); h++)
299  indIon -= temp[h]->PeakAmplitude();
300  temp = fmhc.at(rLook[j].second);
301  for (unsigned int h = 0; h < temp.size(); h++)
302  colIon -= temp[h]->PeakAmplitude();
303  }
304  }
305  }
306  mf::LogInfo("MuonFilter") << "indIon " << indIon * fIonFactor << " colIon " << colIon;
307  if ((indIon * fIonFactor > fMaxIon) && (colIon > fMaxIon))
308  return true;
309  else
310  return false;
311  }
art::PtrVector::reserve
void reserve(size_type n)
Definition: PtrVector.h:337
gar::rec::alg::cluster
void cluster(In first, In last, Out result, Pred *pred)
Definition: NNClusters.h:41
geo::GeometryCore::Plane
PlaneGeo const & Plane(unsigned int const p, unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified wire.
Definition: GeometryCore.h:2994
View
AdcChannelData::View View
Definition: test_AdcChannelData.cxx:23
geo::View_t
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
mf::LogInfo
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
Definition: MessageLogger.h:209
recob::Cluster::StartWire
float StartWire() const
Returns the wire coordinate of the start of the cluster.
Definition: Cluster.h:286
filter::MuonFilter::fCuts
std::vector< double > const fCuts
Definition: MuonFilter_module.cc:50
recob::Cluster::EndTick
float EndTick() const
Returns the tick coordinate of the end of the cluster.
Definition: Cluster.h:342
recob::Cluster::Plane
geo::PlaneID Plane() const
Returns the plane ID this cluster lies on.
Definition: Cluster.h:744
cluster
Cluster finding and building.
Definition: SmallClusterFilter_module.cc:34
reco_momentum_tuples.x2
x2
Definition: reco_momentum_tuples.py:83
cl
QAsciiDict< Entry > cl
Definition: tclscanner.cpp:1148
filter::MuonFilter::fDeltaWire
int const fDeltaWire
allowed differences in wire number between 2 planes
Definition: MuonFilter_module.cc:56
a2
#define a2
Definition: HoughBaseAlg.cxx:64
reco_momentum_tuples.h
h
Definition: reco_momentum_tuples.py:66
geo::PlaneGeo::View
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:184
filter::MuonFilter::fLineModuleLabel
std::string const fLineModuleLabel
Definition: MuonFilter_module.cc:49
art::DataViewImpl::getByLabel
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:633
geo::GeometryCore::Nplanes
unsigned int Nplanes(unsigned int tpc=0, unsigned int cstat=0) const
Returns the total number of wire planes in the specified TPC.
Definition: GeometryCore.h:2848
art::PtrVector::push_back
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:435
art::swap
void swap(Handle< T > &a, Handle< T > &b)
submit_ppfxjobs.temp
temp
Definition: submit_ppfxjobs.py:110
filter::MuonFilter::fClusterModuleLabel
std::string const fClusterModuleLabel
Definition: MuonFilter_module.cc:48
make_input_cov.lines
lines
Definition: make_input_cov.py:15
protoana::distance
double distance(double x1, double y1, double z1, double x2, double y2, double z2)
Definition: truepionXsection_module.cc:217
art::PtrVector::at
reference at(size_type n)
Definition: PtrVector.h:359
geo::GeometryCore::IntersectionPoint
bool IntersectionPoint(geo::WireID const &wid1, geo::WireID const &wid2, double &y, double &z) const
Returns the intersection point of two wires.
Definition: GeometryCore.cxx:1657
art::PtrVector::size
size_type size() const
Definition: PtrVector.h:302
hit
Detector simulation of raw signals on wires.
Definition: NumberOfHitsFilter_module.cc:25
util::quantities::second
second_as<> second
Type of time stored in seconds, in double precision.
Definition: spacetime.h:85
detinfo::sampling_rate
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
Definition: DetectorClocksData.h:643
recob::Cluster::StartTick
float StartTick() const
Returns the tick coordinate of the start of the cluster.
Definition: Cluster.h:297
reco_momentum_tuples.x1
x1
Definition: reco_momentum_tuples.py:82
a1
#define a1
Definition: HoughBaseAlg.cxx:63
recob::Cluster::EndWire
float EndWire() const
Returns the wire coordinate of the end of the cluster.
Definition: Cluster.h:329

Member Data Documentation

std::string const filter::MuonFilter::fClusterModuleLabel
private

Definition at line 48 of file MuonFilter_module.cc.

std::vector<double> const filter::MuonFilter::fCuts
private

Definition at line 50 of file MuonFilter_module.cc.

double const filter::MuonFilter::fDCenter
private

Definition at line 51 of file MuonFilter_module.cc.

double const filter::MuonFilter::fDelay
private

Definition at line 52 of file MuonFilter_module.cc.

int const filter::MuonFilter::fDeltaWire
private

allowed differences in wire number between 2 planes

Definition at line 56 of file MuonFilter_module.cc.

double const filter::MuonFilter::fIonFactor
private

Definition at line 55 of file MuonFilter_module.cc.

std::string const filter::MuonFilter::fLineModuleLabel
private

Definition at line 49 of file MuonFilter_module.cc.

double const filter::MuonFilter::fMaxIon
private

Definition at line 54 of file MuonFilter_module.cc.

double const filter::MuonFilter::fTolerance
private

Definition at line 53 of file MuonFilter_module.cc.

The documentation for this class was generated from the following file:
Generated by   doxygen 1.8.11