24 #include "Utilities/AssociationUtil.h" 118 std::vector<std::tuple<double, double, double>> HitCollectionSingleView,
202 tfs->
make<TTree>(
"VariableTree",
"NNbar Analysis Variable Tree");
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);
240 "TimingCorrelation",
";Timing Correlation Factor;", 100, -1.5, 1.5);
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);
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);
272 double *
x =
new double[pts.size()];
273 double *
y =
new double[pts.size()];
276 for (
auto &
pt : pts) {
278 *(y + i) =
pt.second;
282 TGraph *graph =
new TGraph(pts.size(),
x,
y);
288 double xa = pA.first;
289 double ya = pA.second;
290 double xb = pB.first;
291 double yb = pB.second;
293 return TMath::ATan((ya - yb) / (xa - xb));
301 std::vector<art::Ptr<simb::MCTruth>> mctruths;
305 double trueVertexX = 0, trueVertexY = 0, trueVertexZ = 0, trueVertexT = 0;
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();
318 std::vector<art::Ptr<rb::Cluster>> clusters;
325 double elasticArmVertX;
326 double elasticArmVertY;
327 double elasticArmVertZ;
328 unsigned int clusterIdx = 0;
331 clusterItr != clusters.end(); ++clusterItr) {
332 auto cluster = *clusterItr;
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)
346 elasticArmVertX = vert.at(0)->GetX();
347 elasticArmVertY = vert.at(0)->GetY();
348 elasticArmVertZ = vert.at(0)->GetZ();
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);
366 distFromElasticArmVertToTrueVert);
368 phiOrientationFromElasticArmVertToTrueVert,
369 thetaOrientationFromElasticArmVertToTrueVert);
371 phiOrientationFromElasticArmVertToTrueVert);
373 thetaOrientationFromElasticArmVertToTrueVert);
376 double totalE = (*cluster).TotalGeV();
379 double totalEnergyOfHitsAboveVertexX = 0.;
380 double totalEnergyOfHitsAboveVertexY = 0.;
381 double totalEnergyOfHitsAboveVertexZ = 0.;
384 std::multimap<double, double> distanceToVertex_relTimingToVertex_Map;
385 unsigned int hitIdx = 0;
388 std::multimap<double, double> HitsZ, HitsX, HitsY;
391 itrCell != allCells.
end(); itrCell++) {
400 double cX = recoHit.
X();
401 double cY = recoHit.
Y();
402 double cZ = recoHit.
Z();
404 double cE = recoHit.
GeV();
406 HitsZ.insert(std::pair<double, double>(cZ, cE));
408 HitsX.insert(std::pair<double, double>(cX, cE));
411 HitsY.insert(std::pair<double, double>(cY, cE));
427 if (cX > trueVertexX)
428 totalEnergyOfHitsAboveVertexX += cE;
430 if (cY > trueVertexY)
431 totalEnergyOfHitsAboveVertexY += cE;
433 if (cZ > trueVertexZ)
434 totalEnergyOfHitsAboveVertexZ += cE;
438 HitsXView.push_back(std::make_tuple(cZ, cX, cE));
440 HitsYView.push_back(std::make_tuple(cZ, cY, cE));
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;
458 double reconVertexZ =
459 (lowerBoundOfMedianZVertex + upperBoundOfMedianZVertex) / 2.;
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;
472 double reconVertexX =
473 (lowerBoundOfMedianXVertex + upperBoundOfMedianXVertex) / 2.;
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;
486 double reconVertexY =
487 (lowerBoundOfMedianYVertex + upperBoundOfMedianYVertex) / 2.;
490 itrCell != allCells.
end(); itrCell++) {
498 double cX = recoHit.
X();
499 double cY = recoHit.
Y();
500 double cZ = recoHit.
Z();
501 double cT = recoHit.
T();
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));
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)
523 *(distanceToVertex_array + countElem) = elem.first;
524 *(relTimingToVertex_array + countElem) = elem.second;
530 double reconVertexT = 0.;
531 for (
unsigned int j = 0;
j < 5;
j++) {
532 reconVertexT += *(relTimingToVertex_array +
j);
534 reconVertexT = reconVertexT / 5;
535 for (
unsigned int j = 0;
j < 5;
j++) {
536 *(relTimingToVertex_array + countElem) =
537 *(relTimingToVertex_array + countElem) - reconVertexT;
539 TGraph *timingCorrelation =
540 new TGraph(10, distanceToVertex_array, relTimingToVertex_array);
541 double correlationFactor = timingCorrelation->GetCorrelationFactor();
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);
559 phiOrientationFromReconVertToTrueVert,
560 thetaOrientationFromReconVertToTrueVert);
562 phiOrientationFromReconVertToTrueVert);
564 thetaOrientationFromReconVertToTrueVert);
602 TMath::Abs(diff_true_recon_vertex_z));
604 TMath::Abs(diff_true_elasticArm_vertex_z));
620 if (_tSampleType != 0) {
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;
648 *(x_xv + countHitXView) = std::get<0>(
hit);
649 *(y_xv + countHitXView) = std::get<1>(
hit);
652 TGraph *graph_xv =
new TGraph(nhit_xview, x_xv, y_xv);
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;
660 *(x_yv + countHitYView) = std::get<0>(
hit);
661 *(y_yv + countHitYView) = std::get<1>(
hit);
664 TGraph *graph_yv =
new TGraph(nhit_yview, x_yv, y_yv);
668 new TCanvas(Form(
"c_%i_%i", e.
event(), clusterIdx),
"", 1200, 1200);
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");
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));
703 TCanvas *
canvas =
new TCanvas(
"canvas",
"Canvas", 2400, 800);
704 canvas->Divide(3, 1);
714 canvas->SaveAs(Form(
"EnergyVsTopology.pdf"));
716 TCanvas *canvas_timing =
717 new TCanvas(
"canvas_timing",
"Canvas Timing", 1200, 1200);
723 canvas_timing->SaveAs(Form(
"TimingCorrelation.pdf"));
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);
730 "#deltar_{equale-(v)} (cm)");
735 canvas_vertex_error->cd(2);
737 "#phi_{equale-(v)} (rad)");
742 canvas_vertex_error->cd(3);
744 "#theta_{equale-(v)} (rad)");
749 canvas_vertex_error->SaveAs(Form(
"VertexError.pdf"));
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);
756 "#deltar_{earm-(v)} (cm)");
761 canvas_earm_vertex_error->cd(2);
763 "#phi_{earm-(v)} (rad)");
768 canvas_earm_vertex_error->cd(3);
770 "#theta_{earm-(v)} (rad)");
775 canvas_earm_vertex_error->SaveAs(Form(
"EArmVertexError.pdf"));
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);
782 "#deltaz_{earm-(v)} (cm)");
787 canvas_zvertex_error->cd(2);
789 "#deltaz_{equale-(v)} (cm)");
794 canvas_zvertex_error->SaveAs(Form(
"ZVertexError.pdf"));
RecoAnalysis(fhicl::ParameterSet const &p)
TH1D * ZErrorFromReconVertexToTrueVertex
TH1D * DistanceFromElasticArmVertexToTrueVertex
SubRunNumber_t subRun() const
double HozVertex(std::vector< std::tuple< double, double, double >> HitCollectionSingleView, double vertexZ)
TH1D * AzimuthalAngleFromElasticArmVertexToTrueVertex
double _tTemporalClusterExpand
std::vector< std::tuple< double, double, double > > HitsXView
TH1D * EnergyRatioFromVertexZ
TH1D * PolarAngleFromReconVertexToTrueVertex
Vertical planes which measure X.
std::vector< std::tuple< double, double, double > > HitsYView
double _ExposureTimeInUSec
TH2D * OrientationFromElasticArmVertexToTrueVertex
double _tHitCountXYDifference
TCanvas * canvas(const char *nm, const char *ti, const char *a)
unsigned int _tSubrunNumber
DEFINE_ART_MODULE(TestTMapFile)
unsigned distance(const T &t1, const T &t2)
Horizontal planes which measure Y.
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...
double _tTotalHitCountInYView
TH1D * ZErrorFromElasticArmVertexToTrueVertex
static double hitCountInView(rb::Cluster cluster, geo::View_t aView)
TH1D * DistanceFromReconVertexToTrueVertex
art::InputTag _InputShwrLabel
double _tTotalHitCountInXView
unsigned int _tEventNumber
T get(std::string const &key) const
void evd(art::Event const &e, unsigned int clusterIdx)
TH1D * AzimuthalAngleFromReconVertexToTrueVertex
void analyze(art::Event const &e) override
double _tAverageEnergyPerHitYView
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
TH1D * EnergyRatioFromVertexX
EventNumber_t event() const
EDAnalyzer(Table< Config > const &config)
Vertex location in position and time.
data_t::iterator iterator
double _tTotalVisibleEnergy
EventNumber_t event() const
double slopeAngleOfLineConnect2Points(point2d_t pA, point2d_t pB)
T * make(ARGS...args) const
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Var Sqrt(const Var &v)
Use to take sqrt of a var.
TH1D * EnergyRatioFromVertexY
double _tAverageEnergyPerHitXView
TH2D * OrientationFromReconVertexToTrueVertex
static double totalVisibleEnergyInView(rb::Cluster cluster, geo::View_t aView)
unsigned int _tSampleType
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
TH1D * PolarAngleFromElasticArmVertexToTrueVertex
static double temporalExpand(rb::Cluster cluster)
TGraph * PlotVectorPoint(std::vector< point2d_t > pts)
unsigned int _tSliceNumber
double _tHitVsVertexTimingCorrelation