33 #include "TMathBase.h" 51 std::vector<double>
const fCuts;
66 ,
fCuts{pset.get<std::vector<double>>(
"Cuts")}
67 ,
fDCenter{pset.get<
double>(
"DCenter")}
68 ,
fDelay{pset.get<
double>(
"Delay")}
70 ,
fMaxIon{pset.get<
double>(
"MaxIon")}
84 double drift = detProp.DriftVelocity(detProp.Efield(), detProp.Temperature()) *
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();
100 clusters.
reserve(clustHandle->size());
101 for (
unsigned int i = 0; i < clustHandle->size(); ++i) {
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;
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;
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();
120 if (clusters[
cluster]->View() == uView)
122 else if (clusters[
cluster]->View() == vView)
125 mf::LogInfo(
"MuonFilter") <<
"Ionizations: " << indIon <<
" " << colIon;
130 lineVec.
reserve(lines->size());
131 for (
unsigned int i = 0; i < lines->size(); ++i) {
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]);
146 if (inductionSegments.
size() == 0 || collectionSegments.
size() == 0) {
147 mf::LogInfo(
"MuonFilter") <<
"At least one plane with no track";
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;
159 std::vector<art::Ptr<recob::Hit>> indHits = fmhi.at(i);
160 std::vector<art::Ptr<recob::Hit>> colHits = fmhc.at(j);
165 double trk2End = colSeg->
EndTick();
170 int const vPos2 = colSeg->
EndWire();
171 mf::LogInfo(
"MuonFilter") <<
"I J " << i <<
" " << j;
176 <<
"U's " << uPos1 <<
" " << uPos2 <<
"V's " << vPos1 <<
" " << vPos2 <<
" times " 177 << trk1End <<
" " << trk2End <<
" " << trk1Start <<
" " << trk2Start;
187 std::swap(uPos1, uPos2);
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);
199 <<
"Times: " << trk1Start <<
" " << trk2Start <<
" " << trk1End <<
" " << trk2End;
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)) {
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);
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;
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: ";
243 for (
size_t h = 0; h < indHits.size(); h++) {
244 mf::LogInfo(
"MuonFilter") << indHits[h]->PeakAmplitude() <<
" ";
245 indIon -= indHits[h]->PeakAmplitude();
247 mf::LogInfo(
"MuonFilter") <<
"Removed collection ion: ";
249 for (
size_t h = 0; h < colHits.size(); h++) {
250 mf::LogInfo(
"MuonFilter") << colHits[h]->PeakAmplitude() <<
" ";
251 colIon -= colHits[h]->PeakAmplitude();
254 <<
"Ionization outside track I/C: " << indIon <<
" " << colIon;
256 else if ((x1edge || y1edge || z1edge) && !(x2edge || y2edge || z2edge) &&
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();
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();
270 <<
"Ionization outside track I/C: " << indIon <<
" " << colIon;
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);
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;
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;
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();
void reserve(size_type n)
bool filter(art::Event &evt) override
Float_t y1[n_points_granero]
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
Float_t x1[n_points_granero]
Declaration of signal hit object.
float StartWire() const
Returns the wire coordinate of the start of the cluster.
std::vector< double > const fCuts
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
float EndTick() const
Returns the tick coordinate of the end of the cluster.
geo::PlaneID Plane() const
Returns the plane ID this cluster lies on.
Cluster finding and building.
Float_t y2[n_points_geant4]
int const fDeltaWire
allowed differences in wire number between 2 planes
std::string const fLineModuleLabel
#define DEFINE_ART_MODULE(klass)
MuonFilter(fhicl::ParameterSet const &)
void push_back(Ptr< U > const &p)
std::string const fClusterModuleLabel
The data type to uniquely identify a TPC.
Declaration of cluster object.
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)
Float_t x2[n_points_geant4]
second_as<> second
Type of time stored in seconds, in double precision.
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.
art framework interface to geometry description
float EndWire() const
Returns the wire coordinate of the end of the cluster.