Belle II Software development
BeamBack_arich.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 This macro produces ARICH background plots.
10
11 Run it as: root -l 'BeamBack_arich.cc("path to files")'
12
13 Path to files should be directory containing files produced by "ARICHBkg.py".
14 It is assumed file names contain source of background ...RBB,BHWide...,...
15
16 !! IMPORTANT
17 set the correct time size of each background sample! (in line 42)
18*/
19
20#include <iostream>
21#include <map>
22#include <sstream>
23#include <math.h>
24#include <TH1D.h>
25#include <TH2D.h>
26#include <TGraph2D.h>
27#include <TFile.h>
28#include <TTree.h>
29#include <TChain.h>
30#include <TGaxis.h>
31#include <TBranch.h>
32#include <TObjArray.h>
33#include <TLeaf.h>
34#include <TVector3.h>
35#include <TCanvas.h>
36#include <TPaveText.h>
37
38using namespace std;
39
40// time of simulated sample (total) for all background sources (in ms): "RBB", "BHWide", "Touschek_HER", "Touschek_LER", "Coulomb_HER", "Coulomb_LER","twoPhotons", "brems", "BHWideLargeAngle"
41double sourceTime[9] = {0.1, 50., 8., 8., 80., 80., 5., 800., 250.};
42int sourceType[9] = {0, 1, 2, 3, 4, 5, 6, 7, 1};
43
44// returns "ring number" of HAPD module with id modID
45int mod2row(int modID)
46{
47 if (modID <= 42) return 0;
48 if (modID <= 90) return 1;
49 if (modID <= 144) return 2;
50 if (modID <= 204) return 3;
51 if (modID <= 270) return 4;
52 if (modID <= 342) return 5;
53 if (modID <= 420) return 6;
54 return -1; // -1 if invalid input
55}
56
57TCanvas* make_plot(double data[][8], TString title, int type)
58{
59
60 TH1F* h1 = new TH1F(title + "h1", "RBB", 7, 0.5, 7.5);
61 TH1F* h2 = new TH1F(title + "h2", "BHWide", 7, 0.5, 7.5);
62 TH1F* h3 = new TH1F(title + "h3", "Touschek HER", 7, 0.5, 7.5);
63 TH1F* h4 = new TH1F(title + "h4", "Touschek LER", 7, 0.5, 7.5);
64 TH1F* h5 = new TH1F(title + "h5", "Coulomb HER", 7, 0.5, 7.5);
65 TH1F* h6 = new TH1F(title + "h6", "Coulomb LER", 7, 0.5, 7.5);
66 TH1F* h7 = new TH1F(title + "h7", "2-photon", 7, 0.5, 7.5);
67 TH1F* h8 = new TH1F(title + "h8", "brems", 7, 0.5, 7.5);
68 h1->SetStats(0); h2->SetStats(0); h3->SetStats(0);
69 h4->SetStats(0); h5->SetStats(0); h6->SetStats(0); h7->SetStats(0); h8->SetStats(0);
70 h1->SetBit(TH1::kNoTitle); h2->SetBit(TH1::kNoTitle); h3->SetBit(TH1::kNoTitle); h8->SetBit(TH1::kNoTitle);
71 h4->SetBit(TH1::kNoTitle); h5->SetBit(TH1::kNoTitle); h6->SetBit(TH1::kNoTitle); h7->SetBit(TH1::kNoTitle);
72
73 for (int i = 0; i < 8; i++) {
74 h1->SetBinContent(i + 1, data[i][0]);
75 h2->SetBinContent(i + 1, data[i][0] + data[i][1]);
76 h3->SetBinContent(i + 1, data[i][0] + data[i][1] + data[i][2]);
77 h4->SetBinContent(i + 1, data[i][0] + data[i][1] + data[i][2] + data[i][3]);
78 h5->SetBinContent(i + 1, data[i][0] + data[i][1] + data[i][2] + data[i][3] + data[i][4]);
79 h6->SetBinContent(i + 1, data[i][0] + data[i][1] + data[i][2] + data[i][3] + data[i][4] + data[i][5]);
80 h7->SetBinContent(i + 1, data[i][0] + data[i][1] + data[i][2] + data[i][3] + data[i][4] + data[i][5] + data[i][6]);
81 h8->SetBinContent(i + 1, data[i][0] + data[i][1] + data[i][2] + data[i][3] + data[i][4] + data[i][5] + data[i][6] + data[i][7]);
82 }
83
84 TCanvas* c1 = new TCanvas(title);
85 h8->SetLineColor(8);
86 h8->GetXaxis()->SetTitle("HAPD ring #");
87 TPaveText* pt1 = new TPaveText(0.20, 0.92, 0.8, 0.99, "NDC");
88 if (type == 1) { h8->GetYaxis()->SetTitle("radiation dose [gray/year]"); pt1->AddText("Radiation dose");}
89 if (type == 2) { h8->GetYaxis()->SetTitle("1MeV equiv. neutrons / cm^2 / year"); pt1->AddText("1Mev equiv. neutron flux");}
90 if (type == 3) { h8->GetYaxis()->SetTitle("photons / cm^2 / s"); pt1->AddText("photon flux on HAPDs");}
91
92 h8->SetFillColor(8);
93 h8->Draw();
94 pt1->Draw();
95
96 h7->SetFillColor(7);
97 h7->SetLineColor(7);
98 h7->Draw("same");
99
100 h6->SetFillColor(1);
101 h6->SetLineColor(1);
102 h6->Draw("same");
103
104 h5->SetFillColor(6);
105 h5->SetLineColor(6);
106 h5->Draw("same");
107 h4->SetFillColor(5);
108 h4->SetLineColor(5);
109 h4->Draw("same");
110 h3->SetFillColor(3);
111 h3->SetLineColor(3);
112 h3->Draw("same");
113 h2->SetFillColor(4);
114 h2->SetLineColor(4);
115 h2->Draw("same");
116 h1->SetFillColor(2);
117 h1->SetLineColor(2);
118 h1->Draw("same");
119 c1->BuildLegend(0.6, 0.6, 0.9, 0.9);
120 c1->SetTitle("");
121 return c1;
122}
123
124
125void BeamBack_arich(std::string path = "/gpfs/home/belle/mmanca/basf2/background/16th_campaign/arich_")
126{
127
128 double time = 1000.; // 1 ms
129 gROOT->SetBatch(kTRUE);
130
131 Belle2::ARICHChannelHist* mapsPhotons[9];
132 Belle2::ARICHChannelHist* mapsNeutrons[9];
133 Belle2::ARICHChannelHist* mapsDose[9];
134
135 const TString tHist[9] = {"RBB", "BHWide", "Touschek_HER", "Touschek_LER", "Coulomb_HER", "Coulomb_LER", "twoPhoton", "brems", "total"};
136
137 for (int i = 0; i < 9; i++) {
138 mapsPhotons[i] = new Belle2::ARICHChannelHist("photonFlux_" + tHist[i], "photons / cm^2 / s - " + tHist[i], 1);
139 mapsNeutrons[i] = new Belle2::ARICHChannelHist("neutronFlux_" + tHist[i],
140 "1MeV eq. neutron flux /cm2 / year (x 10^{11}) - " + tHist[i], 1);
141 mapsDose[i] = new Belle2::ARICHChannelHist("radDose_" + tHist[i], "radiation dose [gray/year] - " + tHist[i], 1);
142 }
143
144
145 double edep_board[420][8]; // background contribution from each source
146 double edep_hapd[420][8];
147 double nflux_board[420][8]; // This holds neutron flux on each HAPD module
148 double nflux_hapd[420][8];
149 double flux_phot[420][8];
150 double edep_board_all[420]; // sum of background from all sources
151 double nflux_board_all[420];
152 // double edep_hapd_all[420];
153 // double nflux_hapd_all[420];
154 double phot_all[420];
155
156 // hapd x,y positions
157 double origins[420][2];
158
159 // intialize
160 for (int j = 0; j < 8; j++) {
161 for (int i = 0; i < 420; i++) {
162 edep_board[i][j] = 0;
163 edep_hapd[i][j] = 0;
164 nflux_board[i][j] = 0;
165 nflux_hapd[i][j] = 0;
166 flux_phot[i][j] = 0;
167 origins[i][0] = 0.; origins[i][1] = 0.;
168 edep_board_all[i] = 0.; // edep_hapd_all[i] = 0.;
169 nflux_board_all[i] = 0.; // nflux_hapd_all[i] = 0.;
170 phot_all[i] = 0.;
171 }
172 }
173
174 TChain* tree = new TChain("TrHits");
175
176 // sources
177 //TString tip[7] = {"RBB", "BHWide", "Touschek_HER", "Touschek_LER", "Coulomb_HER", "Coulomb_LER", "2-photon"};
178 // input file names
179 /* for (int j = 0; j < 7; j++) {
180 for (int i = 0; i < 1; i++) {
181 stringstream filename;
182 filename << path.c_str() << "*" << tip[j] << "*" << ".root";
183 printf("%s\n", filename.str().c_str());
184 tree->Add(filename.str().c_str());
185 }
186 }*/
187
188 stringstream filename;
189 filename << path.c_str() << "*.root";
190 tree->Add(filename.str().c_str());
191
192 // this is neutron energy spectrum
193 TH1D* nenergy = new TH1D("nenergy", "neutron spectrum; log10(k.e./MeV)", 2000, -10., 2.);
194 nenergy->SetStats(0);
195 // set variables from input tree files
196 int type = -1;
197 int modID;
198 int pdg = -1;
199 int source;
200 double edep = 0;
201 TVector3* modOrig = 0;
202 TVector3* mom = 0;
203 double phnw = 0; double trlen = 0; double en = 0;
204 tree->SetBranchAddress("phPDG", &pdg);
205 tree->SetBranchAddress("moduleID", &modID);
206 tree->SetBranchAddress("type", &type);
207 tree->SetBranchAddress("edep", &edep);
208 tree->SetBranchAddress("modOrig", &modOrig);
209 tree->SetBranchAddress("phmom", &mom);
210 tree->SetBranchAddress("phnw", &phnw);
211 tree->SetBranchAddress("trlen", &trlen);
212 tree->SetBranchAddress("en", &en);
213 tree->SetBranchAddress("source", &source);
214
215 int nevents = tree->GetEntries();
216
217 // loop over all hits
218 for (int e = 0; e < nevents; e++) {
219
220 if (e % 1000000 == 0) std::cout << "Processed: " << double(e) / double(nevents) << std::endl;
221 edep = 0;
222 type = -1;
223 pdg = -1;
224
225 mom->SetXYZ(0., 0., 0.);
226
227 tree->GetEvent(e);
228
229 if (source < 0) continue;
230
231 // read the corresponding module (HAPD) x,y coordinates
232 if (modID > 0) {
233 origins[modID - 1][0] = modOrig->X();
234 origins[modID - 1][1] = modOrig->Y();
235 } else continue;
236
237 if (sourceTime[source] == 0) continue;
238
239 // this is for photon flux
240 if (type == 2) {
241 flux_phot[modID - 1][sourceType[source]] += 1. / sourceTime[source];
242 }
243
244 if (type != 2) {
245
246 // energy deposit from a hit (two volumes are sensitive el. board and hapd bottom, commonly I show results for board.
247 // Anyhow the rates are approximately equal.)
248 if (pdg != Const::neutron.getPDGCode()) {
249 edep /= sourceTime[source];
250 if (type == 0) edep_board[modID - 1][sourceType[source]] += edep;
251 if (type == 1) edep_hapd[modID - 1][sourceType[source]] += edep;
252 }
253
254 else { // if neutron
255
256 // Calculate neutron kinetic energy in MeV. mom is obtained from BeamBackHit->GetMomentum() and is in GeV
257 double mass = 940.0;
258 double pp = mom->Mag() * 1000.;
259 double ee = sqrt(pp * pp + mass * mass) - mass;
260 double lee = log10(ee);
261
262 // contibution to neutron flux
263 if (type == 0) {
264 phnw /= sourceTime[source];
265 nflux_board[modID - 1][sourceType[source]] += phnw * trlen / 0.2; // trlen/0.2, this is kind of cos(theta) correction to the flux.
266 // trlen = track length in volume, and 0.2cm is thickness of board.
267 nenergy->Fill(lee);
268 }
269 if (type == 1) {
270 phnw *= sourceTime[source];
271 nflux_hapd[modID - 1][sourceType[source]] += phnw * trlen / 0.05;
272 }
273 }
274 }
275
276 } // end of event loop
277
278
279 // values averaged over HAPD rings;
280 double avgeb[7][8];
281 double avgeh[7][8];
282 double avgnb[7][8];
283 double avgnh[7][8];
284 double avgf[7][8];
285
286 for (int i = 0; i < 7; i++) {
287 for (int j = 0; j < 8; j++) {
288 avgeb[i][j] = 0; avgeh[i][j] = 0;
289 avgnb[i][j] = 0; avgnh[i][j] = 0;
290 avgf[i][j] = 0;
291 }
292 }
293
294 const double nhapd[7] = {42., 48., 54., 60., 66., 72., 78.}; // 7 rings
295
296 // This is important. nflux_board holds total flux of neutrons wighted by 1MeV equiv. factor in a given
297 // time "time". Here this is transformed in flux/cm^2/year.
298
299 TH1F* ndist = new TH1F("ndist", "nflux distribution;1MeV equiv. neutrons / cm^2 / year; # of HAPDs", 1000, 0, 8);
300 TH1F* edist = new TH1F("edist", "deposit en. distribution; radiation dose [gray/year]; # of HAPDs", 1000, 0, 10);
301
302 for (int i = 0; i < 420; i++) {
303 int nrow = mod2row(i + 1);
304 for (int j = 0; j < 8; j++) {
305
306 edep_board[i][j] = edep_board[i][j] * 1.6E+6 / 47.25 / time; // 47.25 mass of FEB in grams
307 edep_hapd[i][j] = edep_hapd[i][j] * 1.6E+6 / 5.49 / time; // 5.49 mass of APD chip (volume 62.7x62.7x0.06)
308
309 nflux_board[i][j] = nflux_board[i][j] * 1.E+13 / 56.25 / time; // board surface is 56.25 cm^2
310 nflux_hapd[i][j] = nflux_hapd[i][j] * 1.E+13 / 39.3 / time; // apd surface
311
312 flux_phot[i][j] = flux_phot[i][j] * 1.E+6 / 53.29 / time; // photocathode size = 7.3x7.3cm
313
314 edep_board_all[i] += edep_board[i][j];
315 // edep_hapd_all[i] += edep_hapd[i][j];
316 nflux_board_all[i] += nflux_board[i][j];
317 // nflux_hapd_all[i] += nflux_hapd[i][j];
318 phot_all[i] += flux_phot[i][j];
319
320 avgeb[nrow][j] += edep_board[i][j] / nhapd[nrow];
321 avgeh[nrow][j] += edep_hapd[i][j] / nhapd[nrow];
322 avgnb[nrow][j] += nflux_board[i][j] / nhapd[nrow];
323 avgnh[nrow][j] += nflux_hapd[i][j] / nhapd[nrow];
324 avgf[nrow][j] += flux_phot[i][j] / nhapd[nrow];
325
326 mapsPhotons[j]->setBinContent(i + 1, flux_phot[i][j]);
327 mapsNeutrons[j]->setBinContent(i + 1, nflux_board[i][j] / 1E+11);
328 mapsDose[j]->setBinContent(i + 1, edep_board[i][j]);
329
330 }
331
332 double nn = nflux_board_all[i] / 1E+11;
333 ndist->Fill(nn);
334 edist->Fill(edep_board_all[i]);
335
336 // fill 2D graphs
337 mapsPhotons[8]->setBinContent(i + 1, phot_all[i]);
338 mapsNeutrons[8]->setBinContent(i + 1, nn);
339 mapsDose[8]->setBinContent(i + 1, edep_board_all[i]);
340
341 } // end of HAPD loop
342
343
344 TGaxis::SetMaxDigits(1);
345 TCanvas* c1 = make_plot(avgnb, "n_board", 2);
346
347 TCanvas* c2 = make_plot(avgeb, "e_board", 1);
348
349 TCanvas* c3 = make_plot(avgnh, "n_hapd", 2);
350
351 TCanvas* c4 = make_plot(avgeh, "e_hapd", 1);
352
353 TCanvas* c5 = make_plot(avgf, "photons", 3);
354
355 // write output file
356 TFile* f = new TFile("arich_background_phase3.root", "RECREATE");
357
358 c1->Write(); c2->Write(); c3->Write(); c4->Write(); c5->Write();
359 ndist->Write();
360 edist->Write();
361 nenergy->Write();
362
363 for (int i = 0; i < 8; i++) {
364 mapsPhotons[i]->SetMinimum(0.0);
365 mapsNeutrons[i]->SetMinimum(0.0);
366 mapsDose[i]->SetMinimum(0.0);
367 mapsPhotons[i]->Write();
368 mapsNeutrons[i]->Write();
369 mapsDose[i]->Write();
370 }
371
372 f->Close();
373
374 delete c1;
375 delete c2;
376 delete c3;
377 delete c4;
378 delete c5;
379
380}
E
R E
internal precision of FFTW codelets
Definition: DiscreteCosineTransform_31points.cc:36
Belle2::ARICHChannelHist
ARICH histogram with HAPD plane 3 options for bin segmentation are available type 0 - one bin per HAP...
Definition: ARICHChannelHist.h:27
Belle2::ARICHChannelHist::setBinContent
void setBinContent(unsigned hapdID, unsigned chID, double value)
Set content of bin corresponding to hapd hapdID and channel chID.
Definition: ARICHChannelHist.cc:189
Belle2::sqrt
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28
std
STL namespace.