StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
bsQa.C
1 #include "TCanvas.h"
2 #include "TF1.h"
3 #include "TFile.h"
4 #include "TGraph.h"
5 #include "TGraphErrors.h"
6 #include "TH1.h"
7 #include "TH2.h"
8 #include "TMath.h"
9 #include "TLegend.h"
10 #include "TLine.h"
11 #include "TPaveText.h"
12 #include "TString.h"
13 #include "TSystem.h"
14 #include "TStyle.h"
15 
16 #include <cmath>
17 #include <fstream>
18 #include <iostream>
19 #include <map>
20 using namespace std;
21 
22 //=============================
23 void NormalizeBsHisto(TH2F* H2)
24 {
25  //Normalize bit shift QA histogram for the case of it merged from multiple files
26  const float max = H2->GetMaximum();
27  H2->Scale(1/max);
28 
29  //Remove ambigue bins: only bins with value 1 is valid
30  for (int x=0; x<H2->GetNbinsX(); x++)
31  for (int y=0; y<H2->GetNbinsY(); y++)
32  {
33  const float BC = H2->GetBinContent(x+1, y+1);
34  if (BC < 1.0)
35  {
36  H2->SetBinContent(x+1, y+1, 0);
37  H2->SetBinError (x+1, y+1, 0);
38  }
39  }
40 
41  return;
42 }//NormalizeBsHisto
43 
44 //===================================================================================================================
45 vector<int> GetListBad(TH2F* ADC, TH2F* BS, TH2F* ChMap, int detId, int runNo, bool SHOW = false, bool PRINT = false)
46 {
47  vector<int> List;
48  if (BS->GetMaximum() != 1) { cout <<"WARNING! BS QA histogram is NOT normalized!" <<endl; return List; }
49 
50  //Get average # of hit
51  const float nCh = (ADC->GetNbinsX()>300)?394:238;
52  const float nHitAll = ADC->GetEntries();
53  const float nHitAvg = nHitAll/nCh;
54 
55  //Bad channels decision
56  for (int x=0; x<BS->GetNbinsX(); x++) //Bin # = channel
57  {
58  //0. Good (default)
59  //1. Dead: # of occupied bit (nBitF) = 12 or # of entries < 1 % of average # of hit
60  //2. Bad: ADC RMS < RMS cut (1.0, small RMS = ADCs are concentrated at certain #)
61  int flagStat = 0;
62 
63  //Invalid ch or Dead?
64  bool chValid = ChMap->GetBinContent(x+1, 1);
65  const int nBit = BS->GetNbinsY();
66  int nBitF = 0; //# of bit left filled
67  for (int y=0; y<nBit; y++) { if (BS->GetBinContent(x+1, y+1) == 1) nBitF++; }
68  if (chValid==false || nBitF==nBit)
69  {
70  flagStat = 1;
71  List.push_back(flagStat); //Regard it as dead
72  continue;
73  }
74 
75  //Valid + Some bits are occupied. Now check ADC distribution itself
76  TH1F* ADC_1ch = (TH1F*)ADC->ProjectionY("", x+1, x+1);
77  const float tHit = ADC_1ch->GetEntries();
78  const float tRMS = ADC_1ch->GetRMS();
79  if (tHit < nHitAvg*1.e-2) flagStat = 1;
80  else if (tRMS < 1.) flagStat = 2;
81  ADC_1ch->Delete();
82 
83  List.push_back(flagStat);
84  }//x, ch
85 
86  //Crosscheck functions
87  //-------------------------------------------
88 
89  if (SHOW)
90  {
91  for (unsigned int x=0; x<List.size(); x++)
92  {
93  if (List[x] == false) continue; //Good ch
94  if (ChMap->GetBinContent(x+1, 1) == false) continue; //Invalid ch
95  TH1F* ADC_1ch = (TH1F*)ADC->ProjectionY("", x+1, x+1);
96  const int nHit = ADC_1ch->GetEntries();
97  cout <<Form("Bad: detId = %2i, ch = %3i / Flag = %i", detId, x+1, List[x]) <<endl;
98  ADC_1ch->Delete();
99  }
100  }
101  if (PRINT) //Use this in batch mode
102  {
103  TCanvas *c1 = new TCanvas("c1", "", 800, 800);
104  for (unsigned int x=0; x<List.size(); x++)
105  {
106  if (List[x] == false) continue;
107  if (ChMap->GetBinContent(x+1, 1) == false) continue; //Invalid ch
108  c1->Clear();
109  c1->SetName(Form("Bad_run%i_d%i_ch%i", runNo, detId, x+1));
110  c1->cd();
111  TH1F* ADC_1ch = (TH1F*)ADC->ProjectionY("", x+1, x+1);
112  ADC_1ch->SetTitle(Form("Bad: run = %i detId = %2i, ch = %3i, flag = %i", runNo, detId, x+1, List[x]));
113  ADC_1ch->GetXaxis()->SetRangeUser(0,250);
114  ADC_1ch->DrawCopy("hist e");
115  c1->Print(Form("%s.png", c1->GetName()));
116  }
117  delete c1;
118  }
119 
120  return List;
121 }//GetListBad
122 
123 //========================================================================
124 vector<int> GetListBS(TH2F* BS, TH2F* ChMap, int detId, bool SHOW = false)
125 {
126  vector<int> List;
127  if (BS->GetMaximum() != 1) { cout <<"WARNING! BS QA histogram is NOT normalized!" <<endl; return List; }
128 
129  //A channel can have only 1 type of BS: positive, 0, or negative
130  const int nBit = BS->GetNbinsY();
131  for (int x=0; x<BS->GetNbinsX(); x++) //Bin # (x+1) = channel
132  {
133  int bsValue = 0;
134  bool chValid = ChMap->GetBinContent(x+1, 1);
135  if (!chValid) { List.push_back(bsValue); continue; }
136 
137  int nBitFilled = 0;
138  bool bsPos = false;
139  bool bsNeg = false;
140  for (int y=0; y<nBit; y++)
141  {
142  const int BC = BS->GetBinContent(x+1, y+1);
143  if (BC==1) nBitFilled++;
144 
145  if (y==0 && BC==1) bsPos = true;
146  if (y==nBit-1 && BC==1 && bsPos==false) bsNeg = true;
147  }
148 
149  if (nBitFilled!=0 || nBitFilled!=nBit) //Both valid and not dead
150  {
151  if (bsPos==true)
152  {
153  for (int y=0; y<nBit; y++)
154  {
155  if (BS->GetBinContent(x+1, y+1) == 1) bsValue++;
156  else break;
157  }
158  }
159  else if (bsNeg==true)
160  {
161  for (int y=nBit-1; y>-1; y--)
162  {
163  if (BS->GetBinContent(x+1, y+1) == 1) bsValue--;
164  else break;
165  }
166  }
167  }
168 
169  List.push_back(bsValue);
170  if (SHOW) cout <<Form("BS: %2i, %3i / %2i", detId, x+1, bsValue) <<endl;
171  }//x, ch
172 
173  return List;
174 }//GetListBS
175 
176 //=====================================================================
177 void DrawAdc(TH2F* H2, int runNo, int detId, int ch, bool PRINT = true)
178 {
179  TH1F* H1 = (TH1F*)H2->ProjectionY("", ch, ch);
180  H1->SetTitle(Form("d%i_ch%i_run%i", detId, ch, runNo));
181 
182  int xMax = 1;
183  for (int x=0; x<H1->GetNbinsX(); x++) { if (H1->GetBinContent(x+1) != 0) xMax = x+1; }
184  H1->GetXaxis()->SetRangeUser(0, xMax + 10);
185  H1->GetXaxis()->SetLabelSize(0.05);
186  H1->GetXaxis()->SetTitleOffset(1.25);
187  H1->GetYaxis()->SetLabelSize(0.05);
188 
189  gStyle->SetOptDate(0);
190  gStyle->SetOptStat("emr");
191  TCanvas *c1 = new TCanvas("c1", H1->GetTitle(), 800*1.5, 600*1.5);
192  c1->Divide(1, 2);
193  c1->cd(1)->SetLogy();
194  H1->DrawCopy("hist e");
195  c1->cd(2);
196  H1->GetXaxis()->SetRangeUser(0, 32);
197  H1->DrawCopy("hist e");
198  H1->Delete();
199 
200  if (PRINT) c1->Print(Form("%s.png", c1->GetTitle()));
201  delete c1;
202  return;
203 }//DrawAdc
204 
205 //Main
206 //==================================================================================
207 void bsQa(const char* inList = "inFiles_pptrans.list", int LogLv = 1, int PRINT = 0)
208 {
209  /*
210  LogLv:
211  0: No log file will be created
212  1: Show runs and responsible channels when new list being produced
213  2: Add mismatching (DB <-> data) channels to 1.
214  3: Add bad channels to 2.
215 
216  PRINT:
217  0: No updated list file will be printed out
218  1: BS list will be printed out
219  2: Bad channels list will be printed out in addition to 1.
220  */
221 
222  //Blacklists: skip master list update process for these channels
223  const int ListSkip[][4] =
224  {
225  //Add new row in following format to mark: run# (mark begins), run# (mark ends), detId, and ch
226  //{, , , }, //
227 
228  //RUN15 pp200trans
229  {16073030, 16073037, 8, 483}, //Small # of entries (62) at 16073030 caused accidental update
230  {16068056, 16082011, 8, 566}, //This channel is simply bad: but it keep escapes bad ch judgment routine
231  {16073032, 16073034, 10, 9}, //One or Two fucking entry at ADC ~490
232  {16078028, 16078042, 10, 13}, //Update at 16077021 (0 -> -1) seems reasonable
233  {16073034, 16077021, 10, 47}, //Update at 16073031 (-4 -> -3) seems reasonable
234  {16066035, 16073030, 10, 48}, //Bad channel
235  {16078032, 16078033, 10, 51}, //Adjacent run with only one channel fluctuating... skip it
236  {16066050, 16066050, 10, 65}, //No problem in algorithm, but it seems large ADC (1 entry) is accidental
237  {16077027, 16077040, 10, 73}, //One entry at ADC = 1324 caused update - it continually in and out
238  {16077032, 16093019, 10, 78}, //Heavily jump in and out. Cannot gaurantee quality of list
239  {16073030, 16073033, 10, 151}, //Small RMS (0.85 at run16073030) caused accidental update
240  {16073030, 16073037, 10, 176}, //Small RMS (0.88) at 16073030 caused accidental update
241  {16073037, 16073038, 10, 186}, //One ADC entry ~1,300 at run 16073037 caused update
242  {16073031, 16073040, 10, 282}, //No problem in algorithm but to reduce # of lists
243  {16073035, 16073040, 11, 118}, //Bad channel judgment caused accidental update
244  {16086050, 16086051, 11, 213}, //Adjacent run with only one channel fluctuating... skip it
245  {16067021, 16072024, 11, 287}, //ADC starts from 32 in this period, but from run 16072038 starts from 16
246  {16073031, 16077021, 11, 287}, //Marked as bad channnel: thus old DB value will overwrite. Make it skip
247  //run 16085006 removed from good run list: overall short statistics
248 
249  //RUN15 pAu
250  {16129050, 16129050, 10, 211}, //Typical jump channel, short statistics for this particular run
251  {16129050, 16129051, 10, 221}, //Typical jump channel, short statistics for this particular run
252  //run 16134016 removed from good run list: overall short statistics
253 
254  //RUN15 pAl
255  {16164006, 16167091, 10, 30}, //BS_DB at 16160018 is -1. But at later runs sometimes BS_data judged as 0
256  {16160020, 16169041, 10, 31}, //Actual BS is 2. Suspect in general small statistics cause continual jump
257  {16160023, 16169041, 10, 32}, //Actual BS is 4. Same type of cell to above ch 10_31
258  {16160018, 16162040, 10, 43}, //Actual BS is 2, but BS_DB is -1. Plus low statistics caused update
259  {16160024, 16169041, 10, 177} //Actual BS is 4. Statistics is plenty. Sometimes RMS of ADC is < 1
260  };
261  const int nListSkip = sizeof(ListSkip)/sizeof(ListSkip[0]);
262 
263  //-------------------------------------------
264 
265  //Master lists for continual update/printout
266  vector<int> ListBad[4];
267  vector<int> ListBS[4];
268 
269  TFile *F;
270  ifstream in;
271  in.open(inList);
272  if (!in.is_open()) { cout <<"Cannot open " <<inList <<endl; return; }
273  ofstream out;
274 
275  int iRun = 0;
276  int runNoPrevious = 0;
277  string inFile;
278  while (in.is_open())
279  {
280  //Open file
281  in >> inFile;
282  if (!in.good()) { break; in.close(); }
283  F = TFile::Open(inFile.c_str());
284  if (!F || F->IsZombie()) { cout <<"Cannot open the file " <<inFile.c_str() <<endl; return; }
285  cout <<Form("Processing... %s, %3i", F->GetName(), iRun) <<endl;
286 
287  //Get run number
288  std::size_t strPos = inFile.find(".root");
289  string runNoStr = inFile.substr(strPos-8, 8);
290  const int runNo = std::atoi(runNoStr.c_str());
291  if (LogLv>0)
292  {
293  if (iRun==0) out.open("log_fmsBsQa.txt");
294  else out.open("log_fmsBsQa.txt", ios::app);
295  out <<Form("Run %i", runNo) <<endl;
296  }
297 
298  //-------------------------------------------------
299 
300  //Switch flag if new values found (CAVEAT: this is independent of mismatch btw DB <-> data)
301  bool FlagNew_dead = false;
302  bool FlagNew_bs = false;
303  for (int a=0; a<4; a++) //Loop over detId
304  {
305  TH2F* H2_ADC = (TH2F*)F->Get(Form("Adc_d%i", a+8));
306  TH2F* H2_BS_DB = (TH2F*)F->Get(Form("BitShift_DB_d%i", a+8));
307  TH2F* H2_BS_data = (TH2F*)F->Get(Form("BitShift_data_d%i", a+8));
308  TH2F* H2_ChMap = (TH2F*)F->Get(Form("ChMap_d%i", a+8));
309  if (!H2_ADC || !H2_BS_DB || !H2_BS_data || !H2_ChMap) { cout <<"Cannot open histogram!" <<endl; return; }
310 
311  //Normalize and Remove ambiguity of BS QA histogram
312  NormalizeBsHisto(H2_BS_DB);
313  NormalizeBsHisto(H2_BS_data);
314 
315  //Get dead channels list, Compare with existing value, and Update old one if necessary
316  vector<int> tempBad = GetListBad(H2_ADC, H2_BS_data, H2_ChMap, a+8, runNo, false, false);
317  if (iRun!=0 && ListBad[a].size() != tempBad.size()) { cout <<"Size mismatch! (1)" <<endl; return; }
318  if (iRun==0 || !std::equal(ListBad[a].begin(), ListBad[a].end(), tempBad.begin()))
319  {
320  if (iRun!=0 && LogLv>2)
321  {
322  for (unsigned int x=0; x<ListBad[a].size(); x++)
323  {
324  bool dead_old = ListBad[a][x];
325  bool dead_new = tempBad[x];
326  if (dead_old != dead_new)
327  {
328  out <<Form("Updating dead: %2i, %3i, %i -> %i", a+8,x+1,dead_old,dead_new) <<endl;
329  }
330  }
331  }
332  ListBad[a].clear();
333  ListBad[a] = tempBad;
334  FlagNew_dead = true;
335  }
336 
337  //Get BS lists (DB, data), Compare them, and Update DB list if necessary
338  //Exception 1: skip if the channel is dead
339  //Exception 2: data BS is negative and its absolute value is larger than DB
340  vector<int> tempBS_db = GetListBS(H2_BS_DB, H2_ChMap, a+8, false);
341  vector<int> tempBS_data = GetListBS(H2_BS_data, H2_ChMap, a+8, false);
342  vector<int> tempBS_new;
343  if (tempBS_db.size() != tempBS_data.size()) { cout <<"Size mismatch! (2)" <<endl; return; }
344  for (unsigned int x=0; x<tempBS_db.size(); x++)
345  {
346  int bs_db = tempBS_db[x];
347  int bs_data = tempBS_data[x];
348  if (abs(bs_db)>12 || abs(bs_data)>12) { cout <<Form("Invalid BS found! %2i %2i\n", bs_db,bs_data); }
349  int bs_new = bs_db;
350 
351  if (bs_db!=bs_data && tempBad[x]==false) //Mismatch, NOT dead
352  {
353  //Judge which bs should is correct
354  if (bs_data >= 0) bs_new = bs_data;
355  else //bs_data < 0
356  {
357  if (bs_db >= 0) bs_new = 0; //Statistics issue, in this case bs_data < 0 is NOT real BS
358  else //bs_data < 0 and bs_db < 0
359  {
360  if (abs(bs_data) < abs(bs_db)) bs_new = bs_data;
361  }
362  }
363 
364  if (LogLv>1 && bs_db!=bs_new)
365  {
366  out <<Form("BS mismatch: %2i, %3i | data: %3i, DB: %2i -> %2i",
367  a+8,x+1,bs_data,bs_db,bs_new) <<endl;
368  }
369  }//Mismatch
370 
371  //Check if this channel is blacklisted: skip if true
372  bool FlagSkip = false;
373  for (int i=0; i<nListSkip; i++)
374  {
375  if (runNo>=ListSkip[i][0] && runNo<=ListSkip[i][1] &&
376  a+8==ListSkip[i][2] && x+1==(unsigned int)ListSkip[i][3]) FlagSkip = true;
377  }
378  if (FlagSkip==true) bs_new = ListBS[a][x];
379 
380  tempBS_new.push_back(bs_new);
381  }//Loop over channels
382 
383  //Update old BS (DB) list if difference found
384  if (iRun==0 || !std::equal(ListBS[a].begin(), ListBS[a].end(), tempBS_new.begin()))
385  {
386  if (iRun!=0 && LogLv>0)
387  {
388  for (unsigned int x=0; x<tempBS_new.size(); x++)
389  {
390  int bs_old = ListBS[a][x];
391  int bs_new = tempBS_new[x];
392  if (bs_old != bs_new)
393  {
394  out <<Form("Updating BS: %2i, %3i | %2i -> %2i", a+8,x+1,bs_old,bs_new) <<endl;
395  //DrawAdc(H2_ADC, runNo, a+8, x+1, true); //MUST be used in batch mode
396  }
397  }
398  }
399  ListBS[a].clear();
400  ListBS[a] = tempBS_new;
401  FlagNew_bs = true;
402  }
403 
404  H2_ADC->Delete();
405  H2_BS_DB->Delete();
406  H2_BS_data->Delete();
407  }//a, detId
408  F->Close();
409  iRun++;
410  if (LogLv>0) { out <<endl; out.close(); }
411 
412  //-------------------------------------------------
413 
414  //Printout updated DB values
415  if (PRINT>1 && FlagNew_dead)
416  {
417  cout <<Form("List updated: run %i -> run %i, Dead channels", runNoPrevious, runNo) <<endl;
418  out.open(Form("FmsDead_y2015_run%i.txt", runNo));
419  for (int a=0; a<4; a++)
420  for (unsigned int b=0; b<ListBad[a].size(); b++)
421  {
422  out <<Form("%2i %3i %i", a+8, b+1, ListBad[a][b]) <<endl;
423  }
424  out.close();
425  }
426  if (PRINT>0 && FlagNew_bs)
427  {
428  cout <<Form("List updated: run %i -> run %i, Bit shift", runNoPrevious, runNo) <<endl;
429  out.open(Form("FmsBsGain_y2015_run%i.txt", runNo));
430  for (int a=0; a<4; a++)
431  for (unsigned int b=0; b<ListBS[a].size(); b++)
432  {
433  out <<Form("%2i %3i %2i", a+8, b+1, ListBS[a][b]) <<endl;
434  }
435  out.close();
436  }
437 
438  runNoPrevious = runNo;
439  }//Loop over runs
440 
441  return;
442 }//Main
Generated by   doxygen 1.8.5