AirSlicer_module.cc
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////
2 // \file AirSlicer_module.cc
3 // \brief Modified version of slicer4D for reconstructing high energy
4 // cosmic ray air showers.
5 // \version $Id: AirSlicer_module.cc,v 1.0 2015-04-23 09:30:00 dwwhitti Exp $
6 // \authors grohmc@fnal.gov, dwwhitti@indiana.edu
7 ///////////////////////////////////////////////////////////////////////////
8 
9 #include "art/Framework/Core/EDProducer.h"
10 #include "art/Framework/Core/ModuleMacros.h"
11 #include "art/Framework/Principal/Event.h"
12 #include "art/Framework/Services/Optional/TFileDirectory.h"
13 #include "art/Framework/Services/Optional/TFileService.h"
14 #include "art/Framework/Services/Registry/ServiceHandle.h"
15 
16 #include "Calibrator/Calibrator.h"
17 #include "Geometry/Geometry.h"
18 #include "GeometryObjects/PlaneGeo.h"
19 #include "GeometryObjects/CellGeo.h"
20 #include "NovaDAQConventions/DAQConventions.h"
21 #include "RecoBase/CellHit.h"
22 #include "RecoBase/Cluster.h"
23 #include "RecoBase/HitList.h"
24 #include "RecoBase/HoughResult.h"
25 
26 #include "RawData/RawTrigger.h"
27 
28 #include <cmath>
29 #include <string>
30 #include <vector>
31 
32 #include <TH1D.h>
33 #include <TH2D.h>
34 #include <TSpectrum.h>
35 #include <TFitResult.h>
36 #include <TFile.h>
37 #include "TTree.h"
38 #include "TGraph.h"
39 
40 // AirSlicer header
41 namespace airshower {
42 
43  class AirSlicer : public art::EDProducer {
44  public:
45  explicit AirSlicer(fhicl::ParameterSet const &pset);
46  virtual ~AirSlicer();
47 
48  void produce(art::Event& evt);
49  void beginJob();
50 
51  protected:
52  // internal methods
53  bool GenerateCluster(unsigned int point, int ClID);
54  int CountActiveDCM(const std::vector<rb::CellHit>& hitlist);
55  std::vector<rb::Cluster> TemporalClusters( std::vector<art::Ptr<rb::CellHit > > hitList );
56  void HoughRhoTheta(double x1, double y1, double x2, double y2, double* rho, double* theta, double* sigmaRho, double* sigmaTheta);
57  TH2D* GetHoughMap(rb::HitList hitList);
58  TH1D* GetHoughAngles(rb::HitList hitList);
59  TH1D* GetHoughRhos(rb::HitList hitList, double angle, double angleSig);
60  std::pair<double,double> GetHitPos( art::Ptr<rb::CellHit> hit );
61  void FindNeighbors(rb::Cluster clusterList, geo::View_t view);
62  std::vector<rb::Cluster> MakeTrackSlices( rb::Cluster timeSlice, rb::Cluster& noiseCluster, geo::View_t view, bool& mapItrForward );
63 
64  // configuration settings
65  std::string fCellHitInput; ///< Read CellHits from this input
66  unsigned int fMinPts; ///< minimum number of neighbors for a hit to be in a cluster
67  int fMaxSlices; ///< maximum number of time slices to create
68  int fMaxPeaks; ///< maximum number of rho peaks to count
69  bool fWriteTimeSlices; ///< save time slices to output file
70  std::string fTimeSliceInstance; ///< instance name of time slices
71  std::string fTrackSliceInstance; ///< instance name of track slices
72  int fTimeRangeLow; ///< low edge of time range for temporal slicing (integer microseconds)
73  int fTimeRangeHigh; ///< high edge of time range for temporal slicing (integer microseconds)
74  double fTimePeakCut; ///< cut for temporal slicing
75  double fTimePeakSig; ///< significance for temporal slice peak search
76  double fTimePeakThresh; ///< threshold for temporal slice peak search
77  double fRhoPeakSig; ///< significance for rho peak search
78  double fRhoPeakThresh; ///< threshold for rho peak search
79  double fSigCellSepX; ///< scale for X separation between cells in X-view to include in non-noise slice
80  double fSigCellSepY; ///< scale for Y separation between cells in Y-view to include in non-noise slice
81  double fSigCellSepZ; ///< scale for Z separation between cells to include in non-noise slice
82  double fMaxNeighborScore; ///< neighbor score cutoff
83  double fMaxHitSepX; ///< maximum X separation between hits along track slice in X-view
84  double fMaxHitSepY; ///< maximum Y separation between hits along track slice in Y-view
85  double fMaxHitSepZ; ///< maximum Z separation between hits along track slice
86  double fMinHoughPairs; ///< minimum number of hit pairs in maximum Hough angle or rho histogram to proceed
87  double fMinHoughHitSep; ///< minimum separation of hit pairs for Hough map
88  double fHoughRhoPeakThresh; ///< threshold for Hough map rho peak search
89  double fHoughExtrapPlanes; ///< number of planes required for Hough extrapolation
90  double fMaxHoughMatchVal; ///< maximum match value for Hough extrapolation
91  unsigned int fMinTrackHits; ///< minimum number of hits for a track slice
92  unsigned int fMinTracks; ///< minimum number of tracks in each time slice
93  double fChi2NdofCut; ///< cut on Chi2/Ndof for track fit
94  double fMinTrackLen; ///< minimum track length for good tracks
95  int fMaxSliceHits; ///< maximum number of hits to perform slicing
96  bool fSaveHists; ///< save intermediate histograms
97  std::string fRawDataLabel; ///< Data label for trigger prescale
98 
99  // internal objects
100  std::vector< std::vector<unsigned int> > fNeighborList;
101  std::vector< int > fClusterIDs;
102 
103  TSpectrum *fSpecAna;
104  TSpectrum *fRhoSpecAna;
105 
106  // Global Variables for creating a ROOT tree
107  TTree *tree;
108  int run; // Event run number
109  int subrun; // Event subrun number
110  int evtNum; // Event number
111  int evtTime; // Event time in unix time (s)
112  int prescale; // Trigger prescale
113  int ndcm; // Number of active dcms in the event
114  int houghMultY; // Y-View multiplicity from hough transform
115  int trackMultY; // Y-View multiplicity from slicing
116  int houghMultX; // X-View multiplicity from hough transform
117  int trackMultX; // X-View multiplicity from slicing
118  double houghAngleY; // Y-View peak hough angle
119  double houghStDevY; // Y-View stdev in peak hough angle
120  double houghAngleX; // X-View peak hough angle
121  double houghStDevX; // X-View stdev in peak hough angle
122  int totalHits; // total hits in time slice
123  int totalADC; // total ADC in time slice
124  };
125 
126 }
127 
128 namespace airshower
129 {
130  //----------------------------------------------------------------------
131  AirSlicer::AirSlicer(fhicl::ParameterSet const &pset) :
132  fCellHitInput (pset.get<std::string >("CellHitInput") ),
133  fMinPts (pset.get<unsigned int>("MinPts") ),
134  fMaxSlices (pset.get<int >("MaxSlices") ),
135  fMaxPeaks (pset.get<int >("MaxPeaks") ),
136  fWriteTimeSlices (pset.get<bool >("WriteTimeSlices") ),
137  fTimeSliceInstance (pset.get<std::string >("TimeSliceInstance") ),
138  fTrackSliceInstance(pset.get<std::string >("TrackSliceInstance")),
139  fTimeRangeLow (pset.get<int >("TimeRangeLow") ),
140  fTimeRangeHigh (pset.get<int >("TimeRangeHigh") ),
141  fTimePeakCut (pset.get<double >("TimePeakCut") ),
142  fTimePeakSig (pset.get<double >("TimePeakSig") ),
143  fTimePeakThresh (pset.get<double >("TimePeakThresh") ),
144  fRhoPeakSig (pset.get<double >("RhoPeakSig") ),
145  fRhoPeakThresh (pset.get<double >("RhoPeakThresh") ),
146  fSigCellSepX (pset.get<double >("SigCellSepX") ),
147  fSigCellSepY (pset.get<double >("SigCellSepY") ),
148  fSigCellSepZ (pset.get<double >("SigCellSepZ") ),
149  fMaxNeighborScore (pset.get<double >("MaxNeighborScore") ),
150  fMaxHitSepX (pset.get<double >("MaxHitSepX") ),
151  fMaxHitSepY (pset.get<double >("MaxHitSepY") ),
152  fMaxHitSepZ (pset.get<double >("MaxHitSepZ") ),
153  fMinHoughPairs (pset.get<int >("MinHoughPairs") ),
154  fMinHoughHitSep (pset.get<double >("MinHoughHitSep") ),
155  fHoughExtrapPlanes (pset.get<int >("HoughExtrapPlanes") ),
156  fMaxHoughMatchVal (pset.get<double >("MaxHoughMatchVal") ),
157  fMinTrackHits (pset.get<unsigned int>("MinTrackHits") ),
158  fMinTracks (pset.get<unsigned int>("MinTracks") ),
159  fChi2NdofCut (pset.get<double >("Chi2NdofCut") ),
160  fMinTrackLen (pset.get<double >("MinTrackLen") ),
161  fMaxSliceHits (pset.get<int >("MaxSliceHits") ),
162  fSaveHists (pset.get<bool >("SaveHists") ),
163  fRawDataLabel (pset.get<std::string >("RawDataLabel") )
164  {
165  produces< std::vector<rb::Cluster> >(fTimeSliceInstance);
166  produces< std::vector<rb::Cluster> >(fTrackSliceInstance);
167 
168  //Spectrums for peak counting
169  fSpecAna = new TSpectrum(fMaxSlices,fTimePeakSig);
170  fRhoSpecAna = new TSpectrum(fMaxPeaks,fRhoPeakSig);
171  }
172 
173  //----------------------------------------------------------------------
174  AirSlicer::~AirSlicer()
175  { }
176 
177  //----------------------------------------------------------------------
178  void AirSlicer::beginJob()
179  {
180  art::ServiceHandle<art::TFileService> tfs;
181 
182  // Create the tree branches
183  tree = tfs->make<TTree>("Tree","Tree of info");
184 
185  tree->Branch("run", &run, "run/I" );
186  tree->Branch("subrun", &subrun, "subrun/I" );
187  tree->Branch("evtNum", &evtNum, "evtNum/I" );
188  tree->Branch("eventTime", &evtTime, "eventTime/I" );
189  tree->Branch("ndcm", &ndcm, "ndcm/I" );
190  tree->Branch("prescale", &prescale, "prescale/I" );
191 
192  tree->Branch("houghMultY", &houghMultY, "houghMultY/I" );
193  tree->Branch("trackMultY", &trackMultY, "trackMultY/I" );
194 
195  tree->Branch("houghMultX", &houghMultX, "houghMultX/I" );
196  tree->Branch("trackMultX", &trackMultX, "trackMultX/I" );
197 
198  tree->Branch("houghAngleY",&houghAngleY,"houghAngleY/D");
199  tree->Branch("houghStDevY",&houghStDevY,"houghStDevY/D");
200  tree->Branch("houghAngleX",&houghAngleX,"houghAngleX/D");
201  tree->Branch("houghStDevX",&houghStDevX,"houghStDevX/D");
202 
203  tree->Branch("totalHits", &totalHits, "totalHits/I" );
204  tree->Branch("totalADC" , &totalADC , "totalADC/I" );
205  }
206 
207  //----------------------------------------------------------------------
208  std::vector<rb::Cluster> AirSlicer::TemporalClusters(std::vector<art::Ptr<rb::CellHit > > hitlist)
209  {
210  std::map<int,rb::Cluster> clusterMap;
211  clusterMap[0] = rb::Cluster(); // noise hits (unassociated to peaks)
212  std::vector<art::Ptr<rb::CellHit > > cleanHitlist;
213 
214  // Microsecond-level binning
215  TH1D *fTimeHist = new TH1D("fTimeHist",";time [#mus];hits in view",fTimeRangeHigh-fTimeRangeLow,fTimeRangeLow,fTimeRangeHigh);
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.);
221  if ( hitlist[i]->View() == geo::kY ) maxHitSepXY = fMaxHitSepY;
222  if ( hitlist[i]->View() == geo::kX ) maxHitSepXY = fMaxHitSepX;
223  for ( unsigned int j = 0; j < hitlist.size(); ++j ) {
224  if ( j == i ) continue;
225  double timeDiff = fabs( hitlist[i]->TNS()*1e-3 - hitlist[j]->TNS()*1e-3 );
226  if ( timeDiff > fTimePeakCut ) continue;
227  // Check whether it's isolated in it's own view
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;
234  }
235  }
236 
237  if ( minCellDistXY > maxHitSepXY || minCellDistZ > fMaxHitSepZ ) {
238  clusterMap[0].Add(hitlist[i]);
239  continue;
240  }
241  cleanHitlist.push_back(hitlist[i]);
242  fTimeHist->Fill( hitlist[i]->TNS()*1e-3 );
243  }
244 
245  int nPeaks = fSpecAna->Search( fTimeHist, fTimePeakSig, "goff nodraw", fTimePeakThresh );
246  double *peaks = fSpecAna->GetPositionX();
247 
248  // Loop through hits and associate to clusters
249  for ( unsigned int i = 0; i < cleanHitlist.size(); ++i ) {
250  int matchPeak(0);
251  // Find matching time peak
252  for ( int p = 0; p < nPeaks; ++p ) {
253  if ( fabs( cleanHitlist[i]->TNS()*1e-3 - peaks[p] ) < fTimePeakCut ) {
254  matchPeak = p+1;
255  break;
256  }
257  }
258  // Add in this hit
259  if ( clusterMap.find(matchPeak) == clusterMap.end() ) clusterMap[matchPeak] = rb::Cluster();
260  clusterMap[matchPeak].Add(cleanHitlist[i]);
261  }
262 
263  // Transfer to output cluster vector
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);
268  }
269 
270  delete fTimeHist;
271  return clusterList;
272  }
273 
274  //----------------------------------------------------------------------
275  void AirSlicer::HoughRhoTheta(double x1, double y1,
276  double x2, double y2,
277  double* rho, double* theta,
278  double* sigmarho, double* sigmatheta)
279  {
280  // Modified from MultiHough/MultiHough2P.cxx
281  *theta = atan2(x1-x2,y2-y1); // angle w.r.t. vertical
282  *rho = 0.5*(cos(*theta)*(x1+x2)+sin(*theta)*(y1+y2));
283  //
284  // Keep rho and theta within range [-pi/2,pi/2) (-R,+R)
285  //
286  while (*theta < -M_PI/2.) {
287  *theta += M_PI;
288  *rho = -(*rho);
289  }
290  while (*theta >= M_PI/2.) {
291  *theta -= M_PI;
292  *rho = -(*rho);
293  }
294 
295  // Compute the error estimates
296  double d = sqrt(pow(x2-x1,2)+pow(y2-y1,2));
297  *sigmarho = 3.0/sqrt(12.);
298  *sigmatheta = M_SQRT2*(3.0/sqrt(12.))/d;
299  }
300 
301  //----------------------------------------------------------------------
302  TH2D* AirSlicer::GetHoughMap(rb::HitList hitList)
303  {
304  art::ServiceHandle<geo::Geometry> geoSvc;
305  TH2D *mapHist = new TH2D("HoughMap",";Rho [cm];Angle [rad]",geoSvc->NPlanes()*4,-1.5*geoSvc->DetLength(),2.5*geoSvc->DetLength(),100,-M_PI/2.,M_PI/2.);
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) );
310  if ( hitSep > geoSvc->DetHalfWidth() ) continue;
311  if ( hitSep < fMinHoughHitSep ) continue;
312  double rho, theta, sigmaRho, sigmaTheta;
313  HoughRhoTheta( hitList[i].fZ, hitList[i].fXY,
314  hitList[j].fZ, hitList[j].fXY,
315  &rho, &theta,
316  &sigmaRho, &sigmaTheta );
317  mapHist->Fill(rho,theta);
318  }
319  }
320  return mapHist;
321  }
322 
323  //----------------------------------------------------------------------
324  TH1D* AirSlicer::GetHoughAngles(rb::HitList hitList)
325  {
326  art::ServiceHandle<geo::Geometry> geoSvc;
327  TH1D *angleHist = new TH1D("HoughAngleHist",";Angle [rad];Hit Pairs",200,-M_PI/2.,M_PI/2.);
328 
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) );
333  if ( hitSep > geoSvc->DetHalfWidth() ) continue;
334  if ( hitSep < fMinHoughHitSep ) continue;
335  double rho, theta, sigmaRho, sigmaTheta;
336  HoughRhoTheta( hitList[i].fZ, hitList[i].fXY,
337  hitList[j].fZ, hitList[j].fXY,
338  &rho, &theta,
339  &sigmaRho, &sigmaTheta );
340  angleHist->Fill(theta);
341  }
342  }
343  return angleHist;
344  }
345 
346  //----------------------------------------------------------------------
347  TH1D* AirSlicer::GetHoughRhos(rb::HitList hitList, double angle, double angleSig)
348  {
349  art::ServiceHandle<geo::Geometry> geoSvc;
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) );
355  if ( hitSep > geoSvc->DetHalfWidth() ) continue;
356  if ( hitSep < fMinHoughHitSep ) continue;
357  double rho, theta, sigmaRho, sigmaTheta;
358  HoughRhoTheta( hitList[i].fZ, hitList[i].fXY,
359  hitList[j].fZ, hitList[j].fXY,
360  &rho, &theta,
361  &sigmaRho, &sigmaTheta );
362  if ( fabs( theta - angle ) / angleSig > 2. ) continue;
363  rhoHist->Fill(rho);
364  }
365  }
366  return rhoHist;
367  }
368 
369  //----------------------------------------------------------------------
370  std::pair<double,double> AirSlicer::GetHitPos( art::Ptr<rb::CellHit> hit )
371  {
372  art::ServiceHandle<geo::Geometry> geoSvc;
373 
374  unsigned int plane = hit->Plane();
375  unsigned int cell = hit->Cell();
376  const geo::PlaneGeo* thePlane = geoSvc->Plane(plane);
377  const geo::CellGeo* theCell = thePlane->Cell(cell);
378 
379  double xyz[3];
380  theCell->GetCenter(xyz);
381  double xy(0.);
382  if ( thePlane->View() == geo::kX ) xy = xyz[0];
383  if ( thePlane->View() == geo::kY ) xy = xyz[1];
384  double z = xyz[2];
385 
386  return std::make_pair(xy,z);
387  }
388 
389  //----------------------------------------------------------------------
390  void AirSlicer::FindNeighbors(rb::Cluster slice, geo::View_t view)
391  {
392  fClusterIDs.clear();
393  fNeighborList.clear();
394  fNeighborList.resize( slice.NCell() );
395 
396  double sigCellSepXY(0.);
397  if ( view == geo::kY ) sigCellSepXY = fSigCellSepY;
398  if ( view == geo::kX ) sigCellSepXY = fSigCellSepX;
399 
400  // First add each point to its own list (the DBSCAN algorithm needs this to work correctly.)
401  for ( unsigned int i = 0; i < slice.NCell(); ++i ) {
402  fClusterIDs.push_back(-1);
403  fNeighborList[i].push_back(i);
404  }
405 
406  // Use i+1 < NCell to protect against the case where NCell = 0
407  for ( unsigned int i = 0; i+1 < slice.NCell(); ++i ) {
408  art::Ptr<rb::CellHit> hitI = slice.Cell(i);
409 
410  std::pair<double,double> hitPosI = GetHitPos( hitI );
411  double XYI = hitPosI.first;
412  double ZI = hitPosI.second;
413 
414  for ( unsigned int j = i+1; j < slice.NCell(); ++j ) {
415  art::Ptr<rb::CellHit> hitJ = slice.Cell(j);
416  if ( hitI->View() != hitJ->View() ) continue;
417 
418  std::pair<double,double> hitPosJ = GetHitPos( hitJ );
419  double XYJ = hitPosJ.first;
420  double ZJ = hitPosJ.second;
421 
422  // use physical distance for neighbor score
423  double deltaR = sqrt( pow(XYI-XYJ,2)/pow(sigCellSepXY,2) + pow(ZI-ZJ,2)/pow(fSigCellSepZ,2) );
424 
425  double score = deltaR;
426 
427  if ( score <= fMaxNeighborScore ) {
428  // i and j are neighbors, so they go into each others' lists
429  fNeighborList[i].push_back(j);
430  fNeighborList[j].push_back(i);
431  }
432  }
433  }
434  }
435 
436  //----------------------------------------------------------------------
437  int AirSlicer::CountActiveDCM(const std::vector<rb::CellHit>& hitlist)
438  {
439  unsigned int ndeaddcms = 0;
440  unsigned int ndcmhits[14][12] = {{0}};
441  const daqchannelmap::DAQChannelMap* cmap = art::ServiceHandle<cmap::CMap>()->Map();
442 
443  for(const rb::CellHit& chit: hitlist){
444  const unsigned int db = cmap->getDiBlock(chit.DaqChannel());
445  const unsigned int dcm = cmap->getDCM(chit.DaqChannel());
446 
447  ndcmhits[db-1][dcm-1]++;
448  }
449 
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++;
453 
454  return 168-ndeaddcms;
455  }
456 
457 
458 
459  //----------------------------------------------------------------------
460  void AirSlicer::produce(art::Event& evt)
461  {
462  // Reset tree
463  run = evt.run();
464  subrun = evt.subRun();
465  evtNum = evt.event();
466  evtTime = evt.time().timeHigh();
467  prescale = -999;
468  houghMultY = -999;
469  trackMultY = -999;
470  houghMultX = -999;
471  trackMultX = -999;
472  houghAngleY = -999;
473  houghStDevY = -999;
474  houghAngleX = -999;
475  houghStDevX = -999;
476  totalHits = -999;
477  totalADC = -999;
478 
479  // Modified from DDTPrescaleOffline Package
480  // use rawdata::RawTrigger method to get prescale number
481  art::Handle< std::vector<rawdata::RawTrigger> > rawtrigger;
482  evt.getByLabel(fRawDataLabel, rawtrigger);
483  const rawdata::RawTrigger& pre = (*rawtrigger)[0];
484  prescale = pre.fTriggerMask_Prescale;
485 
486  art::ServiceHandle<calib::Calibrator> cal;
487  art::ServiceHandle<art::TFileService> tfs;
488  art::ServiceHandle<geo::Geometry> geoSvc;
489 
490  // Get the cell hits
491  art::Handle< std::vector<rb::CellHit> > hitcol;
492  evt.getByLabel(fCellHitInput, hitcol);
493 
494  ndcm = CountActiveDCM(*hitcol);
495 
496  // Load the cell hits into a vector for easy use
497  std::vector<art::Ptr<rb::CellHit > > hitlist;
498  for(unsigned int i = 0; i < hitcol->size(); ++i){
499  art::Ptr<rb::CellHit> hit(hitcol, i);
500  hitlist.push_back(hit);
501  }
502 
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>);
505  rb::Cluster noiseCluster;
506  noiseCluster.SetNoise(true);
507 
508  // Avoid slicing on extremely large events
509  if((int)hitlist.size() > fMaxSliceHits){
510  if ( fWriteTimeSlices ) {
511  clustercoltime->push_back(noiseCluster);
512  evt.put(std::move(clustercoltime),fTimeSliceInstance);
513  }
514 
515  trackMultY=1000;
516  houghMultY=1000;
517  trackMultX=1000;
518  houghMultX=1000;
519 
520  tree->Fill();
521  houghMultY=-999;
522  trackMultY=-999;
523  houghMultX=-999;
524  trackMultX=-999;
525  houghAngleY=-999;
526  houghStDevY=-999;
527  houghAngleX=-999;
528  houghStDevX=-999;
529  totalHits=-999;
530  totalADC=-999;
531  }
532  else{
533  // Sort the hits by time to improve the speed.
534  // WARNING: The rest of the code assumes that the hits are time sorted!
535  rb::SortByTime(hitlist);
536 
537  // Slice hits by time (with minimal spatial proximity in each view)
538  std::vector<rb::Cluster> timeSlices = TemporalClusters(hitlist);
539 
540  // Split time-sliced clusters into separate views
541  std::vector<rb::Cluster> timeSlicesY;
542  std::vector<rb::Cluster> timeSlicesX;
543 
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]);
547  }
548  rb::Cluster hitsY(geo::kY);
549  rb::Cluster hitsX(geo::kX);
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] );
553  }
554  hitsY.SetNoise( timeSlices[i].IsNoise() );
555  hitsX.SetNoise( timeSlices[i].IsNoise() );
556  timeSlicesY.push_back(hitsY);
557  timeSlicesX.push_back(hitsX);
558  }
559  noiseCluster.SetNoise(true);
560 
561  if ( fWriteTimeSlices ) {
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]);
564  evt.put(std::move(clustercoltime),fTimeSliceInstance);
565  }
566 
567  // Do track slicing on time slices by view
568  for ( unsigned int i = 0; i < timeSlices.size(); ++i ) {
569  if ( timeSlices[i].IsNoise() ) continue;
570  totalHits = timeSlices[i].NCell();
571  totalADC = timeSlices[i].TotalADC();
572  std::vector<rb::Cluster> trackSlicesY;
573  std::vector<rb::Cluster> trackSlicesX;
574 
575  bool mapItrForward(true); // pick direction based on Y-view Hough angle
576 
577  trackSlicesY = MakeTrackSlices( timeSlicesY[i], noiseCluster, geo::kY, mapItrForward );
578  trackSlicesX = MakeTrackSlices( timeSlicesX[i], noiseCluster, geo::kX, mapItrForward );
579 
580  // Add all of our slices to the event store.
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]);
585 
586  // Fill the ROOT tree and reset the tree variables
587  tree->Fill();
588  houghMultY=-999;
589  trackMultY=-999;
590  houghMultX=-999;
591  trackMultX=-999;
592  houghAngleY=-999;
593  houghStDevY=-999;
594  houghAngleX=-999;
595  houghStDevX=-999;
596  totalHits=-999;
597  totalADC=-999;
598  } // end for
599  }
600 
601  clustercoltracks->push_back(noiseCluster);
602  evt.put(std::move(clustercoltracks),fTrackSliceInstance);
603  }
604 
605  //----------------------------------------------------------------------
606  std::vector<rb::Cluster> AirSlicer::MakeTrackSlices( rb::Cluster timeSlice, rb::Cluster& noiseCluster, geo::View_t view, bool& mapItrForward )
607  {
608  std::vector<rb::Cluster> trackSlices;
609  std::vector<rb::Cluster> finalTrackSlices;
610  if ( timeSlice.IsNoise() ) return trackSlices;
611 
612  art::ServiceHandle<art::TFileService> tfs;
613 
614  art::ServiceHandle<geo::Geometry> geoSvc;
615 
616  double maxHitSepXY(0.);
617  if ( view == geo::kY ) maxHitSepXY = fMaxHitSepY;
618  if ( view == geo::kX ) maxHitSepXY = fMaxHitSepX;
619 
620  rb::HitList thisHitList( timeSlice );
621 
622  // Generate angular Hough map for entire slice
623  std::string viewName("XY");
624  if ( view == geo::kY ) viewName = std::string("Y");
625  if ( view == geo::kX ) viewName = std::string("X");
626 
627  if ( fSaveHists ) {
628  TH2D *mapHist = GetHoughMap( thisHitList );
629  mapHist->SetName(TString::Format("HoughMap_%s_%i_%i_%i",viewName.c_str(),run,subrun,evtNum));
630  tfs->make<TH2D>(*mapHist);
631  delete mapHist;
632  }
633  TH1D *angleHist = GetHoughAngles( thisHitList );
634 
635  angleHist->SetName(TString::Format("HoughAngleHist_%s_%i_%i_%i",viewName.c_str(),run,subrun,evtNum));
636  if ( fSaveHists ) tfs->make<TH1D>(*angleHist);
637  int maxBin = angleHist->GetMaximumBin();
638  if ( angleHist->GetBinContent(maxBin) < fMinHoughPairs ){
639  delete angleHist;
640  return trackSlices; // require minimum number of contributing cell pairs
641  }
642 
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;
650 
651  if ( fitStatus != 0 ){
652  delete agaus;
653  delete angleHist;
654  return trackSlices; // Failed to find maximum Hough angle! Do something?
655  }
656 
657  double maxAngle = fitResult->Parameter(1);
658  double sigAngle = fitResult->Parameter(2);
659  sigAngle=sqrt(sigAngle*sigAngle);
660  delete agaus;
661  delete angleHist;
662 
663  // Generate rho Hough map for entire slice
664  TH1D *rhoHist = GetHoughRhos( thisHitList, maxAngle, sigAngle );
665  rhoHist->SetName(TString::Format("HoughRhoHist_%s_%i_%i_%i",viewName.c_str(),run,subrun,evtNum));
666  if ( fSaveHists ) tfs->make<TH1D>(*rhoHist);
667  int mult=fRhoSpecAna->Search( rhoHist, fRhoPeakSig, "goff nodraw", fRhoPeakThresh );
668 
669  // fill the tree variables for the appropriate view
670  if(view==geo::kY){
671  houghMultY=mult;
672  houghAngleY=maxAngle;
673  houghStDevY=sigAngle;
674  }
675  if(view==geo::kX){
676  houghMultX=mult;
677  houghAngleX=maxAngle;
678  houghStDevX=sigAngle;
679  }
680  delete rhoHist;
681 
682  // Group hits by plane
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; // just to check
686  int plane = timeSlice.AllCells()[hitItr]->Plane();
687  if ( hitsByPlane.find(plane) == hitsByPlane.end() ) hitsByPlane[plane] = rb::Cluster(view);
688  hitsByPlane[plane].Add( timeSlice.AllCells()[hitItr] );
689  }
690 
691  // Loop through plane groupings
692  std::map< int, std::vector<rb::Cluster> > inPlaneClusters;
693  for ( auto planeItr = hitsByPlane.begin(); planeItr != hitsByPlane.end(); ++planeItr ) {
694  int thePlane = planeItr->first;
695  // Cluster hits in this plane
696  int CurrentClusterID = 1;
697  FindNeighbors( planeItr->second, view );
698  for ( unsigned int j = 0; j < planeItr->second.NCell(); ++j ) {
699  // if this cell isn't classified, attempt to expand the cluster
700  if ( fClusterIDs[j] == -1 ) {
701  if ( GenerateCluster( j, CurrentClusterID ) ) CurrentClusterID++;
702  }
703  }
704  // Group hits by cluster in this plane
705  inPlaneClusters[thePlane].resize(CurrentClusterID,rb::Cluster(view));
706  for ( unsigned int j = 0; j < planeItr->second.NCell(); ++j ) {
707  if ( fClusterIDs[j] == 0 ) {
708  noiseCluster.Add( planeItr->second.AllCells()[j] );
709  continue;
710  }
711  inPlaneClusters[thePlane][fClusterIDs[j]].Add( planeItr->second.AllCells()[j] );
712  }
713  } // Now inPlaneClusters contains clusters within each plane.
714 
715  // Loop through planes in this temporal slice
716  if ( view == geo::kY ) mapItrForward = (maxAngle > 0.);
717  // pick the direction for the X- and Y-views from this slice's Y-view, then pass it back for the X-view call
718  auto mapItrStart = inPlaneClusters.begin();
719  auto mapItrEnd = inPlaneClusters.end();
720  if ( !mapItrForward ) {
721  mapItrStart = --(inPlaneClusters.rbegin().base());
722  mapItrEnd = --(inPlaneClusters.rend().base());
723  }
724 
725  // Create first clusters from first plane
726  for ( unsigned int clusItr = 0; clusItr < mapItrStart->second.size(); ++clusItr ) {
727  trackSlices.push_back( mapItrStart->second[clusItr] );
728  }
729 
730  // Loop through remaining planes
731  auto mapItr = mapItrStart;
732  if ( mapItrForward ) ++mapItr;
733  else --mapItr;
734  while ( mapItr != mapItrEnd ) {
735  if ( mapItr->second.size() == 0 ) { // continue
736  if ( mapItrForward ) ++mapItr;
737  else --mapItr;
738  continue;
739  }
740  int thisPlane = mapItr->first;
741 
742  // Match each cluster in this plane with the best seed cluster
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 ) { // cluster 0 is noise
746  if ( mapItr->second[j].NCell() == 0 ) continue;
747  // Find mean position of this cluster in this plane
748  double meanXY(0.);
749  double meanZ(0.);
750  for ( unsigned int c = 0; c < mapItr->second[j].NCell(); ++c ) {
751  art::Ptr<rb::CellHit> hit = mapItr->second[j].Cell(c);
752  std::pair<double,double> hitPos = GetHitPos( hit );
753  meanXY += hitPos.first;
754  meanZ += hitPos.second;
755  }
756  meanXY = meanXY / double(mapItr->second[j].NCell());
757  meanZ = meanZ / double(mapItr->second[j].NCell());
758 
759  // Loop through track clusters
760  for ( unsigned int i = 0; i < trackSlices.size(); ++i ) {
761  // Loop through hits in this track slice that are within the last few planes
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;
768  if ( trackSlices[i].Cell(c)->Plane() > maxPlane ) maxPlane = trackSlices[i].Cell(c)->Plane();
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;
776  nCells++;
777  }
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;
783  posItr++;
784  countItr++;
785  }
786  if ( nCells == 0 ) {
787  continue;
788  }
789  meanXYC = meanXYC / double( nCells );
790  meanZC = meanZC / double( nCells );
791  int nPlanes = maxPlane-minPlane+1;
792  double closestClusterDistXY(1000000.);
793  double closestClusterDistZ(1000000.);
794  double meanHoughAngle(0.), meanHoughAngleSig(0.);
795  double meanHoughRho(0.), meanHoughRhoSig(0.);
796  int nMeanEntries(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;
803  HoughRhoTheta( c->second.second, c->second.first,
804  meanZ, meanXY,
805  &houghRho, &houghAngle, &houghRhoSig, &houghAngleSig );
806  meanHoughAngle += houghAngle;
807  meanHoughRho += houghRho;
808  meanHoughAngleSig += houghAngleSig*houghAngleSig;
809  meanHoughRhoSig += houghAngleSig*houghRhoSig;
810  nMeanEntries++;
811  }
812  if ( closestClusterDistXY > maxHitSepXY || closestClusterDistZ > fMaxHitSepZ ) {
813  continue;
814  }
815  meanHoughAngle = meanHoughAngle / double(nMeanEntries); // should already have excluded nMeanEntries == 0
816  meanHoughRho = meanHoughRho / double(nMeanEntries);
817  meanHoughAngleSig = sqrt(meanHoughAngleSig) / double(nMeanEntries);
818  meanHoughRhoSig = sqrt(meanHoughRhoSig) / double(nMeanEntries);
819  double matchVal;
820  if ( nPlanes <= 2*fHoughExtrapPlanes ) {
821  // Compare Hough angle between mean track slice cluster position and this cluster's position
822  // to the Hough angle of the time slice
823  double houghRho, houghAngle, houghRhoSig, houghAngleSig;
824  HoughRhoTheta( meanZ, meanXY, meanZC, meanXYC, &houghRho, &houghAngle, &houghRhoSig, &houghAngleSig );
825  matchVal = std::min(fabs(houghAngle-maxAngle),M_PI-fabs(houghAngle-maxAngle))/sigAngle;
826  } else {
827  // Compare mean Hough angle between all clusters in the track slice and this cluster's position
828  // to the Hough angle of the time slice
829  matchVal = std::min(fabs(meanHoughAngle-maxAngle),M_PI-fabs(meanHoughAngle-maxAngle))/sigAngle;
830  }
831  if ( matchVal > fMaxHoughMatchVal ) continue;
832  if ( matchVal < bestMatchVal[j] ) {
833  bestMatchInt[j] = i;
834  bestMatchVal[j] = matchVal;
835  }
836  }
837  }
838  for ( unsigned int j = 1; j < mapItr->second.size(); ++j ) { // cluster 0 is noise
839  // If no match found, add this cluster as a new seed
840  if ( bestMatchInt[j] < 0 ) {
841  trackSlices.push_back( mapItr->second[j] );
842  continue;
843  }
844  // Otherwise add the hits from this cluster to the best-matched track cluster
845  for ( unsigned int k = 0; k < mapItr->second[j].NCell(); ++k ) {
846  trackSlices[bestMatchInt[j]].Add( mapItr->second[j].AllCells()[k] );
847  }
848  }
849  if ( mapItrForward ) ++mapItr;
850  else --mapItr;
851  }
852 
853  for ( unsigned int i = 0; i < trackSlices.size(); ++i ) {
854  if ( trackSlices[i].NCell() < fMinTrackHits ) {
855  for ( unsigned int j = 0; j < trackSlices[i].NCell(); ++j ) {
856  noiseCluster.Add( trackSlices[i].AllCells()[j] );
857  }
858  } else {
859  finalTrackSlices.push_back( trackSlices[i] );
860  }
861  }
862  if ( finalTrackSlices.size() < fMinTracks ) finalTrackSlices.clear();
863  if(view==geo::kY) trackMultY=finalTrackSlices.size();
864  if(view==geo::kX) trackMultX=finalTrackSlices.size();
865  return finalTrackSlices;
866  }
867 
868  //----------------------------------------------------------------------
869  bool AirSlicer::GenerateCluster(unsigned int point, int ClID)
870  {
871  // In the original DBSCAN paper, this is the function they called
872  // "ExpandCluster." How this function works is most easily understood
873  // by reading that paper (don't worry, it is pretty short...)
874 
875  // Get the list of neighbors for 'point.'
876  std::vector<unsigned int> seeds = fNeighborList[point];
877 
878  // If 'point' is not a core point, label it as noise and return.
879  if(seeds.size() < fMinPts) {
880  fClusterIDs[point] = 0;
881  return false;
882  }
883  else {
884  // Assign 'point' and the neighbors of 'point' to the cluster ClID.
885  for(unsigned int seedi = 0; seedi < seeds.size(); ++seedi)
886  fClusterIDs[point] = ClID;
887  seeds.erase(seeds.begin());
888 
889  // Loop over the neighbors of 'point' to see if they are also core
890  // points.
891  while(seeds.size() > 0) {
892  // Get the list of neighbors for the first seed point.
893  std::vector<unsigned int> result = fNeighborList[seeds[0]];
894 
895  // If the first seed point is also a core point, add its neighbors
896  // to the list of seeds and assigne them to the cluster ClID.
897  if(result.size() >= fMinPts) {
898  for(unsigned int resulti = 0; resulti < result.size(); ++resulti) {
899  // If the point is unassigned (-1) or previously labeled as
900  // noise (0), asign it ClID.
901  if(fClusterIDs[result[resulti]] == -1 ||
902  fClusterIDs[result[resulti]] == 0) {
903  // Only add this point to the list of seeds if it was
904  // previously unassigned.
905  if(fClusterIDs[result[resulti]] == -1)
906  seeds.push_back(result[resulti]);
907  fClusterIDs[result[resulti]] = ClID;
908  } // end if(resulti is -1 or 0)
909 
910  } // end for resulti
911 
912  } // end if result.size() >= fMinPts
913  seeds.erase(seeds.begin());
914 
915  } // end while
916  return true;
917  } // end else
918 
919  } // end GenerateCluster function
920 
921 } // end namespace airshower
922 /////////////////////////////////////////////////////////////////////////
923 
924 namespace airshower
925 {
926  DEFINE_ART_MODULE(AirSlicer)
927 }
airshower::AirSlicer::fMaxHoughMatchVal
double fMaxHoughMatchVal
maximum match value for Hough extrapolation
Definition: AirSlicer_module.cc:90
art::Event::subRun
SubRunNumber_t subRun() const
Definition: Event.h:72
airshower::AirSlicer::fMinHoughHitSep
double fMinHoughHitSep
minimum separation of hit pairs for Hough map
Definition: AirSlicer_module.cc:87
rawdata::RawTrigger::fTriggerMask_Prescale
uint32_t fTriggerMask_Prescale
Definition: RawTrigger.h:44
angle
Double_t angle
Definition: plot.C:86
rb::Cluster::NCell
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
airshower::AirSlicer::fTimePeakThresh
double fTimePeakThresh
threshold for temporal slice peak search
Definition: AirSlicer_module.cc:76
geo::CellGeo::GetCenter
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
airshower::AirSlicer::fTimeRangeHigh
int fTimeRangeHigh
high edge of time range for temporal slicing (integer microseconds)
Definition: AirSlicer_module.cc:73
airshower::AirSlicer::fMinPts
unsigned int fMinPts
minimum number of neighbors for a hit to be in a cluster
Definition: AirSlicer_module.cc:66
stan::math::fabs
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
gen_flatrecord.size
size
Definition: gen_flatrecord.py:186
y1
Float_t y1[n_points_granero]
Definition: compare.C:5
rb::Cluster::SetNoise
void SetNoise(bool noise)
Declare the cluster to consist of noise hits or not.
Definition: Cluster.h:161
geo::View_t
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
rb::CellHit::Plane
unsigned short Plane() const
Definition: CellHit.h:39
x1
Float_t x1[n_points_granero]
Definition: compare.C:5
airshower::AirSlicer::fSigCellSepX
double fSigCellSepX
scale for X separation between cells in X-view to include in non-noise slice
Definition: AirSlicer_module.cc:79
rb::CellHit::View
geo::View_t View() const
Definition: CellHit.h:41
geo::PlaneGeo::Cell
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
p
const char * p
Definition: xmltok.h:285
std::sqrt
T sqrt(T number)
Definition: d0nt_math.hpp:156
art::Timestamp::timeHigh
constexpr std::uint32_t timeHigh() const
Definition: Timestamp.h:35
cet::pow
constexpr T pow(T x)
Definition: pow.h:75
geo::kX
Vertical planes which measure X.
Definition: PlaneGeo.h:28
airshower::AirSlicer::GetHitPos
std::pair< double, double > GetHitPos(art::Ptr< rb::CellHit > hit)
Definition: AirSlicer_module.cc:370
submit_syst.x2
x2
Definition: submit_syst.py:250
rb::Cluster
A collection of associated CellHits.
Definition: Cluster.h:47
geo::GeometryBase::DetLength
double DetLength() const
Definition: GeometryBase.cxx:481
airshower::AirSlicer::fTimePeakSig
double fTimePeakSig
significance for temporal slice peak search
Definition: AirSlicer_module.cc:75
site_stats_from_log.get
get
Definition: site_stats_from_log.py:189
abs
void abs(TH1 *hist)
Definition: wrong_sign_uncertainty.C:716
submit_syst.y2
y2
Definition: submit_syst.py:250
make_pair
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
rb::SortByTime
void SortByTime(std::vector< art::Ptr< rb::CellHit > > &c)
Sort c in time order (earliest to latest).
Definition: CellHit.cxx:134
geo::GeometryBase::Plane
const PlaneGeo * Plane(unsigned int i) const
Definition: GeometryBase.cxx:749
airshower::AirSlicer::fMinTrackLen
double fMinTrackLen
minimum track length for good tracks
Definition: AirSlicer_module.cc:94
DEFINE_ART_MODULE
DEFINE_ART_MODULE(TestTMapFile)
airshower::AirSlicer::MakeTrackSlices
std::vector< rb::Cluster > MakeTrackSlices(rb::Cluster timeSlice, rb::Cluster &noiseCluster, geo::View_t view, bool &mapItrForward)
Definition: AirSlicer_module.cc:606
M_PI
#define M_PI
Definition: SbMath.h:34
geo::kY
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
rb::Cluster::AllCells
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
airshower::AirSlicer::fSigCellSepZ
double fSigCellSepZ
scale for Z separation between cells to include in non-noise slice
Definition: AirSlicer_module.cc:81
airshower::AirSlicer::fTimeRangeLow
int fTimeRangeLow
low edge of time range for temporal slicing (integer microseconds)
Definition: AirSlicer_module.cc:72
rb::Cluster::Add
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
art::Event::put
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
rb::CellHit::Cell
unsigned short Cell() const
Definition: CellHit.h:40
airshower::AirSlicer::fRawDataLabel
std::string fRawDataLabel
Data label for trigger prescale.
Definition: AirSlicer_module.cc:97
geo::PlaneGeo::View
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
rb::HitList
Just the essential information required for track fitting.
Definition: HitList.h:40
airshower::AirSlicer::fRhoPeakSig
double fRhoPeakSig
significance for rho peak search
Definition: AirSlicer_module.cc:77
airshower::AirSlicer::GetHoughAngles
TH1D * GetHoughAngles(rb::HitList hitList)
Definition: AirSlicer_module.cc:324
airshower::AirSlicer::fTimePeakCut
double fTimePeakCut
cut for temporal slicing
Definition: AirSlicer_module.cc:74
geo::PlaneGeo
Geometry information for a single readout plane.
Definition: PlaneGeo.h:36
evt
int evt
Definition: FillTrainingTrees.h:22
plot_validation_datamc.c
c
Definition: plot_validation_datamc.py:198
d
Float_t d
Definition: plot.C:236
airshower::AirSlicer::HoughRhoTheta
void HoughRhoTheta(double x1, double y1, double x2, double y2, double *rho, double *theta, double *sigmaRho, double *sigmaTheta)
Definition: AirSlicer_module.cc:275
calib::j
const double j
Definition: BetheBloch.cxx:29
airshower::AirSlicer::fMinHoughPairs
double fMinHoughPairs
minimum number of hit pairs in maximum Hough angle or rho histogram to proceed
Definition: AirSlicer_module.cc:86
art::Event::event
EventNumber_t event() const
Definition: Event.h:67
airshower::AirSlicer::fClusterIDs
std::vector< int > fClusterIDs
Definition: AirSlicer_module.cc:101
airshower::AirSlicer::fCellHitInput
std::string fCellHitInput
Read CellHits from this input.
Definition: AirSlicer_module.cc:65
airshower::AirSlicer::fChi2NdofCut
double fChi2NdofCut
cut on Chi2/Ndof for track fit
Definition: AirSlicer_module.cc:93
airshower::AirSlicer::fHoughExtrapPlanes
double fHoughExtrapPlanes
number of planes required for Hough extrapolation
Definition: AirSlicer_module.cc:89
View
Definition: View.py:1
airshower::AirSlicer::fNeighborList
std::vector< std::vector< unsigned int > > fNeighborList
Definition: AirSlicer_module.cc:100
test.z
z
Definition: test.py:28
airshower::AirSlicer::fTrackSliceInstance
std::string fTrackSliceInstance
instance name of track slices
Definition: AirSlicer_module.cc:71
airshower::AirSlicer::CountActiveDCM
int CountActiveDCM(const std::vector< rb::CellHit > &hitlist)
Definition: AirSlicer_module.cc:437
run
Definition: run.py:1
airshower::AirSlicer::fHoughRhoPeakThresh
double fHoughRhoPeakThresh
threshold for Hough map rho peak search
Definition: AirSlicer_module.cc:88
airshower::AirSlicer::fSaveHists
bool fSaveHists
save intermediate histograms
Definition: AirSlicer_module.cc:96
airshower::AirSlicer::TemporalClusters
std::vector< rb::Cluster > TemporalClusters(std::vector< art::Ptr< rb::CellHit > > hitList)
Definition: AirSlicer_module.cc:208
Event.h
febshutoff_auto.dcm
list dcm
Definition: febshutoff_auto.py:63
rb::Cluster::Cell
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
rb::CellHit
A rawdata::RawDigit with channel information decoded.
Definition: CellHit.h:27
airshower::AirSlicer::fMaxSlices
int fMaxSlices
maximum number of time slices to create
Definition: AirSlicer_module.cc:67
art::TFileDirectory::make
T * make(ARGS...args) const
Definition: TFileDirectory.h:82
std::sin
T sin(T number)
Definition: d0nt_math.hpp:132
geo::GeometryBase::DetHalfWidth
double DetHalfWidth() const
Definition: GeometryBase.cxx:467
HoughResult.h
Data resulting from a Hough transform on the cell hit positions.
airshower::AirSlicer::produce
void produce(art::Event &evt)
Definition: AirSlicer_module.cc:460
airshower::AirSlicer::fMaxNeighborScore
double fMaxNeighborScore
neighbor score cutoff
Definition: AirSlicer_module.cc:82
art::DataViewImpl::getByLabel
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
airshower::AirSlicer::fTimeSliceInstance
std::string fTimeSliceInstance
instance name of time slices
Definition: AirSlicer_module.cc:70
hit
Definition: structs.h:12
airshower::AirSlicer::fMinTrackHits
unsigned int fMinTrackHits
minimum number of hits for a track slice
Definition: AirSlicer_module.cc:91
airshower::AirSlicer::GenerateCluster
bool GenerateCluster(unsigned int point, int ClID)
Definition: AirSlicer_module.cc:869
airshower::AirSlicer::fMinTracks
unsigned int fMinTracks
minimum number of tracks in each time slice
Definition: AirSlicer_module.cc:92
airshower::AirSlicer::fMaxHitSepX
double fMaxHitSepX
maximum X separation between hits along track slice in X-view
Definition: AirSlicer_module.cc:83
airshower::AirSlicer::FindNeighbors
void FindNeighbors(rb::Cluster clusterList, geo::View_t view)
Definition: AirSlicer_module.cc:390
std::cos
T cos(T number)
Definition: d0nt_math.hpp:78
art::Event::time
Timestamp time() const
Definition: Event.h:61
airshower::AirSlicer::fSigCellSepY
double fSigCellSepY
scale for Y separation between cells in Y-view to include in non-noise slice
Definition: AirSlicer_module.cc:80
airshower::AirSlicer::fRhoPeakThresh
double fRhoPeakThresh
threshold for rho peak search
Definition: AirSlicer_module.cc:78
cmap
cmap::CMap class source code
Definition: CMap.cxx:17
airshower::AirSlicer::fMaxHitSepZ
double fMaxHitSepZ
maximum Z separation between hits along track slice
Definition: AirSlicer_module.cc:85
airshower::AirSlicer::fWriteTimeSlices
bool fWriteTimeSlices
save time slices to output file
Definition: AirSlicer_module.cc:69
daqchannelmap::DAQChannelMap::getDCM
dcm_id_t getDCM(dchan daqchan) const
Decode the dcm ID from a dchan.
Definition: DAQChannelMap.cpp:115
genie::utils::style::Format
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
geo::GeometryBase::NPlanes
unsigned int NPlanes() const
Definition: GeometryBase.cxx:460
std::min
T min(const caf::Proxy< T > &a, T b)
Definition: BasicTypesProxy.h:404
rb::Cluster::IsNoise
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
e
Float_t e
Definition: plot.C:35
art::Event::run
RunNumber_t run() const
Definition: Event.h:77
airshower::AirSlicer::fMaxHitSepY
double fMaxHitSepY
maximum Y separation between hits along track slice in Y-view
Definition: AirSlicer_module.cc:84
geo::CellGeo
Encapsulate the cell geometry.
Definition: CellGeo.h:25
airshower::AirSlicer::AirSlicer
AirSlicer(fhicl::ParameterSet const &pset)
Definition: AirSlicer_module.cc:131
airshower::AirSlicer::fMaxPeaks
int fMaxPeaks
maximum number of rho peaks to count
Definition: AirSlicer_module.cc:68
airshower::AirSlicer::GetHoughRhos
TH1D * GetHoughRhos(rb::HitList hitList, double angle, double angleSig)
Definition: AirSlicer_module.cc:347
daqchannelmap::DAQChannelMap::getDiBlock
diblock_t getDiBlock(dchan daqchan) const
Decode the diblock ID from a dchan.
Definition: DAQChannelMap.cpp:105
airshower::AirSlicer::GetHoughMap
TH2D * GetHoughMap(rb::HitList hitList)
Definition: AirSlicer_module.cc:302
Geometry.h
Encapsulate the geometry of one entire detector (near, far, ndos)
airshower::AirSlicer::fMaxSliceHits
int fMaxSliceHits
maximum number of hits to perform slicing
Definition: AirSlicer_module.cc:95
std::atan2
T atan2(T number)
Definition: d0nt_math.hpp:72
string
enum BeamMode string
Generated on Thu Dec 3 2020 05:11:29 for NOvA by   doxygen 1.8.11