HorizontalMuonAna_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: HorizontalMuonAna
3 // Plugin Type: filter (art v1_17_06)
4 // File: HorizontalMuonAna_module.cc
5 //
6 // Generated at Fri Dec 8 17:11:54 2017 by Connor Johnson using cetskelgen
7 // from cetlib version v1_15_03.
8 ////////////////////////////////////////////////////////////////////////
9 
19 #include "fhiclcpp/ParameterSet.h"
21 
22 #include <memory>
23 
24 // // ROOT Includes
25 #include <TTree.h>
26 #include <TVector3.h>
27 
35 
36 const uint16_t TOTAL_FD_PLANES = 896;
37 
38 namespace novaddt {
39 
41  public:
42  explicit HorizontalMuonAna(fhicl::ParameterSet const & p);
43  // The compiler-generated destructor is fine for non-base
44  // classes without bare pointers or other resource use.
45 
46  // Plugins should not be copied or assigned.
47  HorizontalMuonAna(HorizontalMuonAna const &) = delete;
51 
52  // Required functions.
53  void analyze(art::Event const & e) override;
54 
55  // Optional functions.
56  void beginJob() override;
57  virtual void endJob() override;
58 
59  private:
60  // FHICL Parameters
63  int16_t fMaxMissing;
64  double fMaxTime;
66  unsigned _prescale;
67 
68  unsigned _trigger_counts = 0;
69  unsigned _after_prescale = 0;
70 
71  // Module variables for endJob
72  int nEvents = 0;
73  int nTracks = 0;
74  int n3DTracks = 0;
75  int nHorzMuons = 0;
76 
77  // Helper functions. Perhaps make public/protected so logic can be reused.
78  bool FrontBackPlanesHits (const std::vector<uint16_t> & missedFrontPlanes, const std::vector<uint16_t> & missedBackPlaness);
79  bool MissingPlanesInBounds(const std::vector<uint16_t> & allMissingPlanes);
80  bool TimingInBounds (const unsigned long long & minTime, const unsigned long long & maxTime);
81  bool IsHorzMuonTrack (const novaddt::HitList & track, unsigned long long & minTime, unsigned long long & maxTime);
82 
83  // Analyzer variables
84  unsigned int fRun; ///< Run number
85  unsigned int fSubRun; ///< Subrun number
86  unsigned int fEvent; ///< Event number
87  TVector3* fStart; ///< Start position of track
88  TVector3* fStop; ///< End position of track
89  double fMissingPlanes; ///< Number of missing planes for passing track.
90  unsigned long long fTimeTDC; ///< Time (in TNS) for track.
91  std::vector<uint16_t> fFrontPlanes;///< Front CellPlanes hit.
92  std::vector<uint16_t> fBackPlanes;
93  std::vector<uint16_t> fHitPlanes;
94  std::vector<unsigned long long> fHitTimes;
95  std::vector<uint8_t> fHitViews;
96 
97  // Art service logic
99  TTree * fTrackTree;
100  };
101 
102  // Constructor from FHICL file
104  :
105  EDAnalyzer(p),
106  fNumFrontBackPlanes( p.get<uint16_t> ("numFrontBackPlanes")),
107  fMinFrontBackPlanes( p.get<int16_t> ("minFrontBackPlanes")),
108  fMaxMissing( p.get<int16_t> ("maxMissingPlanes")),
109  fMaxTime( p.get<double> ("maxTime")),
110  fTracksTag( p.get<std::string> ("tracksTag")),
111  _prescale( p.get<unsigned> ("prescale"))
112 
113  {
114  // Nothing produced into the art event.
115  }
116 
117  // Filter rules, using IsHorzMuonTrack() decision function.
119  {
120  nEvents++;
121 
122  // Grab all Track objects from event
124  // art::Handle<novaddt::HitList> eventHits;
125  e.getByLabel(fTracksTag, eventTracks);
126  if (eventTracks.failedToGet()){
127  mf::LogError("HorizontalMuonAna") << "Error: " << eventTracks.whyFailed()->what() << std::endl;
128  return;
129  }
130 
131  art::FindOneP<novaddt::Track3D> houghTrack(eventTracks, e, fTracksTag);
132 
133  // Loop through tracks. Only one needs to pass to accept the event.
134  for(unsigned long trackIndex = 0; trackIndex < eventTracks->size(); trackIndex++){
135  nTracks++;
136  // Look through each track
137  novaddt::HitList trackHits = eventTracks->at(trackIndex);
138  art::Ptr<novaddt::Track3D> thisTrack3d = houghTrack.at(trackIndex);
139  if (thisTrack3d->Is3D()){
140  n3DTracks++;
141 
142  // Holders for track time (passed by reference)
143  unsigned long long minTime = ULLONG_MAX -1;
144  unsigned long long maxTime = 0;
145  if (IsHorzMuonTrack(trackHits, minTime, maxTime)){
146  nHorzMuons++;
147 
148  fRun = e.run();
149  fSubRun = e.subRun();
150  fEvent = e.event();
151 
152  fStart = new TVector3(thisTrack3d->Start());
153  fStop = new TVector3(thisTrack3d->End());
154  fTimeTDC = maxTime - minTime;
155 
156  fTrackTree->Fill();
157 
158  _trigger_counts++;
159  // if(_trigger_counts%_prescale == _prescale-1){
160  _after_prescale++;
161  // }
162  }
163  }
164  }
165 
166  return;
167  }
168 
169  bool HorizontalMuonAna::FrontBackPlanesHits(const std::vector<uint16_t> & missedFrontPlanes, const std::vector<uint16_t> & missedBackPlanes){
170  // If there is a minimum number of acceptable planes (numFound = Total - numRemaining : numFound >= minFound)
171  if (fMinFrontBackPlanes > 0)
172  return ((int16_t) (fNumFrontBackPlanes - missedFrontPlanes.size()) >= fMinFrontBackPlanes && (int16_t) (fNumFrontBackPlanes - missedBackPlanes.size()) >= fMinFrontBackPlanes);
173 
174  // If no minimum provided, they all must be found.
175  return (missedFrontPlanes.empty() && missedBackPlanes.empty());
176  }
177 
178  bool HorizontalMuonAna::MissingPlanesInBounds(const std::vector<uint16_t> & allMissingPlanes){
179  if (fMaxMissing < 0) return true;
180  return ((int16_t) allMissingPlanes.size() <= fMaxMissing);
181  }
182 
183  bool HorizontalMuonAna::TimingInBounds(const unsigned long long & minTime, const unsigned long long & maxTime){
184  if (fMaxTime < 0) return true;
185  return (maxTime - minTime <= fMaxTime);
186  }
187 
188  bool HorizontalMuonAna::IsHorzMuonTrack(const novaddt::HitList & track, unsigned long long & minTime, unsigned long long & maxTime){
189  // Make list of planes which must have a hit in them.
190  // Different sets for front and back planes, in case requirements differ between them
191  std::vector<uint16_t> RequiredFrontPlanes, RequiredBackPlanes, AllPlanes;
192  for (uint16_t i = 0; i < fNumFrontBackPlanes; i++)
193  RequiredFrontPlanes.push_back(i);
194  for (uint16_t i = 0; i < fNumFrontBackPlanes; i++)
195  RequiredBackPlanes.push_back(TOTAL_FD_PLANES - i - 1);
196  for (uint16_t i = 0; i < TOTAL_FD_PLANES; i++)
197  AllPlanes.push_back(i);
198 
199  fHitPlanes.clear();
200  fHitTimes.clear();
201  fHitViews.clear();
202 
203  // Loop through hits, keeping track of Planes and timing.
204  for (unsigned int i = 0; i < track.size(); i++){
205  const novaddt::DAQHit & hit = track[i];
206  if (hit.TDC().val > maxTime)
207  maxTime = hit.TDC().val;
208  if (hit.TDC().val < minTime)
209  minTime = hit.TDC().val;
210 
211  // If the plane matches any required plane, remove it from the list to search through.
212  RequiredFrontPlanes.erase(
213  std::remove(RequiredFrontPlanes.begin(), RequiredFrontPlanes.end(), hit.Plane().val),
214  RequiredFrontPlanes.end());
215 
216  RequiredBackPlanes.erase(
217  std::remove(RequiredBackPlanes.begin(), RequiredBackPlanes.end(), hit.Plane().val),
218  RequiredBackPlanes.end());
219 
220  AllPlanes.erase(
221  std::remove(AllPlanes.begin(), AllPlanes.end(), hit.Plane().val),
222  AllPlanes.end());
223 
224  fHitPlanes.push_back(hit.Plane().val);
225  fHitTimes.push_back(hit.TDC().val);
226  fHitViews.push_back(hit.View().val);
227  }
228 
229  bool result = TimingInBounds(minTime, maxTime) && MissingPlanesInBounds(AllPlanes) && FrontBackPlanesHits(RequiredFrontPlanes, RequiredBackPlanes);
230 
231  // Fill in analyzer variables if needed.
232  if (result){
233  fMissingPlanes = AllPlanes.size(); // The remianing planes are ones without hits.
234 
235  // Store which planes have hits near the front and back by inverting the "missing" planes that are stored.
236  fFrontPlanes.clear();
237  fBackPlanes.clear();
238  for (uint16_t i = 0; i < fNumFrontBackPlanes; i++){
239  fFrontPlanes.push_back(i);
240  fBackPlanes.push_back(TOTAL_FD_PLANES - i - 1);
241  }
242 
243  // Invert the "missing planes" to only store existing planes in fFrontPlanes and fBackPlanes variables.
244  for (const uint16_t missingPlane : RequiredFrontPlanes)
245  fFrontPlanes.erase(std::remove(fFrontPlanes.begin(), fFrontPlanes.end(), missingPlane), fFrontPlanes.end());
246  for (const uint16_t missingPlane : RequiredBackPlanes)
247  fBackPlanes.erase(std::remove(fBackPlanes.begin(), fBackPlanes.end(), missingPlane), fBackPlanes.end());
248  }
249 
250  return result;
251  }
252 
253  // float HorizontalMuonAna::GetCalibratedTime(double TDC, double dist_to_readout) {
254  // // Taken from UpMuAna code. (Eval/UpMuAna_module.cc)
255  // // This needs to be calibrated and put in a DB table
256  // double speedOfFiberTransport = 15.3; // cm/ns, "first principles" calib.
257  // double TDC_to_ns = 15.625; // 1/64E6=15.625 ns
258 
259  // // Differs from c/n due to zigzag
260  // // paths down fiber. But, this
261  // // value may be way off (10%? 30%?).
262  // return (TDC*TDC_to_ns - dist_to_readout/speedOfFiberTransport);
263  // }
264 
266  fTrackTree = tfs->make<TTree>("HorzMuonTracks", "Horizontal Muon Tracks");
267 
268  // Must instantiate objects before adding to TTree.
269  fStart = new TVector3;
270  fStop = new TVector3;
271 
272  // Define the branches and set the addresses of the TTree
273  fTrackTree->Branch("Run", &fRun);
274  fTrackTree->Branch("SubRun", &fSubRun);
275  fTrackTree->Branch("Event", &fEvent);
276  fTrackTree->Branch("Start", &fStart);
277  fTrackTree->Branch("Stop", &fStop);
278  fTrackTree->Branch("MissingPlanes", &fMissingPlanes);
279  fTrackTree->Branch("Time", &fTimeTDC);
280 
281  fTrackTree->Branch("FrontPlanesHit",&fFrontPlanes);
282  fTrackTree->Branch("BackPlanesHit", &fBackPlanes);
283 
284  fTrackTree->Branch("HitPlanes", &fHitPlanes);
285  fTrackTree->Branch("HitTimes", &fHitTimes);
286  fTrackTree->Branch("HitViews", &fHitViews);
287  }
288 
290  std::cout << "=== novaddt::HorizontalMuonAna endJob" << std::endl;
291  std::cout << "\tNumber of events: " << nEvents << std::endl;
292  std::cout << "\tNumber of tracks: " << nTracks << std::endl;
293  std::cout << "\tNumber of 3D tracks: " << n3DTracks << std::endl;
294  std::cout << "\tNumber of Horizontal Muons: " << nHorzMuons << std::endl;
295  std::cout << "\tNumber of Times Triggered: " << _trigger_counts << std::endl;
296  std::cout << "\tNumber after prescale: " << _after_prescale << std::endl;
297  }
298 
300 
301 } // novaddt namespace
unsigned long long fTimeTDC
Time (in TNS) for track.
TVector3 * fStop
End position of track.
SubRunNumber_t subRun() const
Definition: Event.h:72
art::ServiceHandle< art::TFileService > tfs
std::vector< uint16_t > fBackPlanes
value_type val
Definition: BaseProducts.h:34
double fMissingPlanes
Number of missing planes for passing track.
novaddt::Plane const & Plane() const
Definition: DAQHit.h:70
std::vector< uint16_t > fHitPlanes
std::vector< unsigned long long > fHitTimes
novaddt::TDC const & TDC() const
Definition: DAQHit.h:74
std::vector< DAQHit > HitList
Definition: HitList.h:15
const char * p
Definition: xmltok.h:285
bool TimingInBounds(const unsigned long long &minTime, const unsigned long long &maxTime)
Definition: event.h:19
const uint16_t TOTAL_FD_PLANES
unsigned int fSubRun
Subrun number.
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
DEFINE_ART_MODULE(TestTMapFile)
TVector3 const & End() const
Definition: Track3D.h:46
std::vector< uint16_t > fFrontPlanes
Front CellPlanes hit.
novaddt::View const & View() const
Definition: DAQHit.h:72
bool IsHorzMuonTrack(const novaddt::HitList &track, unsigned long long &minTime, unsigned long long &maxTime)
TVector3 const & Start() const
Definition: Track3D.h:45
value_type val
Definition: BaseProducts.h:84
EventNumber_t event() const
Definition: Event.h:67
bool const & Is3D() const
Definition: Track3D.h:43
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
HorizontalMuonAna(fhicl::ParameterSet const &p)
OStream cout
Definition: OStream.cxx:6
TVector3 * fStart
Start position of track.
bool MissingPlanesInBounds(const std::vector< uint16_t > &allMissingPlanes)
T * make(ARGS...args) const
unsigned int fEvent
Event number.
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
Definition: structs.h:12
value_type val
Definition: BaseProducts.h:137
std::shared_ptr< art::Exception const > whyFailed() const
Definition: Handle.h:219
HorizontalMuonAna & operator=(HorizontalMuonAna const &)=delete
Float_t e
Definition: plot.C:35
RunNumber_t run() const
Definition: Event.h:77
bool FrontBackPlanesHits(const std::vector< uint16_t > &missedFrontPlanes, const std::vector< uint16_t > &missedBackPlaness)
Definition: fwd.h:28
bool failedToGet() const
Definition: Handle.h:196
void analyze(art::Event const &e) override
enum BeamMode string