StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
bsQaShow.C
1 #include "TCanvas.h"
2 #include "TF1.h"
3 #include "TFile.h"
4 #include "TGaxis.h"
5 #include "TGraph.h"
6 #include "TGraphErrors.h"
7 #include "TH1.h"
8 #include "TH2.h"
9 #include "TMath.h"
10 #include "TLegend.h"
11 #include "TLine.h"
12 #include "TPaveText.h"
13 #include "TString.h"
14 #include "TSystem.h"
15 #include "TStyle.h"
16 
17 #include <cmath>
18 #include <fstream>
19 #include <iostream>
20 #include <map>
21 using namespace std;
22 
23 //=============================
24 void NormalizeBsHisto(TH2F* H2)
25 {
26  //Normalize bit shift QA histogram for the case of it merged from multiple files
27  const float max = H2->GetMaximum();
28  H2->Scale(1/max);
29 
30  //Remove ambigue bins: only bins with value 1 is valid
31  for (int x=0; x<H2->GetNbinsX(); x++)
32  for (int y=0; y<H2->GetNbinsY(); y++)
33  {
34  const float BC = H2->GetBinContent(x+1, y+1);
35  if (BC < 1.0)
36  {
37  H2->SetBinContent(x+1, y+1, 0);
38  H2->SetBinError (x+1, y+1, 0);
39  }
40  }
41 
42  return;
43 }//NormalizeBsHisto
44 
45 //======================================================================================
46 void bsQaShow(const char* inFile = "fmsBsQa.root", bool MARK = true, bool PRINT = false)
47 {
48  enum {db, data};
49  const int nSep = 75; //# of channels to be drawn per separated plot
50 
51  TFile* F = TFile::Open(inFile);
52  if (!F || F->IsZombie()) { cout <<"Cannot open " <<inFile <<endl; return; }
53 
54  //-------------------------------------------
55 
56  //Get channel map from histograms
57  map<int, int> nToCh[4];
58  for (int a=0; a<4; a++) //detId
59  {
60  const int detId = a+8;
61  TH2F* tempChMap = (TH2F*)F->Get(Form("ChMap_d%i", detId));
62  if (!tempChMap) { cout <<"Cannot open channel map histogram!" <<endl; return; }
63 
64  int index = 0;
65  for (int x=0; x<tempChMap->GetNbinsX(); x++)
66  {
67  const int ch = x+1;
68  if (tempChMap->GetBinContent(x+1, 1) == 0) continue;
69  else
70  {
71  nToCh[a].insert(std::pair<int, int>(index, ch));
72  index++;
73  }
74  }
75  tempChMap->Delete();
76  //for (unsigned int x=0; x<nToCh[a].size(); x++) cout <<Form("%2i %3i %3i", detId, x, nToCh[a][x]) <<endl;
77  }//a, detId
78 
79  //-------------------------------------------
80 
81  #if 1
82  //Show QA: remove invalid channels and draw x channels per plot
83  for (int a=0; a<4; a++) //detId
84  {
85  const int nSepHisto = nToCh[a].size()/nSep + 1;
86  TH2F* H2Sep[2][nSepHisto];
87  for (int b=0; b<2; b++) //DB or data
88  {
89  TH2F* H2_BS = (TH2F*)F->Get(Form("BitShift_%s_d%i", b==0?"DB":"data", a+8));
90  TH2F* H2_ADC = (TH2F*)F->Get(Form("Adc_d%i", a+8));
91  if (!H2_BS || !H2_ADC) { cout <<"Cannot open histograms! (1)" <<endl; return; }
92  NormalizeBsHisto(H2_BS);
93 
94  //Get average # of hit for this detId
95  const float nCh = (H2_ADC->GetNbinsX()>300)?394:238;
96  const float nHitAll = H2_ADC->GetEntries();
97  const float nHitAvg = nHitAll/nCh;
98 
99  //Remove invalid channels + Group nSep channels per plot
100  int Seq = 0;
101  int Count = 0;
102  const int nBit = H2_BS->GetNbinsY();
103  for (unsigned int x=0; x<nToCh[a].size(); x++)
104  {
105  //Declare new sub histogram
106  if (Count == 0)
107  {
108  if (x != 0) Seq++;
109  const char* tempN = Form("%s_%i", H2_BS->GetName(), Seq);
110  const char* tempT = Form("%s, from ch = %i;ch;Bit (BS)", H2_BS->GetTitle(), nToCh[a][nSep*Seq]);
111  H2Sep[b][Seq] = new TH2F(tempN, tempT, nSep,0,nSep, nBit,0.5,nBit+0.5);
112  for (int i=0; i<nSep; i++)
113  {
114  const int chTarget = i + (nSep * Seq);
115  if (nToCh[a][chTarget] == 0) continue;
116  H2Sep[b][Seq]->GetXaxis()->SetBinLabel(i+1, Form("%i", nToCh[a][chTarget]));
117  }
118  H2Sep[b][Seq]->LabelsOption("v", "X");
119  for (int y=0; y<nBit; y++)
120  {
121  if (y+1 <= 6) H2Sep[b][Seq]->GetYaxis()->SetBinLabel(y+1, Form("%2i (#plus%i)", y+1, y+1));
122  else H2Sep[b][Seq]->GetYaxis()->SetBinLabel(y+1, Form("%2i (#minus%i)", y+1, abs(y-12)));
123  }
124  H2Sep[b][Seq]->GetYaxis()->SetTitleOffset(1.2);
125  }
126 
127  //Copy bin content to new sub histogram
128  const int xBinTarget = H2_BS->GetXaxis()->FindBin(nToCh[a][x] + 1.e-3);
129  const int xBinToFill = x - (nSep * Seq) + 1;
130  int nBitFilled = 0;
131  for (int y=0; y<nBit; y++)
132  {
133  float BC = H2_BS->GetBinContent(xBinTarget, y+1); if (BC==0) continue;
134  float BE = H2_BS->GetBinError (xBinTarget, y+1);
135  H2Sep[b][Seq]->SetBinContent(xBinToFill, y+1, BC);
136  H2Sep[b][Seq]->SetBinError (xBinToFill, y+1, BE);
137  nBitFilled++;
138  }
139 
140  //Check if this channel is good or not: 0 for good, 1 for dead, and 2 for bad
141  int chQuality = 0;
142  if (nBitFilled == nBit) chQuality = 1;
143  else if (abs(nBitFilled) > 5)
144  {
145  TH1F* H1_ADC_1ch = (TH1F*)H2_ADC->ProjectionY("", nToCh[a][x], nToCh[a][x]);
146  if (H1_ADC_1ch->GetEntries() < nHitAvg*1.e-2) chQuality = 1;
147  else if (H1_ADC_1ch->GetRMS() < 1.) chQuality = 2;
148  H1_ADC_1ch->Delete();
149  }
150 
151  //Mark bad channels
152  if (MARK==true && chQuality!=0)
153  {
154  for (int y=0; y<nBit; y++)
155  {
156  if (chQuality==1) H2Sep[b][Seq]->SetBinContent(xBinToFill, y+1, 0.1);
157  if (chQuality==2) H2Sep[b][Seq]->SetBinContent(xBinToFill, y+1, 0.3);
158  }
159  }
160 
161  Count++;
162  if (Count == nSep) Count = 0;
163  }//x, loop over channel index (valid channels only)
164 
165  H2_BS ->Delete();
166  H2_ADC->Delete();
167  }//b, DB or data
168 
169  //Draw
170  gStyle->SetOptDate(0);
171  gStyle->SetOptStat(0);
172  TCanvas *c1;
173  int iCVS = 0;
174  int iPAD = 1;
175  for (int c=0; c<nSepHisto; c++)
176  {
177  if (iPAD==1)
178  {
179  c1 = new TCanvas(Form("FmsBsQa_d%i_%i", a+8,iCVS), Form("detId=%i, %i",a+8,iCVS), 1600, 900);
180  c1->Divide(2, 2);
181  }
182  c1->cd(iPAD+0); H2Sep[0][c]->DrawCopy("col");
183  c1->cd(iPAD+2); H2Sep[1][c]->DrawCopy("col");
184  iPAD++;
185  if (iPAD==3)
186  {
187  if (PRINT) c1->Print(Form("%s.png", c1->GetName()));
188  iCVS++;
189  iPAD=1;
190  }
191  }
192  }//a, detId
193  #endif
194 
195  //-------------------------------------------
196 
197  #if 10
198  //Make and Show mismatching (DB <-> data) summary plot
199  vector<int> ListMM1[4]; //Mismatch
200  vector<int> ListMM2[4]; //Mismatch, only due to dead channels
201  for (int a=0; a<4; a++)
202  {
203  //Get ADC histogram
204  TH2F* H2_ADC = (TH2F*)F->Get(Form("Adc_d%i", a+8));
205  const float nCh = (H2_ADC->GetNbinsX()>300)?394:238;
206  const float nHitAll = H2_ADC->GetEntries();
207  const float nHitAvg = nHitAll/nCh;
208 
209  //Get, Normalize, and Remove ambiguity of histograms
210  TH2F* H2[2]; //DB or data
211  for (int b=0; b<2; b++)
212  {
213  H2[b] = (TH2F*)F->Get(Form("BitShift_%s_d%i", b==0?"DB":"data", a+8));
214  if (!H2[b]) { cout <<"Cannot open histograms! (2)" <<endl; return; }
215  NormalizeBsHisto(H2[b]);
216  }//b, DB or data
217 
218  //Listup mismatching channels
219  const int nBit = H2[db]->GetNbinsY();
220  for (unsigned int x=0; x<nToCh[a].size(); x++)
221  {
222  const int ch = nToCh[a][x];
223  bool FlagMM = false;
224 
225  //Check in "DB-wise" point of view
226  for (int y=0; y<nBit; y++)
227  {
228  const int BC_DB = H2[db] ->GetBinContent(ch, y+1); if (BC_DB == 0) continue;
229  const int BC_data = H2[data]->GetBinContent(ch, y+1);
230  if (BC_DB != BC_data)
231  {
232  ListMM1[a].push_back(ch);
233  FlagMM = true;
234  break;
235  }
236  }
237  if (FlagMM == true) continue;
238 
239  int nBitFilled = 0;
240  for (int y=0; y<nBit; y++) { if (H2[data]->GetBinContent(ch, y+1) != 0) nBitFilled++; }
241 
242  //Check if this channel is good or not: 0 for good, 1 for dead, and 2 for bad
243  int chQuality = 0;
244  if (nBitFilled == nBit) chQuality = 1;
245  else if (abs(nBitFilled) > 5)
246  {
247  TH1F* H1_ADC_1ch = (TH1F*)H2_ADC->ProjectionY("", nToCh[a][x], nToCh[a][x]);
248  if (H1_ADC_1ch->GetEntries() < nHitAvg*1.e-2) chQuality = 1;
249  else if (H1_ADC_1ch->GetRMS() < 1.) chQuality = 2;
250  H1_ADC_1ch->Delete();
251  }
252 
253  //Check in data-wise point of view
254  for (int y=0; y<nBit; y++) //Ascending order: targets positive bit shifted case
255  {
256  const int BC_data = H2[data]->GetBinContent(ch, y+1);
257  if (BC_data != 0)
258  {
259  const int BC_DB = H2[db]->GetBinContent(ch, y+1);
260  if (BC_data != BC_DB)
261  {
262  if (chQuality == 0) ListMM1[a].push_back(ch);
263  else ListMM2[a].push_back(ch);
264  break;
265  }
266  }
267  else break; //Met empty bin (corresponding ADc exist)
268  }
269  }//x, Loop over ch index (valid channels only): listup finished
270 
271  H2_ADC->Delete();
272  H2[db]->Delete();
273  H2[data]->Delete();
274  }//a, detId
275 
276  //Fill QA histograms and Draw
277  TH2F* H2_Mismatch[2][2];
278  for (int a=0; a<2; a++) //DB or data
279  for (int b=0; b<2; b++) //Pure mismatch or by dead channel
280  {
281  const char* tempN = Form("H2_MisMatch_%s_%i", a==0?"DB":"data", b);
282  const char* tempT = Form("%s, %s;;Bit (BS)", a==0?"DB":"data",b==0?"Mismatch":"Mismatch due to bad channel");
283  unsigned int nMM;
284  if (b==0) nMM = ListMM1[0].size() + ListMM1[1].size() + ListMM1[2].size() + ListMM1[3].size();
285  else nMM = ListMM2[0].size() + ListMM2[1].size() + ListMM2[2].size() + ListMM2[3].size();
286  H2_Mismatch[a][b] = new TH2F(tempN, tempT, (nMM>30)?nMM:30,0,(nMM>30)?nMM:30, 12,0.5,12+0.5);
287  for (int x=0; x<H2_Mismatch[a][b]->GetNbinsX(); x++) H2_Mismatch[a][b]->GetXaxis()->SetBinLabel(x+1, "");
288  for (int y=0; y<H2_Mismatch[a][b]->GetNbinsY(); y++)
289  {
290  if (y+1 <= 6) H2_Mismatch[a][b]->GetYaxis()->SetBinLabel(y+1, Form("%2i (#plus%i)", y+1, y+1));
291  else H2_Mismatch[a][b]->GetYaxis()->SetBinLabel(y+1, Form("%2i (#minus%i)", y+1, abs(y-12)));
292  }
293  H2_Mismatch[a][b]->LabelsOption("v", "X");
294  H2_Mismatch[a][b]->GetYaxis()->SetLabelSize(0.045);
295  H2_Mismatch[a][b]->GetYaxis()->SetTitleOffset(1.2);
296  }//a, b
297 
298  int iMM1 = 1;
299  int iMM2 = 1;
300  for (int a=0; a<4; a++)
301  {
302  //Get, Normalize, and Remove ambiguity of histograms
303  TH2F* H2[2]; //DB or data
304  for (int b=0; b<2; b++)
305  {
306  H2[b] = (TH2F*)F->Get(Form("BitShift_%s_d%i", b==0?"DB":"data", a+8));
307  if (!H2[b]) { cout <<"Cannot open histograms! (2)" <<endl; return; }
308 
309  //Nomarlize + Remove ambiguity
310  const float tempMax = H2[b]->GetMaximum();
311  H2[b]->Scale(1/tempMax);
312  for (int x=0; x<H2[b]->GetNbinsX(); x++)
313  for (int y=0; y<H2[b]->GetNbinsY(); y++)
314  {
315  const float BC = H2[b]->GetBinContent(x+1, y+1);
316  if (BC < 1.0)
317  {
318  H2[b]->SetBinContent(x+1, y+1, 0);
319  H2[b]->SetBinError (x+1, y+1, 0);
320  }
321  }
322  }//b, DB or data
323 
324  //Fill
325  const int nBit = H2[db]->GetNbinsY();
326  for (unsigned int x=0; x<ListMM1[a].size(); x++)
327  {
328  const int ch = ListMM1[a][x];
329  for (int y=0; y<nBit; y++)
330  {
331  for (int b=0; b<2; b++) //DB or data
332  {
333  const int BC = H2[b]->GetBinContent(ch, y+1);
334  H2_Mismatch[b][0]->SetBinContent(iMM1, y+1, BC);
335  H2_Mismatch[b][0]->GetXaxis()->SetBinLabel(iMM1, Form("%i_%i", a+8, ch));
336  }
337  }
338  iMM1++;
339  }//x, ch
340  for (unsigned int x=0; x<ListMM2[a].size(); x++)
341  {
342  const int ch = ListMM2[a][x];
343  for (int y=0; y<nBit; y++)
344  {
345  for (int b=0; b<2; b++) //DB or data
346  {
347  const int BC = H2[b]->GetBinContent(ch, y+1);
348  H2_Mismatch[b][1]->SetBinContent(iMM2, y+1, BC);
349  H2_Mismatch[b][1]->GetXaxis()->SetBinLabel(iMM2, Form("%i_%i", a+8, ch));
350  }
351  }
352  iMM2++;
353  }//x, ch
354 
355  H2[db] ->Delete();
356  H2[data]->Delete();
357  }//a, detId
358 
359  gStyle->SetOptDate(0);
360  gStyle->SetOptStat(0);
361  TCanvas *c2 = new TCanvas("FmsBsQa_Summary", "Summary", 1600, 900);
362  c2->Divide(2, 2);
363  for (int a=0; a<2; a++) //DB or data
364  for (int b=0; b<2; b++)
365  {
366  c2->cd(2*a + b + 1);
367  H2_Mismatch[a][b]->DrawCopy("col");
368  }
369  if (PRINT) c2->Print(Form("%s.png", c2->GetName()));
370  #endif
371 
372  return;
373 }//Main
Generated by   doxygen 1.8.5