LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
CheckProtonProduction_module.cc
Go to the documentation of this file.
1 //
2 // __ __ __ __ __
3 // ____ ______/ /_____ _/ // / / /_/ /__
4 // / __ `/ ___/ __/ __ `/ // /_/ __/ //_/
5 // / /_/ / / / /_/ /_/ /__ __/ /_/ ,<
6 // \__,_/_/ \__/\__, / /_/ \__/_/|_|
7 // /____/
8 //
9 // artg4tk: art based Geant 4 Toolkit
10 //
11 //=============================================================================
12 // CheckProtonProduction_module.cc: Analyzer module that demonstrates access to
13 // and makes some histograms
14 //
15 // Author: Hans Wenzel (Fermilab)
16 //=============================================================================
17 
18 // C++ includes.
19 #include <cmath>
20 #include <iostream>
21 #include <sstream>
22 #include <string>
23 #include <vector>
24 
25 // art Framework includes.
26 #include "art/Framework/Core/EDAnalyzer.h"
27 #include "art/Framework/Core/ModuleMacros.h"
28 #include "art/Framework/Principal/Event.h"
29 #include "art/Framework/Principal/Handle.h"
30 #include "art/Framework/Principal/Run.h"
31 #include "art/Framework/Services/Registry/ServiceHandle.h"
32 #include "art_root_io/TFileService.h"
33 
34 // artg4tk includes:
35 #include "artg4tk/pluginDetectors/gdml/myInteractionArtHitData.hh"
36 
37 // Root includes.
38 #include "TH1F.h"
39 
40 using namespace std;
41 namespace artg4tk {
42  class CheckProtonProduction;
43 }
44 
45 class artg4tk::CheckProtonProduction : public art::EDAnalyzer {
46 public:
47  explicit CheckProtonProduction(fhicl::ParameterSet const& p);
48  void beginJob() override;
49  void beginRun(const art::Run& Run) override;
50  void endJob() override;
51  void analyze(const art::Event& event) override;
52 
53 private:
54  int fNThetaBinsFW;
55  double fThetaMinFW;
56  double fDeltaThetaFW;
57  int fNThetaBinsLA;
58  double fThetaMinLA;
59  double fDeltaThetaLA;
60  //
61  TH1F* _fHistoNSec; // number of secondaries
62  //
63  TH1F* _fHistoSecPMom; // momentum of secondary protons
64  TH1F* _fHistoSecPkinE54; // kinetic Energy of secondary protons
65  TH1F* _fHistoSecPkinE68; // kinetic Energy of secondary protons
66  TH1F* _fHistoSecPkinE90; // kinetic Energy of secondary protons
67  TH1F* _fHistoSecPkinE121; // kinetic Energy of secondary protons
68  TH1F* _fHistoSecPkinE164; // kinetic Energy of secondary protons
69  TH1F* _fHistoSecPkinE; // kinetic Energy of secondary protons
70  TH1F* _fHistoSecPTheta; // theta of secondary protons
71  //
72  TH1F* _fHistoSecDMom; // momentum of secondary deuterons
73  TH1F* _fHistoSecDkinE54; // kinetic Energy of secondary deuterons
74  TH1F* _fHistoSecDkinE68; // kinetic Energy of secondary deuterons
75  TH1F* _fHistoSecDkinE90; // kinetic Energy of secondary deuterons
76  TH1F* _fHistoSecDkinE121; // kinetic Energy of secondary deuterons
77  TH1F* _fHistoSecDkinE164; // kinetic Energy of secondary deuterons
78  TH1F* _fHistoSecDkinE; // kinetic Energy of secondary deuterons
79  TH1F* _fHistoSecDTheta; // theta of secondary deuterons
80  //
81  TH1F* _fHistoSecTMom; // momentum of secondary tritons
82  TH1F* _fHistoSecTkinE54; // kinetic Energy of secondary tritons
83  TH1F* _fHistoSecTkinE68; // kinetic Energy of secondary tritons
84  TH1F* _fHistoSecTkinE90; // kinetic Energy of secondary tritons
85  TH1F* _fHistoSecTkinE121; // kinetic Energy of secondary tritons
86  TH1F* _fHistoSecTkinE164; // kinetic Energy of secondary tritons
87  TH1F* _fHistoSecTkinE; // kinetic Energy of secondary tritons
88  TH1F* _fHistoSecTTheta; // theta of secondary tritons
89  std::vector<TH1D*> fHistoSecProtonFW;
90 };
91 
92 artg4tk::CheckProtonProduction::CheckProtonProduction(fhicl::ParameterSet const& p)
93  : art::EDAnalyzer(p)
94  , fNThetaBinsFW(p.get<int>("ThetaBinsFW", 4))
95  , fThetaMinFW(p.get<double>("ThetaMinFW", 0.05))
96  , fDeltaThetaFW(p.get<double>("DeltaThetaFW", 0.05))
97  , fNThetaBinsLA(p.get<int>("ThetaBinsLA", 9))
98  , fThetaMinLA(p.get<double>("ThetaMinLA", 0.35))
99  , fDeltaThetaLA(p.get<double>("DeltaThetaLA", 0.2))
100  , _fHistoNSec(0)
101  , _fHistoSecPMom(0)
102  , _fHistoSecPkinE54(0)
103  , _fHistoSecPkinE68(0)
104  , _fHistoSecPkinE90(0)
105  , _fHistoSecPkinE121(0)
106  , _fHistoSecPkinE164(0)
107  , _fHistoSecPkinE(0)
108  , _fHistoSecPTheta(0)
109  , _fHistoSecDMom(0)
110  , _fHistoSecDkinE54(0)
111  , _fHistoSecDkinE68(0)
112  , _fHistoSecDkinE90(0)
113  , _fHistoSecDkinE121(0)
114  , _fHistoSecDkinE164(0)
115  , _fHistoSecDkinE(0)
116  , _fHistoSecDTheta(0)
117  , _fHistoSecTMom(0)
118  , _fHistoSecTkinE54(0)
119  , _fHistoSecTkinE68(0)
120  , _fHistoSecTkinE90(0)
121  , _fHistoSecTkinE121(0)
122  , _fHistoSecTkinE164(0)
123  , _fHistoSecTkinE(0)
124  , _fHistoSecTTheta(0)
125 {}
126 
127 void
128 artg4tk::CheckProtonProduction::beginRun(const art::Run& thisRun)
129 {
130  thisRun.beginTime();
131 }
132 
133 void
134 artg4tk::CheckProtonProduction::beginJob()
135 {
136  art::ServiceHandle<art::TFileService> tfs;
137  _fHistoNSec = tfs->make<TH1F>("NSec", "neutron + Cu", 100, 0., 100.);
138  //
139  _fHistoSecPMom =
140  tfs->make<TH1F>("pmomentum", "neutron + Cu momentum of secondary protons", 100, 0., 600);
141  _fHistoSecPkinE54 = tfs->make<TH1F>(
142  "pkinE54", "neutron + Cu kinetic Energy of secondary protons (54)", 100, 0., 400.);
143  _fHistoSecPkinE68 = tfs->make<TH1F>(
144  "pkinE68", "neutron + Cu kinetis Energy of secondary protons (68)", 100, 0., 400.);
145  _fHistoSecPkinE90 = tfs->make<TH1F>(
146  "pkinE90", "neutron + Cu kinetic Energy of secondary protons (90)", 100, 0., 400.);
147  _fHistoSecPkinE121 = tfs->make<TH1F>(
148  "pkinE121", "neutron + Cu kinetic Energy of secondary protons (121)", 100, 0., 400);
149  _fHistoSecPkinE164 = tfs->make<TH1F>(
150  "pkinE164", "neutron + Cu kinetic Energy of secondary protons (168)", 100, 0., 400.);
151  _fHistoSecPkinE =
152  tfs->make<TH1F>("pkinE", "neutron + Cu kinetic Energy of secondary protons", 100, 0., 400.);
153  _fHistoSecPTheta =
154  tfs->make<TH1F>("ptheta", "neutron + Cu kin. theta of secondary protons", 100, 0., 3.2);
155  //
156  _fHistoSecDMom =
157  tfs->make<TH1F>("dmomentum", "neutron + Cu momentum of secondary deuterons", 100, 0., 600);
158  _fHistoSecDkinE54 = tfs->make<TH1F>(
159  "dkinE54", "neutron + Cu kinetic Energy of secondary deuterons (54)", 100, 0., 400.);
160  _fHistoSecDkinE68 = tfs->make<TH1F>(
161  "dkinE68", "neutron + Cu kinetis Energy of secondary deuterons (68)", 100, 0., 400.);
162  _fHistoSecDkinE90 = tfs->make<TH1F>(
163  "dkinE90", "neutron + Cu kinetic Energy of secondary deuterons (90)", 100, 0., 400.);
164  _fHistoSecDkinE121 = tfs->make<TH1F>(
165  "dkinE121", "neutron + Cu kinetic Energy of secondary deuterons (121)", 100, 0., 400);
166  _fHistoSecDkinE164 = tfs->make<TH1F>(
167  "dkinE164", "neutron + Cu kinetic Energy of secondary deuterons (168)", 100, 0., 400.);
168  _fHistoSecDkinE =
169  tfs->make<TH1F>("dkinE", "neutron + Cu kinetic Energy of secondary deuterons", 100, 0., 400.);
170  _fHistoSecDTheta =
171  tfs->make<TH1F>("dtheta", "neutron + Cu kin. theta of secondary deuterons", 100, 0., 3.2);
172  //
173  _fHistoSecTMom =
174  tfs->make<TH1F>("tmomentum", "neutron + Cu momentum of secondary tritons", 100, 0., 600);
175  _fHistoSecTkinE54 = tfs->make<TH1F>(
176  "tkinE54", "neutron + Cu kinetic Energy of secondary tritons (54)", 100, 0., 400.);
177  _fHistoSecTkinE68 = tfs->make<TH1F>(
178  "tkinE68", "neutron + Cu kinetis Energy of secondary tritons (68)", 100, 0., 400.);
179  _fHistoSecTkinE90 = tfs->make<TH1F>(
180  "tkinE90", "neutron + Cu kinetic Energy of secondary tritons (90)", 100, 0., 400.);
181  _fHistoSecTkinE121 = tfs->make<TH1F>(
182  "tkinE121", "neutron + Cu kinetic Energy of secondary tritons (121)", 100, 0., 400);
183  _fHistoSecTkinE164 = tfs->make<TH1F>(
184  "tkinE164", "neutron + Cu kinetic Energy of secondary tritons (168)", 100, 0., 400.);
185  _fHistoSecTkinE =
186  tfs->make<TH1F>("tkinE", "neutron + Cu kinetic Energy of secondary tritons", 100, 0., 400.);
187  _fHistoSecTTheta =
188  tfs->make<TH1F>("ttheta", "neutron + Cu kin. theta of secondary tritons", 100, 0., 3.2);
189  //
190  string htitle;
191  string hname;
192  string ht = "neutron + Cu";
193 
194  double thetaMin = 0.;
195  double thetaMax = 0.;
196  std::string theta_bin_fw;
197  std::string theta_bin_la;
198 
199  for (int i = 0; i < fNThetaBinsFW; i++) {
200  thetaMin = fThetaMinFW + fDeltaThetaFW * i;
201  thetaMax = thetaMin + fDeltaThetaFW;
202 
203  std::ostringstream osTitle1;
204  std::ostringstream osTitle2;
205  std::ostringstream osTitle3;
206 
207  osTitle1.clear();
208  osTitle1 << thetaMin;
209  theta_bin_fw = osTitle1.str() + " < theta < ";
210  osTitle2.clear();
211  osTitle2 << thetaMax;
212  theta_bin_fw += osTitle2.str();
213  theta_bin_fw += "(rad)";
214 
215  osTitle3.clear();
216  osTitle3 << i;
217 
218  htitle = ht + " -> X + proton, " + theta_bin_fw;
219  hname = "proton_fw_" + osTitle3.str();
220  TH1D* histo = tfs->make<TH1D>(hname.c_str(), htitle.c_str(), 80, 0., 8.0);
221  fHistoSecProtonFW.push_back(histo);
222  }
223 
224 } // end beginJob
225 
226 void
227 artg4tk::CheckProtonProduction::analyze(const art::Event& event)
228 {
229  typedef std::vector<art::Handle<myInteractionArtHitDataCollection>> HandleVector;
230  auto allSims = event.getMany<myInteractionArtHitDataCollection>();
231 
232  for (HandleVector::const_iterator i = allSims.begin(); i != allSims.end(); ++i) {
233  const myInteractionArtHitDataCollection& sims(**i);
234  if (sims.size() > 0) {
235  _fHistoNSec->Fill(sims.size());
236  }
237  for (myInteractionArtHitDataCollection::const_iterator j = sims.begin(); j != sims.end(); ++j) {
238  const myInteractionArtHitData& hit = *j;
239  if (hit.pname == "proton") {
240  _fHistoSecPMom->Fill(hit.pmom);
241  double angle = hit.theta * 57.29577951;
242  if (fabs(angle - 54.) < 2.)
243  _fHistoSecPkinE54->Fill(hit.Ekin);
244  if (fabs(angle - 68.) < 2.)
245  _fHistoSecPkinE68->Fill(hit.Ekin);
246  if (fabs(angle - 90.) < 2.)
247  _fHistoSecPkinE90->Fill(hit.Ekin);
248  if (fabs(angle - 121.) < 2.)
249  _fHistoSecPkinE121->Fill(hit.Ekin);
250  if (fabs(angle - 164.) < 2.)
251  _fHistoSecPkinE164->Fill(hit.Ekin);
252  _fHistoSecPkinE->Fill(hit.Ekin);
253  _fHistoSecPTheta->Fill(hit.theta);
254  } else if (hit.pname == "deuteron") {
255  _fHistoSecDMom->Fill(hit.pmom);
256  double angle = hit.theta * 57.29577951;
257  if (fabs(angle - 54.) < 2.)
258  _fHistoSecDkinE54->Fill(hit.Ekin);
259  if (fabs(angle - 68.) < 2.)
260  _fHistoSecDkinE68->Fill(hit.Ekin);
261  if (fabs(angle - 90.) < 2.)
262  _fHistoSecDkinE90->Fill(hit.Ekin);
263  if (fabs(angle - 121.) < 2.)
264  _fHistoSecDkinE121->Fill(hit.Ekin);
265  if (fabs(angle - 164.) < 2.)
266  _fHistoSecDkinE164->Fill(hit.Ekin);
267  _fHistoSecDkinE->Fill(hit.Ekin);
268  _fHistoSecDTheta->Fill(hit.theta);
269  } else if (hit.pname == "triton") {
270  _fHistoSecTMom->Fill(hit.pmom);
271  double angle = hit.theta * 57.29577951;
272  if (fabs(angle - 54.) < 2.)
273  _fHistoSecTkinE54->Fill(hit.Ekin);
274  if (fabs(angle - 68.) < 2.)
275  _fHistoSecTkinE68->Fill(hit.Ekin);
276  if (fabs(angle - 90.) < 2.)
277  _fHistoSecTkinE90->Fill(hit.Ekin);
278  if (fabs(angle - 121.) < 2.)
279  _fHistoSecTkinE121->Fill(hit.Ekin);
280  if (fabs(angle - 164.) < 2.)
281  _fHistoSecTkinE164->Fill(hit.Ekin);
282  _fHistoSecTkinE->Fill(hit.Ekin);
283  _fHistoSecTTheta->Fill(hit.theta);
284  }
285  }
286  }
287 
288 } // end analyze
289 
290 void
291 artg4tk::CheckProtonProduction::endJob()
292 {
293  cout << " ********************************CheckProtonProduction: now normalizing the histos "
294  << endl;
295  double norm = _fHistoNSec->Integral();
296  double xbin = _fHistoNSec->GetBinWidth(1);
297  double scale = 1. / (xbin * norm);
298  _fHistoNSec->Scale(scale);
299 } // end endJob
300 
301 DEFINE_ART_MODULE(artg4tk::CheckProtonProduction)
artg4tk::CheckProtonProduction::analyze
void analyze(const art::Event &event) override
Definition: CheckProtonProduction_module.cc:227
std
STL namespace.
const_iterator
intermediate_table::const_iterator const_iterator
Definition: intermediate_table.cc:28
art::Run
Definition: Run.h:37
scale
Double_t scale
Definition: plot.C:24
artg4tk::CheckProtonProduction::CheckProtonProduction
CheckProtonProduction(fhicl::ParameterSet const &p)
Definition: CheckProtonProduction_module.cc:92
DEFINE_ART_MODULE
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
breakpoints::beginJob
void beginJob()
Definition: Breakpoints.cc:14
artg4tk::CheckProtonProduction::beginRun
void beginRun(const art::Run &Run) override
Definition: CheckProtonProduction_module.cc:128
Event.h
hit
Detector simulation of raw signals on wires.
Definition: DumpHits_module.cc:25
util::get
decltype(auto) get(T &&obj)
ADL-aware version of std::to_string.
Definition: StdUtils.h:120
norm
Float_t norm
Definition: comparison_ascii.C:71
art
Definition: MVAAlg.h:12
art::Run::beginTime
Timestamp const & beginTime() const
Definition: Run.cc:33
event
Event finding and building.
Definition: EventCheater_module.cc:32