SigVsBkgSimCompare_module.cc
Go to the documentation of this file.
14 #include "fhiclcpp/ParameterSet.h"
16 
17 #include "HoughTransform.h"
18 #include "NNbarUtilities.h"
19 #include "RecoBase/CellHit.h"
20 #include "RecoBase/Cluster.h"
21 #include "RecoBase/Prong.h"
22 #include "RecoBase/RecoHit.h"
23 #include "RecoBase/Vertex.h"
24 #include "Utilities/AssociationUtil.h"
28 
29 #include <math.h>
30 #include <cmath>
31 #include <iterator>
32 #include <string>
33 #include <vector>
34 
35 #include <TCanvas.h>
36 #include <TColor.h>
37 #include <TGraph.h>
38 #include <TH1.h>
39 #include <TH2.h>
40 #include <TLatex.h>
41 #include <TLegend.h>
42 #include <TMath.h>
43 #include <TROOT.h>
44 #include <TStyle.h>
45 #include <TTree.h>
46 
47 namespace nnbar {
48 class SigVsBkgSimCompare;
49 }
50 
52  public:
53  explicit SigVsBkgSimCompare(fhicl::ParameterSet const &p);
54  void analyze(art::Event const &e) override;
55  void beginJob() override;
56  void endJob() override;
57 
58  private:
60  unsigned int _SampleType;
61 
62  // TTree *_VariableTree;
63  // unsigned int _tRunNumber;
64  // unsigned int _tSubrunNumber;
65  // unsigned int _tEventNumber;
66  // unsigned int _tSliceNumber;
67  // unsigned int _tSampleType;
68  // double _tTotalVisibleEnergy;
69  // double _tTemporalClusterExpand;
70  // double _tTotalHitCountInXView;
71  // double _tTotalHitCountInYView;
72  // double _tHitCountXYDifference;
73  // double _tAverageEnergyPerHitXView;
74  // double _tAverageEnergyPerHitYView;
75  // double _tHitVsVertexTimingCorrelation;
76  // double _tTrueVertexX;
77  // double _tTrueVertexY;
78  // double _tTrueVertexZ;
79  // double _tTrueVertexT;
80  // double _tReconVertexX;
81  // double _tReconVertexY;
82  // double _tReconVertexZ;
83  // double _tReconVertexT;
84 
88 
90 
96 
97  std::vector<std::tuple<double, double, double>> HitsXView, HitsYView;
98 };
99 
101  : EDAnalyzer(p) {
102  _InputShwrLabel = "shower";
103  _SampleType = p.get<unsigned int>("SampleType");
104 }
105 
108 
109  // _VariableTree =
110  // tfs->make<TTree>("VariableTree", "NNbar Analysis Variable Tree");
111  // _VariableTree->Branch("RunNumber", &_tRunNumber);
112  // _VariableTree->Branch("SubrunNumber", &_tSubrunNumber);
113  // _VariableTree->Branch("EventNumber", &_tEventNumber);
114  // _VariableTree->Branch("SliceNumber", &_tSliceNumber);
115  // _VariableTree->Branch("SampleType", &_tSampleType);
116  // _VariableTree->Branch("TotalVisibleEnergy", &_tTotalVisibleEnergy);
117  // _VariableTree->Branch("TotalHitCountInXView", &_tTotalHitCountInXView);
118  // _VariableTree->Branch("TotalHitCountInYView", &_tTotalHitCountInYView);
119  // _VariableTree->Branch("HitCountXYDifference", &_tHitCountXYDifference);
120  // _VariableTree->Branch("TemporalClusterExpand", &_tTemporalClusterExpand);
121  // _VariableTree->Branch("AverageEnergyPerHitXView",
122  // &_tAverageEnergyPerHitXView);
123  // _VariableTree->Branch("AverageEnergyPerHitYView",
124  // &_tAverageEnergyPerHitYView);
125  // _VariableTree->Branch("HitVsVertexTimingCorrelation",
126  // &_tHitVsVertexTimingCorrelation);
127  // _VariableTree->Branch("TrueVertexX", &_tTrueVertexX);
128  // _VariableTree->Branch("TrueVertexY", &_tTrueVertexY);
129  // _VariableTree->Branch("TrueVertexZ", &_tTrueVertexZ);
130  // _VariableTree->Branch("TrueVertexT", &_tTrueVertexT);
131  // _VariableTree->Branch("ReconVertexX", &_tReconVertexX);
132  // _VariableTree->Branch("ReconVertexY", &_tReconVertexY);
133  // _VariableTree->Branch("ReconVertexZ", &_tReconVertexZ);
134  // _VariableTree->Branch("ReconVertexT", &_tReconVertexT);
135 
137  tfs->make<TH1D>("EnergyRatioFromVertexX", ";Energy Ratio;", 100, 0, 1);
139  tfs->make<TH1D>("EnergyRatioFromVertexY", ";Energy Ratio;", 100, 0, 1);
141  tfs->make<TH1D>("EnergyRatioFromVertexZ", ";Energy Ratio;", 100, 0, 1);
142 
143  TimingCorrelation = tfs->make<TH1D>(
144  "TimingCorrelation", ";Timing Correlation Factor;", 100, -1.5, 1.5);
145 
147  tfs->make<TH1D>("DistanceFromReconVertexToTrueVertex",
148  ";Vertex Error (cm);", 100, 0, 200);
150  tfs->make<TH1D>("ZErrorFromReconVertexToTrueVertex",
151  ";Vertex Z Error (cm);", 100, 0, 200);
153  tfs->make<TH2D>("OrientationFromReconVertexToTrueVertex", ";#phi;#theta",
154  50, -TMath::Pi(), TMath::Pi(), 50, 0, TMath::Pi());
156  "PolarAngleFromReconVertexToTrueVertex", ";#phi;", 50, -1, 7);
158  "AzimuthalAngleFromReconVertexToTrueVertex", ";#theta;", 50, -2, 2);
159 }
160 
163  std::vector<art::Ptr<simb::MCTruth>> mctruths;
164  if (e.getByLabel("generator", MCTruthHandle))
165  art::fill_ptr_vector(mctruths, MCTruthHandle);
166 
167  /*
168  double trueVertexX, trueVertexY, trueVertexZ, trueVertexT;
169  for (std::vector<art::Ptr<simb::MCTruth>>::iterator mctruthItr =
170  mctruths.begin();
171  mctruthItr != mctruths.end(); ++mctruthItr) {
172  auto mctruth = *mctruthItr;
173  auto particle = (*mctruth).GetParticle(0);
174  trueVertexX = particle.Vx();
175  trueVertexY = particle.Vy();
176  trueVertexZ = particle.Vz();
177  trueVertexT = particle.T();
178  }
179  */
180 
182  std::vector<art::Ptr<rb::Cluster>> clusters;
183  if (e.getByLabel(_InputShwrLabel, clusterHandle))
184  art::fill_ptr_vector(clusters, clusterHandle);
185 
186  unsigned int clusterIdx = 0;
187  for (std::vector<art::Ptr<rb::Cluster>>::iterator clusterItr =
188  clusters.begin();
189  clusterItr != clusters.end(); ++clusterItr) {
190  auto cluster = *clusterItr;
191 
192  art::PtrVector<rb::CellHit> allCells = (*cluster).AllCells();
193 
194  unsigned int countHits = 0;
195  for (art::PtrVector<rb::CellHit>::iterator itrCell = allCells.begin();
196  itrCell != allCells.end(); itrCell++) {
197  countHits++;
198  }
199  if (countHits == 0) {
200  std::cout << "countHits: " << countHits << std::endl << std::endl << std::endl << std::endl << std::endl << std::endl << std::endl << std::endl;
201  }
202 /*
203  // CONSIDERING THE ENERGY DISTRIBUTION TOPOLOGY
204  double totalE = (*cluster).TotalGeV();
205  double totalEX = 0.;
206  double totalEY = 0.;
207  double totalEnergyOfHitsAboveVertexX = 0.;
208  double totalEnergyOfHitsAboveVertexY = 0.;
209  double totalEnergyOfHitsAboveVertexZ = 0.;
210 
211  // Consider timing correlation
212  std::multimap<double, double> distanceToVertex_relTimingToVertex_Map;
213  unsigned int hitIdx = 0;
214 
215  // Compare distance between true vertex and reconstructed vertex
216  std::multimap<double, double> HitsZ, HitsX, HitsY;
217 
218  for (art::PtrVector<rb::CellHit>::iterator itrCell = allCells.begin();
219  itrCell != allCells.end(); itrCell++) {
220  if ((*itrCell)->View() == geo::View_t::kXorY) continue;
221 
222  bool cIsYView = (*itrCell)->View() == geo::View_t::kY ? true : false;
223  rb::RecoHit recoHit = (*cluster).RecoHit(*itrCell);
224  if (!recoHit.IsCalibrated()) continue;
225 
226  double cX = recoHit.X();
227  double cY = recoHit.Y();
228  double cZ = recoHit.Z();
229  double cT = recoHit.T();
230  double cE = recoHit.GeV();
231 
232  HitsZ.insert(std::pair<double, double>(cZ, cE));
233  if (!cIsYView) {
234  HitsX.insert(std::pair<double, double>(cX, cE));
235  totalEX += cE;
236  } else {
237  HitsY.insert(std::pair<double, double>(cY, cE));
238  totalEY += cE;
239  }
240 
241  if (!cIsYView) {
242  if (cX > trueVertexX) totalEnergyOfHitsAboveVertexX += cE;
243  } else {
244  if (cY > trueVertexY) totalEnergyOfHitsAboveVertexY += cE;
245  }
246  if (cZ > trueVertexZ) totalEnergyOfHitsAboveVertexZ += cE;
247 
248  // COUNTING THE NUMBER OF PRONGS
249  if (!cIsYView) {
250  HitsXView.push_back(std::make_tuple(cZ, cX, cE));
251  } else {
252  HitsYView.push_back(std::make_tuple(cZ, cY, cE));
253  }
254 
255  hitIdx++;
256  }
257 
258  // Compare distance between true vertex and reconstructed vertex
259  double lowerBoundOfMedianZVertex = 0.;
260  double upperBoundOfMedianZVertex = 0.;
261  double halfOfTotalE = 0.;
262  for (auto hit = HitsZ.begin(); hit != HitsZ.end(); hit++) {
263  halfOfTotalE += (*hit).second;
264  if (halfOfTotalE < totalE / 2) {
265  lowerBoundOfMedianZVertex = (*hit).first;
266  upperBoundOfMedianZVertex = (*(hit++)).first;
267  hit--;
268  }
269  }
270  double reconVertexZ =
271  (lowerBoundOfMedianZVertex + upperBoundOfMedianZVertex) / 2.;
272 
273  double lowerBoundOfMedianXVertex = 0.;
274  double upperBoundOfMedianXVertex = 0.;
275  double halfOfTotalEX = 0.;
276  for (auto hit = HitsX.begin(); hit != HitsX.end(); hit++) {
277  halfOfTotalEX += (*hit).second;
278  if (halfOfTotalEX < totalEX / 2) {
279  lowerBoundOfMedianXVertex = (*hit).first;
280  upperBoundOfMedianXVertex = (*(hit++)).first;
281  hit--;
282  }
283  }
284  double reconVertexX =
285  (lowerBoundOfMedianXVertex + upperBoundOfMedianXVertex) / 2.;
286 
287  double lowerBoundOfMedianYVertex = 0.;
288  double upperBoundOfMedianYVertex = 0.;
289  double halfOfTotalEY = 0.;
290  for (auto hit = HitsY.begin(); hit != HitsY.end(); hit++) {
291  halfOfTotalEY += (*hit).second;
292  if (halfOfTotalEY < totalEY / 2) {
293  lowerBoundOfMedianYVertex = (*hit).first;
294  upperBoundOfMedianYVertex = (*(hit++)).first;
295  hit--;
296  }
297  }
298  double reconVertexY =
299  (lowerBoundOfMedianYVertex + upperBoundOfMedianYVertex) / 2.;
300 
301  for (art::PtrVector<rb::CellHit>::iterator itrCell = allCells.begin();
302  itrCell != allCells.end(); itrCell++) {
303  if ((*itrCell)->View() == geo::View_t::kXorY) continue;
304 
305  rb::RecoHit recoHit = (*cluster).RecoHit(*itrCell);
306  if (!recoHit.IsCalibrated()) continue;
307 
308  double cX = recoHit.X();
309  double cY = recoHit.Y();
310  double cZ = recoHit.Z();
311  double cT = recoHit.T();
312 
313  // Consider timing correlation
314  // Timing correlation with recon vertex
315  double distanceToVertex = TMath::Sqrt(TMath::Power(cX - reconVertexX, 2) +
316  TMath::Power(cY - reconVertexY, 2) +
317  TMath::Power(cZ - reconVertexZ, 2));
318  double relTimingToVertex = cT;
319  distanceToVertex_relTimingToVertex_Map.insert(
320  std::pair<double, double>(distanceToVertex, relTimingToVertex));
321  }
322 
323  //--------------------------------------------------------------------------
324  // Consider timing correlation
325  // Recon Vertex
326  const unsigned int mapSize = distanceToVertex_relTimingToVertex_Map.size();
327  unsigned int countElem = 0;
328  double *distanceToVertex_array = new double[10];
329  double *relTimingToVertex_array = new double[10];
330  for (auto &elem : distanceToVertex_relTimingToVertex_Map) {
331  if (4 < countElem && countElem < mapSize - 5) continue;
332  *(distanceToVertex_array + countElem) = elem.first;
333  *(relTimingToVertex_array + countElem) = elem.second;
334  countElem++;
335  }
336  // The timing of reconstructed vertex is calculated as the mean timing of
337  // the 5 hits nearest to the vertex.
338  // So this loop find that average timing.
339  double reconVertexT = 0.;
340  for (unsigned int j = 0; j < 5; j++) {
341  reconVertexT += *(relTimingToVertex_array + j);
342  }
343  reconVertexT = reconVertexT / 5;
344  for (unsigned int j = 0; j < 5; j++) {
345  *(relTimingToVertex_array + countElem) =
346  *(relTimingToVertex_array + countElem) - reconVertexT;
347  }
348  TGraph *timingCorrelation =
349  new TGraph(10, distanceToVertex_array, relTimingToVertex_array);
350  double correlationFactor = timingCorrelation->GetCorrelationFactor();
351  TimingCorrelation->Fill(correlationFactor);
352 
353  //--------------------------------------------------------------------------
354 
355  double diff_true_recon_vertex_x = reconVertexX - trueVertexX;
356  double diff_true_recon_vertex_y = reconVertexY - trueVertexY;
357  double diff_true_recon_vertex_z = reconVertexZ - trueVertexZ;
358  double distFromReconVertToTrueVert =
359  TMath::Sqrt(TMath::Power(diff_true_recon_vertex_x, 2) +
360  TMath::Power(diff_true_recon_vertex_y, 2) +
361  TMath::Power(diff_true_recon_vertex_z, 2));
362  double thetaOrientationFromReconVertToTrueVert =
363  TMath::ATan(diff_true_recon_vertex_y / diff_true_recon_vertex_x);
364  double phiOrientationFromReconVertToTrueVert =
365  TMath::ACos(diff_true_recon_vertex_z / distFromReconVertToTrueVert);
366  DistanceFromReconVertexToTrueVertex->Fill(distFromReconVertToTrueVert);
367  OrientationFromReconVertexToTrueVertex->Fill(
368  phiOrientationFromReconVertToTrueVert,
369  thetaOrientationFromReconVertToTrueVert);
370  PolarAngleFromReconVertexToTrueVertex->Fill(
371  phiOrientationFromReconVertToTrueVert);
372  AzimuthalAngleFromReconVertexToTrueVertex->Fill(
373  thetaOrientationFromReconVertToTrueVert);
374 
375  // CONSIDERING THE ENERGY DISTRIBUTION TOPOLOGY
376  EnergyRatioFromVertexX->Fill(totalEnergyOfHitsAboveVertexX / totalE);
377  EnergyRatioFromVertexY->Fill(totalEnergyOfHitsAboveVertexY / totalE);
378  EnergyRatioFromVertexZ->Fill(totalEnergyOfHitsAboveVertexZ / totalE);
379 
380  // Filling analysis tree
381  // _tRunNumber = e.run();
382  // _tSubrunNumber = e.subRun();
383  // _tEventNumber = e.id().event();
384  // _tSliceNumber = std::distance(clusters.begin(), clusterItr);
385  // _tSampleType = _SampleType;
386  // double totVisEX =
387  // nnbar::NNbarUtilities::totalVisibleEnergyInView(*cluster, geo::kX);
388  // double totVisEY =
389  // nnbar::NNbarUtilities::totalVisibleEnergyInView(*cluster, geo::kY);
390  // _tTotalVisibleEnergy = totVisEX + totVisEY;
391 
392  ZErrorFromReconVertexToTrueVertex->Fill(
393  TMath::Abs(diff_true_recon_vertex_z));
394 
395  // _tTemporalClusterExpand = nnbar::NNbarUtilities::temporalExpand(*cluster);
396  // _tTotalHitCountInXView =
397  // nnbar::NNbarUtilities::hitCountInView(*cluster, geo::kX);
398  // _tTotalHitCountInYView =
399  // nnbar::NNbarUtilities::hitCountInView(*cluster, geo::kY);
400  // _tHitCountXYDifference = _tTotalHitCountInYView - _tTotalHitCountInXView;
401  // _tAverageEnergyPerHitXView = totVisEX / _tTotalHitCountInXView;
402  // _tAverageEnergyPerHitYView = totVisEY / _tTotalHitCountInYView;
403  // _tHitVsVertexTimingCorrelation = correlationFactor;
404  // _tReconVertexX = reconVertexX;
405  // _tReconVertexY = reconVertexY;
406  // _tReconVertexZ = reconVertexZ;
407  // _tReconVertexT = reconVertexT;
408  // if (_tSampleType != 0) {
409  // _tTrueVertexX = trueVertexX;
410  // _tTrueVertexY = trueVertexY;
411  // _tTrueVertexZ = trueVertexZ;
412  // _tTrueVertexT = trueVertexT;
413  // } else {
414  // _tTrueVertexX = -1.0E9;
415  // _tTrueVertexY = -1.0E9;
416  // _tTrueVertexZ = -1.0E9;
417  // _tTrueVertexT = -1.0E9;
418  // }
419  // _VariableTree->Fill();
420 */
421  clusterIdx++;
422  }
423 
424  HitsXView.clear();
425  HitsYView.clear();
426 }
427 
429  TCanvas *canvas = new TCanvas("canvas", "Canvas", 2400, 800);
430  canvas->Divide(3, 1);
431  canvas->cd(1);
432  EnergyRatioFromVertexX->SetTitle("x");
433  EnergyRatioFromVertexX->Draw();
434  canvas->cd(2);
435  EnergyRatioFromVertexY->SetTitle("y");
436  EnergyRatioFromVertexY->Draw();
437  canvas->cd(3);
438  EnergyRatioFromVertexZ->SetTitle("z");
439  EnergyRatioFromVertexZ->Draw();
440  canvas->SaveAs(Form("EnergyVsTopology.pdf"));
441 
442  TCanvas *canvas_timing =
443  new TCanvas("canvas_timing", "Canvas Timing", 1200, 1200);
444  TimingCorrelation->GetXaxis()->SetTitle("Timing Correlation Factor");
445  TimingCorrelation->GetXaxis()->CenterTitle();
446  TimingCorrelation->GetYaxis()->SetTitle("");
447  TimingCorrelation->GetYaxis()->CenterTitle();
448  TimingCorrelation->Draw("HIST");
449  canvas_timing->SaveAs(Form("TimingCorrelation.pdf"));
450 
451  TCanvas *canvas_vertex_error =
452  new TCanvas("canvas_vertex_error", "Canvas Vertex Error", 2400, 800);
453  canvas_vertex_error->Divide(3, 1);
454  canvas_vertex_error->cd(1);
455  DistanceFromReconVertexToTrueVertex->GetXaxis()->SetTitle(
456  "#deltar_{equale-(v)} (cm)");
457  DistanceFromReconVertexToTrueVertex->GetXaxis()->CenterTitle();
458  DistanceFromReconVertexToTrueVertex->GetYaxis()->SetTitle("");
459  DistanceFromReconVertexToTrueVertex->GetYaxis()->CenterTitle();
461  canvas_vertex_error->cd(2);
462  PolarAngleFromReconVertexToTrueVertex->GetXaxis()->SetTitle(
463  "#phi_{equale-(v)} (rad)");
464  PolarAngleFromReconVertexToTrueVertex->GetXaxis()->CenterTitle();
465  PolarAngleFromReconVertexToTrueVertex->GetYaxis()->SetTitle("");
466  PolarAngleFromReconVertexToTrueVertex->GetYaxis()->CenterTitle();
468  canvas_vertex_error->cd(3);
469  AzimuthalAngleFromReconVertexToTrueVertex->GetXaxis()->SetTitle(
470  "#theta_{equale-(v)} (rad)");
471  AzimuthalAngleFromReconVertexToTrueVertex->GetXaxis()->CenterTitle();
472  AzimuthalAngleFromReconVertexToTrueVertex->GetYaxis()->SetTitle("");
473  AzimuthalAngleFromReconVertexToTrueVertex->GetYaxis()->CenterTitle();
475  canvas_vertex_error->SaveAs(Form("VertexError.pdf"));
476 
477  TCanvas *canvas_zvertex_error =
478  new TCanvas("canvas_zvertex_error", "Canvas Z Vertex Error", 2400, 1200);
479  ZErrorFromReconVertexToTrueVertex->GetXaxis()->SetTitle(
480  "#deltaz_{equale-(v)} (cm)");
481  ZErrorFromReconVertexToTrueVertex->GetXaxis()->CenterTitle();
482  ZErrorFromReconVertexToTrueVertex->GetYaxis()->SetTitle("");
483  ZErrorFromReconVertexToTrueVertex->GetYaxis()->CenterTitle();
484  ZErrorFromReconVertexToTrueVertex->Draw("HIST");
485  canvas_zvertex_error->SaveAs(Form("ZVertexError.pdf"));
486 }
487 
void analyze(art::Event const &e) override
iterator begin()
Definition: PtrVector.h:223
const char * p
Definition: xmltok.h:285
std::vector< std::tuple< double, double, double > > HitsYView
DEFINE_ART_MODULE(TestTMapFile)
Particle class.
std::vector< std::tuple< double, double, double > > HitsXView
T get(std::string const &key) const
Definition: ParameterSet.h:231
iterator end()
Definition: PtrVector.h:237
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
Vertex location in position and time.
data_t::iterator iterator
Definition: PtrVector.h:60
OStream cout
Definition: OStream.cxx:6
T * make(ARGS...args) const
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
SigVsBkgSimCompare(fhicl::ParameterSet const &p)
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:464
Float_t e
Definition: plot.C:35