BPFEnergyAna_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file BPFEnergyAna_module.cc
3 // \brief Module to create the TTree and histos used to examine and
4 // validate the performance of the BPFEnergyEstimator module.
5 // \version
6 // \author $Author: mbaird42 $
7 ////////////////////////////////////////////////////////////////////////
8 
9 // C/C++ includes
10 #include <cmath>
11 #include <iostream>
12 #include <vector>
13 
14 // ROOT includes
15 #include "TH1F.h"
16 #include "TH2F.h"
17 #include "TTree.h"
18 
19 // Framework includes
29 #include "fhiclcpp/ParameterSet.h"
31 
32 // NOvASoft includes
33 #include "RecoBase/Cluster.h"
35 
36 #include "Utilities/AssociationUtil.h"
37 
38 #include "MCCheater/BackTracker.h"
41 
42 
43 
44 // BPFEnergyAna header
45 namespace bpfit {
46  class BPFEnergyAna : public art::EDAnalyzer {
47  public:
48  explicit BPFEnergyAna(fhicl::ParameterSet const& pset);
49  ~BPFEnergyAna();
50 
51  void analyze(const art::Event& evt);
52  void reconfigure(const fhicl::ParameterSet& pset);
53  void beginJob();
54  void resetVars();
55  void setTruthVars(const std::vector<cheat::NeutrinoEffPur>& nuEP);
56 
57  private:
59  std::string fEnergyLabel; ///< label for module that made the energy objects
60  bool fIsMC; ///< are you running over MC?
61 
62  TTree *fBPFEnergyTree; ///< TTree with slice level info (for energy estimator training)
63 
64 
65 
66  //
67  // Variables to go into the Energy TTree.
68  //
69 
70  // variables from the BPFEnergy object
71  float fEreco1; ///< first reconstructed event energy
72  float fEres1; ///< first estimated event energy resolution
73  float fEreco2; ///< second reconstructed event energy
74  float fEres2; ///< second estimated event energy resolution
75  float fEreco3; ///< third reconstructed event energy
76  float fEres3; ///< third estimated event energy resolution
77  float fPmu; ///< reconstructed muon momentum from the track selected to be the muon
78  float fDirZmu; ///< track.dirZ() from the track selected to be the muon
79  float fN3Dprongs; ///< number of fuzzyk 3D prongs
80  float fE3Dprongs; ///< sum of hit energies in 3D prongs (excluding the "muon" prong) - this does NOT account for dead material
81  float fEremain; ///< sum of remaining slice energy not in 3D prongs - this does NOT account for dead material
82  float fSumPE; ///< sum of PE from all hits not on the muon track
83  float fEventID; ///< best BPF muon PID value for the event
84 
85  // MCTruth and slice variables used for cutting etc.
86  float fE; ///< true energy
87  float fMinX; ///< minimum hit X value in the slice
88  float fMinY; ///< minimum hit Y value in the slice
89  float fMinZ; ///< minimum hit Z value in the slice
90  float fMaxX; ///< maximum hit X value in the slice
91  float fMaxY; ///< maximum hit Y value in the slice
92  float fMaxZ; ///< maximum hit Z value in the slice
93  int fNHitX; ///< number of X-view hits in the slice
94  int fNHitY; ///< number of Y-view hits in the slice
95  int fCCNC; ///< is this CC or NC (from truth)
96  int fIntType; ///< neutrino interaction type (from truth)
97  int fnuPDG; ///< PDG code for the best matched nu (from truth)
98 
99  };
100 }
101 
102 
103 
104 namespace bpfit
105 {
106  //.......................................................................
108  : EDAnalyzer(pset),
109  fSlicerToken(consumes<std::vector<rb::Cluster>>(pset.get<std::string>("SlicerLabel")))
110  {
111  this->reconfigure(pset);
112  }
113 
114  //......................................................................
116  {
117  //======================================================================
118  // Clean up any memory allocated by your module... ...if you dare...
119  //======================================================================
120  }
121 
122  //......................................................................
124  {
125  fEnergyLabel = pset.get<std::string>("EnergyLabel");
126  fIsMC = pset.get<bool> ("IsMC");
127  }
128 
129  //......................................................................
131  {
133 
134  //
135  // Book the Energy TTree
136  //
137  fBPFEnergyTree = tfs->make<TTree>("BPFEnergyTree","Tree containing BPFEnergy object info");
138 
139  // attach the variables:
140  fBPFEnergyTree->Branch("Ereco1",&fEreco1,"Ereco1/F");
141  fBPFEnergyTree->Branch("Eres1",&fEres1,"Eres1/F");
142  fBPFEnergyTree->Branch("Ereco2",&fEreco2,"Ereco2/F");
143  fBPFEnergyTree->Branch("Eres2",&fEres2,"Eres2/F");
144  fBPFEnergyTree->Branch("Ereco3",&fEreco3,"Ereco3/F");
145  fBPFEnergyTree->Branch("Eres3",&fEres3,"Eres3/F");
146  fBPFEnergyTree->Branch("Pmu",&fPmu,"Pmu/F");
147  fBPFEnergyTree->Branch("DirZmu",&fDirZmu,"DirZmu/F");
148  fBPFEnergyTree->Branch("N3Dprongs",&fN3Dprongs,"N3Dprongs/F");
149  fBPFEnergyTree->Branch("E3Dprongs",&fE3Dprongs,"E3Dprongs/F");
150  fBPFEnergyTree->Branch("Eremain",&fEremain,"Eremain/F");
151  fBPFEnergyTree->Branch("SumPE",&fSumPE,"SumPE/F");
152  fBPFEnergyTree->Branch("EventID",&fEventID,"EventID/F");
153 
154  fBPFEnergyTree->Branch("E",&fE,"E/F");
155  fBPFEnergyTree->Branch("MinX",&fMinX,"MinX/F");
156  fBPFEnergyTree->Branch("MinY",&fMinY,"MinY/F");
157  fBPFEnergyTree->Branch("MinZ",&fMinZ,"MinZ/F");
158  fBPFEnergyTree->Branch("MaxX",&fMaxX,"MaxX/F");
159  fBPFEnergyTree->Branch("MaxY",&fMaxY,"MaxY/F");
160  fBPFEnergyTree->Branch("MaxZ",&fMaxZ,"MaxZ/F");
161  fBPFEnergyTree->Branch("NHitX",&fNHitX,"NHitX/I");
162  fBPFEnergyTree->Branch("NHitY",&fNHitY,"NHitY/I");
163  fBPFEnergyTree->Branch("CCNC",&fCCNC,"CCNC/I");
164  fBPFEnergyTree->Branch("IntType",&fIntType,"IntType/I");
165  fBPFEnergyTree->Branch("nuPDG",&fnuPDG,"nuPDG/I");
166 
167  }
168 
169  //......................................................................
171  {
172  fEreco1 = -5.0;
173  fEres1 = -5.0;
174  fEreco2 = -5.0;
175  fEres2 = -5.0;
176  fEreco3 = -5.0;
177  fEres3 = -5.0;
178  fPmu = -5.0;
179  fDirZmu = -5.0;
180  fN3Dprongs = -5.0;
181  fE3Dprongs = -5.0;
182  fEremain = -5.0;
183  fSumPE = -5.0;
184  fEventID = -5.0;
185 
186  fE = -5.0;
187  fMinX = 1.0e9;
188  fMinY = 1.0e9;
189  fMinZ = 1.0e9;
190  fMaxX = -1.0e9;
191  fMaxY = -1.0e9;
192  fMaxZ = -1.0e9;
193  fNHitX = -5;
194  fNHitY = -5;
195  fCCNC = -5;
196  fIntType = -5;
197  fnuPDG = 0;
198  }
199 
200  //......................................................................
201  void BPFEnergyAna::setTruthVars(const std::vector<cheat::NeutrinoEffPur>& nuEP)
202  {
203  //
204  // This function sets any variables for the TTree that come from truth.
205  // Currently, the most pure neutrino is assumed to be the best
206  // match.
207  //
208 
209  float mostPure = 0.0;
210  unsigned int match = 0;
211  for(unsigned int i = 0; i < nuEP.size(); ++i) {
212  if(nuEP[i].purity > mostPure && nuEP[i].neutrinoInt->NeutrinoSet()) {
213  mostPure = nuEP[i].purity;
214  match = i;
215  }
216  }
217 
218  if(mostPure > 0 && nuEP[match].neutrinoInt->NeutrinoSet()) {
219  simb::MCNeutrino nu = nuEP[match].neutrinoInt->GetNeutrino();
220  simb::MCParticle Nu = nu.Nu();
221  fE = Nu.E();
222  fCCNC = nu.CCNC();
223  fIntType = nu.InteractionType();
224  fnuPDG = Nu.PdgCode();
225  }
226 
227  }
228 
229  //......................................................................
231  {
232 
233  // get the slices
235  evt.getByToken(fSlicerToken, slices);
236 
237  // create the FindMany object for getting the energy objects
238  art::FindManyP<bpfit::BPFEnergy> energyFMP(slices, evt, fEnergyLabel);
239 
240  // Make the required service handles...
242 
243  // Put the slices into a PtrVector and get the list neutrinos matched
244  // to each slice from backtracker.
245  art::PtrVector<rb::Cluster> ARTslices;
246  for(unsigned int i = 0; i < slices->size(); ++i) {
247  ARTslices.push_back(art::Ptr<rb::Cluster>(slices,i));
248  }
249  std::vector< std::vector< cheat::NeutrinoEffPur > > ALLnuEffPur;
250  if(fIsMC) {
251  ALLnuEffPur = BT->SlicesToMCTruthsTable(ARTslices);
252  }
253 
254 
255  // loop over all slices and get the fuzzyk 3D prongs
256  for(unsigned int islice = 0; islice < slices->size(); ++islice) {
257 
258  // As usual, skip the noise slice...
259  if((*slices)[islice].IsNoise()) continue;
260 
261  // NOTE: Currently, the fundamental object of analysis is the slice.
262  // If in the future we change things so that the fundamental object
263  // of analysis becomes the vertex, then we will have to get prongs
264  // associated with the vertex instead.
265 
266  // Reset all TTree variables.
267  this->resetVars();
268 
269  // Set slice level and truth variables for the TTrees.
270  fMinX = (*slices)[islice].MinX();
271  fMinY = (*slices)[islice].MinY();
272  fMinZ = (*slices)[islice].MinZ();
273  fMaxX = (*slices)[islice].MaxX();
274  fMaxY = (*slices)[islice].MaxY();
275  fMaxZ = (*slices)[islice].MaxZ();
276  fNHitX = (*slices)[islice].NXCell();
277  fNHitY = (*slices)[islice].NYCell();
278  if(fIsMC) {
279  this->setTruthVars(ALLnuEffPur[islice]);
280  }
281 
282  // get the energy object associated with this slice
283  std::vector< art::Ptr< bpfit::BPFEnergy > > energies;
284  if(energyFMP.isValid()) energies = energyFMP.at(islice);
285 
286  // There can be only one!!!
287  // Something went wrong if there is more than one energy object
288  // per slice.
289  if(energies.size() != 1) continue;
290 
291  // Set remaining branch values and fill the TTree.
292  fEreco1 = energies[0]->E1();
293  fEres1 = energies[0]->Eres1();
294  fEreco2 = energies[0]->E2();
295  fEres2 = energies[0]->Eres2();
296  fEreco3 = energies[0]->E3();
297  fEres3 = energies[0]->Eres3();
298  fPmu = energies[0]->PMuon();
299  fDirZmu = energies[0]->DirZMuon();
300  fN3Dprongs = energies[0]->N3DProngs();
301  fE3Dprongs = energies[0]->EFuzzyK3D();
302  fEremain = energies[0]->ERemain();
303  fSumPE = energies[0]->SumPE();
304  fEventID = energies[0]->EventID();
305 
306  fBPFEnergyTree->Fill();
307 
308  } // end loop over slices (islice)
309 
310  return;
311 
312  }
313 
314 } // end namespace bpfit
315 ////////////////////////////////////////////////////////////////////////
316 
317 
318 
320 
321 namespace bpfit
322 {
324 }
double E(const int i=0) const
Definition: MCParticle.h:232
back track the reconstruction to the simulation
int PdgCode() const
Definition: MCParticle.h:211
int CCNC() const
Definition: MCNeutrino.h:148
float fEreco2
second reconstructed event energy
float fMaxY
maximum hit Y value in the slice
float fE3Dprongs
sum of hit energies in 3D prongs (excluding the "muon" prong) - this does NOT account for dead materi...
int fCCNC
is this CC or NC (from truth)
float fN3Dprongs
number of fuzzyk 3D prongs
int fNHitX
number of X-view hits in the slice
float fEres3
third estimated event energy resolution
float fEreco3
third reconstructed event energy
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
float fEremain
sum of remaining slice energy not in 3D prongs - this does NOT account for dead material ...
void analyze(const art::Event &evt)
DEFINE_ART_MODULE(TestTMapFile)
float fEventID
best BPF muon PID value for the event
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 fEnergyLabel
label for module that made the energy objects
float fEres2
second estimated event energy resolution
bool fIsMC
are you running over MC?
float fE
true energy
int InteractionType() const
Definition: MCNeutrino.h:150
float fSumPE
sum of PE from all hits not on the muon track
void setTruthVars(const std::vector< cheat::NeutrinoEffPur > &nuEP)
float fPmu
reconstructed muon momentum from the track selected to be the muon
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
float fMaxX
maximum hit X value in the slice
float fEres1
first estimated event energy resolution
T get(std::string const &key) const
Definition: ParameterSet.h:231
int evt
float fMinX
minimum hit X value in the slice
int fIntType
neutrino interaction type (from truth)
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
Perform a "2 point" Hough transform on a collection of hits.
float fMaxZ
maximum hit Z value in the slice
float fDirZmu
track.dirZ() from the track selected to be the muon
BPFEnergyAna(fhicl::ParameterSet const &pset)
float fMinY
minimum hit Y value in the slice
float fEreco1
first reconstructed event energy
void reconfigure(const fhicl::ParameterSet &pset)
T * make(ARGS...args) const
float fMinZ
minimum hit Z value in the slice
int fnuPDG
PDG code for the best matched nu (from truth)
const art::ProductToken< std::vector< rb::Cluster > > fSlicerToken
TTree * fBPFEnergyTree
TTree with slice level info (for energy estimator training)
bool getByToken(ProductToken< PROD > const &token, Handle< PROD > &result) const
Definition: DataViewImpl.h:387
Event generator information.
Definition: MCNeutrino.h:18
ProductToken< T > consumes(InputTag const &)
int fNHitY
number of Y-view hits in the slice
enum BeamMode string