LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
MuonFilter_module.cc
Go to the documentation of this file.
1 //
3 // MuonFilter
4 //
5 // This event filter can act to identify events with only through-going tracks
6 //
7 // pagebri3@msu.edu
8 //
10 
11 // LArSoft includes
19 
20 // Framework includes
29 #include "fhiclcpp/ParameterSet.h"
31 
32 // ROOT includes
33 #include "TMathBase.h"
34 
35 // C++ includes
36 #include <map>
37 #include <memory>
38 #include <string>
39 
40 namespace filter {
41 
42  class MuonFilter : public art::EDFilter {
43  public:
44  explicit MuonFilter(fhicl::ParameterSet const&);
45 
46  private:
47  bool filter(art::Event& evt) override;
48 
49  std::string const fClusterModuleLabel;
50  std::string const fLineModuleLabel;
51  std::vector<double> const fCuts;
52  double const fDCenter;
53  double const fDelay;
54  double const fTolerance;
55  double const fMaxIon;
56  double const fIonFactor;
57  int const fDeltaWire;
58 
59  }; // class MuonFilter
60 
61  //-------------------------------------------------
63  : EDFilter{pset}
64  , fClusterModuleLabel{pset.get<std::string>("ClusterModuleLabel")}
65  , fLineModuleLabel{pset.get<std::string>("LineModuleLabel")}
66  , fCuts{pset.get<std::vector<double>>("Cuts")}
67  , fDCenter{pset.get<double>("DCenter")}
68  , fDelay{pset.get<double>("Delay")}
69  , fTolerance{pset.get<double>("Tolerance")}
70  , fMaxIon{pset.get<double>("MaxIon")}
71  , fIonFactor{pset.get<double>("IonFactor")}
72  , fDeltaWire{pset.get<int>("DeltaWire")}
73  {}
74 
75  //-------------------------------------------------
77  {
78  auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout const>()->Get();
79  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
80  auto const detProp =
82 
83  // Drift Velocity in cm/us Sampling rate in ns
84  double drift = detProp.DriftVelocity(detProp.Efield(), detProp.Temperature()) *
85  sampling_rate(clockData) / 1000.0;
86 
87  // This code only works comparing 2 planes so for now these are the
88  // last induction plane and collection plane
89  geo::TPCID const tpcid{0, 0};
90  unsigned int vPlane = wireReadoutGeom.Nplanes(tpcid) - 1;
91  geo::View_t vView = wireReadoutGeom.Plane({tpcid, vPlane}).View();
92  unsigned int uPlane = vPlane - 1;
93  geo::View_t uView = wireReadoutGeom.Plane({tpcid, uPlane}).View();
95  evt.getByLabel(fClusterModuleLabel, clustHandle);
96 
97  art::FindManyP<recob::Hit> fmh(clustHandle, evt, fClusterModuleLabel);
98 
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;
127  art::PtrVector<recob::Cluster> inductionSegments, collectionSegments;
128  evt.getByLabel(fLineModuleLabel, lines);
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  std::vector<art::Ptr<recob::Hit>> colHits = fmhc.at(j);
161 
162  double trk1Start = indSeg->StartTick() + fDelay;
163  double trk1End = indSeg->EndTick() + fDelay;
164  double trk2Start = colSeg->StartTick();
165  double trk2End = colSeg->EndTick();
166 
167  int uPos1 = indSeg->StartWire();
168  int uPos2 = indSeg->EndWire();
169  int const vPos1 = colSeg->StartWire();
170  int const vPos2 = colSeg->EndWire();
171  mf::LogInfo("MuonFilter") << "I J " << i << " " << j;
172  mf::LogInfo("MuonFilter")
173  << "Start/end " << indSeg->StartWire() << " " << colSeg->StartWire() << " "
174  << indSeg->EndWire() << " " << colSeg->EndWire();
175  mf::LogInfo("MuonFilter")
176  << "U's " << uPos1 << " " << uPos2 << "V's " << vPos1 << " " << vPos2 << " times "
177  << trk1End << " " << trk2End << " " << trk1Start << " " << trk2Start;
178  // need to have the corresponding endpoints matched
179  // check if they match in this order else switch
180  // really should use the crossing function and then have limits
181  // on distance outide tpc, or some other way of dealing with
182  // imperfect matches
183  if ((TMath::Abs(uPos1 - vPos1) > fDeltaWire || TMath::Abs(uPos2 - vPos2) > fDeltaWire) &&
184  (TMath::Abs(uPos1 - vPos2) <= fDeltaWire &&
185  TMath::Abs(uPos2 - vPos1) <= fDeltaWire)) {
186  mf::LogInfo("MuonFilter") << "Swapped1";
187  std::swap(uPos1, uPos2);
188  }
189  // check for time tolerance
190  if ((TMath::Abs(trk1Start - trk2Start) > fTolerance &&
191  TMath::Abs(trk1End - trk2End) > fTolerance) &&
192  (TMath::Abs(trk1Start - trk2End) < fTolerance &&
193  TMath::Abs(trk1End - trk2Start) < fTolerance)) {
194  std::swap(trk1Start, trk1End);
195  std::swap(uPos1, uPos2);
196  mf::LogInfo("MuonFilter") << "Swapped2";
197  }
198  mf::LogInfo("MuonFilter")
199  << "Times: " << trk1Start << " " << trk2Start << " " << trk1End << " " << trk2End;
200  // again needs to be fixed
204  if ((TMath::Abs(trk1Start - trk2Start) < fTolerance &&
205  TMath::Abs(trk1End - trk2End) < fTolerance) &&
206  (TMath::Abs(uPos1 - vPos1) <= fDeltaWire + 2 &&
207  TMath::Abs(uPos2 - vPos2) <= fDeltaWire + 2)) {
208  geo::WireID u_wID1(indSeg->Plane(), uPos1);
209  geo::WireID u_wID2(indSeg->Plane(), uPos2);
210  geo::WireID v_wID1(colSeg->Plane(), vPos1);
211  geo::WireID v_wID2(colSeg->Plane(), vPos2);
212 
213  auto const intersection1 = wireReadoutGeom.WireIDsIntersect(u_wID1, v_wID1);
214  auto const intersection2 = wireReadoutGeom.WireIDsIntersect(u_wID2, v_wID2);
215  if (!intersection1 || !intersection2) { continue; }
216  auto const [y1, z1] = std::make_pair(intersection1->y, intersection1->z);
217  auto const [y2, z2] = std::make_pair(intersection2->y, intersection2->z);
218 
219  double const x1 = (trk1Start + trk2Start) / 2.0 * drift - fDCenter;
220  double const x2 = (trk1End + trk2End) / 2.0 * drift - fDCenter;
221  mf::LogInfo("MuonFilter") << "Match " << matchNum << " " << x1 << " " << y1 << " " << z1
222  << " " << x2 << " " << y2 << " " << z2;
223  bool x1edge, x2edge, y1edge, y2edge, z1edge, z2edge;
224  indMap[i] = matchNum;
225  colMap[j] = matchNum;
226  matchNum++;
227  pointTemp[0] = x1;
228  pointTemp[1] = y1;
229  pointTemp[2] = z1;
230  pointTemp[3] = x2;
231  pointTemp[4] = y2;
232  pointTemp[5] = z2;
233  x1edge = (TMath::Abs(x1) - fCuts[0] > 0);
234  x2edge = (TMath::Abs(x2) - fCuts[0] > 0);
235  y1edge = (TMath::Abs(y1) - fCuts[1] > 0);
236  y2edge = (TMath::Abs(y2) - fCuts[1] > 0);
237  z1edge = (TMath::Abs(z1) - fCuts[2] > 0);
238  z2edge = (TMath::Abs(z2) - fCuts[2] > 0);
239  if ((x1edge || y1edge || z1edge) && (x2edge || y2edge || z2edge)) {
240  tGoing.push_back(pointTemp);
241  mf::LogInfo("MuonFilter") << "outside Removed induction ion: ";
242 
243  for (size_t h = 0; h < indHits.size(); h++) {
244  mf::LogInfo("MuonFilter") << indHits[h]->PeakAmplitude() << " ";
245  indIon -= indHits[h]->PeakAmplitude();
246  }
247  mf::LogInfo("MuonFilter") << "Removed collection ion: ";
248 
249  for (size_t h = 0; h < colHits.size(); h++) {
250  mf::LogInfo("MuonFilter") << colHits[h]->PeakAmplitude() << " ";
251  colIon -= colHits[h]->PeakAmplitude();
252  }
253  mf::LogInfo("MuonFilter")
254  << "Ionization outside track I/C: " << indIon << " " << colIon;
255  }
256  else if ((x1edge || y1edge || z1edge) && !(x2edge || y2edge || z2edge) &&
257  (z2 - z1) > 1.2) {
258  tGoing.push_back(pointTemp);
259  mf::LogInfo("MuonFilter") << "stopping Removed induction ion: ";
260  for (size_t h = 0; h < indHits.size(); h++) {
261  mf::LogInfo("MuonFilter") << indHits[h]->PeakAmplitude() << " ";
262  indIon -= indHits[h]->PeakAmplitude();
263  }
264  mf::LogInfo("MuonFilter") << "Removed collection ion: ";
265  for (size_t h = 0; h < colHits.size(); h++) {
266  mf::LogInfo("MuonFilter") << colHits[h]->PeakAmplitude() << " ";
267  colIon -= colHits[h]->PeakAmplitude();
268  }
269  mf::LogInfo("MuonFilter")
270  << "Ionization outside track I/C: " << indIon << " " << colIon;
271  }
272  else {
273  pairTemp = std::make_pair(i, j);
274  mf::LogInfo("MuonFilter") << "rLook matchnum " << matchNum << " " << i << " " << j;
275  rLook.push_back(pairTemp);
276  matched.push_back(pointTemp);
277  }
278  break; // advances i, makes j=0;
279  }
280  }
281  }
282  }
283  // after all matches are made, remove deltas
284  double distance = 0;
285  for (unsigned int i = 0; i < tGoing.size(); i++)
286  for (unsigned int j = 0; j < matched.size(); j++) {
287  mf::LogInfo("MuonFilter") << tGoing.size() << " " << matched.size() << " " << i << " " << j;
288  // test if one is contained within the other in the z-direction
289  if ((tGoing[i][2] <= matched[j][2]) && (tGoing[i][5] >= matched[j][5])) {
290  TVector3 a1(&tGoing[i][0]);
291  TVector3 a2(&tGoing[i][3]);
292  TVector3 b1(&matched[j][0]);
293  distance = TMath::Abs((((a1 - a2).Cross((a1 - a2).Cross(a1 - b1))).Unit()).Dot(a1 - b1));
294  mf::LogInfo("MuonFilter") << "distance " << distance;
295  if (distance < 6) {
296  mf::LogInfo("MuonFilter")
297  << "Removing delta ion " << rLook.size() << " " << rLook[j].first << " " << matchNum;
298  std::vector<art::Ptr<recob::Hit>> temp = fmhi.at(rLook[j].first);
299  for (unsigned int h = 0; h < temp.size(); h++)
300  indIon -= temp[h]->PeakAmplitude();
301  temp = fmhc.at(rLook[j].second);
302  for (unsigned int h = 0; h < temp.size(); h++)
303  colIon -= temp[h]->PeakAmplitude();
304  }
305  }
306  }
307  mf::LogInfo("MuonFilter") << "indIon " << indIon * fIonFactor << " colIon " << colIon;
308  if ((indIon * fIonFactor > fMaxIon) && (colIon > fMaxIon))
309  return true;
310  else
311  return false;
312  }
313 
315 
316 } // namespace filt
void reserve(size_type n)
Definition: PtrVector.h:337
bool filter(art::Event &evt) override
Float_t y1[n_points_granero]
Definition: compare.C:5
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
Float_t x1[n_points_granero]
Definition: compare.C:5
Declaration of signal hit object.
float StartWire() const
Returns the wire coordinate of the start of the cluster.
Definition: Cluster.h:276
std::vector< double > const fCuts
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
float EndTick() const
Returns the tick coordinate of the end of the cluster.
Definition: Cluster.h:331
geo::PlaneID Plane() const
Returns the plane ID this cluster lies on.
Definition: Cluster.h:717
Cluster finding and building.
Framework includes.
Float_t y2[n_points_geant4]
Definition: compare.C:26
int const fDeltaWire
allowed differences in wire number between 2 planes
#define a2
std::string const fLineModuleLabel
void hits()
Definition: readHits.C:15
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
MuonFilter(fhicl::ParameterSet const &)
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:435
std::string const fClusterModuleLabel
The data type to uniquely identify a TPC.
Definition: geo_types.h:306
Declaration of cluster object.
size_type size() const
Definition: PtrVector.h:302
Detector simulation of raw signals on wires.
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Encapsulate the construction of a single detector plane .
EDFilter(fhicl::ParameterSet const &pset)
Definition: EDFilter.cc:6
TCEvent evt
Definition: DataStructs.cxx:8
Float_t x2[n_points_geant4]
Definition: compare.C:26
second_as<> second
Type of time stored in seconds, in double precision.
Definition: spacetime.h:82
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
float StartTick() const
Returns the tick coordinate of the start of the cluster.
Definition: Cluster.h:287
art framework interface to geometry description
#define a1
float EndWire() const
Returns the wire coordinate of the end of the cluster.
Definition: Cluster.h:318