28 #include <TObjArray.h>
32 #include <TPaveText.h>
37 double scaling[8] = {100., 1, 10., 10., 10., 10., 10., 0.001};
38 int sourceType[8] = {0, 1, 2, 3, 4, 5, 6, 1};
41 int mod2row(
int modID)
43 if (modID <= 42)
return 0;
44 if (modID <= 90)
return 1;
45 if (modID <= 144)
return 2;
46 if (modID <= 204)
return 3;
47 if (modID <= 270)
return 4;
48 if (modID <= 342)
return 5;
49 if (modID <= 420)
return 6;
53 TCanvas* make_plot(
double data[][7], TString title,
int type)
56 TH1F* h1 =
new TH1F(title +
"h1",
"RBB", 7, 0.5, 7.5);
57 TH1F* h2 =
new TH1F(title +
"h2",
"BHWide", 7, 0.5, 7.5);
58 TH1F* h3 =
new TH1F(title +
"h3",
"Touschek HER", 7, 0.5, 7.5);
59 TH1F* h4 =
new TH1F(title +
"h4",
"Touschek LER", 7, 0.5, 7.5);
60 TH1F* h5 =
new TH1F(title +
"h5",
"Coulomb HER", 7, 0.5, 7.5);
61 TH1F* h6 =
new TH1F(title +
"h6",
"Coulomb LER", 7, 0.5, 7.5);
62 TH1F* h7 =
new TH1F(title +
"h7",
"2-photon", 7, 0.5, 7.5);
63 h1->SetStats(0); h2->SetStats(0); h3->SetStats(0);
64 h4->SetStats(0); h5->SetStats(0); h6->SetStats(0); h7->SetStats(0);
65 h1->SetBit(TH1::kNoTitle); h2->SetBit(TH1::kNoTitle); h3->SetBit(TH1::kNoTitle);
66 h4->SetBit(TH1::kNoTitle); h5->SetBit(TH1::kNoTitle); h6->SetBit(TH1::kNoTitle); h7->SetBit(TH1::kNoTitle);
68 for (
int i = 0; i < 7; i++) {
69 h1->SetBinContent(i + 1, data[i][0]);
70 h2->SetBinContent(i + 1, data[i][0] + data[i][1]);
71 h3->SetBinContent(i + 1, data[i][0] + data[i][1] + data[i][2]);
72 h4->SetBinContent(i + 1, data[i][0] + data[i][1] + data[i][2] + data[i][3]);
73 h5->SetBinContent(i + 1, data[i][0] + data[i][1] + data[i][2] + data[i][3] + data[i][4]);
74 h6->SetBinContent(i + 1, data[i][0] + data[i][1] + data[i][2] + data[i][3] + data[i][4] + data[i][5]);
75 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]);
77 TCanvas* c1 =
new TCanvas(title);
79 h7->GetXaxis()->SetTitle(
"HAPD ring #");
80 TPaveText* pt1 =
new TPaveText(0.20, 0.92, 0.8, 0.99,
"NDC");
81 if (type == 1) { h7->GetYaxis()->SetTitle(
"radiation dose [gray/year]"); pt1->AddText(
"Radiation dose");}
82 if (type == 2) { h7->GetYaxis()->SetTitle(
"1MeV equiv. neutrons / cm^2 / year"); pt1->AddText(
"1Mev equiv. neutron flux");}
83 if (type == 3) { h7->GetYaxis()->SetTitle(
"photons / cm^2 / s"); pt1->AddText(
"photon flux on HAPDs");}
108 c1->BuildLegend(0.6, 0.6, 0.9, 0.9);
115 void BeamBack_arich(
double time = 1000, std::string path =
"/gpfs/home/belle/mmanca/basf2/background/16th_campaign/arich_")
118 gROOT->SetBatch(kTRUE);
119 TGraph2D* dteb =
new TGraph2D(420);
120 TGraph2D* dtnb =
new TGraph2D(420);
124 TGraph2D* dteh =
new TGraph2D(420);
125 TGraph2D* dtnh =
new TGraph2D(420);
126 TGraph2D* dtp =
new TGraph2D(420);
127 dteb->SetName(
"dteb"); dteb->SetTitle(
"Deposited energy [gray/year] in HAPD board");
128 dteh->SetName(
"dteh"); dteh->SetTitle(
"Deposited energy [gray/year] in HAPD bottom");
129 dtnb->SetName(
"dtnb"); dtnb->SetTitle(
"1MeV equiv. neutron flux / cm2 / year on HAPD board [x10^{11}]");
130 dtnh->SetName(
"dtnh"); dtnh->SetTitle(
"1MeV equiv. neutron flux / cm2 / year on HAPD bottom [x10^{11}]");
131 dtp->SetName(
"dtp"); dtp->SetTitle(
"Optical photon flux on hapd [photons / s]");
133 double edep_board[420][7];
134 double edep_hapd[420][7];
135 double nflux_board[420][7];
136 double nflux_hapd[420][7];
137 double flux_phot[420][7];
138 double edep_board_all[420];
139 double nflux_board_all[420];
140 double edep_hapd_all[420];
141 double nflux_hapd_all[420];
142 double phot_all[420];
145 double origins[420][2];
148 for (
int j = 0; j < 7; j++) {
149 for (
int i = 0; i < 420; i++) {
150 edep_board[i][j] = 0;
152 nflux_board[i][j] = 0;
153 nflux_hapd[i][j] = 0;
155 origins[i][0] = 0.; origins[i][1] = 0.;
156 edep_board_all[i] = 0.; edep_hapd_all[i] = 0.;
157 nflux_board_all[i] = 0.; nflux_hapd_all[i] = 0.;
162 TChain* tree =
new TChain(
"TrHits");
179 filename << path.c_str() <<
"*.root";
183 TH1D* nenergy =
new TH1D(
"nenergy",
"neutron spectrum; log10(k.e./MeV)", 2000, -10., 2.);
184 nenergy->SetStats(0);
191 TVector3* modOrig = 0;
193 double phnw = 0;
double trlen = 0;
double en = 0;
194 tree->SetBranchAddress(
"phPDG", &pdg);
195 tree->SetBranchAddress(
"moduleID", &modID);
196 tree->SetBranchAddress(
"type", &type);
197 tree->SetBranchAddress(
"edep", &edep);
198 tree->SetBranchAddress(
"modOrig", &modOrig);
199 tree->SetBranchAddress(
"phmom", &mom);
200 tree->SetBranchAddress(
"phnw", &phnw);
201 tree->SetBranchAddress(
"trlen", &trlen);
202 tree->SetBranchAddress(
"en", &en);
203 tree->SetBranchAddress(
"source", &source);
205 int nevents = tree->GetEntries();
208 for (
int e = 0;
e < nevents;
e++) {
214 mom->SetXYZ(0., 0., 0.);
218 if (source < 0)
continue;
219 if (modID < 1)
continue;
223 origins[modID - 1][0] = modOrig->X();
224 origins[modID - 1][1] = modOrig->Y();
229 flux_phot[modID - 1][sourceType[source]] += scaling[source];
237 edep *= scaling[source];
238 if (type == 0) edep_board[modID - 1][sourceType[source]] += edep;
239 if (type == 1) edep_hapd[modID - 1][sourceType[source]] += edep;
246 double pp = mom->Mag() * 1000.;
247 double ee = sqrt(pp * pp + mass * mass) - mass;
248 double lee = log10(ee);
252 if (modID == 10)
continue;
253 phnw *= scaling[source];
254 nflux_board[modID - 1][sourceType[source]] += phnw * trlen / 0.2;
259 phnw *= scaling[source];
260 nflux_hapd[modID - 1][sourceType[source]] += phnw * trlen / 0.05;
275 for (
int i = 0; i < 7; i++) {
276 for (
int j = 0; j < 7; j++) {
277 avgeb[i][j] = 0; avgeh[i][j] = 0;
278 avgnb[i][j] = 0; avgnh[i][j] = 0;
283 double nhapd[7] = {42., 48., 54., 60., 66., 72., 78.};
288 TH1F* ndist =
new TH1F(
"ndist",
"nflux distribution;1MeV equiv. neutrons / cm^2 / year; # of HAPDs", 1000, 0, 8);
289 TH1F* edist =
new TH1F(
"edist",
"deposit en. distribution; radiation dose [gray/year]; # of HAPDs", 1000, 0, 10);
291 for (
int i = 0; i < 420; i++) {
292 int nrow = mod2row(i + 1);
293 for (
int j = 0; j < 7; j++) {
295 edep_board[i][j] = edep_board[i][j] * 1.6E+6 / 47.25 / time;
296 edep_hapd[i][j] = edep_hapd[i][j] * 1.6E+6 / 5.49 / time;
298 nflux_board[i][j] = nflux_board[i][j] * 1.E+13 / 56.25 / time;
299 nflux_hapd[i][j] = nflux_hapd[i][j] * 1.E+13 / 39.3 / time;
301 flux_phot[i][j] = flux_phot[i][j] * 1.E+6 / 53.29 / time;
303 edep_board_all[i] += edep_board[i][j];
304 edep_hapd_all[i] += edep_hapd[i][j];
305 nflux_board_all[i] += nflux_board[i][j];
306 nflux_hapd_all[i] += nflux_hapd[i][j];
307 phot_all[i] += flux_phot[i][j];
309 avgeb[nrow][j] += edep_board[i][j] / nhapd[nrow];
310 avgeh[nrow][j] += edep_hapd[i][j] / nhapd[nrow];
311 avgnb[nrow][j] += nflux_board[i][j] / nhapd[nrow];
312 avgnh[nrow][j] += nflux_hapd[i][j] / nhapd[nrow];
313 avgf[nrow][j] += flux_phot[i][j] / nhapd[nrow];
316 double nn = nflux_board_all[i] / 1E+11;
318 edist->Fill(edep_board_all[i]);
320 dteb->SetPoint(i, origins[i][0], origins[i][1], edep_board_all[i]);
321 dtnb->SetPoint(i, origins[i][0], origins[i][1], nn);
322 dteh->SetPoint(i, origins[i][0], origins[i][1], edep_hapd_all[i]);
323 dtnh->SetPoint(i, origins[i][0], origins[i][1], nflux_hapd_all[i]);
324 dtp->SetPoint(i, origins[i][0], origins[i][1], phot_all[i]);
329 dteb->SetMarkerStyle(21);
330 dteh->SetMarkerStyle(21);
331 dtnb->SetMarkerStyle(21);
332 dteh->SetMarkerStyle(21);
333 dtp->SetMarkerStyle(21);
337 TGaxis::SetMaxDigits(1);
338 TCanvas* c1 = make_plot(avgnb,
"n_board", 2);
340 TCanvas* c2 = make_plot(avgeb,
"e_board", 1);
342 TCanvas* c3 = make_plot(avgnh,
"n_hapd", 2);
344 TCanvas* c4 = make_plot(avgeh,
"e_hapd", 1);
346 TCanvas* c5 = make_plot(avgf,
"photons", 3);
349 TFile* f =
new TFile(
"arich_background_phase3.root",
"RECREATE");
351 c1->Write(); c2->Write(); c3->Write(); c4->Write(); c5->Write();
352 dteb->Write(); dteh->Write();
353 dtnb->Write(); dtnh->Write();