BreakPointAna_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: BreakPointAna
3 // Module Type: analyzer
4 // File: BreakPointAna_module.cc
5 //
6 // Author: Michael Baird (mbaird42@fnal.gov)
7 ////////////////////////////////////////////////////////////////////////
8 
9 #include <cmath>
10 #include <iostream>
11 #include <algorithm>
12 
13 // Histogram includes
14 #include "TH1.h"
15 #include "TH2.h"
16 #include "TMath.h"
17 
18 // ART includes
32 #include "fhiclcpp/ParameterSet.h"
34 
36 
37 // Required for CAFAna Cuts
39 
40 #include "RecoBase/Cluster.h"
41 #include "RecoBase/Track.h"
42 #include "RecoBase/Vertex.h"
43 #include "RecoBase/FitSum.h"
44 #include "RecoBase/RecoHit.h"
45 
47 
48 #include "Utilities/AssociationUtil.h"
49 #include "Utilities/func/MathUtil.h"
50 #include "RecoBase/FilterList.h"
51 
52 namespace bpfit {
53  class BreakPointAna;
54 }
55 
57 public:
58  explicit BreakPointAna(fhicl::ParameterSet const & p);
60 
61  // Required functions.
62  void analyze(art::Event const & e) override;
63  void beginJob();
64  void endJob();
65 
66 private:
67 
68  // fcl parameters
69  std::string fSlicerLabel; ///< slicer label
70  std::string fVertexLabel; ///< vertex label
71  std::string fProngLabel; ///< prong label
72  std::string fProngInstance; ///< prong instance label
73  std::string fTrackLabel; ///< tracker label
74  std::string fFitSumLabel; ///< fitsum maker label
75  std::string fBPFPIdLabel; ///< BPF PID object label
76  int fTrackPDG; ///< which BPF tracking assumption to make plots for
77 
78  std::string fCAFLabel; ///< standard record maer label
79  bool fApplyCAFCuts; ///< apply the cuts?
80  int fCutType; ///< which cuts to apply
81 
82  recovalid::CAFCutter fCAFCutter; ///< the cafcutter class...
83 
84 
85 
86  //
87  // histos:
88  //
89 
90  // slice/prong-level histos
91  TH1F *fNTracks; ///< number of tracks per slice
92  TH1F *fNTracksPerProng; ///< number of tracks per prong
93  TH1F *fNPDGs; ///< number of times a track was made for each tracking assumption
94 
95  // track-level histos
96  TH1F *fLength; ///< track length
97  TH1F *fNHits; ///< number of hits per track
98  TH1F *fTheta; ///< track theta angle
99  TH1F *fPhi; ///< track phi angle
100  TH1F *fStartX; ///< track start X
101  TH1F *fStartY; ///< track start Y
102  TH1F *fStartZ; ///< track start Z
103  TH1F *fStopX; ///< track stop X
104  TH1F *fStopY; ///< track stop Y
105  TH1F *fStopZ; ///< track stop Z
106 
107  TH1F *fPID; ///< track PID value
108  TH1F *fChi2T; ///< track total chi^2 (input to the PID)
109  TH1F *fdEdxLL; ///< track dE/dx log liklihood (input to the PID)
110  TH1F *fHitRatio; ///< track/prong hit ratio (input to the PID)
111 
112  // best-muon histos
113  TH1F *fBestMu_PID; ///< PID value for highest scoring muon track
114  TH1F *fBestMu_Length; ///< length for highest scoring muon track
115  TH1F *fBestMu_Chi2T; ///< total chi^2 for highest scoring muon track
116  TH1F *fBestMu_dEdxLL; ///< dE/dx log likelihood for highest scoring muon track
117  TH1F *fBestMu_HitRatio; ///< track/prong hit ratio for highest scoring muon track
118 
119 };
120 
121 
123  :
124  EDAnalyzer(p),
125  fSlicerLabel (p.get< std::string >("SlicerLabel" )),
126  fVertexLabel (p.get< std::string >("VertexLabel" )),
127  fProngLabel (p.get< std::string >("ProngLabel" )),
128  fProngInstance (p.get< std::string >("ProngInstance" )),
129  fTrackLabel (p.get< std::string >("TrackLabel" )),
130  fFitSumLabel (p.get< std::string >("FitSumLabel" )),
131  fBPFPIdLabel (p.get< std::string >("BPFPIdLabel" )),
132  fTrackPDG (p.get< int >("TrackPDG" )),
133  fCAFLabel (p.get< std::string >("CAFLabel" )),
134  fApplyCAFCuts (p.get< bool >("ApplyCAFCuts" )),
135  fCutType (p.get< int >("CutType" ))
136 {}
137 
138 //......................................................................
140 {
141  //======================================================================
142  // Strike me down and I shall become more powerful than you could
143  // possibly imagine...
144  //======================================================================
145 }
146 
147 
149 {
151 
152  // book histos:
153  fNTracks = tfs->make<TH1F>("NTracks",
154  "Number of tracks per slice;N Tracks;count",
155  31, -0.5, 30.5);
156 
157  fNTracksPerProng = tfs->make<TH1F>("NTracksPerProng",
158  "Number of tracks per prong;N Tracks;count",
159  6, -0.5, 5.5);
160 
161  fNPDGs = tfs->make<TH1F>("NPDGs",
162  "Number of times a track was made for each particle assumption;0 = mu / 1 = pi / 2 = p;count",
163  3, -0.5, 2.5);
164 
165  fLength = tfs->make<TH1F>("Length",
166  "Track length;length [cm];",
167  1200, 0.0, 6000.0);
168 
169  fNHits = tfs->make<TH1F>("NHits",
170  "Number of hits per track;NHits;",
171  1001, -0.5, 1000.5);
172 
173  fTheta = tfs->make<TH1F>("Theta",
174  "Track #theta;#theta [deg.];",
175  180, 0.0, 180.0);
176 
177  fPhi = tfs->make<TH1F>("Phi",
178  "Track #phi;#phi [deg.];",
179  360, 0.0, 360.0);
180 
181  fStartX = tfs->make<TH1F>("StartX",
182  "Track start X;Start X [cm];",
183  800, -800.0, 800.0);
184 
185  fStartY = tfs->make<TH1F>("StartY",
186  "Track start Y;Start Y [cm];",
187  800, -800.0, 800.0);
188 
189  fStartZ = tfs->make<TH1F>("StartZ",
190  "Track start Z;Start Z [cm];",
191  3200, -200.0, 6200.0);
192 
193  fStopX = tfs->make<TH1F>("StopX",
194  "Track stop X;Stop X [cm];",
195  800, -800.0, 800.0);
196 
197  fStopY = tfs->make<TH1F>("StopY",
198  "Track stop Y;Stop Y [cm];",
199  800, -800.0, 800.0);
200 
201  fStopZ = tfs->make<TH1F>("StopZ",
202  "Track stop Z;Stop Z [cm];",
203  3200, -200.0, 6200.0);
204 
205  fPID = tfs->make<TH1F>("PID",
206  "Muon PID value;PID;",
207  220, -1.1, 1.1);
208 
209  fChi2T = tfs->make<TH1F>("Chi2T",
210  "Track Fit total Chi-squared per D.O.F.;#chi^{2}/Ndof;",
211  150, 0.0, 15.0);
212 
213  fdEdxLL = tfs->make<TH1F>("dEdxLL",
214  "dE/dx Log-Likelihood;LL;",
215  700, -9.0, -2.0);
216 
217  fHitRatio = tfs->make<TH1F>("HitRatio",
218  "Track/Prong hit ratio;hit ratio;",
219  100, 0.0, 1.0);
220 
221 
222 
223  fBestMu_PID = tfs->make<TH1F>("BestMu_PID",
224  "Muon PID value - Best Muon;PID;",
225  220, -1.1, 1.1);
226 
227  fBestMu_Length = tfs->make<TH1F>("BestMu_Length",
228  "Track length - Best Muon;length [cm];",
229  1200, 0.0, 6000.0);
230 
231  fBestMu_Chi2T = tfs->make<TH1F>("BestMu_Chi2T",
232  "Track Fit total Chi-squared per D.O.F. - Best Muon;#chi^{2}/Ndof;",
233  150, 0.0, 15.0);
234 
235  fBestMu_dEdxLL = tfs->make<TH1F>("BestMu_dEdxLL",
236  "dE/dx Log-Likelihood - Best Muon;LL;",
237  700, -9.0, -2.0);
238 
239  fBestMu_HitRatio = tfs->make<TH1F>("BestMu_HitRatio",
240  "Track/Prong hit ratio - Best Muon;hit ratio;",
241  100, 0.0, 1.0);
242 
243 }
244 
246 {
247 }
248 
250 {
251 
252  // Get the slices and put them into a more convenient container.
254  evt.getByLabel(fSlicerLabel, slhandle);
256  for(unsigned int i=0; i < slhandle->size(); ++i) {
257  slices.push_back(art::Ptr<rb::Cluster>(slhandle,i));
258  }
259 
260  // Get the Standard Records
261  art::FindManyP<caf::StandardRecord> recordFMP(slices, evt, fCAFLabel);
262 
263  // Create the FindMany object for getting prongs.
265 
266 
267 
268  //
269  // Loop over all slices and get the fuzzyk 3D prongs:
270  //
271  for(unsigned int islice = 0; islice < slhandle->size(); ++islice) {
272 
273  unsigned int NTracks = 0;
274  float BestMu_PID = -10.0;
275  float BestMu_Length = -10.0;
276  float BestMu_Chi2T = -10.0;
277  float BestMu_dEdxLL = -10.0;
278  float BestMu_HitRatio = -10.0;
279 
280  // As usual, skip the noise slice...
281  if((*slhandle)[islice].IsNoise()) continue;
282 
283  // Apply the CAF-level cuts
284  bool pass = true;
285  if( fApplyCAFCuts ) {
286  // get record associated with the slice
287  std::vector< art::Ptr<caf::StandardRecord> > records = recordFMP.at(islice);
288  if( records.size() > 0 && recordFMP.isValid()) {
289  pass = fCAFCutter.passCuts(fCutType, records[0].get());
290  }
291  else { pass = false; }
292  }
293 
294  if(!pass) continue;
295 
296 
297 
298  // NOTE: Currently, the fundamental object of analysis is the slice.
299  // If in the future we change things so that the fundamental object
300  // of analysis becomes the vertex, or if we start putting in more
301  // than one vertex per slice, then we will have to get prongs
302  // associated with the vertex instead.
303 
304  // get the 3D prongs associated with this slice
305  std::vector< art::Ptr< rb::Prong > > prongs;
306  if(prongFMP.isValid()) prongs = prongFMP.at(islice);
307 
308  // create the FindMany object for getting tracks
309  art::FindManyP<rb::Track> trackFMP(prongs, evt, fTrackLabel);
310 
311 
312 
313  //
314  // Loop over all prongs and get the tracks:
315  //
316  for(unsigned int iprong = 0; iprong < prongs.size(); ++iprong) {
317 
318  // get the tracks associated with this prong
319  std::vector< art::Ptr< rb::Track > > tracks;
320  if(trackFMP.isValid()) tracks = trackFMP.at(iprong);
321 
322  // create the FindMany object for getting FitSum objects
323  art::FindManyP<rb::FitSum> fitsumFMP(tracks, evt, fTrackLabel);
324 
325  // create the FindMany object for getting the bpf PID objects
326  art::FindManyP<bpfit::BPFPId> bpfpidFMP(tracks, evt, fBPFPIdLabel);
327 
328 
329 
330  //
331  // Loop over all tracks...
332  //
333  for(unsigned int itrack = 0; itrack < tracks.size(); ++itrack) {
334 
335  NTracks++;
336 
337  // get the FitSum object associated with this track
338  std::vector< art::Ptr< rb::FitSum > > fitsums;
339  if(fitsumFMP.isValid()) fitsums = fitsumFMP.at(itrack);
340 
341  // get the BPFPId object associated with this track
342  std::vector< art::Ptr< bpfit::BPFPId > > bpfpids;
343  if(bpfpidFMP.isValid()) bpfpids = bpfpidFMP.at(itrack);
344 
345  if(fitsums.size() == 0) continue; // this should never happen, but just in case...
346 
347 
348 
349  //
350  // Fill histos...
351  //
352 
353  if(fitsums[0]->PDG() == 13) fNPDGs->Fill(0);
354  if(fitsums[0]->PDG() == 211) fNPDGs->Fill(1);
355  if(fitsums[0]->PDG() == 2212) fNPDGs->Fill(2);
356 
357  // only fill track plots for the requested tracking particle assumption
358  if(fitsums[0]->PDG() == fTrackPDG) {
359 
360  fLength->Fill(tracks[itrack]->TotalLength());
361  fNHits ->Fill(tracks[itrack]->NCell());
362 
363  float dirX = tracks[itrack]->Dir().X();
364  float dirY = tracks[itrack]->Dir().Y();
365  float dirZ = tracks[itrack]->Dir().Z();
366 
367  float theta = acos(dirZ)*180.0/3.14159;
368  fTheta->Fill(theta);
369 
370  float phi = (atan2(dirY,dirX))*180.0/3.14159;
371  if(phi < 0.0) phi += 360.0;
372  fPhi->Fill(phi);
373 
374  fStartX->Fill(tracks[itrack]->Start().X());
375  fStartY->Fill(tracks[itrack]->Start().Y());
376  fStartZ->Fill(tracks[itrack]->Start().Z());
377 
378  fStopX->Fill(tracks[itrack]->Stop().X());
379  fStopY->Fill(tracks[itrack]->Stop().Y());
380  fStopZ->Fill(tracks[itrack]->Stop().Z());
381 
382  // NOTE: Currently, we only make the muon PID for the muon assumption track.
383  if(bpfpids.size() > 0) {
384  float value = (float)bpfpids[0]->Value();
385  float Chi2T = bpfpids[0]->GetChi2T();
386  float dEdxLL = bpfpids[0]->GetdEdXLL();
387  float HitRatio = bpfpids[0]->GetHitRatio();
388 
389  fPID->Fill(value);
390  fChi2T->Fill(Chi2T);
391  fdEdxLL->Fill(dEdxLL);
392  fHitRatio->Fill(HitRatio);
393 
394  if(value > BestMu_PID && fitsums[0]->PDG() == 13) {
395  BestMu_PID = value;
396  BestMu_Length = tracks[itrack]->TotalLength();
397  BestMu_Chi2T = Chi2T;
398  BestMu_dEdxLL = dEdxLL;
399  BestMu_HitRatio = HitRatio;
400  }
401  }
402 
403  }
404 
405  } // end loop over tracks
406 
407 
408 
409  //
410  // Fill prong-level plots
411  //
412  fNTracksPerProng->Fill(tracks.size());
413 
414  } // end loop over prongs
415 
416 
417 
418  //
419  // Fill slice-level plots
420  //
421  fNTracks->Fill(NTracks);
422 
423  if(BestMu_PID > -10.0) {
424  fBestMu_PID ->Fill(BestMu_PID);
425  fBestMu_Length ->Fill(BestMu_Length);
426  fBestMu_Chi2T ->Fill(BestMu_Chi2T);
427  fBestMu_dEdxLL ->Fill(BestMu_dEdxLL);
428  fBestMu_HitRatio->Fill(BestMu_HitRatio);
429  }
430 
431  } // end loop over slices
432 
433  return;
434 
435 } // end analyze module
436 
438 
439 
TH1F * fStopZ
track stop Z
TH1F * fBestMu_Length
length for highest scoring muon track
std::string fProngLabel
prong label
TH1F * fPID
track PID value
BreakPointAna(fhicl::ParameterSet const &p)
const char * p
Definition: xmltok.h:285
TH1F * fStopX
track stop X
TH1F * fBestMu_Chi2T
total chi^2 for highest scoring muon track
TH1F * fChi2T
track total chi^2 (input to the PID)
recovalid::CAFCutter fCAFCutter
the cafcutter class...
TH1F * fStartZ
track start Z
TH1F * fNPDGs
number of times a track was made for each tracking assumption
TH1F * fTheta
track theta angle
DEFINE_ART_MODULE(TestTMapFile)
T acos(T number)
Definition: d0nt_math.hpp:54
TH1F * fNHits
number of hits per track
TH1F * fBestMu_HitRatio
track/prong hit ratio for highest scoring muon track
TH1F * fPhi
track phi angle
Float_t Y
Definition: plot.C:38
TH1F * fStartX
track start X
Float_t Z
Definition: plot.C:38
bool fApplyCAFCuts
apply the cuts?
bool passCuts(int cut, const caf::StandardRecord *sr)
Definition: CAFCutter.cxx:37
TH1F * fBestMu_PID
PID value for highest scoring muon track.
void analyze(art::Event const &e) override
const XML_Char int const XML_Char * value
Definition: expat.h:331
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
int evt
TH1F * fNTracks
number of tracks per slice
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
std::string fSlicerLabel
slicer label
Vertex location in position and time.
std::string fCAFLabel
standard record maer label
int fTrackPDG
which BPF tracking assumption to make plots for
TH1F * fBestMu_dEdxLL
dE/dx log likelihood for highest scoring muon track
TH1F * fNTracksPerProng
number of tracks per prong
TH1F * fLength
track length
TH1F * fdEdxLL
track dE/dx log liklihood (input to the PID)
std::string fFitSumLabel
fitsum maker label
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
std::string fProngInstance
prong instance label
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
int fCutType
which cuts to apply
std::string fBPFPIdLabel
BPF PID object label.
std::string fTrackLabel
tracker label
TH1F * fStartY
track start Y
Helper class for Reco Validation modules.
Definition: CAFCutter.h:58
std::string fVertexLabel
vertex label
TH1F * fStopY
track stop Y
Float_t e
Definition: plot.C:35
Float_t X
Definition: plot.C:38
T atan2(T number)
Definition: d0nt_math.hpp:72
TH1F * fHitRatio
track/prong hit ratio (input to the PID)