SNMichelAnalyzer_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \brief Cluster Michels and output TTrees
3 // \author Justin Vasel <justin.vasel@gmail.com>
4 // \date September 2019
5 ////////////////////////////////////////////////////////////////////////
6 
7 // STL includes
8 #include <deque>
9 #include <iostream>
10 
11 // Framework includes
19 #include "fhiclcpp/ParameterSet.h"
20 
21 // NOvASoft includes
22 #include "Geometry/Geometry.h"
25 #include "NovaTimingUtilities/TimingUtilities.h"
26 #include "RawData/DAQHeader.h"
27 #include "RawData/RawTrigger.h"
28 #include "RecoBase/CellHit.h"
29 #include "RecoBase/Cluster.h"
30 #include "RecoBase/Track.h"
31 
32 // ROOT includes
33 #include "TTree.h"
34 #include "TVector3.h"
35 
36 
37 namespace sn {
38  class SNMichelAnalyzer;
39  class MichelBackgroundManager;
40 
42  bool SortTrackEndsByTime(const std::pair<int, std::pair<double, TVector3>>& t1, const std::pair<int, std::pair<double, TVector3>>& t2);
43  bool SortClustersByTime(const rb::Cluster& c1, const rb::Cluster& c2);
44 
45  enum TrackType {
46  kUnknown, // Track type is none of the below
47  kStopping, // Track enters the detector and stops inside
48  kThroughGoing, // Track enters and exits the detector
49  kExiting, // Track starts inside the detector and exits
50  kContained // Track starts and ends entirely inside the detector
51  };
52 }
53 
54 
55 ///////////////////////////////////////////////////////////////////////////////
57  public:
58  MichelBackgroundManager(unsigned int size) {this->fSize = size;};
60 
61  void AddEntry(std::vector<std::pair<int, std::pair<double, TVector3>>> entry);
62  void SetSize(unsigned int size) {this->fSize = size;}
63 
64  std::vector<std::pair<int, std::pair<double, TVector3>>> GetLastEntry();
65  unsigned int Size() {return fSize;};
66  size_t NEntries() {return fEntries.size();};
67  bool IsReady() {return fEntries.size() == fSize;};
68  std::deque<std::vector<std::pair<int, std::pair<double, TVector3>>>> GetEntries() const {return fEntries;};
69 
70  friend std::ostream& operator<<(std::ostream& os, const sn::MichelBackgroundManager& m);
71 
72  private:
73  unsigned int fSize;
74  std::deque<std::vector<std::pair<int, std::pair<double, TVector3>>>> fEntries;
75 };
76 
77 namespace sn {
78  std::ostream& operator<<(std::ostream& os, const sn::MichelBackgroundManager& m)
79  {
80  os << "[ ";
81  for (std::vector<std::pair<int, std::pair<double, TVector3>>> entry : m.GetEntries()) {
82  os << entry.size() << " ";
83  }
84  os << "]";
85  return os;
86  }
87 }
88 
89 // ............................................................................
90 void sn::MichelBackgroundManager::AddEntry(std::vector<std::pair<int, std::pair<double, TVector3>>> entry)
91 {
92  if (this->fEntries.size() == fSize) this->fEntries.pop_back();
93  fEntries.push_front(entry);
94 
95  assert(fEntries.size() <= fSize);
96 
97  return;
98 }
99 
100 
101 // ............................................................................
102 std::vector<std::pair<int, std::pair<double, TVector3>>> sn::MichelBackgroundManager::GetLastEntry()
103 {
104  return this->fEntries.back();
105 }
106 
107 
108 
109 ///////////////////////////////////////////////////////////////////////////////
111 public:
112  explicit SNMichelAnalyzer(fhicl::ParameterSet const & p);
113 
114  SNMichelAnalyzer(SNMichelAnalyzer const &) = delete;
115  SNMichelAnalyzer(SNMichelAnalyzer &&) = delete;
116  SNMichelAnalyzer & operator = (SNMichelAnalyzer const &) = delete;
117  SNMichelAnalyzer & operator = (SNMichelAnalyzer &&) = delete;
118 
119  // Required functions
120  void beginJob();
121  void analyze(art::Event const & e) override;
122 
123  void ClusterAroundTrackEnds(std::vector<std::pair<int, std::pair<double, TVector3>>> trackEnds, art::Event const & e, TTree*& tree, bool kIsBackground=false);
124  double DistanceHitToPoint(const rb::CellHit hit, const TVector3 point);
125  bool PointIsContained(TVector3 point);
126  sn::TrackType GetTrackType(rb::Track track);
127 
128 private:
129  std::string fRawDataLabel; ///< ART module label for raw data
130  std::string fCellHitLabel; ///< ART module label for cell hits
131  std::string fClusterLabel; ///< ART module label for clusters
132  std::string fTrackLabel; ///< ART module label for tracks
133 
134  bool fBackgroundSampleEnable; ///< Switch to enable collecting a background sample
135  int fBackgroundSampleDepth; ///< Number of events to look back for background sample
136 
137  int fMuonEndSphereRadius; ///< Distance from end of muon track to veto
138  int fMuonEndTimeCut; ///< Time from end of muon track to veto
139 
140  int fEvt; ///< Event number
141  int fRun; ///< Run number
142  int fSubrun; ///< Subrun number
143  int fMichelSumADC; ///< Total Michel cluster ADC
144  float fMichelSumGeV; ///< Total Michel cluster GeV
145  float fMichelMeanTNS; ///< Mean Michel cluster time (ns)
146  int fMichelNHits; ///< Number of hits in Michel cluster
147  int fMichelTrackId; ///< Identifier for stopping tracks
148  float fMichelDistance; ///< Distance from Michel cluster to track
149  float fMichelTimeDiff; ///< Time diff from Michel cluster to track
150  int fTrigNMicroslices; ///< Number of microslices in trigger
151  unsigned long long fTrigUnixTimeStart; ///< Unix time of trigger (in ns)
152  uint32_t fTrigLength; ///< Duration of trigger (in ns)
153 
154  TTree* fTreeMichels;
156 
159 };
160 
161 
162 // ............................................................................
164 EDAnalyzer(p),
165 fRawDataLabel(p.get<std::string>("RawDataLabel")),
166 fCellHitLabel(p.get<std::string>("CellHitLabel")),
167 fClusterLabel(p.get<std::string>("ClusterLabel")),
168 fTrackLabel(p.get<std::string>("TrackLabel")),
169 fBackgroundSampleEnable(p.get<bool>("BackgroundSampleEnable")),
170 fBackgroundSampleDepth(p.get<int>("BackgroundSampleDepth")),
171 fMuonEndSphereRadius(p.get<int>("MuonEndSphereRadius")),
172 fMuonEndTimeCut(p.get<int>("MuonEndTimeCut")),
173 fEvt(0),
174 fRun(0),
175 fSubrun(0),
176 fTrigNMicroslices(0),
177 fTrigUnixTimeStart(0),
178 fTrigLength(0)
179 {
181 }
182 
183 
184 // ............................................................................
186 {
188 
189  fTreeMichels = tfs->make<TTree>("Michels", "Michels");
190  fTreeMichels->Branch("Run", &fRun);
191  fTreeMichels->Branch("Subrun", &fSubrun);
192  fTreeMichels->Branch("Event", &fEvt);
193  fTreeMichels->Branch("TrigNMicroslices", &fTrigNMicroslices);
194  fTreeMichels->Branch("TrigUnixTimeStart", &fTrigUnixTimeStart);
195  fTreeMichels->Branch("TrigLength", &fTrigLength);
196  fTreeMichels->Branch("MichelSumADC", &fMichelSumADC);
197  fTreeMichels->Branch("MichelSumGeV", &fMichelSumGeV);
198  fTreeMichels->Branch("MichelMeanTNS", &fMichelMeanTNS);
199  fTreeMichels->Branch("MichelNHits", &fMichelNHits);
200  fTreeMichels->Branch("MichelTrackId", &fMichelTrackId);
201  fTreeMichels->Branch("MichelDistance", &fMichelDistance);
202  fTreeMichels->Branch("MichelTimeDiff", &fMichelTimeDiff);
203 
205  fTreeBackground = tfs->make<TTree>("Backgrounds", "Backgrounds");
206  fTreeBackground->Branch("MichelSumGeV", &fMichelSumGeV);
207  fTreeBackground->Branch("MichelMeanTNS", &fMichelMeanTNS);
208  fTreeBackground->Branch("MichelNHits", &fMichelNHits);
209  fTreeBackground->Branch("MichelTrackId", &fMichelTrackId);
210  fTreeBackground->Branch("MichelDistance", &fMichelDistance);
211  fTreeBackground->Branch("MichelTimeDiff", &fMichelTimeDiff);
212  }
213 
214  return;
215 }
216 
217 
218 // ............................................................................
219 void sn::SNMichelAnalyzer::ClusterAroundTrackEnds(std::vector<std::pair<int, std::pair<double, TVector3>>> trackEnds, art::Event const & e, TTree*& tree, bool kIsBackground)
220 {
221  /* Get CellHits from event, put them into a vector, and sort them by time */
223  e.getByLabel(fCellHitLabel, hdlHits);
224  std::vector<art::Ptr<rb::CellHit>> hits;
225  hits.reserve(hdlHits->size());
226  for (unsigned int idx=0; idx<hdlHits->size(); ++idx){
227  art::Ptr<rb::CellHit> hit(hdlHits, idx);
228  hits.push_back(hit);
229  }
230  rb::SortByTime(hits);
231 
232  std::sort(trackEnds.begin(), trackEnds.end(), sn::SortTrackEndsByTime);
233 
234  // This will allow us to extract the hits from a track as a cluster
236  e.getByLabel(fTrackLabel, hdlTracks);
237 
238  art::FindOneP<rb::Cluster> clusterFromTrack(hdlTracks, e, fTrackLabel);
239 
240  /* Loop over track ends and find nearby hits */
241  unsigned int idxStartCellHit = 0;
242  for (std::pair<int, std::pair<double, TVector3>> trackEnd : trackEnds) {
243  rb::Cluster* michelCluster = new rb::Cluster();
244  unsigned int trackId = trackEnd.first;
245  double trackTNS = trackEnd.second.first;
246  TVector3 xyz = trackEnd.second.second;
247 
248  art::PtrVector<rb::CellHit> trackHits;
249  if (!kIsBackground) {
250  trackHits = clusterFromTrack.at(trackId)->AllCells();
251  }
252 
253  for (size_t idx=idxStartCellHit; idx<hits.size(); ++idx) {
254  art::Ptr<rb::CellHit> hit = hits.at(idx);
255  if (hit->TNS() < trackTNS) {
256  idxStartCellHit = idx;
257  continue;
258  }
259 
260  if (!kIsBackground) {
261  if (std::find(trackHits.begin(), trackHits.end(), hit) != trackHits.end()) {
262  idxStartCellHit = idx;
263  continue;
264  }
265  }
266 
267  if (hit->TNS() > (trackTNS + fMuonEndTimeCut)) {
268  break;
269  }
270 
271  TVector3 cellCenter;
272  fGeom->Plane(hit->Plane())->Cell(hit->Cell())->GetCenter(cellCenter);
273 
274  float fDistance = this->DistanceHitToPoint(*hit, xyz);
275  if (fDistance < fMuonEndSphereRadius) {
276  michelCluster->Add(hit);
277  michelCluster->SetID(trackId);
278  }
279  }
280 
281  if (michelCluster->NXCell() > 0 && michelCluster->NYCell() > 0) {
282  fMichelNHits = michelCluster->NCell();
283  fMichelTrackId = michelCluster->ID();
284 
285  size_t cSize = michelCluster->NCell();
286  fMichelSumADC = (cSize==0) ? 0 : michelCluster->TotalADC();
287  fMichelSumGeV = (cSize==0) ? 0 : michelCluster->TotalGeV();
288  fMichelMeanTNS = (cSize==0) ? 0 : michelCluster->MeanTNS();
289  fMichelDistance = (cSize==0) ? 0 : (michelCluster->MeanXYZ() - xyz).Mag();
290  fMichelTimeDiff = (cSize==0) ? 0 : michelCluster->MeanTNS() - trackTNS;
291  tree->Fill();
292  }
293 
294  delete michelCluster;
295  }
296 
297  return;
298 }
299 
300 
301 // ............................................................................
303 {
304  fRun = e.run();
305  fSubrun = e.subRun();
306  ++fEvt;
307 
308  /* Get DAQHeader from event */
310  e.getByLabel(fRawDataLabel, daqheader);
311  fTrigNMicroslices = daqheader->TotalMicroSlices();
312 
313  /* Get RawTrigger info from event */
315  e.getByLabel(fRawDataLabel, trigs);
316  const rawdata::RawTrigger trig = trigs->at(0);
317 
318  /* Get trigger length and start time in Unix epoch (ms) */
319  struct timespec unixtime;
321 
322  fTrigUnixTimeStart = unixtime.tv_sec*1e3 + unixtime.tv_nsec*1e-6;
323  fTrigLength = trig.fTriggerRange_TriggerLength * 500; // units: ns
324 
325  /* Get Tracks from event, put them in a ptrVector, and sort them by time */
327  e.getByLabel(fTrackLabel, hdlTracks);
328  std::vector<art::Ptr<rb::Track>> tracks;
329  art::fill_ptr_vector(tracks, hdlTracks);
330  std::sort(tracks.begin(), tracks.end(), sn::SortTracksByTime);
331 
332  /* Get all the stoppoing track end points */
333  std::vector<std::pair<int, std::pair<double, TVector3>>> trackEnds;
334  for (size_t trackId=0; trackId<tracks.size(); ++trackId) {
335  rb::Track track = hdlTracks->at(trackId);
336  sn::TrackType trackType = this->GetTrackType(track);
337  double trackTNS = track.MeanTNS();
338 
339  // Only add the track end point if it's a stopping track
340  if (trackType == sn::TrackType::kStopping){
341  trackEnds.push_back(
342  std::pair<int, std::pair<double, TVector3>>(
343  trackId,
344  std::pair<double, TVector3>(trackTNS, track.Stop())
345  )
346  );
347  }
348  }
349 
351  if (fBkgManager->IsReady()) {
352  this->ClusterAroundTrackEnds(trackEnds, e, fTreeMichels, false);
354  }
355  fBkgManager->AddEntry(trackEnds);
356  } else {
357  this->ClusterAroundTrackEnds(trackEnds, e, fTreeMichels, false);
358  }
359 
360  return;
361 }
362 
363 
364 // ............................................................................
366 {
367  if (!PointIsContained(track.Start()) && PointIsContained(track.Stop())) return sn::TrackType::kStopping;
368  else if (!PointIsContained(track.Start()) && !PointIsContained(track.Stop())) return sn::TrackType::kThroughGoing;
369  else if ( PointIsContained(track.Start()) && PointIsContained(track.Stop())) return sn::TrackType::kContained;
370  else if ( PointIsContained(track.Start()) && !PointIsContained(track.Stop())) return sn::TrackType::kExiting;
371  else return sn::TrackType::kUnknown;
372 }
373 
374 
375 // ............................................................................
376 /**
377  * Return the distance between a cell hit and a point in space.
378  *
379  * @param[in] hit <code>rb::CellHit</code>
380  * @param[in] point TVector3
381  * @return distance between the hit and the point
382  */
383 double sn::SNMichelAnalyzer::DistanceHitToPoint(const rb::CellHit hit, const TVector3 point)
384 {
385  // Get hit and track vectors in XYZ-space
386  TVector3 x0;
387  TVector3 x1 = point;
388 
389  // Note: x0 is passed by reference here, and thus assigned at this line
390  fGeom->Plane(hit.Plane())->Cell(hit.Cell())->GetCenter(x0);
391 
392  // Define ZW coordinate system based on the view of the hit
393  TVector3 w(0,0,0);
394  if (hit.View() == geo::kX) w.SetX(1);
395  if (hit.View() == geo::kY) w.SetY(1);
396 
397  // Transform from XYZ-space into ZW-space
398  TVector3 w0(x0.Z(), x0.Dot(w), 0);
399  TVector3 w1(x1.Z(), x1.Dot(w), 0);
400 
401  return (w1-w0).Mag();
402 }
403 
404 
405 // ............................................................................
407 {
408  return (point.X() > -720) && (point.X() < 720) &&
409  (point.Y() > -720) && (point.Y() < 700) &&
410  (point.Z() > 80) && (point.Z() < 5920);
411 }
412 
413 // ............................................................................
415 {
416  return t1->MinTNS() < t2->MinTNS();
417 }
418 
419 // ............................................................................
420 bool sn::SortTrackEndsByTime(const std::pair<int, std::pair<double, TVector3>>& t1, const std::pair<int, std::pair<double, TVector3>>& t2)
421 {
422  return t1.second.first < t2.second.first;
423 }
424 
425 // ............................................................................
427 {
428  return c1.MinTNS() < c2.MinTNS();
429 }
430 
int fBackgroundSampleDepth
Number of events to look back for background sample.
float TNS() const
Definition: CellHit.h:46
sn::MichelBackgroundManager * fBkgManager
SubRunNumber_t subRun() const
Definition: Event.h:72
int fTrigNMicroslices
Number of microslices in trigger.
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
float fMichelSumGeV
Total Michel cluster GeV.
void AddEntry(std::vector< std::pair< int, std::pair< double, TVector3 >>> entry)
unsigned short Plane() const
Definition: CellHit.h:39
Float_t x1[n_points_granero]
Definition: compare.C:5
iterator begin()
Definition: PtrVector.h:223
geo::View_t View() const
Definition: CellHit.h:41
uint32_t fTrigLength
Duration of trigger (in ns)
int fMichelSumADC
Total Michel cluster ADC.
const char * p
Definition: xmltok.h:285
Vertical planes which measure X.
Definition: PlaneGeo.h:28
int fMuonEndSphereRadius
Distance from end of muon track to veto.
Definition: event.h:19
A collection of associated CellHits.
Definition: Cluster.h:47
Track enters the detector and stops inside.
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
void SortByTime(std::vector< art::Ptr< rb::CellHit > > &c)
Sort c in time order (earliest to latest).
Definition: CellHit.cxx:134
const PlaneGeo * Plane(unsigned int i) const
TVector3 MeanXYZ(rb::AveragingScheme=kDefaultScheme) const
Definition: Cluster.cxx:538
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
int TotalMicroSlices() const
Definition: DAQHeader.h:27
bool convertNovaTimeToUnixTime(uint64_t const &inputNovaTime, struct timespec &outputUnixTime)
Track type is none of the below.
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
float fMichelMeanTNS
Mean Michel cluster time (ns)
bool PointIsContained(TVector3 point)
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
int fMichelTrackId
Identifier for stopping tracks.
std::deque< std::vector< std::pair< int, std::pair< double, TVector3 > > > > fEntries
c2
Definition: demo5.py:33
unsigned short Cell() const
Definition: CellHit.h:40
double MinTNS() const
Definition: Cluster.cxx:482
unsigned long long fTrigUnixTimeStart
Unix time of trigger (in ns)
double TotalADC() const
Sum of the ADC of all the contained hits.
Definition: Cluster.cxx:360
Track starts inside the detector and exits.
void hits()
Definition: readHits.C:15
Track enters and exits the detector.
SNMichelAnalyzer(fhicl::ParameterSet const &p)
void beginJob()
std::string fClusterLabel
ART module label for clusters.
iterator end()
Definition: PtrVector.h:237
double DistanceHitToPoint(const rb::CellHit hit, const TVector3 point)
void analyze(art::Event const &e) override
bool SortTracksByTime(const art::Ptr< rb::Track > &t1, const art::Ptr< rb::Track > &t2)
Remove hits from hot and cold channels.
float fMichelDistance
Distance from Michel cluster to track.
void SetID(int id)
Definition: Cluster.h:74
std::vector< std::pair< int, std::pair< double, TVector3 > > > GetLastEntry()
double t2
std::string fCellHitLabel
ART module label for cell hits.
void ClusterAroundTrackEnds(std::vector< std::pair< int, std::pair< double, TVector3 >>> trackEnds, art::Event const &e, TTree *&tree, bool kIsBackground=false)
std::string fRawDataLabel
ART module label for raw data.
std::deque< std::vector< std::pair< int, std::pair< double, TVector3 > > > > GetEntries() const
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
bool SortClustersByTime(const rb::Cluster &c1, const rb::Cluster &c2)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
A rawdata::RawDigit with channel information decoded.
Definition: CellHit.h:27
T * make(ARGS...args) const
Track starts and ends entirely inside the detector.
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
unsigned long long fTriggerTimingMarker_TimeStart
Definition: RawTrigger.h:38
bool fBackgroundSampleEnable
Switch to enable collecting a background sample.
Definition: event.h:1
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
double MeanTNS(rb::AveragingScheme scheme=kDefaultScheme) const
Definition: Cluster.cxx:554
art::ServiceHandle< geo::Geometry > fGeom
bool SortTrackEndsByTime(const std::pair< int, std::pair< double, TVector3 >> &t1, const std::pair< int, std::pair< double, TVector3 >> &t2)
float Mag() const
int fMichelNHits
Number of hits in Michel cluster.
assert(nhit_max >=nhit_nbins)
c1
Definition: demo5.py:24
double TotalGeV(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple sum of the estimated GeV of all the hits.
Definition: Cluster.cxx:378
uint32_t fTriggerRange_TriggerLength
Definition: RawTrigger.h:40
const int ID() const
Definition: Cluster.h:75
sn::TrackType GetTrackType(rb::Track track)
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:464
std::string fTrackLabel
ART module label for tracks.
Float_t e
Definition: plot.C:35
RunNumber_t run() const
Definition: Event.h:77
Float_t w
Definition: plot.C:20
def entry(str)
Definition: HTMLTools.py:26
int fMuonEndTimeCut
Time from end of muon track to veto.
Encapsulate the geometry of one entire detector (near, far, ndos)
friend std::ostream & operator<<(std::ostream &os, const sn::MichelBackgroundManager &m)
MichelBackgroundManager(unsigned int size)
float fMichelTimeDiff
Time diff from Michel cluster to track.