12 #include <dqm/analysis/modules/DQMHistAnalysisCDCMonObj.h>
15 #include <dqm/analysis/modules/DQMHistAnalysis.h>
18 #include <cdc/geometry/CDCGeometryPar.h>
48 DQMHistAnalysisCDCMonObjModule::DQMHistAnalysisCDCMonObjModule()
52 setDescription(
"Modify and analyze the data quality histograms of CDCMonObj");
54 for (
int i = 0; i < 300; i++) {
59 for (
int i = 0; i < 56; i++)
m_hHits[i] =
nullptr;
71 if (!(*m_channelMapFromDB).isValid()) {
72 B2FATAL(
"Channel map is not valid");
77 B2FATAL(
"CDCGeometryp is not valid");
84 gStyle->SetOptStat(0);
85 gStyle->SetPalette(kViridis);
86 gStyle->SetPadTopMargin(0.1);
87 gStyle->SetPadRightMargin(0.05);
88 gStyle->SetPadBottomMargin(0.1);
89 gStyle->SetPadLeftMargin(0.15);
91 m_cMain =
new TCanvas(
"cdc_main",
"cdc_main", 1500, 1000);
94 m_cADC =
new TCanvas(
"cdc_adc",
"cdc_adc", 2000, 10000);
96 m_cTDC =
new TCanvas(
"cdc_tdc",
"cdc_tdc", 2000, 10000);
98 m_cHit =
new TCanvas(
"cdc_hit",
"cdc_hit", 1500, 6000);
101 B2DEBUG(20,
"DQMHistAnalysisCDCMonObj: initialized.");
108 const int isl = cm.getISuperLayer();
109 const int il = cm.getILayer();
110 const int iw = cm.getIWire();
111 const int iBoard = cm.getBoardID();
112 const int iCh = cm.getBoardChannel();
113 const WireID wireId(isl, il, iw);
114 m_chMap.insert(std::make_pair(wireId, std::make_pair(iBoard, iCh)));
120 for (
const auto& sense : geom.getSenseLayers()) {
121 int i = sense.getId();
122 if (i < 0 || i > 55) {
123 B2FATAL(
"no such sense layer");
130 for (
const auto& field : geom.getFieldLayers()) {
131 int i = field.getId();
132 if (i < 0 || i > 54) {
133 B2FATAL(
"no such sense layer");
137 m_fieldR[56] = geom.getOuterWall(0).getRmin();
138 m_fieldR[0] = geom.getInnerWall(0).getRmax();
148 for (
int ilayer = 0; ilayer < 56; ++ilayer) {
152 for (
int iwire = 0; iwire <
m_nSenseWires[ilayer]; ++iwire) {
153 float phi = dPhi * (iwire +
m_offset[ilayer]);
154 float phi1 = phi - dPhi * 0.5;
155 float phi2 = phi + dPhi * 0.5;
156 Double_t x_pos[] = {r1* (sin(phi)*tan(phi - phi1) + cos(phi)),
159 r1* (sin(phi)*tan(phi - phi2) + cos(phi))
161 Double_t y_pos[] = {r1* (-cos(phi)*tan(phi - phi1) + sin(phi)),
164 r1* (-cos(phi)*tan(phi - phi2) + sin(phi))
166 h->AddBin(4, x_pos, y_pos);
175 for (
int il = 0; il < 56; ++il) {
177 const int y =
m_hHits[il]->GetBinContent(iw + 1);
188 TH1D* hist = (TH1D*)h->Clone();
189 hist->SetBinContent(1, 0.0);
190 float m = hist->GetMean();
197 const WireID w(layer, wire);
202 B2ERROR(
"no corresponding board/channel found layer " << layer <<
" wire " << wire);
203 return std::make_pair(-1, -1);
210 B2DEBUG(20,
"end run");
216 TF1* fitFunc[300] = {};
217 for (
int i = 0; i < 300; ++i) {
218 fitFunc[i] =
new TF1(Form(
"f%d", i),
"[0]+[6]*x+[1]*(exp([2]*(x-[3]))/(1+exp(-([4]-x)/[5])))",
219 4921 - 100, 4921 + 100);
220 fitFunc[i]->SetParLimits(0, 0, 100);
221 fitFunc[i]->SetParLimits(6, 0, 0.1);
222 fitFunc[i]->SetParLimits(4, 4850., 5000.0);
223 fitFunc[i]->SetParLimits(5, 0, 50.0);
227 if (neve == 0)neve = 1;
229 B2DEBUG(20,
"adc related");
232 TH1F* hADCMean =
new TH1F(
"hADCMean",
"ADC mean;board;adc mean", 300, 0, 300);
233 TH1F* hADC1000 =
new TH1F(
"ADC1000",
"ADC1000", 300, 0, 300);
234 TH1F* hADC0 =
new TH1F(
"ADC0",
"ADC0", 300, 0, 300);
236 std::vector<float> means = {};
237 for (
int i = 0; i < 300; ++i) {
240 m_hADCs[i] =
m_hADC->ProjectionY(Form(
"hADC%d", i), i + 1, i + 1,
"");
241 m_hADCs[i]->SetTitle(Form(
"hADC%d", i));
242 float n =
static_cast<float>(
m_hADCs[i]->GetEntries());
243 if (
m_hADCs[i]->GetEntries() == 0) {
245 hADC0->SetBinContent(i + 1, -0.1);
247 float n0 =
static_cast<float>(
m_hADCs[i]->GetBinContent(1));
249 B2DEBUG(21,
"bad adc bid " << i <<
" " << n0 <<
" " << n);
252 float bin1 =
m_hADCs[i]->GetBinContent(1);
255 hADCMean->SetBinContent(i + 1, m);
256 hADCMean->SetBinError(i + 1, 0);
257 double overflow =
m_hADCs[i]->GetBinContent(1001);
258 hADC1000->SetBinContent(i + 1, overflow / (overflow + n));
259 hADC0->SetBinContent(i + 1, bin1 / (overflow + n));
264 B2DEBUG(20,
"tdc related");
266 TH1F* hTDCEdge =
new TH1F(
"hTDCEdge",
"TDC edge;board;tdc edge [nsec]", 300, 0, 300);
267 TH1F* hTDCSlope =
new TH1F(
"hTDCSlope",
"TDC slope;board;tdc slope [nsec]", 300, 0, 300);
268 std::vector<float> tdcEdges = {};
269 std::vector<float> tdcSlopes = {};
270 for (
int i = 0; i < 300; ++i) {
271 m_hTDCs[i] =
m_hTDC->ProjectionY(Form(
"hTDC%d", i), i + 1, i + 1);
272 m_hTDCs[i]->SetTitle(Form(
"hTDC%d", i));
276 fitFunc[i]->SetParameters(0, 100, 0.01, 4700, 4900, 2, 0.01);
277 fitFunc[i]->SetParameter(0, 10);
278 fitFunc[i]->SetParameter(6, 0.02);
281 xxx =
m_hTDCs[i]->Fit(fitFunc[i],
"qM0",
"", 4850, 5000);
283 xxx =
m_hTDCs[i]->Fit(fitFunc[i],
"qM0",
"", 4800, 5000);
285 float p4 = fitFunc[i]->GetParameter(4);
286 float p5 = fitFunc[i]->GetParameter(5);
287 if (xxx != -1 && 4850 < p4 && p4 < 5000) {
288 hTDCEdge->SetBinContent(i + 1, p4);
289 hTDCEdge->SetBinError(i + 1, 0);
290 hTDCSlope->SetBinContent(i + 1, p5);
291 hTDCSlope->SetBinError(i + 1, 0);
292 tdcEdges.push_back(p4);
293 tdcSlopes.push_back(p5);
301 B2DEBUG(20,
"hit related");
302 TH1F* hHitPerLayer =
new TH1F(
"hHitPerLayer",
"hit/Layer;layer", 56, 0, 56);
304 for (
int i = 0; i < 56; ++i) {
305 m_hHits[i] =
m_hHit->ProjectionY(Form(
"hHit%d", i), i + 1, i + 1);
306 m_hHits[i]->SetTitle(Form(
"hHit%d", i));
310 for (
int j = 0; j < nBins; ++j) {
311 nhitSumL +=
m_hHits[i]->GetBinContent(j + 1);
314 hHitPerLayer->SetBinContent(i + 1,
static_cast<float>(nhitSumL / neve));
315 }
else hHitPerLayer->SetBinContent(i + 1,
static_cast<float>(nhitSumL));
316 hHitPerLayer->SetBinError(i + 1, 0);
322 B2DEBUG(20,
"bad wire related");
323 TH2F* hBadChannel =
new TH2F(
"hbadch",
"bad channel map;wire;layer", 400, 0, 400, 56, 0, 56);
324 for (
int i = 0; i < 400; ++i) {
325 for (
int j = 0; j < 56; ++j) {
326 hBadChannel->Fill(i, j, -1);
330 TH2F* hBadChannelBC =
new TH2F(
"hbadchBC",
"bad channel map per board/channel;board;channel", 300, 0, 300, 48, 0, 48);
331 for (
int i = 0; i < 300; ++i) {
332 for (
int j = 0; j < 48; ++j) {
333 hBadChannelBC->Fill(i, j, -1);
337 TH2Poly* h2p =
new TH2Poly();
339 h2p->SetTitle(
"bad wires in xy view");
340 h2p->GetXaxis()->SetTitle(
"X [cm]");
341 h2p->GetYaxis()->SetTitle(
"Y [cm]");
344 const int l = lw.first;
345 const int w = lw.second;
346 B2DEBUG(21,
"l " << l <<
" w " << w);
347 hBadChannel->Fill(w, l);
349 hBadChannelBC->Fill(bc.first, bc.second);
351 float dPhi =
static_cast<float>(2.0 * M_PI /
m_nSenseWires[l]);
352 float phi = dPhi * (w +
m_offset[l]);
353 float x = r * cos(phi);
354 float y = r * sin(phi);
355 h2p->Fill(x, y, 1.1);
358 B2DEBUG(20,
"writing");
362 hADCMean->SetMinimum(0);
363 hADCMean->SetMaximum(300);
367 hTDCEdge->SetMinimum(4800);
368 hTDCEdge->SetMaximum(5000);
372 hTDCSlope->SetMinimum(0);
373 hTDCSlope->SetMaximum(50);
376 hBadChannel->Draw(
"col");
379 hBadChannelBC->Draw(
"col");
381 hHitPerLayer->Draw();
388 for (
int i = 0; i < 56; i++) {
394 m_cADC->Divide(6, 50, 0.0002, 0.0002);
395 m_cTDC->Divide(6, 50, 0.0002, 0.0002);
397 for (
int i = 0; i < 300; i++) {
401 m_hADCs[i]->GetYaxis()->SetRangeUser(0, 3 * max);
407 fitFunc[i]->SetLineColor(kRed);
408 fitFunc[i]->Draw(
"same");
409 max =
m_hTDCs[i]->GetMaximum();
410 TLine* l1 =
new TLine(tdcEdges[i], 0, tdcEdges[i], max * 1.05);
411 l1->SetLineColor(kRed);
413 TLine* l0 =
new TLine(4900, 0, 4900, max * 1.05);
419 float superLayerR[10] = {16.3, 24.3, 35.66, 46.63, 57.55, 68.47,
420 79.39, 90.31, 101.23, 112.05
424 for (
int i = 0; i < 10; ++i) {
425 circs[i] =
new TEllipse(0, 0, superLayerR[i], superLayerR[i]);
426 circs[i]->SetFillStyle(4000);
427 circs[i]->SetLineStyle(kDashed);
428 circs[i]->SetLineColor(0);
429 circs[i]->Draw(
"same");
435 m_monObj->
setVariable(
"adcMean", std::accumulate(means.begin(), means.end(), 0.0) / means.size());
438 m_monObj->
setVariable(
"tdcEdge", std::accumulate(tdcEdges.begin(), tdcEdges.end(), 0.0) / tdcEdges.size());
440 m_monObj->
setVariable(
"tdcSlope", std::accumulate(tdcSlopes.begin(), tdcSlopes.end(), 0.0) / tdcSlopes.size());
449 delete hBadChannelBC;
457 B2DEBUG(20,
"terminate called");