BuildMetricsTree_OnMon.C
Go to the documentation of this file.
1 #include <iostream>
2 #include <fstream>
3 #include <iomanip>
4 #include <sstream>
5 #include <string>
6 #include <streambuf>
7 #include <cstring>
8 #include <algorithm>
9 
10 #include <TStyle.h>
11 #include <TSystem.h>
12 #include <TFile.h>
13 #include <TString.h>
14 #include <TTree.h>
15 #include <TH2F.h>
16 #include <TH1F.h>
17 #include <TMath.h>
18 
19 #include <TSQLServer.h>
20 #include <TSQLResult.h>
21 #include <TSQLStatement.h>
22 
23 using namespace std;
24 
25 struct RunInfo{
26 
27  int run;
28  bool isLast;
30  vector< double > livetime;
31  vector< int > subrun;
32  vector< TString > fnames;
33  vector< TString > rnames;
34  vector< bool > isNewRelom;
35  vector< bool > isNewRelana;
36 
37 };
38 
39 double GetMedian(TH1D* hist);
40 
41 float GetMedian(vector<float> scores);
42 
43 bool GetDBPartition(int run, int &par);
44 
45 int GetFileName(TString& filename, int run, int subrun, bool isOnMon = true);
46 
48 
49 vector< RunInfo > GetSubrunSets(RunInfo wholerun, double run_len);
50 
52 
53  gStyle->SetTimeOffset(0);
54 
55  Int_t apply_runlen = 1; //1 if want to apply live time cut; 0 if do not
56  Int_t apply_empty = 1; //1 if want to apply empty cut; 0 if do not
57  Int_t apply_time = 1; //1 if want to apply negative time cut; 0 if do not (start time before end time)
58  Int_t apply_db = 1; //1 if want to apply diblock cut; 0 if do not
59  Int_t apply_hits = 1; //1 if want to apply hit rate cut; 0 if do not
60  Int_t apply_reco = 1; //1 if want to apply reco cut; 0 if do not
61 
62  //Cut Parameter Values
63 
64  Double_t run_len = 1.00; //Minimum live time of run in seconds.
65  Double_t hits_up = 23.0; //Median MIP rate smaller than this parameter
66  Double_t hits_down = 13.0; //Median MIP rate larger than this parameter
67  Double_t PixelsPerFEB = 26; //Minimum # of good pixels per FEB
68  Double_t FEBsPerDCM = 56; //Minimum # of good FEBs per DCM
69  Double_t DCMsPerDB = 12; //Minimum # of good DCMs per DiBlock
70  Double_t MaxTopSideAsym = 0.5; //Maximum Top/Side MIP hit rate Asymmetry
71  Double_t MinTopSideAsym = 0.1; //Minimum Top/Side MIP hit rate Asymmetry
72  Double_t MinGoodDB = 2; //Minimum # of good consecutive DiBlocks
73  Double_t empty_cut = 100.0; //The percentage of empty events less than this value. This only cuts out files with ONLY empty events.
74  Double_t MaxTrkFrac2D = 0.15; //Maximum 2D track fraction
75  Double_t MinNumSlc = 1.2; //Minimum # of slices / trigger / 10^4 channels
76  Double_t MaxNumSlc = 3.2; //Maximum # of slices / trigger / 10^4 channels
77 
78  // MIP rate cuts. Entries are for DCMs {01 - 06 (Top), 07, 08, 09, 10, 11, 12}
79 
80  float lowpixmip[7] = {13, 5, 4, 4, 2, 2, 1.5}; //Minimum MIP rate for pixels in different DCMs
81  float highpixmip[7] = {31, 45, 36, 30, 26, 23, 20.}; //Maximum MIP rate for pixels in different DCMs
82 
83  float lowfebmip[7] = {400, 300, 200, 200, 100, 100, 100}; //Minimum MIP rate for FEBs in different DCMs
84  float highfebmip[7] = {1000, 1200, 1000, 850, 700, 600, 550}; //Maximum MIP rate for FEBs in different DCMs
85 
86 
87  // **************END OF TYPICAL USER'S NEED TO ALTER THIS SCRIPT*******************************************
88 
89  // Open file to keep track of processed subruns
90  ofstream procfile;
91  procfile.open("listProcessedTemp.txt");
92 
93  // Variables to fill Good Run metrics
97 
100 
102  Double_t numslc, trkfrac2D;
103 
104  Int_t preliminary;
105 
106  // Create root file to store tree
107  TFile *outfile = new TFile("MetricsNew.root","recreate");
108 
109  // Create tree to store Good Run metrics
110  TTree *mytree = new TTree("evts","");
111  mytree->Branch("run",&run,"run/I");
112  mytree->Branch("subrun",&subrun,"subrun/I");
113  mytree->Branch("par",&par,"par/I");
114  mytree->Branch("nevents",&nevents,"nevents/I");
115  mytree->Branch("firstsec",&firstsec,"firstsec/I");
116  mytree->Branch("lastsec",&lastsec,"lastsec/I");
117  mytree->Branch("ngoodpix",&ngoodpix,"ngoodpix/I");
118  mytree->Branch("ngoodfeb",&ngoodfeb,"ngoodfeb/I");
119  mytree->Branch("ngooddcm",&ngooddcm,"ngooddcm/I");
120  mytree->Branch("ngooddb",&ngooddb,"ngooddb/I");
121  mytree->Branch("ngoodmip",&ngoodmip,"ngoodmip/I");
122  mytree->Branch("ngoodcdb",&ngoodcdb,"ngoodcdb/I");
123  mytree->Branch("nactivefeb",&nactivefeb,"nactivefeb/I");
124  mytree->Branch("nactivedcm",&nactivedcm,"nactivedcm/I");
125  mytree->Branch("nactivedb",&nactivedb,"nactivedb/I");
126  mytree->Branch("pass_runlen",&pass_runlen,"pass_runlen/I");
127  mytree->Branch("pass_empty",&pass_empty,"pass_empty/I");
128  mytree->Branch("pass_hits",&pass_hits,"pass_hits/I");
129  mytree->Branch("pass_db",&pass_db,"pass_db/I");
130  mytree->Branch("pass_time",&pass_time,"pass_time/I");
131  mytree->Branch("pass_all",&pass_all,"pass_all/I");
132  mytree->Branch("setsize",&setsize,"setsize/I");
133  mytree->Branch("rectimesec",&rectimesec,"rectimesec/D");
134  mytree->Branch("setlivetime",&setlivetime,"setlivetime/D");
135  mytree->Branch("hitrate",&hitrate,"hitrate/D");
136  mytree->Branch("midhitrate",&midhitrate,"midhitrate/D");
137  mytree->Branch("miprate",&miprate,"miprate/D");
138  mytree->Branch("midmiprate",&midmiprate,"midmiprate/D");
139  mytree->Branch("mipratio",&mipratio,"mipratio/D");
140  mytree->Branch("mipasym",&mipasym,"mipasym/D");
141  mytree->Branch("nactivechannels",&nactivechannels,"nactivechannels/D");
142  mytree->Branch("emptypercentage",&emptypercentage,"emptypercentage/D");
143 
144  mytree->Branch("dbencoded",&dbencoded,"dbencoded/I");
145  mytree->Branch("dbaencoded",&dbaencoded,"dbaencoded/I");
146 
147  mytree->Branch("pass_reco",&pass_reco,"pass_reco/I");
148  mytree->Branch("pass_slc",&pass_slc,"pass_slc/I");
149  mytree->Branch("pass_trk",&pass_trk,"pass_trk/I");
150  mytree->Branch("corrupted",&corrupted,"corrupted/I");
151  mytree->Branch("procsec",&procsec,"procsec/I");
152  mytree->Branch("numslc",&numslc,"numslc/D");
153  mytree->Branch("trkfrac2D",&trkfrac2D,"trkfrac2D/D");
154 
155  mytree->Branch("preliminary",&preliminary,"preliminary/I");
156 
157  // Open list of subruns to process
158  ifstream fileListItr(filelist.Data());
159 
160  // Loop over subruns
161  while(fileListItr >> run){
162 
163  // Get a vector of subrun sets
164  vector< RunInfo > setvec = GetSubrunSets( GetRunInfo(run), run_len );
165 
166  // If vector is empty, fill tree with empty set
167  if(setvec.size()==0){
168 
169  subrun = 0; par = 4; nevents = 0; pass_runlen = 0; pass_hits = 0;
170  pass_all = 0; pass_db = 0; pass_time = 0; pass_empty = 0;
171  firstsec = 0; lastsec = 0; ngoodpix = 0; ngoodfeb = 0;
172  ngooddcm = 0; ngooddb = 0; ngoodmip = 0; nactivefeb = 0;
173  nactivedcm = 0; nactivedb = 0; ngoodcdb = 0; dbencoded = 0;
174  dbaencoded = 0; setsize = 0;
175 
176  hitrate = 0; nactivechannels = 0; rectimesec = 0; emptypercentage = 0;
177  setlivetime = 0; mipratio = 0; midhitrate = 0; midmiprate = 0;
178  mipasym = 0; miprate =0;
179 
180  pass_reco = 0; pass_slc = 0; pass_trk = 0; corrupted = 0; procsec = 0;
181  numslc = 0; trkfrac2D = 0;
182 
183  preliminary = 0;
184 
185  mytree->Fill();
186 
187  }
188 
189  // Loop over subrun sets
190  for(int seti=0; seti<int(setvec.size()); seti++){
191 
192  // 2D vector of Pixel MIP hits maps (DB x DCM)
193  vector< vector< TH2F* > > TotPixMapMip(14, vector< TH2F* >(12,0));
194 
195  // Reset live time of set
196  setlivetime = 0;
197 
198  // First loop over subruns in set to add Mip hits
199  for(int subruni=0; subruni<int(setvec[seti].subrun.size()); subruni++){
200 
201  // Open OnMon file
202  TFile omFile(setvec[seti].fnames[subruni],"read");
203 
204  // Add MIP hits for whole run
205  for(int db=1;db<=14;db++){
206  for(int dcm=1;dcm<=12;dcm++){
207 
208  // Temporary histogram to store maps
209  TH2F *h2temp = 0;
210 
211  // Get map from OnMon file
212  h2temp = (TH2F*)omFile.Get(TString::Format("MipADCPixelsDCM_%02d_%02d",db,dcm));
213 
214  // Check that map exists in file, otherwise continue
215  if(!h2temp){
216  delete h2temp;
217  continue;
218  }
219 
220  // If it's the first map, clone it, else add histograms
221  if(!TotPixMapMip[db-1][dcm-1]){
222  TotPixMapMip[db-1][dcm-1] = (TH2F*)h2temp->Clone("");
223  }
224  else{
225  TotPixMapMip[db-1][dcm-1]->Add(h2temp);
226  }
227 
228  // Delete temporary histogram
229  delete h2temp;
230 
231  }}
232 
233  // Close OnMon File
234  omFile.Clear();
235  omFile.Close();
236 
237  // Add the live time for this subrun
238  setlivetime += setvec[seti].livetime[subruni];
239 
240  }// End of first subruns loop
241 
242  // Fill run and preliminary status
243  run = setvec[seti].run;
244  preliminary = setvec[seti].isIncomplete;
245 
246  // Now loop over subruns again to apply cuts
247  for(int subruni=0; subruni<int(setvec[seti].subrun.size()); subruni++){
248 
249  // Fill subrun number and set size
250  subrun = setvec[seti].subrun[subruni];
251  setsize = setvec[seti].subrun.size();
252 
253  // Fill single subrun live time
254  rectimesec = setvec[seti].livetime[subruni];
255 
256  bool isNewRelom = setvec[seti].isNewRelom[subruni];
257 
258  // Get Header information
259  double startHour, endHour;
260  Int_t Nevents;
261  UInt_t Partition;
262  UInt_t startYear, endYear;
263  UShort_t startMonth, endMonth;
264  UShort_t startDay, endDay;
265 
266  TFile omFile(setvec[seti].fnames[subruni],"read");
267 
268  TTree *header = (TTree*)omFile.Get("Header");
269  header->SetBranchAddress("Nevents",&Nevents);
270  header->SetBranchAddress("StartHour",&startHour);
271  header->SetBranchAddress("StartYear",&startYear);
272  header->SetBranchAddress("StartMonth",&startMonth);
273  header->SetBranchAddress("StartDay",&startDay);
274  header->SetBranchAddress("EndHour",&endHour);
275  header->SetBranchAddress("EndYear",&endYear);
276  header->SetBranchAddress("EndMonth",&endMonth);
277  header->SetBranchAddress("EndDay",&endDay);
278  // Partition info only exists in newer releases
279  if(isNewRelom) header->SetBranchAddress("Partition",&Partition);
280 
281  header->GetEntry(0);
282 
283  // Number of triggers
284  nevents = Nevents;
285 
286  // Get times in seconds. Need to think about time zones.
287  double hour = startHour;
288  int hourInt = int(hour);
289  int minutes = int((hour - hourInt)*60);
290  int seconds = int(((hour-hourInt)*60 - minutes)*60);
291  TDatime startDate;
292  if(isNewRelom) startDate.Set(startYear,startMonth,startDay,hourInt,minutes,seconds);
293  else startDate.Set(startYear+1900,startMonth+1,startDay,hourInt,minutes,seconds);
294 
295  hour = endHour;
296  hourInt = int(hour);
297  minutes = int((hour - hourInt)*60);
298  seconds = int(((hour-hourInt)*60 - minutes)*60);
299  TDatime endDate;
300  if(isNewRelom) endDate.Set(endYear,endMonth,endDay,hourInt,minutes,seconds);
301  else endDate.Set(endYear+1900,endMonth+1,endDay,hourInt,minutes,seconds);
302 
303  firstsec = startDate.Convert();
304  lastsec = endDate.Convert();
305 
306  delete header;
307 
308  bool connected = false;
309 
310  // Get partition number
311  if(isNewRelom){
312  par = Partition;
313  connected = true;
314  }
315  // Old releases need to use the database
316  else{
317  connected = GetDBPartition(run,par);
318  }
319 
320  // Skip if database connection failed
321  if(!connected) continue;
322 
323  // Get number of empty events
324  TH2F *NhitVsHour = (TH2F*)omFile.Get("NhitVsHour");
325 
326  int nemptyevents = NhitVsHour->Integral(61,1500,1,1);
327  int nallevents = NhitVsHour->Integral(61,1500);
328  if(nallevents!=0) emptypercentage = 100.0*nemptyevents/nallevents;
329  else emptypercentage = 0;
330 
331  delete NhitVsHour;
332 
333  // Initialize variables
334  nactivechannels = 0;
335  hitrate = 0;
336 
337  nactivefeb = 0; nactivedb = 0; nactivedcm = 0; ngoodcdb = 0;
338  ngoodpix = 0; ngoodfeb = 0; ngooddb = 0; ngooddcm = 0; ngoodmip = 0;
339  dbencoded = 0; dbaencoded = 0;
340 
341  // Vectors to store all pixel hits
342  vector<float> allMipHitsTop;
343  vector<float> allMipHitsSide;
344  vector<float> allHitRates;
345 
346  // Nested loops to determine good diblock cut
347 
348  // This is to keep track of good consecutive diblocks
349  int tempdb = 0;
350 
351  // Loop over diblocks
352  for(int db=1;db<=14;db++){
353 
354  // Count good dcms in this db
355  int nGoodDCMInDB = 0;
356 
357  // Vectors to store pixel hits for this db
358  vector<float> pixMipHitsTop;
359  vector<float> pixMipHitsSide;
360 
361  // Keep track of active dcms
362  int tempdcm = nactivedcm;
363 
364  // Loop over DCMs
365  for(int dcm=1;dcm<=12;dcm++){
366 
367  // Check if MIP hit map exists for this DCM
368  if(!TotPixMapMip[db-1][dcm-1]) continue;
369 
370  // Get overall hit rate map for this DCM
371  TH2F *PixelRate = (TH2F*)omFile.Get(TString::Format("PixelsRateDCM_%02d_%02d",db,dcm));
372  if(!PixelRate) continue; // Skip if it doesn't exist
373 
374  // Count good FEBs in this DCM
375  int nGoodFEBInDCM = 0;
376 
377  // Keep track of active FEBs
378  int tempfeb = nactivefeb;
379 
380  // Loop over FEBs
381  for(int feb=0;feb<64;feb++){
382 
383  // Count good Pixels in this FEB
384  int nGoodPixInFEB = 0;
385 
386  // Keep track of active Pixels
387  int temppix = nactivechannels;
388 
389  // Sum of MIP hits in this FEB
390  double febmiphits = 0;
391 
392  // Loop over Pixels
393  for(int pix=0;pix<32;pix++){
394 
395  // Figure out bins x and y for this pixel
396  int ix = (feb%16)*4 + pix%4 + 1;
397  int iy = int(feb/16)*8 + int(pix/4) + 1;
398 
399  // Check this pixel has MIP hits, i.e. it is active
400  if( TotPixMapMip[db-1][dcm-1]->GetBinContent(ix,iy)==0 ) continue;
401 
402  // Increment active channels
403  nactivechannels++;
404 
405  // Store overall hit rate and MIP hits in this Pixel
406  double pixhitrate = PixelRate->GetBinContent(ix,iy);
407  double pixmiphits = TotPixMapMip[db-1][dcm-1]->GetBinContent(ix,iy);
408 
409  // Increment MIP hits in this FEB with a max threshold of 4x the pixel high MIP rate cut
410  // The threshold mitigates the effect of few pixels with very large rates
411  febmiphits += TMath::Min(pixmiphits, highpixmip[TMath::Max(dcm-6,0)] * setlivetime * 4);
412 
413  // Fill total hit rate variable and vectors
414  hitrate += pixhitrate;
415  allHitRates.push_back(pixhitrate);
416  if(dcm<7) pixMipHitsTop.push_back(pixmiphits);
417  else pixMipHitsSide.push_back(pixmiphits);
418 
419  // Determine if pixel is good, i.e. neither cold nor hot
420  if(pixhitrate<pow(10,3.5) && pixhitrate>pow(10,0.5) &&
421  TMath::Gamma(lowpixmip[TMath::Max(dcm-6,0)] * setlivetime , pixmiphits) > 0.1 &&
422  TMath::Gamma(highpixmip[TMath::Max(dcm-6,0)] * setlivetime , pixmiphits) < 0.9 ){
423 
424  // Increment good pixels in total and in this FEB
425  ngoodpix++;
426  nGoodPixInFEB++;
427 
428  }
429 
430  }// End of Pixel loop
431 
432  // If active channels has increased, this FEB is active.
433  // Increment active FEBs
434  if(temppix < nactivechannels) nactivefeb++;
435 
436  // Determine if FEB is good, i.e. good MIP rate and enough good pixels
437  if(nGoodPixInFEB >= PixelsPerFEB &&
438  TMath::Gamma(lowfebmip[TMath::Max(dcm-6,0)] * setlivetime , febmiphits) > 0.1 &&
439  TMath::Gamma(highfebmip[TMath::Max(dcm-6,0)] * setlivetime , febmiphits) < 0.9 ){
440 
441  // Increment good FEBs in total and in this DCM
442  ngoodfeb++;
443  nGoodFEBInDCM++;
444 
445  }
446 
447  }// End of FEB loop
448 
449  // If active FEBs has increased, this DCM is active.
450  // Increment active DCMs
451  if(tempfeb < nactivefeb) nactivedcm++;
452 
453  // Determine if DCM is good, i.e. enough good FEBs
454  if(nGoodFEBInDCM >= FEBsPerDCM){
455  // Increment good DCMs in total and in this diblock
456  ngooddcm++;
457  nGoodDCMInDB++;
458  }
459 
460  // Delete pixel rate map
461  delete PixelRate;
462 
463  }// End of DCM loop
464 
465  // If active DCMs has increased, this diblock is active.
466  // Increment active diblocks
467  if(tempdcm < nactivedcm){
468  nactivedb++;
469  // Also, fill active diblock mask
470  dbaencoded += pow(2,db-1);
471  }
472 
473  // Fill vectors for whole detector with hits from this diblock
474  allMipHitsTop.insert(allMipHitsTop.end(), pixMipHitsTop.begin(), pixMipHitsTop.end());
475  allMipHitsSide.insert(allMipHitsSide.end(), pixMipHitsSide.begin(), pixMipHitsSide.end());
476 
477  // Find median of MIP hits for top and side of this diblock
478  double medianTop = GetMedian(pixMipHitsTop);
479  double medianSide = GetMedian(pixMipHitsSide);
480 
481  // Calculate top-side asymmetry if there are MIP hits
482  if(medianTop + medianSide != 0) medianTop = (medianTop - medianSide) / (medianTop + medianSide);
483  else medianTop = 0;
484 
485  // Determine if diblock is good
486  if(nGoodDCMInDB >= DCMsPerDB){
487  // Increment good diblocks (DCMs statuses only)
488  ngooddb++;
489  // Check for bad top-side asymmetry
490  if(medianTop < MaxTopSideAsym &&
491  medianTop > MinTopSideAsym ){
492  // Increment good diblocks (including asymmetry)
493  ngoodmip++;
494  // Counting consecutive good diblocks
495  tempdb++;
496  // Fill good diblock mask
497  dbencoded += pow(2,db-1);
498  }
499  // Else, reset consecutive diblocks
500  else tempdb = 0;
501  }
502  // Else, reset consecutive diblocks
503  else tempdb = 0;
504 
505  // Fill largest number of consecutive diblocks
506  if(tempdb > ngoodcdb) ngoodcdb = tempdb;
507 
508  }// End of diblock loop
509 
510 
511  // Compute median of all MIP hits for top and side
512  double medianAllTop = GetMedian(allMipHitsTop);
513  double medianAllSide = GetMedian(allMipHitsSide);
514 
515  // Calculate overall top-side asymmetry
516  if(medianAllTop + medianAllSide != 0) mipasym = (medianAllTop - medianAllSide) / (medianAllTop + medianAllSide);
517  else mipasym = 0;
518 
519  // Calculate top-side ratio
520  if(medianAllSide!=0) mipratio = medianAllTop / medianAllSide;
521  else mipratio = -1;
522 
523  // Normalise hit rate per active channel (mean)
524  if(nactivechannels!=0) hitrate /= nactivechannels;
525  else hitrate = 0;
526 
527  // Get overall median of hit rates
528  midhitrate = GetMedian(allHitRates);
529 
530  // Add top and side MIP hits and compute median
531  allMipHitsTop.insert(allMipHitsTop.end(), allMipHitsSide.begin(), allMipHitsSide.end());
532  midmiprate = GetMedian(allMipHitsTop);
533  if(setlivetime != 0) midmiprate /= setlivetime;
534  else midmiprate = 0;
535 
536  // Reset mean MIP rate
537  miprate = 0;
538 
539  // Add all MIP hits
540  for(int i=0;i<int(allMipHitsTop.size());i++){
541  miprate += allMipHitsTop[i];
542  }
543 
544  // Compute mean MIP rate
545  double chtime = allMipHitsTop.size() * setlivetime;
546  if(chtime!=0) miprate /= chtime;
547  else miprate = 0;
548 
549  // Close OnMon file
550  omFile.Clear();
551  omFile.Close();
552 
553 
554  // Start looking at reconstructed variables
555 
556  // Initialize variables. If no nearline file exists, subrun passes reco cut.
557  pass_reco = 1;
558  pass_trk = 1;
559  pass_slc = 1;
560  corrupted = 0;
561  numslc = -1;
562  trkfrac2D = -1;
563 
564  // Get file name and release type
565  TString anaFileName = setvec[seti].rnames[subruni];
566  bool isNewRelana = setvec[seti].isNewRelana[subruni];
567 
568  // Check if file exists
569  bool hasReco = true;
570  if(anaFileName == "") hasReco = false;
571 
572  // Fill timestamp in which this tree was filled
573  TDatime procDate;
574  procsec = procDate.Convert(kTRUE);
575 
576  // If file exists, get reco info
577  if(hasReco){
578 
579  // Open Ana file
580  TFile recoFile(anaFileName,"read");
581 
582  // Get number of slices / trigger
583  TH1D* fNumSlices = (TH1D*)recoFile.Get("nearlineana/fNumSlices");
584  // Check if histogram exists
585  if(fNumSlices){
586  // Make sure there are active channels
587  if(nactivechannels!=0){
588  // Fill slicing rate normalised to 10k active channels
589  numslc = 1e4 * fNumSlices->GetMean() / nactivechannels;
590  // Fill pass slicing rate cut
591  pass_slc = ((MinNumSlc < numslc) && (numslc < MaxNumSlc));
592  }
593  }
594  // If histogram failed, file is corrupted
595  else corrupted = 1;
596  //Delete histogram
597  delete fNumSlices;
598 
599  // If this is a new release, also get fraction of tracks that are 2D
600  if(isNewRelana){
601  // Get 2D track fraction histogram
602  TH1D* fTrackFractionAll2D = (TH1D*)recoFile.Get("nearlineana/fTrackFractionAll2D");
603  // Check histogram exists
604  if(fTrackFractionAll2D){
605  // Overall fraction is just the mean
606  trkfrac2D = fTrackFractionAll2D->GetMean();
607  // Fill pass tracking rate cut
608  pass_trk = trkfrac2D < MaxTrkFrac2D;
609  }
610  // If histogram failed, file is corrupted
611  else corrupted = 1;
612  //Delete histogram
613  delete fTrackFractionAll2D;
614  }
615 
616  // Close Ana file
617  recoFile.Clear();
618  recoFile.Close();
619 
620  }// End of reco if
621 
622 
623  // Apply all cuts
624 
625  // Pass reco cuts if both metrics are good
626  if(apply_reco) pass_reco = ( pass_trk && pass_slc );
627  else pass_reco = 1;
628 
629  if (apply_runlen) pass_runlen = ( setlivetime > run_len );
630  else pass_runlen = 1;
631 
632  if (apply_empty) pass_empty = ( emptypercentage < empty_cut );
633  else pass_empty = 1;
634 
635  if (apply_hits) pass_hits = ( midmiprate < hits_up && midmiprate > hits_down );
636  else pass_hits = 1;
637 
638  // Partition 2 requires only 1 good diblock
639  if (apply_db) pass_db = ( (ngoodcdb >= MinGoodDB) || (par == 2 && ngoodcdb > 0) );
640  else pass_db = 1;
641 
642  if (apply_time) pass_time = ( firstsec < lastsec && startDate.GetYear() >= 2013 && endDate.GetYear() >= 2013);
643  else pass_time = 1;
644 
645  // Fill variable for passing all cuts
646  if ( pass_runlen && pass_hits && pass_db && pass_time && pass_empty && pass_reco ){
647 
648  pass_all = 1;
649 
650  }
651  else pass_all = 0;
652 
653  // If subrun is bad, cout more info
654  if ( !pass_all ){
655 
656  cout << "Run " << run << "." << subrun << " has failed: ";
657  if(!pass_runlen) cout << "Run Live Time = " << setlivetime << " seconds; ";
658  if(!pass_empty) cout << "Empty Percentage = " << emptypercentage << "%; ";
659  if(!pass_hits) cout << "Median MIP/Hit Rate = " << midmiprate << "/" << midhitrate << "; ";
660  if(!pass_db) cout << "Good MIP/DiBlocks = " << ngoodmip << "/" << ngooddb << "; ";
661  if(!pass_slc) cout << "Num Slices / 1e4 Pixels = " << numslc << "; ";
662  if(!pass_trk) cout << "Fraction of 2D tracks = " << trkfrac2D << "; ";
663  if(!pass_time){
664  cout << "Bad timing => Diff = " << lastsec - firstsec << ";" << endl;
665  startDate.Print();
666  endDate.Print();
667  }
668  cout << endl;
669 
670  }
671 
672  // Fill tree
673  mytree->Fill();
674 
675  } //End of subrun loop
676 
677  } //End of set loop
678 
679  // Write this run to processed list
680  procfile << TString::Format("%08d",run) << endl;
681 
682  } //End of while loop over run list
683 
684  // Close run list file
685  fileListItr.close();
686 
687  // Close processed runs file
688  procfile.close();
689 
690  // Write tree to root file
691  outfile->cd();
692  mytree->Write("mytree");
693  outfile->Close();
694 
695 
696 } //End of Good Runs
697 
698 ///////////////////////////////////////////////////////////////
699 
700 int GetFileName(TString& filename, int run, int subrun, bool isOnMon){
701 
702  TString fDir;
703  vector<TString> fRel;
704  vector<TString> fTrigDir;
705  vector<TString> fForm;
706  vector<TString> fName;
707  vector<int> fRelType;
708 
709  TString fSubDir = TString::Format("/%06d/%08d/",int(run)/100,run);
710 
711  if(isOnMon){
712 
713  fDir = "/nova/data/nearline-OnMon/FarDet/";
714 
715  fRel.push_back("S15-02-05");
716  fRel.push_back("FA14-09-23");
717  fRel.push_back("S14-08-01");
718  fRel.push_back("S14-02-05");
719  fRel.push_back("S13-09-17");
720 
721  fTrigDir.push_back("");
722  fTrigDir.push_back("t02-t00/");
723 
724  fForm.push_back( TString::Format("FDNL-onmon-endsubrun-%08d-%03d.root",run,subrun) );
725  fForm.push_back( TString::Format("fardet_nearline_r%08d_s%03d.onmon.root",run,subrun) );
726 
727  fName.push_back(fDir + fRel[0] + fSubDir + fTrigDir[1] + fForm[1]);
728  fName.push_back(fDir + fRel[1] + fSubDir + fTrigDir[1] + fForm[1]);
729  fName.push_back(fDir + fRel[2] + fSubDir + fTrigDir[0] + fForm[1]);
730  fName.push_back(fDir + fRel[3] + fSubDir + fTrigDir[0] + fForm[1]);
731  fName.push_back(fDir + fRel[4] + fSubDir + fTrigDir[0] + fForm[0]);
732 
733  fRelType.push_back(2);
734  fRelType.push_back(2);
735  fRelType.push_back(2);
736  fRelType.push_back(2);
737  fRelType.push_back(1);
738 
739  }
740  else{
741 
742  fDir = "/nova/data/nearline-Ana/FarDet/";
743 
744  fRel.push_back("S15-02-05");
745  fRel.push_back("FA14-09-23");
746  fRel.push_back("S14-08-01");
747  fRel.push_back("S14-02-24");
748  fRel.push_back("S14-01-20");
749 
750  fTrigDir.push_back("");
751  fTrigDir.push_back("t02/");
752 
753  fForm.push_back( TString::Format("fardet_r%08d_s%02d_reco_hist_t02.root",run,subrun) );
754 
755  fName.push_back(fDir + fRel[0] + fSubDir + fTrigDir[1] + fForm[0]);
756  fName.push_back(fDir + fRel[1] + fSubDir + fTrigDir[1] + fForm[0]);
757  fName.push_back(fDir + fRel[2] + fSubDir + fTrigDir[0] + fForm[0]);
758  fName.push_back(fDir + fRel[3] + fSubDir + fTrigDir[0] + fForm[0]);
759  fName.push_back(fDir + fRel[4] + fSubDir + fTrigDir[0] + fForm[0]);
760 
761  fRelType.push_back(2);
762  fRelType.push_back(2);
763  fRelType.push_back(2);
764  fRelType.push_back(2);
765  fRelType.push_back(1);
766 
767  }
768 
769  int ntries = 0;
770  int filetype = 0;
771 
772  while(!filetype && ntries < int(fRel.size())){
773 
774  filename = fName[ntries];
775 
776  if(!gSystem->AccessPathName(filename)){
777 
778  filetype = fRelType[ntries];
779 
780  }
781 
782  ntries++;
783 
784  }
785 
786  return filetype;
787 
788 }
789 
790 ///////////////////////////////////////////////////////////////
791 
793 
794  RunInfo myrun;
795 
796  myrun.run = run;
797 
798  bool incomplete = false;
799 
800  for(int subrun=0; subrun<64; subrun++){
801 
802  // Find OnMon file for subrun. Use older release if needed
803 
804  TString omFileName;
805 
806  int omfiletype = GetFileName(omFileName, run, subrun, true);
807 
808  bool isNewRelom = false;
809 
810  if(omfiletype == 0) continue;
811  else if(omfiletype > 1) isNewRelom = true;
812 
813  ////////////////////////////////////////////////////////////////
814  // HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK//
815  // HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK//
816  // HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK//
817  /* HACK */ if((run == 17647 && subrun > 25) // HACK//
818  /* HACK */ || (run == 18531 && subrun == 13)) // HACK//
819  /* HACK */ continue; // HACK//
820  // HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK//
821  // HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK//
822  // HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK//
823  ////////////////////////////////////////////////////////////////
824 
825  TFile omFile(omFileName,"read");
826 
827  TDatime tnow;
828  TDatime tfile = omFile.GetCreationDate();
829 
830  if(tnow.Convert() < tfile.Convert() + 4500){
831  incomplete = true;
832  }
833 
834  // Get live time
835  TH1F *RecordedTime = (TH1F*)omFile.Get("RecordedTime");
836 
837  if(RecordedTime) myrun.livetime.push_back(RecordedTime->Integral()*500e-9);
838  else{
839  incomplete = true;
840  delete RecordedTime;
841  continue;
842  }
843 
844  delete RecordedTime;
845 
846  omFile.Clear();
847  omFile.Close();
848 
849  TString anaFileName;
850 
851  int anafiletype = GetFileName(anaFileName, run, subrun, false);
852 
853  bool isNewRelana = false;
854 
855  if(anafiletype == 0){
856  anaFileName = "";
857  if(tnow.Convert() < tfile.Convert() + 3600*24){
858  incomplete = true;
859  }
860  }
861  else if(anafiletype > 1) isNewRelana = true;
862 
863  myrun.subrun.push_back(subrun);
864  myrun.fnames.push_back(omFileName);
865  myrun.rnames.push_back(anaFileName);
866  myrun.isNewRelom.push_back(isNewRelom);
867  myrun.isNewRelana.push_back(isNewRelana);
868 
869  }// End of subrun loop
870 
871  if(incomplete){
872  cout << "Run " << run << " is incomplete. Marking as preliminary..." << endl;
873  }
874 
875  myrun.isIncomplete = incomplete;
876 
877  return myrun;
878 
879 }
880 
881 ///////////////////////////////////////////////////////////////
882 
883 vector< RunInfo > GetSubrunSets(RunInfo wholerun, double run_len){
884 
885  int nsubruns = wholerun.subrun.size();
886  double remlt = 0;
887  for(int subruni=0; subruni<nsubruns; subruni++){
888  remlt += wholerun.livetime[subruni];
889  }
890 
891  double setlt = 0;
892  int nsets = 0;
893 
894  vector< RunInfo > setvec;
895 
896  RunInfo emptystruct;
897  RunInfo thisset;
898 
899  for(int subruni=0; subruni<nsubruns; subruni++){
900 
901  thisset.subrun.push_back(wholerun.subrun[subruni]);
902  thisset.livetime.push_back(wholerun.livetime[subruni]);
903  thisset.fnames.push_back(wholerun.fnames[subruni]);
904  thisset.rnames.push_back(wholerun.rnames[subruni]);
905  thisset.isNewRelom.push_back(wholerun.isNewRelom[subruni]);
906  thisset.isNewRelana.push_back(wholerun.isNewRelana[subruni]);
907 
908  setlt += thisset.livetime.back();
909  remlt -= thisset.livetime.back();
910 
911  if((setlt > run_len && remlt > run_len) || subruni == nsubruns-1){
912  if(subruni<nsubruns-1) thisset.isLast = false;
913  else thisset.isLast = true;
914  thisset.run = wholerun.run;
915  thisset.isIncomplete = wholerun.isIncomplete;
916  setvec.push_back(thisset);
917  thisset = emptystruct;
918  setlt = 0;
919  nsets++;
920  }
921 
922  }
923 
924  cout << TString::Format("Run %d; %d subruns; %d sets",wholerun.run,nsubruns,nsets) << endl;
925 
926  return setvec;
927 
928 }
929 
930 ///////////////////////////////////////////////////////////////
931 
932 float GetMedian(vector<float> scores)
933 {
934  float median;
935  size_t size = scores.size();
936 
937  if(size==0) return 0;
938 
939  sort(scores.begin(), scores.end());
940 
941  if (size % 2 == 0)
942  {
943  median = (scores[size / 2 - 1] + scores[size / 2]) / 2;
944  }
945  else
946  {
947  median = scores[size / 2];
948  }
949 
950  return median;
951 }
952 
953 double GetMedian(TH1D* hist){
954 
955  int nbins = hist->GetNbinsX();
956 
957  double *x = new double[nbins];
958  hist->GetXaxis()->GetCenter(x);
959 
960  double *y = hist->GetArray();
961 
962  double out = TMath::Median(nbins,x,y);
963 
964  delete x;
965 
966  return out;
967 
968 }
969 
970 ///////////////////////////////////////////////////////////////
971 
972 bool GetDBPartition(int run, int &par){
973 
974  par = 4;
975 
976  std::ifstream t;
977  t.open(gSystem->Getenv("NOVADBPWDFILE")); // open input file
979  getline(t,buffer);
980  //cout<<buffer<<endl;
981  char * dbpass;
982  dbpass = new char [8];
983  strcpy(dbpass, buffer.c_str());
984  //cout<<dbpass<<endl;
985 
986  std::stringstream dbstream;
987  std::string dbstring;
988  dbstream << "select run, partition from fardet.runs where run = "
989  << run
990  << " order by run asc";
991  dbstring = dbstream.str();
992  // cout<<dbstring<<endl;
993  char * dbcstr;
994  dbcstr = new char [dbstring.size()+1];
995  strcpy(dbcstr, dbstring.c_str());
996  TSQLServer *dbserv = TSQLServer::Connect("pgsql://ifdbprod.fnal.gov:5433/nova_prod", "nova_reader", dbpass);
997 
998  int attempts = 1;
999 
1000  while ( (dbserv==0 || !dbserv->IsConnected()) && attempts<100 ) {
1001 
1002  cout << "Try to connect again..." << endl;
1003 
1004  dbserv = TSQLServer::Connect("pgsql://ifdbprod.fnal.gov:5433/nova_prod", "nova_reader", dbpass);
1005 
1006  attempts++;
1007 
1008  }
1009 
1010  bool connected = false;
1011 
1012  if ( dbserv!=0 && dbserv->IsConnected() ) {
1013 
1014  connected = true;
1015 
1016  TSQLStatement *stmt = dbserv->Statement(dbcstr,1000);
1017 
1018  if(stmt->Process()) {
1019  // store result of statement in buffer
1020  stmt->StoreResult();
1021  while (stmt->NextResultRow()) {
1022  // Extract rows one after another
1023  par = stmt->GetUInt(1);
1024  }
1025  }
1026  // std::cout << "\n" << run << " " << subrun << " " << par;
1027 
1028  delete stmt;
1029  }
1030 
1031  delete dbserv;
1032  delete dbcstr;
1033  delete dbpass;
1034 
1035  if(attempts>1 && connected) cout << "Connected after " << attempts << " attempts." << endl;
1036  if(!connected) cout << "Failed after " << attempts << " attempts." << endl;
1037 
1038  return connected;
1039 
1040 }
vector< int > subrun
Int_t pass_trk
Definition: SimpleIterate.C:58
Int_t lastsec
Definition: SimpleIterate.C:27
Double_t numslc
Definition: SimpleIterate.C:61
int GetFileName(TString &filename, int run, int subrun, bool isOnMon=true)
vector< TString > fnames
void BuildMetricsTree_OnMon(TString filelist)
Int_t ngooddcm
Definition: SimpleIterate.C:30
Int_t pass_all
Definition: SimpleIterate.C:42
Double_t trkfrac2D
Definition: SimpleIterate.C:62
double GetMedian(TH1D *hist)
Int_t ngooddb
Definition: SimpleIterate.C:31
constexpr T pow(T x)
Definition: pow.h:75
Int_t par
Definition: SimpleIterate.C:24
Double_t nactivechannels
Definition: SimpleIterate.C:52
string filename
Definition: shutoffs.py:106
Int_t dbencoded
Definition: SimpleIterate.C:54
float Gamma() const
bool GetDBPartition(int run, int &par)
Double_t miprate
Definition: SimpleIterate.C:48
Double_t rectimesec
Definition: SimpleIterate.C:44
Int_t nevents
Definition: SimpleIterate.C:25
RunInfo GetRunInfo(int run)
::xsd::cxx::tree::buffer< char > buffer
Definition: Database.h:179
vector< RunInfo > GetSubrunSets(RunInfo wholerun, double run_len)
const int nbins
Definition: cellShifts.C:15
Int_t pass_empty
Definition: SimpleIterate.C:38
Int_t nactivefeb
Definition: SimpleIterate.C:34
Int_t dbaencoded
Definition: SimpleIterate.C:55
Double_t midhitrate
Definition: SimpleIterate.C:47
Int_t nactivedcm
Definition: SimpleIterate.C:35
Int_t ngoodmip
Definition: SimpleIterate.C:32
Double_t emptypercentage
Definition: SimpleIterate.C:53
Int_t preliminary
Definition: SimpleIterate.C:63
Int_t ngoodcdb
Definition: SimpleIterate.C:33
TTree * mytree
Definition: SimpleIterate.C:18
double median(TH1D *h)
Definition: absCal.cxx:524
Definition: run.py:1
Double_t midmiprate
Definition: SimpleIterate.C:49
Int_t nactivedb
Definition: SimpleIterate.C:36
OStream cout
Definition: OStream.cxx:6
Int_t procsec
Definition: SimpleIterate.C:60
Int_t ngoodpix
Definition: SimpleIterate.C:28
vector< bool > isNewRelana
Int_t firstsec
Definition: SimpleIterate.C:26
Double_t mipratio
Definition: SimpleIterate.C:50
ifstream fileListItr("ManBadSubruns.txt")
vector< double > livetime
Int_t pass_db
Definition: SimpleIterate.C:40
vector< bool > isNewRelom
ofstream procfile
Double_t hitrate
Definition: SimpleIterate.C:46
Int_t pass_slc
Definition: SimpleIterate.C:57
vector< TString > rnames
Int_t setsize
Definition: SimpleIterate.C:43
Double_t setlivetime
Definition: SimpleIterate.C:45
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
Float_t e
Definition: plot.C:35
Int_t pass_reco
Definition: SimpleIterate.C:56
Double_t mipasym
Definition: SimpleIterate.C:51
Int_t ngoodfeb
Definition: SimpleIterate.C:29
Int_t pass_runlen
Definition: SimpleIterate.C:37
Int_t corrupted
Definition: SimpleIterate.C:59
FILE * outfile
Definition: dump_event.C:13
Int_t pass_time
Definition: SimpleIterate.C:41
Int_t pass_hits
Definition: SimpleIterate.C:39
enum BeamMode string