All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions | List of all members
phot::CreateHybridLibrary Class Reference
Inheritance diagram for phot::CreateHybridLibrary:

Public Member Functions

 CreateHybridLibrary (const fhicl::ParameterSet &p)
 
 CreateHybridLibrary (const CreateHybridLibrary &)=delete
 
 CreateHybridLibrary (CreateHybridLibrary &&)=delete
 
CreateHybridLibrary & operator= (const CreateHybridLibrary &)=delete
 
CreateHybridLibrary & operator= (CreateHybridLibrary &&)=delete
 
void analyze (const art::Event &e) override
 

Detailed Description

Definition at line 32 of file CreateHybridLibrary_module.cc.

Constructor & Destructor Documentation

phot::CreateHybridLibrary::CreateHybridLibrary ( const fhicl::ParameterSet &  p)
explicit

Definition at line 47 of file CreateHybridLibrary_module.cc.

48  : EDAnalyzer(p)
49  {
50 
51  art::ServiceHandle<geo::Geometry const> geom;
52 
53  art::ServiceHandle<phot::PhotonVisibilityService const> pvs;
54  sim::PhotonVoxelDef voxdef = pvs->GetVoxelDef();
55 
56  TFile* fout_full = new TFile("full.root", "RECREATE");
57  TFile* fout_fit = new TFile("fit.root", "RECREATE");
58 
59  std::cout << voxdef.GetNVoxels() << " voxels for each of " << geom->NOpDets() << " OpDets" << 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", &psi); //not needed to parameterize the visibilities, useful for tests in SP
81  tr_full->Branch("theta", &theta); //not needed to parameterize the visibilities, useful for tests in SP
82  tr_full->Branch("xpos", &xpos); //not needed to parameterize the visibilities, useful for tests in SP
83 
84  const geo::OpDetGeo& opdet = geom->OpDetGeoFromOpDet(opdetIdx);
85  auto const opdetCenter = opdet.GetCenter();
86  //std::cout << "OpDet position " << opdetIdx << ": " << opdetCenter << std::endl;
87 
88  struct Visibility{
89  Visibility(int vx, int t, float d, float v, float p, float th, float xp) : vox(vx), taken(t), dist(d), vis(v), psi(p), theta(th), xpos(xp) {}
90  int vox;
91  int taken;
92  float dist;
93  float vis;
94  float psi;
95  float theta;
96  float xpos;
97  };
98  TCanvas *c1=new TCanvas("c1", "c1");
99  //c1->SetCanvasSize(1500, 1500);
100  c1->SetWindowSize(600, 600);
101  //c1->Divide(1,2);
102  TGraph g;
103  TGraph g2;
104 
105  std::vector<Visibility> viss;
106  viss.reserve(voxdef.GetNVoxels());
107 
108  for(unsigned int voxIdx = 0U; voxIdx < voxdef.GetNVoxels(); ++voxIdx){
109 
110  const auto voxvec = voxdef.GetPhotonVoxel(voxIdx).GetCenter();
111  const double xyzvox[] = {voxvec.X(), voxvec.Y(), voxvec.Z()};
112  const double fc_y = 600; //624cm is below the center of the first voxel outside the Field Cage
113  const double fc_z = 1394;
114  const double fc_x = 350;
115  taken = 0;
116  //DP does not need variable "taken" because all voxels are inside the Field Cage for the Photon Library created in LightSim.
117  //DP taken = 1;
118  dist = opdet.DistanceToPoint(voxvec);
119  vis = pvs->GetVisibility(xyzvox, opdetIdx); // TODO use Point_t when available
120  // all voxels outside the Field Cage would be assigned these values of psi and theta
121  psi = 100;
122  theta = 200;
123  xpos = voxvec.X();
124 
125  if((voxvec.X() - opdetCenter.X())<0){
126  psi = atan((voxvec.Y() - opdetCenter.Y())/(-voxvec.X() +opdetCenter.X()));
127  }
128 
129  if(voxvec.X()<fc_x && voxvec.X()>-fc_x && voxvec.Y()<fc_y && voxvec.Y()>-fc_y && voxvec.Z()> -9 && voxvec.Z()<fc_z){
130  g.SetPoint(g.GetN(), dist, vis*dist*dist);
131  taken = 1;
132  psi = atan((voxvec.Y() - opdetCenter.Y())/(voxvec.X() - opdetCenter.X()))* 180.0/PI ; // psi takes values within (-PI/2, PI/2)
133  theta = acos((voxvec.Z() - opdetCenter.Z())/dist) * 180.0/PI; // theta takes values within (0 (beam direction, z), PI (-beam direction, -z))
134 
135  if((voxvec.X() - opdetCenter.X())<0){
136  psi = atan((voxvec.Y() - opdetCenter.Y())/(-voxvec.X() +opdetCenter.X()));
137  }
138 
139  }
140  tr_full->Fill();
141  viss.emplace_back(voxIdx, taken, dist, vis, psi, theta, xpos);
142  } // end for voxIdx
143 
144  d_full->cd();
145  tr_full->Write();
146  delete tr_full;
147 
148  g.SetMarkerStyle(7);
149  c1->cd();
150  g.Draw("ap");
151  g.GetXaxis()->SetTitle("Distance (cm)");
152  g.GetYaxis()->SetTitle("Visibility #times r^{2}");
153  g.Fit("expo");
154  TF1* fit = g.GetFunction("expo");
155 
156  d_fit->cd();
157  TVectorD fitres(2);
158  fitres[0] = fit->GetParameter(0);
159  fitres[1] = fit->GetParameter(1);
160  fitres.Write("fit");
161 
162  gPad->SetLogy();
163 
164  TH1F h("", "", 200, 0, 20);
165 
166  d_fit->cd();
167  TTree* tr_fit = new TTree("tr", "tr");
168  tr_fit->Branch("vox", &vox);
169  tr_fit->Branch("dist", &dist);
170  tr_fit->Branch("vis", &vis);
171  tr_fit->Branch("taken", &taken);
172  tr_fit->Branch("psi", &psi);
173  tr_fit->Branch("theta", &theta);
174  tr_fit->Branch("xpos", &xpos);
175 
176 
177  for(const Visibility& v: viss){
178  // 2e-5 is the magic scaling factor to get back to integer photon
179  // counts. TODO this will differ for new libraries, should work out a
180  // way to communicate it or derive it.
181  const double obs = v.vis / 2e-5; //taken from the Photon Library
182  const double pred = fit->Eval(v.dist) / (v.dist*v.dist) / 2e-5; //calculated with parameterization
183 
184  //DP const double obs = v.vis / 1e-7; //magic scaling factor for DP library created in LightSim
185  //Minimal amount of detected photons is 50, bc of Landau dustribution
186  //Those voxels with detected photons < 50 were set to 0
187  //DP const double pred = fit->Eval(v.dist) / (v.dist*v.dist) / 1e-7; //calculated with parametrisation
188  //std::cout << "observed = "<<obs<<" predicted (by parametrization) = "<<pred <<std::endl;
189 
190 
191  // Log-likelihood ratio for poisson statistics
192  double chisq = 2*(pred-obs);
193  if(obs) chisq += 2*obs*log(obs/pred);
194 
195  vox = v.vox;
196  dist = v.dist;
197  vis = pred *2e-5;
198  psi = v.psi;
199  theta = v.theta;
200  xpos = v.xpos;
201  //DP vis = pred *1e-7;
202 
203 
204  if (v.taken==1){
205  h.Fill(chisq);
206  }
207 
208 
209  if(chisq > 9){ //equivalent to more than 9 chisquare = 3 sigma //maybe play around with this cutoff
210  g2.SetPoint(g2.GetN(), v.dist, v.vis*v.dist*v.dist);
211  vis = obs *2e-5;
212  //DP vis = obs *1e-7;
213  tr_fit->Fill();
214  }
215 
216  }
217 
218  d_fit->cd();
219  tr_fit->Write();
220  delete tr_fit;
221  g2.SetMarkerSize(1);
222  g2.SetLineWidth(3);
223  g2.SetMarkerStyle(7);
224  g2.SetMarkerColor(kRed);
225  gStyle->SetLabelSize(.04, "XY");
226  gStyle->SetTitleSize(.04,"XY");
227  c1->cd();
228  g2.Draw("p same");
229  c1->SaveAs(TString::Format("plots/Chris_vis_vs_dist_%d.png", opdetIdx).Data());
230 
231 
232  h.GetXaxis()->SetTitle("#chi^{2}");
233  h.GetYaxis()->SetTitle("Counts");
234  gStyle->SetLabelSize(.04, "XY");
235  gStyle->SetTitleSize(.04,"XY");
236  h.Draw("hist");
237  std::cout <<"Integral(90,-1)/Integral(all) = "<< h.Integral(90, -1) / h.Integral(0, -1) << std::endl;
238  std::cout <<"*****************************" << std::endl;
239  totExceptions += h.Integral(90, -1);
240  totPts += h.Integral(0, -1);
241  gPad->Print(TString::Format("plots/chisq_opdet_%d.eps", opdetIdx).Data());
242 
243  delete c1;
244  } // end for opdetIdx
245 
246  std::cout << totExceptions << " exceptions from " << totPts << " points = " << (100.*totExceptions)/totPts << "%" << std::endl;
247 
248  delete fout_full;
249  delete fout_fit;
250  exit(0); // We're done :)
251  }
exit
exit
Definition: generate_icarus_geometry_splitwire.pl:430
PI
#define PI
Definition: CreateHybridLibrary_module.cc:26
p
pdgs p
Definition: selectors.fcl:22
sim::PhotonVoxelDef
Representation of a region of space diced into voxels.
Definition: larsim/larsim/Simulation/PhotonVoxels.h:58
geo::OpDetGeo::GetCenter
void GetCenter(double *xyz, double localz=0.0) const
Definition: OpDetGeo.cxx:40
g
BEGIN_PROLOG g
Definition: decoderTools_icarus.fcl:7
sim::PhotonVoxelDef::GetNVoxels
unsigned int GetNVoxels() const
Returns the total number of voxels in the volume.
Definition: larsim/larsim/Simulation/PhotonVoxels.cxx:58
h
while getopts h
Definition: getStatistics.sh:13
sim::PhotonVoxel::GetCenter
Point GetCenter() const
Returns the center of the voxel (type Point).
Definition: larsim/larsim/Simulation/PhotonVoxels.h:199
geo::OpDetGeo::DistanceToPoint
double DistanceToPoint(geo::Point_t const &point) const
Returns the distance of the specified point from detector center [cm].
Definition: OpDetGeo.cxx:123
larg4::dist
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
Definition: OpFastScintillation.hh:563
e
do i e
Definition: updateVersionList.sh:68
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10
phot::CreateHybridLibrary::CreateHybridLibrary ( const CreateHybridLibrary &  )
delete
phot::CreateHybridLibrary::CreateHybridLibrary ( CreateHybridLibrary &&  )
delete

Member Function Documentation

void phot::CreateHybridLibrary::analyze ( const art::Event &  e)
override

Definition at line 254 of file CreateHybridLibrary_module.cc.

255  {
256  }
CreateHybridLibrary& phot::CreateHybridLibrary::operator= ( const CreateHybridLibrary &  )
delete
CreateHybridLibrary& phot::CreateHybridLibrary::operator= ( CreateHybridLibrary &&  )
delete
The documentation for this class was generated from the following file:
Generated by   doxygen 1.8.5