Belle II Software  release-05-01-25
DQMHistAnalysisPXDFits.cc
1 //+
2 // File : DQMHistAnalysisPXDFits.cc
3 // Description : DQM module, which fits many PXD histograms and writes out fit parameters in new histograms
4 //
5 // Author : Bjoern Spruck, University Mainz
6 // Date : 5 - May - 2017
7 //-
8 
9 
10 #include <dqm/analysis/modules/DQMHistAnalysisPXDFits.h>
11 #include <TROOT.h>
12 
13 #include <boost/format.hpp>
14 
15 using namespace std;
16 using namespace Belle2;
17 
18 using boost::format;
19 
20 //-----------------------------------------------------------------
21 // Register the Module
22 //-----------------------------------------------------------------
23 REG_MODULE(DQMHistAnalysisPXDFits)
24 
25 //-----------------------------------------------------------------
26 // Implementation
27 //-----------------------------------------------------------------
28 
31 {
32  //Parameter definition
33  addParam("histogramDirectoryName", m_histogramDirectoryName, "Name of the directory where histograms were placed",
34  std::string("PXDRAW"));
35 // addParam("HistoName", m_histoname, "Name of Histogram (incl dir)", std::string(""));
36  for (auto i = 0, j = 0; i < 64; i++) {
37  auto layer = (((i >> 5) & 0x1) + 1);
38  auto ladder = ((i >> 1) & 0xF);
39  if (ladder == 0) continue; // numbering starts at 1
40  if (layer == 1 && ladder > 8) continue; // 8 inner ladders
41  if (layer == 2 && ladder > 12) continue; // 12 outer ladders
42  m_id_to_inx[i] = j; // i = id , j - index
43  m_inx_to_id[j] = i;
44  j++;
45  if (j == NUM_MODULES) break;
46  }
47  for (auto i = 0; i < NUM_MODULES; i++) {
48  m_hSignal[i] = nullptr;
49  m_hCommon[i] = nullptr;
50  m_hCounts[i] = nullptr;
51  m_cSignal[i] = nullptr;
52  m_cCommon[i] = nullptr;
53  m_cCounts[i] = nullptr;
54  }
55  B2DEBUG(1, "DQMHistAnalysisPXDFits: Constructor done.");
56 }
57 
58 void DQMHistAnalysisPXDFitsModule::initialize()
59 {
60  B2DEBUG(1, "DQMHistAnalysisPXDFits: initialized.");
61 
62  gROOT->cd(); // this seems to be important, or strange things happen
63 
64  TString a;
65 
66  a = "pxdraw/hSignalAll";
67  a.ReplaceAll("/", "_");
68  m_cSignalAll = new TCanvas((m_histogramDirectoryName + "/c_").data() + a);
69  m_hSignalAll = new TH1F(a, a, NUM_MODULES, 0, NUM_MODULES);
70  m_hSignalAll->SetDirectory(0);// dont mess with it, this is MY histogram
71  m_hSignalAll->SetStats(false);
72 
73  a = "pxdraw/hCommonAll";
74  a.ReplaceAll("/", "_");
75  m_cCommonAll = new TCanvas((m_histogramDirectoryName + "/c_").data() + a);
76  m_hCommonAll = new TH1F(a, a, NUM_MODULES, 0, NUM_MODULES);
77  m_hCommonAll->SetDirectory(0);// dont mess with it, this is MY histogram
78  m_hCommonAll->SetStats(false);
79 
80  a = "pxdraw/hCountsAll";
81  a.ReplaceAll("/", "_");
82  m_cCountsAll = new TCanvas((m_histogramDirectoryName + "/c_").data() + a);
83  m_hCountsAll = new TH1F(a, a, NUM_MODULES, 0, NUM_MODULES);
84  m_hCountsAll->SetDirectory(0);// dont mess with it, this is MY histogram
85  m_hCountsAll->SetStats(false);
86 
87  a = "pxdraw/hOccupancyAll";
88  a.ReplaceAll("/", "_");
89  m_cOccupancyAll = new TCanvas((m_histogramDirectoryName + "/c_").data() + a);
90  m_hOccupancyAll = new TH1F(a, a, NUM_MODULES, 0, NUM_MODULES);
91  m_hOccupancyAll->SetDirectory(0);// dont mess with it, this is MY histogram
92  m_hOccupancyAll->SetStats(false);
93 
94  for (auto i = 0; i < NUM_MODULES; i++) {
95  auto id = m_inx_to_id[i];
96  auto layer = (((id >> 5) & 0x1) + 1);
97  auto ladder = ((id >> 1) & 0xF);
98  auto sensor = ((id & 0x1) + 1);
99  string s2 = str(format("_%d.%d.%d") % layer % ladder % sensor);
100 
101  m_hSignalAll->GetXaxis()->SetBinLabel(i + 1, TString(s2));
102  m_hCommonAll->GetXaxis()->SetBinLabel(i + 1, TString(s2));
103  m_hCountsAll->GetXaxis()->SetBinLabel(i + 1, TString(s2));
104  m_hOccupancyAll->GetXaxis()->SetBinLabel(i + 1, TString(s2));
105 
106  a = "pxdraw/hSignal";
107  a.ReplaceAll("/", "_");
108  a += s2;
109 
110  m_cSignal[i] = new TCanvas((m_histogramDirectoryName + "/c_").data() + a);
111  m_hSignal[i] = new TH2F(a, a, 6, 0, 6, 4, 0, 4);
112  m_hSignal[i]->SetDirectory(0);// dont mess with it, this is MY histogram
113  m_hSignal[i]->SetStats(false);
114  m_hSignal[i]->SetMinimum(0);
115  m_hSignal[i]->SetMaximum(64);
116 
117  a = "pxdraw/hCommon";
118  a.ReplaceAll("/", "_");
119  a += s2;
120  m_cCommon[i] = new TCanvas((m_histogramDirectoryName + "/c_").data() + a);
121  m_hCommon[i] = new TH2F(a, a, 6, 0, 6, 4, 0, 4);
122  m_hCommon[i]->SetDirectory(0);// dont mess with it, this is MY histogram
123  m_hCommon[i]->SetStats(false);
124  m_hCommon[i]->SetMinimum(0);
125  m_hCommon[i]->SetMaximum(256);
126 
127  a = "pxdraw/hCounts";
128  a.ReplaceAll("/", "_");
129  a += s2;
130  m_cCounts[i] = new TCanvas((m_histogramDirectoryName + "/c_").data() + a);
131  m_hCounts[i] = new TH2F(a, a, 6, 0, 6, 4, 0, 4);
132  m_hCounts[i]->SetDirectory(0);// dont mess with it, this is MY histogram
133  m_hCounts[i]->SetStats(false);
134  }
135 
136  m_fLandau = new TF1("f_Landau", "landau", 0, 256);
137  m_fLandau->SetParameter(0, 1000);
138  m_fLandau->SetParameter(1, 0);
139  m_fLandau->SetParameter(2, 10);
140  m_fLandau->SetLineColor(4);
141  m_fLandau->SetNpx(256);
142  m_fLandau->SetNumberFitPoints(256);
143 
144  m_fGaus = new TF1("f_Gaus", "gaus", 0, 8096);
145  m_fGaus->SetParameter(0, 1000);
146  m_fGaus->SetParameter(1, 0);
147  m_fGaus->SetParameter(2, 10);
148  m_fGaus->SetLineColor(4);
149  m_fGaus->SetNpx(256);
150  m_fGaus->SetNumberFitPoints(256);
151 }
152 
153 
154 void DQMHistAnalysisPXDFitsModule::beginRun()
155 {
156  B2DEBUG(1, "DQMHistAnalysisPXDFits: beginRun called.");
157 
158  // not much we can do here ... as we have created everything already
159  for (auto i = 0; i < NUM_MODULES; i++) {
160  m_cSignal[i]->Clear();
161  m_cCommon[i]->Clear();
162  m_cCounts[i]->Clear();
163  // no need to Cd and Draw yet
164  }
165 }
166 
167 void DQMHistAnalysisPXDFitsModule::event()
168 {
169 // bool flag = false;
170 
171  m_hSignalAll->Reset(); // dont sum up!!!
172  m_hCommonAll->Reset(); // dont sum up!!!
173  m_hCountsAll->Reset(); // dont sum up!!!
174  m_hOccupancyAll->Reset(); // dont sum up!!!
175 
176  for (auto i = 0; i < NUM_MODULES; i++) {
177  auto id = m_inx_to_id[i];
178  auto layer = (((id >> 5) & 0x1) + 1);
179  auto ladder = ((id >> 1) & 0xF);
180  auto sensor = ((id & 0x1) + 1);
181 
182  m_hSignal[i]->Reset(); // dont sum up!!!
183  m_hCommon[i]->Reset(); // dont sum up!!!
184  m_hCounts[i]->Reset(); // dont sum up!!!
185 
186  for (auto j = 0; j < 6; j++) {
187  for (auto k = 0; k < 4; k++) {
188  //TH1* hh1 = NULL;
189  string s2 = str(format("_%d.%d.%d_%d_%d") % layer % ladder % sensor % j % k);
190 
191  std::string name = "hrawPxdHitsCharge" + s2;
192  TH1* hh1 = findHist(name);
193  if (hh1 == NULL) {
194  hh1 = findHist(m_histogramDirectoryName, name);
195  }
196 
197  if (hh1 != NULL) {
198 // cout << "do da fit " << endl;
199 // m_fLandau->SetParameter(0, 1000);
200 // m_fLandau->SetParameter(1, 0);
201 // m_fLandau->SetParameter(2, 10);
202 // hh1->Fit(m_fLandau, "0");
203 // m_hSignal[i]->Fill(j, k, m_fLandau->GetParameter(1));
204 // cout << m_fLandau->GetParameter(0) << " " << m_fLandau->GetParameter(1) << " " << m_fLandau->GetParameter(2) << endl;
205  m_hSignal[i]->Fill(j, k, hh1->GetMean());
206  m_hSignalAll->Fill(i, hh1->GetMean());
207  } else {
208  B2INFO("Histo " << name << " not found");
209  }
210 
211  name = "hrawPxdHitsCommonMode" + s2;
212  hh1 = findHist(name);
213  if (hh1 == NULL) {
214  hh1 = findHist(m_histogramDirectoryName, name);
215  }
216 
217  if (hh1 != NULL) {
218 // cout << "do da fit " << endl;
219 // m_fGaus->SetParameter(0, 1000);
220 // m_fGaus->SetParameter(1, 10);
221 // m_fGaus->SetParameter(2, 10);
222 // hh1->Fit(m_fGaus, "0");
223 // m_hCommon[i]->Fill(j, k, m_fGaus->GetParameter(1));
224 // cout << m_fGaus->GetParameter(0) << " " << m_fGaus->GetParameter(1) << " " << m_fGaus->GetParameter(2) << endl;
225  m_hCommon[i]->Fill(j, k, hh1->GetMean());
226  m_hCommonAll->Fill(i, hh1->GetMean());
227  } else {
228  B2INFO("Histo " << name << " not found");
229  }
230 
231  name = "hrawPxdCount" + s2;
232  hh1 = findHist(name);
233  if (hh1 == NULL) {
234  hh1 = findHist(m_histogramDirectoryName, name);
235  }
236 
237  if (hh1 != NULL) {
238 // cout << "do da fit " << endl;
239 // m_fGaus->SetParameter(0, 1000);
240 // m_fGaus->SetParameter(1, 100);
241 // m_fGaus->SetParameter(2, 10);
242 // hh1->Fit(m_fGaus, "0");
243 // m_hCounts[i]->Fill(j, k, m_fGaus->GetParameter(1));
244 // cout << m_fGaus->GetParameter(0) << " " << m_fGaus->GetParameter(1) << " " << m_fGaus->GetParameter(2) << endl;
245  m_hCounts[i]->Fill(j, k, hh1->GetMean());
246  m_hCountsAll->Fill(i, hh1->GetMean());
247  m_hOccupancyAll->Fill(i, hh1->GetMean() / (250 * 768 / 24)); // Occupancy in percent
248  } else {
249  B2INFO("Histo " << name << " not found");
250  }
251  }
252  }
253  if (m_cSignal[i]) {
254  m_cSignal[i]->cd();
255  if (m_hSignal[i]) m_hSignal[i]->Draw("colz");
256  m_cSignal[i]->Modified();
257  m_cSignal[i]->Update();
258  }
259  if (m_cCommon[i]) {
260  m_cCommon[i]->cd();
261  if (m_hCommon[i]) m_hCommon[i]->Draw("colz");
262  m_cCommon[i]->Modified();
263  m_cCommon[i]->Update();
264  }
265  if (m_cCounts[i]) {
266  m_cCounts[i]->cd();
267  if (m_hCounts[i]) m_hCounts[i]->Draw("colz");
268  m_cCounts[i]->Modified();
269  m_cCounts[i]->Update();
270  }
271  }
272  if (m_cSignalAll) {
273  m_cSignalAll->cd();
274  if (m_hSignalAll) {
275  m_hSignalAll->Scale(1.0 / 24.0); // need to scale
276  m_hSignalAll->Draw();
277  }
278  m_cSignalAll->Modified();
279  m_cSignalAll->Update();
280  }
281  if (m_cCommonAll) {
282  m_cCommonAll->cd();
283  if (m_hCommonAll) {
284  m_hCommonAll->Scale(1.0 / 24.0); // need to scale
285  m_hCommonAll->Draw();
286  }
287  m_cCommonAll->Modified();
288  m_cCommonAll->Update();
289  }
290  if (m_cCountsAll) {
291  m_cCountsAll->cd();
292  if (m_hCountsAll) {
293  // m_hCountsAll->Scale(1.0/24.0); // dont scale counts, we want to have sum
294  // but we would need to scale if we directly calculate occupancy!
295  m_hCountsAll->Draw();
296  }
297  m_cCountsAll->Modified();
298  m_cCountsAll->Update();
299  }
300  if (m_cOccupancyAll) {
301  m_cOccupancyAll->cd();
302  if (m_hOccupancyAll) {
303  m_hOccupancyAll->Scale(1.0 / 24.0); // need to scale
304  m_hOccupancyAll->Draw();
305  }
306  m_cOccupancyAll->Modified();
307  m_cOccupancyAll->Update();
308  }
309 }
310 
311 void DQMHistAnalysisPXDFitsModule::endRun()
312 {
313  B2DEBUG(1, "DQMHistAnalysisPXDFits : endRun called");
314 }
315 
316 
317 void DQMHistAnalysisPXDFitsModule::terminate()
318 {
319  B2DEBUG(1, "DQMHistAnalysisPXDFits: terminate called");
320 }
321 
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::DQMHistAnalysisPXDFitsModule
Class definition for the output module of Sequential ROOT I/O.
Definition: DQMHistAnalysisPXDFits.h:24
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::DQMHistAnalysisModule
The base class for the histogram analysis module.
Definition: DQMHistAnalysis.h:27