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;
29  bool isIncomplete;
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 = 28.0; //Median MIP rate smaller than this parameter
66  Double_t hits_down = 18.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.02; //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] = {17.3, 9.5, 9.3, 7.5, 3.8, 3.8, 1.9}; //Minimum MIP rate for pixels in different DCMs
81  float highpixmip[7] = {36, 85, 52, 44, 36, 32, 28}; //Maximum MIP rate for pixels in different DCMs
82 
83  float lowfebmip[7] = {571, 492, 432, 360, 191, 212, 190}; //Minimum MIP rate for FEBs in different DCMs
84  float highfebmip[7] = {1150, 2100, 1450, 1200, 1050, 950, 850}; //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) < .95 ){
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  // Determine if DCM is good, i.e. enough good FEBs
453  if(nGoodFEBInDCM >= FEBsPerDCM){
454  // Increment good DCMs in total and in this diblock
455  ngooddcm++;
456  nGoodDCMInDB++;
457  }
458 
459  // Delete pixel rate map
460  delete PixelRate;
461 
462  }// End of DCM loop
463 
464  // If active DCMs has increased, this diblock is active.
465  // Increment active diblocks
466  if(tempdcm < nactivedcm){
467  nactivedb++;
468  // Also, fill active diblock mask
469  dbaencoded += pow(2,db-1);
470  }
471 
472  // Fill vectors for whole detector with hits from this diblock
473  allMipHitsTop.insert(allMipHitsTop.end(), pixMipHitsTop.begin(), pixMipHitsTop.end());
474  allMipHitsSide.insert(allMipHitsSide.end(), pixMipHitsSide.begin(), pixMipHitsSide.end());
475 
476  // Find median of MIP hits for top and side of this diblock
477  double medianTop = GetMedian(pixMipHitsTop);
478  double medianSide = GetMedian(pixMipHitsSide);
479 
480  // Calculate top-side asymmetry if there are MIP hits
481  if(medianTop + medianSide != 0) medianTop = (medianTop - medianSide) / (medianTop + medianSide);
482  else medianTop = 0;
483  // Determine if diblock is good
484  if(nGoodDCMInDB >= DCMsPerDB){
485  // Increment good diblocks (DCMs statuses only)
486  ngooddb++;
487  // Check for bad top-side asymmetry
488  if(medianTop < MaxTopSideAsym &&
489  medianTop > MinTopSideAsym ){
490  // Increment good diblocks (including asymmetry)
491  ngoodmip++;
492  // Counting consecutive good diblocks
493  tempdb++;
494  // Fill good diblock mask
495  dbencoded += pow(2,db-1);
496  }
497  // Else, reset consecutive diblocks
498  else tempdb = 0;
499  }
500  // Else, reset consecutive diblocks
501  else tempdb = 0;
502 
503  // Fill largest number of consecutive diblocks
504  if(tempdb > ngoodcdb) ngoodcdb = tempdb;
505 
506  }// End of diblock loop
507 
508 
509  // Compute median of all MIP hits for top and side
510  double medianAllTop = GetMedian(allMipHitsTop);
511  double medianAllSide = GetMedian(allMipHitsSide);
512 
513  // Calculate overall top-side asymmetry
514  if(medianAllTop + medianAllSide != 0) mipasym = (medianAllTop - medianAllSide) / (medianAllTop + medianAllSide);
515  else mipasym = 0;
516 
517  // Calculate top-side ratio
518  if(medianAllSide!=0) mipratio = medianAllTop / medianAllSide;
519  else mipratio = -1;
520 
521  // Normalise hit rate per active channel (mean)
522  if(nactivechannels!=0) hitrate /= nactivechannels;
523  else hitrate = 0;
524 
525  // Get overall median of hit rates
526  midhitrate = GetMedian(allHitRates);
527 
528  // Add top and side MIP hits and compute median
529  allMipHitsTop.insert(allMipHitsTop.end(), allMipHitsSide.begin(), allMipHitsSide.end());
530  midmiprate = GetMedian(allMipHitsTop);
531  if(setlivetime != 0) midmiprate /= setlivetime;
532  else midmiprate = 0;
533 
534  // Reset mean MIP rate
535  miprate = 0;
536 
537  // Add all MIP hits
538  for(int i=0;i<int(allMipHitsTop.size());i++){
539  miprate += allMipHitsTop[i];
540  }
541 
542  // Compute mean MIP rate
543  double chtime = allMipHitsTop.size() * setlivetime;
544  if(chtime!=0) miprate /= chtime;
545  else miprate = 0;
546 
547  // Close OnMon file
548  omFile.Clear();
549  omFile.Close();
550 
551 
552  // Start looking at reconstructed variables
553 
554  // Initialize variables. If no nearline file exists, subrun passes reco cut.
555  pass_reco = 1;
556  pass_trk = 1;
557  pass_slc = 1;
558  corrupted = 0;
559  numslc = -1;
560  trkfrac2D = -1;
561 
562  // Get file name and release type
563  TString anaFileName = setvec[seti].rnames[subruni];
564  bool isNewRelana = setvec[seti].isNewRelana[subruni];
565 
566  // Check if file exists
567  bool hasReco = true;
568  if(anaFileName == "") hasReco = false;
569 
570  // Fill timestamp in which this tree was filled
571  TDatime procDate;
572  procsec = procDate.Convert(kTRUE);
573 
574  // If file exists, get reco info
575  if(hasReco){
576 
577  // Open Ana file
578  TFile recoFile(anaFileName,"read");
579 
580  // Get number of slices / trigger
581  TH1D* fNumSlices = (TH1D*)recoFile.Get("nearlineana/fNumSlices");
582  // Check if histogram exists
583  if(fNumSlices){
584  // Make sure there are active channels
585  if(nactivechannels!=0){
586  // Fill slicing rate normalised to 10k active channels
587  numslc = 1e4 * fNumSlices->GetMean() / nactivechannels;
588  // Fill pass slicing rate cut
589  pass_slc = ((MinNumSlc < numslc) && (numslc < MaxNumSlc));
590  }
591  }
592  // If histogram failed, file is corrupted
593  else corrupted = 1;
594  //Delete histogram
595  delete fNumSlices;
596 
597  // If this is a new release, also get fraction of tracks that are 2D
598  if(isNewRelana){
599  // Get 2D track fraction histogram
600  TH1D* fTrackFractionAll2D = (TH1D*)recoFile.Get("nearlineana/fTrackFractionAll2D");
601  // Check histogram exists
602  if(fTrackFractionAll2D){
603  // Overall fraction is just the mean
604  trkfrac2D = fTrackFractionAll2D->GetMean();
605  // Fill pass tracking rate cut
606  pass_trk = trkfrac2D < MaxTrkFrac2D;
607  }
608  // If histogram failed, file is corrupted
609  else corrupted = 1;
610  //Delete histogram
611  delete fTrackFractionAll2D;
612  }
613 
614  // Close Ana file
615  recoFile.Clear();
616  recoFile.Close();
617 
618  }// End of reco if
619 
620 
621  // Apply all cuts
622 
623  // Pass reco cuts if both metrics are good
624  // pass reco
625  if(apply_reco) pass_reco = ( pass_trk && pass_slc );
626  else pass_reco = 1;
627  //pass live time
628  if (apply_runlen) pass_runlen = ( setlivetime > run_len );
629  else pass_runlen = 1;
630  //pass diblock
631  if (apply_empty) pass_empty = ( emptypercentage < empty_cut );
632  else pass_empty = 1;
633  //pass hit rate
634  if (apply_hits) pass_hits = ( midmiprate < hits_up && midmiprate > hits_down );
635  else pass_hits = 1;
636 
637  // Partition 2 requires only 1 good diblock
638  //pass partial detector
639  if (apply_db) pass_db = ( (ngoodcdb >= MinGoodDB) || (par == 2 && ngoodcdb > 0) );
640  else pass_db = 1;
641  //pass other
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 
716 // fRel.push_back("development");
717  fRel.push_back("S19-10-22");
718  fRel.push_back("S18-09-02");
719  fRel.push_back("S18-07-03");
720  fRel.push_back("S18-01-05");
721  fRel.push_back("S17-10-10");
722  fRel.push_back("S17-07-18");
723  fRel.push_back("S16-10-26");
724  fRel.push_back("S16-10-07");
725  fRel.push_back("S16-08-23");
726  fRel.push_back("S15-11-06");
727  fRel.push_back("S15-02-05");
728  fRel.push_back("FA14-09-23");
729  fRel.push_back("S14-08-01");
730  fRel.push_back("S14-02-05");
731  fRel.push_back("S13-09-17");
732 
733  fTrigDir.push_back("");
734  fTrigDir.push_back("t02-t00/");
735 
736  fForm.push_back( TString::Format("FDNL-onmon-endsubrun-%08d-%03d.root",run,subrun) );
737  fForm.push_back( TString::Format("fardet_nearline_r%08d_s%03d.onmon.root",run,subrun) );
738 
739  fName.push_back(fDir + fRel[0] + fSubDir + fTrigDir[1] + fForm[1]);
740  fName.push_back(fDir + fRel[1] + fSubDir + fTrigDir[1] + fForm[1]);
741  fName.push_back(fDir + fRel[2] + fSubDir + fTrigDir[1] + fForm[1]);
742  fName.push_back(fDir + fRel[3] + fSubDir + fTrigDir[1] + fForm[1]);
743  fName.push_back(fDir + fRel[4] + fSubDir + fTrigDir[1] + fForm[1]);
744  fName.push_back(fDir + fRel[5] + fSubDir + fTrigDir[1] + fForm[1]);
745  fName.push_back(fDir + fRel[6] + fSubDir + fTrigDir[1] + fForm[1]);
746  fName.push_back(fDir + fRel[7] + fSubDir + fTrigDir[1] + fForm[1]);
747  fName.push_back(fDir + fRel[8] + fSubDir + fTrigDir[1] + fForm[1]);
748  fName.push_back(fDir + fRel[9] + fSubDir + fTrigDir[1] + fForm[1]);
749  fName.push_back(fDir + fRel[10] + fSubDir + fTrigDir[1] + fForm[1]);
750  fName.push_back(fDir + fRel[11] + fSubDir + fTrigDir[1] + fForm[1]);
751  fName.push_back(fDir + fRel[12] + fSubDir + fTrigDir[0] + fForm[1]);
752  fName.push_back(fDir + fRel[13] + fSubDir + fTrigDir[0] + fForm[0]);
753  fName.push_back(fDir + fRel[14] + fSubDir + fTrigDir[0] + fForm[0]);
754 
755 
756 
757  fRelType.push_back(2);
758  fRelType.push_back(2);
759  fRelType.push_back(2);
760  fRelType.push_back(2);
761  fRelType.push_back(2);
762  fRelType.push_back(2);
763  fRelType.push_back(2);
764  fRelType.push_back(2);
765  fRelType.push_back(2);
766  fRelType.push_back(2);
767  fRelType.push_back(2);
768  fRelType.push_back(2);
769  fRelType.push_back(2);
770  fRelType.push_back(2);
771  fRelType.push_back(1);
772 
773  }
774  else{
775 
776  fDir = "/nova/data/nearline-Ana/FarDet/";
777 
778  fRel.push_back("S19-10-22");
779  fRel.push_back("S18-09-02");
780  fRel.push_back("S18-07-03");
781  fRel.push_back("S18-01-05");
782  fRel.push_back("S17-10-10");
783  fRel.push_back("S17-07-18");
784  fRel.push_back("S16-10-26");
785  fRel.push_back("S16-10-07");
786  fRel.push_back("S16-08-23");
787  fRel.push_back("S15-11-06");
788  fRel.push_back("S15-02-05");
789  fRel.push_back("FA14-09-23");
790  fRel.push_back("S14-08-01");
791  fRel.push_back("S14-02-24");
792  fRel.push_back("S14-01-20");
793 
794  fTrigDir.push_back("");
795  fTrigDir.push_back("t02/");
796 
797  fForm.push_back( TString::Format("fardet_r%08d_s%02d_reco_hist_t02.root",run,subrun) );
798 
799  fName.push_back(fDir + fRel[0] + fSubDir + fTrigDir[1] + fForm[0]);
800  fName.push_back(fDir + fRel[1] + fSubDir + fTrigDir[1] + fForm[0]);
801  fName.push_back(fDir + fRel[2] + fSubDir + fTrigDir[1] + fForm[0]);
802  fName.push_back(fDir + fRel[3] + fSubDir + fTrigDir[1] + fForm[0]);
803  fName.push_back(fDir + fRel[4] + fSubDir + fTrigDir[1] + fForm[0]);
804  fName.push_back(fDir + fRel[5] + fSubDir + fTrigDir[1] + fForm[0]);
805  fName.push_back(fDir + fRel[6] + fSubDir + fTrigDir[1] + fForm[0]);
806  fName.push_back(fDir + fRel[7] + fSubDir + fTrigDir[1] + fForm[0]);
807  fName.push_back(fDir + fRel[8] + fSubDir + fTrigDir[1] + fForm[0]);
808  fName.push_back(fDir + fRel[9] + fSubDir + fTrigDir[1] + fForm[0]);
809  fName.push_back(fDir + fRel[10] + fSubDir + fTrigDir[1] + fForm[0]);
810  fName.push_back(fDir + fRel[11] + fSubDir + fTrigDir[0] + fForm[0]);
811  fName.push_back(fDir + fRel[12] + fSubDir + fTrigDir[0] + fForm[0]);
812  fName.push_back(fDir + fRel[13] + fSubDir + fTrigDir[0] + fForm[0]);
813  fName.push_back(fDir + fRel[14] + fSubDir + fTrigDir[0] + fForm[0]);
814 
815 
816 
817 
818  fRelType.push_back(2);
819  fRelType.push_back(2);
820  fRelType.push_back(2);
821  fRelType.push_back(2);
822  fRelType.push_back(2);
823  fRelType.push_back(2);
824  fRelType.push_back(2);
825  fRelType.push_back(2);
826  fRelType.push_back(2);
827  fRelType.push_back(2);
828  fRelType.push_back(2);
829  fRelType.push_back(2);
830  fRelType.push_back(2);
831  fRelType.push_back(2);
832  fRelType.push_back(1);
833 
834 
835  std::cout << fName[0] << std::endl;
836 
837  }
838 
839  int ntries = 0;
840  int filetype = 0;
841 
842  while(!filetype && ntries < int(fRel.size())){
843 
844  filename = fName[ntries];
845 
846  if(!gSystem->AccessPathName(filename)){
847 
848  filetype = fRelType[ntries];
849 
850  }
851 
852  ntries++;
853 
854  }
855 
856  return filetype;
857 
858 }
859 
860 ///////////////////////////////////////////////////////////////
861 
863 
864  RunInfo myrun;
865 
866  myrun.run = run;
867 
868  bool incomplete = false;
869 
870  for(int subrun=0; subrun<64; subrun++){
871 
872  // Find OnMon file for subrun. Use older release if needed
873 
874  TString omFileName;
875 
876  int omfiletype = GetFileName(omFileName, run, subrun, true);
877 
878  bool isNewRelom = false;
879 
880  if(omfiletype == 0) continue;
881  else if(omfiletype > 1) isNewRelom = true;
882 
883  ////////////////////////////////////////////////////////////////
884  // HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK//
885  // HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK//
886  // HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK//
887  /* HACK */ if((run == 17647 && subrun > 25) // HACK//
888  /* HACK */ || (run == 18531 && subrun == 13)) // HACK//
889  /* HACK */ continue; // HACK//
890  // HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK//
891  // HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK//
892  // HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK//
893  ////////////////////////////////////////////////////////////////
894 
895  TFile omFile(omFileName,"read");
896 
897  TDatime tnow;
898  TDatime tfile = omFile.GetCreationDate();
899  TDatime ttest = omFile.GetModificationDate();
900 
901  if(tnow.Convert() < tfile.Convert() + 4500){
902  incomplete = true;
903  std::cout << "File too new" << std::endl;
904  }
905 
906  // Get live time
907  TH1F *RecordedTime = (TH1F*)omFile.Get("RecordedTime");
908 
909  if(RecordedTime) myrun.livetime.push_back(RecordedTime->Integral()*500e-9);
910  else{
911  incomplete = true;
912  std::cout << "No RecordedTime" << std::endl;
913  delete RecordedTime;
914  continue;
915  }
916 
917  delete RecordedTime;
918 
919  omFile.Clear();
920  omFile.Close();
921 
922  TString anaFileName;
923 
924  int anafiletype = GetFileName(anaFileName, run, subrun, false);
925 
926  bool isNewRelana = false;
927 
928  if(anafiletype == 0){
929  anaFileName = "";
930  if(tnow.Convert() < tfile.Convert() + 3600*24){
931  incomplete = true;
932  std::cout << "No Ana File" << std::endl;
933  }
934  }
935  else if(anafiletype > 1) isNewRelana = true;
936 
937  myrun.subrun.push_back(subrun);
938  myrun.fnames.push_back(omFileName);
939  myrun.rnames.push_back(anaFileName);
940  myrun.isNewRelom.push_back(isNewRelom);
941  myrun.isNewRelana.push_back(isNewRelana);
942 
943  }// End of subrun loop
944 
945  if(incomplete){
946  cout << "Run " << run << " is incomplete. Marking as preliminary..." << endl;
947  }
948 
949  myrun.isIncomplete = incomplete;
950 
951  return myrun;
952 
953 }
954 
955 ///////////////////////////////////////////////////////////////
956 
957 vector< RunInfo > GetSubrunSets(RunInfo wholerun, double run_len){
958 
959  int nsubruns = wholerun.subrun.size();
960  double remlt = 0;
961  for(int subruni=0; subruni<nsubruns; subruni++){
962  remlt += wholerun.livetime[subruni];
963  }
964 
965  double setlt = 0;
966  int nsets = 0;
967 
968  vector< RunInfo > setvec;
969 
970  RunInfo emptystruct;
971  RunInfo thisset;
972 
973  for(int subruni=0; subruni<nsubruns; subruni++){
974 
975  thisset.subrun.push_back(wholerun.subrun[subruni]);
976  thisset.livetime.push_back(wholerun.livetime[subruni]);
977  thisset.fnames.push_back(wholerun.fnames[subruni]);
978  thisset.rnames.push_back(wholerun.rnames[subruni]);
979  thisset.isNewRelom.push_back(wholerun.isNewRelom[subruni]);
980  thisset.isNewRelana.push_back(wholerun.isNewRelana[subruni]);
981 
982  setlt += thisset.livetime.back();
983  remlt -= thisset.livetime.back();
984 
985  if((setlt > run_len && remlt > run_len) || subruni == nsubruns-1){
986  if(subruni<nsubruns-1) thisset.isLast = false;
987  else thisset.isLast = true;
988  thisset.run = wholerun.run;
989  thisset.isIncomplete = wholerun.isIncomplete;
990  setvec.push_back(thisset);
991  thisset = emptystruct;
992  setlt = 0;
993  nsets++;
994  }
995 
996  }
997 
998  cout << TString::Format("Run %d; %d subruns; %d sets",wholerun.run,nsubruns,nsets) << endl;
999 
1000  return setvec;
1001 
1002 }
1003 
1004 ///////////////////////////////////////////////////////////////
1005 
1006 float GetMedian(vector<float> scores)
1007 {
1008  float median;
1009  size_t size = scores.size();
1010 
1011  if(size==0) return 0;
1012 
1013  sort(scores.begin(), scores.end());
1014 
1015  if (size % 2 == 0)
1016  {
1017  median = (scores[size / 2 - 1] + scores[size / 2]) / 2;
1018  }
1019  else
1020  {
1021  median = scores[size / 2];
1022  }
1023 
1024  return median;
1025 }
1026 
1027 double GetMedian(TH1D* hist){
1028 
1029  int nbins = hist->GetNbinsX();
1030 
1031  double *x = new double[nbins];
1032  hist->GetXaxis()->GetCenter(x);
1033 
1034  double *y = hist->GetArray();
1035 
1036  double out = TMath::Median(nbins,x,y);
1037 
1038  delete x;
1039 
1040  return out;
1041 
1042 }
1043 
1044 ///////////////////////////////////////////////////////////////
1045 
1046 bool GetDBPartition(int run, int &par){
1047 
1048  par = 4;
1049 
1050  std::ifstream t;
1051  t.open(gSystem->Getenv("NOVADBPWDFILE")); // open input file
1053  getline(t,buffer);
1054  //cout<<buffer<<endl;
1055  char * dbpass;
1056  dbpass = new char [8];
1057  strcpy(dbpass, buffer.c_str());
1058  //cout<<dbpass<<endl;
1059 
1060  std::stringstream dbstream;
1061  std::string dbstring;
1062  dbstream << "select run, partition from fardet.runs where run = "
1063  << run
1064  << " order by run asc";
1065  dbstring = dbstream.str();
1066  // cout<<dbstring<<endl;
1067  char * dbcstr;
1068  dbcstr = new char [dbstring.size()+1];
1069  strcpy(dbcstr, dbstring.c_str());
1070  TSQLServer *dbserv = TSQLServer::Connect("pgsql://ifdbprod.fnal.gov:5433/nova_prod", "nova_reader", dbpass);
1071 
1072  int attempts = 1;
1073 
1074  while ( (dbserv==0 || !dbserv->IsConnected()) && attempts<100 ) {
1075 
1076  cout << "Try to connect again..." << endl;
1077 
1078  dbserv = TSQLServer::Connect("pgsql://ifdbprod.fnal.gov:5433/nova_prod", "nova_reader", dbpass);
1079 
1080  attempts++;
1081 
1082  }
1083 
1084  bool connected = false;
1085 
1086  if ( dbserv!=0 && dbserv->IsConnected() ) {
1087 
1088  connected = true;
1089 
1090  TSQLStatement *stmt = dbserv->Statement(dbcstr,1000);
1091 
1092  if(stmt->Process()) {
1093  // store result of statement in buffer
1094  stmt->StoreResult();
1095  while (stmt->NextResultRow()) {
1096  // Extract rows one after another
1097  par = stmt->GetUInt(1);
1098  }
1099  }
1100  // std::cout << "\n" << run << " " << subrun << " " << par;
1101 
1102  delete stmt;
1103  }
1104 
1105  delete dbserv;
1106  delete dbcstr;
1107  delete dbpass;
1108 
1109  if(attempts>1 && connected) cout << "Connected after " << attempts << " attempts." << endl;
1110  if(!connected) cout << "Failed after " << attempts << " attempts." << endl;
1111 
1112  return connected;
1113 
1114 
1115 }
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
double livetime
Definition: saveFDMCHists.C:21
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