ThroughgoingSelection_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file ThroughgoingSelection.cxx
3 // \brief Selects throughgoing cosmics
4 // \author Christopher Backhouse - bckhouse@caltech.edu
5 ////////////////////////////////////////////////////////////////////////
6 
13 
14 // NOvASoft includes
15 #include "Geometry/Geometry.h"
16 #include "MCCheater/BackTracker.h"
17 #include "MEFinder/MEClusters.h"
18 #include "NovaDAQConventions/DAQConventions.h"
19 #include "RecoBase/CellHit.h"
20 #include "RecoBase/Track.h"
21 
22 namespace calib
23 {
24  /// Selects throughgoing cosmics
26  {
27  public:
28  explicit ThroughgoingSelection(const fhicl::ParameterSet& pset);
30 
31  virtual bool filter(art::Event& evt);
32 
33  void reconfigure(const fhicl::ParameterSet& pset);
34 
35  void beginJob();
36  void endJob();
37 
38  protected:
39  bool IsUncontained(const rb::Track& trk) const;
40  bool IsTrulyUncontained(const rb::Track& trk) const;
41 
42  int fTable[2][2];
43 
44  };
45 
46  //......................................................................
48  {
49  produces< std::vector<rb::Track> >();
50 
51  reconfigure(pset);
52 
53  for(int i = 0; i < 2; ++i) for(int j = 0; j < 2; ++j) fTable[i][j] = 0;
54  }
55 
56  //......................................................................
58  {
59  }
60 
61  //......................................................................
63  {
64  }
65 
66  //......................................................................
68  {
69  }
70 
71  //......................................................................
73  {
74  std::unique_ptr<std::vector<rb::Track> > outTracks(new std::vector<rb::Track>);
75 
77  evt.getByLabel("cosmictrack", tracks);
78 
79  art::FindManyP<me::TrkME> fmp(tracks, evt, "mefinder");
80 
81  for(size_t trIdx = 0; trIdx < tracks->size(); ++trIdx){
82  // If there's any Michel associated, don't take this track
83  if(!fmp.at(trIdx).empty()) continue;
84 
85  const rb::Track& track = (*tracks)[trIdx];
86 
87  const bool thru = IsUncontained(track);
88 
89  if(thru) outTracks->push_back(track);
90 
91  if(!evt.isRealData()){
92  const bool trueThru = IsTrulyUncontained(track);
93  ++fTable[thru][trueThru];
94  }
95  } // end for trIdx
96 
97  const bool ret = !outTracks->empty();
98  evt.put(std::move(outTracks));
99  return ret;
100  }
101 
102  //......................................................................
104  {
105  if(fTable[0][0]+fTable[1][0] == 0 || fTable[0][1]+fTable[1][1] == 0) return;
106  std::cout << "Throughgoing efficiency table. x-axis selection, y-axis truth:" << std::endl;
107  std::cout << fTable[0][0] << "\t" << fTable[1][0] << "\t" << .1*(1000*fTable[0][0]/(fTable[0][0]+fTable[1][0])) << "%" << std::endl;
108  std::cout << fTable[0][1] << "\t" << fTable[1][1] << "\t" << .1*(1000*fTable[1][1]/(fTable[0][1]+fTable[1][1])) << "%" << std::endl;
109  std::cout << .1*(1000*fTable[0][0]/(fTable[0][0]+fTable[0][1])) << "%\t" << .1*(1000*fTable[1][1]/(fTable[1][0]+fTable[1][1])) << "%" << std::endl;
110  std::cout << "Purity = " << .1*(1000*fTable[1][1]/(fTable[1][0]+fTable[1][1])) << "%" << std::endl;
111  std::cout << "Efficiency = " << .1*(1000*fTable[1][1]/(fTable[0][1]+fTable[1][1])) << "%" << std::endl;
112  }
113 
114  //......................................................................
116  {
117  // Start and end don't mean much, pick whichever is lowest in y.
118  const TVector3 v1 = trk.Start();
119  const TVector3 v2 = trk.Stop();
120  const TVector3 stop = (v1.Y() < v2.Y()) ? v1 : v2;
121 
122  if(v1.X() == v2.X() ||
123  v1.Y() == v2.Y()) return false; // dx or dy = 0 means bad fit
124 
125  const double x = stop.X();
126  const double y = stop.Y();
127  const double z = stop.Z();
128 
130 
131  if(fabs(x) > geom->DetHalfWidth()-20) return true;
132  if(fabs(y) > geom->DetHalfHeight()-20) return true;
133  if(z < 20 || z > geom->DetLength()-20) return true;
134 
135  return false;
136  }
137 
138  //......................................................................
141  {
143 
144  std::vector<const sim::Particle*> parts = bt->HitsToParticle(trk.AllCells());
145  if(parts.empty()) return false;
146 
147  const sim::Particle* part = parts.front();
148  if(abs(part->PdgCode()) != 13) return false;
149 
150  const TVector3 stop = part->EndPosition().Vect();
151 
153 
154  if(fabs(stop.X()) > geom->DetHalfWidth()) return true;
155  if(fabs(stop.Y()) > geom->DetHalfHeight()) return true;
156  if(stop.Z() < 0) return true;
157  if(stop.Z() > geom->DetLength()) return true;
158 
159  return false;
160  }
161 
163 
164 } // end namespace
165 ////////////////////////////////////////////////////////////////////////
back track the reconstruction to the simulation
int PdgCode() const
Definition: MCParticle.h:211
bool IsTrulyUncontained(const rb::Track &trk) const
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
const TLorentzVector & EndPosition() const
Definition: MCParticle.h:224
Definition: event.h:19
double DetLength() const
void abs(TH1 *hist)
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
bool isRealData() const
Definition: Event.h:83
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
bool IsUncontained(const rb::Track &trk) const
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
Track finder for cosmic rays.
TString part[npart]
Definition: Style.C:32
CDPStorage service.
int evt
const double j
Definition: BetheBloch.cxx:29
Selects throughgoing cosmics.
double DetHalfHeight() const
z
Definition: test.py:28
ThroughgoingSelection(const fhicl::ParameterSet &pset)
OStream cout
Definition: OStream.cxx:6
double DetHalfWidth() const
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
void reconfigure(const fhicl::ParameterSet &pset)
virtual bool filter(art::Event &evt)
void geom(int which=0)
Definition: geom.C:163
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
Encapsulate the geometry of one entire detector (near, far, ndos)
std::vector< const sim::Particle * > HitsToParticle(const std::vector< const rb::CellHit * > &hits) const
Returns vector of sim::Particle objects contributing to the given collection of hits.