Belle II Software  release-08-01-10
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 
38 using 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"
41 double sourceTime[9] = {0.1, 50., 8., 8., 80., 80., 5., 800., 250.};
42 int sourceType[9] = {0, 1, 2, 3, 4, 5, 6, 7, 1};
43 
44 // returns "ring number" of HAPD module with id modID
45 int 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 
57 TCanvas* 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 
125 void 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:185
Belle2::sqrt
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28