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();
92 Visibility(
int vx,
int t,
float d,
float v,
float p,
float th,
float xp)
93 : vox(vx), taken(t), dist(d), vis(v), psi(p), theta(th), xpos(xp)
103 TCanvas*
c1 =
new TCanvas(
"c1",
"c1");
105 c1->SetWindowSize(600, 600);
110 std::vector<Visibility> viss;
113 for (
unsigned int voxIdx = 0U; voxIdx < voxdef.
GetNVoxels(); ++voxIdx) {
116 const double xyzvox[] = {voxvec.X(), voxvec.Y(), voxvec.Z()};
119 const double fc_z = 1394;
120 const double fc_x = 350;
131 if ((voxvec.X() - opdetCenter.X()) < 0) {
132 psi = atan((voxvec.Y() - opdetCenter.Y()) / (-voxvec.X() + opdetCenter.X()));
135 if (voxvec.X() < fc_x && voxvec.X() > -fc_x && voxvec.Y() < fc_y && voxvec.Y() > -fc_y &&
136 voxvec.Z() > -9 && voxvec.Z() < fc_z) {
137 g.SetPoint(g.GetN(),
dist, vis * dist *
dist);
139 psi = atan((voxvec.Y() - opdetCenter.Y()) / (voxvec.X() - opdetCenter.X())) * 180.0 /
141 theta = acos((voxvec.Z() - opdetCenter.Z()) / dist) * 180.0 /
144 if ((voxvec.X() - opdetCenter.X()) < 0) {
145 psi = atan((voxvec.Y() - opdetCenter.Y()) / (-voxvec.X() + opdetCenter.X()));
149 viss.emplace_back(voxIdx, taken, dist, vis, psi, theta, xpos);
159 g.GetXaxis()->SetTitle(
"Distance (cm)");
160 g.GetYaxis()->SetTitle(
"Visibility #times r^{2}");
162 TF1* fit = g.GetFunction(
"expo");
166 fitres[0] = fit->GetParameter(0);
167 fitres[1] = fit->GetParameter(1);
172 TH1F h(
"",
"", 200, 0, 20);
175 TTree* tr_fit =
new TTree(
"tr",
"tr");
176 tr_fit->Branch(
"vox", &vox);
177 tr_fit->Branch(
"dist", &dist);
178 tr_fit->Branch(
"vis", &vis);
179 tr_fit->Branch(
"taken", &taken);
180 tr_fit->Branch(
"psi", &psi);
181 tr_fit->Branch(
"theta", &theta);
182 tr_fit->Branch(
"xpos", &xpos);
184 for (
const Visibility& v : viss) {
188 const double obs = v.vis / 2
e-5;
190 fit->Eval(v.dist) / (v.dist * v.dist) / 2
e-5;
199 double chisq = 2 * (pred - obs);
200 if (obs) chisq += 2 * obs * log(obs / pred);
210 if (v.taken == 1) { h.Fill(chisq); }
215 g2.SetPoint(g2.GetN(), v.dist, v.vis * v.dist * v.dist);
227 g2.SetMarkerStyle(7);
228 g2.SetMarkerColor(kRed);
229 gStyle->SetLabelSize(.04,
"XY");
230 gStyle->SetTitleSize(.04,
"XY");
233 c1->SaveAs(TString::Format(
"plots/Chris_vis_vs_dist_%d.png", opdetIdx).Data());
235 h.GetXaxis()->SetTitle(
"#chi^{2}");
236 h.GetYaxis()->SetTitle(
"Counts");
237 gStyle->SetLabelSize(.04,
"XY");
238 gStyle->SetTitleSize(.04,
"XY");
240 std::cout <<
"Integral(90,-1)/Integral(all) = " << h.Integral(90, -1) / h.Integral(0, -1)
242 std::cout <<
"*****************************" << std::endl;
243 totExceptions += h.Integral(90, -1);
244 totPts += h.Integral(0, -1);
245 gPad->Print(TString::Format(
"plots/chisq_opdet_%d.eps", opdetIdx).Data());
250 std::cout << totExceptions <<
" exceptions from " << totPts
251 <<
" points = " << (100. * totExceptions) / totPts <<
"%" << std::endl;
OpDetGeo const & OpDetGeoFromOpDet(unsigned int OpDet) const
Returns the geo::OpDetGeo object for the given detector number.
Representation of a region of space diced into voxels.
EDAnalyzer(fhicl::ParameterSet const &pset)
unsigned int NOpDets() const
Number of OpDets in the whole detector.
geo::Point_t const & GetCenter() const
const sim::PhotonVoxelDef & GetVoxelDef() const
double DistanceToPoint(geo::Point_t const &point) const
Returns the distance of the specified point from detector center [cm].
constexpr 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
PhotonVoxel GetPhotonVoxel(int ID) const