test_blipsource.cxx
Go to the documentation of this file.
1 #include "WireCellUtil/ExecMon.h"
6 
7 #include "TFile.h"
8 #include "TStyle.h"
9 #include "TCanvas.h"
10 #include "TGraph.h"
11 #include "TH1D.h"
12 #include "TH2F.h"
13 
14 #include <iostream>
15 
16 
17 #include "ar39.h"
18 
19 using namespace WireCell;
20 using namespace WireCell::Gen::Test;
21 
22 
23 int main(int argc, char* argv[])
24 {
25  ExecMon em("start");
27  pm.add("WireCellGen");
28 
29  const int nebins=100;
30  const int specsize = ar39_mev.size();
31  const double enemin = ar39_mev.front();
32  const double enemax = ar39_mev.back();
33  const double enebin = (enemax-enemin)/specsize;
34 
35  // note: this all assumes something close to a MIP dE/dX
36  // for the full story see: http://lar.bnl.gov/properties/#particle-pass
37  // 0.7 * 1 MeV / 23.6 eV = 29661 electrons at 500 V/cm
38  // 0.64* 1 MeV / 23.6 eV = 27118 electrons at 300 V/cm
39  const double ele_per_mev = 28000; // ballpark
40 
41  TGraph spectrum;
42  double spectot = 0.0;
43 
44  // Normal user component code should never do this
45  {
46  auto rng_cfg = Factory::lookup<IConfigurable>("Random");
47  {
48  auto cfg = rng_cfg->default_configuration();
49  rng_cfg->configure(cfg);
50  }
51  auto bs_cfg = Factory::lookup<IConfigurable>("BlipSource");
52  {
53  auto cfg = bs_cfg->default_configuration();
54  Configuration edges, pdf, nele;
55  for (int ind=0; ind<specsize; ++ind) {
56  const double e = ar39_mev[ind];
57  const double p = ar39_pdf[ind];
58  edges[ind] = e;
59  nele[ind] = e*ele_per_mev;
60  pdf[ind] = p;
61  spectot += p;
62  }
63  for (int ind=0; ind<specsize; ++ind) {
64  const double prob = ar39_pdf[ind]/(spectot*enebin);
65  const double ene = ar39_mev[ind];
66  spectrum.SetPoint(ind, ene, prob);
67 
68  }
69  Configuration ene;
70  ene["type"] = "pdf";
71  ene["edges"] = edges;
72  ene["pdf"] = pdf;
73  cfg["charge"] = ene;
74  //std::cout << cfg << std::endl;
75  bs_cfg->configure(cfg);
76  ene["edges"] = nele;
77 
78  // pos config is (2m)^2 box centered on origin
79 
80  std::cout << "Configuration for Ar39 blips:\n" << ene << std::endl;
81  }
82  }
83 
84 
85  // Normal user component code might do this but should make the
86  // name a configurable parameter.
87  auto deposrc = Factory::lookup<IDepoSource>("BlipSource");
88 
89  em("initialized");
90 
91  TFile* rootfile = TFile::Open(Form("%s.root", argv[0]), "recreate");
92  gStyle->SetOptStat(111111);
93  gStyle->SetOptFit(111111);
94  TCanvas canvas("canvas","canvas", 1000, 800);
95 
96  auto hxy = new TH2F("hxy", "Decay location (Y vs X)",
97  200, -1.0, 1.0, 200, -1.0, 1.0);
98  auto hxz = new TH2F("hxz", "Decay location (Z vs X)",
99  200, -1.0, 1.0, 200, -1.0, 1.0);
100 
101  auto ht = new TH1D("ht", "Ar39 decay time", 1000, 0, 100);
102  ht->GetXaxis()->SetTitle("time between decays [us]");
103  auto he = new TH1D("he", "Ar39 decay spectrum", nebins, 0, 1); // fixme: this needs to change to electrons!
104  he->GetXaxis()->SetTitle("beta kinetic energy (MeV)");
105  auto hde = new TH1D("hde", "Energy bin size (round-off errors)", 200, 0.01-0.0001, 0.01+0.0001);
106  { // make sure we only see round-off errors
107  double last_e = 0.0;
108  for (const double e : ar39_mev) {
109  hde->Fill(e-last_e+0.00001);
110  last_e = e;
111  }
112  }
113 
114  double last_time = 0.0;
115  for (int ind=0; ind != 10000000; ++ind) {
116  IDepo::pointer depo;
117  bool ok = (*deposrc)(depo);
118  if (!ok) {
119  std::cerr << "BlipSource failed!\n";
120  return 1;
121  }
122  if (!depo) {
123  std::cerr << "BlipSource: EOS\n";
124  break;
125  }
126  const double charge = depo->charge();
127  if (charge < 0) {
128  std::cerr << "Got negative charge\n";
129  return 1;
130  }
131  he->Fill(charge);
132  const double time = depo->time();
133  ht->Fill((time - last_time)/units::us);
134  last_time = time;
135 
136  auto& pos = depo->pos();
137  hxy->Fill(pos.x()/units::m, pos.y()/units::m, charge);
138  hxz->Fill(pos.x()/units::m, pos.z()/units::m, charge);
139 
140  }
141  he->Scale(1.0/(he->Integral() / nebins));
142 
143  em("filled");
144 
145  std::string fname = Form("%s.pdf", argv[0]);
146  canvas.Print((fname+"[").c_str(), "pdf");
147 
148  ht->Fit("expo");
149  ht->Draw();
150  canvas.Print(fname.c_str(), "pdf");
151  em("fitted");
152 
153  spectrum.SetLineColor(2);
154  he->Draw();
155  spectrum.Draw("L");
156  canvas.Print(fname.c_str(), "pdf");
157 
158  hde->Draw();
159  canvas.Print(fname.c_str(), "pdf");
160 
161  hxy->Draw("colz");
162  canvas.Print(fname.c_str(), "pdf");
163 
164  hxz->Draw("colz");
165  canvas.Print(fname.c_str(), "pdf");
166 
167  canvas.Print((fname+"]").c_str(), "pdf");
168  em("plotted");
169 
170  ht->Write();
171  he->Write();
172  hde->Write();
173  rootfile->Close();
174 
175  em("done");
176  std::cerr << em.summary() << std::endl;
177  return 0;
178 }
std::shared_ptr< const IDepo > pointer
Definition: IData.h:19
static const double m
Definition: Units.h:63
std::string string
Definition: nybbler.cc:12
const std::vector< double > ar39_pdf
Definition: ar39.h:50
cfg
Definition: dbjson.py:29
static PluginManager & instance()
const double e
Definition: Main.h:22
p
Definition: test.py:223
const std::vector< double > ar39_mev
Definition: ar39.h:35
Json::Value Configuration
Definition: Configuration.h:50
static const double us
Definition: Units.h:105
int main(int argc, char *argv[])
Plugin * add(const std::string &plugin_name, const std::string &libname="")
Add a plugin. If libname is not given, try to derive it.
QTextStream & endl(QTextStream &s)