55 TFile* fout_full =
new TFile(
"full.root",
"RECREATE");
56 TFile* fout_fit =
new TFile(
"fit.root",
"RECREATE");
58 std::cout << voxdef.
GetNVoxels() <<
" voxels for each of " << geom->
NOpDets() <<
" OpDets" 60 std::cout << std::endl;
63 long totExceptions = 0;
66 for (
unsigned int opdetIdx = 0; opdetIdx < geom->
NOpDets(); ++opdetIdx) {
67 std::cout << opdetIdx <<
" / " << geom->
NOpDets() << std::endl;
69 TDirectory* d_full = fout_full->mkdir(TString::Format(
"opdet_%d", opdetIdx).Data());
70 TDirectory* d_fit = fout_fit->mkdir(TString::Format(
"opdet_%d", opdetIdx).Data());
73 TTree* tr_full =
new TTree(
"tr",
"tr");
75 float dist, vis, psi, theta, xpos;
76 tr_full->Branch(
"vox", &vox);
77 tr_full->Branch(
"dist", &dist);
78 tr_full->Branch(
"vis", &vis);
79 tr_full->Branch(
"taken", &taken);
80 tr_full->Branch(
"psi",
82 tr_full->Branch(
"theta",
84 tr_full->Branch(
"xpos",
88 auto const opdetCenter = opdet.
GetCenter();
91 Visibility(
int vx,
int t,
float d,
float v,
float p,
float th,
float xp)
92 :
vox(vx), taken(t),
dist(d), vis(v), psi(p), theta(th), xpos(xp)
102 TCanvas*
c1 =
new TCanvas(
"c1",
"c1");
103 c1->SetWindowSize(600, 600);
107 std::vector<Visibility> viss;
110 for (
unsigned int voxIdx = 0U; voxIdx < voxdef.
GetNVoxels(); ++voxIdx) {
113 const double xyzvox[] = {voxvec.X(), voxvec.Y(), voxvec.Z()};
116 const double fc_z = 1394;
117 const double fc_x = 350;
128 if ((voxvec.X() - opdetCenter.X()) < 0) {
129 psi = atan((voxvec.Y() - opdetCenter.Y()) / (-voxvec.X() + opdetCenter.X()));
132 if (voxvec.X() < fc_x && voxvec.X() > -fc_x && voxvec.Y() < fc_y && voxvec.Y() > -fc_y &&
133 voxvec.Z() > -9 && voxvec.Z() < fc_z) {
134 g.SetPoint(g.GetN(),
dist, vis * dist *
dist);
136 psi = atan((voxvec.Y() - opdetCenter.Y()) / (voxvec.X() - opdetCenter.X())) * 180.0 /
138 theta = acos((voxvec.Z() - opdetCenter.Z()) / dist) * 180.0 /
141 if ((voxvec.X() - opdetCenter.X()) < 0) {
142 psi = atan((voxvec.Y() - opdetCenter.Y()) / (-voxvec.X() + opdetCenter.X()));
146 viss.emplace_back(voxIdx, taken, dist, vis, psi, theta, xpos);
156 g.GetXaxis()->SetTitle(
"Distance (cm)");
157 g.GetYaxis()->SetTitle(
"Visibility #times r^{2}");
159 TF1* fit = g.GetFunction(
"expo");
163 fitres[0] = fit->GetParameter(0);
164 fitres[1] = fit->GetParameter(1);
169 TH1F h(
"",
"", 200, 0, 20);
172 TTree* tr_fit =
new TTree(
"tr",
"tr");
173 tr_fit->Branch(
"vox", &vox);
174 tr_fit->Branch(
"dist", &dist);
175 tr_fit->Branch(
"vis", &vis);
176 tr_fit->Branch(
"taken", &taken);
177 tr_fit->Branch(
"psi", &psi);
178 tr_fit->Branch(
"theta", &theta);
179 tr_fit->Branch(
"xpos", &xpos);
181 for (
const Visibility& v : viss) {
185 const double obs = v.vis / 2
e-5;
187 fit->Eval(v.dist) / (v.dist * v.dist) / 2
e-5;
193 double chisq = 2 * (pred - obs);
194 if (obs) chisq += 2 * obs * log(obs / pred);
203 if (v.taken == 1) { h.Fill(chisq); }
208 g2.SetPoint(g2.GetN(), v.dist, v.vis * v.dist * v.dist);
219 g2.SetMarkerStyle(7);
220 g2.SetMarkerColor(kRed);
221 gStyle->SetLabelSize(.04,
"XY");
222 gStyle->SetTitleSize(.04,
"XY");
225 c1->SaveAs(TString::Format(
"plots/Chris_vis_vs_dist_%d.png", opdetIdx).Data());
227 h.GetXaxis()->SetTitle(
"#chi^{2}");
228 h.GetYaxis()->SetTitle(
"Counts");
229 gStyle->SetLabelSize(.04,
"XY");
230 gStyle->SetTitleSize(.04,
"XY");
232 std::cout <<
"Integral(90,-1)/Integral(all) = " << h.Integral(90, -1) / h.Integral(0, -1)
234 std::cout <<
"*****************************" << std::endl;
235 totExceptions += h.Integral(90, -1);
236 totPts += h.Integral(0, -1);
237 gPad->Print(TString::Format(
"plots/chisq_opdet_%d.eps", opdetIdx).Data());
242 std::cout << totExceptions <<
" exceptions from " << totPts
243 <<
" points = " << (100. * totExceptions) / totPts <<
"%" << std::endl;
Definitions of voxel data structures.
Point_t const & GetCenter() const
Representation of a region of space diced into voxels.
void analyze(const art::Event &e) override
double DistanceToPoint(Point_t const &point) const
Returns the distance of the specified point from detector center [cm].
#define DEFINE_ART_MODULE(klass)
unsigned int NOpDets() const
Number of OpDets in the whole detector.
Encapsulate the geometry of an optical detector.
const sim::PhotonVoxelDef & GetVoxelDef() const
General LArSoft Utilities.
double dist(const TReal *x, const TReal *y, const unsigned int dimension)
Point GetCenter() const
Returns the center of the voxel (type Point).
unsigned int GetNVoxels() const
Returns the total number of voxels in the volume.
float GetVisibility(Point const &p, unsigned int OpChannel, bool wantReflected=false) const
art framework interface to geometry description
OpDetGeo const & OpDetGeoFromOpDet(unsigned int OpDet) const
Returns the geo::OpDetGeo object for the given detector number.
PhotonVoxel GetPhotonVoxel(int ID) const