Belle II Software  release-08-01-10
DQMHistAnalysisARICHMonObj.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 // Own header.
10 #include <dqm/analysis/modules/DQMHistAnalysisARICHMonObj.h>
11 
12 #include <TF1.h>
13 #include <TH1F.h>
14 #include <TH3.h>
15 #include <TH2.h>
16 #include <TLine.h>
17 #include <TH2Poly.h>
18 #include <TText.h>
19 #include <TROOT.h>
20 #include <TLegend.h>
21 #include <iostream>
22 #include <TStyle.h>
23 #include <TGaxis.h>
24 
25 using namespace std;
26 using namespace Belle2;
27 
28 //-----------------------------------------------------------------
29 // Register module
30 //-----------------------------------------------------------------
31 
32 REG_MODULE(DQMHistAnalysisARICHMonObj);
33 
34 DQMHistAnalysisARICHMonObjModule::DQMHistAnalysisARICHMonObjModule()
35  : DQMHistAnalysisModule()
36 {
37  // set module description (e.g. insert text)
38  setDescription("Example module for making MonitoringObject in DQMHistAnalysis module");
39  setPropertyFlags(c_ParallelProcessingCertified);
40 }
41 
42 DQMHistAnalysisARICHMonObjModule::~DQMHistAnalysisARICHMonObjModule()
43 {
44 }
45 
46 void DQMHistAnalysisARICHMonObjModule::initialize()
47 {
48  // make monitoring object related to this module (use arich as an example)
49  // if monitoring object already exists this will return pointer to it
50  m_monObj = getMonitoringObject("arich");
51 
52  // make canvases to be added to MonitoringObject
53  m_c_main = new TCanvas("arich_main", "main", 1500, 800);
54  m_c_mask = new TCanvas("arich_mask", "mask", 750, 1600);
55  m_c_mirror = new TCanvas("arich_mirror", "mirror", 1000, 1000);
56  m_c_tracks = new TCanvas("arich_tracks", "tracks", 500, 500);
57 
58  m_apdHist = new ARICHChannelHist("tmpApdHist", "tmpApdHist", 2);
59  //m_chHist = new ARICHChannelHist("tmpChHist", "tmpChHist", 0); /**<ARICH TObject to draw hit map for each channel*/
60  m_hapdHist = new ARICHChannelHist("tmpHapdHist", "tmpHapdHist", 1);
62  // add canvases to MonitoringObject
63  m_monObj->addCanvas(m_c_main);
64  m_monObj->addCanvas(m_c_mask);
65  m_monObj->addCanvas(m_c_mirror);
66  m_monObj->addCanvas(m_c_tracks);
67 
68 }
69 
70 
71 void DQMHistAnalysisARICHMonObjModule::endRun()
72 {
73 
74 
75  TGaxis::SetMaxDigits(3);
76  gStyle->SetPalette(1);
77  gStyle->SetOptStat(0);
78 
79 
80  // get existing histograms produced by DQM modules
81  TH1* chHit = findHist("ARICH/chHit");
82  TH1* bits = findHist("ARICH/bits");
83  TH1* hitsPerTrack = findHist("ARICH/hitsPerTrack");
84  TH1* theta = findHist("ARICH/theta");
85  TH2* tracks2D = (TH2*)findHist("ARICH/tracks2D");
86  TH1* hitsPerEvent = findHist("ARICH/hitsPerEvent");
87  TH2* hapdHitPerEvent = (TH2*)findHist("ARICH/hapdHitPerEvent");
88  TH2* thetaPhi = (TH2*)findHist("ARICH/thetaPhi");
89  TH3* mirrorThetaPhi = (TH3*)findHist("ARICH/mirrorThetaPhi");
90  TH1* chDigit = findHist("ARICH/chDigit");
91  TH1* hapdDigit = findHist("ARICH/hapdDigit");
92 
93  if (chHit == NULL) {m_monObj->setVariable("comment", "No arich histograms in file"); B2INFO("Histogram named chHit is not found."); return;}
94  if (chHit->GetEntries() == 0) {m_monObj->setVariable("comment", "No arich hits in histograms"); B2INFO("No arich hits in histograms."); return;}
95 
96  // set the content of main canvas
97  m_c_main->Clear(); // clear existing content
98  m_c_main->Divide(3, 2);
99  m_c_mirror->Clear(); // clear existing content
100  m_c_mirror->Divide(3, 3);
101  m_c_mask->Clear(); // clear existing content
102  m_c_mask->Divide(2, 4);
103  m_c_tracks->Clear(); // clear existing content
104 
105  //Draw 2D hit map of channels and APDs
106  pp1 = new TH2Poly();
107  pp1->SetTitle("Number of hits / APD / event");
108  pp1->SetOption("colz");
109 
110  pp2 = new TH2Poly();
111  pp2->SetMaximum(0.1 / 36.);
112  pp2->SetMinimum(0.0001 / 36.);
113  pp2->SetTitle("Number of hits / channel / event");
114 
115  pflash = new TH2Poly();
116  pflash->SetTitle("Number of flash (>40 hits) / event");
117 
118  int nevt = 0;
119  if (hitsPerEvent) nevt = hitsPerEvent->GetEntries();
120  m_apdHist->fillFromTH1(chHit);
121  // m_chHist->fillFromTH1(chHit);
122  if (nevt) {
123  m_apdHist->Scale(1. / float(nevt));
124  // m_chHist->Scale(1. / float(nevt));
125  }
126 
127  pp1->SetMaximum(0.1);
128  m_apdHist->setPoly(pp1);
129  pp1->SetMinimum(0.0001);
130  pp2->SetMaximum(0.1 / 36.);
131  // m_chHist->setPoly(pp2);
132  // pp2->SetMinimum(0.0001 / 36.);
133 
134 
135  //TCanvas main
136  m_c_main->cd(1);
137  pp1->Draw("colz");
138  pp1->GetXaxis()->SetTickLength(0);
139  pp1->GetYaxis()->SetTickLength(0);
140  gPad->SetLogz();
141  //m_c_main->Update();
142 
143  TH1F* flash = (TH1F*)hapdHitPerEvent->ProjectionX("flash", 40, 144);
144  m_hapdHist->fillFromTH1(flash);
145  if (nevt) m_hapdHist->Scale(1. / float(nevt));
146 
147  m_hapdHist->setPoly(pflash);
148  // draw sector lines
149  double rlin = 40;
150  double rlout = 113;
151  for (int isec = 0; isec < 6; isec++) {
152  double x1 = rlin * cos(M_PI / 3.*isec);
153  double x2 = rlout * cos(M_PI / 3.*isec);
154  double y1 = rlin * sin(M_PI / 3.*isec);
155  double y2 = rlout * sin(M_PI / 3.*isec);
156  TLine* line = new TLine(x1, y1, x2, y2);
157  line->Draw();
158  x1 = rlin * cos(M_PI / 3.*isec + M_PI / 6.);
159  y1 = rlin * sin(M_PI / 3.*isec + M_PI / 6.);
160  TText* lab = new TText(x1, y1, TString::Format("S-%d", isec + 1));
161  lab->SetTextAlign(22);
162  lab->SetTextSize(0.03);
163  lab->Draw();
164  }
165 
166 
167  //TCanvas tracks
168  m_c_tracks->cd();
169  tracks2D->RebinX(4);
170  tracks2D->RebinY(4);
171  double trkevt = nevt > 0 ? tracks2D->GetEntries() / nevt : 0;
172  int ntracks = tracks2D->GetEntries();
173  m_monObj->setVariable("ntracks", ntracks ? ntracks : 0);
174  m_monObj->setVariable("ntracksPerEvent", trkevt ? trkevt : 0);
175  string comment = "";
176  if (ntracks == 0) comment = "No arich tracks in file.";
177  if (theta->GetEntries() == 0) comment.append(" No cherenkov photons available");
178  m_monObj->setVariable("comment", comment);
179 
180  tracks2D->SetTitle(TString::Format("Track distribution (avg %f trk/evt)", trkevt));
181  tracks2D->Draw("colz");
182 
183 
184  if (ntracks) {
185  double sigbkg[8] = {0};
186  //fit two gauses
187  if (theta->GetEntries() == 0) return;
188  TF1* f1 = new TF1("thcFit", "gaus(0)+gaus(3)", 0.2, 0.4);
189  f1->SetParameters(0.8 * theta->GetMaximum(), 0.323, 0.016, 0.2 * theta->GetMaximum(), 0.323, 0.13);
190  f1->FixParameter(5, 0.13);
191  f1->SetParName(0, "C");
192  f1->SetParName(1, "mean");
193  f1->SetParName(2, "sigma");
194  f1->SetParName(3, "p0");
195  f1->SetParName(4, "p1");
196  int status = theta->Fit(f1, "R");
197  double xmin = f1->GetParameter(1) - 2.*f1->GetParameter(2);
198  double xmax = f1->GetParameter(1) + 2.*f1->GetParameter(2);
199  double tmp = f1->GetParameter(3);
200  f1->SetParameter(3, 0);
201  double nphot = f1->Integral(xmin, xmax);
202  f1->SetParameter(3, tmp);
203  tmp = f1->GetParameter(0);
204  f1->SetParameter(0, 0);
205  double nbkg = f1->Integral(xmin, xmax);
206  f1->SetParameter(0, tmp);
207  if (status) return;
208 
209  sigbkg[0] = nphot;
210  sigbkg[1] = f1->GetParameter(0) > 0 ? nphot * f1->GetParError(0) / f1->GetParameter(0) : 0.;
211  sigbkg[2] = nbkg;
212  sigbkg[3] = f1->GetParameter(3) > 0 ? nbkg * f1->GetParError(3) / f1->GetParameter(3) : 0.;
213  sigbkg[4] = f1->GetParameter(1);
214  sigbkg[5] = f1->GetParError(1);
215  sigbkg[6] = f1->GetParameter(2);
216  sigbkg[7] = f1->GetParError(2);
217  if (theta->GetBinWidth(1) == 0) return;
218  m_monObj->setVariable("nsig", sigbkg[0] / float(ntracks) / theta->GetBinWidth(1),
219  sigbkg[1] / float(ntracks) / theta->GetBinWidth(1));
220  m_monObj->setVariable("nbgr", sigbkg[2] / float(ntracks) / theta->GetBinWidth(1),
221  sigbkg[3] / float(ntracks) / theta->GetBinWidth(1));
222  m_monObj->setVariable("theta", sigbkg[4], sigbkg[5]);
223  m_monObj->setVariable("sigma", sigbkg[6], sigbkg[7]);
224  std::cout << sigbkg[0] << " " << sigbkg[1] / float(ntracks) / theta->GetBinWidth(1) << std::endl;
225  }
226 
227 
228  //TCanvas mirror
229  TH1F* thetaCl = (TH1F*)theta->Clone("thetaCl");
230  thetaCl->SetLineColor(16);
231  thetaCl->SetLineWidth(2);
232  thetaCl->SetTitle("");
233  TLegend* leg[9];
234  gStyle->SetOptTitle(0);
235 
236  if (mirrorThetaPhi) {
237  for (int i = 1; i < 18 + 1; i++) {
238  TH1F* hmir = (TH1F*)mirrorThetaPhi->ProjectionZ(TString::Format("hmir_%d", i), i, i, 1, 10000);
239  hmir->SetTitle(TString::Format("mirror %d", i));
240  if (hmir->GetEntries() > 0) hmir->Scale(theta->GetEntries() / hmir->GetEntries());
241  hmir->Rebin(2);
242  hmir->SetLineWidth(2);
243  int iplot = (i - 1) / 2;
244  if ((i - 1) % 2 == 0) { m_c_mirror->cd(iplot + 1); thetaCl->Draw("hist"); hmir->SetLineColor(1); hmir->Draw("sames hist"); leg[iplot] = new TLegend(0.1, 0.75, 0.4, 0.9); leg[iplot]->AddEntry(hmir, TString::Format("mirror %d", i));}
245  else {
246  hmir->SetLineColor(2); hmir->Draw("sames hist");
247  leg[iplot]->AddEntry(hmir, TString::Format("mirror %d", i));
248  leg[iplot]->Draw();
249  }
250  }
251  }
252 
253  gStyle->SetOptTitle(1);
254 
255 
256  //chDigit
257  if (chDigit != NULL && nevt) chDigit->Scale(1. / nevt);
258  if (nevt) {
259  chHit->Scale(1. / nevt);
260  flash->Scale(1. / nevt);
261  }
262  int nhot = 0, ndead = 0;
263  TH2F* hotCh = new TH2F("arich_hot", "Number of channels in APD with >0.5% occ.", 42, 0.5, 42.5, 40, 0.5, 40.5);
264  TH2F* hotCh1 = new TH2F("arich_hot1", "Number of channels in APD with >0.5% occ. after mask", 42, 0.5, 42.5, 40, 0.5, 40.5);
265  TH2F* deadCh = new TH2F("arich_dead", "Number of channels in APD with no hits", 42, 0.5, 42.5, 40, 0.5, 40.5);
266  TH2F* falseCh = new TH2F("arich_false", "Number of wrongly masked channels in APD (masked but not dead/hot)", 42, 0.5, 42.5, 40,
267  0.5,
268  40.5);
269  TH1F* occ = new TH1F("arich_occ", "nhits / nevt for all channels; nhits/nevt;# of chn", 500, 0, 0.005);
270 
271  int ndeadHapd = 0;
272  if (chDigit != NULL) {
273  double hotlim = 0.005;
274  for (int i = 0; i < chDigit->GetNbinsX(); i++) {
275  int hapd = i / 144 + 1; int chip = (i % 144) / 36 + 1;
276  int binx = (hapd - 1) % 42 + 1; int biny = ((hapd - 1) / 42) * 4 + chip;
277  if (chDigit->GetBinContent(i + 1) > hotlim) { hotCh->Fill(binx, biny); nhot++;}
278  if (chDigit->GetBinContent(i + 1) == 0) {
279  deadCh->Fill(binx, biny);
280  if (hapdDigit->GetBinContent(hapd) != 0) ndead++;
281  }
282  if (chHit->GetBinContent(i + 1) == 0 && chDigit->GetBinContent(i + 1) < hotlim
283  && chDigit->GetBinContent(i + 1) > 0.00002) falseCh->Fill(binx, biny);
284  if (chHit->GetBinContent(i + 1) > hotlim) hotCh1->Fill(binx, biny);
285  occ->Fill(chHit->GetBinContent(i + 1));
286  }
287 
288  for (int i = 0; i < hapdDigit->GetNbinsX(); i++) {
289  if (hapdDigit->GetBinContent(i + 1) == 0) ndeadHapd++;
290  }
291  }
292  m_monObj->setVariable("nhot", nhot);
293  m_monObj->setVariable("ndead", ndead);
294  m_monObj->setVariable("ndeadHapd", ndeadHapd);
295 
296 
297  // TCanvas mask
298  m_c_mask->cd(1);
299  occ->Draw();
300  m_c_mask->cd(2);
301  pflash->Draw("colz");
302  m_c_mask->cd(3);
303  // pp2->Draw("colz");
304  //pp2->GetXaxis()->SetTickLength(0);
305  //pp2->GetYaxis()->SetTickLength(0);
306  gPad->SetLogz();
307  m_c_mask->cd(4);
308  hotCh->Draw("colz text");
309  m_c_mask->cd(5);
310  deadCh->Draw("colz text");
311  m_c_mask->cd(6);
312  hotCh1->Draw("colz text");
313  m_c_mask->cd(7);
314  falseCh->Draw("colz text");
315 
316  for (int i = 0; i < 42; i++) {
317  double x1 = i + 0.5 + 1;
318  TLine* line = new TLine(x1, 0.5, x1, 40.5);
319  m_c_mask->cd(5); line->Draw();
320  m_c_mask->cd(6); line->Draw();
321  m_c_mask->cd(7); line->Draw();
322  m_c_mask->cd(4); line->Draw();
323  }
324  for (int i = 0; i < 40; i++) {
325  double y1 = i + 0.5 + 1;
326  TLine* line = new TLine(0.5, y1, 42.5, y1);
327  if ((i + 1) % 4 == 0) line->SetLineColor(2);
328  else line->SetLineColor(15);
329  m_c_mask->cd(5); line->Draw();
330  m_c_mask->cd(6); line->Draw();
331  m_c_mask->cd(7); line->Draw();
332  m_c_mask->cd(4); line->Draw();
333  }
334 
335  if (bits) {
336  bits->SetLineWidth(2);
337  bits->SetLineColor(2);
338  bits->SetOption("hist");
339  bits->SetFillStyle(3010);
340  bits->SetFillColor(3);
341  }
342 
343  if (hitsPerTrack) {
344  hitsPerTrack->SetLineWidth(2);
345  hitsPerTrack->SetLineColor(2);
346  hitsPerTrack->SetOption("hist");
347  }
348  if (hitsPerEvent) {
349  hitsPerEvent->SetLineWidth(2);
350  hitsPerEvent->SetLineColor(2);
351  hitsPerEvent->SetOption("hist");
352  }
353  //Tcanvas main
354  m_c_main->cd(2);
355  if (bits) { bits->Draw(); bits->GetYaxis()->SetTitleOffset(0.5); }
356  m_c_main->cd(3);
357  if (theta) theta->Draw();
358  m_c_main->cd(4);
359  if (hitsPerTrack) hitsPerTrack->Draw();
360  m_c_main->cd(5);
361  if (hitsPerEvent) hitsPerEvent->Draw();
362  m_c_main->cd(6);
363  if (thetaPhi != NULL) thetaPhi->Draw("colz");
364 
365  // set values of monitoring variables (if variable already exists this will change its value, otherwise it will insert new variable)
366  // with error (also asymmetric error can be used as m_monObj->setVariable(name, value, upError, downError))
367  m_monObj->setVariable("hitsPerEvent", hitsPerEvent ? hitsPerEvent->GetMean() : 0, hitsPerEvent ? hitsPerEvent->GetMeanError() : -1);
368  // without error
369  m_monObj->setVariable("bitsMean", bits ? bits->GetMean() : 0);
370  B2DEBUG(20, "DQMHistAnalysisARICHMonObj : endRun called");
371 }
372 
373 void DQMHistAnalysisARICHMonObjModule::terminate()
374 {
375 
376  B2DEBUG(20, "terminate called");
377 }
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::fillFromTH1
void fillFromTH1(TH1 *hist)
Fill the channelHist from the histogram Type 0 channelHist has to be filled with 420*144bin TH1 (each...
Definition: ARICHChannelHist.cc:203
Belle2::ARICHChannelHist::setPoly
void setPoly(TH2Poly *poly)
Fill pure TH2Poly from ARICHChannelHist, makes bins and fills content.
Definition: ARICHChannelHist.cc:236
Belle2::DQMHistAnalysisARICHMonObjModule::initialize
void initialize() override final
Initialize the Module.
Definition: DQMHistAnalysisARICHMonObj.cc:46
Belle2::DQMHistAnalysisARICHMonObjModule::m_apdHist
Belle2::ARICHChannelHist * m_apdHist
ARICH TObject to draw hit map for each APD.
Definition: DQMHistAnalysisARICHMonObj.h:66
Belle2::DQMHistAnalysisARICHMonObjModule::m_monObj
MonitoringObject * m_monObj
MonitoringObject to be produced by this module.
Definition: DQMHistAnalysisARICHMonObj.h:69
Belle2::DQMHistAnalysisARICHMonObjModule::m_hapdHist
Belle2::ARICHChannelHist * m_hapdHist
ARICH TObject to draw flash map for each hapd.
Definition: DQMHistAnalysisARICHMonObj.h:68
Belle2::DQMHistAnalysisARICHMonObjModule::terminate
void terminate() override final
Termination action.
Definition: DQMHistAnalysisARICHMonObj.cc:373
Belle2::DQMHistAnalysisARICHMonObjModule::m_c_mask
TCanvas * m_c_mask
Canvas with histograms related to channel masking.
Definition: DQMHistAnalysisARICHMonObj.h:59
Belle2::DQMHistAnalysisARICHMonObjModule::m_c_main
TCanvas * m_c_main
Canvas with main run summary histograms.
Definition: DQMHistAnalysisARICHMonObj.h:58
Belle2::DQMHistAnalysisARICHMonObjModule::m_c_mirror
TCanvas * m_c_mirror
Canvas with histograms related to mirrors.
Definition: DQMHistAnalysisARICHMonObj.h:60
Belle2::DQMHistAnalysisARICHMonObjModule::endRun
void endRun() override final
End-of-run action.
Definition: DQMHistAnalysisARICHMonObj.cc:71
Belle2::DQMHistAnalysisARICHMonObjModule::pp2
TH2Poly * pp2
2D hitmap of number of hits per channel per event
Definition: DQMHistAnalysisARICHMonObj.h:63
Belle2::DQMHistAnalysisARICHMonObjModule::m_c_tracks
TCanvas * m_c_tracks
Canvas with histograms related to tracks.
Definition: DQMHistAnalysisARICHMonObj.h:61
Belle2::DQMHistAnalysisARICHMonObjModule::pflash
TH2Poly * pflash
2D hitmap of number of flash (>40 hits) per event
Definition: DQMHistAnalysisARICHMonObj.h:64
Belle2::DQMHistAnalysisARICHMonObjModule::pp1
TH2Poly * pp1
2D hitmap of number of hits per APD per event
Definition: DQMHistAnalysisARICHMonObj.h:62
Belle2::DQMHistAnalysisModule
The base class for the histogram analysis module.
Definition: DQMHistAnalysis.h:39
Belle2::DQMHistAnalysisModule::findHist
static TH1 * findHist(const std::string &histname, bool onlyIfUpdated=false)
Get histogram from list (no other search).
Definition: DQMHistAnalysis.cc:142
Belle2::DQMHistAnalysisModule::getMonitoringObject
static MonitoringObject * getMonitoringObject(const std::string &histname)
Get MonitoringObject with given name (new object is created if non-existing)
Definition: DQMHistAnalysis.cc:112
Belle2::Module::setDescription
void setDescription(const std::string &description)
Sets the description of the module.
Definition: Module.cc:214
Belle2::Module::setPropertyFlags
void setPropertyFlags(unsigned int propertyFlags)
Sets the flags for the module properties.
Definition: Module.cc:208
Belle2::Module::c_ParallelProcessingCertified
@ c_ParallelProcessingCertified
This module can be run in parallel processing mode safely (All I/O must be done through the data stor...
Definition: Module.h:80
Belle2::MonitoringObject::setVariable
void setVariable(const std::string &var, float val, float upErr=-1., float dwErr=-1)
set value to float variable (new variable is made if not yet existing)
Definition: MonitoringObject.h:124
Belle2::MonitoringObject::addCanvas
void addCanvas(TCanvas *canv)
Add Canvas to monitoring object.
Definition: MonitoringObject.h:48
Belle2::hapd
TH1F * hapd[6]
histogram of hits for each hapd
Definition: arichBtestModule.cc:81
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:650
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17