Belle II Software prerelease-10-00-00a
TOPBackgroundModule.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9// Own header.
10#include <top/modules/TOPBackground/TOPBackgroundModule.h>
11
12// TOP headers.
13#include <top/dataobjects/TOPBarHit.h>
14#include <top/dataobjects/TOPDigit.h>
15#include <top/dataobjects/TOPSimHit.h>
16#include <top/geometry/TOPGeometryPar.h>
17
18#include <mdst/dataobjects/MCParticle.h>
19#include <simulation/dataobjects/BeamBackHit.h>
20
21// framework - DataStore
22#include <framework/datastore/DataStore.h>
23#include <framework/datastore/StoreArray.h>
24
25// framework aux
26#include <framework/logging/Logger.h>
27#include <framework/gearbox/Const.h>
28
29#include <TCanvas.h>
30#include <TLegend.h>
31#include <TLine.h>
32#include <TPad.h>
33#include <TROOT.h>
34#include <TStyle.h>
35#include <Math/VectorUtil.h>
36#include <TMath.h>
37
38#include <geometry/Materials.h>
39#include <G4Material.hh>
40
41using namespace std;
42
43namespace Belle2 {
48 //-----------------------------------------------------------------
50 //-----------------------------------------------------------------
51
52 REG_MODULE(TOPBackground);
53
54
55 //-----------------------------------------------------------------
56 // Implementation
57 //-----------------------------------------------------------------
58
60 m_rootFile(0),
61 origingamma(0),
62 originpe(0),
63 peflux(0),
64 nflux(0),
65 rdose(0),
66 zdist(0),
67 genergy(0),
68 genergy2(0),
69 zdistg(0),
70 originpt(0),
71 nflux_bar(0),
72 gflux_bar(0),
73 cflux_bar(0),
74 norigin(0),
75 gorigin(0),
76 corigin(0),
77 nprim(0),
78 gprim(0),
79 cprim(0),
80 origin_zx(0),
81 origin_zy(0),
82 prim_zx(0),
83 prim_zy(0),
85 PCBmass(0),
86 PCBarea(0),
87 yearns(0),
88 evtoJ(0),
89 mtoc(0),
90 count(0),
91 count_occ(0),
95 originpe_x(0),
96 originpe_y(0),
97 originpe_z(0)
98 {
99 // Set description()
100 setDescription("A module to analyze beam background simulations regarding TOP");
101
102 // Add parameters
103 addParam("Type", m_BkgType, "Backgound type", string("Backgound"));
104 addParam("Output", m_OutputFileName, "Name of the output file",
105 string("Backgound.root"));
106 addParam("TimeOfSimulation", m_TimeOfSimulation,
107 "Real time in micro seconds that corresponds to simulated data", 5.);
108 }
109
114
116 {
117
118 // CPU time start
119
120 // Initializing the output root file
121 m_rootFile = TFile::Open(m_OutputFileName.c_str(), "RECREATE");
122 origingamma = new TTree("origingamma", "tree");
123 originpe = new TTree("originpe", "tree2");
124
125 origingamma->Branch("x", &origingamma_x);
126 origingamma->Branch("y", &origingamma_y);
127 origingamma->Branch("z", &origingamma_z);
128 originpe->Branch("x", &originpe_x);
129 originpe->Branch("y", &originpe_y);
130 originpe->Branch("z", &originpe_z);
131
132 peflux = new TH1F("Photoelectron flux", "Photoelectron flux", 33, -11.25, 360);
133 nflux = new TH1F("Neutron flux", "Neutron flux", 33, -11.25, 360);
134 rdose = new TH1F("Radiation dose", "Radiation dose", 33, -11.25, 360);
135 zdist = new TH1F("Photoelectron origin", "Photoelectron origin", 200, -400, 400);
136 genergy = new TH1F("Energy distribution of photons", "Energy distribution of photons", 150, 0, 5);
137 genergy2 = new TH1F("Energy distribution of gammas", "Energy distribution of gammas", 500, 0, 5);
138
139 zdistg = new TH1F("Photoelectron flux z", "Photoelectron flux Z projection", 800, -400, 400);
140 originpt = new TH1F("P_t of origin electron", "P_t of origin electron", 300, 0.06, 0.14);
141
142 nflux_bar = new TH2F("Neutron flux on bar", "Neutron flux on bar", 32, -114.8, 211.5, 16, -0, 360);
143 gflux_bar = new TH2F("Gamma flux on bar", "Gamma flux on bar MHz/cm^{2}", 32, -114.8, 211.5, 16, -0, 360);
144 cflux_bar = new TH2F("Charged flux on bar", "Charged flux on bar MHz/cm^{2}", 32, -114.8, 211.5, 16, -0, 360);
145
146 norigin = new TH1F("neutron(BAR) origin", "neutron(BAR) origin", 200, -400, 400);
147 gorigin = new TH1F("gamma(BAR) origin", "gamma(BAR) origin", 200, -400, 400);
148 corigin = new TH1F("charged(BAR) origin", "charged(BAR) origin", 200, -400, 400);
149
150 nprim = new TH1F("neutron(BAR) primary", "neutron(BAR) primary", 200, -400, 400);
151 gprim = new TH1F("gamma(BAR) primary", "gamma(BAR) primary", 200, -400, 400);
152 cprim = new TH1F("charged(BAR) primary", "charged(BAR) primary", 200, -400, 400);
153
154 origin_zx = new TGraph();
155 origin_zy = new TGraph();
156
157 prim_zx = new TGraph();
158 prim_zy = new TGraph();
159 module_occupancy = new TGraph();
160
161 origin_zy->SetName("originZY");
162 origin_zx->SetName("originZX");
163 module_occupancy->SetName("occupancy");
164
165 const auto& geo = TOP::TOPGeometryPar::Instance()->getGeometry()->getFrontEnd();
166 double S1 = geo.getFrontBoardWidth() * geo.getFrontBoardHeight();
167 double S2 = geo.getHVBoardWidth() * geo.getHVBoardLength();
168 double V1 = S1 * geo.getFrontBoardThickness();
169 double V2 = S2 * geo.getHVBoardThickness();
170 G4Material* material1 = geometry::Materials::get(geo.getFrontBoardMaterial());
171 if (!material1) B2FATAL("Material '" << geo.getFrontBoardMaterial() << "' not found");
172 G4Material* material2 = geometry::Materials::get(geo.getHVBoardMaterial());
173 if (!material2) B2FATAL("Material '" << geo.getHVBoardMaterial() << "' not found");
174 double density1 = material1->GetDensity() / CLHEP::kg * CLHEP::cm * CLHEP::cm * CLHEP::cm;
175 double density2 = material2->GetDensity() / CLHEP::kg * CLHEP::cm * CLHEP::cm * CLHEP::cm;
176
177 PCBarea = S1 + S2; // [cm^2], old value was: 496.725
178 PCBmass = V1 * density1 + V2 * density2; // [kg], old value was: 0.417249
179
180 yearns = 1.e13;
181 evtoJ = 1.60217653 * 1e-10;
182 mtoc = 1.97530864197531;
183 count = 0;
184 count_occ = 0;
185
186 origingamma_x = 0;
187 origingamma_y = 0;
188 origingamma_z = 0;
189 originpe_x = 0;
190 originpe_y = 0;
191 originpe_z = 0;
192 }
193
195 {
196 // Print run number
197 B2INFO("TOPBackground: Processing:");
198
199 }
200
202 {
203 StoreArray<TOPSimHit> topSimhits;
204 StoreArray<MCParticle> mcParticles;
205 StoreArray<TOPDigit> topDigits;
206 StoreArray<TOPBarHit> topTracks;
207
208 int nHits = topDigits.getEntries();
209 for (int i = 0; i < nHits; i++) {
210 TOPDigit* aDigit = topDigits[i];
211 int barID = aDigit->getModuleID();
212
213 peflux->AddBinContent(barID * 2, 1. / m_TimeOfSimulation / 32.0);
214
215 const TOPSimHit* simHit = DataStore::getRelated<TOPSimHit>(aDigit);
216 if (!simHit) continue;
217 int PMTID = simHit->getPmtID();
218
219 module_occupancy->SetPoint(count_occ, PMTID, barID);
220 count_occ++;
221
222 genergy->Fill(simHit->getEnergy());
223
224 const MCParticle* particle = DataStore::getRelated<MCParticle>(simHit);
225
226 if (particle) {
227 const MCParticle* currParticle = particle;
228
229 const MCParticle* mother = currParticle->getMother();
230
231 while (mother) {
232 const MCParticle* pommother = mother->getMother();
233 if (!pommother) {
234
235 zdist->Fill(mother->getVertex().Z());
236 zdistg->Fill(mother->getVertex().Z(), 1. / m_TimeOfSimulation / 32.0 / 16.0);
237 originpe_x = mother->getVertex().X();
238 originpe_y = mother->getVertex().Y();
239 originpe_z = mother->getVertex().Z();
240 originpe->Fill();
241 ROOT::Math::XYZVector momentum = mother->getMomentum();
242 if (m_BkgType.at(m_BkgType.size() - 3) == 'L') ROOT::Math::VectorUtil::RotateY(momentum, 0.0415);
243 else if (m_BkgType.at(m_BkgType.size() - 3) == 'H') ROOT::Math::VectorUtil::RotateY(momentum, -0.0415);
244 double px = momentum.X();
245 double py = momentum.Y();
246 double pt = sqrt(px * px + py * py);
247 originpt->Fill(pt);
248 }
249 mother = pommother;
250 }
251 }
252 }
253
254
255 StoreArray<BeamBackHit> beamBackHits;
256 nHits = beamBackHits.getEntries();
257
258 for (int iHit = 0; iHit < nHits; ++iHit) {
259 BeamBackHit* tophit = beamBackHits[iHit];
260 int subdet = tophit->getSubDet();
261 if (subdet != 5) continue;
262
263 auto pos = tophit->getPosition();
264 double phi = ROOT::Math::VectorUtil::Phi_0_2pi(pos.Phi()) * TMath::RadToDeg();
265 int barID = int (phi / 22.5 + 0.5);
266 if (barID == 16) {
267 barID = 0;
268 }
269 barID++;
270
271
272 if (tophit->getPDG() == Const::neutron.getPDGCode()) {
273 double w = tophit->getNeutronWeight();
274 double tlen = tophit->getTrackLength();
275
276 nflux->AddBinContent(barID * 2, w / m_TimeOfSimulation / PCBarea * yearns * tlen / 0.2);
277
278 } else {
279 double edep = tophit->getEnergyDeposit();
280 rdose->AddBinContent(barID * 2, edep / m_TimeOfSimulation * yearns / PCBmass * evtoJ);
281 }
282 }
283
284 nHits = topTracks.getEntries();
285 for (int iHit = 0; iHit < nHits; ++iHit) {
286 TOPBarHit* toptrk = topTracks[iHit];
287
288 int PDG = toptrk->getPDG();
289 int barID = toptrk->getModuleID();
290
291 if (PDG == Const::neutron.getPDGCode()) {
292 nflux_bar->Fill(toptrk->getPosition().Z(), (barID - 1) * 22.5,
293 1. / 917.65 / m_TimeOfSimulation * yearns * 2.);
294 norigin->Fill(toptrk->getProductionPoint().Z());
295 } else {
296 if (PDG == Const::photon.getPDGCode()) {
297 gflux_bar->Fill(toptrk->getPosition().Z(), (barID - 1) * 22.5,
298 1. / 917.65 / m_TimeOfSimulation * 2.);
299 gorigin->Fill(toptrk->getProductionPoint().Z());
300 genergy2->Fill(toptrk->getMomentum().R() * 1000);
301 origin_zx->SetPoint(count, toptrk->getProductionPoint().Z(),
302 toptrk->getProductionPoint().X());
303 origin_zy->SetPoint(count, toptrk->getProductionPoint().Z() / 0.999143,
304 toptrk->getProductionPoint().Y());
305 origingamma_x = toptrk->getProductionPoint().X();
306 origingamma_y = toptrk->getProductionPoint().Y();
307 origingamma_z = toptrk->getProductionPoint().Z();
308 origingamma->Fill();
309 count++;
310
311 } else {
312 cflux_bar->Fill(toptrk->getPosition().Z(), (barID - 1) * 22.5,
313 1. / 917.65 / m_TimeOfSimulation * 2.);
314 corigin->Fill(toptrk->getProductionPoint().Z());
315 }
316 }
317
318 }
319 }
320
321
322 void TOPBackgroundModule::myprint(TH1F* histo, const char* path, const char* xtit = "", const char* ytit = "", double tresh = 0)
323 {
324
325 gROOT->Reset();
326 gStyle->SetOptStat("");
327 gStyle->SetOptFit(1111);
328
329 gStyle->SetCanvasColor(-1);
330 gStyle->SetPadColor(-1);
331 gStyle->SetFrameFillColor(-1);
332 gStyle->SetHistFillColor(-1);
333 gStyle->SetTitleFillColor(-1);
334 gStyle->SetFillColor(-1);
335 gStyle->SetFillStyle(4000);
336 gStyle->SetStatStyle(0);
337 gStyle->SetTitleStyle(0);
338 gStyle->SetCanvasBorderSize(0);
339 gStyle->SetCanvasBorderMode(0);
340 gStyle->SetPadBorderMode(0);
341 // gStyle->SetTitleMode(0);
342 gStyle->SetFrameBorderSize(0);
343 gStyle->SetLegendBorderSize(0);
344 gStyle->SetStatBorderSize(0);
345 gStyle->SetTitleBorderSize(0);
346 //gROOT->ForceStyle();*/
347
348
349
350 TCanvas* c1 = new TCanvas("c1", "", 1920, 1200);
351
352 double x1 = histo->GetBinLowEdge(1);
353 double nb = histo->GetNbinsX();
354 double bin = histo->GetBinWidth(1);
355 double x2 = x1 + bin * nb;
356
357 double max = histo->GetBinContent(histo->GetMaximumBin());
358
359 if (max < tresh) {
360 histo->GetYaxis()->SetRangeUser(0, tresh * 1.1);
361 }
362
363 TLine* line = new TLine(x1, tresh, x2, tresh);
364 line->SetLineColor(1);
365 line->SetLineWidth(3);
366 line->SetLineStyle(2);
367
368
369 histo->SetFillColor(2);
370 histo->SetLineColor(1);
371
372 gPad->SetTopMargin(0.08);
373 gPad->SetBottomMargin(0.15);
374 gPad->SetGridy();
375
376 histo->GetXaxis()->SetLabelSize(0.06);
377 histo->GetYaxis()->SetLabelSize(0.06);
378 histo->GetXaxis()->SetTitleSize(0.06);
379 histo->GetYaxis()->SetTitleSize(0.06);
380 histo->GetXaxis()->SetTitle(xtit);
381 histo->GetYaxis()->SetTitle(ytit);
382 histo->GetXaxis()->SetTitleOffset(0.9);
383 histo->GetYaxis()->SetTitleOffset(0.7);
384
385 histo->Draw();
386
387 TLegend* leg = new TLegend(0.75, 0.95, 0.90, 1.00);
388 leg->AddEntry(histo, m_BkgType.c_str(), "pf");
389 leg->Draw("SAME");
390 if (tresh > 0.01) {
391 line->Draw("SAME");
392 }
393
394 c1->Print(path);
395 }
396
397
398
400 {
401 B2INFO("TOPBackground: Finished:");
402 }
403
405 {
406 /*
407 myprint(peflux, ("peflux_" + m_BkgType + ".pdf").c_str(), "#phi", "MHz / PMT", 1.);
408 myprint(zdist, ("zdist_" + m_BkgType + ".pdf").c_str(), "z[cm]", "", 0.0);
409 myprint(nflux, ("nflux_" + m_BkgType + ".pdf").c_str(), "#phi", "neutrons / cm^{2} / year", 0.0);
410 myprint(rdose, ("rdose_" + m_BkgType + ".pdf").c_str(), "#phi", "Gy/year", 0.0);
411 */
412
413 m_rootFile->cd();
414 origingamma->Write();
415 originpe->Write();
416 peflux->Write();
417 zdist->Write();
418 nflux->Write();
419 rdose->Write();
420 genergy->Write();
421 genergy2->Write();
422 nflux_bar->Write();
423 gflux_bar->Write();
424 origin_zx->Write();
425 origin_zy->Write();
426 module_occupancy->Write();
427 gorigin->Write();
428 norigin->Write();
429 zdistg->Write();
430 originpt->Write();
431 m_rootFile->Close();
432
433 // Announce
434 B2INFO("TOPBackground finished");
435
436
437 }
438
440 {
441
442 }
443
445} // end Belle2 namespace
Class BeamBackHit - Stores hits from beam backgound simulation.
Definition BeamBackHit.h:28
double getNeutronWeight() const
get the effective neutron weigth
double getEnergyDeposit() const
Get particle energy deposit in sensitive volume.
int getPDG() const
Get the lund code of the particle that hit the sensitive area.
Definition BeamBackHit.h:89
double getTrackLength() const
the length of the track in the volume
ROOT::Math::XYZVector getPosition() const
Get global position of the particle hit.
Definition BeamBackHit.h:95
int getSubDet() const
Det the index of subdetector in which hit occured.
Definition BeamBackHit.h:86
static const ParticleType neutron
neutron particle
Definition Const.h:675
static const ParticleType photon
photon particle
Definition Const.h:673
static T * getRelated(const TObject *object, const std::string &name="", const std::string &namedRelation="")
Get the object to or from which another object has a relation.
Definition DataStore.h:432
A Class to store the Monte Carlo particle information.
Definition MCParticle.h:32
ROOT::Math::XYZVector getVertex() const
Return production vertex position, shorthand for getProductionVertex().
Definition MCParticle.h:172
ROOT::Math::XYZVector getMomentum() const
Return momentum.
Definition MCParticle.h:187
void setDescription(const std::string &description)
Sets the description of the module.
Definition Module.cc:214
Module()
Constructor.
Definition Module.cc:30
Accessor to arrays stored in the data store.
Definition StoreArray.h:113
int getEntries() const
Get the number of objects in the array.
Definition StoreArray.h:216
double origingamma_z
points from origin_zx and _zy graphs: z
TTree * originpe
tree for saving originpe points
TFile * m_rootFile
root file for saving histograms
TH2F * cflux_bar
charged flux on bar
double origingamma_y
points from origin_zx and _zy graphs: y
double originpe_x
points for origin of mother of photoelectrons: x
TH1F * zdistg
z distribution of the photoelectron flux
TH2F * gflux_bar
gamma flux on bar
TH1F * genergy2
energy of gamma that hits the bar
TH2F * nflux_bar
neutron flux on bar
double origingamma_x
points from origin_zx and _zy graphs: x
double originpe_z
points for origin of mother of photoelectrons: z (=zdist)
std::string m_BkgType
Type of background.
std::string m_OutputFileName
Output filename.
double originpe_y
points for origin of mother of photoelectrons: y
TTree * origingamma
tree for saving origingamma points
double m_TimeOfSimulation
Time of the simulation of the input file.
TH1F * zdist
z distribution of primaries
TH1F * genergy
energy of gamma that hits the bar
Class to store track parameters of incoming MC particles relation to MCParticle filled in top/simulat...
Definition TOPBarHit.h:27
ROOT::Math::XYZPoint getPosition() const
Returns impact point.
Definition TOPBarHit.h:104
int getModuleID() const
Returns module ID.
Definition TOPBarHit.h:86
ROOT::Math::XYZPoint getProductionPoint() const
Returns production point.
Definition TOPBarHit.h:98
int getPDG() const
Returns PDG code of particle.
Definition TOPBarHit.h:92
ROOT::Math::XYZVector getMomentum() const
Returns impact momentum.
Definition TOPBarHit.h:116
Class to store TOP digitized hits (output of TOPDigitizer or raw data unpacker) relations to TOPSimHi...
Definition TOPDigit.h:24
int getModuleID() const
Returns module ID.
Definition TOPDigit.h:270
double getFrontBoardWidth() const
Returns front board width.
const TOPGeoFrontEnd & getFrontEnd() const
Returns front-end.
Class to store simulated hits of Cherenkov photons on PMT's input for digitization module (TOPDigitiz...
Definition TOPSimHit.h:27
double getEnergy() const
Returns photon energy.
Definition TOPSimHit.h:108
int getPmtID() const
Returns PMT ID.
Definition TOPSimHit.h:72
const TOPGeometry * getGeometry() const
Returns pointer to geometry object using basf2 units.
static TOPGeometryPar * Instance()
Static method to obtain the pointer to its instance.
static G4Material * get(const std::string &name)
Find given material.
Definition Materials.h:63
void addParam(const std::string &name, T &paramVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
Definition Module.h:559
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition Module.h:649
double sqrt(double a)
sqrt for double
Definition beamHelpers.h:28
MCParticle * getMother() const
Returns a pointer to the mother particle.
Definition MCParticle.h:590
virtual void initialize() override
Initialize the Module.
virtual void event() override
Event processor.
void myprint(TH1F *histo, const char *path, const char *xtit, const char *ytit, double tresh)
Print histogram 1D, helper function.
virtual void endRun() override
End-of-run action.
virtual void terminate() override
Termination action.
virtual ~TOPBackgroundModule()
Destructor.
virtual void beginRun() override
Called when entering a new run.
void printModuleParams() const
Prints module parameters.
Abstract base class for different kinds of events.
STL namespace.