RecoValidationTutorial_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file
3 // \brief Simple demonstration of grabbing validation from art::Event
4 // \author Fernanda Psihas - psihas@fnal.gov
5 //
6 ////////////////////////////////////////////////////////////////////////
7 
8 // ART includes
16 #include "fhiclcpp/ParameterSet.h"
18 
19 // NOvA includes
20 #include "MCCheater/BackTracker.h"
21 #include "RecoBaseHit/CellHit.h"
22 #include "RecoBaseHit/RecoHit.h"
23 #include "RecoBase/Cluster.h"
24 #include "RecoBase/Track.h"
25 #include "RecoBase/Vertex.h"
29 
30 #include "TH1F.h"
31 #include "TVector3.h"
32 
33 namespace tut
34 {
36  {
37 
38  public:
41 
42  // Required functions.
43  void analyze(const art::Event& evt);
44 
45  // Selected optional functions.
46  void reconfigure(const fhicl::ParameterSet& p);
47 
48  void beginJob();
49 
50  protected:
51 
52  std::string fCellHitLabel; ///< module label for CellHit
53  std::string fClusterLabel; ///< module label for clusters
54  std::string fVertexLabel; ///< module label for vertices
55 
57 
58  TH1* fVertexX;
59  TH1* fVertexY;
60  TH1* fVertexZ;
61 
62  };
63 }
64 
65 namespace tut {
66 
67  //.......................................................................
69  : EDAnalyzer(pset)
70  {
71  reconfigure(pset);
72  }
73 
74  //.......................................................................
76  {
77  }
78 
79  //.......................................................................
81  {
82  fClusterLabel = pset.get<std::string>("ClusterLabel");
83  fVertexLabel = pset.get<std::string>("VertexLabel");
84  }
85 
86  //.......................................................................
88  {
90  fVertexX = tfs->make<TH1F>("vertexXdiff", ";Reco - True Vertex X (cm)",80, -200, 200);
91  fVertexY = tfs->make<TH1F>("vertexYdiff", ";Reco - True Vertex Y (cm)",80, -200, 200);
92  fVertexZ = tfs->make<TH1F>("vertexZdiff", ";Reco - True Vertex Z (cm)",80, -200, 200);
93  }
94 
95  //-----------------------------------------------------------------------
97  {
98 
99  // ---------- Get data products ----------
100 
101  // Get all Slices in the event
103  evt.getByLabel(fClusterLabel, sliceHandle);
104 
106  for ( unsigned int i = 0; i < sliceHandle->size(); ++i ) {
107  art::Ptr<rb::Cluster> sliceptr(sliceHandle, i);
108  slices.push_back(sliceptr);
109  }
110 
111 
112  // Get vertices associated with the slices
113  art::FindManyP<rb::Vertex> fmVert(sliceHandle, evt, fVertexLabel);
114 
115 
116  // Get all truths related to this slice
117  std::vector<std::vector<cheat::NeutrinoEffPur>> slTruth = bt->SlicesToMCTruthsTable(slices);
118  // Get all neutrinos from the event
119  std::vector<cheat::NeutrinoWithIndex> nus = bt->allMCTruth();
120  std::vector<std::vector<cheat::NeutrinoEffPur>> sEffPur = bt->SlicesToMCTruthsTable(slices);
121  // and order them by purity
122  std::vector<int> sliceNuId = bt->SliceToOrderedNuIds(nus, slTruth,
125 
126 
127  // ---------- Loop Over Slices ----------
128  for ( unsigned int iSlice = 0; iSlice < slices.size(); ++iSlice ){
129 
130  // Simple example. Print number of cells in each slice.
131  const art::Ptr<rb::Cluster> slice = slices[iSlice];
132  unsigned int nCells = slice->NCell();
133 
134  mf::LogInfo("RecoValidationTutorial") << " Slice # " << iSlice
135  << " has " << nCells
136  << " cells." << std::endl;
137 
138  if ( slice->IsNoise() ) continue;
139 
140 
141  // ---------- Reco Vertex Info ----------
142 
143  // Always check the assns are ok
144  if( !fmVert.isValid() ) continue;
145 
146  const std::vector< art::Ptr<rb::Vertex> > vtx = fmVert.at(iSlice);
147  if ( vtx.size() != 1) continue;
148 
149  TVector3 vxyz = vtx[0]->GetXYZ();
150  mf::LogInfo("RecoValidationTutorial") << " It has Reco Vertex at "
151  << "( "<< vxyz.X()
152  << ", "<< vxyz.Y()
153  << ", "<< vxyz.Z()
154  << " ) " <<std::endl;
155 
156 
157  // ---------- True Vertex Info ----------
158 
159  // Check there are related nus
160  if ( sliceNuId[iSlice] == -1 ) return;
161 
162  // Get this slice's highest eff. neutrino
163  art::Ptr<simb::MCTruth> truth = slTruth[iSlice][sliceNuId[iSlice]].neutrinoInt;
164  simb::MCNeutrino thisNu = truth->GetNeutrino();
165 
166  //now get vertex info and make comparisons
167  fVertexX->Fill( vxyz.X() - thisNu.Nu().Vx() );
168  fVertexY->Fill( vxyz.Y() - thisNu.Nu().Vy() );
169  fVertexZ->Fill( vxyz.Z() - thisNu.Nu().Vz() );
170 
171  } //iSlice
172  }// analyze
173  //-----------------------------------------------------------------------
174 
175 } // end namespce tut
176 
177 // Magic. Do not touch this
back track the reconstruction to the simulation
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:76
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
const char * p
Definition: xmltok.h:285
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
std::vector< NeutrinoWithIndex > allMCTruth() const
art::ServiceHandle< cheat::BackTracker > bt
DEFINE_ART_MODULE(TestTMapFile)
Particle class.
std::vector< std::vector< cheat::NeutrinoEffPur > > SlicesToMCTruthsTable(const std::vector< const rb::Cluster * > &sliceList) const
Given ALL the slices in an event, including the noise slice, returns a vector of vector of structures...
std::string fCellHitLabel
module label for CellHit
std::vector< int > SliceToOrderedNuIds(const std::vector< cheat::NeutrinoWithIndex > &nusWithIdx, const std::vector< std::vector< cheat::NeutrinoEffPur >> &slTruthTable, std::function< double(const cheat::NeutrinoEffPur &)> slMetric, std::function< double(const cheat::NeutrinoEffPur &)> nuMetric) const
Given a vector of indexed neutrino interaction and a vector of slice truth tables, returns a vector of neutrino interaction indices ordered by best match to the corresponding slice. Here, best match is determined according to the given cheat::NeutrinoEffPur functions for the slice and the nu.
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
std::string fClusterLabel
module label for clusters
T get(std::string const &key) const
Definition: ParameterSet.h:231
int evt
TODO.
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
Vertex location in position and time.
size_type size() const
Definition: PtrVector.h:308
RecoValidationTutorial(const fhicl::ParameterSet &p)
double Vx(const int i=0) const
Definition: MCParticle.h:220
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
std::string fVertexLabel
module label for vertices
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
double Vz(const int i=0) const
Definition: MCParticle.h:222
void reconfigure(const fhicl::ParameterSet &p)
double EffMetric(const cheat::NeutrinoEffPur &ep)
Function for NeutrinoEffPur&#39;s nu interaction to slice efficiency.
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
Event generator information.
Definition: MCNeutrino.h:18
double Vy(const int i=0) const
Definition: MCParticle.h:221