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