LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
CreateHybridLibrary_module.cc
Go to the documentation of this file.
1 // Chris Backhouse, UCL, Nov 2017
4 
9 
13 
14 #include "TBranch.h"
15 #include "TCanvas.h"
16 #include "TF1.h"
17 #include "TFile.h"
18 #include "TGraph.h"
19 #include "TH1.h"
20 #include "TStyle.h"
21 #include "TTree.h"
22 #include "TVectorD.h"
23 
24 #include <iostream>
25 
26 #define PI 3.14159265
27 
28 namespace fhicl {
29  class ParameterSet;
30 }
31 
32 namespace phot {
34  public:
35  explicit CreateHybridLibrary(const fhicl::ParameterSet& p);
36 
37  // Plugins should not be copied or assigned.
40  CreateHybridLibrary& operator=(const CreateHybridLibrary&) = delete;
41  CreateHybridLibrary& operator=(CreateHybridLibrary&&) = delete;
42 
43  void analyze(const art::Event& e) override;
44  };
45 
46  //--------------------------------------------------------------------
47  CreateHybridLibrary::CreateHybridLibrary(const fhicl::ParameterSet& p) : EDAnalyzer(p)
48  {
49 
51 
53  sim::PhotonVoxelDef voxdef = pvs->GetVoxelDef();
54 
55  TFile* fout_full = new TFile("full.root", "RECREATE");
56  TFile* fout_fit = new TFile("fit.root", "RECREATE");
57 
58  std::cout << voxdef.GetNVoxels() << " voxels for each of " << geom->NOpDets() << " OpDets"
59  << std::endl;
60  std::cout << std::endl;
61 
62  //EP = Exception points: the parameterization is not a good description of the visibility and the value of the Photon Library is kept.
63  long totExceptions = 0;
64  long totPts = 0;
65 
66  for (unsigned int opdetIdx = 0; opdetIdx < geom->NOpDets(); ++opdetIdx) {
67  std::cout << opdetIdx << " / " << geom->NOpDets() << std::endl;
68 
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());
71 
72  d_full->cd();
73  TTree* tr_full = new TTree("tr", "tr");
74  int vox, taken;
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",
81  &psi); //not needed to parameterize the visibilities, useful for tests in SP
82  tr_full->Branch("theta",
83  &theta); //not needed to parameterize the visibilities, useful for tests in SP
84  tr_full->Branch("xpos",
85  &xpos); //not needed to parameterize the visibilities, useful for tests in SP
86 
87  const geo::OpDetGeo& opdet = geom->OpDetGeoFromOpDet(opdetIdx);
88  auto const opdetCenter = opdet.GetCenter();
89 
90  struct Visibility {
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)
93  {}
94  int vox;
95  int taken;
96  float dist;
97  float vis;
98  float psi;
99  float theta;
100  float xpos;
101  };
102  TCanvas* c1 = new TCanvas("c1", "c1");
103  c1->SetWindowSize(600, 600);
104  TGraph g;
105  TGraph g2;
106 
107  std::vector<Visibility> viss;
108  viss.reserve(voxdef.GetNVoxels());
109 
110  for (unsigned int voxIdx = 0U; voxIdx < voxdef.GetNVoxels(); ++voxIdx) {
111 
112  const auto voxvec = voxdef.GetPhotonVoxel(voxIdx).GetCenter();
113  const double xyzvox[] = {voxvec.X(), voxvec.Y(), voxvec.Z()};
114  const double fc_y =
115  600; //624cm is below the center of the first voxel outside the Field Cage
116  const double fc_z = 1394;
117  const double fc_x = 350;
118  taken = 0;
119  //DP does not need variable "taken" because all voxels are inside the Field Cage for the Photon Library created in LightSim.
120  //DP taken = 1;
121  dist = opdet.DistanceToPoint(voxvec);
122  vis = pvs->GetVisibility(xyzvox, opdetIdx); // TODO use Point_t when available
123  // all voxels outside the Field Cage would be assigned these values of psi and theta
124  psi = 100;
125  theta = 200;
126  xpos = voxvec.X();
127 
128  if ((voxvec.X() - opdetCenter.X()) < 0) {
129  psi = atan((voxvec.Y() - opdetCenter.Y()) / (-voxvec.X() + opdetCenter.X()));
130  }
131 
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);
135  taken = 1;
136  psi = atan((voxvec.Y() - opdetCenter.Y()) / (voxvec.X() - opdetCenter.X())) * 180.0 /
137  PI; // psi takes values within (-PI/2, PI/2)
138  theta = acos((voxvec.Z() - opdetCenter.Z()) / dist) * 180.0 /
139  PI; // theta takes values within (0 (beam direction, z), PI (-beam direction, -z))
140 
141  if ((voxvec.X() - opdetCenter.X()) < 0) {
142  psi = atan((voxvec.Y() - opdetCenter.Y()) / (-voxvec.X() + opdetCenter.X()));
143  }
144  }
145  tr_full->Fill();
146  viss.emplace_back(voxIdx, taken, dist, vis, psi, theta, xpos);
147  } // end for voxIdx
148 
149  d_full->cd();
150  tr_full->Write();
151  delete tr_full;
152 
153  g.SetMarkerStyle(7);
154  c1->cd();
155  g.Draw("ap");
156  g.GetXaxis()->SetTitle("Distance (cm)");
157  g.GetYaxis()->SetTitle("Visibility #times r^{2}");
158  g.Fit("expo");
159  TF1* fit = g.GetFunction("expo");
160 
161  d_fit->cd();
162  TVectorD fitres(2);
163  fitres[0] = fit->GetParameter(0);
164  fitres[1] = fit->GetParameter(1);
165  fitres.Write("fit");
166 
167  gPad->SetLogy();
168 
169  TH1F h("", "", 200, 0, 20);
170 
171  d_fit->cd();
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);
180 
181  for (const Visibility& v : viss) {
182  // 2e-5 is the magic scaling factor to get back to integer photon
183  // counts. TODO this will differ for new libraries, should work out a
184  // way to communicate it or derive it.
185  const double obs = v.vis / 2e-5; //taken from the Photon Library
186  const double pred =
187  fit->Eval(v.dist) / (v.dist * v.dist) / 2e-5; //calculated with parameterization
188 
189  //Minimal amount of detected photons is 50, bc of Landau dustribution
190  //Those voxels with detected photons < 50 were set to 0
191 
192  // Log-likelihood ratio for poisson statistics
193  double chisq = 2 * (pred - obs);
194  if (obs) chisq += 2 * obs * log(obs / pred);
195 
196  vox = v.vox;
197  dist = v.dist;
198  vis = pred * 2e-5;
199  psi = v.psi;
200  theta = v.theta;
201  xpos = v.xpos;
202 
203  if (v.taken == 1) { h.Fill(chisq); }
204 
205  if (
206  chisq >
207  9) { //equivalent to more than 9 chisquare = 3 sigma //maybe play around with this cutoff
208  g2.SetPoint(g2.GetN(), v.dist, v.vis * v.dist * v.dist);
209  vis = obs * 2e-5;
210  tr_fit->Fill();
211  }
212  }
213 
214  d_fit->cd();
215  tr_fit->Write();
216  delete tr_fit;
217  g2.SetMarkerSize(1);
218  g2.SetLineWidth(3);
219  g2.SetMarkerStyle(7);
220  g2.SetMarkerColor(kRed);
221  gStyle->SetLabelSize(.04, "XY");
222  gStyle->SetTitleSize(.04, "XY");
223  c1->cd();
224  g2.Draw("p same");
225  c1->SaveAs(TString::Format("plots/Chris_vis_vs_dist_%d.png", opdetIdx).Data());
226 
227  h.GetXaxis()->SetTitle("#chi^{2}");
228  h.GetYaxis()->SetTitle("Counts");
229  gStyle->SetLabelSize(.04, "XY");
230  gStyle->SetTitleSize(.04, "XY");
231  h.Draw("hist");
232  std::cout << "Integral(90,-1)/Integral(all) = " << h.Integral(90, -1) / h.Integral(0, -1)
233  << std::endl;
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());
238 
239  delete c1;
240  } // end for opdetIdx
241 
242  std::cout << totExceptions << " exceptions from " << totPts
243  << " points = " << (100. * totExceptions) / totPts << "%" << std::endl;
244 
245  delete fout_full;
246  delete fout_fit;
247  exit(0); // We're done :)
248  }
249 
250  //--------------------------------------------------------------------
252 
254 } // namespace
Definitions of voxel data structures.
Float_t vox
Definition: plot.C:38
Point_t const & GetCenter() const
Definition: OpDetGeo.h:71
Representation of a region of space diced into voxels.
Definition: PhotonVoxels.h:58
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].
Definition: OpDetGeo.cxx:94
TCanvas * c1
Definition: plotHisto.C:7
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
parameter set interface
Float_t d
Definition: plot.C:235
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).
Definition: PhotonVoxels.h:186
unsigned int GetNVoxels() const
Returns the total number of voxels in the volume.
Float_t e
Definition: plot.C:35
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