extract_resolution.C
Go to the documentation of this file.
1 // Extract the vertex resolution from ElasticArms and HoughVertex
2 // and save to "outname" the 1d and 2d resolution and reco-true
3 // distributions for comparing the two trackers.
4 // Run as: cafe -bq extract_resoltion.C <output_filename>
5 //
6 // Authors: C. Backhouse and G. Davies
7 /////////////////////////////////////////////////////////////////
8 
9 #include "CAFAna/Core/Binning.h"
10 #include "CAFAna/Core/Spectrum.h"
13 #include "CAFAna/Vars/Vars.h"
14 
15 using namespace ana;
16 
18 
19 #include "TFile.h"
20 
21 #include <cmath>
22 
23 const Cut kVertices([](const caf::SRProxy* sr)
24  {
25  return sr->mc.nnu > 0 && sr->vtx.nhough > 0 && sr->vtx.elasitc.IsValid;
26  });
27 
28 const Cut kMuCatch([](const caf::SRProxy* sr)
29  {
30  return sr->hdr.det == 1 && sr->mc.nnu > 0 && sr->mc.nu[0].vtx.Z() > 1270;
31  });
32 
33 class DiffHough
34 {
35 public:
36  DiffHough(int a) : axis(a) {}
37  double operator()(const caf::SRProxy* sr){
38  return sr->vtx.hough[0].vtx(axis) - sr->mc.nu[0].vtx(axis);
39  }
40  int axis;
41 };
42 
44 {
45 public:
46  DiffElastic(int a) : axis(a) {}
47  double operator()(const caf::SRProxy* sr){
48  return sr->vtx.elastic.vtx(axis) - sr->mc.nu[0].vtx(axis);
49  }
50  int axis;
51 };
52 
53 const Var kXDiffHough(DiffHough(0));
54 const Var kYDiffHough(DiffHough(1));
55 const Var kZDiffHough(DiffHough(2));
56 
57 const Var kXDiffElastic(DiffElastic(0));
58 const Var kYDiffElastic(DiffElastic(1));
59 const Var kZDiffElastic(DiffElastic(2));
60 
61 const Var kDistHough([](const caf::SRProxy* sr)
62  {
63  return (sr->mc.nu[0].vtx-sr->vtx.hough[0].vtx).Mag();
64  });
65 
66 const Var kVolumeNormHough([](const caf::SRProxy* sr)
67  {
68  const double r = (sr->mc.nu[0].vtx-sr->vtx.hough[0].vtx).Mag();
69  return 1/(4.*M_PI*r*r);
70  });
71 
72 const Var kDistElastic([](const caf::SRProxy* sr)
73  {
74  return (sr->mc.nu[0].vtx-sr->vtx.elastic.vtx).Mag();
75  });
76 
77 const Var kVolumeNormElastic([](const caf::SRProxy* sr)
78  {
79  const double r = (sr->mc.nu[0].vtx-sr->vtx.elastic.vtx).Mag();
80  return 1/(4.*M_PI*r*r);
81  });
82 
83 // const Var kDistDiff = kDistHough - kDistElastic;
84 const Var kDistDiff( [](const caf::SRProxy* sr)
85  {
86  return (sr->mc.nu[0].vtx-sr->vtx.hough[0].vtx).Mag() - (sr->mc.nu[0].vtx-sr->vtx.elastic.vtx).Mag();
87  });
88 
90 {
91  SpectrumLoader loader("prod_decaf_R16-03-03-prod2reco.f_fd_genie_nonswap_fhc_nova_v08_period1_nue_or_numu_contain_v1_prod2-snapshot");
92 
93  const Binning diffBins = Binning::Simple(500, -100, +100);
94  Spectrum xdh("Reco - True Vertex X (cm)", diffBins, loader, kXDiffHough, kVertices && !kMuCatch);
95  Spectrum ydh("Reco - True Vertex Y (cm)", diffBins, loader, kYDiffHough, kVertices && !kMuCatch);
96  Spectrum zdh("Reco - True Vertex Z (cm)", diffBins, loader, kZDiffHough, kVertices && !kMuCatch);
97 
98  Spectrum xde("Reco - True Vertex X (cm)", diffBins, loader, kXDiffElastic, kVertices && !kMuCatch);
99  Spectrum yde("Reco - True Vertex Y (cm)", diffBins, loader, kYDiffElastic, kVertices && !kMuCatch);
100  Spectrum zde("Reco - True Vertex Z (cm)", diffBins, loader, kZDiffElastic, kVertices && !kMuCatch);
101 
102 
103  const Binning diff2DBins = Binning::Simple(100, -100, +100);
104  Spectrum xd2d(";HoughVertex Reco - True X (cm);ElasticArms Reco - True X (cm)", loader, diff2DBins, kXDiffHough, diff2DBins, kXDiffElastic, kVertices && !kMuCatch);
105  Spectrum yd2d(";HoughVertex Reco - True Y (cm);ElasticArms Reco - True Y (cm)", loader, diff2DBins, kYDiffHough, diff2DBins, kXDiffElastic, kVertices && !kMuCatch);
106  Spectrum zd2d(";HoughVertex Reco - True Z (cm);ElasticArms Reco - True Z (cm)", loader, diff2DBins, kZDiffHough, diff2DBins, kZDiffElastic, kVertices && !kMuCatch);
107 
108 
109  const Binning distBins = Binning::Simple(500, 0, +100);
110  const Binning hitBins = Binning::Simple(250, 0, +250);
111  const Binning calBins = Binning::Simple(50, 0, +5);
112  Spectrum disth("Distance from true vertex (cm)", distBins, loader, kDistHough, kVertices && !kMuCatch);
113  Spectrum diste("Distance from true vertex (cm)", distBins, loader, kDistElastic, kVertices && !kMuCatch);
114 
115  Spectrum disth2d(";HV: Distance from true vertex (cm); No. of hits", loader, distBins, kDistHough, hitBins, kNHit, kVertices && !kMuCatch);
116  Spectrum diste2d(";EA: Distance from true vertex (cm); No. of hits", loader, distBins, kDistElastic, hitBins, kNHit, kVertices && !kMuCatch);
117  Spectrum calh2d(";HV: Distance from true vertex (cm); Calorimetric energy (GeV)", loader, distBins, kDistHough, calBins, kCaloE, kVertices && !kMuCatch);
118  Spectrum cale2d(";EA: Distance from true vertex (cm); Calorimetric energy (GeV)", loader, distBins, kDistElastic, calBins, kCaloE, kVertices && !kMuCatch);
119 
120  const Binning distNormBins = Binning::Simple(2500, 0, +100);
121  Spectrum distnormh("Distance from true vertex (cm)", distNormBins, loader, kDistHough, kVertices && !kMuCatch, kNoShift, kVolumeNormHough);
122  Spectrum distnorme("Distance from true vertex (cm)", distNormBins, loader, kDistElastic, kVertices && !kMuCatch, kNoShift, kVolumeNormElastic);
123 
124  Spectrum distdiff("HoughVertex - ElasticArms distance from true vertex (cm)", diffBins, loader, kDistDiff, kVertices && !kMuCatch);
125 
126 
127  loader.Go();
128 
129 
130  TFile* fout = new TFile(outname.c_str(), "RECREATE");
131 
132  SaveTo(xdh, fout, "xdh");
133  SaveTo(ydh, fout, "ydh");
134  SaveTo(zdh, fout, "zdh");
135  SaveTo(xde, fout, "xde");
136  SaveTo(yde, fout, "yde");
137  SaveTo(zde, fout, "zde");
138  SaveTo(xd2d, fout, "xd2d");
139  SaveTo(yd2d, fout, "yd2d");
140  SaveTo(zd2d, fout, "zd2d");
141  SaveTo(disth, fout, "disth");
142  SaveTo(diste, fout, "diste");
143  SaveTo(disth2d, fout, "disth2d");
144  SaveTo(diste2d, fout, "diste2d");
145  SaveTo(calh2d, fout, "calh2d");
146  SaveTo(cale2d, fout, "cale2d");
147  SaveTo(distnormh, fout, "distnormh");
148  SaveTo(distnorme, fout, "distnorme");
149  SaveTo(distdiff, fout, "distdiff");
150 }
const Var kYDiffHough(DiffHough(1))
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const Var kXDiffHough(DiffHough(0))
void SaveTo(const osc::IOscCalc &x, TDirectory *dir, const std::string &name)
caf::Proxy< caf::SRHeader > hdr
Definition: SRProxy.h:2136
const Var kYDiffElastic(DiffElastic(1))
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2125
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:617
const Cut kVertices([](const caf::SRProxy *sr){return sr->mc.nnu > 0 &&sr->vtx.nhough > 0 &&sr->vtx.elasitc.IsValid;})
const Var kZDiffHough(DiffHough(2))
const Var kXDiffElastic(DiffElastic(0))
caf::Proxy< short int > nnu
Definition: SRProxy.h:616
const Cut kMuCatch([](const caf::SRProxy *sr){return sr->hdr.det==1 &&sr->mc.nnu > 0 &&sr->mc.nu[0].vtx.Z() > 1270;})
#define M_PI
Definition: SbMath.h:34
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
caf::Proxy< caf::SRElastic > elastic
Definition: SRProxy.h:2117
caf::Proxy< std::vector< caf::SRHoughVertex > > hough
Definition: SRProxy.h:2118
caf::Proxy< size_t > nhough
Definition: SRProxy.h:2119
const Var kCaloE
Summed calorimetric energy of all hits in slice, uncorrected.
Definition: Vars.cxx:52
void extract_resolution(std::string outname)
const double a
const Var kNHit
Definition: Vars.cxx:71
caf::StandardRecord * sr
virtual void Go() override
Load all the registered spectra.
loader
Definition: demo0.py:10
const Var kZDiffElastic(DiffElastic(2))
const SystShifts kNoShift
Definition: SystShifts.cxx:21
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2137
const Var kVolumeNormHough([](const caf::SRProxy *sr){const double r=(sr->mc.nu[0].vtx-sr->vtx.hough[0].vtx).Mag();return 1/(4.*M_PI *r *r);})
const Var kDistDiff([](const caf::SRProxy *sr){return(sr->mc.nu[0].vtx-sr->vtx.hough[0].vtx).Mag()-(sr->mc.nu[0].vtx-sr->vtx.elastic.vtx).Mag();})
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const Var kDistHough([](const caf::SRProxy *sr){return(sr->mc.nu[0].vtx-sr->vtx.hough[0].vtx).Mag();})
float Mag() const
caf::Proxy< caf::SRVector3D > vtx
Definition: SRProxy.h:2072
TRandom3 r(0)
const Var kVolumeNormElastic([](const caf::SRProxy *sr){const double r=(sr->mc.nu[0].vtx-sr->vtx.elastic.vtx).Mag();return 1/(4.*M_PI *r *r);})
caf::Proxy< caf::SRVertexBranch > vtx
Definition: SRProxy.h:2145
double operator()(const caf::SRProxy *sr)
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
const Var kDistElastic([](const caf::SRProxy *sr){return(sr->mc.nu[0].vtx-sr->vtx.elastic.vtx).Mag();})
double operator()(const caf::SRProxy *sr)
caf::Proxy< caf::Det_t > det
Definition: SRProxy.h:231
enum BeamMode string