StopperSelection_module.cc
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////
2 // \file StopperSelection.cxx
3 // \brief Selects stopping cosmics
4 // \author bckhouse@caltech.edu, bays@caltech.edu
5 ////////////////////////////////////////////////////////////////////////
6 
13 
14 // NOvASoft includes
15 #include "MCCheater/BackTracker.h"
16 #include "MEFinder/MEClusters.h"
17 #include "RecoBase/CellHit.h"
18 #include "RecoBase/Track.h"
19 #include "Geometry/LiveGeometry.h"
20 #include "Geometry/Geometry.h"
21 #include "NovaDAQConventions/DAQConventions.h"
22 #include "Utilities/AssociationUtil.h"
23 
24 namespace calib
25 {
26  /// Selects stopping cosmics
28  {
29  public:
30  explicit StopperSelection(const fhicl::ParameterSet& pset);
32 
33  virtual bool filter(art::Event& evt);
34 
35  void reconfigure(const fhicl::ParameterSet& pset);
36 
37  void beginJob();
38  void endJob();
39 
40  protected:
41  bool IsContainedStopper(const rb::Track* trk, const bool isData) const;
42  bool IsTrulyContainedStopper(const rb::Track* trk) const;
43 
47  double fMinDistToEdge; ///< minimum distance in cm to edge of detector
48  double fMaxElectronEnergy; ///< maximum energy for the electron
49  int fTable[2][2];
50 
53 
54  };
55 
56  //......................................................................
58  {
59  produces< std::vector<rb::Track> >();
60  produces< art::Assns<rb::Track, rb::Cluster> >();
61 
62  reconfigure(pset);
63 
64  for(int i = 0; i < 2; ++i) for(int j = 0; j < 2; ++j) fTable[i][j] = 0;
65  }
66 
67  //......................................................................
69  {
70  }
71 
72  //......................................................................
74  {
75  fMichelLabel = pset.get< std::string >("MichelLabel" );
76  fTrackLabel = pset.get< std::string >("TrackLabel" );
77  fApplyFilter = pset.get< bool >("ApplyFilter" );
78  fMinDistToEdge = pset.get< double >("MinDistToEdge", 50. );
79  fMaxElectronEnergy = pset.get< double >("MaxElectronEnergy", 0.035);
80 
81  return;
82  }
83 
84  //......................................................................
86  {
87  }
88 
89  //......................................................................
91  {
92  std::unique_ptr<std::vector<rb::Track> > outTracks(new std::vector<rb::Track>);
93  std::unique_ptr<art::Assns<rb::Track, rb::Cluster> > assns(new art::Assns<rb::Track, rb::Cluster>);
94 
96  evt.getByLabel(fTrackLabel, trkInfo);
97 
98  art::FindManyP<rb::Cluster> fmpslc(trkInfo, evt, fTrackLabel);
99  art::FindManyP<me::TrkME> fmpMEs(trkInfo, evt, fMichelLabel);
100 
101  for (unsigned int trkIdx = 0; trkIdx < trkInfo->size(); ++trkIdx){
102  const rb::Track* track = &(trkInfo->at(trkIdx));
103  bool isStopper = IsContainedStopper(track, evt.isRealData());
104  bool hasME = !fmpMEs.at(trkIdx).empty();
105  if (!isStopper) continue;
106  if (hasME)
107  outTracks->push_back(*track);
108  if (!evt.isRealData()){
109  const bool trueStopper = IsTrulyContainedStopper(track);
110  ++fTable[hasME][trueStopper];
111  }
112  // Find the slice the original track was associated to
113  if (hasME){
114  const std::vector<art::Ptr<rb::Cluster> > slices = fmpslc.at(trkIdx);
115  assert(slices.size() == 1);
116  util::CreateAssn(*this, evt, *outTracks, slices[0], *assns);
117  }
118  }
119 
120  const bool ret = !outTracks->empty();
121  evt.put(std::move(outTracks));
122  evt.put(std::move(assns));
123  if(fApplyFilter) return ret; else return true;
124  }
125 
126  //......................................................................
128  {
129  // Ensure all the denominators are non-zero
130  if(fTable[0][0]+fTable[1][0] == 0 ||
131  fTable[0][1]+fTable[1][1] == 0 ||
132  fTable[0][0]+fTable[0][1] == 0 ||
133  fTable[1][0]+fTable[1][1] == 0) return;
134 
135  std::cout << "Stopper efficiency table. x-axis selection, y-axis truth:" << std::endl
136  << fTable[0][0] << "\t" << fTable[1][0] << "\t" << .1*(1000*fTable[0][0]/(fTable[0][0]+fTable[1][0])) << "%" << std::endl
137  << fTable[0][1] << "\t" << fTable[1][1] << "\t" << .1*(1000*fTable[1][1]/(fTable[0][1]+fTable[1][1])) << "%" << std::endl
138  << .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
139  << "Purity = " << .1*(1000*fTable[1][1]/(fTable[1][0]+fTable[1][1])) << "%" << std::endl
140  << "Efficieny = " << .1*(1000*fTable[1][1]/(fTable[0][1]+fTable[1][1])) << "%" << std::endl;
141  }
142 
143  //......................................................................
144  bool StopperSelection::
145  IsContainedStopper(const rb::Track* trk, const bool isData) const
146  {
147  // Start and end don't mean much, pick whichever is lowest in y.
148  const TVector3 v1 = trk->Start();
149  const TVector3 v2 = trk->Stop();
150  const TVector3 stop = (v1.Y() < v2.Y()) ? v1 : v2;
151 
152 
153  //For TB manual stopping selection, half instrumented det
154  if(fGeom->DetId()==novadaq::cnv::kTESTBEAM && isData){
155  if( fLiveGeom->DistToEdgeXY(stop) < fMinDistToEdge ||
156  ((fLiveGeom->DistToFront(stop) < fMinDistToEdge) ||
157  (fLiveGeom->DistToBack(stop) < ((fLiveGeom->InstrumentedDetLength())*(32/63) + fMinDistToEdge)))){
158  return false;
159  }
160  }
161  else if(fLiveGeom->DistToEdgeXY(stop) < fMinDistToEdge ||
162  fLiveGeom->DistToEdgeZ(stop) < fMinDistToEdge ) return false;
163  if(v1.X() == v2.X() || v1.Y() == v2.Y()) return false; // dx or dy = 0 means bad fit
164 
165  return true;
166  }
167 
168  //......................................................................
170  {
172 
173  std::vector<const sim::Particle*> parts = bt->HitsToParticle(trk->AllCells());
174  if(parts.empty()) return false;
175 
176  const sim::Particle* part = parts.front();
177 
178  const TVector3 stop = part->EndPosition().Vect();
179 
180  if(fLiveGeom->DistToEdgeXY(stop) < 0.1 ||
181  fLiveGeom->DistToEdgeZ(stop) < 0.1 ) return false;
182  if(std::abs(part->PdgCode())!=13) return false;
183 
184  return true;
185  }
186 
188 
189 } // end namespace
190 ////////////////////////////////////////////////////////////////////////
void reconfigure(const fhicl::ParameterSet &pset)
back track the reconstruction to the simulation
static bool CreateAssn(art::EDProducer const &prod, art::Event &evt, std::vector< T > &a, art::Ptr< U > b, art::Assns< T, U > &assn, size_t indx=UINT_MAX, std::string const &instance=std::string())
Create a 1 to 1 association between a new product and one already in the event.
int PdgCode() const
Definition: MCParticle.h:211
const TLorentzVector & EndPosition() const
Definition: MCParticle.h:224
double fMaxElectronEnergy
maximum energy for the electron
double DistToEdgeXY(TVector3 vertex)
double fMinDistToEdge
minimum distance in cm to edge of detector
Definition: event.h:19
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
double DistToFront(TVector3 vertex)
double DistToBack(TVector3 vertex)
bool isRealData() const
Definition: Event.h:83
bool IsTrulyContainedStopper(const rb::Track *trk) const
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
virtual bool filter(art::Event &evt)
bool IsContainedStopper(const rb::Track *trk, const bool isData) const
float abs(float number)
Definition: d0nt_math.hpp:39
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.
double InstrumentedDetLength()
get instrumented detector length of downstream part
TString part[npart]
Definition: Style.C:32
CDPStorage service.
double DistToEdgeZ(TVector3 vertex)
Selects stopping cosmics.
T get(std::string const &key) const
Definition: ParameterSet.h:231
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
int evt
StopperSelection(const fhicl::ParameterSet &pset)
const double j
Definition: BetheBloch.cxx:29
art::ServiceHandle< geo::Geometry > fGeom
OStream cout
Definition: OStream.cxx:6
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
assert(nhit_max >=nhit_nbins)
art::ServiceHandle< geo::LiveGeometry > fLiveGeom
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.
enum BeamMode string