20 #include "NovaDAQConventions/DAQConventions.h" 34 #include <TSpectrum.h> 35 #include <TFitResult.h> 56 void HoughRhoTheta(
double x1,
double y1,
double x2,
double y2,
double* rho,
double*
theta,
double* sigmaRho,
double* sigmaTheta);
183 tree = tfs->
make<TTree>(
"Tree",
"Tree of info");
185 tree->Branch(
"run", &
run,
"run/I" );
188 tree->Branch(
"eventTime", &
evtTime,
"eventTime/I" );
189 tree->Branch(
"ndcm", &
ndcm,
"ndcm/I" );
210 std::map<int,rb::Cluster> clusterMap;
212 std::vector<art::Ptr<rb::CellHit > > cleanHitlist;
216 for (
unsigned int i = 0;
i < hitlist.size(); ++
i ) {
217 double minCellDistXY(1000.);
218 double minCellDistZ(1000.);
219 std::pair<double,double> hitPosI =
GetHitPos( hitlist[
i] );
220 double maxHitSepXY(0.);
223 for (
unsigned int j = 0;
j < hitlist.size(); ++
j ) {
224 if (
j == i )
continue;
225 double timeDiff =
fabs( hitlist[i]->TNS()*1
e-3 - hitlist[
j]->TNS()*1
e-3 );
228 if ( hitlist[
j]->
View() == hitlist[i]->
View() ) {
229 std::pair<double,double> hitPosJ =
GetHitPos( hitlist[
j] );
230 double cellDistXY =
fabs(hitPosI.first-hitPosJ.first);
231 if ( cellDistXY < minCellDistXY ) minCellDistXY = cellDistXY;
232 double cellDistZ =
fabs(hitPosI.second-hitPosJ.second);
233 if ( cellDistZ < minCellDistZ ) minCellDistZ = cellDistZ;
237 if ( minCellDistXY > maxHitSepXY || minCellDistZ >
fMaxHitSepZ ) {
238 clusterMap[0].Add(hitlist[i]);
241 cleanHitlist.push_back(hitlist[i]);
242 fTimeHist->Fill( hitlist[i]->TNS()*1
e-3 );
246 double *peaks =
fSpecAna->GetPositionX();
249 for (
unsigned int i = 0;
i < cleanHitlist.size(); ++
i ) {
252 for (
int p = 0;
p < nPeaks; ++
p ) {
259 if ( clusterMap.find(matchPeak) == clusterMap.end() ) clusterMap[matchPeak] =
rb::Cluster();
260 clusterMap[matchPeak].
Add(cleanHitlist[
i]);
264 std::vector<rb::Cluster> clusterList;
265 for (
auto itr = clusterMap.begin(); itr != clusterMap.end(); ++itr ) {
266 if ( itr->first == 0 ) itr->second.SetNoise(
true);
267 clusterList.push_back(itr->second);
276 double x2,
double y2,
277 double* rho,
double*
theta,
278 double* sigmarho,
double* sigmatheta)
281 *theta =
atan2(x1-x2,y2-y1);
282 *rho = 0.5*(
cos(*theta)*(x1+
x2)+
sin(*theta)*(y1+
y2));
286 while (*theta < -
M_PI/2.) {
290 while (*theta >=
M_PI/2.) {
297 *sigmarho = 3.0/
sqrt(12.);
298 *sigmatheta = M_SQRT2*(3.0/
sqrt(12.))/d;
306 for (
unsigned int i = 0;
i < hitList.size(); ++
i ) {
307 for (
unsigned int j =
i+1;
j < hitList.size(); ++
j ) {
308 if ( hitList[
j].fView != hitList[
i].fView )
continue;
309 double hitSep =
sqrt(
pow(hitList[
i].fZ-hitList[
j].fZ,2) +
pow(hitList[
i].fXY-hitList[
j].fXY,2) );
312 double rho,
theta, sigmaRho, sigmaTheta;
314 hitList[
j].fZ, hitList[
j].fXY,
316 &sigmaRho, &sigmaTheta );
317 mapHist->Fill(rho,theta);
327 TH1D *angleHist =
new TH1D(
"HoughAngleHist",
";Angle [rad];Hit Pairs",200,-
M_PI/2.,
M_PI/2.);
329 for (
unsigned int i = 0;
i < hitList.size(); ++
i ) {
330 for (
unsigned int j =
i+1;
j < hitList.size(); ++
j ) {
331 if ( hitList[
j].fView != hitList[
i].fView )
continue;
332 double hitSep =
sqrt(
pow(hitList[
i].fZ-hitList[
j].fZ,2) +
pow(hitList[
i].fXY-hitList[
j].fXY,2) );
335 double rho,
theta, sigmaRho, sigmaTheta;
337 hitList[
j].fZ, hitList[
j].fXY,
339 &sigmaRho, &sigmaTheta );
340 angleHist->Fill(theta);
350 TH1D *rhoHist =
new TH1D(
"HoughRhoHist",
";Rho [cm];Hit Pairs",geoSvc->
NPlanes()*3,-0.25*geoSvc->
DetLength(),1.25*geoSvc->
DetLength());
351 for (
unsigned int i = 0;
i < hitList.size(); ++
i ) {
352 for (
unsigned int j =
i+1;
j < hitList.size(); ++
j ) {
353 if ( hitList[
j].fView != hitList[
i].fView )
continue;
354 double hitSep =
sqrt(
pow(hitList[
i].fZ-hitList[
j].fZ,2) +
pow(hitList[
i].fXY-hitList[
j].fXY,2) );
357 double rho,
theta, sigmaRho, sigmaTheta;
359 hitList[
j].fZ, hitList[
j].fXY,
361 &sigmaRho, &sigmaTheta );
362 if (
fabs( theta - angle ) / angleSig > 2. )
continue;
396 double sigCellSepXY(0.);
401 for (
unsigned int i = 0;
i < slice.
NCell(); ++
i ) {
407 for (
unsigned int i = 0;
i+1 < slice.
NCell(); ++
i ) {
410 std::pair<double,double> hitPosI =
GetHitPos( hitI );
411 double XYI = hitPosI.first;
412 double ZI = hitPosI.second;
414 for (
unsigned int j =
i+1;
j < slice.
NCell(); ++
j ) {
416 if ( hitI->
View() != hitJ->
View() )
continue;
418 std::pair<double,double> hitPosJ =
GetHitPos( hitJ );
419 double XYJ = hitPosJ.first;
420 double ZJ = hitPosJ.second;
425 double score = deltaR;
439 unsigned int ndeaddcms = 0;
440 unsigned int ndcmhits[14][12] = {{0}};
444 const unsigned int db = cmap->
getDiBlock(chit.DaqChannel());
445 const unsigned int dcm = cmap->
getDCM(chit.DaqChannel());
447 ndcmhits[db-1][dcm-1]++;
450 for(
unsigned int i = 0;
i < 14; ++
i)
451 for(
unsigned int j = 0;
j < 12; ++
j)
452 if(ndcmhits[
i][
j]==0) ndeaddcms++;
454 return 168-ndeaddcms;
497 std::vector<art::Ptr<rb::CellHit > > hitlist;
498 for(
unsigned int i = 0;
i < hitcol->size(); ++
i){
500 hitlist.push_back(hit);
503 std::unique_ptr< std::vector<rb::Cluster> > clustercoltime(
new std::vector<rb::Cluster>);
504 std::unique_ptr< std::vector<rb::Cluster> > clustercoltracks(
new std::vector<rb::Cluster>);
511 clustercoltime->push_back(noiseCluster);
541 std::vector<rb::Cluster> timeSlicesY;
542 std::vector<rb::Cluster> timeSlicesX;
544 for (
unsigned int i = 0;
i < timeSlices.size(); ++
i ) {
545 if ( timeSlices[
i].IsNoise() ) {
546 for (
unsigned int j = 0;
j < timeSlices[
i].NCell(); ++
j ) noiseCluster.
Add(timeSlices[
i].AllCells()[
j]);
550 for (
unsigned int j = 0;
j < timeSlices[
i].NCell(); ++
j ) {
551 if ( timeSlices[
i].AllCells()[
j]->View() ==
geo::kY ) hitsY.
Add( timeSlices[
i].AllCells()[
j] );
552 if ( timeSlices[
i].AllCells()[
j]->View() ==
geo::kX ) hitsX.
Add( timeSlices[
i].AllCells()[
j] );
554 hitsY.
SetNoise( timeSlices[
i].IsNoise() );
555 hitsX.
SetNoise( timeSlices[
i].IsNoise() );
556 timeSlicesY.push_back(hitsY);
557 timeSlicesX.push_back(hitsX);
562 clustercoltime->push_back(noiseCluster);
563 for (
unsigned int i = 0;
i < timeSlices.size(); ++
i )
if ( !timeSlices[
i].IsNoise() && timeSlices[
i].NCell() > 0 ) clustercoltime->push_back(timeSlices[
i]);
568 for (
unsigned int i = 0;
i < timeSlices.size(); ++
i ) {
569 if ( timeSlices[
i].IsNoise() )
continue;
572 std::vector<rb::Cluster> trackSlicesY;
573 std::vector<rb::Cluster> trackSlicesX;
575 bool mapItrForward(
true);
581 for (
unsigned int j = 0;
j < trackSlicesY.size(); ++
j )
582 if ( trackSlicesY[
j].NCell() > 0 ) clustercoltracks->push_back(trackSlicesY[
j]);
583 for (
unsigned int j = 0; j < trackSlicesX.size(); ++
j )
584 if ( trackSlicesX[j].NCell() > 0 ) clustercoltracks->push_back(trackSlicesX[j]);
601 clustercoltracks->push_back(noiseCluster);
608 std::vector<rb::Cluster> trackSlices;
609 std::vector<rb::Cluster> finalTrackSlices;
610 if ( timeSlice.
IsNoise() )
return trackSlices;
616 double maxHitSepXY(0.);
630 tfs->
make<TH2D>(*mapHist);
637 int maxBin = angleHist->GetMaximumBin();
643 TF1 *agaus=
new TF1(
"agaus",
"gaus",angleHist->GetBinCenter(maxBin-5),angleHist->GetBinCenter(maxBin+5));
644 agaus->SetParameter(0,angleHist->GetBinContent(maxBin));
645 agaus->SetParameter(1,angleHist->GetBinCenter(maxBin));
646 agaus->SetParameter(2,.01);
647 agaus->SetParLimits(1,angleHist->GetBinCenter(maxBin-1),angleHist->GetBinCenter(maxBin+1));
648 TFitResultPtr fitResult = angleHist->Fit(
"agaus",
"QSNB",
"",angleHist->GetBinCenter(maxBin-10),angleHist->GetBinCenter(maxBin+10));
649 int fitStatus = fitResult;
651 if ( fitStatus != 0 ){
657 double maxAngle = fitResult->Parameter(1);
658 double sigAngle = fitResult->Parameter(2);
659 sigAngle=
sqrt(sigAngle*sigAngle);
664 TH1D *rhoHist =
GetHoughRhos( thisHitList, maxAngle, sigAngle );
683 std::map< int, rb::Cluster > hitsByPlane;
684 for (
unsigned int hitItr = 0; hitItr < timeSlice.
NCell(); ++hitItr ) {
685 if ( timeSlice.
AllCells()[hitItr]->View() !=
view )
continue;
687 if ( hitsByPlane.find(plane) == hitsByPlane.end() ) hitsByPlane[plane] =
rb::Cluster(view);
692 std::map< int, std::vector<rb::Cluster> > inPlaneClusters;
693 for (
auto planeItr = hitsByPlane.begin(); planeItr != hitsByPlane.end(); ++planeItr ) {
694 int thePlane = planeItr->first;
696 int CurrentClusterID = 1;
698 for (
unsigned int j = 0;
j < planeItr->second.NCell(); ++
j ) {
705 inPlaneClusters[thePlane].resize(CurrentClusterID,
rb::Cluster(view));
706 for (
unsigned int j = 0;
j < planeItr->second.NCell(); ++
j ) {
708 noiseCluster.
Add( planeItr->second.AllCells()[
j] );
711 inPlaneClusters[thePlane][
fClusterIDs[
j]].Add( planeItr->second.AllCells()[
j] );
716 if ( view ==
geo::kY ) mapItrForward = (maxAngle > 0.);
718 auto mapItrStart = inPlaneClusters.begin();
719 auto mapItrEnd = inPlaneClusters.end();
720 if ( !mapItrForward ) {
721 mapItrStart = --(inPlaneClusters.rbegin().base());
722 mapItrEnd = --(inPlaneClusters.rend().base());
726 for (
unsigned int clusItr = 0; clusItr < mapItrStart->second.size(); ++clusItr ) {
727 trackSlices.push_back( mapItrStart->second[clusItr] );
731 auto mapItr = mapItrStart;
732 if ( mapItrForward ) ++mapItr;
734 while ( mapItr != mapItrEnd ) {
735 if ( mapItr->second.size() == 0 ) {
736 if ( mapItrForward ) ++mapItr;
740 int thisPlane = mapItr->first;
743 std::vector<int> bestMatchInt(inPlaneClusters[thisPlane].
size(),-1);
744 std::vector<double> bestMatchVal(inPlaneClusters[thisPlane].
size(),100000.);
745 for (
unsigned int j = 1;
j < mapItr->second.size(); ++
j ) {
746 if ( mapItr->second[
j].NCell() == 0 )
continue;
750 for (
unsigned int c = 0;
c < mapItr->second[
j].NCell(); ++
c ) {
752 std::pair<double,double> hitPos =
GetHitPos( hit );
753 meanXY += hitPos.first;
754 meanZ += hitPos.second;
756 meanXY = meanXY / double(mapItr->second[
j].NCell());
757 meanZ = meanZ / double(mapItr->second[
j].NCell());
760 for (
unsigned int i = 0;
i < trackSlices.size(); ++
i ) {
762 std::map< int, std::pair<double,double> > trackSliceClusterPositions;
763 std::map< int, int > trackSliceClusterCounts;
764 double meanXYC(0.), meanZC(0.);
765 int nCells(0),
maxPlane(0), minPlane(9999);
766 for (
unsigned int c = 0;
c < trackSlices[
i].NCell(); ++
c ) {
767 if (
abs(
int(trackSlices[
i].Cell(
c)->Plane()) -
int(thisPlane) ) > 10 )
continue;
769 if ( trackSlices[
i].Cell(
c)->Plane() < minPlane ) minPlane = trackSlices[
i].Cell(
c)->Plane();
770 std::pair<double,double> hitPosC =
GetHitPos( trackSlices[
i].Cell(
c) );
771 trackSliceClusterPositions[
int(trackSlices[
i].Cell(
c)->Plane())].first += hitPosC.first;
772 trackSliceClusterPositions[
int(trackSlices[
i].Cell(
c)->Plane())].second += hitPosC.second;
773 trackSliceClusterCounts[
int(trackSlices[
i].Cell(
c)->Plane())]++;
774 meanXYC += hitPosC.first;
775 meanZC += hitPosC.second;
778 auto posItr = trackSliceClusterPositions.begin();
779 auto countItr = trackSliceClusterCounts.begin();
780 for (
unsigned int c = 0;
c < trackSliceClusterPositions.size(); ++
c ) {
781 posItr->second.first = posItr->second.first / countItr->second;
782 posItr->second.second = posItr->second.second / countItr->second;
789 meanXYC = meanXYC / double( nCells );
790 meanZC = meanZC / double( nCells );
792 double closestClusterDistXY(1000000.);
793 double closestClusterDistZ(1000000.);
794 double meanHoughAngle(0.), meanHoughAngleSig(0.);
795 double meanHoughRho(0.), meanHoughRhoSig(0.);
797 for (
auto c = trackSliceClusterPositions.begin();
c != trackSliceClusterPositions.end(); ++
c ) {
798 double deltaXY =
fabs(
c->second.first-meanXY);
799 double deltaZ =
fabs(
c->second.second-meanZ);
800 if ( deltaXY < closestClusterDistXY ) closestClusterDistXY = deltaXY;
801 if ( deltaZ < closestClusterDistZ ) closestClusterDistZ = deltaZ;
802 double houghRho, houghAngle, houghRhoSig, houghAngleSig;
805 &houghRho, &houghAngle, &houghRhoSig, &houghAngleSig );
806 meanHoughAngle += houghAngle;
807 meanHoughRho += houghRho;
808 meanHoughAngleSig += houghAngleSig*houghAngleSig;
809 meanHoughRhoSig += houghAngleSig*houghRhoSig;
812 if ( closestClusterDistXY > maxHitSepXY || closestClusterDistZ >
fMaxHitSepZ ) {
815 meanHoughAngle = meanHoughAngle / double(nMeanEntries);
816 meanHoughRho = meanHoughRho / double(nMeanEntries);
817 meanHoughAngleSig =
sqrt(meanHoughAngleSig) / double(nMeanEntries);
818 meanHoughRhoSig =
sqrt(meanHoughRhoSig) / double(nMeanEntries);
823 double houghRho, houghAngle, houghRhoSig, houghAngleSig;
824 HoughRhoTheta( meanZ, meanXY, meanZC, meanXYC, &houghRho, &houghAngle, &houghRhoSig, &houghAngleSig );
832 if ( matchVal < bestMatchVal[
j] ) {
834 bestMatchVal[
j] = matchVal;
838 for (
unsigned int j = 1;
j < mapItr->second.size(); ++
j ) {
840 if ( bestMatchInt[
j] < 0 ) {
841 trackSlices.push_back( mapItr->second[
j] );
845 for (
unsigned int k = 0; k < mapItr->second[
j].NCell(); ++k ) {
846 trackSlices[bestMatchInt[
j]].Add( mapItr->second[
j].AllCells()[k] );
849 if ( mapItrForward ) ++mapItr;
853 for (
unsigned int i = 0;
i < trackSlices.size(); ++
i ) {
855 for (
unsigned int j = 0;
j < trackSlices[
i].NCell(); ++
j ) {
856 noiseCluster.
Add( trackSlices[
i].AllCells()[
j] );
859 finalTrackSlices.push_back( trackSlices[
i] );
862 if ( finalTrackSlices.size() <
fMinTracks ) finalTrackSlices.clear();
865 return finalTrackSlices;
885 for(
unsigned int seedi = 0; seedi < seeds.size(); ++seedi)
887 seeds.erase(seeds.begin());
891 while(seeds.size() > 0) {
898 for(
unsigned int resulti = 0; resulti < result.size(); ++resulti) {
906 seeds.push_back(result[resulti]);
913 seeds.erase(seeds.begin());
double fMaxHoughMatchVal
maximum match value for Hough extrapolation
SubRunNumber_t subRun() const
double fMinHoughHitSep
minimum separation of hit pairs for Hough map
uint32_t fTriggerMask_Prescale
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
double fTimePeakThresh
threshold for temporal slice peak search
void GetCenter(double *xyz, double localz=0.0) const
int fTimeRangeHigh
high edge of time range for temporal slicing (integer microseconds)
unsigned int fMinPts
minimum number of neighbors for a hit to be in a cluster
fvar< T > fabs(const fvar< T > &x)
Float_t y1[n_points_granero]
void SetNoise(bool noise)
Declare the cluster to consist of noise hits or not.
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
unsigned short Plane() const
Float_t x1[n_points_granero]
double fSigCellSepX
scale for X separation between cells in X-view to include in non-noise slice
const CellGeo * Cell(int icell) const
constexpr std::uint32_t timeHigh() const
Vertical planes which measure X.
std::pair< double, double > GetHitPos(art::Ptr< rb::CellHit > hit)
A collection of associated CellHits.
double fTimePeakSig
significance for temporal slice peak search
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
void SortByTime(std::vector< art::Ptr< rb::CellHit > > &c)
Sort c in time order (earliest to latest).
const PlaneGeo * Plane(unsigned int i) const
double fMinTrackLen
minimum track length for good tracks
DEFINE_ART_MODULE(TestTMapFile)
std::vector< rb::Cluster > MakeTrackSlices(rb::Cluster timeSlice, rb::Cluster &noiseCluster, geo::View_t view, bool &mapItrForward)
Horizontal planes which measure Y.
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
double fSigCellSepZ
scale for Z separation between cells to include in non-noise slice
int fTimeRangeLow
low edge of time range for temporal slicing (integer microseconds)
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
ProductID put(std::unique_ptr< PROD > &&product)
unsigned short Cell() const
std::string fRawDataLabel
Data label for trigger prescale.
View_t View() const
Which coordinate does this plane measure.
Just the essential information required for track fitting.
double fRhoPeakSig
significance for rho peak search
TH1D * GetHoughAngles(rb::HitList hitList)
double fTimePeakCut
cut for temporal slicing
Geometry information for a single readout plane.
void HoughRhoTheta(double x1, double y1, double x2, double y2, double *rho, double *theta, double *sigmaRho, double *sigmaTheta)
double fMinHoughPairs
minimum number of hit pairs in maximum Hough angle or rho histogram to proceed
EventNumber_t event() const
std::vector< int > fClusterIDs
std::string fCellHitInput
Read CellHits from this input.
double fChi2NdofCut
cut on Chi2/Ndof for track fit
double fHoughExtrapPlanes
number of planes required for Hough extrapolation
std::vector< std::vector< unsigned int > > fNeighborList
std::string fTrackSliceInstance
instance name of track slices
int CountActiveDCM(const std::vector< rb::CellHit > &hitlist)
double fHoughRhoPeakThresh
threshold for Hough map rho peak search
bool fSaveHists
save intermediate histograms
std::vector< rb::Cluster > TemporalClusters(std::vector< art::Ptr< rb::CellHit > > hitList)
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
A rawdata::RawDigit with channel information decoded.
int fMaxSlices
maximum number of time slices to create
T * make(ARGS...args) const
double DetHalfWidth() const
Data resulting from a Hough transform on the cell hit positions.
void produce(art::Event &evt)
double fMaxNeighborScore
neighbor score cutoff
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
std::string fTimeSliceInstance
instance name of time slices
unsigned int fMinTrackHits
minimum number of hits for a track slice
bool GenerateCluster(unsigned int point, int ClID)
unsigned int fMinTracks
minimum number of tracks in each time slice
double fMaxHitSepX
maximum X separation between hits along track slice in X-view
void FindNeighbors(rb::Cluster clusterList, geo::View_t view)
double fSigCellSepY
scale for Y separation between cells in Y-view to include in non-noise slice
double fRhoPeakThresh
threshold for rho peak search
cmap::CMap class source code
double fMaxHitSepZ
maximum Z separation between hits along track slice
bool fWriteTimeSlices
save time slices to output file
dcm_id_t getDCM(dchan daqchan) const
Decode the dcm ID from a dchan.
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
unsigned int NPlanes() const
T min(const caf::Proxy< T > &a, T b)
bool IsNoise() const
Is the noise flag set?
double fMaxHitSepY
maximum Y separation between hits along track slice in Y-view
Encapsulate the cell geometry.
AirSlicer(fhicl::ParameterSet const &pset)
int fMaxPeaks
maximum number of rho peaks to count
TH1D * GetHoughRhos(rb::HitList hitList, double angle, double angleSig)
diblock_t getDiBlock(dchan daqchan) const
Decode the diblock ID from a dchan.
TH2D * GetHoughMap(rb::HitList hitList)
Encapsulate the geometry of one entire detector (near, far, ndos)
int fMaxSliceHits
maximum number of hits to perform slicing