38 #include "Utilities/func/MathUtil.h" 39 #include "Utilities/AssociationUtil.h" 128 produces<std::vector<rb::Cluster> >();
154 assert(0 &&
"Unknown detector");
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 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 188 std::unique_ptr< std::vector<rb::Cluster> > showercol(
new std::vector<rb::Cluster>);
190 showercol->emplace_back();
191 showercol->back().SetNoise(
true);
195 for(
unsigned int iSlice =0; iSlice < slicecol->size(); ++iSlice){
198 if(slc->
NCell()<=0)
continue;
200 bool goodSlice =
true;
202 if (slc->
MeanTNS () <25000 || slc->
MeanTNS () > 475000) {goodSlice =
false;}
204 std::vector <art::Ptr<rb::Track> > tracks = fmt.at(iSlice);
212 bool isShower =
false;
213 bool isCandTrack =
false;
215 if(!slc->
IsNoise() && tracks.size()==1 && tracks[0]->Is3D() && goodSlice){
217 bool startContained =
false;
218 bool endContained =
false;
231 else nPlaneTrack = tracks[0]->MaxPlane() - tracks[0]->MinPlane();
238 if(startContained && !endContained) isCandTrack =
false;
264 if(startPlane>0 && endPlane>0) {
266 int parameter =
std::abs(
int(tracks[0]->MaxPlane()) - endPlane);
269 if(!
eparm(slc,tracks[0],startPlane, endPlane)){
289 muonstub(slc, tracks[0],startPlane );
293 for(
unsigned int iCell =0; iCell<slc->
NCell(); ++iCell){
296 if(chit->
Plane()<startPlane || chit->
Plane()>endPlane)
continue;
305 showercol->push_back(shower);
311 showercol->push_back(*tracks[0]);
315 showercol->push_back(*slc);
321 evt.
put(std::move(showercol));
329 bool isContained =
false;
331 double xyzMax[3], xyzMin[3];
362 bool isContained =
false;
365 for(
int iTrajPt=0; iTrajPt<nTrajPt; iTrajPt++){
373 isContained =
inFiducial(xTrajpt, yTrajpt, zTrajpt);
381 isContained =
inFiducial(xTrajpt, yTrajpt, zTrajpt);
401 std::map<unsigned int, double> ePlane;
405 for(
unsigned int iCell=0; iCell<trk->
NCell(); ++iCell){
410 ePlane[cHit->
Plane()] += rHit.
GeV();
414 for(
int i=pn;
i >=minPlane;
i--){
415 double mipDeDx = 0.00157;
418 double avg = (ePlane[
i-2] + ePlane[
i] + ePlane[
i-1])/3;
422 if(avg<=0.9*mipPlane)
439 std::map<unsigned int, double> ePlane;
443 for(
unsigned int iCell=0; iCell<trk->
NCell(); ++iCell){
448 ePlane[cHit->
Plane()] += rHit.
GeV();
452 for(
int i=pln;
i <=maxPlane-2;
i++){
453 double mipDeDx = 0.00157;
456 double avg = (ePlane[
i+2] + ePlane[
i] + ePlane[
i+1])/3;
459 if(avg<=0.9*mipPlane)
464 if(pln+3 <=
int(trk->
MaxPlane()))
return pln+3;
477 std::map<unsigned int, double> ePlane;
483 for(
unsigned int iCell=0; iCell<trk->
NCell(); ++iCell){
488 ePlane[cHit->
Plane()] += rHit.
GeV();
495 for(
int i=endPlane;
i >endPlane-7;
i--){
496 double mipDeDx = 0.00157;
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);
519 std::map<unsigned int, double> ePlanec;
520 std::map<unsigned int, double> ePlane;
521 std::map<unsigned int, double> ePlaneAve;
526 for(
unsigned int iCell=0; iCell<slc->
NCell(); ++iCell){
533 ePlanec[cHit->
Plane()] += rHit.
GeV();
543 for(
unsigned int iCell=0; iCell<trk->
NCell(); ++iCell){
550 ePlane[cHit->
Plane()] += rHit.
GeV();
558 for(
int iPlane=minPlane; iPlane<=
maxPlane; iPlane++){
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;
575 bool isShower =
false;
577 bool showerEnd =
false;
579 const int nPlaneStart =5;
580 const int nPlaneEnd = 2;
583 for(
int iPlane=minPlane; iPlane<=
maxPlane; iPlane++){
584 double mipDeDx = 0.00157;
587 if(iPlane<=(maxPlane-nPlaneStart))
588 for(
int k=0; k<nPlaneStart; ++k){
591 if(k==0) showerStart= (ePlane[iPlane+k]>2.5*mipPlane);
592 else showerStart *= (ePlane[iPlane+k]>2.5*mipPlane);
593 if(!showerStart)
break;
596 else showerStart =
false;
607 for(
int k=0; k<nPlaneEnd; k++)
608 if((iPlane+k)<=maxPlane){
610 if(k==0) showerEnd = (ePlaneAve[iPlane+k]<1.5*mipPlane && ePlaneAve[iPlane+k]>0.01*mipPlane);
611 else showerEnd *= (ePlaneAve[iPlane+k]<1.5*mipPlane && ePlaneAve[iPlane+k]>0.01*mipPlane);
613 if(!showerEnd)
break;
627 for(
int iPlane = startPlane; iPlane <= endPlane; iPlane++) {
630 eshower += ePlanec[iPlane] ;
646 TVector3
dir = trk->
Dir();
650 const TVector3 dxyz = (cHit->
View() ==
geo::kX) ? TVector3(0, 1, 0): TVector3(1, 0, 0);
667 if(spp>maxPlane) return ;
668 std::map<unsigned int, double> ePlane;
672 for(
unsigned int iCell=0; iCell<slc->
NCell(); ++iCell){
676 ePlane[cHit->
Plane()] += rHit.
GeV();
684 for(
int i=sp;
i <sp+7;
i++){
685 double mipDeDx = 0.00157;
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);
695 if(eparm){ startPlane=
i+4;}
size_t NTrajectoryPoints() const
SubRunNumber_t subRun() const
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.
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
void GetCenter(double *xyz, double localz=0.0) const
TVector3 TrajectoryPoint(unsigned int i) const
The ith point on the trajectory, a 3-vector in cm.
fhicl::ParameterSet fPSetND
fvar< T > fabs(const fvar< T > &x)
int adjustPlane(art::Ptr< rb::Track > &trk, int p)
unsigned short Plane() const
double DistToEdgeY(TVector3 vertex)
const CellGeo * Cell(int icell) const
bool trackStartContained(art::Ptr< rb::Track > &trk)
Vertical planes which measure X.
virtual void produce(art::Event &evt) override
A collection of associated CellHits.
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
void muonstub(art::Ptr< rb::Cluster > &slc, art::Ptr< rb::Track > &trk, int &startplane)
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
bool inFiducial(double x, double y, double z)
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
ProductID put(std::unique_ptr< PROD > &&product)
unsigned short Cell() const
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)
bool fWriteParentSliceCluster
T get(std::string const &key) const
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
virtual TVector3 Dir() const
Unit vector describing prong direction.
rb::RecoHit RecoHit(const art::Ptr< rb::CellHit > &chit) const
Return calibrated hit based on assumed W coordinate.
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.
EventNumber_t event() const
bool fWriteParentTrackCluster
double DetHalfHeight() const
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
virtual bool Is3D() const
bool eparm(art::Ptr< rb::Cluster > &slc, art::Ptr< rb::Track > &trk, int startplane, int endPlane)
A rawdata::RawDigit with channel information decoded.
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.
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
virtual ~DiFShowerFinder()
bool trackEndContained(art::Ptr< rb::Track > &trk)
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
double MeanTNS(rb::AveragingScheme scheme=kDefaultScheme) const
void beginRun(art::Run &)
std::string fClusterLabel
fhicl::ParameterSet fPSetTB
assert(nhit_max >=nhit_nbins)
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
fhicl::ParameterSet fPSetFD
bool IsNoise() const
Is the noise flag set?
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