FuzzyKValidate_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file FuzzyKValidate.cxx
3 // \brief Analysis module for performance of FuzzyK
4 // \version $Id: FuzzyKValidate.cxx,v 1.7 2012-12-12 16:24:35 edniner Exp $
5 // \author $Author: edniner $
6 ////////////////////////////////////////////////////////////////////////
7 
8 // C/C++ includes
9 #include <cmath>
10 #include <iostream>
11 #include <map>
12 #include <string>
13 #include <vector>
14 
15 // ROOT includes
16 #include "TH1F.h"
17 #include "TH2F.h"
18 #include "TTree.h"
19 
20 // Framework includes
31 #include "fhiclcpp/ParameterSet.h"
33 
34 // NOvASoft includes
36 #include "Simulation/Particle.h"
40 
41 #include "Calibrator/Calibrator.h"
42 #include "MCCheater/BackTracker.h"
43 #include "RecoBase/CellHit.h"
44 #include "RecoBase/RecoHit.h"
45 #include "RecoBase/Prong.h"
46 #include "RecoBase/Vertex.h"
47 #include "Utilities/func/MathUtil.h"
48 #include "Utilities/AssociationUtil.h"
49 
50 // Required for CAFAna Cuts
52 #include "CAFAna/Cuts/Cuts.h"
53 #include "CAFAna/Cuts/TruthCuts.h"
54 #include "CAFAna/Cuts/TimingCuts.h"
55 #include "CAFAna/Vars/Vars.h"
56 
59 
60 
61 
62 namespace fuzz {
64  public:
65  explicit FuzzyKValidate(fhicl::ParameterSet const& pset);
66  ~FuzzyKValidate();
67 
68  void analyze(const art::Event& evt);
69  void reconfigure(const fhicl::ParameterSet& pset);
70  void beginJob();
71  void ClearVectors();
73  std::set<int> FindVisibleProngs3D(const art::PtrVector<rb::CellHit>& allhits, std::set<int> ids);
74 
75  private:
76 
77  /// MC information
78  std::string fSliceLabel; ///< label for the process that made the slices
79  std::string fCellHitLabel; ///< label for the process that made cell hits
80  std::string fTrackLabel; ///< label for the process that made tracks
81  std::string fVertexLabel; ///< label for the process that made the vertices
82  std::string fAssn3DLabel; ///< association label for 3D prongs
83  std::string fAssn2DLabel; ///< association label for 2D prongs
84  int fCompMCNeutrino; ///< ???
85 
86  std::string fCAFLabel; ///< label for the process that made standard records
87  bool fApplyCAFCuts; ///< should we apply the CAF-level cuts?
88  int fCutType; ///< what cuts to apply (see CAFCutter.h)
89 
90  recovalid::CAFCutter fCAFCutter; ///< the CAFCutter helper class
91 
92  TTree* fOutTree;
93  //event summary info
94  int run;
95  int subrun;
96  int event;
97 
98  //reco Prong info and direction
99  int nt3D; //number 3D Prongs
100  int nt2D; //number 2D Prongs
101  std::vector<int> view; //Prong view
102  std::vector<float> sX; //start position x
103  std::vector<float> sY;//start position y
104  std::vector<float> sZ;//start position z
105  std::vector<float> dX; //direction x
106  std::vector<float> dY;//direction y
107  std::vector<float> dZ;//direction z
108  std::vector<int> ncells; //num cells total
109  std::vector<int> ncellsX;//x cells
110  std::vector<int> ncellsY;//y cells
111  std::vector<float> gevTot; //total GeV
112  std::vector<float> gevX;//x GeV
113  std::vector<float> gevY;//y GeV
114  std::vector<float> totGeV;
115  std::vector<float> len; //Prong length
116  //Prong comparisons these are only for 3D Prongs
117  //below vectors are for true visible primaries that were matched to a reco Prong
118  std::vector<int> matchedPDG;//pdg of visible particle being matched
119  std::vector<int> matchedProngIndx;//index of reco Prong being matched
120  std::vector<float> matchedVET;//visible energy of true paricle in entire slice
121  std::vector<float> matchedVER;//visible energy of true particle captured in Prong
122  std::vector<float> matchedET;//total energy of true particle
123  std::vector<float> matchedXStartT;//x start of true particle
124  std::vector<float> matchedYStartT;//y start of true particle
125  std::vector<float> matchedZStartT;//z start of true particle
126  std::vector<float> matchedPXT;//x momentum of true particle
127  std::vector<float> matchedPYT;//y momentum of true particle
128  std::vector<float> matchedPZT;//z momentum of true particle
129  std::vector<float> matchedAngleDiff;//angle difference between reco Prong and truth
130  std::vector<float> matchedEff;//efficiency of match
131  std::vector<float> matchedPur;//purity of match
132  std::vector<int> matchedMother;//pdg of mother particle
133  //the below are for true visible primaries that were NOT matched to a reco Prong
134  std::vector<int> unmatchedMother;//mother pdg code
135  std::vector<int> unmatchedPDG;//unmatched primary pdg code
136  std::vector<float> unmatchedVET;//visible energy of unmatched in the slice
137  std::vector<float> unmatchedET;//total energy of unmatched
138  std::vector<float> unmatchedXStartT;//x start of true particle
139  std::vector<float> unmatchedYStartT;//y start of true particle
140  std::vector<float> unmatchedZStartT;//z start of true particle
141  std::vector<float> unmatchedPXT;//x momentum of true particle
142  std::vector<float> unmatchedPYT;//y momentum of true particle
143  std::vector<float> unmatchedPZT;//z momentum of true particle
144  //below vectors are for reco Prongs not matched as first choice to a visible primary. The matches made are too the best option
145  std::vector<int> danglerPDG;//pdg of visible particle being matched
146  std::vector<int> danglerProngIndx;//index of reco Prong being matched
147  std::vector<float> danglerVET;//visible energy of true paricle in entire slice
148  std::vector<float> danglerVER;//visible energy of true particle captured in Prong
149  std::vector<float> danglerET;//total energy of true particle
150  std::vector<float> danglerXStartT;//x start of true particle
151  std::vector<float> danglerYStartT;//y start of true particle
152  std::vector<float> danglerZStartT;//z start of true particle
153  std::vector<float> danglerPXT;//x momentum of true particle
154  std::vector<float> danglerPYT;//y momentum of true particle
155  std::vector<float> danglerPZT;//z momentum of true particle
156  std::vector<float> danglerAngleDiff;//angle difference between reco Prong and truth
157  std::vector<float> danglerEff;//efficiency of match
158  std::vector<float> danglerPur;//purity of match
159  std::vector<int> danglerMother;//pdg of mother particle
160  int matchedProngs;//number of successfully matched primaries
161  int unmatchedProngs;//number of unmatched primaries
162  int danglerProngs;//number of reco Prongs not matched
163 
164 
165  //overall optimization score for the event
166  double eventScore;
167 
168 
169  //Prong histograms
170  TH1F* fHistoNt;
171  TH1F* fHisto3D;
172  TH1F* fHisto2D;
174  TH1F* fHistoLen;
175  //TH1F* fHistoNCell;
176  TH1F* fHistoSX;
177  TH1F* fHistoSY;
178  TH1F* fHistoSZ;
179  TH1F* fHistoDX;
180  TH1F* fHistoDY;
181  TH1F* fHistoDZ;
182  TH1F* fHistoGeV;
188  //Prong comparisons
189  TH1F* fHistoEff;
190  TH1F* fHistoPur;
191  TH1F* fHistoAngle;
194  TH1F* fHistoLeft;
198  //TH2F* fHistoUvL;
199  //TH1F* fHistoUnVisE;
200  //TH1F* fHistoMVisE;
201  //TH2F* fHistoEffvEvis;
202  //TH2F* fHistoPurvEvis;
203  //TH2F* fHistoAnglevEvis;
204 
205  //efficiency vs purity for all visible primaries matched to their best reco Prong
206  //TH2F* fHistoEffvPur;
207  //TH2F* fHistoDEffvPur;
208  //TH1F* fHistoMatchMetric; //store value of match quality between reco Prong and visible particle
209 
210  //store event optimization score
211  //TH1F* fHistoEventScore;
212 
213  //neutrino information
214  float Vx;//start x
215  float Vy;//start y
216  float Vz;//start z
217  float NuEff;//efficiency of neutrino for slice
218  float NuPur;//purity of neutrino for slice
219  float NuEVSl;//visible energy of neutrino in slice
220  float NuEVTot;//visible energy of neutrino in all slices
221  float x;
222  float y;
223  float q2;
224  float W2;
225  int NuPdg;//neutrino pdg
226  int ccnc;//is cc or nc
227  int mode;//interaction mode
228  float NuE;//total energy of neutrino
229 
230  float rVx;//start x
231  float rVy;//start y
232  float rVz;//start z
233 
234 
235 
236 
237 
238  };
239 }
240 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
241 
242 namespace fuzz
243 {
244  //.......................................................................
246  : EDAnalyzer(pset)
247  {
248  this->reconfigure(pset);
249  }
250 
251  //......................................................................
253  {
254  //======================================================================
255  // Clean up any memory allocated by your module
256  //======================================================================
257  }
258 
259  //......................................................................
261  {
262 
263  fSliceLabel = pset.get< std::string >("SliceLabel");
264  fCellHitLabel = pset.get< std::string >("CellHitLabel");
265  fVertexLabel = pset.get< std::string >("VertexLabel");
266  fTrackLabel = pset.get< std::string >("TrackLabel");
267  fAssn3DLabel = pset.get< std::string >("Assn3DLabel");
268  fAssn2DLabel = pset.get< std::string >("Assn2DLabel");
269  fCompMCNeutrino = pset.get< int >("CompMCNeutrino");
270  fCAFLabel = pset.get< std::string >("CAFLabel");
271  fApplyCAFCuts = pset.get< bool >("ApplyCAFCuts");
272  fCutType = pset.get< int >("CutType");
273 
274  }
275 
276  //......................................................................
278  {
279 
281 
282  // define the tree
283  fOutTree = tfs->make<TTree>("FuzzyKValidate","FuzzyK Validation");
284 
285  //store event summary information
286  fOutTree->Branch("run", &run, "run/I");
287  fOutTree->Branch("subrun", &subrun, "subrun/I");
288  fOutTree->Branch("event", &event, "event/I");
289 
290  //store Prong info
291  fOutTree->Branch("nt3D", &nt3D);
292  fOutTree->Branch("nt2D", &nt2D);
293  fOutTree->Branch("view", &view);
294  fOutTree->Branch("sX", &sX);
295  fOutTree->Branch("sY", &sY);
296  fOutTree->Branch("sZ", &sZ);
297  fOutTree->Branch("dX", &dX);
298  fOutTree->Branch("dY", &dY);
299  fOutTree->Branch("dZ", &dZ);
300  fOutTree->Branch("ncells", &ncells);
301  fOutTree->Branch("ncellsX", &ncellsX);
302  fOutTree->Branch("ncellsY", &ncellsY);
303  fOutTree->Branch("gevTot", &gevTot);
304  fOutTree->Branch("gevX", &gevX);
305  fOutTree->Branch("gevY", &gevY);
306  fOutTree->Branch("totGeV", &totGeV);
307  fOutTree->Branch("len", &len);
308  //store Prong comparisons
309  fOutTree->Branch("matchedPDG", &matchedPDG);
310  fOutTree->Branch("matchedProngIndx", &matchedProngIndx);
311  fOutTree->Branch("matchedVET", &matchedVET);
312  fOutTree->Branch("matchedVER", &matchedVER);
313  fOutTree->Branch("matchedET", &matchedET);
314  fOutTree->Branch("matchedXStartT", &matchedXStartT);
315  fOutTree->Branch("matchedYStartT", &matchedYStartT);
316  fOutTree->Branch("matchedZStartT", &matchedZStartT);
317  fOutTree->Branch("matchedPXT", &matchedPXT);
318  fOutTree->Branch("matchedPYT", &matchedPYT);
319  fOutTree->Branch("matchedPZT", &matchedPZT);
320  fOutTree->Branch("matchedAngleDiff", &matchedAngleDiff);
321  fOutTree->Branch("matchedEff", &matchedEff);
322  fOutTree->Branch("matchedPur", &matchedPur);
323  fOutTree->Branch("matchedMother", &matchedMother);
324  fOutTree->Branch("unmatchedMother", &unmatchedMother);
325  fOutTree->Branch("unmatchedPDG", &unmatchedPDG);
326  fOutTree->Branch("unmatchedVET", &unmatchedVET);
327  fOutTree->Branch("unmatchedXStartT", &unmatchedXStartT);
328  fOutTree->Branch("unmatchedYStartT", &unmatchedYStartT);
329  fOutTree->Branch("unmatchedZStartT", &unmatchedZStartT);
330  fOutTree->Branch("unmatchedPXT", &unmatchedPXT);
331  fOutTree->Branch("unmatchedPYT", &unmatchedPYT);
332  fOutTree->Branch("unmatchedPZT", &unmatchedPZT);
333  fOutTree->Branch("matchedProngs", &matchedProngs);
334  fOutTree->Branch("unmatchedProngs", &unmatchedProngs);
335  fOutTree->Branch("danglerProngs", &danglerProngs);
336  fOutTree->Branch("danglerPDG", &danglerPDG);
337  fOutTree->Branch("danglerProngIndx", &danglerProngIndx);
338  fOutTree->Branch("danglerVET", &danglerVET);
339  fOutTree->Branch("danglerVER", &danglerVER);
340  fOutTree->Branch("danglerET", &danglerET);
341  fOutTree->Branch("danglerXStartT", &danglerXStartT);
342  fOutTree->Branch("danglerYStartT", &danglerYStartT);
343  fOutTree->Branch("danglerZStartT", &danglerZStartT);
344  fOutTree->Branch("danglerPXT", &danglerPXT);
345  fOutTree->Branch("danglerPYT", &danglerPYT);
346  fOutTree->Branch("danglerPZT", &danglerPZT);
347  fOutTree->Branch("danglerAngleDiff", &danglerAngleDiff);
348  fOutTree->Branch("danglerEff", &danglerEff);
349  fOutTree->Branch("danglerPur", &danglerPur);
350  fOutTree->Branch("danglerMother", &danglerMother);
351 
352  //store overall event score for optimization
353  //score = (1 - E_dangle/E_total)sqrt(sum[ (E_i/E_tot)^2*(eff^2 + pur^2)*0.5] )
354  //the sum is over all the visible particles
355  //E_i is the energy of the true particle used in the match and E_tot is the total visible energy in the slice
356  //a visible primary or a reco Prong with no match goes into the score as a 0,0
357  //The factor (1- E_dangle/E_total) is a weight factor that accounts for unmatched reco Prongs.
358  //E_dangle/E_total is the fraction of fls energy put into reco Prongs that were not associated with true particles
359  //this score will have a value between 0 and 1, the higher the score the better the event performance
360  fOutTree->Branch("eventScore", &eventScore);
361 
362 
363  //Prong histograms
364  fHistoNt = tfs->make<TH1F>("TotalProngs","Number of Prongs (3D+2D) per vertex;num 2D+3D Prongs",17,-0.5,16.5);
365  fHisto3D = tfs->make<TH1F>("Num3DProngs","Number of 3D Prongs per vertex;num 3D Prongs",9,-0.5,8.5);
366  fHisto2D = tfs->make<TH1F>("Num2DProngs","Number of 2D Prongs per vertex;num 2D Prongs",9,-0.5,8.5);
367  fHisto3Dvs2D = tfs->make<TH2F>("NumProngs3Dvs2D","3D Prongs vs 2D Prongs;num 2D Prongs;num 3D Prongs",9,-0.5,8.5,9,-0.5,8.5);
368  fHistoLen = tfs->make<TH1F>("ProngLength","Prong length;prong length (cm)",350,0.0,7000.0);
369  fHistoSX = tfs->make<TH1F>("ProngStartX","Start X coordinate ;prong start x (cm)",160,-800.0,800.0);
370  fHistoSY = tfs->make<TH1F>("ProngStartY","Start Y coordinate ;prong start y (cm)",160,-800.0,800.0);
371  fHistoSZ = tfs->make<TH1F>("ProngStartZ","Start Z coordinate ;prong start z (cm)",120,0.0,6000.0);
372  fHistoDX = tfs->make<TH1F>("ProngDirX","X direction ;prong x dir",101,-1.01,1.01);
373  fHistoDY = tfs->make<TH1F>("ProngDirY","Y direction ;prong y dir",101,-1.01,1.01);
374  fHistoDZ = tfs->make<TH1F>("ProngDirZ","Z direction ;prong z dir",101,-1.01,1.01);
375  fHistoNCells = tfs->make<TH1F>("ProngCells","Number of cells in Prong; cells",200,0.0,400.0);
376  fHistoNCellsX = tfs->make<TH1F>("ProngCellsX","Number of cells in Prong X view;prong x cells",150,0.0,300.0);
377  fHistoNCellsY = tfs->make<TH1F>("ProngCellsY","Number of cells in Prong Y view;prong y cells",150,0.0,300.0);
378  fHistoNCellsXY = tfs->make<TH2F>("ProngCellsYvsX","Y cells vs. X cells;prong x cells;prong y cells"
379  ,150,0.0,300.0,150,0.0,300.0);
380  fHistoGeV = tfs->make<TH1F>("ProngGeV","GeV of Prong ;prong energy (GeV)",500,0.0,50.0);
381  fHistoGeVYvX = tfs->make<TH2F>("ProngGeVYvX","GeVY vs GeVX ;prong energy X (GeV);prong energy Y (GeV)",
382  250,0.0,25.0,250,0.0,25.0);
383  fHistoEff = tfs->make<TH1F>("ProngEfficiency","Efficiency of matched visible particles;prong efficiency",101,0.0,1.01);
384  fHistoPur = tfs->make<TH1F>("ProngPurity","Purity of matched visible particles;prong purity",101,0.0,1.01);
385  fHistoAngle = tfs->make<TH1F>("ProngMatchAngle","Angle between true particle and reco Prong;angle between true and matched reco prong (theta)",91,0.0,182.0);
386  fHistoMatched = tfs->make<TH1F>("NumMatchedProngs","Number of visible primaries matched to reco Prongs;number of reco prongs matched to true particle"
387  ,13,-0.5,12.5);
388  fHistoUnMatched = tfs->make<TH1F>("UnMatchedPrimaries","Number of visible primaries unmatched to reco Prongs;num visible primaries unmatched to reco prongs"
389  ,13,-0.5,12.5);
390  fHistoLeft = tfs->make<TH1F>("UnMatchedRecoProngs","Number of reco Prongs not matched;num reco prongs unmatched to true particles",9,-0.5,8.5);
391  fHistoLeptonEff = tfs->make<TH1F>("ProngPrimaryLeptonEfficiency","Completeness of matched primary leptons;prong efficiency for primary lepton",101,0.0,1.01);
392  fHistoLeptonPur = tfs->make<TH1F>("ProngPrimaryLeptonPurity","Purity of matched primary leptons;prong purity for primary lepton",101,0.0,1.01);
393  fHistoLeptonAngle = tfs->make<TH1F>("PrimaryLeptonAngle","Angle between primary lepton and reco Prong;angle between primary lepton and reco prong (theta)",91,0.0,182.0);
394  //fHistoUvL = tfs->make<TH2F>("fHistoUvL","Unmatched primaries vs Unmatched reco Prongs;reco Prongs;primaries",9,-0.5,8.5,9,-0.5,8.5);
395  //fHistoUnVisE = tfs->make<TH1F>("fHistoUnVisE","Visible energy of unmatched primaries;GeV",250,0.0,5.0);
396  //fHistoMVisE = tfs->make<TH1F>("fHistoMVisE","Visible energy of matched primaries;GeV",250,0.0,5.0);
397  //fHistoEffvEvis = tfs->make<TH2F>("fHistoEffvEvis","Completeness vs Visible energy for match;visE(GeV);eff",250,0.0,5.0,101,0.0,1.01);
398  //fHistoPurvEvis = tfs->make<TH2F>("fHistoPurvEvis","Purity vs Visible energy for match;visE(GeV);purity",250,0.0,5.0,101,0.0,1.01);
399  //fHistoAnglevEvis = tfs->make<TH2F>("fHistoAnglevEvis","Angle diff vs Visible energy for match;visE(GeV);theta",250,0.0,5.0,91,0.0,182.0);
400  //fHistoEffvPur = tfs->make<TH2F>("fHistoEffvPur","Completeness vs Purity for visible primaries;purity;eff",101,0.0,1.01,101,0.0,1.01);
401 
402  //fHistoDEffvPur = tfs->make<TH2F>("fHistoDEffvPur","Completeness vs Purity for dangling reco Prongs;purity;eff",101,0.0,1.01,101,0.0,1.01);
403 
404  //fHistoMatchMetric = tfs->make<TH1F>("fHistoMatchMetric","|E_truepart - E_recoProng|/E_truepart",200,0.0,100.0);
405 
406  //fHistoEventScore = tfs->make<TH1F>("fHistoEventScore","Event Score for optimization;score;events",101,0.0,1.01);
407 
408  //store neutrino info
409  fOutTree->Branch("Vx", &Vx);//neutrino position info
410  fOutTree->Branch("Vy", &Vy);
411  fOutTree->Branch("Vz", &Vz);
412  fOutTree->Branch("NuEff", &NuEff);//info on quality of neutrino in slice
413  fOutTree->Branch("NuPur", &NuPur);
414  fOutTree->Branch("NuEVSl", &NuEVSl);
415  fOutTree->Branch("NuEVTot", &NuEVTot);
416  fOutTree->Branch("NuPdg", &NuPdg);
417  fOutTree->Branch("ccnc", &ccnc);
418  fOutTree->Branch("mode", &mode);
419  fOutTree->Branch("NuE", &NuE);
420  fOutTree->Branch("y", &y);
421  fOutTree->Branch("x", &x);
422  fOutTree->Branch("q2", &q2);
423  fOutTree->Branch("W2", &W2);
424 
425  fOutTree->Branch("rVx", &rVx);
426  fOutTree->Branch("rVy", &rVy);
427  fOutTree->Branch("rVz", &rVz);
428 
429 
430  }
431 
432  //......................................................................
434  {
436 
438  evt.getByLabel(fSliceLabel, sHandle);
439  if(sHandle.failedToGet()){
440  std::cout<<"No slices found, skipping event" << std::endl;
441  }
443  for(unsigned int i=0; i<sHandle->size(); ++i){
444  slices.push_back(art::Ptr<rb::Cluster>(sHandle, i));
445  }
446 
447  //get all the hits to calculate neutrino efficiency/purity for slices
449  evt.getByLabel(fCellHitLabel, htHandle);
451  for(unsigned int i=0; i<htHandle->size(); ++i){
452  art::Ptr<rb::CellHit> hit(htHandle,i);
453  allHits.push_back(hit);
454  }
455 
456  //get associations between slices and vertices
457  art::FindManyP<rb::Vertex> fmVert(sHandle, evt, fVertexLabel);
458  //try to get associations between Prongs and slices
461 
462  // Get the FMP for the Standard Records
463  art::FindManyP<caf::StandardRecord> recordFMP(sHandle, evt, fCAFLabel);
464 
465  //if these are events with neutrino information to compare too, get it
466  std::vector<int> bestNuId;
467  std::vector<cheat::NeutrinoWithIndex> nus;
468  std::vector<std::vector<cheat::NeutrinoEffPur>> sEffPur;
469  if(fCompMCNeutrino == 1){
470  nus = bt->allMCTruth();
471  //std::cout<<allNus.size()<<std::endl;
472  sEffPur = bt->SlicesToMCTruthsTable(slices);
473  //std::cout<<"Done matching"<<std::endl;
474 
475  bestNuId = bt->SliceToOrderedNuIdsByEff(nus, sEffPur);
476  }
477 
478 
479  //loop over slices
480  for( unsigned int i=0; i<slices.size(); ++i){
481  this->ClearVectors();
482 
483 
484  // Apply the CAF-level cuts
485  bool pass = true;
486  if( fApplyCAFCuts ) {
487  // get record associated with the slice
488  std::vector< art::Ptr<caf::StandardRecord> > records = recordFMP.at(i);
489  if (records.size() > 0 && recordFMP.isValid()) {
490  pass = fCAFCutter.passCuts(fCutType, records[0].get());
491  }
492  else { pass = false; }
493  }
494 
495  if(!pass) continue;
496 
497 
498 
499  //skip noise slice
500  if (slices[i]->IsNoise()) continue;
501  //if we are comparing to truth, make sure there is a neutrino with this slice
502  std::vector<art::Ptr<rb::Prong>> fk3D;
503  std::vector<art::Ptr<rb::Prong>> fk2D;
504  //if there are associations between Prongs and slices use those
505  if(fmSlProng3D.isValid()){
506  fk3D = fmSlProng3D.at(i);
507  fk2D = fmSlProng2D.at(i);
508  }
509  //else see if there are Prongs connected to vertex
510  else if(fmVert.isValid()){
511  std::vector<art::Ptr<rb::Vertex>> v = fmVert.at(i);
512  //should only be one vertex
513  if (v.size() != 1) continue;
516  //skip if there are no prongs
517  if(!(fmVertProng3D.isValid())) continue;
518  //only one vertex per slice so only need index 0
519  fk3D = fmVertProng3D.at(0);
520  fk2D = fmVertProng2D.at(0);
521 
522  rVx = v[0]->GetX();
523  rVy = v[0]->GetY();
524  rVz = v[0]->GetZ();
525  }
526  else{
527  continue;
528  }
529 
530  run = evt.run();
531  subrun = evt.subRun();
532  event = evt.id().event();
533 
534  //fill Prong statistics
535  nt3D = fk3D.size();
536  nt2D = fk2D.size();
537  for(unsigned int it=0; it<fk3D.size(); ++it){
538  view.push_back(fk3D[it]->View());
539  sX.push_back(fk3D[it]->Start().X());
540  sY.push_back(fk3D[it]->Start().Y());
541  sZ.push_back(fk3D[it]->Start().Z());
542  dX.push_back(fk3D[it]->Dir().X());
543  dY.push_back(fk3D[it]->Dir().Y());
544  dZ.push_back(fk3D[it]->Dir().Z());
545  ncells.push_back(fk3D[it]->NCell());
546  ncellsX.push_back(fk3D[it]->NXCell());
547  ncellsY.push_back(fk3D[it]->NYCell());
548  gevTot.push_back(fk3D[it]->TotalGeV());
549  totGeV.push_back(fk3D[it]->TotalGeV());
550  float gevCX = GeVCalc(fk3D[it], geo::kX);
551  float gevCY = GeVCalc(fk3D[it], geo::kY);
552  gevX.push_back(gevCX);
553  gevY.push_back(gevCY);
554  len.push_back(fk3D[it]->TotalLength());
555 
556  fHistoLen->Fill(len.back());
557  fHistoSX->Fill(sX.back());
558  fHistoSY->Fill(sY.back());
559  fHistoSZ->Fill(sZ.back());
560  fHistoDX->Fill(dX.back());
561  fHistoDY->Fill(dY.back());
562  fHistoDZ->Fill(dZ.back());
563  fHistoGeV->Fill(gevTot.back());
564  fHistoGeVYvX->Fill(gevX.back(),gevY.back());
565  fHistoNCells->Fill(ncells.back());
566  fHistoNCellsX->Fill(ncellsX.back());
567  fHistoNCellsY->Fill(ncellsY.back());
568  fHistoNCellsXY->Fill(ncellsX.back(),ncellsY.back());
569  }
570  fHistoNt->Fill(nt3D+nt2D);
571  fHisto3D->Fill(nt3D);
572  fHisto2D->Fill(nt2D);
573  fHisto3Dvs2D->Fill(nt2D,nt3D);
574  //now fill 2D Prong info
575  for(unsigned int it=0; it<fk2D.size(); ++it){
576  view.push_back(fk2D[it]->View());
577  sX.push_back(fk2D[it]->Start().X());
578  sY.push_back(fk2D[it]->Start().Y());
579  sZ.push_back(fk2D[it]->Start().Z());
580  dX.push_back(fk2D[it]->Dir().X());
581  dY.push_back(fk2D[it]->Dir().Y());
582  dZ.push_back(fk2D[it]->Dir().Z());
583  ncells.push_back(fk2D[it]->NCell());
584  ncellsX.push_back(fk2D[it]->NXCell());
585  ncellsY.push_back(fk2D[it]->NYCell());
586  gevTot.push_back(fk2D[it]->TotalGeV());
587  totGeV.push_back(fk2D[it]->TotalGeV());
588  float gevCX = GeVCalc(fk2D[it], geo::kX);
589  float gevCY = GeVCalc(fk2D[it], geo::kY);
590  gevX.push_back(gevCX);
591  gevY.push_back(gevCY);
592  len.push_back(fk2D[it]->TotalLength());
593  fHistoLen->Fill(len.back());
594  if (fk2D[it]->View() == geo::kX){
595  fHistoSX->Fill(sX.back());
596  fHistoDX->Fill(dX.back());
597  fHistoNCellsX->Fill(ncellsX.back());
598  }
599  else if(fk2D[it]->View() == geo::kY){
600  fHistoSY->Fill(sY.back());
601  fHistoDY->Fill(dY.back());
602  fHistoNCellsY->Fill(ncellsY.back());
603  }
604  fHistoNCells->Fill(ncells.back());
605  fHistoSZ->Fill(sZ.back());
606  fHistoDZ->Fill(dZ.back());
607  fHistoGeV->Fill(gevTot.back());
608  }
609 
610  //if we do not want to compare to truth, dump into to tree and move to next slice
611  if(fCompMCNeutrino != 1){
612  fOutTree->Fill();
613  continue;
614  }
615  else if (bestNuId[i] == -1){
616  fOutTree->Fill();
617  continue;
618  }
619 
620  //now access true neutrino info for further validation
621  art::Ptr<simb::MCTruth> nuTru = sEffPur[i][bestNuId[i]].neutrinoInt;
622 
623  NuEff = sEffPur[i][bestNuId[i]].efficiency;
624  NuPur = sEffPur[i][bestNuId[i]].purity;
625  NuEVSl = sEffPur[i][bestNuId[i]].energySlice;
626  NuEVTot = sEffPur[i][bestNuId[i]].energyTotal;
627 
628  //if this a neutrino object, store this information, if this is a cosmic,
629  //then skip
630  if (nuTru->NeutrinoSet()){
631  NuPdg = nuTru->GetNeutrino().Nu().PdgCode();
632  ccnc = nuTru->GetNeutrino().CCNC();
633  mode = nuTru->GetNeutrino().Mode();
634  NuE = nuTru->GetNeutrino().Nu().Momentum().E();
635 
636  q2 = nuTru->GetNeutrino().QSqr();
637  x = nuTru->GetNeutrino().X();
638  y = nuTru->GetNeutrino().Y();
639  W2 = nuTru->GetNeutrino().W()*nuTru->GetNeutrino().W();
640 
641  //now get vertex info and make comparisons
642  Vx = nuTru->GetNeutrino().Nu().Vx();
643  Vy = nuTru->GetNeutrino().Nu().Vy();
644  Vz = nuTru->GetNeutrino().Nu().Vz();
645  }
646 
647 
648  //now go about finding list of primary visible daughters sorted from largest to smallest visible energy
649  std::vector<cheat::TrackIDE> sliceProngs = bt->HitsToTrackIDE(slices[i]->AllCells());
650  std::set<int> sT;
651  for(unsigned int j=0; j<sliceProngs.size(); ++j){
652  sT.insert(sliceProngs[j].trackID);
653  }
654  std::set<int> visIds = FindVisibleProngs3D(slices[i]->AllCells(), sT);
655 
656  //std::cout<<"There are: "<<visIds.size()<<" visible particles."<<std::endl;
657 
658  std::vector<const sim::Particle*> particles = bt->MCTruthToParticles(nuTru);
659  //std::cout<<"Checking against the list of: "<<particles.size()<<" particles."<<std::endl;
660 
661  //now keep particles that meet the visible criteria and are primary.
662  //NOTE: testing using all visible partilces and not just the primaries. The variable names will
663  //be changed to reflect this if this change works. For now visPrimaries really means visParticles
664  std::vector<const sim::Particle*> visPrimaries;
665  for(std::set<int>::iterator itr = visIds.begin(); itr!=visIds.end(); itr++){
666  //std::cout<<"Looking for primary that matches visible Prong id: "<<(*itr)<<std::endl;
667  for(unsigned int k=0; k<particles.size(); ++k){
668  //std::cout<<"Checking particle: "<<k<<" with id: "<<particles[k]->ProngId()<<" and pdg: "<<particles[k]->PdgCode()<<std::endl;
669  if (visIds.find(particles[k]->TrackId()) == visIds.end()) continue;
670  if (*(visIds.find(particles[k]->TrackId())) == (*itr)){
671  visPrimaries.push_back(particles[k]);
672  break;
673  }
674  /*if ((*(visIds.find(particles[k]->TrackId())) == (*itr)) && (particles[k]->Mother() == 0)){
675  visPrimaries.push_back(particles[k]);
676  std::cout<<"Found a visible primary: "<<particles[k]->PdgCode()<<" with track id "<<particles[k]->TrackId()<<std::endl;
677  break;
678  }
679  //in the case that these are pi0 daughters make sure they get counted
680  if ((particles[k]->Mother() != 0) && (*(visIds.find(particles[k]->TrackId())) == (*itr))){
681  const sim::Particle* mTest = bt->TrackIDToMotherParticle(particles[k]->TrackId());
682  if ((mTest->PdgCode() == 111) && (mTest->Mother() == 0)){
683  visPrimaries.push_back(particles[k]);
684  std::cout<<"Found a visible primary: "<<particles[k]->PdgCode()<<std::endl;
685  break;
686  }
687  }*/
688  }
689  }
690 
691  //initialize score for slice
692  double score = 0.0;
693 
694  //now loop through visible primaries, starting with highest energy.
695  //for each one find the 3D track for which it is most efficient.
696  //visible primaries are allowed to match to the same reco track if that is what is desired
697  int matched = 0;
698  int unmatched = 0;
699  std::set<int> usedProngs;
700  std::map<int, double> purMap, effMap;
701  std::map<int, int> parents;
702  //std::cout<<"There are: "<<visPrimaries.size()<<" primary particles. Making matches."<<std::endl;
703  for(unsigned int j=0; j<visPrimaries.size(); ++j){
704  //std::cout<<"Matching particle: "<<j<<std::endl;
705  std::set<int> temp;
706  temp.insert(visPrimaries[j]->TrackId());
707  double eff, pur;//, matchMetric;
708  double maxEff = -1.0;
709  double maxPur = -1.0;
710  int maxTr = -1;
711  double dE, tE, dEMax, tEMax;
712  dEMax = -1.0;
713  tEMax = -1.0;
714  //match a truth particle to reco track with metric |E_truepart - E_recotrack|/E_truepart
715  //a good match should have the smallest metric possible
716  // double maxMatchMetric = 9999.0;
717  for(unsigned int k=0; k<fk3D.size(); ++k){
718  //if (usedTracks.find(k) != usedTracks.end()) continue;
719  eff = bt->HitCollectionEfficiency(temp, fk3D[k]->AllCells(), slices[i]->AllCells(), geo::kXorY, &effMap,true,&dE,&tE);
720  pur = bt->HitCollectionPurity(temp, fk3D[k]->AllCells(), &purMap, &parents,true);
721  // if (pur > 0.0) matchMetric = fabs(tE - dE/pur)/tE;
722  // else matchMetric = 9999.0;
723  //std::cout<<"***************With reco track: "<<k<<" efficiency: "<<eff<<" and metric: "<<matchMetric<<std::endl;
724  if(eff > maxEff){
725  //if(matchMetric < maxMatchMetric){
726  // maxMatchMetric = matchMetric;
727  maxEff = eff;
728  maxPur = pur;
729  maxTr = k;
730  dEMax = dE;
731  tEMax = tE;
732  }
733  //if efficieny is the same, tiebreak on purity
734  else if ((eff == maxEff) && (pur > maxPur)){
735  // maxMatchMetric = matchMetric;
736  maxEff = eff;
737  maxPur = pur;
738  maxTr = k;
739  dEMax = dE;
740  tEMax = tE;
741  }
742  }
743  if (maxEff > 0.0){
744  //if (maxMatchMetric < 9999.0){
745  usedProngs.insert(maxTr);
746  const sim::Particle* mTest = bt->TrackIDToMotherParticle(visPrimaries[j]->TrackId());
747  matchedMother.push_back(mTest->PdgCode());
748  matched++;
749  matchedPDG.push_back(visPrimaries[j]->PdgCode());
750  matchedProngIndx.push_back(maxTr);
751  matchedVET.push_back(tEMax);
752  matchedVER.push_back(dEMax);
753  matchedET.push_back(visPrimaries[j]->E());
754  matchedXStartT.push_back(visPrimaries[j]->Vx());
755  matchedYStartT.push_back(visPrimaries[j]->Vy());
756  matchedZStartT.push_back(visPrimaries[j]->Vz());
757  matchedPXT.push_back(visPrimaries[j]->Px());
758  matchedPYT.push_back(visPrimaries[j]->Py());
759  matchedPZT.push_back(visPrimaries[j]->Pz());
760  float costheta = (dX[maxTr]*matchedPXT.back() + dY[maxTr]*matchedPYT.back() + dZ[maxTr]*matchedPZT.back())/
761  (util::pythag(dX[maxTr],dY[maxTr],dZ[maxTr])*util::pythag(matchedPXT.back(),matchedPYT.back()
762  ,matchedPZT.back()));
763  matchedAngleDiff.push_back(acos(costheta)*180.0/M_PI);
764  matchedEff.push_back(maxEff);
765  matchedPur.push_back(maxPur);
766 
767  fHistoEff->Fill(maxEff);
768  fHistoPur->Fill(maxPur);
769  fHistoAngle->Fill(matchedAngleDiff.back());
770  if ((abs(matchedPDG.back()) == 13)){
771  if (matchedMother.back() != 111){
772  fHistoLeptonEff->Fill(maxEff);
773  fHistoLeptonPur->Fill(maxPur);
774  fHistoLeptonAngle->Fill(matchedAngleDiff.back());
775  }
776  }
777  //fHistoMVisE->Fill(tEMax);
778  //fHistoEffvEvis->Fill(tEMax,maxEff);
779  //fHistoPurvEvis->Fill(tEMax,maxPur);
780  //fHistoAnglevEvis->Fill(tEMax,matchedAngleDiff.back());
781  //fHistoEffvPur->Fill(maxPur,maxEff);
782 
783  score += (tEMax*NuEff/NuEVSl)*(tEMax*NuEff/NuEVSl)*(maxEff*maxEff + maxPur*maxPur)*0.5;
784  //fHistoMatchMetric->Fill(maxMatchMetric);
785 
786  //std::cout<<"Made a match for particle: "<<j<<" to reco track: "<<maxTr<<" with eff: "<<maxEff<<" and pur: "<< maxPur<<" score is now: "<<score<<std::endl;
787  //std::cout<<"The true particle had total energy: "<<tEMax<<std::endl;
788  }
789  else{
790  //std::cout<<"No match made for particle: "<<j<<std::endl;
791  unmatched++;
792  const sim::Particle* mTest = bt->TrackIDToMotherParticle(visPrimaries[j]->TrackId());
793  unmatchedMother.push_back(mTest->PdgCode());
794  unmatchedPDG.push_back(visPrimaries[j]->PdgCode());
795  unmatchedVET.push_back(tEMax);
796  unmatchedET.push_back(visPrimaries[j]->E());
797  unmatchedXStartT.push_back(visPrimaries[j]->Vx());
798  unmatchedYStartT.push_back(visPrimaries[j]->Vy());
799  unmatchedZStartT.push_back(visPrimaries[j]->Vz());
800  unmatchedPXT.push_back(visPrimaries[j]->Px());
801  unmatchedPYT.push_back(visPrimaries[j]->Py());
802  unmatchedPZT.push_back(visPrimaries[j]->Pz());
803  //fHistoUnVisE->Fill(tEMax);
804  //fHistoEffvPur->Fill(0.0,0.0);
805  }
806  }
807 
808  //now if there are unmatched 3D tracks, look at them and get the info for the most efficient visible primary
809  int dangler = 0.0;
810  std::set<art::Ptr<rb::CellHit>> dCells; //store energy from cells from all danglers to calculate energy weight
811  for (unsigned int j=0; j<fk3D.size(); j++){
812  if (usedProngs.find(j) != usedProngs.end()) continue;
813  //std::cout<<"Reco track: "<<j<<" with "<<fk3D[j]->NCell()<<" cells has no match, looking now."<<std::endl;
814  for (unsigned int k=0; k<fk3D[j]->NCell(); k++){
815  dCells.insert(fk3D[j]->Cell(k));
816  }
817  //std::cout<<"There are now: "<<dCells.size()<<" cells in dangler."<<std::endl;
818  double eff, pur;//, matchMetric;
819  double maxEff = -1.0;
820  double maxPur = -1.0;
821  // double maxMatchMetric = 9999.0;
822  int maxPr = -1;
823  double dE, tE, dEMax, tEMax;
824  dEMax = -1.0;
825  tEMax = -1.0;
826  for (unsigned int k=0; k<visPrimaries.size(); k++){
827  std::set<int> temp;
828  temp.insert(visPrimaries[k]->TrackId());
829 
830  eff = bt->HitCollectionEfficiency(temp, fk3D[j]->AllCells(), slices[i]->AllCells(), geo::kXorY, &effMap,true,&dE,&tE);
831  pur = bt->HitCollectionPurity(temp, fk3D[j]->AllCells(), &purMap, &parents,true);
832  // if (pur > 0.0) matchMetric = fabs(tE - dE/pur)/tE;
833  // else matchMetric = 9999.0;
834  //std::cout<<"**********This track matches primary: "<<k<<" with efficiency: "<<eff<<" and metric: "<<matchMetric<<std::endl;
835  if(eff > maxEff){
836  //if(matchMetric < maxMatchMetric){
837  // maxMatchMetric = matchMetric;
838  maxEff = eff;
839  maxPur = pur;
840  maxPr = k;
841  dEMax = dE;
842  tEMax = tE;
843  }
844  else if ((eff == maxEff) && (pur > maxPur)){
845  // maxMatchMetric = matchMetric;
846  maxEff = eff;
847  maxPur = pur;
848  maxPr = k;
849  dEMax = dE;
850  tEMax = tE;
851  }
852  }
853  if (maxEff > 0.0){
854  //if (maxMatchMetric < 9999.0){
855  std::set<int> temp;
856  temp.insert(visPrimaries[maxPr]->TrackId());
857  const sim::Particle* mTest = bt->TrackIDToMotherParticle(visPrimaries[maxPr]->TrackId());
858  danglerMother.push_back(mTest->PdgCode());
859  dangler++;
860  danglerPDG.push_back(visPrimaries[maxPr]->PdgCode());
861  danglerProngIndx.push_back(j);
862  danglerVET.push_back(tEMax);
863  danglerVER.push_back(dEMax);
864  danglerET.push_back(visPrimaries[maxPr]->E());
865  danglerXStartT.push_back(visPrimaries[maxPr]->Vx());
866  danglerYStartT.push_back(visPrimaries[maxPr]->Vy());
867  danglerZStartT.push_back(visPrimaries[maxPr]->Vz());
868  danglerPXT.push_back(visPrimaries[maxPr]->Px());
869  danglerPYT.push_back(visPrimaries[maxPr]->Py());
870  danglerPZT.push_back(visPrimaries[maxPr]->Pz());
871  float costheta = (dX[j]*matchedPXT.back() + dY[j]*matchedPYT.back() + dZ[j]*matchedPZT.back())/
872  (util::pythag(dX[j],dY[j],dZ[j])*util::pythag(matchedPXT.back(),matchedPYT.back()
873  ,matchedPZT.back()));
874  danglerAngleDiff.push_back(acos(costheta)*180.0/M_PI);
875  danglerEff.push_back(maxEff);
876  danglerPur.push_back(maxPur);
877 
878  //fHistoDEffvPur->Fill(maxPur,maxEff);
879 
880  //score += (tEMax*NuEff/NuEVSl)*(tEMax*NuEff/NuEVSl)*(maxEff*maxEff + maxPur*maxPur)*0.5;
881  //std::cout<<"Made a match for reco track: "<<j<<" to particle: "<<maxPr<<" with eff: "<<maxEff<<" and pur: "<< maxPur<<" score is now: "<<score<<std::endl;
882  //fHistoMatchMetric->Fill(maxMatchMetric);
883 
884  }
885  else {
886  //std::cout<<"No match made for reco track: "<<j<<std::endl;
887  maxPur = 0.0;
888  maxEff = 0.0;
889  danglerMother.push_back(-99);
890  dangler++;
891  danglerPDG.push_back(-99);
892  danglerProngIndx.push_back(j);
893  danglerVET.push_back(-1.0);
894  danglerVER.push_back(-1.0);
895  danglerET.push_back(-1.0);
896  danglerXStartT.push_back(-999.0);
897  danglerYStartT.push_back(-999.0);
898  danglerZStartT.push_back(-999.0);
899  danglerPXT.push_back(-999.0);
900  danglerPYT.push_back(-999.0);
901  danglerPZT.push_back(-999.0);
902  danglerAngleDiff.push_back(-1.0);
903  danglerEff.push_back(maxEff);
904  danglerPur.push_back(maxPur);
905 
906  //fHistoDEffvPur->Fill(maxPur,maxEff);
907  }
908  }
909 
910 
911  matchedProngs = matched;
912  unmatchedProngs = unmatched;
913  danglerProngs = dangler;
914  fHistoMatched->Fill(matched);
915  fHistoUnMatched->Fill(unmatched);
916  fHistoLeft->Fill(dangler);
917  //fHistoUvL->Fill(dangler,unmatched);
918 
919  double dEnergy = 0.0;
920  for(std::set<art::Ptr<rb::CellHit>>::iterator dItr = dCells.begin(); dItr!=dCells.end(); dItr++){
921  std::vector<sim::FLSHit> tHits = bt->HitToFLSHit((*dItr));
922  for (unsigned int f=0; f<tHits.size(); ++f){
923  dEnergy += tHits[f].GetEdep();
924  }
925  }
926  score = sqrt(score);
927  //add weight factor to penalize for unmatched reco tracks
928  //std::cout<<dEnergy<<" "<<NuEVSl/NuEff<<std::endl;
929  score = (1.0 - dEnergy*NuEff/NuEVSl)*score;
930  eventScore = score;
931  //fHistoEventScore->Fill(score);
932 
933  //std::cout<<"matched: "<<matched<<" unmatched: "<<unmatched<<" dangler: "<<dangler<<" score: "<<score<<std::endl;
934 
935  //fill tree
936  fOutTree->Fill();
937  }//end loop over slices
938 
939 
940 
941  }
942 
943  //.........................................................................
945  {
946  view.clear();
947  sX.clear();
948  sY.clear();
949  sZ.clear();
950  dX.clear();
951  dY.clear();
952  dZ.clear();
953  ncells.clear();
954  ncellsX.clear();
955  ncellsY.clear();
956  gevTot.clear();
957  totGeV.clear();
958  gevX.clear();
959  gevY.clear();
960  len.clear();
961  matchedPDG.clear();
962  matchedProngIndx.clear();
963  matchedVET.clear();
964  matchedVER.clear();
965  matchedET.clear();
966  matchedXStartT.clear();
967  matchedYStartT.clear();
968  matchedZStartT.clear();
969  matchedPXT.clear();
970  matchedPYT.clear();
971  matchedPZT.clear();
972  matchedAngleDiff.clear();
973  matchedEff.clear();
974  matchedPur.clear();
975  matchedMother.clear();
976  unmatchedMother.clear();
977  unmatchedPDG.clear();
978  unmatchedVET.clear();
979  unmatchedET.clear();
980  unmatchedXStartT.clear();
981  unmatchedYStartT.clear();
982  unmatchedZStartT.clear();
983  unmatchedPXT.clear();
984  unmatchedPYT.clear();
985  unmatchedPZT.clear();
986  matchedProngs = 0;
987  unmatchedProngs = 0;
988  ccnc = -1;
989  NuE = -1.0;
990  NuEff = -1.0;
991  NuPur = -1.0;
992  NuEVSl = -1.0;
993  NuEVTot = -1.0;
994  Vx = -999.0;
995  Vy = -999.0;
996  Vz = -999.0;
997  NuPdg = 0;
998  mode = -1;
999  }
1000 
1001  //.........................................................................
1003  {
1005  int ncells = clus->NCell(v);
1006  float gevTot = 0.0;
1007  for(int j = 0; j < ncells; j++){
1008  rb::RecoHit rhit = clus->RecoHit(v, j);
1009  if(!rhit.IsCalibrated()){
1010  std::cerr << "Not calibrated?!" << std::endl;
1011  continue;
1012  }
1013  gevTot += rhit.GeV()*clus->Weight(v, j);
1014  }
1015  return gevTot;
1016  }
1017 
1018  //.........................................................................
1019  struct Counts
1020  {
1022  {
1023  nHits[0] = nHits[1] = 0;
1024  nOnlyHits[0] = nOnlyHits[1] = 0;
1025  }
1026  int nHits[2]; // Indexed by view
1027  int nOnlyHits[2]; // Hit by only one particle
1028  };
1029 
1030  //......................................................................
1032  std::set<int> ids) {
1033 
1034  std::map<int, Counts> counts;
1036 
1037  // loop over all hits and check if this track id contributed to them
1038  //std::cout<<"Finding visible prongs, preparing to loop over hits"<<std::endl;
1039  for(unsigned int i = 0; i < allhits.size(); ++i){
1040 
1041  //get vector of track id's that contribute to this hit
1042  const art::Ptr<rb::CellHit> hit = allhits[i];
1043  //std::cout<<"Getting hit "<<i<<std::endl;
1044  const std::vector<cheat::TrackIDE> ProngIDs = bt->HitToTrackIDE(hit);
1045 
1046  // loop over all the returned track ids
1047  int nProngs = 0; // number of Prongs that contributed energy to the hit
1048  int loneID = -1; // Which ID it was if there was only one
1049 
1050  for(unsigned int j = 0; j < ProngIDs.size(); ++j){
1051 
1052  const int ProngID = ProngIDs[j].trackID;
1053  const float eFrac = ProngIDs[j].energyFrac;
1054  //detector energy resolution is around 10%
1055  //to be a visible track we are looking for at least 1 hit in each view
1056  //where a true particle deposited 20% (2 sigma) of the energy in that cell
1057  //this is a very minimum definition but leaves the most room for optimization
1058  if(eFrac >= 0.2){
1059  ++nProngs;
1060  loneID = ProngID;
1061 
1062  // check if the current id matches the input Prong id
1063  if(ids.find(ProngID) != ids.end()){
1064  ++counts[ProngID].nHits[hit->View()];
1065  }
1066  }
1067  } // end for j
1068 
1069  if(nProngs == 1){
1070  ++counts[loneID].nOnlyHits[hit->View()];
1071  }
1072 
1073  } // end for i
1074 
1075  const int kMinHits = 1; // minimum number of hits to be considered visible
1076  const int kMinOnlyHit = 0; // minimum number of hits in which it is only particle contributing energy to be visible
1077  std::set<int> ret;
1078 
1079  for(std::set<int>::iterator it = ids.begin(); it != ids.end(); ++it){
1080 
1081  const Counts c = counts[*it];
1082  if((c.nOnlyHits[geo::kX] >= kMinOnlyHit) && (c.nOnlyHits[geo::kY] >= kMinOnlyHit)){
1083  if((c.nHits[geo::kX] >= kMinHits) && (c.nHits[geo::kY] >= kMinHits)){
1084  //std::cout<<"For track: "<<(*it)<<" xcells: "<<c.nHits[geo::kX]<<" ycells "<<c.nHits[geo::kY]<<std::endl;
1085  ret.insert(*it);
1086  }
1087  }
1088 
1089  }
1090  return ret;
1091  } // end FindVisibleProngs
1092 
1093 
1094 } // end namespace numuCC
1095 ////////////////////////////////////////////////////////////////////////
1096 
1097 namespace fuzz
1098 {
1100 }
1101 ////////////////////////////////////////////////////////////////////////
const XML_Char int len
Definition: expat.h:262
std::vector< float > matchedPZT
int fCutType
what cuts to apply (see CAFCutter.h)
std::vector< float > unmatchedYStartT
SubRunNumber_t subRun() const
Definition: Event.h:72
back track the reconstruction to the simulation
int PdgCode() const
Definition: MCParticle.h:211
int CCNC() const
Definition: MCNeutrino.h:148
std::vector< float > danglerXStartT
std::vector< float > dY
std::vector< float > totGeV
double QSqr() const
Definition: MCNeutrino.h:157
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
set< int >::iterator it
const sim::Particle * TrackIDToMotherParticle(int const &id) const
std::vector< TrackIDE > HitToTrackIDE(const rb::CellHit &hit, bool useBirksE=false) const
Convenience function. HitsToTrackIDE but for a single hit.
X or Y views.
Definition: PlaneGeo.h:30
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
std::vector< float > matchedXStartT
std::vector< float > matchedYStartT
std::vector< float > matchedAngleDiff
std::string fVertexLabel
label for the process that made the vertices
geo::View_t View() const
Definition: CellHit.h:41
std::vector< float > danglerPur
std::vector< float > unmatchedPZT
std::vector< float > matchedPur
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
recovalid::CAFCutter fCAFCutter
the CAFCutter helper class
T sqrt(T number)
Definition: d0nt_math.hpp:156
Vertical planes which measure X.
Definition: PlaneGeo.h:28
std::vector< float > unmatchedVET
std::vector< float > unmatchedET
std::vector< float > danglerZStartT
std::vector< float > danglerVET
float GeVCalc(art::Ptr< rb::Prong > clus, geo::View_t v)
OStream cerr
Definition: OStream.cxx:7
void abs(TH1 *hist)
std::vector< float > matchedVET
std::vector< float > len
std::string fSliceLabel
MC information.
std::vector< float > danglerPXT
DEFINE_ART_MODULE(TestTMapFile)
Particle class.
T acos(T number)
Definition: d0nt_math.hpp:54
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...
double dE
std::vector< int > matchedMother
#define M_PI
Definition: SbMath.h:34
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
std::vector< float > danglerYStartT
double HitCollectionEfficiency(const std::set< int > &trackIDs, const std::vector< rb::WeightedHit > &whits, const std::vector< const rb::CellHit * > &allhits, const geo::View_t &view, std::map< int, double > *effMap=0, bool energyEff=false, double *desiredEnergy=0, double *totalEnergy=0, int *desiredHits=0, int *totalHits=0) const
Returns the fraction of all energy in an event from a specific set of Geant4 track IDs that are repre...
Float_t Y
Definition: plot.C:38
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
Float_t Z
Definition: plot.C:38
std::string fTrackLabel
label for the process that made tracks
std::vector< int > ncellsY
std::vector< float > unmatchedPXT
std::vector< float > unmatchedXStartT
std::vector< float > unmatchedZStartT
std::vector< int > matchedPDG
std::vector< float > danglerET
void reconfigure(const fhicl::ParameterSet &pset)
FuzzyKValidate(fhicl::ParameterSet const &pset)
bool passCuts(int cut, const caf::StandardRecord *sr)
Definition: CAFCutter.cxx:37
std::vector< float > matchedVER
int PdgCode(Resonance_t res, int Q)
(resonance id, charge) -> PDG code
double W() const
Definition: MCNeutrino.h:154
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
double Y() const
Definition: MCNeutrino.h:156
Float_t E
Definition: plot.C:20
std::vector< sim::FLSHit > HitToFLSHit(const rb::CellHit &hit) const
All the FLSHits that contributed to this hit, sorted from most to least light.
std::vector< float > gevY
std::string fCAFLabel
label for the process that made standard records
T get(std::string const &key) const
Definition: ParameterSet.h:231
std::vector< float > sZ
rb::RecoHit RecoHit(const art::Ptr< rb::CellHit > &chit) const
Return calibrated hit based on assumed W coordinate.
Definition: Cluster.cxx:259
std::vector< int > danglerMother
void analyze(const art::Event &evt)
std::vector< float > dX
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
Definition: RecoHit.cxx:35
std::vector< int > danglerProngIndx
double X() const
Definition: MCNeutrino.h:155
const double j
Definition: BetheBloch.cxx:29
Fuzzy k-Means prong-finding algorithm.
std::vector< float > sY
std::vector< int > unmatchedMother
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
std::string fAssn3DLabel
association label for 3D prongs
Vertex location in position and time.
std::vector< int > unmatchedPDG
std::vector< float > danglerAngleDiff
Definition: View.py:1
Definition: run.py:1
std::vector< float > sX
size_type size() const
Definition: PtrVector.h:308
double HitCollectionPurity(const std::set< int > &trackIDs, const std::vector< rb::WeightedHit > &whits, std::map< int, double > *purMap=0, std::map< int, int > *parents=0, bool energyPur=false) const
Returns the fraction of hits in a collection that come from the specified Geant4 track ids...
OStream cout
Definition: OStream.cxx:6
std::vector< float > gevTot
std::vector< float > matchedPXT
EventNumber_t event() const
Definition: EventID.h:116
std::vector< TrackIDE > HitsToTrackIDE(const std::vector< const rb::CellHit * > &hits, bool useBirksE=false) const
Returns vector of TrackIDE structs contributing to the given collection of hits.
double Vx(const int i=0) const
Definition: MCParticle.h:220
std::vector< int > danglerPDG
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
std::string fAssn2DLabel
association label for 2D prongs
float GeV() const
Definition: RecoHit.cxx:69
std::vector< float > danglerPYT
bool fApplyCAFCuts
should we apply the CAF-level cuts?
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
std::vector< float > danglerVER
Definition: event.h:1
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:219
std::set< int > FindVisibleProngs3D(const art::PtrVector< rb::CellHit > &allhits, std::set< int > ids)
std::string fCellHitLabel
label for the process that made cell hits
double Vz(const int i=0) const
Definition: MCParticle.h:222
std::vector< int > ncellsX
std::vector< float > matchedPYT
std::vector< float > danglerEff
Helper class for Reco Validation modules.
Definition: CAFCutter.h:58
std::vector< float > unmatchedPYT
bool NeutrinoSet() const
Definition: MCTruth.h:77
std::vector< float > matchedZStartT
Int_t trackID
Definition: plot.C:84
std::vector< float > danglerPZT
std::vector< float > matchedEff
RunNumber_t run() const
Definition: Event.h:77
Float_t X
Definition: plot.C:38
std::vector< int > matchedProngIndx
std::vector< float > matchedET
int Mode() const
Definition: MCNeutrino.h:149
std::vector< float > dZ
Definition: fwd.h:28
EventID id() const
Definition: Event.h:56
double Vy(const int i=0) const
Definition: MCParticle.h:221
double Weight(unsigned int globalIdx) const
Weight assigned to the cell.
Definition: Cluster.cxx:209
std::vector< float > gevX
bool failedToGet() const
Definition: Handle.h:196
std::vector< const sim::Particle * > MCTruthToParticles(art::Ptr< simb::MCTruth > const &mct) const