DiFShowerFinder_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: DiFShowerFinder
3 // Module Type: producer
4 // File: DiFShowerFinder_module.cc
5 //
6 // Generated at Thu Aug 13 13:48:09 2015 by Nitin Yadav using artmod
7 // from cetpkgsupport v1_08_06.
8 //
9 //
10 // Updated to include use at ND 10/30/2019
11 ////////////////////////////////////////////////////////////////////////
12 
20 
21 #include "fhiclcpp/ParameterSet.h"
27 
28 #include "Geometry/Geometry.h"
29 #include "Geometry/LiveGeometry.h"
30 #include "GeometryObjects/Geo.h"
32 
33 #include "RecoBase/Track.h"
34 #include "RecoBase/Cluster.h"
35 #include "RecoBase/CellHit.h"
36 #include "RecoBase/RecoHit.h"
37 
38 #include "Utilities/func/MathUtil.h"
39 #include "Utilities/AssociationUtil.h"
40 
41 #include <vector>
42 #include <algorithm>
43 
44 //-----------------------------------------------------------------------------
45 /////////////////////////////////////////////////////////////////////////////////
46 
47 namespace dif
48 {
49 
51 {
52  public:
53  explicit DiFShowerFinder(fhicl::ParameterSet const & pset);
54  virtual ~DiFShowerFinder();
55 
56  //void reconfigure(const fhicl::ParameterSet& pset); //redid hwo parameters are read in
57  void beginJob() override;
58  void beginRun(art::Run&);
59  virtual void produce(art::Event & evt) override;
60 
61  void findShowerByReco(art::Ptr<rb::Cluster> & slc,art::Ptr<rb::Track> &trk, int &startPlane, int &endPlane);
62  bool inFiducial(double x, double y, double z);
65  int adjustPlane(art::Ptr<rb::Track> & trk, int p);
66  int adjustPlane_end(art::Ptr<rb::Track> & trk, int p);
67  bool eparm(art::Ptr<rb::Cluster> & slc,art::Ptr<rb::Track> & trk, int startplane, int endPlane);
69  void muonstub(art::Ptr<rb::Cluster> & slc,art::Ptr<rb::Track> & trk, int & startplane);
70 
71  void endJob() override;
72 
73 
74  private:
77 
81 
83  bool fTrackForwardGoing; //Assume the track end is the farthest in z instead of lowest in y?
84  float fAngleCut;
85  float fEnergyCut;
86  float fFiducialX;
87  float fFiducialY;
88  float fFiducialZ;
90  int fPlaneMin;
91  int fPlaneMax;
92  int fMinTrackLength;//"Optimized" cut on track quality
93  int fMaxTrackTail; //"Optimized" cut on Brem background
94 
96 
99 };
100 
101 //-----------------------------------------------------------------------------
102 /////////////////////////////////////////////////////////////////////////////////
103 
105 
106 
107  fClusterLabel ( pset.get< std::string >("ClusterLabel")),
108  fTrackLabel ( pset.get< std::string >("TrackLabel")),
109  fWriteDiFCluster ( pset.get< bool >("WriteDiFCluster")),
110  fWriteParentTrackCluster ( pset.get< bool >("WriteParentTrackCluster")),
111  fWriteParentSliceCluster ( pset.get< bool >("WriteParentSliceCluster")),
112  fShowerByReco ( pset.get< bool >("ShowerByReco")),
113  fTrackForwardGoing (false),
114  fAngleCut ( pset.get< float >("AngleCut")),
115  fEnergyCut ( pset.get< float >("EnergyCut")),
116  fFiducialX (0),
117  fFiducialY (0),
118  fFiducialZ (0),
119  fMinDistToEdge (0),
120  fPlaneMin (0),
121  fPlaneMax (0),
122  fMinTrackLength (0),
123  fMaxTrackTail (0),
124  fPSetND (pset.get<fhicl::ParameterSet>("nd")),
125  fPSetFD (pset.get<fhicl::ParameterSet>("fd")),
126  fPSetTB (pset.get<fhicl::ParameterSet>("tb"))
127  {
128  produces<std::vector<rb::Cluster> >();
129  if(fWriteDiFCluster) produces< art::Assns<rb::Cluster, rb::Track> >();
130  }
131 //-----------------------------------------------------------------------------
132 /////////////////////////////////////////////////////////////////////////////////
133 
135 
137  {
138  }
139 
141  {
142  fhicl::ParameterSet pset;
143  switch(geom->DetId()){
145  pset = fPSetND;
146  break;
148  pset = fPSetFD;
149  break;
151  pset = fPSetTB;
152  break;
153  default:
154  assert(0 && "Unknown detector");
155  }
156 
157  fTrackForwardGoing = pset.get< bool >("TrackForwardGoing");
158  fFiducialX = pset.get< float >("FiducialX");
159  fFiducialY = pset.get< float >("FiducialY");
160  fFiducialZ = pset.get< float >("FiducialZ");
161  fMinDistToEdge = pset.get< float >("MinDistToEdge");
162  fPlaneMin = pset.get< int >("PlaneMin");
163  fPlaneMax = pset.get< int >("PlaneMax");
164  fMinTrackLength = pset.get< int >("MinTrackLength");
165  fMaxTrackTail = pset.get< int >("MaxTrackTail");
166 
168  mf::LogWarning("DiFShowerFinder") << "DifShowerFinder: WARNING! your configured options will include some hits in more than one output cluster. Recommend WriteParentSliceCluster is not used in combination with the other options" << std::endl
169  << pset.to_string();
170  }
171 
173  mf::LogWarning("DiFShowerFinder") << "DifShowerFinder: WARNING! you haven't configured any output. You probably want one or more of the Write*Cluster settings." << std::endl
174  << pset.to_string();
175  }
176 
177  }
178 
179 
181 {
183  evt.getByLabel(fClusterLabel, slicecol);
184 
185  art::FindManyP<rb::Track> fmt(slicecol, evt, fTrackLabel);
186 
187  std::unique_ptr< art::Assns<rb::Cluster, rb::Track> > trackAssassin(new art::Assns<rb::Cluster, rb::Track>);
188  std::unique_ptr< std::vector<rb::Cluster> > showercol(new std::vector<rb::Cluster>);
189  // Downstream modules love their noise slices
190  showercol->emplace_back();
191  showercol->back().SetNoise(true);
192 
193 
194  // Loop over slices
195  for( unsigned int iSlice =0; iSlice < slicecol->size(); ++iSlice){
196  art::Ptr<rb::Cluster> slc(slicecol, iSlice);
197 
198  if(slc->NCell()<=0) continue;
199 
200  bool goodSlice = true;
201 
202  if (slc->MeanTNS () <25000 || slc->MeanTNS () > 475000) {goodSlice =false;}
203 
204  std::vector <art::Ptr<rb::Track> > tracks = fmt.at(iSlice);
205 
206  int startPlane = -1;
207  int endPlane = -1;
208 
209  double disToBadX=-1;
210  double disToBadY=-1;
211  double disToBadZ=-1;
212  bool isShower = false;
213  bool isCandTrack = false;
214 
215  if(!slc->IsNoise() && tracks.size()==1 && tracks[0]->Is3D() && goodSlice){
216 
217  bool startContained = false;
218  bool endContained =false;
219 
220  //Decide if this track is a candidat track?
221  if(int(tracks[0]->MaxPlane() ) >fPlaneMin && int(tracks[0]->MinPlane() ) <fPlaneMax){
222 
223  startContained = trackStartContained(tracks[0]);
224  endContained = trackEndContained(tracks[0]);
225 
226  //part of "optimized" cuts
227  int nPlaneTrack = 0;
228 
229  if(int(tracks[0]->MaxPlane())>fPlaneMax) nPlaneTrack = fPlaneMax - int(tracks[0]->MinPlane());
230  else if(int(tracks[0]->MinPlane())<fPlaneMin) nPlaneTrack = int(tracks[0]->MaxPlane()) - fPlaneMin;
231  else nPlaneTrack = tracks[0]->MaxPlane() - tracks[0]->MinPlane();
232 
233  disToBadX=livegeom->DistToEdgeX(tracks[0]->Stop());
234  disToBadY=livegeom->DistToEdgeY(tracks[0]->Stop());
235  disToBadZ=livegeom->DistToEdgeZ(tracks[0]->Stop());
236 
237 
238  if(startContained && !endContained) isCandTrack =false;
239  else{
240 
241  if(tracks[0]->Dir().Z()<fAngleCut) isCandTrack = false;
242  else {
243 
244  if(nPlaneTrack<fMinTrackLength) isCandTrack =false; //"optimized" cut
245  else {
246 
247  if(abs(disToBadX)<fMinDistToEdge || abs(disToBadY)<fMinDistToEdge || abs(disToBadZ)<fMinDistToEdge) isCandTrack = false;
248  else {
249  isCandTrack = true;
250  }
251  }
252  }
253 
254 
255  }
256  } // candidate track ends
257 
258 
259  if(isCandTrack){
260 
261  if(fShowerByReco && tracks[0]->Dir().Z()>0){
262 
263  findShowerByReco(slc, tracks[0], startPlane, endPlane);
264  if(startPlane>0 && endPlane>0) {
265 
266  int parameter = std::abs(int(tracks[0]->MaxPlane()) - endPlane);
267 
268  if(parameter<fMaxTrackTail){ //"optimized" cut
269  if(!eparm(slc,tracks[0],startPlane, endPlane)){
270 
271  isShower=true;
272  // endPlane = tracks[0]->MaxPlane();
273 
274  } //eparm
275  } //param
276  } // startplane and endplane >0
277 
278  } // fShowerByReco
279 
280  } // isCandTrack
281  } //good slice ends end if track size is 1
282 
283  if(isShower){
284 
285  std::cout<<"DiF Shower Found! Run: "<<evt.run()<<" Subrun: "<<evt.subRun()<<" Event: "<<evt.event()<<"\n";
286 
287  if(fWriteDiFCluster){
288  // removing the muon stub
289  muonstub(slc, tracks[0],startPlane );
290 
291  rb::Cluster shower;
292 
293  for(unsigned int iCell =0; iCell<slc->NCell(); ++iCell){
294  art::Ptr<rb::CellHit> chit = slc->Cell(iCell);
295 
296  if(chit->Plane()<startPlane || chit->Plane()>endPlane) continue;
297 
298  if(distancefromtrack(chit, tracks[0])<=80){
299  shower.Add(chit);
300  }
301  }//looping cells
302 
303  //non 3D clusters aren't useful and break things down the line
304  if(shower.Is3D()){
305  showercol->push_back(shower);
306  if(fWriteDiFCluster) util::CreateAssn( *this, evt, *showercol, tracks[0], *trackAssassin );
307  }
308  }
309 
311  showercol->push_back(*tracks[0]);
312  }
313 
315  showercol->push_back(*slc);
316  }
317  } //isShower
318 
319 } // looping over slice ends
320 
321  evt.put(std::move(showercol));
322  if(fWriteDiFCluster) evt.put(std::move(trackAssassin));
323 }//end produce
324 
325 //-----------------------------------------------------------------------------
326 ///////////////////////////////////////////////////////////////////////////////
327 
328 bool DiFShowerFinder::inFiducial(double x, double y, double z){
329  bool isContained = false;
330 
331  double xyzMax[3], xyzMin[3];
332  geom->Plane(fPlaneMax)->Cell(0)->GetCenter(xyzMax);
333  geom->Plane(fPlaneMin)->Cell(0)->GetCenter(xyzMin);
334  double zMax = xyzMax[2] + geom->Plane(fPlaneMax)->Cell(0)->HalfD();
335  double zMin = xyzMin[2] - geom->Plane(fPlaneMin)->Cell(0)->HalfD();
336 
337  if(x < (geom->DetHalfWidth() - fFiducialX) && x > (-geom->DetHalfWidth() + fFiducialX) &&
338  y < (geom->DetHalfHeight() - fFiducialY) && y > (-geom->DetHalfHeight() + fFiducialY) &&
339  z < (zMax - fFiducialZ) && z > zMin + fFiducialZ
340  )
341 
342  isContained = true;
343  return isContained;
344 
345 }
346 //-----------------------------------------------------------------------------
347 /////////////////////////////////////////////////////////////////////////////////
348 
350 
351  return inFiducial(trk->Start().X(), trk->Start().Y(), trk->Start().Z());
352 }
353 
354 //-----------------------------------------------------------------------------
355 /////////////////////////////////////////////////////////////////////////////////
356 
358 
359  double ymin = 10000;
360  double zmax = 0;
361 
362  bool isContained = false;
363  int nTrajPt = trk->NTrajectoryPoints();
364 
365  for(int iTrajPt=0; iTrajPt<nTrajPt; iTrajPt++){
366  double xTrajpt = trk->TrajectoryPoint(iTrajPt).X();
367  double yTrajpt = trk->TrajectoryPoint(iTrajPt).Y();
368  double zTrajpt = trk->TrajectoryPoint(iTrajPt).Z();
369 
370  //assume forward track if specified
371  if(fTrackForwardGoing == true){
372  if(zTrajpt>zmax){
373  isContained = inFiducial(xTrajpt, yTrajpt, zTrajpt);
374  zmax = zTrajpt;
375  }
376  }
377 
378  //assume downward track otherwise
379  else{
380  if(yTrajpt<ymin){
381  isContained = inFiducial(xTrajpt, yTrajpt, zTrajpt);
382  ymin = yTrajpt;
383  }
384  }
385 
386  }
387 
388  return isContained;
389 }
390 //-----------------------------------------------------------------------------
391 /////////////////////////////////////////////////////////////////////////////////
392 
394 
395  int pn=p;
396 
397  int minPlane = trk->MinPlane();
398  int maxPlane = trk->MaxPlane();
399  if(minPlane<fPlaneMin) minPlane = fPlaneMin;
400  if(maxPlane>fPlaneMax) maxPlane = fPlaneMax;
401  std::map<unsigned int, double> ePlane;
402 
403  ePlane.clear();
404 
405  for(unsigned int iCell=0; iCell<trk->NCell(); ++iCell){
406  const rb::CellHit* cHit = trk->Cell(iCell).get();
407  const rb::RecoHit rHit = trk->RecoHit(trk->Cell(iCell));
408  if(cHit->Plane()<fPlaneMin || cHit->Plane()>fPlaneMax) continue;
409  if( rHit.IsCalibrated() ){
410  ePlane[cHit->Plane()] += rHit.GeV();
411  }
412  }//end loop over cells
413 
414  for(int i=pn; i >=minPlane; i--){
415  double mipDeDx = 0.00157;
416  double mipPlane = 2*mipDeDx*geom->Plane(i)->Cell(0)->HalfD()/fabs(trk->Dir().Z());
417  pn=i;
418  double avg = (ePlane[i-2] + ePlane[i] + ePlane[i-1])/3;//the average of 3 planes
419  //double avg = (ePlane[i-3] + ePlane[i-2] + ePlane[i-1]+ ePlane[i])/4;//the average of 3 planes
420  //if(ePlane[i]<0.5*mipPlane) break;
421  //if(avg<=1.0*mipPlane)
422  if(avg<=0.9*mipPlane)
423 
424  break;
425  }
426 
427  return pn;
428 }
429 //-----------------------------------------------------------------------------
430 /////////////////////////////////////////////////////////////////////////////////
431 
433  int pln=p;
434 
435  int minPlane = trk->MinPlane();
436  int maxPlane = trk->MaxPlane();
437  if(minPlane<fPlaneMin) minPlane = fPlaneMin;
438  if(maxPlane>fPlaneMax) maxPlane = fPlaneMax;
439  std::map<unsigned int, double> ePlane;
440 
441  ePlane.clear();
442 
443  for(unsigned int iCell=0; iCell<trk->NCell(); ++iCell){
444  const rb::CellHit* cHit = trk->Cell(iCell).get();
445  const rb::RecoHit rHit = trk->RecoHit(trk->Cell(iCell));
446  if(cHit->Plane()<fPlaneMin || cHit->Plane()>fPlaneMax) continue;
447  if( rHit.IsCalibrated() ){
448  ePlane[cHit->Plane()] += rHit.GeV();
449  }
450  }//end loop over cells
451 
452  for(int i=pln; i <=maxPlane-2; i++){
453  double mipDeDx = 0.00157;
454  double mipPlane = 2*mipDeDx*geom->Plane(i)->Cell(0)->HalfD()/fabs(trk->Dir().Z());
455  pln=i;
456  double avg = (ePlane[i+2] + ePlane[i] + ePlane[i+1])/3;//the average of 3 planes
457  //double avg = (ePlane[i-3] + ePlane[i-2] + ePlane[i-1]+ ePlane[i])/4;//the average of 3 planes
458  //if(ePlane[i]<0.5*mipPlane) break;
459  if(avg<=0.9*mipPlane)
460 
461  break;
462  }
463 
464  if(pln+3 <= int(trk->MaxPlane())) return pln+3;
465  else return pln;
466 }
467 
468 //-----------------------------------------------------------------------------
469 /////////////////////////////////////////////////////////////////////////////////
470 
471 bool DiFShowerFinder::eparm(art::Ptr<rb::Cluster> & slc,art::Ptr<rb::Track> & trk, int startPlane , int endPlane){
472 
473  int minPlane = trk->MinPlane();
474  int maxPlane = trk->MaxPlane();
475  if(minPlane<fPlaneMin) minPlane = fPlaneMin;
476  if(maxPlane>fPlaneMax) maxPlane = fPlaneMax;
477  std::map<unsigned int, double> ePlane;
478  // std::map<unsigned int, double> ePlanec;
479 
480  ePlane.clear();
481  //ePlanec.clear();
482 
483  for(unsigned int iCell=0; iCell<trk->NCell(); ++iCell){
484  const rb::CellHit* cHit = trk->Cell(iCell).get();
485  const rb::RecoHit rHit = trk->RecoHit(trk->Cell(iCell));
486  if(cHit->Plane()<fPlaneMin || cHit->Plane()>fPlaneMax) continue;
487  if( rHit.IsCalibrated() ){
488  ePlane[cHit->Plane()] += rHit.GeV();
489  }
490  }//end loop over cells
491 
492  bool eparm=false;
493  //int eparmcount=0;
494 
495  for(int i=endPlane; i >endPlane-7; i--){
496  double mipDeDx = 0.00157;
497  double mipPlane = 2*mipDeDx*geom->Plane(i)->Cell(0)->HalfD()/fabs(trk->Dir().Z());
498 
499  for(int k=0;k<=3;k++){
500  if(k==0) eparm=(ePlane[i-k] >0.5*mipPlane && ePlane[i-k]<=1.5*mipPlane);
501  else eparm*= (ePlane[i-k] >0.5*mipPlane && ePlane[i-k]<=1.5*mipPlane);
502  if(!eparm) break;
503  }
504  if(eparm) break;
505  }
506 
507  return eparm;
508 }
509 //-----------------------------------------------------------------------------
510 /////////////////////////////////////////////////////////////////////////////////
511 
513 
514  int minPlane = trk->MinPlane();
515  int maxPlane = trk->MaxPlane();
516  if(minPlane<fPlaneMin) minPlane = fPlaneMin;
517  if(maxPlane>fPlaneMax) maxPlane = fPlaneMax;
518 
519  std::map<unsigned int, double> ePlanec;
520  std::map<unsigned int, double> ePlane;
521  std::map<unsigned int, double> ePlaneAve;
522 
523  ePlane.clear();
524  ePlanec.clear();
525 
526 for(unsigned int iCell=0; iCell<slc->NCell(); ++iCell){
527  const rb::CellHit* cHit = slc->Cell(iCell).get();
528  const rb::RecoHit rHit = slc->RecoHit(slc->Cell(iCell));
529 
530  if(cHit->Plane()<fPlaneMin || cHit->Plane()>fPlaneMax) continue;
531 
532  if( rHit.IsCalibrated() ){
533  ePlanec[cHit->Plane()] += rHit.GeV();
534  // fHEdep->Fill(rHit.GeV());
535  // fHADC->Fill(cHit->ADC(0));
536  // fHPE->Fill(cHit->PE());
537  // fHPECorr->Fill(rHit.PECorr());
538  }
539  }//end loop over cells
540 
541 
542 
543  for(unsigned int iCell=0; iCell<trk->NCell(); ++iCell){
544  const rb::CellHit* cHit = trk->Cell(iCell).get();
545  const rb::RecoHit rHit = trk->RecoHit(trk->Cell(iCell));
546 
547  if(cHit->Plane()<fPlaneMin || cHit->Plane()>fPlaneMax) continue;
548 
549  if( rHit.IsCalibrated() ){
550  ePlane[cHit->Plane()] += rHit.GeV();
551  // fHEdep->Fill(rHit.GeV());
552  // fHADC->Fill(cHit->ADC(0));
553  // fHPE->Fill(cHit->PE());
554  // fHPECorr->Fill(rHit.PECorr());
555  }
556  }//end loop over cells
557 
558  for(int iPlane=minPlane; iPlane<=maxPlane; iPlane++){
559  //opt stuff
560  //if(iPlane == minPlane) ePlaneAve[iPlane] = 0.5*(ePlane[iPlane] + ePlane[iPlane+1]);
561  //else if(iPlane == maxPlane) ePlaneAve[iPlane] = 0.5*(ePlane[iPlane-1] + ePlane[iPlane]);
562  // else ePlaneAve[iPlane] = (ePlane[iPlane-1] + ePlane[iPlane] + ePlane[iPlane+1])/3;//the average of 3 planes
563 
564  //ePlaneAve[iPlane] = ePlane[iPlane];
565 
566 
567  if(iPlane == minPlane) ePlaneAve[iPlane] = (ePlane[iPlane] + ePlane[iPlane+1] + ePlane[iPlane+2])/3;
568  else if(iPlane == minPlane+1) ePlaneAve[iPlane] = ( ePlane[iPlane-1] + ePlane[iPlane] + ePlane[iPlane+1]+ePlane[iPlane+2] )/4;
569  else if(iPlane == maxPlane) ePlaneAve[iPlane] = (ePlane[iPlane-2] + ePlane[iPlane-1] + ePlane[iPlane])/3;
570  else if(iPlane == maxPlane-1) ePlaneAve[iPlane] = (ePlane[iPlane-2]+ePlane[iPlane-1] + ePlane[iPlane] + ePlane[iPlane+1])/4;
571  else ePlaneAve[iPlane] = (ePlane[iPlane-2] + ePlane[iPlane-1] + ePlane[iPlane] + ePlane[iPlane+1] + ePlane[iPlane+2])/5;//the average of 5 planes
572 
573  }
574 
575  bool isShower = false;
576  bool showerStart = false;
577  bool showerEnd = false;
578 
579  const int nPlaneStart =5;// it is 2
580  const int nPlaneEnd = 2;//opt stuff change from 2 to 1
581 
582  double eshower = 0;
583  for(int iPlane=minPlane; iPlane<=maxPlane; iPlane++){
584  double mipDeDx = 0.00157;
585  double mipPlane = 2*mipDeDx*geom->Plane(iPlane)->Cell(0)->HalfD()/fabs(trk->Dir().Z());
586  if(!isShower){//if it's not already in a shower, try to find one
587  if(iPlane<=(maxPlane-nPlaneStart))
588  for(int k=0; k<nPlaneStart; ++k){
589  // if(k==0) showerStart= (ePlane[iPlane+k]>3.5*mipPlane);
590  // else showerStart *= (ePlane[iPlane+k]>3.5*mipPlane);
591  if(k==0) showerStart= (ePlane[iPlane+k]>2.5*mipPlane);//opt stuff 3.0 to 1.5
592  else showerStart *= (ePlane[iPlane+k]>2.5*mipPlane); //opt stuff
593  if(!showerStart) break;
594  }
595 
596  else showerStart = false;
597 
598 
599  if(showerStart){
600  isShower = true;
601  startPlane=adjustPlane(trk, iPlane ) ;
602  // startPlane = iPlane;
603  }//end if shower start
604  }//end if not shower
605 
606  else{
607  for(int k=0; k<nPlaneEnd; k++)
608  if((iPlane+k)<=maxPlane){
609 
610  if(k==0) showerEnd = (ePlaneAve[iPlane+k]<1.5*mipPlane && ePlaneAve[iPlane+k]>0.01*mipPlane); //opt stuff change from 1.5 to 2.5
611  else showerEnd *= (ePlaneAve[iPlane+k]<1.5*mipPlane && ePlaneAve[iPlane+k]>0.01*mipPlane);
612 
613  if(!showerEnd) break;
614  }
615  if(showerEnd){
616  isShower = false;
617  endPlane= adjustPlane_end(trk,iPlane);
618 // endPlane = iPlane;
619  }
620 
621  }//end else
622 
623  }//end loop over track planes
624 
625 
626  if(startPlane>0 && endPlane>0 && startPlane>fPlaneMin && endPlane<fPlaneMax){
627  for(int iPlane = startPlane; iPlane <= endPlane; iPlane++) {
628  //double mipDeDx = 0.00157;
629  //double mipPlane = 2*mipDeDx*geom->Plane(iPlane)->Cell(0)->HalfD()/fabs(trk->Dir().Z());
630  eshower += ePlanec[iPlane] ;//- mipPlane;
631  }
632 
633  // fHEShower->Fill(eshower);
634  }
635 
636  if(eshower < fEnergyCut) startPlane = -1;
637  return;
638 }//end findShowerByReco
639 
640 //-----------------------------------------------------------------------------
641 /////////////////////////////////////////////////////////////////////////////////
642 
644 {
645  TVector3 start = trk->Start();
646  TVector3 dir = trk->Dir();
647 
648  TVector3 xyz;
649  geom->Plane(cHit->Plane())->Cell(cHit->Cell())->GetCenter(xyz);
650  const TVector3 dxyz = (cHit->View() == geo::kX) ? TVector3(0, 1, 0): TVector3(1, 0, 0);
651 
652  double di;
653  return sqrt(geo::ClosestApproach(start, start+dir, xyz, xyz+dxyz, &di));
654 
655 } // end of DistanceToCore
656 
657 //-----------------------------------------------------------------------------
658 /////////////////////////////////////////////////////////////////////////////////
659 
661 
662  int sp=startPlane;
663  //int minPlane = slc->MinPlane();
664  int maxPlane = slc->MaxPlane();
665 
666  int spp=sp+7;
667  if(spp>maxPlane) return ;
668  std::map<unsigned int, double> ePlane;
669 
670  ePlane.clear();
671 
672  for(unsigned int iCell=0; iCell<slc->NCell(); ++iCell){
673  const rb::CellHit* cHit = slc->Cell(iCell).get();
674  const rb::RecoHit rHit = slc->RecoHit(slc->Cell(iCell));
675  if( rHit.IsCalibrated() ){
676  ePlane[cHit->Plane()] += rHit.GeV();
677  }
678  }//end loop over cells
679 
680  bool eparm=false;
681  //int eparmcount=0;
682 
683 
684  for(int i=sp; i <sp+7; i++){
685  double mipDeDx = 0.00157;
686  double mipPlane = 2*mipDeDx*geom->Plane(i)->Cell(0)->HalfD()/fabs(trk->Dir().Z());
687 
688  for(int k=0;k<=4;k++){
689  if(i+k>maxPlane) continue;
690  if(k==0) eparm=(ePlane[i+k] >0.2*mipPlane && ePlane[i+k]<=2.5*mipPlane);
691  else eparm*= (ePlane[i+k] >0.2*mipPlane && ePlane[i+k]<=2.5*mipPlane);
692  if(!eparm) break;
693  }
694 
695  if(eparm){ startPlane=i+4;}
696  }
697 }
698 
699 //-----------------------------------------------------------------------------
700 /////////////////////////////////////////////////////////////////////////////////
701 
703 
704 //-----------------------------------------------------------------------------
705 /////////////////////////////////////////////////////////////////////////////////
706 
708 
709 } // end namespace
size_t NTrajectoryPoints() const
Definition: Track.h:83
SubRunNumber_t subRun() const
Definition: Event.h:72
static bool CreateAssn(art::EDProducer const &prod, art::Event &evt, std::vector< T > &a, art::Ptr< U > b, art::Assns< T, U > &assn, size_t indx=UINT_MAX, std::string const &instance=std::string())
Create a 1 to 1 association between a new product and one already in the event.
double HalfD() const
Definition: CellGeo.cxx:205
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
TVector3 TrajectoryPoint(unsigned int i) const
The ith point on the trajectory, a 3-vector in cm.
Definition: Track.cxx:158
fhicl::ParameterSet fPSetND
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
int adjustPlane(art::Ptr< rb::Track > &trk, int p)
unsigned short Plane() const
Definition: CellHit.h:39
double DistToEdgeY(TVector3 vertex)
geo::View_t View() const
Definition: CellHit.h:41
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
bool trackStartContained(art::Ptr< rb::Track > &trk)
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
Vertical planes which measure X.
Definition: PlaneGeo.h:28
virtual void produce(art::Event &evt) override
A collection of associated CellHits.
Definition: Cluster.h:47
void abs(TH1 *hist)
double DistToEdgeX(TVector3 vertex)
const PlaneGeo * Plane(unsigned int i) const
art::ServiceHandle< geo::LiveGeometry > livegeom
int adjustPlane_end(art::Ptr< rb::Track > &trk, int p)
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
void muonstub(art::Ptr< rb::Cluster > &slc, art::Ptr< rb::Track > &trk, int &startplane)
float abs(float number)
Definition: d0nt_math.hpp:39
Definition: Run.h:31
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
bool inFiducial(double x, double y, double z)
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
Float_t Z
Definition: plot.C:38
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
unsigned short Cell() const
Definition: CellHit.h:40
Track finder for cosmic rays.
double distancefromtrack(art::Ptr< rb::CellHit > cHit, art::Ptr< rb::Track > &trk)
Far Detector at Ash River, MN.
double DistToEdgeZ(TVector3 vertex)
T get(std::string const &key) const
Definition: ParameterSet.h:231
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
rb::RecoHit RecoHit(const art::Ptr< rb::CellHit > &chit) const
Return calibrated hit based on assumed W coordinate.
Definition: Cluster.cxx:259
Near Detector in the NuMI cavern.
Collect Geo headers and supply basic geometry functions.
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
Definition: RecoHit.cxx:35
EventNumber_t event() const
Definition: Event.h:67
double DetHalfHeight() const
z
Definition: test.py:28
OStream cout
Definition: OStream.cxx:6
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:462
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
virtual bool Is3D() const
Definition: Cluster.h:95
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
bool eparm(art::Ptr< rb::Cluster > &slc, art::Ptr< rb::Track > &trk, int startplane, int endPlane)
A rawdata::RawDigit with channel information decoded.
Definition: CellHit.h:27
float GeV() const
Definition: RecoHit.cxx:69
double DetHalfWidth() const
double ClosestApproach(const double point[], const double intercept[], const double slopes[], double closest[])
Find the distance of closest approach between point and line.
Definition: Geo.cxx:222
T const * get() const
Definition: Ptr.h:321
TDirectory * dir
Definition: macro.C:5
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
bool trackEndContained(art::Ptr< rb::Track > &trk)
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
double MeanTNS(rb::AveragingScheme scheme=kDefaultScheme) const
Definition: Cluster.cxx:554
fhicl::ParameterSet fPSetTB
assert(nhit_max >=nhit_nbins)
Double_t ymin
Definition: plot.C:24
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:508
double showerStart[3]
fhicl::ParameterSet fPSetFD
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
RunNumber_t run() const
Definition: Event.h:77
art::ServiceHandle< geo::Geometry > geom
Encapsulate the geometry of one entire detector (near, far, ndos)
void findShowerByReco(art::Ptr< rb::Cluster > &slc, art::Ptr< rb::Track > &trk, int &startPlane, int &endPlane)
DiFShowerFinder(fhicl::ParameterSet const &pset)
std::string to_string() const
Definition: ParameterSet.h:137