test_impactresponse.cxx
Go to the documentation of this file.
3 
4 #include "WireCellUtil/ExecMon.h"
5 #include "WireCellUtil/Testing.h"
7 
11 
12 #include "MultiPdf.h" // local helper shared by a few tests
13 #include "TH2F.h"
14 #include "TLine.h"
15 #include "TStyle.h"
16 #include "TFile.h"
17 
18 #include <iostream>
19 #include <vector>
20 #include <algorithm>
21 
22 #include "Utils.h"
23 
24 using namespace WireCell;
25 using namespace WireCell::Test;
26 using namespace std;
27 
29  int iplane, Binning tbins,
30  const std::string& name, const std::string& title)
31 {
32  // only show bins where we think the response is
33  const double tmin = tbins.min();
34  const double tspan = 100*units::us;
35  const int ntbins = tbins.bin(tmin+tspan);
36  const double tmax = tbins.edge(ntbins);
37 
38  const int nwires = pir->nwires();
39  cerr << "Plane " << iplane << " with " << nwires << " wires\n";
40 
41  const char *uvw = "UVW";
42 
43  const double half_pitch = 0.5*pir->pitch_range();
44  const double impact_dist = pir->impact();
45 
46  const double pmin = -36*units::mm, pmax=36*units::mm;
47  const int npbins = (pmax-pmin)/impact_dist;
48 
49  // dr:
50  std::string zunit = "negative microvolt";
51  double zunitval = -units::microvolt;
52  vector<double> zextent{1.0, 1.0, 2.0};
53  if (name=="fr") {
54  zunit = "induced electrons";
55  zunitval = -units::eplus;
56  zextent = vector<double>{0.3, 0.15, 0.6};
57  }
58  std::cerr <<"zunits: " << zunit
59  << " tbinsize: " << tbins.binsize()/units::us
60  << " us\n";
61 
62 
63  // they all suck. black body sucks the least.
64  set_palette(kBlackBody);
65  //set_palette(kLightTemperature);
66  //set_palette(kRedBlue);
67  //set_palette(kTemperatureMap);
68  //set_palette(kThermometer);
69  // set_palette(kVisibleSpectrum);
70  //set_palette();
71  gStyle->SetOptStat(0);
72  TH2F* hist = new TH2F(Form("h%s_%c", name.c_str(), uvw[iplane]),
73  Form("%s, 1e-/impact %c-plane", title.c_str(), uvw[iplane]),
74  ntbins, tmin/units::us, tmax/units::us,
75  npbins, pmin/units::mm, pmax/units::mm);
76  hist->SetXTitle("time (us)");
77  hist->SetYTitle("pitch (mm)");
78  hist->SetZTitle(zunit.c_str());
79 
80  hist->GetZaxis()->SetRangeUser(-zextent[iplane], +zextent[iplane]);
81 
82  TH1F* htot = new TH1F(Form("htot%s_%c", name.c_str(), uvw[iplane]),
83  Form("%s total, 1e-/impact %c-plane", title.c_str(), uvw[iplane]),
84  npbins, pmin/units::mm, pmax/units::mm);
85  htot->SetXTitle("pitch (mm)");
86  htot->SetYTitle(Form("impact total [%s]", zunit.c_str()));
87 
88  for (double pitch = -half_pitch; pitch <= half_pitch; pitch += impact_dist) {
89  auto ir = pir->closest(pitch);
90  if (!ir) {
91  std::cerr << "No closest for pitch " << pitch << endl;
92  continue;
93  }
94  auto spec = ir->spectrum();
95  auto wave = Waveform::idft(spec);
96  pitch += 0.001*impact_dist;
97  for (int ind=0; ind < ntbins; ++ind) {
98  const double time = tbins.center(ind);
99  hist->Fill(time/units::us, pitch/units::mm, wave[ind]/zunitval);
100  htot->Fill(pitch/units::mm, wave[ind]/zunitval);
101  }
102  }
103  hist->Write();
104  hist->Draw("colz");
105 
106  mpdf.canvas.SetRightMargin(0.15);
107  mpdf.canvas.SetLeftMargin(0.15);
108 
109  TLine wline, hline;
110  wline.SetLineColorAlpha(1, 0.5);
111  wline.SetLineStyle(1);
112  hline.SetLineColorAlpha(2, 0.5);
113  hline.SetLineStyle(2);
114  for (int iwire=0; iwire<nwires/2; ++iwire) {
115  double wpitch = iwire * pir->pitch();
116  if (wpitch < pmax) {
117  wline.DrawLine(tmin/units::us, wpitch, tmax/units::us, wpitch);
118  wline.DrawLine(tmin/units::us, -wpitch, tmax/units::us, -wpitch);
119  }
120  wpitch += 0.5*pir->pitch();
121  if (wpitch < pmax) {
122  hline.DrawLine(tmin/units::us, wpitch, tmax/units::us, wpitch);
123  hline.DrawLine(tmin/units::us, -wpitch, tmax/units::us, -wpitch);
124  }
125  }
126 
127  mpdf();
128  htot->Draw("hist");
129  mpdf();
130 
131 }
132 
133 
134 
135 int main(int argc, const char* argv[])
136 {
138  pm.add("WireCellGen");
139  pm.add("WireCellSigProc");
140 
141  const int nticks = 9595;
142  const double tick = 0.5*units::us;
143  const double gain = 14.0*units::mV/units::fC;
144  const double shaping = 2.0*units::us;
145 
146  const double t0 = 0.0*units::s;
147  const double readout_time = nticks*tick;
148 
149  Binning tbins(nticks, t0, t0 + readout_time);
150 
151  string out_basename = argv[0];
152  string response_file = "ub-10-half.json.bz2";
153  if (argc > 1) {
154  response_file = argv[1];
155  };
156  if (argc > 2) {
157  out_basename = argv[2];
158  }
159  cerr << "Using response file: " << response_file << endl;
160  cerr << "Writing to " << out_basename << endl;
161 
162 
163  const std::string er_tn = "ElecResponse", rc_tn = "RCResponse";
164 
165  { // configure elecresponse
166  auto icfg = Factory::lookup_tn<IConfigurable>(er_tn);
167  auto cfg = icfg->default_configuration();
168  cfg["gain"] = gain;
169  cfg["shaping"] = shaping;
170  cfg["nticks"] = nticks;
171  cerr << "Setting: " << cfg["nticks"].asInt() << " ticks\n";
172  cfg["tick"] = tick;
173  cfg["start"] = t0;
174  icfg->configure(cfg);
175  }
176  { // configure rc response
177  auto icfg = Factory::lookup_tn<IConfigurable>(rc_tn);
178  auto cfg = icfg->default_configuration();
179  cfg["nticks"] = nticks;
180  cfg["tick"] = tick;
181  cfg["start"] = t0;
182  icfg->configure(cfg);
183  }
184  {
185  auto icfg = Factory::lookup<IConfigurable>("FieldResponse");
186  auto cfg = icfg->default_configuration();
187  cfg["filename"] = response_file;
188  icfg->configure(cfg);
189  }
190 
191  std::vector<std::string> pir_tns{ "PlaneImpactResponse:frU", "PlaneImpactResponse:frV", "PlaneImpactResponse:frW"};
192  { // configure pirs, just FR
193  for (int iplane=0; iplane<3; ++iplane) {
194  auto icfg = Factory::lookup_tn<IConfigurable>(pir_tns[iplane]);
195  auto cfg = icfg->default_configuration();
196  cfg["plane"] = iplane;
197  cfg["nticks"] = nticks;
198  cfg["tick"] = tick;
199  icfg->configure(cfg);
200  }
201  }
202  std::vector<std::string> pir_ele_tns{ "PlaneImpactResponse:frerU", "PlaneImpactResponse:frerV", "PlaneImpactResponse:frerW"};
203  { // configure pirs, FR + ER + RCRC
204  for (int iplane=0; iplane<3; ++iplane) {
205  auto icfg = Factory::lookup_tn<IConfigurable>(pir_ele_tns[iplane]);
206  auto cfg = icfg->default_configuration();
207  cfg["plane"] = iplane;
208  cfg["nticks"] = nticks;
209  cfg["tick"] = tick;
210  cfg["other_responses"][0] = er_tn;
211  cfg["other_responses"][1] = rc_tn; // double it so
212  cfg["other_responses"][2] = rc_tn; // we get RC^2
213  icfg->configure(cfg);
214  }
215  }
216 
217  //auto ifr = Factory::find_tn<IFieldResponse>("FieldResponse");
218  //const auto& fr = ifr->field_response();
219 
220  TFile* rootfile = TFile::Open(Form("%s.root", out_basename.c_str()), "recreate");
221 
222 
223  MultiPdf mpdf(out_basename.c_str());
224  for (int iplane=0; iplane<3; ++iplane) {
225  auto pir = Factory::find_tn<IPlaneImpactResponse>(pir_tns[iplane]);
226  plot_time(mpdf, pir, iplane, tbins, "fr", "Field Response");
227 
228  auto pir_ele = Factory::find_tn<IPlaneImpactResponse>(pir_ele_tns[iplane]);
229  plot_time(mpdf, pir_ele, iplane, tbins, "dr", "Detector Response");
230  }
231 
232  mpdf.close();
233 
234  cerr << "Closing ROOT file: " << rootfile->GetName() << endl;
235  rootfile->Close();
236 
237  return 0;
238 
239 }
240 
241 
static QCString name
Definition: declinfo.cpp:673
code to link reconstructed objects back to the MC truth information
const int nwires
static const double eplus
Definition: Units.h:110
void set_palette(int pal=-1)
Definition: Utils.h:8
std::shared_ptr< IPlaneImpactResponse > pointer
Access subclass facet by pointer.
Definition: IComponent.h:33
double center(int ind) const
Definition: Binning.h:86
std::string string
Definition: nybbler.cc:12
STL namespace.
static const double mV
Definition: Units.h:180
cfg
Definition: dbjson.py:29
static PluginManager & instance()
const double tick
const int nticks
double binsize() const
Definition: Binning.h:73
Binning tbins(nticks, t0, t0+readout_time)
static const double microvolt
Definition: Units.h:179
const double t0
int bin(double val) const
Definition: Binning.h:80
double min() const
Definition: Binning.h:47
static const double fC
Definition: Units.h:113
static const double mm
Definition: Units.h:55
Definition: Main.h:22
realseq_t idft(compseq_t spec)
Definition: Waveform.cxx:149
static const double s
Definition: Units.h:103
static const double us
Definition: Units.h:105
double edge(int ind) const
Definition: Binning.h:99
Plugin * add(const std::string &plugin_name, const std::string &libname="")
Add a plugin. If libname is not given, try to derive it.
void plot_time(MultiPdf &mpdf, IPlaneImpactResponse::pointer pir, int iplane, Binning tbins, const std::string &name, const std::string &title)
int main(int argc, const char *argv[])
QTextStream & endl(QTextStream &s)
const double readout_time