RecoAnalysis_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 <cmath>
30 #include <iterator>
31 #include <math.h>
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 // struct Point2D {
48 // point2d_t point;
49 //
50 // Point2D(double xa, double ya) {
51 // point.first = xa;
52 // point.second = ya;
53 // }
54 //
55 // struct Point2D &operator+=(const Point2D &rhs) {
56 // point.first += rhs.point.first;
57 // point.second += rhs.point.second;
58 // return *this;
59 // }
60 //
61 // struct Point2D &operator-=(const Point2D &rhs) {
62 // point.first -= rhs.point.first;
63 // point.second -= rhs.point.second;
64 // return *this;
65 // }
66 //
67 // struct Point2D &operator*=(const double &k) {
68 // point.first *= k;
69 // point.second *= k;
70 // return *this;
71 // }
72 //
73 // double operator*=(const Point2D &rhs) {
74 // return point.first * rhs.point.first + point.second * rhs.point.second;
75 // }
76 //
77 // struct Point2D &operator=(const Point2D &rhs) {
78 // point.first = rhs.point.first;
79 // point.second = rhs.point.second;
80 // return *this;
81 // }
82 // };
83 //
84 // Point2D operator+(Point2D lhs, const Point2D &rhs) { return lhs += rhs; }
85 // Point2D operator-(Point2D lhs, const Point2D &rhs) { return lhs -= rhs; }
86 // double operator*(Point2D lhs, const Point2D &rhs) { return lhs *= rhs; }
87 // Point2D operator*(const double k, Point2D rhs) { return rhs *= k; }
88 // Point2D operator*(Point2D lhs, const double k) { return lhs *= k; }
89 //
90 // double length(Point2D A) {
91 // return TMath::Sqrt(TMath::Power(A.point.first, 2) + TMath::Power(A.point.second, 2));
92 // }
93 //
94 // unsigned int check_isolation(double* array, unsigned int maxIndex, unsigned int hitIndex) {
95 // if (TMath::Abs(*(array + hitIndex + 1) - *(array + hitIndex)) > 5.) {
96 // return hitIndex;
97 // } else {
98 // if (hitIndex + 1 == maxIndex) {
99 // return hitIndex;
100 // } else {
101 // return check_isolation(array, maxIndex, hitIndex + 1);
102 // }
103 // }
104 // }
105 
106 namespace nnbar {
107 class RecoAnalysis;
108 }
109 
111 public:
112  explicit RecoAnalysis(fhicl::ParameterSet const &p);
113  void beginJob() override;
114  void endJob() override;
115  void analyze(art::Event const &e) override;
116  TGraph *PlotVectorPoint(std::vector<point2d_t> pts);
117  double HozVertex(
118  std::vector<std::tuple<double, double, double>> HitCollectionSingleView,
119  double vertexZ);
121  void evd(art::Event const &e, unsigned int clusterIdx);
122 
123 private:
126  unsigned int _SampleType;
127 
129  unsigned int _tRunNumber;
130  unsigned int _tSubrunNumber;
131  unsigned int _tEventNumber;
132  unsigned int _tSliceNumber;
133  unsigned int _tSampleType;
154 
158 
160 
166 
172 
173  std::vector<std::tuple<double, double, double>> HitsXView, HitsYView;
174 };
175 
177  : EDAnalyzer(p) {
178  _InputShwrLabel = "shower:showers";
179  _SampleType = p.get<unsigned int>("SampleType");
180 }
181 
182 // double nnbar::RecoAnalysis::HozVertex(
183 // std::vector<std::tuple<double, double, double>> HitCollectionSingleView,
184 // double vertexZ) {
185 // std::vector<Point2D> Hits;
186 // for (auto &hit : HitCollectionSingleView) {
187 // double x = std::get<0>(hit);
188 // double y = std::get<1>(hit);
189 // Point2D aHit(x, y);
190 // Hits.push_back(aHit);
191 // }
192 //
193 // return 0;
194 // }
195 
197  _ExposureTimeInUSec = 0.;
198 
200 
201  _VariableTree =
202  tfs->make<TTree>("VariableTree", "NNbar Analysis Variable Tree");
203  _VariableTree->Branch("RunNumber", &_tRunNumber);
204  _VariableTree->Branch("SubrunNumber", &_tSubrunNumber);
205  _VariableTree->Branch("EventNumber", &_tEventNumber);
206  _VariableTree->Branch("SliceNumber", &_tSliceNumber);
207  _VariableTree->Branch("SampleType", &_tSampleType);
208  _VariableTree->Branch("TotalVisibleEnergy", &_tTotalVisibleEnergy);
209  _VariableTree->Branch("TotalHitCountInXView", &_tTotalHitCountInXView);
210  _VariableTree->Branch("TotalHitCountInYView", &_tTotalHitCountInYView);
211  _VariableTree->Branch("HitCountXYDifference", &_tHitCountXYDifference);
212  _VariableTree->Branch("TemporalClusterExpand", &_tTemporalClusterExpand);
213  _VariableTree->Branch("AverageEnergyPerHitXView",
215  _VariableTree->Branch("AverageEnergyPerHitYView",
217  _VariableTree->Branch("HitVsVertexTimingCorrelation",
219  _VariableTree->Branch("TrueVertexX", &_tTrueVertexX);
220  _VariableTree->Branch("TrueVertexY", &_tTrueVertexY);
221  _VariableTree->Branch("TrueVertexZ", &_tTrueVertexZ);
222  _VariableTree->Branch("TrueVertexT", &_tTrueVertexT);
223  _VariableTree->Branch("ReconVertexX", &_tReconVertexX);
224  _VariableTree->Branch("ReconVertexY", &_tReconVertexY);
225  _VariableTree->Branch("ReconVertexZ", &_tReconVertexZ);
226  _VariableTree->Branch("ReconVertexT", &_tReconVertexT);
227  _VariableTree->Branch("EArmVertexX", &_tEArmVertexX);
228  _VariableTree->Branch("EArmVertexY", &_tEArmVertexY);
229  _VariableTree->Branch("EArmVertexZ", &_tEArmVertexZ);
230  _VariableTree->Branch("EArmVertexT", &_tEArmVertexT);
231 
233  tfs->make<TH1D>("EnergyRatioFromVertexX", ";Energy Ratio;", 100, 0, 1);
235  tfs->make<TH1D>("EnergyRatioFromVertexY", ";Energy Ratio;", 100, 0, 1);
237  tfs->make<TH1D>("EnergyRatioFromVertexZ", ";Energy Ratio;", 100, 0, 1);
238 
239  TimingCorrelation = tfs->make<TH1D>(
240  "TimingCorrelation", ";Timing Correlation Factor;", 100, -1.5, 1.5);
241 
243  tfs->make<TH1D>("DistanceFromReconVertexToTrueVertex",
244  ";Vertex Error (cm);", 100, 0, 200);
246  tfs->make<TH1D>("ZErrorFromReconVertexToTrueVertex",
247  ";Vertex Z Error (cm);", 100, 0, 200);
249  tfs->make<TH2D>("OrientationFromReconVertexToTrueVertex", ";#phi;#theta",
250  50, -TMath::Pi(), TMath::Pi(), 50, 0, TMath::Pi());
252  "PolarAngleFromReconVertexToTrueVertex", ";#phi;", 50, -1, 7);
254  "AzimuthalAngleFromReconVertexToTrueVertex", ";#theta;", 50, -2, 2);
255 
257  tfs->make<TH1D>("DistanceFromElasticArmVertexToTrueVertex",
258  ";EArm Vertex Error (cm);", 100, 0, 200);
260  tfs->make<TH1D>("ZErrorFromElasticArmVertexToTrueVertex",
261  ";Vertex Z Error (cm);", 100, 0, 200);
263  "OrientationFromElasticArmVertexToTrueVertex", ";#phi;#theta", 50,
264  -TMath::Pi(), TMath::Pi(), 50, 0, TMath::Pi());
266  "PolarAngleFromElasticArmVertexToTrueVertex", ";#phi;", 50, -1, 7);
268  "AzimuthalAngleFromElasticArmVertexToTrueVertex", ";#theta;", 50, -2, 2);
269 }
270 
271 TGraph *nnbar::RecoAnalysis::PlotVectorPoint(std::vector<point2d_t> pts) {
272  double *x = new double[pts.size()];
273  double *y = new double[pts.size()];
274 
275  size_t i = 0;
276  for (auto &pt : pts) {
277  *(x + i) = pt.first;
278  *(y + i) = pt.second;
279  i++;
280  }
281 
282  TGraph *graph = new TGraph(pts.size(), x, y);
283  return graph;
284 }
285 
287  point2d_t pB) {
288  double xa = pA.first;
289  double ya = pA.second;
290  double xb = pB.first;
291  double yb = pB.second;
292 
293  return TMath::ATan((ya - yb) / (xa - xb));
294 }
295 
297 
298  _ExposureTimeInUSec += 50;
299 
301  std::vector<art::Ptr<simb::MCTruth>> mctruths;
302  if (e.getByLabel("generator", MCTruthHandle))
303  art::fill_ptr_vector(mctruths, MCTruthHandle);
304 
305  double trueVertexX = 0, trueVertexY = 0, trueVertexZ = 0, trueVertexT = 0;
306  for (std::vector<art::Ptr<simb::MCTruth>>::iterator mctruthItr =
307  mctruths.begin();
308  mctruthItr != mctruths.end(); ++mctruthItr) {
309  auto mctruth = *mctruthItr;
310  auto particle = (*mctruth).GetParticle(0);
311  trueVertexX = particle.Vx();
312  trueVertexY = particle.Vy();
313  trueVertexZ = particle.Vz();
314  trueVertexT = particle.T();
315  }
316 
318  std::vector<art::Ptr<rb::Cluster>> clusters;
319  if (e.getByLabel(_InputShwrLabel, clusterHandle))
320  art::fill_ptr_vector(clusters, clusterHandle);
321 
322  art::FindManyP<rb::Vertex> assnClusterVert(clusterHandle, e, "elasticarmshs");
323 
324  // double elasticArmVertT;
325  double elasticArmVertX;
326  double elasticArmVertY;
327  double elasticArmVertZ;
328  unsigned int clusterIdx = 0;
329  for (std::vector<art::Ptr<rb::Cluster>>::iterator clusterItr =
330  clusters.begin();
331  clusterItr != clusters.end(); ++clusterItr) {
332  auto cluster = *clusterItr;
333 
334  art::PtrVector<rb::CellHit> allCells = (*cluster).AllCells();
335 
336  // Conventional Vertex Reconstruction
337  // elasticArmVertT = -1E9;
338  elasticArmVertX = -1E9;
339  elasticArmVertY = -1E9;
340  elasticArmVertZ = -1E9;
341  if (assnClusterVert.isValid()) {
342  std::vector<art::Ptr<rb::Vertex>> vert = assnClusterVert.at(clusterIdx);
343  if (vert.size() != 1)
344  continue;
345  // elasticArmVertT = vert.at(0)->GetT();
346  elasticArmVertX = vert.at(0)->GetX();
347  elasticArmVertY = vert.at(0)->GetY();
348  elasticArmVertZ = vert.at(0)->GetZ();
349  }
350 
351  //--------------------------------------------------------------------------
352 
353  double diff_true_elasticArm_vertex_x = elasticArmVertX - trueVertexX;
354  double diff_true_elasticArm_vertex_y = elasticArmVertY - trueVertexY;
355  double diff_true_elasticArm_vertex_z = elasticArmVertZ - trueVertexZ;
356  double distFromElasticArmVertToTrueVert =
357  TMath::Sqrt(TMath::Power(diff_true_elasticArm_vertex_x, 2) +
358  TMath::Power(diff_true_elasticArm_vertex_y, 2) +
359  TMath::Power(diff_true_elasticArm_vertex_z, 2));
360  double thetaOrientationFromElasticArmVertToTrueVert = TMath::ATan(
361  diff_true_elasticArm_vertex_y / diff_true_elasticArm_vertex_x);
362  double phiOrientationFromElasticArmVertToTrueVert = TMath::ACos(
363  diff_true_elasticArm_vertex_z / distFromElasticArmVertToTrueVert);
364 
366  distFromElasticArmVertToTrueVert);
368  phiOrientationFromElasticArmVertToTrueVert,
369  thetaOrientationFromElasticArmVertToTrueVert);
371  phiOrientationFromElasticArmVertToTrueVert);
373  thetaOrientationFromElasticArmVertToTrueVert);
374 
375  // CONSIDERING THE ENERGY DISTRIBUTION TOPOLOGY
376  double totalE = (*cluster).TotalGeV();
377  double totalEX = 0.;
378  double totalEY = 0.;
379  double totalEnergyOfHitsAboveVertexX = 0.;
380  double totalEnergyOfHitsAboveVertexY = 0.;
381  double totalEnergyOfHitsAboveVertexZ = 0.;
382 
383  // Consider timing correlation
384  std::multimap<double, double> distanceToVertex_relTimingToVertex_Map;
385  unsigned int hitIdx = 0;
386 
387  // Compare distance between true vertex and reconstructed vertex
388  std::multimap<double, double> HitsZ, HitsX, HitsY;
389 
390  for (art::PtrVector<rb::CellHit>::iterator itrCell = allCells.begin();
391  itrCell != allCells.end(); itrCell++) {
392  if ((*itrCell)->View() == geo::View_t::kXorY)
393  continue;
394 
395  bool cIsYView = (*itrCell)->View() == geo::View_t::kY ? true : false;
396  rb::RecoHit recoHit = (*cluster).RecoHit(*itrCell);
397  if (!recoHit.IsCalibrated())
398  continue;
399 
400  double cX = recoHit.X();
401  double cY = recoHit.Y();
402  double cZ = recoHit.Z();
403  // double cT = recoHit.T();
404  double cE = recoHit.GeV();
405 
406  HitsZ.insert(std::pair<double, double>(cZ, cE));
407  if (!cIsYView) {
408  HitsX.insert(std::pair<double, double>(cX, cE));
409  totalEX += cE;
410  } else {
411  HitsY.insert(std::pair<double, double>(cY, cE));
412  totalEY += cE;
413  }
414 
415  // // Consider timing correlation
416  // // Timing correlation with true vertex
417  // double distanceToVertex = TMath::Sqrt(TMath::Power(cX - trueVertexX, 2)
418  // +
419  // TMath::Power(cY - trueVertexY, 2)
420  // + TMath::Power(cZ - trueVertexY,
421  // 2));
422  // double relTimingToVertex = cT - trueVertexT;
423  // distanceToVertex_relTimingToVertex_Map.insert(
424  // std::pair<double, double>(distanceToVertex, relTimingToVertex));
425 
426  if (!cIsYView) {
427  if (cX > trueVertexX)
428  totalEnergyOfHitsAboveVertexX += cE;
429  } else {
430  if (cY > trueVertexY)
431  totalEnergyOfHitsAboveVertexY += cE;
432  }
433  if (cZ > trueVertexZ)
434  totalEnergyOfHitsAboveVertexZ += cE;
435 
436  // COUNTING THE NUMBER OF PRONGS
437  if (!cIsYView) {
438  HitsXView.push_back(std::make_tuple(cZ, cX, cE));
439  } else {
440  HitsYView.push_back(std::make_tuple(cZ, cY, cE));
441  }
442 
443  hitIdx++;
444  }
445 
446  // Compare distance between true vertex and reconstructed vertex
447  double lowerBoundOfMedianZVertex = 0.;
448  double upperBoundOfMedianZVertex = 0.;
449  double halfOfTotalE = 0.;
450  for (auto hit = HitsZ.begin(); hit != HitsZ.end(); hit++) {
451  halfOfTotalE += (*hit).second;
452  if (halfOfTotalE < totalE / 2) {
453  lowerBoundOfMedianZVertex = (*hit).first;
454  upperBoundOfMedianZVertex = (*(hit++)).first;
455  hit--;
456  }
457  }
458  double reconVertexZ =
459  (lowerBoundOfMedianZVertex + upperBoundOfMedianZVertex) / 2.;
460 
461  double lowerBoundOfMedianXVertex = 0.;
462  double upperBoundOfMedianXVertex = 0.;
463  double halfOfTotalEX = 0.;
464  for (auto hit = HitsX.begin(); hit != HitsX.end(); hit++) {
465  halfOfTotalEX += (*hit).second;
466  if (halfOfTotalEX < totalEX / 2) {
467  lowerBoundOfMedianXVertex = (*hit).first;
468  upperBoundOfMedianXVertex = (*(hit++)).first;
469  hit--;
470  }
471  }
472  double reconVertexX =
473  (lowerBoundOfMedianXVertex + upperBoundOfMedianXVertex) / 2.;
474 
475  double lowerBoundOfMedianYVertex = 0.;
476  double upperBoundOfMedianYVertex = 0.;
477  double halfOfTotalEY = 0.;
478  for (auto hit = HitsY.begin(); hit != HitsY.end(); hit++) {
479  halfOfTotalEY += (*hit).second;
480  if (halfOfTotalEY < totalEY / 2) {
481  lowerBoundOfMedianYVertex = (*hit).first;
482  upperBoundOfMedianYVertex = (*(hit++)).first;
483  hit--;
484  }
485  }
486  double reconVertexY =
487  (lowerBoundOfMedianYVertex + upperBoundOfMedianYVertex) / 2.;
488 
489  for (art::PtrVector<rb::CellHit>::iterator itrCell = allCells.begin();
490  itrCell != allCells.end(); itrCell++) {
491  if ((*itrCell)->View() == geo::View_t::kXorY)
492  continue;
493 
494  rb::RecoHit recoHit = (*cluster).RecoHit(*itrCell);
495  if (!recoHit.IsCalibrated())
496  continue;
497 
498  double cX = recoHit.X();
499  double cY = recoHit.Y();
500  double cZ = recoHit.Z();
501  double cT = recoHit.T();
502 
503  // Consider timing correlation
504  // Timing correlation with recon vertex
505  double distanceToVertex = TMath::Sqrt(TMath::Power(cX - reconVertexX, 2) +
506  TMath::Power(cY - reconVertexY, 2) +
507  TMath::Power(cZ - reconVertexZ, 2));
508  double relTimingToVertex = cT;
509  distanceToVertex_relTimingToVertex_Map.insert(
510  std::pair<double, double>(distanceToVertex, relTimingToVertex));
511  }
512 
513  //--------------------------------------------------------------------------
514  // Consider timing correlation
515  // Recon Vertex
516  const unsigned int mapSize = distanceToVertex_relTimingToVertex_Map.size();
517  unsigned int countElem = 0;
518  double *distanceToVertex_array = new double[10];
519  double *relTimingToVertex_array = new double[10];
520  for (auto &elem : distanceToVertex_relTimingToVertex_Map) {
521  if (4 < countElem && countElem < mapSize - 5)
522  continue;
523  *(distanceToVertex_array + countElem) = elem.first;
524  *(relTimingToVertex_array + countElem) = elem.second;
525  countElem++;
526  }
527  // The timing of reconstructed vertex is calculated as the mean timing of
528  // the 5 hits nearest to the vertex.
529  // So this loop find that average timing.
530  double reconVertexT = 0.;
531  for (unsigned int j = 0; j < 5; j++) {
532  reconVertexT += *(relTimingToVertex_array + j);
533  }
534  reconVertexT = reconVertexT / 5;
535  for (unsigned int j = 0; j < 5; j++) {
536  *(relTimingToVertex_array + countElem) =
537  *(relTimingToVertex_array + countElem) - reconVertexT;
538  }
539  TGraph *timingCorrelation =
540  new TGraph(10, distanceToVertex_array, relTimingToVertex_array);
541  double correlationFactor = timingCorrelation->GetCorrelationFactor();
542  TimingCorrelation->Fill(correlationFactor);
543 
544  //--------------------------------------------------------------------------
545 
546  double diff_true_recon_vertex_x = reconVertexX - trueVertexX;
547  double diff_true_recon_vertex_y = reconVertexY - trueVertexY;
548  double diff_true_recon_vertex_z = reconVertexZ - trueVertexZ;
549  double distFromReconVertToTrueVert =
550  TMath::Sqrt(TMath::Power(diff_true_recon_vertex_x, 2) +
551  TMath::Power(diff_true_recon_vertex_y, 2) +
552  TMath::Power(diff_true_recon_vertex_z, 2));
553  double thetaOrientationFromReconVertToTrueVert =
554  TMath::ATan(diff_true_recon_vertex_y / diff_true_recon_vertex_x);
555  double phiOrientationFromReconVertToTrueVert =
556  TMath::ACos(diff_true_recon_vertex_z / distFromReconVertToTrueVert);
557  DistanceFromReconVertexToTrueVertex->Fill(distFromReconVertToTrueVert);
559  phiOrientationFromReconVertToTrueVert,
560  thetaOrientationFromReconVertToTrueVert);
562  phiOrientationFromReconVertToTrueVert);
564  thetaOrientationFromReconVertToTrueVert);
565 
566  // // Consider timing correlation
567  // // True Vertex
568  // const unsigned int mapSize =
569  // distanceToVertex_relTimingToVertex_Map.size(); unsigned int countElem =
570  // 0; double *distanceToVertex_array = new double[10]; double
571  // *relTimingToVertex_array = new double[10]; for (auto &elem :
572  // distanceToVertex_relTimingToVertex_Map) {
573  // if (4 < countElem && countElem < mapSize - 5)
574  // continue;
575  // *(distanceToVertex_array + countElem) = elem.first;
576  // *(relTimingToVertex_array + countElem) = elem.second;
577  // countElem++;
578  // }
579  // TGraph *timingCorrelation =
580  // new TGraph(10, distanceToVertex_array, relTimingToVertex_array);
581  // double correlationFactor = timingCorrelation->GetCorrelationFactor();
582  // TimingCorrelation->Fill(correlationFactor);
583 
584  // CONSIDERING THE ENERGY DISTRIBUTION TOPOLOGY
585  EnergyRatioFromVertexX->Fill(totalEnergyOfHitsAboveVertexX / totalE);
586  EnergyRatioFromVertexY->Fill(totalEnergyOfHitsAboveVertexY / totalE);
587  EnergyRatioFromVertexZ->Fill(totalEnergyOfHitsAboveVertexZ / totalE);
588 
589  // Filling analysis tree
590  _tRunNumber = e.run();
591  _tSubrunNumber = e.subRun();
592  _tEventNumber = e.id().event();
593  _tSliceNumber = std::distance(clusters.begin(), clusterItr);
595  double totVisEX =
597  double totVisEY =
599  _tTotalVisibleEnergy = totVisEX + totVisEY;
600  if (_tTotalVisibleEnergy > 0.8) {
602  TMath::Abs(diff_true_recon_vertex_z));
604  TMath::Abs(diff_true_elasticArm_vertex_z));
605  }
606 
615  _tHitVsVertexTimingCorrelation = correlationFactor;
616  _tReconVertexX = reconVertexX;
617  _tReconVertexY = reconVertexY;
618  _tReconVertexZ = reconVertexZ;
619  _tReconVertexT = reconVertexT;
620  if (_tSampleType != 0) {
621  _tTrueVertexX = trueVertexX;
622  _tTrueVertexY = trueVertexY;
623  _tTrueVertexZ = trueVertexZ;
624  _tTrueVertexT = trueVertexT;
625  } else {
626  _tTrueVertexX = -1.0E9;
627  _tTrueVertexY = -1.0E9;
628  _tTrueVertexZ = -1.0E9;
629  _tTrueVertexT = -1.0E9;
630  }
631  _VariableTree->Fill();
632 
633  clusterIdx++;
634  }
635 
636  HitsXView.clear();
637  HitsYView.clear();
638 }
639 
640 void nnbar::RecoAnalysis::evd(art::Event const &e, unsigned int clusterIdx) {
641 
642  // TGRAPH XVIEW
643  const unsigned int nhit_xview = HitsXView.size();
644  double *x_xv = new double[nhit_xview];
645  double *y_xv = new double[nhit_xview];
646  unsigned int countHitXView = 0;
647  for (auto &hit : HitsXView) {
648  *(x_xv + countHitXView) = std::get<0>(hit);
649  *(y_xv + countHitXView) = std::get<1>(hit);
650  countHitXView++;
651  }
652  TGraph *graph_xv = new TGraph(nhit_xview, x_xv, y_xv);
653 
654  // TGRAPH YVIEW
655  const unsigned int nhit_yview = HitsYView.size();
656  double *x_yv = new double[nhit_yview];
657  double *y_yv = new double[nhit_yview];
658  unsigned int countHitYView = 0;
659  for (auto &hit : HitsYView) {
660  *(x_yv + countHitYView) = std::get<0>(hit);
661  *(y_yv + countHitYView) = std::get<1>(hit);
662  countHitYView++;
663  }
664  TGraph *graph_yv = new TGraph(nhit_yview, x_yv, y_yv);
665 
666  // DRAW
667  TCanvas *c =
668  new TCanvas(Form("c_%i_%i", e.event(), clusterIdx), "", 1200, 1200);
669  c->Divide(1, 2);
670  c->cd(1);
671  graph_xv->SetTitle("");
672  graph_xv->GetXaxis()->SetLimits(0, 6400);
673  graph_xv->GetYaxis()->SetRangeUser(-800, 800);
674  graph_xv->GetXaxis()->SetTitle("z");
675  graph_xv->GetYaxis()->SetTitle("x");
676  graph_xv->GetXaxis()->CenterTitle();
677  graph_xv->GetYaxis()->CenterTitle();
678  graph_xv->SetMarkerColor(kRed);
679  graph_xv->SetMarkerStyle(20);
680  graph_xv->SetMarkerSize(1);
681  graph_xv->Draw("AP");
682  c->cd(2);
683  graph_yv->SetTitle("");
684  graph_yv->GetXaxis()->SetLimits(0, 6400);
685  graph_yv->GetYaxis()->SetRangeUser(-800, 800);
686  graph_yv->GetXaxis()->SetTitle("z");
687  graph_yv->GetYaxis()->SetTitle("y");
688  graph_yv->GetXaxis()->CenterTitle();
689  graph_yv->GetYaxis()->CenterTitle();
690  graph_yv->SetMarkerColor(kRed);
691  graph_yv->SetMarkerStyle(20);
692  graph_yv->SetMarkerSize(1);
693  graph_yv->Draw("AP");
694  c->SaveAs(Form("Plots/EVD_%i_%i.pdf", e.event(), clusterIdx));
695 }
696 
698  std::cout << "Exposure Time In USec: " << _ExposureTimeInUSec << std::endl;
699  // gStyle->SetOptStat(0);
700 
701  //----------------------------------------------------------------------------
702 
703  TCanvas *canvas = new TCanvas("canvas", "Canvas", 2400, 800);
704  canvas->Divide(3, 1);
705  canvas->cd(1);
706  EnergyRatioFromVertexX->SetTitle("x");
707  EnergyRatioFromVertexX->Draw();
708  canvas->cd(2);
709  EnergyRatioFromVertexY->SetTitle("y");
710  EnergyRatioFromVertexY->Draw();
711  canvas->cd(3);
712  EnergyRatioFromVertexZ->SetTitle("z");
713  EnergyRatioFromVertexZ->Draw();
714  canvas->SaveAs(Form("EnergyVsTopology.pdf"));
715 
716  TCanvas *canvas_timing =
717  new TCanvas("canvas_timing", "Canvas Timing", 1200, 1200);
718  TimingCorrelation->GetXaxis()->SetTitle("Timing Correlation Factor");
719  TimingCorrelation->GetXaxis()->CenterTitle();
720  TimingCorrelation->GetYaxis()->SetTitle("");
721  TimingCorrelation->GetYaxis()->CenterTitle();
722  TimingCorrelation->Draw("HIST");
723  canvas_timing->SaveAs(Form("TimingCorrelation.pdf"));
724 
725  TCanvas *canvas_vertex_error =
726  new TCanvas("canvas_vertex_error", "Canvas Vertex Error", 2400, 800);
727  canvas_vertex_error->Divide(3, 1);
728  canvas_vertex_error->cd(1);
729  DistanceFromReconVertexToTrueVertex->GetXaxis()->SetTitle(
730  "#deltar_{equale-(v)} (cm)");
731  DistanceFromReconVertexToTrueVertex->GetXaxis()->CenterTitle();
732  DistanceFromReconVertexToTrueVertex->GetYaxis()->SetTitle("");
733  DistanceFromReconVertexToTrueVertex->GetYaxis()->CenterTitle();
735  canvas_vertex_error->cd(2);
736  PolarAngleFromReconVertexToTrueVertex->GetXaxis()->SetTitle(
737  "#phi_{equale-(v)} (rad)");
738  PolarAngleFromReconVertexToTrueVertex->GetXaxis()->CenterTitle();
739  PolarAngleFromReconVertexToTrueVertex->GetYaxis()->SetTitle("");
740  PolarAngleFromReconVertexToTrueVertex->GetYaxis()->CenterTitle();
742  canvas_vertex_error->cd(3);
743  AzimuthalAngleFromReconVertexToTrueVertex->GetXaxis()->SetTitle(
744  "#theta_{equale-(v)} (rad)");
745  AzimuthalAngleFromReconVertexToTrueVertex->GetXaxis()->CenterTitle();
746  AzimuthalAngleFromReconVertexToTrueVertex->GetYaxis()->SetTitle("");
747  AzimuthalAngleFromReconVertexToTrueVertex->GetYaxis()->CenterTitle();
749  canvas_vertex_error->SaveAs(Form("VertexError.pdf"));
750 
751  TCanvas *canvas_earm_vertex_error = new TCanvas(
752  "canvas_earm_vertex_error", "Canvas EArm Vertex Error", 2400, 800);
753  canvas_earm_vertex_error->Divide(3, 1);
754  canvas_earm_vertex_error->cd(1);
755  DistanceFromElasticArmVertexToTrueVertex->GetXaxis()->SetTitle(
756  "#deltar_{earm-(v)} (cm)");
757  DistanceFromElasticArmVertexToTrueVertex->GetXaxis()->CenterTitle();
758  DistanceFromElasticArmVertexToTrueVertex->GetYaxis()->SetTitle("");
759  DistanceFromElasticArmVertexToTrueVertex->GetYaxis()->CenterTitle();
761  canvas_earm_vertex_error->cd(2);
762  PolarAngleFromElasticArmVertexToTrueVertex->GetXaxis()->SetTitle(
763  "#phi_{earm-(v)} (rad)");
764  PolarAngleFromElasticArmVertexToTrueVertex->GetXaxis()->CenterTitle();
765  PolarAngleFromElasticArmVertexToTrueVertex->GetYaxis()->SetTitle("");
766  PolarAngleFromElasticArmVertexToTrueVertex->GetYaxis()->CenterTitle();
768  canvas_earm_vertex_error->cd(3);
770  "#theta_{earm-(v)} (rad)");
771  AzimuthalAngleFromElasticArmVertexToTrueVertex->GetXaxis()->CenterTitle();
772  AzimuthalAngleFromElasticArmVertexToTrueVertex->GetYaxis()->SetTitle("");
773  AzimuthalAngleFromElasticArmVertexToTrueVertex->GetYaxis()->CenterTitle();
775  canvas_earm_vertex_error->SaveAs(Form("EArmVertexError.pdf"));
776 
777  TCanvas *canvas_zvertex_error =
778  new TCanvas("canvas_zvertex_error", "Canvas Z Vertex Error", 2400, 1200);
779  canvas_zvertex_error->Divide(2, 1);
780  canvas_zvertex_error->cd(1);
781  ZErrorFromElasticArmVertexToTrueVertex->GetXaxis()->SetTitle(
782  "#deltaz_{earm-(v)} (cm)");
783  ZErrorFromElasticArmVertexToTrueVertex->GetXaxis()->CenterTitle();
784  ZErrorFromElasticArmVertexToTrueVertex->GetYaxis()->SetTitle("");
785  ZErrorFromElasticArmVertexToTrueVertex->GetYaxis()->CenterTitle();
787  canvas_zvertex_error->cd(2);
788  ZErrorFromReconVertexToTrueVertex->GetXaxis()->SetTitle(
789  "#deltaz_{equale-(v)} (cm)");
790  ZErrorFromReconVertexToTrueVertex->GetXaxis()->CenterTitle();
791  ZErrorFromReconVertexToTrueVertex->GetYaxis()->SetTitle("");
792  ZErrorFromReconVertexToTrueVertex->GetYaxis()->CenterTitle();
793  ZErrorFromReconVertexToTrueVertex->Draw("HIST");
794  canvas_zvertex_error->SaveAs(Form("ZVertexError.pdf"));
795 }
796 
RecoAnalysis(fhicl::ParameterSet const &p)
float T() const
Definition: RecoHit.h:39
TH1D * DistanceFromElasticArmVertexToTrueVertex
SubRunNumber_t subRun() const
Definition: Event.h:72
double HozVertex(std::vector< std::tuple< double, double, double >> HitCollectionSingleView, double vertexZ)
TH1D * AzimuthalAngleFromElasticArmVertexToTrueVertex
std::vector< std::tuple< double, double, double > > HitsXView
iterator begin()
Definition: PtrVector.h:223
Definition: geo.h:1
TH1D * PolarAngleFromReconVertexToTrueVertex
const char * p
Definition: xmltok.h:285
Vertical planes which measure X.
Definition: PlaneGeo.h:28
float Z() const
Definition: RecoHit.h:38
std::vector< std::tuple< double, double, double > > HitsYView
TH2D * OrientationFromElasticArmVertexToTrueVertex
DEFINE_ART_MODULE(TestTMapFile)
unsigned distance(const T &t1, const T &t2)
Particle class.
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
const Var kY([](const caf::SRProxy *sr){float tmp=0.f;if(sr->mc.nu.empty()) return tmp;tmp=sr->mc.nu[0].y;return tmp;})
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
TH1D * ZErrorFromElasticArmVertexToTrueVertex
static double hitCountInView(rb::Cluster cluster, geo::View_t aView)
T get(std::string const &key) const
Definition: ParameterSet.h:231
iterator end()
Definition: PtrVector.h:237
void evd(art::Event const &e, unsigned int clusterIdx)
TH1D * AzimuthalAngleFromReconVertexToTrueVertex
void analyze(art::Event const &e) override
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
Definition: RecoHit.cxx:35
const double j
Definition: BetheBloch.cxx:29
EventNumber_t event() const
Definition: Event.h:67
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
EventNumber_t event() const
Definition: EventID.h:116
double slopeAngleOfLineConnect2Points(point2d_t pA, point2d_t pB)
T * make(ARGS...args) const
float GeV() const
Definition: RecoHit.cxx:69
std::pair< double, double > point2d_t
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
Definition: structs.h:12
float X() const
Definition: RecoHit.h:36
TH2D * OrientationFromReconVertexToTrueVertex
static double totalVisibleEnergyInView(rb::Cluster cluster, geo::View_t aView)
float Y() const
Definition: RecoHit.h:37
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:464
Float_t e
Definition: plot.C:35
RunNumber_t run() const
Definition: Event.h:77
TH1D * PolarAngleFromElasticArmVertexToTrueVertex
EventID id() const
Definition: Event.h:56
static double temporalExpand(rb::Cluster cluster)
TGraph * PlotVectorPoint(std::vector< point2d_t > pts)