ElasticArmsValidate_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file ElasticArmsValidate.cxx
3 // \brief Analysis module for performance of FuzzyK
4 // \version $Id: ElasticArmsValidate.cxx,v 1.2 2012-12-12 15:55:23 edniner Exp $
5 // \author $Author: edniner $
6 ////////////////////////////////////////////////////////////////////////
7 
8 // C/C++ includes
9 #include <cmath>
10 #include <iostream>
11 #include <map>
12 #include <vector>
13 #include <string>
14 
15 // ROOT includes
16 #include "TTree.h"
17 #include "TH1F.h"
18 #include "TH2F.h"
19 
20 // Framework includes
30 #include "fhiclcpp/ParameterSet.h"
32 
33 // NOvASoft includes
37 #include "Utilities/func/MathUtil.h"
38 #include "MCCheater/BackTracker.h"
39 #include "Geometry/Geometry.h"
41 #include "GeometryObjects/Geo.h"
42 #include "RecoBase/Cluster.h"
43 #include "RecoBase/Vertex.h"
44 
45 // Required for CAFAna Cuts
47 #include "CAFAna/Cuts/Cuts.h"
48 #include "CAFAna/Cuts/TruthCuts.h"
49 #include "CAFAna/Cuts/TimingCuts.h"
50 #include "CAFAna/Vars/Vars.h"
51 
54 
55 namespace earms {
57  public:
58  explicit ElasticArmsValidate(fhicl::ParameterSet const& pset);
60 
61  void analyze(const art::Event& evt);
62  void reconfigure(const fhicl::ParameterSet& pset);
63  void beginJob();
64  void ClearVectors();
65 
66  private:
67 
68  /// MC information
69  std::string fSliceLabel; ///< label for the process that made the slices
70  std::string fCellHitLabel; ///< label for the process that made the cell hits
71  std::string fVertexLabel; ///< label for the process that made the vertices
72  int fCompMCNeutrino; ///< don't know what this one does... :(
73 
74  std::string fCAFLabel; ///< label for the process that made the standard records
75  bool fApplyCAFCuts; ///< should we apply the CAF-level cuts?
76  int fCutType; ///< what cuts to apply (see CAFCutter.h)
77 
78  recovalid::CAFCutter fCAFCutter; ///< the CAFCutter helper class
79 
80  TTree* fOutTree;
81  //event summary info
82  int run;
83  int subrun;
84  int event;
85 
86  //vertex position, reconstructed
87  float rVx;
88  float rVy;
89  float rVz;
90 
91  //vertex histograms
92  TH1F* fHistoRVX;//reco x position
93  TH1F* fHistoRVY;//reco y position
94  TH1F* fHistoRVZ;//reco z position
95  TH1F* fHistoDVX;//delta x
96  TH1F* fHistoDVY;//delta y
97  TH1F* fHistoDVZ;//delta z
98  TH1F* fHistoDVT;//total diff between reco and true vertex
99 
100  //for cosmics make histograms for endpoint as well
101  //TH1F* fHistoDEX;//delta x
102  //TH1F* fHistoDEY;//delta y
103  //TH1F* fHistoDEZ;//delta z
104  //TH1F* fHistoDET;//total diff between reco and true vertex
105 
106  //neutrino information
107  float Vx;//start x
108  float Vy;//start y
109  float Vz;//start z
110  float NuEff;//efficiency of neutrino for slice
111  float NuPur;//purity of neutrino for slice
112  float NuEVSl;//visible energy of neutrino in slice
113  float NuEVTot;//visible energy of neutrino in all slices
114  int NuPdg;//neutrino pdg
115  int ccnc;//is cc or nc
116  int mode;//interaction mode
117  float NuE;//total energy of neutrino
118 
119  //for cosmics also store endpoint since elastic arms can get confused
120  float Ex;//end x
121  float Ey;//end y
122  float Ez;//end z
123  };
124 }
125 
126 namespace earms
127 {
128  //.......................................................................
130  : EDAnalyzer(pset)
131  {
132  this->reconfigure(pset);
133  }
134 
135  //......................................................................
137  {
138  //======================================================================
139  // Clean up any memory allocated by your module
140  //======================================================================
141  }
142 
143  //......................................................................
145  {
146 
147  fSliceLabel = pset.get< std::string >("SliceLabel");
148  fCellHitLabel = pset.get< std::string >("CellHitLabel");
149  fVertexLabel = pset.get< std::string >("VertexLabel");
150  fCompMCNeutrino = pset.get< int >("CompMCNeutrino");
151 
152  fCAFLabel = pset.get< std::string >("CAFLabel");
153  fApplyCAFCuts = pset.get< bool > ("ApplyCAFCuts");
154  fCutType = pset.get< int > ("CutType");
155  }
156 
157  //......................................................................
159  {
160 
162 
163  // define the tree
164  fOutTree = tfs->make<TTree>("ElasticArmsValidate","Elastic Arms Validation");
165 
166  //store event summary information
167  fOutTree->Branch("run", &run, "run/I");
168  fOutTree->Branch("subrun", &subrun, "subrun/I");
169  fOutTree->Branch("event", &event, "event/I");
170 
171  //store vertex info
172  fOutTree->Branch("rVx", &rVx);
173  fOutTree->Branch("rVy", &rVy);
174  fOutTree->Branch("rVz", &rVz);
175 
176  //vertex histograms
177  fHistoRVX = tfs->make<TH1F>("VertexRecoX","Reconstructed X vertex coordinate;vertex x (cm)",160,-800.0,800.0);
178  fHistoRVY = tfs->make<TH1F>("VertexRecoY","Reconstructed Y vertex coordinate;vertex y (cm)",160,-800.0,800.0);
179  fHistoRVZ = tfs->make<TH1F>("VertexRecoZ","Reconstructed Z vertex coordinate;vertex z (cm)",320,0.0,6400.0);
180  fHistoDVX = tfs->make<TH1F>("VertexDeltaX","Vertex True X - Reco X;True Vertex X - Reco Vertex X (cm)",201,-100.5,100.5);
181  fHistoDVY = tfs->make<TH1F>("VertexDeltaY","Vertex True Y - Reco Y;True Vertex Y - Reco Vertex Y (cm)",201,-100.5,100.5);
182  fHistoDVZ = tfs->make<TH1F>("VertexDeltaZ","Vertex True Z - Reco Z;True Vertex Z - Reco Vertex Z (cm)",201,-100.5,100.5);
183  fHistoDVT = tfs->make<TH1F>("VertexDelta3D","Distance from Reco to True vertex;|True Vertex - Reco Vertex| (cm)",150,0.0,300.0);
184  //fHistoDEX = tfs->make<TH1F>("fHistoDEX","Endpoint True X - Reco X;(cm)",201,-100.5,100.5);
185  //fHistoDEY = tfs->make<TH1F>("fHistoDEY","Endpoint True Y - Reco Y;(cm)",201,-100.5,100.5);
186  //fHistoDEZ = tfs->make<TH1F>("fHistoDEZ","Endpoint True Z - Reco Z;(cm)",201,-100.5,100.5);
187  //fHistoDET = tfs->make<TH1F>("fHistoDET","Distance from Reco to Endpoint;(cm)",150,0.0,300.0);
188 
189  //store neutrino info
190  fOutTree->Branch("Vx", &Vx);//neutrino position info
191  fOutTree->Branch("Vy", &Vy);
192  fOutTree->Branch("Vz", &Vz);
193  fOutTree->Branch("Ex", &Ex);
194  fOutTree->Branch("Ey", &Ey);
195  fOutTree->Branch("Ez", &Ez);
196  fOutTree->Branch("NuEff", &NuEff);//info on quality of neutrino in slice
197  fOutTree->Branch("NuPur", &NuPur);
198  fOutTree->Branch("NuEVSl", &NuEVSl);
199  fOutTree->Branch("NuEVTot", &NuEVTot);
200  fOutTree->Branch("NuPdg", &NuPdg);
201  fOutTree->Branch("ccnc", &ccnc);
202  fOutTree->Branch("mode", &mode);
203  fOutTree->Branch("NuE", &NuE);
204 
205  }
206 
207  //......................................................................
209  {
212 
214  evt.getByLabel(fSliceLabel, sHandle);
215  if(sHandle.failedToGet()){
216  std::cout<<"No slices found, skipping event" << std::endl;
217  }
219  for(unsigned int i=0; i<sHandle->size(); ++i){
220  slices.push_back(art::Ptr<rb::Cluster>(sHandle, i));
221  }
222 
223  // Get the FMP for the Standard Records
224  art::FindManyP<caf::StandardRecord> recordFMP(slices, evt, fCAFLabel);
225 
226  //get associations between slices and vertices
227  art::FindManyP<rb::Vertex> fmVert(sHandle, evt, fVertexLabel);
228 
229  //if these are events with neutrino information to compare too, get it
230  std::vector<int> bestNuId;
231  std::vector<cheat::NeutrinoWithIndex> nus;
232  std::vector<std::vector<cheat::NeutrinoEffPur>> sEffPur;
233  if(fCompMCNeutrino == 1){
234  nus = bt->allMCTruth();
235  //std::cout<<allNus.size()<<std::endl;
236  sEffPur = bt->SlicesToMCTruthsTable(slices);
237  //std::cout<<"Done matching"<<std::endl;
238 
239  bestNuId = bt->SliceToOrderedNuIdsByEff(nus, sEffPur);
240  }
241 
242 
243  //loop over slices
244  for( unsigned int i=0; i<slices.size(); ++i){
245  this->ClearVectors();
246 
247  // Apply the CAF-level cuts
248  bool pass = true;
249  if( fApplyCAFCuts ) {
250  // get record associated with the slice
251  std::vector< art::Ptr<caf::StandardRecord> > records = recordFMP.at(i);
252  if (records.size() > 0 && recordFMP.isValid()) {
253  pass = fCAFCutter.passCuts(fCutType, records[0].get());
254  }
255  else { pass = false; }
256  }
257 
258  if(!pass) continue;
259 
260 
261 
262  //skip noise slice
263  if (slices[i]->IsNoise()) continue;
264  //if we are comparing to truth, make sure there is a neutrino with this slice
265  //make sure there is a vertex associated with this slice
266  if(!(fmVert.isValid())) continue;
267 
268  std::vector<art::Ptr<rb::Vertex>> v = fmVert.at(i);
269  //continue if there is no vertex or more then 1 vertex
270  if (v.size() != 1) continue;
271 
272 
273  //get reco vertex position
274  rVx = v[0]->GetX();
275  rVy = v[0]->GetY();
276  rVz = v[0]->GetZ();
277  //fill position histograms
278  fHistoRVX->Fill(rVx);
279  fHistoRVY->Fill(rVy);
280  fHistoRVZ->Fill(rVz);
281 
282  //if we do not want to compare to truth, dump into to tree and move to next slice
283  if(fCompMCNeutrino != 1){
284  fOutTree->Fill();
285  continue;
286  }
287  else if (bestNuId[i] == -1){
288  fOutTree->Fill();
289  continue;
290  }
291 
292  //now access true neutrino info for further validation
293  art::Ptr<simb::MCTruth> nuTru = sEffPur[i][bestNuId[i]].neutrinoInt;
294 
295  NuEff = sEffPur[i][0].efficiency;
296  NuPur = sEffPur[i][0].purity;
297  NuEVSl = sEffPur[i][0].energySlice;
298  NuEVTot = sEffPur[i][0].energyTotal;
299 
300  //if we have a neutrino in this mc truth, proceed
301  if (nuTru->NeutrinoSet()){
302  NuPdg = nuTru->GetNeutrino().Nu().PdgCode();
303  ccnc = nuTru->GetNeutrino().CCNC();
304  mode = nuTru->GetNeutrino().Mode();
305  NuE = nuTru->GetNeutrino().Nu().Momentum().E();
306 
307  //now get vertex info and make comparisons
308  Vx = nuTru->GetNeutrino().Nu().Vx();
309  Vy = nuTru->GetNeutrino().Nu().Vy();
310  Vz = nuTru->GetNeutrino().Nu().Vz();
311  fHistoDVX->Fill(Vx-rVx);
312  fHistoDVY->Fill(Vy-rVy);
313  fHistoDVZ->Fill(Vz-rVz);
314  float distdiff = util::pythag(Vx-rVx,Vy-rVy,Vz-rVz);
315  fHistoDVT->Fill(distdiff);
316  }
317  else {
318  TVector3 enter;
319  TVector3 exit;
320  TVector3 stop;
321  //if this is a cosmic ray then we need to find the start and end points inside the detector
322  std::vector<const sim::Particle*> particles = bt->MCTruthToParticles(nuTru);
323  const simb::MCTrajectory& traj = particles[0]->Trajectory();
324 
325 
326  bool entered = 0;
327  bool exited = 0;
328 
329  for(size_t pt = 0; pt < traj.size(); ++pt)
330  {
331  TVector3 point = traj.Position(pt).Vect();
332 
333  // Check to see if trajectory point is in the detector
334  bool inDet = fabs(point.X()) <= geom->DetHalfWidth() // Check X inside detector
335  && fabs(point.Y()) <= geom->DetHalfHeight() /// Then check Y
336  && point.Z() >= 0 && point.Z() <= geom->DetLength(); // then check Z
337 
338  if(inDet && !entered)
339  {
340  // then we just entered the detector
341  entered = true;
342 
343  if(pt > 0)
344  {
345  // Put xyz from last point into an array for geo function call
346  TVector3 lastPoint = traj.Position(pt - 1).Vect();
347  double pointArray[3] = {point.X(), point.Y(), point.Z()};
348 
349  // Put xyz from direction vector into an array.
350  TVector3 dir = (lastPoint - point).Unit();
351  double dirArray[3] = {dir.X(), dir.Y(), dir.Z()};
352 
353  double enterArray[3] = {0,0,0};
354 
355  geo::ProjectToBoxEdge(pointArray, dirArray,
356  - geom->DetHalfWidth() , geom->DetHalfWidth(),
357  - geom->DetHalfHeight() , geom->DetHalfHeight(),
358  0 , geom->DetLength(),
359  enterArray);
360 
361  enter.SetXYZ(enterArray[0], enterArray[1], enterArray[2]);
362 
363  }
364  else
365  {// if this is the case, we started inside the detector.
366  enter = point;
367  }
368  }
369 
370  if(!inDet && entered && !exited){
371  // then we just exited the detector
372  exited= true;
373 
374  // Put xyz from last point into an array for geo function call
375  TVector3 lastPoint = traj.Position(pt - 1).Vect();
376  double lastPointArray[3] = {lastPoint.X(), lastPoint.Y(), lastPoint.Z()};
377 
378  // Put xyz from direction vector into an array.
379  TVector3 dir = (point - lastPoint).Unit();
380  double dirArray[3] = {dir.X(), dir.Y(), dir.Z()};
381 
382  double exitArray[3] = {0,0,0};
383 
384  geo::ProjectToBoxEdge(lastPointArray, dirArray,
385  - geom->DetHalfWidth() , geom->DetHalfWidth(),
386  - geom->DetHalfHeight() , geom->DetHalfHeight(),
387  0 , geom->DetLength(),
388  exitArray);
389 
390  exit.SetXYZ(exitArray[0], exitArray[1], exitArray[2]);
391 
392  }
393 
394 
395  }
396  stop = traj.Position(traj.size() - 1).Vect();
397 
398  // make sure the particle actually entered the detector
399  if(!entered){
400  /// if the primary never entered the detector, set the coordinates to -5*10^9
401 
402  enter.SetXYZ(-5e9, -5e9, -5e9);
403  exit.SetXYZ(-5e9, -5e9, -5e9);
404 
405  }
406 
407  /// If the particle entered and did not exit, set the exit position to stop
408  if(entered && !exited){
409  exit = stop;
410  }
411 
412  Vx = enter.X();
413  Vy = enter.Y();
414  Vz = enter.Z();
415  Ex = exit.X();
416  Ey = exit.Y();
417  Ez = exit.Z();
418 
419  fHistoDVX->Fill(Vx-rVx);
420  fHistoDVY->Fill(Vy-rVy);
421  fHistoDVZ->Fill(Vz-rVz);
422  float distdiff = util::pythag(Vx-rVx,Vy-rVy,Vz-rVz);
423  fHistoDVT->Fill(distdiff);
424 
425  //fHistoDEX->Fill(Ex-rVx);
426  //fHistoDEY->Fill(Ey-rVy);
427  //fHistoDEZ->Fill(Ez-rVz);
428  //float distdiff2 = util::pythag(Ex-rVx,Ey-rVy,Ez-rVz);
429  //fHistoDET->Fill(distdiff2);
430 
431 
432  }
433 
434  //fill tree
435  fOutTree->Fill();
436  }//end loop over slices
437 
438 
439 
440  }
441 
442  //.........................................................................
444  {
445  ccnc = -1;
446  NuE = -1.0;
447  NuEff = -1.0;
448  NuPur = -1.0;
449  NuEVSl = -1.0;
450  NuEVTot = -1.0;
451  Vx = -999.0;
452  Vy = -999.0;
453  Vz = -999.0;
454  NuPdg = 0;
455  mode = -1;
456  }
457 
458 } // end namespace
459 ////////////////////////////////////////////////////////////////////////
460 namespace earms
461 {
463 }
464 ////////////////////////////////////////////////////////////////////////
std::string fSliceLabel
MC information.
back track the reconstruction to the simulation
int PdgCode() const
Definition: MCParticle.h:211
int CCNC() const
Definition: MCNeutrino.h:148
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:76
bool fApplyCAFCuts
should we apply the CAF-level cuts?
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
void analyze(const art::Event &evt)
std::vector< int > SliceToOrderedNuIdsByEff(const std::vector< cheat::NeutrinoWithIndex > &nusWithIdx, const std::vector< std::vector< cheat::NeutrinoEffPur >> &slTruthTable) const
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
std::vector< NeutrinoWithIndex > allMCTruth() const
double DetLength() const
void reconfigure(const fhicl::ParameterSet &pset)
DEFINE_ART_MODULE(TestTMapFile)
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...
Elastic Arms vertex finding algorithm.
recovalid::CAFCutter fCAFCutter
the CAFCutter helper class
bool passCuts(int cut, const caf::StandardRecord *sr)
Definition: CAFCutter.cxx:37
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
T get(std::string const &key) const
Definition: ParameterSet.h:231
int evt
Collect Geo headers and supply basic geometry functions.
std::string fCAFLabel
label for the process that made the standard records
double DetHalfHeight() const
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
TVector3 Unit() const
Vertex location in position and time.
const TLorentzVector & Position(const size_type) const
The accessor methods described above.
Definition: run.py:1
size_type size() const
Definition: PtrVector.h:308
std::string fCellHitLabel
label for the process that made the cell hits
OStream cout
Definition: OStream.cxx:6
int fCutType
what cuts to apply (see CAFCutter.h)
double Vx(const int i=0) const
Definition: MCParticle.h:220
std::string fVertexLabel
label for the process that made the vertices
T * make(ARGS...args) const
ElasticArmsValidate(fhicl::ParameterSet const &pset)
double DetHalfWidth() const
size_type size() const
Definition: MCTrajectory.h:165
TDirectory * dir
Definition: macro.C:5
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:219
void geom(int which=0)
Definition: geom.C:163
exit(0)
double Vz(const int i=0) const
Definition: MCParticle.h:222
Helper class for Reco Validation modules.
Definition: CAFCutter.h:58
void ProjectToBoxEdge(const double xyz[], const double dxyz[], double xlo, double xhi, double ylo, double yhi, double zlo, double zhi, double xyzout[])
Project along a direction from a particular starting point to the edge of a box.
Definition: Geo.cxx:38
bool NeutrinoSet() const
Definition: MCTruth.h:77
int Mode() const
Definition: MCNeutrino.h:149
int fCompMCNeutrino
don&#39;t know what this one does... :(
double Vy(const int i=0) const
Definition: MCParticle.h:221
Encapsulate the geometry of one entire detector (near, far, ndos)
bool failedToGet() const
Definition: Handle.h:196
TVector3 Vect() const
std::vector< const sim::Particle * > MCTruthToParticles(art::Ptr< simb::MCTruth > const &mct) const
enum BeamMode string