Belle II Software  release-05-01-25
BeamBack_arich.cc
1 /*
2  This macro produces ARICH background plots.
3 
4  Run it as: root -l 'BeamBack_arich.cc(simulated time in us,"path to files")'
5 
6  Path to files should be directory containing files produced by "ARICHBkg.py".
7  It is assumed file names contain source of background ...RBB,BHWide...,...
8 
9  !! IMPORTANT
10  set the correct scalings in case of using samples of different time lenghts (in line 36)
11 
12  Author: Luka Santelj
13  Contributor: Manca Mrvar, Leonid Burmistrov
14 */
15 
16 #include <iostream>
17 #include <map>
18 #include <sstream>
19 #include <math.h>
20 #include <TH1D.h>
21 #include <TH2D.h>
22 #include <TGraph2D.h>
23 #include <TFile.h>
24 #include <TTree.h>
25 #include <TChain.h>
26 #include <TGaxis.h>
27 #include <TBranch.h>
28 #include <TObjArray.h>
29 #include <TLeaf.h>
30 #include <TVector3.h>
31 #include <TCanvas.h>
32 #include <TPaveText.h>
33 
34 using namespace std;
35 
36 // scalings for "RBB", "BHWide", "Touschek_HER", "Touschek_LER", "Coulomb_HER", "Coulomb_LER","2-photon", "BHWideLargeAngle"
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};
39 
40 // returns "ring number" of HAPD module with id modID
41 int mod2row(int modID)
42 {
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;
50  return -1; // -1 if invalid input
51 }
52 
53 TCanvas* make_plot(double data[][7], TString title, int type)
54 {
55 
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);
67 
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]);
76  }
77  TCanvas* c1 = new TCanvas(title);
78  h7->SetLineColor(7);
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");}
84 
85  h7->SetFillColor(7);
86  h7->Draw();
87  pt1->Draw();
88 
89  h6->SetFillColor(1);
90  h6->SetLineColor(1);
91  h6->Draw("same");
92 
93  h5->SetFillColor(6);
94  h5->SetLineColor(6);
95  h5->Draw("same");
96  h4->SetFillColor(5);
97  h4->SetLineColor(5);
98  h4->Draw("same");
99  h3->SetFillColor(3);
100  h3->SetLineColor(3);
101  h3->Draw("same");
102  h2->SetFillColor(4);
103  h2->SetLineColor(4);
104  h2->Draw("same");
105  h1->SetFillColor(2);
106  h1->SetLineColor(2);
107  h1->Draw("same");
108  c1->BuildLegend(0.6, 0.6, 0.9, 0.9);
109  c1->SetTitle("");
110  return c1;
111 }
112 
113 
114 //int BeamBack_arich(double time = 1000, std::string path = "/gpfs/home/belle/mmanca/basf2/background/16th_campaign/arich_")
115 void BeamBack_arich(double time = 1000, std::string path = "/gpfs/home/belle/mmanca/basf2/background/16th_campaign/arich_")
116 {
117 
118  gROOT->SetBatch(kTRUE);
119  TGraph2D* dteb = new TGraph2D(420);
120  TGraph2D* dtnb = new TGraph2D(420);
121  Belle2::ARICHChannelHist* chHits = new Belle2::ARICHChannelHist("neutronFlux", "1MeV eq. neutron flux /cm2 / year on HAPD board",
122  1);
123 
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]");
132 
133  double edep_board[420][7]; // background contribution from each source
134  double edep_hapd[420][7];
135  double nflux_board[420][7]; // This holds neutron flux on each HAPD module
136  double nflux_hapd[420][7];
137  double flux_phot[420][7];
138  double edep_board_all[420]; // sum of background from all sources
139  double nflux_board_all[420];
140  double edep_hapd_all[420];
141  double nflux_hapd_all[420];
142  double phot_all[420];
143 
144  // hapd x,y positions
145  double origins[420][2];
146 
147  // intialize
148  for (int j = 0; j < 7; j++) {
149  for (int i = 0; i < 420; i++) {
150  edep_board[i][j] = 0;
151  edep_hapd[i][j] = 0;
152  nflux_board[i][j] = 0;
153  nflux_hapd[i][j] = 0;
154  flux_phot[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.;
158  phot_all[i] = 0.;
159  }
160  }
161 
162  TChain* tree = new TChain("TrHits");
163 
164  // sources
165  //TString tip[7] = {"RBB", "BHWide", "Touschek_HER", "Touschek_LER", "Coulomb_HER", "Coulomb_LER", "2-photon"};
166 
167 
168  // input file names
169  /* for (int j = 0; j < 7; j++) {
170  for (int i = 0; i < 1; i++) {
171  stringstream filename;
172  filename << path.c_str() << "*" << tip[j] << "*" << ".root";
173  printf("%s\n", filename.str().c_str());
174  tree->Add(filename.str().c_str());
175  }
176  }*/
177 
178  stringstream filename;
179  filename << path.c_str() << "*.root";
180  tree->Add(filename.str().c_str());
181 
182  // this is neutron energy spectrum
183  TH1D* nenergy = new TH1D("nenergy", "neutron spectrum; log10(k.e./MeV)", 2000, -10., 2.);
184  nenergy->SetStats(0);
185  // set variables from input tree files
186  int type = -1;
187  int modID;
188  int pdg = -1;
189  int source;
190  double edep = 0;
191  TVector3* modOrig = 0;
192  TVector3* mom = 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);
204 
205  int nevents = tree->GetEntries();
206  //double nmax = 0;
207  // loop over all hits
208  for (int e = 0; e < nevents; e++) {
209 
210  edep = 0;
211  type = -1;
212  pdg = -1;
213  //modID = 1;
214  mom->SetXYZ(0., 0., 0.);
215 
216  tree->GetEvent(e);
217 
218  if (source < 0) continue;
219  if (modID < 1) continue;
220 
221  // read the corresponding module (HAPD) x,y coordinates
222  if (modID > 0) {
223  origins[modID - 1][0] = modOrig->X();
224  origins[modID - 1][1] = modOrig->Y();
225  } else continue;
226 
227  // this is for photon flux
228  if (type == 2) {
229  flux_phot[modID - 1][sourceType[source]] += scaling[source];
230  }
231 
232  if (type != 2) {
233 
234  // energy deposit from a hit (two volumes are sensitive el. board and hapd bottom, commonly I show results for board.
235  // Anyhow the rates are approximately equal.)
236  if (pdg != 2112) {
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;
240  }
241 
242  else { // if neutron
243 
244  // Calculate neutron kinetic energy in MeV. mom is obtained from BeamBackHit->GetMomentum() and is in GeV
245  double mass = 940.0;
246  double pp = mom->Mag() * 1000.;
247  double ee = sqrt(pp * pp + mass * mass) - mass;
248  double lee = log10(ee);
249 
250  // contibution to neutron flux
251  if (type == 0) {
252  if (modID == 10) continue;
253  phnw *= scaling[source];
254  nflux_board[modID - 1][sourceType[source]] += phnw * trlen / 0.2; // trlen/0.2, this is kind of cos(theta) correction to the flux.
255  // trlen = track length in volume, and 0.2cm is thickness of board.
256  nenergy->Fill(lee);
257  }
258  if (type == 1) {
259  phnw *= scaling[source];
260  nflux_hapd[modID - 1][sourceType[source]] += phnw * trlen / 0.05;
261  }
262  }
263  }
264 
265  } // end of event loop
266 
267 
268  // values averaged over HAPD rings;
269  double avgeb[7][7];
270  double avgeh[7][7];
271  double avgnb[7][7];
272  double avgnh[7][7];
273  double avgf[7][7];
274 
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;
279  avgf[i][j] = 0;
280  }
281  }
282 
283  double nhapd[7] = {42., 48., 54., 60., 66., 72., 78.}; // 7 rings
284 
285  // This is important. nflux_board holds total flux of neutrons wighted by 1MeV equiv. factor in a given
286  // time "time". Here this is transformed in flux/cm^2/year.
287 
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);
290 
291  for (int i = 0; i < 420; i++) {
292  int nrow = mod2row(i + 1);
293  for (int j = 0; j < 7; j++) {
294 
295  edep_board[i][j] = edep_board[i][j] * 1.6E+6 / 47.25 / time; // 47.25 mass of FEB in grams
296  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)
297 
298  nflux_board[i][j] = nflux_board[i][j] * 1.E+13 / 56.25 / time; // board surface is 56.25 cm^2
299  nflux_hapd[i][j] = nflux_hapd[i][j] * 1.E+13 / 39.3 / time; // apd surface
300 
301  flux_phot[i][j] = flux_phot[i][j] * 1.E+6 / 53.29 / time; // photocathode size = 7.3x7.3cm
302 
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];
308 
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];
314  }
315 
316  double nn = nflux_board_all[i] / 1E+11;
317  ndist->Fill(nn);
318  edist->Fill(edep_board_all[i]);
319  // fill 2D graphs
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]);
325  chHits->setBinContent(i + 1, nn);
326 
327  } // end of HAPD loop
328 
329  dteb->SetMarkerStyle(21);
330  dteh->SetMarkerStyle(21);
331  dtnb->SetMarkerStyle(21);
332  dteh->SetMarkerStyle(21);
333  dtp->SetMarkerStyle(21);
334  //chHits->SetMarkerStyle(21);
335 
336 
337  TGaxis::SetMaxDigits(1);
338  TCanvas* c1 = make_plot(avgnb, "n_board", 2);
339 
340  TCanvas* c2 = make_plot(avgeb, "e_board", 1);
341 
342  TCanvas* c3 = make_plot(avgnh, "n_hapd", 2);
343 
344  TCanvas* c4 = make_plot(avgeh, "e_hapd", 1);
345 
346  TCanvas* c5 = make_plot(avgf, "photons", 3);
347 
348  // write output file
349  TFile* f = new TFile("arich_background_phase3.root", "RECREATE");
350 
351  c1->Write(); c2->Write(); c3->Write(); c4->Write(); c5->Write();
352  dteb->Write(); dteh->Write();
353  dtnb->Write(); dtnh->Write();
354  dtp->Write();
355  chHits->Write();
356  ndist->Write();
357  edist->Write();
358  nenergy->Write();
359  f->Close();
360 
361  delete c1;
362  delete c2;
363  delete c3;
364  delete c4;
365  delete c5;
366 
367  //return 0;
368 }
Belle2::ARICHChannelHist
ARICH histogram with HAPD plane 3 options for bin segmentation are available type 0 - one bin per HAP...
Definition: ARICHChannelHist.h:37
prepareAsicCrosstalkSimDB.e
e
aux.
Definition: prepareAsicCrosstalkSimDB.py:53
prepareAsicCrosstalkSimDB.ee
ee
aux.
Definition: prepareAsicCrosstalkSimDB.py:55
alignment.constraints_generator.filename
filename
File name.
Definition: constraints_generator.py:224
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:187