21 #include "Utilities/func/MathUtil.h" 23 #include "TGeoManager.h" 76 double* vertexRegionDeDxCutOff)
95 bool xViewCleaned =
false;
96 bool yViewCleaned =
false;
105 std::vector<rb::WeightedHit > trackHits;
109 trackHits.push_back( newHit );
112 for(
int plidx = 0; plidx < nPlanes; plidx++ ){
113 if ( xViewCleaned && yViewCleaned )
120 if (iPlane > vertexRegion[0] ){
126 if ( iPlane > vertexRegion[1] ){
149 for(
int iHit = 0; iHit< nPlaneHits; iHit++){
166 for(
int k = iHit+1; k<nPlaneHits; k++){
199 unsigned int cell[ nPlaneHits ];
200 std::vector< art::Ptr<rb::CellHit> > slicePlaneHits;
202 int nCells = slice.
NCell();
204 for(
int h = 0;
h < nPlaneHits;
h++)
208 for(
int iHit = 0; iHit < nPlaneHits; iHit++){
211 slicePlaneHits.push_back(slice.
Cell(
hit));
218 int nSlicePlaneHits = slicePlaneHits.size();
219 for(
int hit = 0;
hit< nSlicePlaneHits;
hit++){
221 bool closeEnough =
false;
222 for(
int c = 0;
c < nPlaneHits;
c++ ){
223 int diff = slicePlaneHits[
hit]->Cell() - cell[
c];
224 if ( diff < 5 && diff > -5 )
236 trialE += rHit.
MIP();
244 trackHits.push_back(newHit);
253 trackHits.push_back(newHit);
260 else if (cellFlag == 2){
262 float tempEnergy = eDep + rHit.
MIP() - trialE;
267 trackHits.push_back( tempHit );
325 int minPlane = muonTrack.
MinPlane();
333 if(minPlane == maxPlane){
345 for(
auto const& iHit : trkHits){
358 for(
auto const& iHit : trkHits){
359 int iPlane = iHit->Plane();
381 for(
int iPl = 0, nPlanes =
fPlanes.size(); iPl < nPlanes; ++iPl){
388 if( muonTrack.
Dir().Z() < 0.2 && muonTrack.
Dir().Z() > -0.2 ){
401 double cellCenter[3];
403 double xbound1, ybound1, xbound2, ybound2;
404 muonTrack.
InterpolateXY( cellCenter[2]-planeHalfWidth, xbound1, ybound1);
405 muonTrack.
InterpolateXY( cellCenter[2]+planeHalfWidth, xbound2, ybound2);
406 TVector3 bound1( xbound1, ybound1, cellCenter[2]-planeHalfWidth);
407 TVector3 bound2( xbound2, ybound2, cellCenter[2]+planeHalfWidth);
431 int nTrackHits = trackHits.size();
432 for(
int iHit = 0; iHit < nTrackHits ; iHit++){
433 if(*(trackHits[iHit].hit) == *(hit))
434 trackHits[iHit].weight =
weight;
449 double* vertexRegionDeDxCutOff)
451 if(vertexRegionDeDxCutOff)
454 if( !(muonTrack.
NCell() > 0) )
460 vertexRegion[0] = -1;
461 vertexRegion[1] = -1;
465 std::vector<int> planes;
472 int nPlanes = planes.size();
474 for(
int ipl = 0; ipl < nPlanes - 2 ; ipl++ ){
476 for(
int jpl = ipl; jpl < ipl+3; jpl++){
480 fAvgDeDx[ planes[ipl] ] = avgDeDx/3;
484 for(
int ipl = 0; ipl < nPlanes - 3; ipl++ ){
490 vertexRegion[
view] = planes[ipl];
494 if( vertexRegion[
view] == -1 && nPlanes>0 )
495 vertexRegion[
view] = planes[ nPlanes - 1 ];
498 for(
int ipl = 0; ipl < nPlanes - 3; ipl++ ){
503 vertexRegion[
view] = planes[ipl];
507 if( vertexRegion[
view] == -1 && nPlanes>0 )
508 vertexRegion[
view] = planes[ nPlanes - 1 ];
523 int nHits =
hits.size();
525 for(
int iHit = 0; iHit < nHits ; iHit++){
527 TVector3 cellCenter, trkCenter;
535 TVector3 trkPoint( x, y, trkCenter[2]);
537 float iCellDist = ( cellCenter - trkPoint ).
Mag();
539 for(
int jHit = iHit+1; jHit < nHits; jHit++){
542 float jCellDist = ( cellCenter - trkPoint ).
Mag();
544 if( jCellDist < iCellDist ){
547 hits[jHit] = tmpCellHit;
548 iCellDist = jCellDist;
562 double* mipRangeHigh,
565 double* vertexRegionDeDxCutOff)
567 CleanUpTrack( track, slice, mipRangeHigh, mipRangeLow, mipValue, vertexRegionDeDxCutOff);
573 double* mipRangeHigh,
576 double* vertexRegionDeDxCutOff)
587 if(vertexRegionDeDxCutOff)
592 double mipInMIPs = 1.0;
595 CleanUpTrack( track, slice, &highInMIPs, &lowInMIPs, &mipInMIPs, &vertexRegionInMIPs);
604 double* mipRangeHigh,
607 double* vertexRegionDeDxCutOff)
610 CleanUpTrack( track, slice, mipRangeHigh, mipRangeLow, mipValue, vertexRegionDeDxCutOff);
616 double* mipRangeHigh,
619 double* vertexRegionDeDxCutOff)
630 if(vertexRegionDeDxCutOff)
635 double mipInMIPs = 1.0;
638 CleanUpTrack( track, slice, &highInMIPs, &lowInMIPs, &mipInMIPs, &vertexRegionInMIPs);
639 std::map<int, float> fExtraEOnTrackPlanesInGeV;
641 fExtraEOnTrackPlanesInGeV.insert ( std::pair<int, float>(itr->first, itr->second / 94.62) );
643 return fExtraEOnTrackPlanesInGeV;
650 double* vertexRegionDeDxCutOff)
654 vertexPlanes[0] = vertexRegion[0];
655 vertexPlanes[1] = vertexRegion[1] ;
662 double* mipRangeHigh,
665 double* vertexRegionDeDxCutOff)
668 CleanUpTrack( track, slice, mipRangeHigh, mipRangeLow, mipValue, vertexRegionDeDxCutOff);
677 for(
unsigned int i = 0;
i < track.
NCell(); ++
i){
std::map< int, float > fExtraEOnTrackPlanes
A map of planes in the vertex region of a track and the extra energy found on them.
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
TrackCleanUpAlg(fhicl::ParameterSet const &pset)
void GetCenter(double *xyz, double localz=0.0) const
std::map< int, std::vector< art::Ptr< rb::CellHit > > > fPlaneHits
A map to store the cell numbers that are in each plane of the track.
std::map< geo::OfflineChan, double > fCellEnergy
A map of cell and energy in Plane left by the track.
void reconfigure(const fhicl::ParameterSet &pset)
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
unsigned short Plane() const
const CellGeo * Cell(int icell) const
Vertical planes which measure X.
std::map< unsigned int, double > fPlaneDeDx
A map of plane and dE/dX of the track in the plane.
A collection of associated CellHits.
std::vector< rb::WeightedHit > fCleanUpWeightedHits
A vector of rb::WeightedHit that are impacted by CleanUp procedure.
float fMipRangeHigh
Default upper bound of allowed MIP range (FHiCL parameter)
art::ServiceHandle< geo::Geometry > fgeom
Handle to Geometry service.
A rb::Prong with full reconstructed trajectory.
std::vector< rb::WeightedHit > ResetHitWeight(std::vector< rb::WeightedHit > &trackHits, art::Ptr< rb::CellHit > &hit, float &weight)
A function to reset the weight of a rb::WeightedHit in a given vector of WeightedHits.
const PlaneGeo * Plane(unsigned int i) const
float fVertexRegionDeDxCutOffInGev
the same as fVertexRegionDeDxCutOff but in GeV/cm units (FHiCL parameter)
void ComputeVertexRegion(const rb::Track &muonTrack, int *vertexRegion, double *vertexRegionDeDxCutOff=NULL)
A function to decide how large a vertex region should be considered for clean up. vertexRegion is an ...
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
std::map< int, float > ExtraEOnTrackPlanes(const rb::Track &track, const rb::Cluster &slice, double *mipRangeHigh=NULL, double *mipRangeLow=NULL, double *mipValue=NULL, double *vertexRegionDeDxCutOff=NULL)
Calculates the total hadronic energy on the provided track based on the assumption of MIP for each pl...
virtual double TotalLength() const
Length (cm) of all the track segments.
int IsMIP(double dEdX) const
A function to check if the given dE/dx is consistent with MIP.
std::map< int, float > ExtraEOnTrackPlanesInGeV(const rb::Track &track, const rb::Cluster &slice, double *mipRangeHigh=NULL, double *mipRangeLow=NULL, double *mipValue=NULL, double *vertexRegionDeDxCutOff=NULL)
float fVertexRegionDeDxCutOff
Default vertex region dedx cut off value (FHiCL parameter)
View_t View() const
Which coordinate does this plane measure.
float ExtraEOnTrack(const rb::Track &track, const rb::Cluster &slice, double *mipRangeHigh=NULL, double *mipRangeLow=NULL, double *mipValue=NULL, double *vertexRegionDeDxCutOff=NULL)
Calculates the total hadronic energy on the provided track based on the assumption of MIP...
float fMipValueInGev
the same as fMipValue but in GeV/cm units (FHiCL parameter)
T get(std::string const &key) const
virtual TVector3 Dir() const
Unit vector describing prong direction.
float fMipValue
Default MIP value in GeV/cm (FHiCL parameter)
float TrackEinMIP(const rb::Track &track)
A function for total energy calculation in terms of MIP.
std::map< unsigned int, double > DeDxInPlane(const rb::Track &muonTrack)
A function to calculate the energy, dE/dx and tracklength of a track through a plane.
rb::RecoHit RecoHit(const art::Ptr< rb::CellHit > &chit) const
Return calibrated hit based on assumed W coordinate.
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
std::vector< int > fPlanes
A vector of planes on muon track with hits in them.
std::vector< rb::WeightedHit > CleanUpTrack(const rb::Track &muonTrack, const rb::Cluster &slice, double *mipRangeHigh=NULL, double *mipRangeLow=NULL, double *mipValue=NULL, double *vertexRegionDeDxCutOff=NULL)
A function to disentangle the muon track from the hadronic energy close to the vertex. Returns a vector of rb::WeightedHits where the hit is the cellhits that were on the given track and weight is the fraction of the energy in the hit that came from the track. mipRangeLow, mipRangeHigh, mipValue and vertexRegionDeDxCutOff are pointers that point to the value of mip range limits, peak values andthe dedx cut off used to find the vertex region. If these pointers are not provided, the default values are used. The default values are good for muons, not for other particles.
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
float fExtraEOnTrack
FHiCL Parameters.
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
std::vector< rb::WeightedHit > CleanUpWeightedHits(const rb::Track &muonTrack, const rb::Cluster &slice, double *mipRangeHigh=NULL, double *mipRangeLow=NULL, double *mipValue=NULL, double *vertexRegionDeDxCutOff=NULL)
Returns a vector of rb::WeightedHit that have been impacted by the CleanUp procedure. hits are all the hits on the muonTrack within the vertex region and weight indicates the fraction of energy deposited in the cell that belonged to the track. mipRangeLow, mipRangeHigh, mipValue and vertexRegionDeDxCutOff are pointers that point to the value of mip range limits, peak values and the dedx cut off used to find the vertex region. If these pointers are not provided, the default values are used. The default values are good for muons, not for other particles.
std::vector< int > fPlanesX
A vector of X-view planes on muon track with hits in them.
std::vector< int > fPlanesY
A vector of Y-view planes on muon track with hits in them.
Simple hit+weight pair, returned from rb::Cluster::WeightedHits.
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
virtual void InterpolateXY(double z, double &x, double &y) const
void SortByDistFromTrack(std::vector< art::Ptr< rb::CellHit > > &hits, const rb::Track &track, int plane)
Sorts a vector of cellhits in ascending order by distance from the trajPt on track.
void SortByPlaneAndCell(std::vector< rb::WeightedHit > &whits)
std::map< unsigned int, double > fPlaneTrkLength
A map of plane and track length through plane.
float fMipRangeLowInGev
the same as fMipRangeLow but in GeV/cm units (FHiCL parameter)
void LastVertexPlanes(const rb::Track &muonTrack, int *vertexPlanes, double *vertexRegionDeDxCutOff=NULL)
A function to decide how large a vertex region should be considered for clean up. vertexRegion is an ...
float ExtraEOnTrackInGeV(const rb::Track &track, const rb::Cluster &slice, double *mipRangeHigh=NULL, double *mipRangeLow=NULL, double *mipValue=NULL, double *vertexRegionDeDxCutOff=NULL)
double Weight(unsigned int globalIdx) const
Weight assigned to the cell.
Encapsulate the geometry of one entire detector (near, far, ndos)
std::map< unsigned int, double > fPlaneEnergy
A map of plane and energy in plane left by the track.
float fMipRangeHighInGev
the same as fMipRangeHigh but in GeV/cm units (FHiCL parameter)
std::map< unsigned int, double > fAvgDeDx
A map of plane and average dedx over this and the next 2 planes.
float fMipRangeLow
Default lower bound of allowed MIP range (FHiCL parameter)