123 auto const * detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();
126 double drift = detprop->DriftVelocity(detprop->Efield(), detprop->Temperature())*detprop->SamplingRate()/1000.0;
130 int vPlane = geom->
Nplanes() - 1;
132 int uPlane = vPlane-1;
140 clusters.
reserve(clustHandle->size());
141 for(
unsigned int i = 0; i < clustHandle->size(); ++i) {
144 double indIon(0),colIon(0);
145 std::map<int,int> indMap;
146 std::map<int,int> colMap;
147 std::vector<std::pair<int,int> > rLook;
149 std::vector<std::vector<double> > tGoing;
150 std::vector<std::vector<double> > matched;
151 std::vector<double> pointTemp(6);
152 std::pair<int,int> pairTemp;
156 std::vector< art::Ptr<recob::Hit> >
hits = fmh.at(
cluster);
157 for(
unsigned int hit = 0;
hit < hits.size();
hit++) {
158 ionSum+=hits[
hit]->PeakAmplitude();
160 if(clusters[
cluster]->View() == uView) indIon+=ionSum;
161 else if(clusters[
cluster]->View() == vView) colIon+=ionSum;
163 mf::LogInfo(
"MuonFilter") <<
"Ionizations: " << indIon <<
" " << colIon ;
168 lineVec.
reserve(lines->size());
169 for(
unsigned int i = 0; i < lines->size(); ++i) {
173 for(
size_t cl = 0;cl < clusters.
size(); cl++) {
174 std::vector< art::Ptr<recob::Hit> > hits = fmh.at(cl);
175 if (hits.size()>0 && clusters[cl]->View()==uView)
176 inductionSegments.
push_back(clusters[cl]);
177 else if(hits.size()>0 && clusters[cl]->View() == vView) collectionSegments.
push_back(clusters[cl]);
183 if(inductionSegments.
size() == 0 || collectionSegments.
size() == 0) {
184 mf::LogInfo(
"MuonFilter") <<
"At least one plane with no track";
188 for(
unsigned int i = 0; i < inductionSegments.
size(); i++) {
189 if(indMap[i])
continue;
190 for(
unsigned int j = 0; j < collectionSegments.
size(); j++) {
191 if(colMap[j])
continue;
196 std::vector< art::Ptr<recob::Hit> > indHits = fmhi.
at(i);
198 std::vector< art::Ptr<recob::Hit> > colHits = fmhc.at(j);
203 double trk2End =colSeg->
EndTick();
208 int const vPos2 = colSeg->
EndWire();
209 mf::LogInfo(
"MuonFilter") <<
"I J " << i <<
" " << j ;
214 mf::LogInfo(
"MuonFilter")<<
"U's "<< uPos1 <<
" " << uPos2
215 <<
"V's "<< vPos1 <<
" " << vPos2
216 <<
" times " << trk1End <<
" "<< trk2End
217 <<
" "<< trk1Start <<
" "<< trk2Start ;
229 if((TMath::Abs(trk1Start-trk2Start) >
fTolerance && TMath::Abs(trk1End-trk2End) >
fTolerance) &&
235 mf::LogInfo(
"MuonFilter") <<
"Times: " << trk1Start <<
" "<< trk2Start <<
" "<<trk1End <<
" "<<trk2End;
239 if((TMath::Abs(trk1Start-trk2Start) <
fTolerance && TMath::Abs(trk1End-trk2End) <
fTolerance) &&
247 double y1,
y2, z1, z2;
251 double const x1 = (trk1Start+trk2Start)/2.0*drift-
fDCenter;
252 double const x2 = (trk1End+trk2End)/2.0*drift-
fDCenter;
254 <<
" " << x1 <<
" " << y1 <<
" " << z1
255 <<
" " << x2 <<
" " << y2 <<
" " << z2;
256 bool x1edge,x2edge,y1edge, y2edge,z1edge,z2edge;
266 x1edge =(TMath::Abs(x1) -
fCuts[0] > 0);
267 x2edge =(TMath::Abs(x2) -
fCuts[0] > 0);
268 y1edge =(TMath::Abs(y1) -
fCuts[1] > 0);
269 y2edge =(TMath::Abs(y2) -
fCuts[1] > 0);
270 z1edge =(TMath::Abs(z1) -
fCuts[2] > 0);
271 z2edge =(TMath::Abs(z2) -
fCuts[2] > 0);
272 if((x1edge||y1edge||z1edge) && (x2edge||y2edge||z2edge)) {
273 tGoing.push_back(pointTemp);
274 mf::LogInfo(
"MuonFilter") <<
"outside Removed induction ion: ";
276 for(
size_t h = 0; h < indHits.size(); h++){
277 mf::LogInfo(
"MuonFilter") << indHits[h]->PeakAmplitude() <<
" ";
278 indIon -= indHits[h]->PeakAmplitude();
280 mf::LogInfo(
"MuonFilter") <<
"Removed collection ion: ";
282 for(
size_t h = 0; h < colHits.size(); h++){
283 mf::LogInfo(
"MuonFilter") << colHits[h]->PeakAmplitude() <<
" ";
284 colIon -= colHits[h]->PeakAmplitude();
286 mf::LogInfo(
"MuonFilter")<<
"Ionization outside track I/C: " << indIon <<
" "<<colIon;
288 else if((x1edge || y1edge || z1edge) &&
289 !(x2edge || y2edge || z2edge) &&
291 tGoing.push_back(pointTemp);
292 mf::LogInfo(
"MuonFilter") <<
"stopping Removed induction ion: ";
293 for(
size_t h = 0; h < indHits.size(); h++){
294 mf::LogInfo(
"MuonFilter") <<indHits[h]->PeakAmplitude() <<
" ";
295 indIon -= indHits[h]->PeakAmplitude();
297 mf::LogInfo(
"MuonFilter") <<
"Removed collection ion: ";
298 for(
size_t h = 0; h < colHits.size(); h++){
299 mf::LogInfo(
"MuonFilter") << colHits[h]->PeakAmplitude()<<
" ";
300 colIon -= colHits[h]->PeakAmplitude();
302 mf::LogInfo(
"MuonFilter") <<
"Ionization outside track I/C: " << indIon <<
" "<<colIon;
305 pairTemp = std::make_pair(i,j);
306 mf::LogInfo(
"MuonFilter") <<
"rLook matchnum " <<matchNum <<
" "<<i <<
" "<<j ;
307 rLook.push_back(pairTemp);
308 matched.push_back(pointTemp);
317 for(
unsigned int i = 0; i < tGoing.size(); i++)
318 for(
unsigned int j = 0; j < matched.size();j++){
319 mf::LogInfo(
"MuonFilter") << tGoing.size() <<
" " << matched.size() <<
" " << i <<
" " << j;
321 if((tGoing[i][2] <= matched[j][2]) &&
322 (tGoing[i][5] >= matched[j][5])) {
323 TVector3
a1(&tGoing[i][0]);
324 TVector3
a2(&tGoing[i][3]);
325 TVector3 b1(&matched[j][0]);
326 distance = TMath::Abs((((
a1-
a2).Cross((
a1-
a2).Cross(
a1-b1))).Unit()).Dot(
a1-b1));
327 mf::LogInfo(
"MuonFilter") <<
"distance "<< distance;
329 mf::LogInfo(
"MuonFilter") <<
"Removing delta ion "<<rLook.size()<<
" "<<rLook[j].first<<
" "<<matchNum;
330 std::vector< art::Ptr<recob::Hit> > temp = fmhi.at(rLook[j].first);
331 for(
unsigned int h = 0; h < temp.size();h++)
332 indIon -= temp[h]->PeakAmplitude();
333 temp = fmhc.at(rLook[j].second);
334 for(
unsigned int h = 0; h < temp.size();h++)
335 colIon -= temp[h]->PeakAmplitude();
void reserve(size_type n)
PlaneGeo const & Plane(unsigned int const p, unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified wire.
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]
std::string fLineModuleLabel
float StartWire() const
Returns the wire coordinate of the start of the cluster.
float EndTick() const
Returns the tick coordinate of the end of the cluster.
std::vector< double > fCuts
geo::PlaneID Plane() const
Returns the plane ID this cluster lies on.
Cluster finding and building.
Float_t y2[n_points_geant4]
View_t View() const
Which coordinate does this plane measure.
unsigned int Nplanes(unsigned int tpc=0, unsigned int cstat=0) const
Returns the total number of wire planes in the specified TPC.
std::string fClusterModuleLabel
void push_back(Ptr< U > const &p)
reference at(size_type n)
bool IntersectionPoint(geo::WireID const &wid1, geo::WireID const &wid2, double &y, double &z) const
Returns the intersection point of two wires.
int fDeltaWire
allowed differences in wire number between 2 planes
Detector simulation of raw signals on wires.
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Float_t x2[n_points_geant4]
float StartTick() const
Returns the tick coordinate of the start of the cluster.
float EndWire() const
Returns the wire coordinate of the end of the cluster.