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 #include <iostream>
24 
25 // // ROOT Includes
26 #include <TTree.h>
27 #include <TVector3.h>
28 #include <TGraph.h>
29 #include <TF1.h>
30 
31 #include "RecoBase/Cluster.h"
32 #include "RecoBase/Track.h"
33 #include "RecoBase/CellHit.h"
34 
35 const uint16_t TOTAL_FD_PLANES = 896;
36 
37 namespace novaddt {
38 
39  class HorizontalMuonAna : public art::EDAnalyzer {
40  public:
41  explicit HorizontalMuonAna(fhicl::ParameterSet const & p);
42  // The compiler-generated destructor is fine for non-base
43  // classes without bare pointers or other resource use.
44 
45  // Plugins should not be copied or assigned.
46  HorizontalMuonAna(HorizontalMuonAna const &) = delete;
50 
51  // Required functions.
52  void analyze(art::Event const & e) override;
53 
54  // Optional functions.
55  void beginJob() override;
56  virtual void endJob() override;
57 
58  private:
59  // FHICL Parameters
60  unsigned short fNumFrontBackPlanes;
62  short fMaxMissing;
63  double fMaxTime;
65  unsigned short fStoreAllHits;
66 
67  // Module variables for endJob
68  int nEvents = 0;
69  int nTracks = 0;
70  int n3DTracks = 0;
71  int nHorzMuons = 0;
72 
73  // Helper functions. Perhaps make public/protected so logic can be reused.
74  bool FrontBackPlanesHits (const std::vector<unsigned short> & missedFrontPlanes, const std::vector<unsigned short> & missedBackPlaness);
75  bool MissingPlanesInBounds(const std::vector<unsigned short> & allMissingPlanes);
76  bool TimingInBounds (const float & minTime, const float & maxTime);
77  bool IsHorzMuonTrack (const rb::Track & track, float & minTime, float & maxTime);
78 
79  // bool FillStartEndInfo (const novaddt::HitList & track, unsigned long long offsetTime);
80 
81  // Analyzer variables
82  unsigned int fRun; ///< Run number
83  unsigned int fSubRun; ///< Subrun number
84  unsigned int fEvent; ///< Event number
85  TVector3* fStart; ///< Start position of track
86  TVector3* fStop; ///< End position of track
87  unsigned int fMissingPlanes; ///< Number of missing planes for passing track.
88  float fTimeTNS; ///< Time (in TNS) for track.
89  float fTotalGeV; ///< TotalGeV deposited in track
90  float fTotalPE; ///< TotalPE deposited in track
91 
92 
93  std::vector<unsigned short> fFrontPlanes;///< Front CellPlanes hit.
94  std::vector<unsigned short> fBackPlanes;
95  std::vector<unsigned short> fHitPlanes;
96  std::vector<float> fHitTimes;
97  std::vector<unsigned short> fHitViews;
98  std::vector<unsigned short> fHitCells;
99 
100  // Art service logic
102  TTree * fTrackTree;
103  };
104 
105  // Constructor from FHICL file
107  :
108  EDAnalyzer(p),
109  fNumFrontBackPlanes( p.get<unsigned short> ("numFrontBackPlanes")),
110  fMinFrontBackPlanes( p.get<short> ("minFrontBackPlanes")),
111  fMaxMissing( p.get<short> ("maxMissingPlanes")),
112  fMaxTime( p.get<double> ("maxTime")),
113  fTracksTag( p.get<std::string> ("tracksTag")),
114  fStoreAllHits( p.get<unsigned short> ("storeAllHits"))
115  {
116  // Nothing produced into the art event.
117  }
118 
119  // Filter rules, using IsHorzMuonTrack() decision function.
121  {
122  nEvents++;
123 
124  // Grab all Track objects from event
126  // art::Handle<novaddt::HitList> eventHits;
127  e.getByLabel(fTracksTag, eventTracks);
128  if (eventTracks.failedToGet()){
129  mf::LogError("HorizontalMuonAna") << "Error: " << eventTracks.whyFailed()->what() << std::endl;
130  return;
131  }
132 
133  // art::FindOneP<novaddt::Track3D> windowtrack(eventTracks, e, fTracksTag);
134 
135  // Loop through tracks. Only one needs to pass to accept the event.
136  for(unsigned long trackIndex = 0; trackIndex < eventTracks->size(); trackIndex++){
137  nTracks++;
138  // Look through each track
139  rb::Track track = eventTracks->at(trackIndex);
140  // art::Ptr<novaddt::Track3D> thisTrack3d = houghTrack.at(trackIndex);
141  if (track.Is3D()){
142  n3DTracks++;
143 
144  // Holders for track time (passed by reference)
145  float minTime = ULLONG_MAX -1;
146  float maxTime = 0;
147  if (IsHorzMuonTrack(track, minTime, maxTime)){
148  nHorzMuons++;
149 
150  fRun = e.run();
151  fSubRun = e.subRun();
152  fEvent = e.event();
153 
154  // this->FillStartEndInfo(trackHits, minTime);
155  delete fStart;
156  delete fStop;
157  fStart = new TVector3(track.Start());
158  fStop = new TVector3(track.Stop());
159  fTimeTNS = maxTime - minTime;
160  fTotalGeV = track.TotalGeV();
161  fTotalPE = track.TotalPE();
162 
163  fTrackTree->Fill();
164  }
165  }
166  }
167 
168  return;
169  }
170 
171  bool HorizontalMuonAna::FrontBackPlanesHits(const std::vector<uint16_t> & missedFrontPlanes, const std::vector<uint16_t> & missedBackPlanes){
172  // If there is a minimum number of acceptable planes (numFound = Total - numRemaining : numFound >= minFound)
173  if (fMinFrontBackPlanes > 0)
174  return ((int16_t) (fNumFrontBackPlanes - missedFrontPlanes.size()) >= fMinFrontBackPlanes && (int16_t) (fNumFrontBackPlanes - missedBackPlanes.size()) >= fMinFrontBackPlanes);
175 
176  // If no minimum provided, they all must be found.
177  return (missedFrontPlanes.empty() && missedBackPlanes.empty());
178  }
179 
180  bool HorizontalMuonAna::MissingPlanesInBounds(const std::vector<uint16_t> & allMissingPlanes){
181  if (fMaxMissing < 0) return true;
182  return ((int16_t) allMissingPlanes.size() <= fMaxMissing);
183  }
184 
185  bool HorizontalMuonAna::TimingInBounds(const float & minTime, const float & maxTime){
186  if (fMaxTime < 0) return true;
187  return (maxTime - minTime <= fMaxTime);
188  }
189 
190  bool HorizontalMuonAna::IsHorzMuonTrack(const rb::Track & track, float & minTime, float & maxTime){
191  // Initial obvious tests, to make sure we're not doing something stupid.
192  if (track.TotalLength() < 4000 || track.IsNoise())
193  return false;
194 
195  // At this point, study the track
196  const art::PtrVector<rb::CellHit> trackCells = track.AllCells();
197  // Make list of planes which must have a hit in them.
198  // Different sets for front and back planes, in case requirements differ between them
199  std::vector<uint16_t> RequiredFrontPlanes, RequiredBackPlanes, AllPlanes;
200  for (uint16_t i = 0; i < fNumFrontBackPlanes; i++)
201  RequiredFrontPlanes.push_back(i);
202  for (uint16_t i = 0; i < fNumFrontBackPlanes; i++)
203  RequiredBackPlanes.push_back(TOTAL_FD_PLANES - i - 1);
204  for (uint16_t i = 0; i < TOTAL_FD_PLANES; i++)
205  AllPlanes.push_back(i);
206 
207  fHitPlanes.clear();
208  fHitTimes.clear();
209  fHitViews.clear();
210  fHitCells.clear();
211 
212  // Loop through hits, keeping track of Planes and timing.
213  for (unsigned int i = 0; i < trackCells.size(); i++){
214  const art::Ptr<rb::CellHit> hit = trackCells[i];
215  if (hit->TNS() > maxTime)
216  maxTime = hit->TNS();
217  if (hit->TNS() < minTime)
218  minTime = hit->TNS();
219 
220  // If the plane matches any required plane, remove it from the list to search through.
221  RequiredFrontPlanes.erase(
222  std::remove(RequiredFrontPlanes.begin(), RequiredFrontPlanes.end(), hit->Plane()),
223  RequiredFrontPlanes.end());
224 
225  RequiredBackPlanes.erase(
226  std::remove(RequiredBackPlanes.begin(), RequiredBackPlanes.end(), hit->Plane()),
227  RequiredBackPlanes.end());
228 
229  AllPlanes.erase(
230  std::remove(AllPlanes.begin(), AllPlanes.end(), hit->Plane()),
231  AllPlanes.end());
232 
233  fHitPlanes.push_back(hit->Plane());
234  fHitTimes.push_back(hit->TNS());
235  fHitViews.push_back(hit->View());
236  fHitCells.push_back(hit->Cell());
237  }
238 
239  bool result = TimingInBounds(minTime, maxTime) && MissingPlanesInBounds(AllPlanes) && FrontBackPlanesHits(RequiredFrontPlanes, RequiredBackPlanes);
240 
241  // Fill in analyzer variables if needed.
242  if (result){
243  fMissingPlanes = AllPlanes.size(); // The remianing planes are ones without hits.
244 
245  // Store which planes have hits near the front and back by inverting the "missing" planes that are stored.
246  fFrontPlanes.clear();
247  fBackPlanes.clear();
248  for (uint16_t i = 0; i < fNumFrontBackPlanes; i++){
249  fFrontPlanes.push_back(i);
250  fBackPlanes.push_back(TOTAL_FD_PLANES - i - 1);
251  }
252 
253  // Invert the "missing planes" to only store existing planes in fFrontPlanes and fBackPlanes variables.
254  for (const uint16_t missingPlane : RequiredFrontPlanes)
255  fFrontPlanes.erase(std::remove(fFrontPlanes.begin(), fFrontPlanes.end(), missingPlane), fFrontPlanes.end());
256  for (const uint16_t missingPlane : RequiredBackPlanes)
257  fBackPlanes.erase(std::remove(fBackPlanes.begin(), fBackPlanes.end(), missingPlane), fBackPlanes.end());
258  }
259 
260  return result;
261  }
262 
264  fTrackTree = tfs->make<TTree>("HorzMuonTracks", "Horizontal Muon Tracks");
265 
266  // Must instantiate objects before adding to TTree.
267  fStart = new TVector3;
268  fStop = new TVector3;
269 
270  // Define the branches and set the addresses of the TTree
271  fTrackTree->Branch("Run", &fRun);
272  fTrackTree->Branch("SubRun", &fSubRun);
273  fTrackTree->Branch("Event", &fEvent);
274  fTrackTree->Branch("Start", &fStart);
275  fTrackTree->Branch("Stop", &fStop);
276  fTrackTree->Branch("MissingPlanes", &fMissingPlanes);
277  // fTrackTree->Branch("StartTime", &fStartTimeTNS);
278  // fTrackTree->Branch("EndTime", &fEndTimeTNS);
279  fTrackTree->Branch("Time", &fTimeTNS);
280 
281  fTrackTree->Branch("FrontPlanesHit",&fFrontPlanes);
282  fTrackTree->Branch("BackPlanesHit", &fBackPlanes);
283 
284  if (fStoreAllHits){
285  fTrackTree->Branch("HitPlanes", &fHitPlanes);
286  fTrackTree->Branch("HitTimes", &fHitTimes);
287  fTrackTree->Branch("HitViews", &fHitViews);
288  fTrackTree->Branch("HitCells", &fHitCells);
289  }
290 
291  fTrackTree->Branch("TotalGeV", &fTotalGeV);
292  fTrackTree->Branch("TotalPE", &fTotalPE);
293  }
294 
296  std::cout << "=== novaddt::HorizontalMuonAna endJob" << std::endl;
297  std::cout << "\tNumber of events: " << nEvents << std::endl;
298  std::cout << "\tNumber of tracks: " << nTracks << std::endl;
299  std::cout << "\tNumber of 3D tracks: " << n3DTracks << std::endl;
300  std::cout << "\tNumber of Horizontal Muons: " << nHorzMuons << std::endl;
301  }
302 
304 
305 } // novaddt namespace
TVector3 * fStop
End position of track.
float TNS() const
Definition: CellHit.h:46
double TotalPE() const
Sum of the PE value of all the contained hits.
Definition: Cluster.cxx:369
SubRunNumber_t subRun() const
Definition: Event.h:72
art::ServiceHandle< art::TFileService > tfs
std::vector< uint16_t > fBackPlanes
double fMissingPlanes
Number of missing planes for passing track.
std::vector< uint16_t > fHitPlanes
unsigned short Plane() const
Definition: CellHit.h:39
std::vector< unsigned long long > fHitTimes
geo::View_t View() const
Definition: CellHit.h:41
const char * p
Definition: xmltok.h:285
std::vector< unsigned short > fFrontPlanes
Front CellPlanes hit.
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.
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
std::vector< unsigned short > fHitViews
virtual double TotalLength() const
Length (cm) of all the track segments.
Definition: Track.cxx:213
std::vector< uint16_t > fFrontPlanes
Front CellPlanes hit.
unsigned short Cell() const
Definition: CellHit.h:40
std::vector< unsigned short > fHitPlanes
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
bool IsHorzMuonTrack(const novaddt::HitList &track, unsigned long long &minTime, unsigned long long &maxTime)
EventNumber_t event() const
Definition: Event.h:67
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
HorizontalMuonAna(fhicl::ParameterSet const &p)
size_type size() const
Definition: PtrVector.h:308
OStream cout
Definition: OStream.cxx:6
TVector3 * fStart
Start position of track.
bool MissingPlanesInBounds(const std::vector< uint16_t > &allMissingPlanes)
float fTotalPE
TotalPE deposited in track.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
float fTimeTNS
Time (in TNS) for track.
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
double TotalGeV(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple sum of the estimated GeV of all the hits.
Definition: Cluster.cxx:378
std::vector< unsigned short > fBackPlanes
std::shared_ptr< art::Exception const > whyFailed() const
Definition: Handle.h:219
HorizontalMuonAna & operator=(HorizontalMuonAna const &)=delete
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
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)
float fTotalGeV
TotalGeV deposited in track.
unsigned int fMissingPlanes
Number of missing planes for passing track.
std::vector< unsigned short > fHitCells
virtual bool Is3D() const
Definition: Prong.h:71
bool failedToGet() const
Definition: Handle.h:196
void analyze(art::Event const &e) override