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