BuildMetricsTree_OnMon.C
Go to the documentation of this file.
1 R__LOAD_LIBRARY(libDAQChannelMap.so)
2 
3 #include <iostream>
4 #include <fstream>
5 #include <iomanip>
6 #include <sstream>
7 #include <string>
8 #include <streambuf>
9 #include <cstring>
10 #include <algorithm>
11 
12 #include <TStyle.h>
13 #include <TSystem.h>
14 #include <TFile.h>
15 #include <TString.h>
16 #include <TTree.h>
17 #include <TH2F.h>
18 #include <TH1F.h>
19 #include <TMath.h>
20 
21 #include "DAQChannelMap/HardwareDisplay.h"
22 #include "DAQChannelMap/DAQChannelMap.h"
23 #include "NovaDAQConventions/DAQConventions.h"
24 
25 using namespace std;
26 
27 float GetMedian(vector<float> scores);
28 
29 bool IsAntiNeutrinoRun(int run);
30 
31 int GetFileName(TString& filename, int run, int subrun, bool isOnMon);
32 
34 
35  gStyle->SetTimeOffset(0);
36 
37  int apply_time = 1; //1 if want to apply negative time cut; 0 if do not (start time before end time)
38  int apply_db = 1; //1 if want to apply diblock cut; 0 if do not
39  int apply_hits = 1; //1 if want to apply hit rate cut; 0 if do not
40  int apply_reco = 1; //1 if want to apply reco cut; 0 if do not
41  int apply_slc = 1; //1 if want to apply reco cut; 0 if do not
42  int apply_trk = 0; //1 if want to apply reco cut; 0 if do not
43 
44  int apply_emptyspill = 1;
45  int apply_timingpeak = 1;
46  int apply_duration = 1;
47 
48  //Cut Parameter Values
49  double NomtimingPeakStart = 218e3; // start time of timing peak (ns)
50  double NomtimingPeakEnd = 228e3; // end time of timing peak (ns)
51  double TimingPeakWindow = 1e3; // allowable error in timing peak edges (ns)
52  double MinTrack3D = 0.95; // min 3d track fraction
53  double numslc_up = 5.5; // maximum average number of slices after goodbeam
54  double numslc_down = 3.5;
55  double emptyspillfract_up = 0.03;
56  double minDuration = 1000; //
57  double hits_up = 20.0; //Median MIP rate smaller than this parameter
58  double hits_down = 12.0; //Median MIP rate larger than this parameter
59 
60 
61  //Values that changes for antineutrino (RHC) mode
62  //PLEASE NOTE: pixhitratepotscaled MUST ALSO BE CHANGED!
63  double rhc_hits_up = 10.0;
64  double rhc_hits_down = 6.0;
65  double rhc_emptyspillfract_up = 0.05;
66  double rhc_numslc_up = 3.5;
67  double rhc_numslc_down = 1.5;
68  //
69 
70  double MaxTopSideAsymActive = 0.1; //Maximum Top/Side MIP hit rate Asymmetry in active DBs
71  double MinTopSideAsymActive = -0.1; //Minimum Top/Side MIP hit rate Asymmetry in active DBs
72  double MaxTopSideAsymMuonC = 0.1; //Maximum Top/Side MIP hit rate Asymmetry in muon catcher
73  double MinTopSideAsymMuonC = -0.3; //Minimum Top/Side MIP hit rate Asymmetry in muon catcher
74 
75  double MinGoodDB = 4; //Minimum # of good consecutive DiBlocks
76  double MinGoodFEBFract = 0.80;
77  double MinGoodPixPerFEB = 32 - 7; //Minimum # of good pixels that an FEB must have
78 
79  // **************END OF TYPICAL USER'S NEED TO ALTER THIS SCRIPT*******************************************
80 
81  // Open file to keep track of processed subruns
82  ofstream procfile;
83  procfile.open("listProcessedTemp.txt");
84 
85  // Create tree to store Good Run metrics
86  int run, subrun, par, nevents;
87  int firstsec, lastsec;
88 
90  int pass_reco, pass_slc, pass_trk3D;
91  int pass_timingpeak, pass_duration;
92  int passemptyspill;
93 
96 
99 
100  double MIPTopSideSym[4];
101 
102  double goodbeam, pot, emptyspillfract;
103  double hitratepotscaled, midmipratepotscaled;
104 
105  double timingPeakStart, timingPeakEnd;
106 
107  int corrupted, procsec;
108  int hasReco, preliminary;
109 
110  double numslc, trkfrac3D;
111  double numslcpotscaled;
112 
113  TFile *outfile = new TFile("MetricsNew.root","recreate");
114 
115  TTree *mytree = new TTree("evts","");
116  mytree->Branch("run", &run, "run/I");
117  mytree->Branch("subrun", &subrun, "subrun/I");
118  mytree->Branch("par", &par, "par/I");
119  mytree->Branch("nevents", &nevents, "nevents/I");
120  mytree->Branch("firstsec", &firstsec, "firstsec/I");
121  mytree->Branch("lastsec", &lastsec, "lastsec/I");
122 
123  mytree->Branch("pass_all", &pass_all, "pass_all/I");
124  mytree->Branch("pass_hits", &pass_hits, "pass_hits/I");
125  mytree->Branch("pass_db", &pass_db, "pass_db/I");
126  mytree->Branch("pass_time", &pass_time, "pass_time/I");
127  mytree->Branch("pass_reco", &pass_reco, "pass_reco/I");
128  mytree->Branch("pass_slc", &pass_slc, "pass_slc/I");
129  mytree->Branch("pass_trk3D", &pass_trk3D, "pass_trk3D/I");
130 
131  mytree->Branch("pass_timingpeak", &pass_timingpeak, "pass_timingpeak/I");
132  mytree->Branch("pass_duration", &pass_duration, "pass_duration/I");
133  mytree->Branch("passemptyspill", &passemptyspill, "passemptyspill/I");
134 
135  mytree->Branch("ngoodpix", &ngoodpix, "ngoodpix/I");
136  mytree->Branch("ngoodfeb", &ngoodfeb, "ngoodfeb/I");
137  mytree->Branch("ngooddcm", &ngooddcm, "ngooddcm/I");
138  mytree->Branch("ngooddb", &ngooddb, "ngooddb/I");
139  mytree->Branch("ngoodmip", &ngoodmip, "ngoodmip/I");
140  mytree->Branch("ngoodcdb", &ngoodcdb, "ngoodcdb/I");
141  mytree->Branch("nactivefeb", &nactivefeb, "nactivefeb/I");
142  mytree->Branch("nactivedcm", &nactivedcm, "nactivedcm/I");
143  mytree->Branch("nactivedb", &nactivedb, "nactivedb/I");
144 
145  mytree->Branch("duration", &duration, "duration/D");
146  mytree->Branch("rectimesec", &rectimesec, "rectimesec/D");
147  mytree->Branch("nactivechannels", &nactivechannels, "nactivechannels/D");
148 
149  mytree->Branch("hitrate", &hitrate, "hitrate/D");
150  mytree->Branch("miprate", &miprate, "miprate/D");
151  mytree->Branch("mipasym", &mipasym, "mipasym/D");
152  mytree->Branch("midhitrate", &midhitrate, "midhitrate/D");
153  mytree->Branch("midmiprate", &midmiprate, "midmiprate/D");
154 
155  mytree->Branch("mipasymdb1", &MIPTopSideSym[0], "mipasymdb1/D");
156  mytree->Branch("mipasymdb2", &MIPTopSideSym[1], "mipasymdb2/D");
157  mytree->Branch("mipasymdb3", &MIPTopSideSym[2], "mipasymdb3/D");
158  mytree->Branch("mipasymdb4", &MIPTopSideSym[3], "mipasymdb4/D");
159 
160  mytree->Branch("pot", &pot, "pot/D");
161  mytree->Branch("goodbeam", &goodbeam, "goodbeam/D");
162 
163  mytree->Branch("emptyspillfract", &emptyspillfract, "emptyspillfract/D");
164  mytree->Branch("hitratepotscaled", &hitratepotscaled, "hitratepotscaled/D");
165  mytree->Branch("midmipratePoTscaled", &midmipratepotscaled, "midmipratePoTscaled/D");
166 
167  mytree->Branch("timingPeakStart", &timingPeakStart, "timingPeakStart/D");
168  mytree->Branch("timingPeakEnd", &timingPeakEnd, "timingPeakEnd/D");
169 
170  mytree->Branch("hasReco", &hasReco, "hasReco/I");
171  mytree->Branch("preliminary", &preliminary, "preliminary/I");
172  mytree->Branch("corrupted", &corrupted, "corrupted/I");
173  mytree->Branch("procsec", &procsec, "procsec/I");
174 
175  mytree->Branch("numslc", &numslc, "numslc/D");
176  mytree->Branch("trkfrac3D", &trkfrac3D, "trkfrac3D/D");
177  mytree->Branch("numslcpotscaled", &numslcpotscaled, "numslcpotscaled/D");
178 
179  // Open list of subruns to process
180  ifstream fileListItr(filelist.Data());
181  TString omFileName;
182  TString omFileNameReco;
183 
184  // Use HardwareDisplay object to read OnMon maps properly
185  daqchannelmap::HardwareDisplay fHardwareMapping;
186  fHardwareMapping.SetupDet(novadaq::cnv::kNEARDET);
187 
188  // Loop over subruns
189  while(fileListItr >> run >> subrun){
190 
191  if( IsAntiNeutrinoRun(run) ){
192  std::cout << "Running in antineutrino mode!" << std::endl;
193  hits_up = rhc_hits_up;
194  hits_down = rhc_hits_down;
195  numslc_up = rhc_numslc_up;
196  numslc_down = rhc_numslc_down;
197  emptyspillfract_up = rhc_emptyspillfract_up;
198  }
199  //cout << "Processing run: " << run << "/" << subrun << endl;
200 
201  // Find OnMon file for subrun. Use older release if needed
202  int hasOnMonFile = GetFileName(omFileName, run, subrun, true);
203 
204  if(hasOnMonFile == 0){
205  cout << "No file with run number " << run << " and sub run number " << subrun << endl;
206  continue;
207  }
208 
209  // Open OnMon file
210  TFile omFile(omFileName,"read");
211 
212  // Get Header information
213  double startHour, endHour;
214  int Nevents_header;
215  unsigned int Partition;
216  unsigned int startYear, endYear;
217  UShort_t startMonth, endMonth;
218  UShort_t startDay, endDay;
219 
220  TTree *header = (TTree*)omFile.Get("Header");
221  header->SetBranchAddress("Nevents", &Nevents_header);
222  header->SetBranchAddress("StartHour", &startHour);
223  header->SetBranchAddress("StartYear", &startYear);
224  header->SetBranchAddress("StartMonth", &startMonth);
225  header->SetBranchAddress("StartDay", &startDay);
226  header->SetBranchAddress("EndHour", &endHour);
227  header->SetBranchAddress("EndYear", &endYear);
228  header->SetBranchAddress("EndMonth", &endMonth);
229  header->SetBranchAddress("EndDay", &endDay);
230  header->SetBranchAddress("Partition", &Partition);
231 
232  header->GetEntry();
233 
234  // Number of triggers
235  nevents = Nevents_header; // this value is wrong! don't know why...
236 
237  // Get partition number
238  par = Partition;
239 
240  // Get times in seconds. Need to think about time zones.
241  double hour = startHour;
242  int hourInt = int(hour);
243  int minutes = int((hour - hourInt)*60);
244  int seconds = int(((hour-hourInt)*60 - minutes)*60);
245  TDatime startDate;
246  startDate.Set(startYear,startMonth,startDay,hourInt,minutes,seconds);
247 
248  hour = endHour;
249  hourInt = int(hour);
250  minutes = int((hour - hourInt)*60);
251  seconds = int(((hour-hourInt)*60 - minutes)*60);
252  TDatime endDate;
253  endDate.Set(endYear,endMonth,endDay,hourInt,minutes,seconds);
254 
255  firstsec = startDate.Convert();
256  lastsec = endDate.Convert();
257  duration = lastsec-firstsec;
258 
259  delete header;
260 
261  // Find nearline-Ana files.
262  hasReco = true;
263 
264  int hasAnaFile = GetFileName(omFileNameReco, run, subrun, false);
265 
266  if(hasAnaFile == 0){
267  hasReco = false;
268  cout << " No Reco File Found at " << omFileNameReco << endl;
269  }
270 
271  // Initialize PoT
272  pot = -1;
273  preliminary = false;
274 
275  if(hasReco){
276 
277  TFile recoFile(omFileNameReco,"read");
278 
279  // get pot
280  TH1D *fpotsum = (TH1D*)recoFile.Get("nearlineana/fPOTSum");
281  if(fpotsum) {
282  pot = fpotsum->GetBinContent(1);
283  }
284 
285  }
286  else{
287  TDatime tnow;
288  TDatime omtime = omFile.GetCreationDate();
289  if(tnow.Convert() < omtime.Convert() + 24*3600){
290  preliminary = true;
291  cout << "Subrun " << run << "." << subrun
292  << " is less than 24h old and has no Reco File."
293  << " Marking as preliminary..."
294  << endl;
295  }
296  }
297 
298  // Get live time
299  TH1F *RecordedTime = (TH1F*)omFile.Get("RecordedTime");
300 
301  rectimesec = RecordedTime->Integral()*500e-9;
302 
303  delete RecordedTime;
304 
305  // Initialize variables
306  nactivechannels = 0;
307  hitrate = 0;
308 
309  nactivefeb = 0; nactivedb = 0; nactivedcm = 0; ngoodcdb = 0;
310  ngoodpix = 0; ngoodfeb = 0; ngooddb = 0; ngooddcm = 0; ngoodmip = 0;
311  MIPTopSideSym[0] = -9999;
312  MIPTopSideSym[1] = -9999;
313  MIPTopSideSym[2] = -9999;
314  MIPTopSideSym[3] = -9999;
315 
316  vector<float> allMipHitsTop;
317  vector<float> allMipHitsSide;
318  vector<float> allHitRates;
319 
320  int tempdb = 0;
321 
322  // Nested loops to determine good diblock cut
323  for(int db=1; db<=4; db++){
324 
325  int nGoodDCMInDB = 0;
326 
327  vector<float> pixMipHitsTop;
328  vector<float> pixMipHitsSide;
329 
330  int tempdcm = nactivedcm;
331  int nGoodFEBInDCM;
332 
333  int DCMsPerDB;
334 
335  for(int dcm=1; dcm<=4; dcm++){
336 
337  int FEBsPerDB;
338 
339  if(db != 4){
340  DCMsPerDB = 4;
341  if(dcm == 1 || dcm == 4) FEBsPerDB = 32;
342  else FEBsPerDB = 64;
343  }
344  else{
345  DCMsPerDB = 2;
346  if(dcm == 1) FEBsPerDB = 33;
347  else if(dcm == 2) FEBsPerDB = 22;
348  else continue;
349  }
350 
351  TH2F *PixelHitMapMip = (TH2F*)omFile.Get(TString::Format("MipADCPixelsDCM_%02d_%02d",db,dcm));
352  if(!PixelHitMapMip) continue;
353 
354  TH2F *PixelRate = (TH2F*)omFile.Get(TString::Format("PixelsRateDCM_%02d_%02d",db,dcm));
355  if(!PixelRate) continue;
356 
357  nGoodFEBInDCM = 0;
358 
359  int tempfeb = nactivefeb;
360  int nGoodPixInFEB;
361  for(int feb=0; feb<FEBsPerDB; feb++){
362 
363  nGoodPixInFEB = 0;
364 
365  int temppix = nactivechannels;
366 
367  for(int pix=0; pix<32; pix++){
368 
369  // The real coordinates on the onmon plots for a given FEB and Pix
370  // are now (x,y)=(xfeb+xpix,yfeb+ypix)
371  unsigned int xfeb, yfeb, xpix, ypix;
372 
373  nactivechannels++;
374 
375  // Set correct values for xfeb, yfeb, xpix, ypix
376  fHardwareMapping.FEBXY(feb,&xfeb,&yfeb);
377  fHardwareMapping.PixXY(pix,&xpix,&ypix);
378 
379  // Get Bin Content for whichever bin is at (x,y)=(xfeb+xpix,yfeb+ypix)
380  double pixhitrate = PixelRate->GetBinContent(PixelRate->FindFixBin(xfeb+xpix,yfeb+ypix));
381  double pixmiphits = PixelHitMapMip->GetBinContent(PixelRate->FindFixBin(xfeb+xpix,yfeb+ypix));
382  double pixhitratepotscaled = pixhitrate;
383 
384  if(pot >= 0 && nevents > 0){
385  //pixhitratepotscaled += 40*(1 - pot/(4e13*nevents));
386  pixhitratepotscaled += 40*(1 - pot/(4e13*nevents)); // we are now using the same number for FHC and RHC
387  }
388 
389  hitrate += pixhitrate;
390  hitratepotscaled += pixhitratepotscaled;
391 
392  allHitRates.push_back(pixhitratepotscaled);
393 
394  if(db != 4){
395  if(dcm < 3) pixMipHitsTop.push_back(pixmiphits);
396  else pixMipHitsSide.push_back(pixmiphits);
397  }
398  else {
399  if(dcm == 1) pixMipHitsTop.push_back(pixmiphits);
400  else pixMipHitsSide.push_back(pixmiphits);
401  }
402 
403  if(pixhitratepotscaled<pow(10,3.5) && pixhitratepotscaled>pow(10,0.5)){
404  ngoodpix++;
405  nGoodPixInFEB++;
406  }
407  else{
408  std::cout << "FAILED PIXEL HIT OF " << pixhitratepotscaled << " RATE FOR RUN " << run << " SUBRUN " << subrun << " DB " << db << " DCM " << dcm << " FEB " << feb << " PIXEL " << pix << std::endl;
409 
410  if(pixhitratepotscaled < 1 ) std::cout << "NEGATIVE HIT RATE PixHitRate: " << pixhitrate << " POT " << pot << " nevents " << nevents<< std::endl;
411  }
412  }
413 
414  if(temppix < nactivechannels) nactivefeb++;
415 
416  if(nGoodPixInFEB >= MinGoodPixPerFEB){
417  ngoodfeb++;
418  nGoodFEBInDCM++;
419  }
420 
421  }
422 
423  if(tempfeb < nactivefeb) nactivedcm++;
424 
425  // 20% gives you maximum of 12 bad in 64 case and 6 in 32 case.
426  if (nGoodFEBInDCM >= MinGoodFEBFract * FEBsPerDB){
427  ngooddcm++;
428  nGoodDCMInDB++;
429  }
430  else{
431  std::cout << "FAILED NUMBER OF FEBS WITH " << nGoodFEBInDCM <<" IN RUN " << run << " SUBRUN " << subrun << "DCM" << dcm << std::endl;
432  }
433  delete PixelRate;
434  delete PixelHitMapMip;
435 
436  }
437 
438  if(tempdcm < nactivedcm) nactivedb++;
439 
440  allMipHitsTop.insert(allMipHitsTop.end(), pixMipHitsTop.begin(), pixMipHitsTop.end());
441  allMipHitsSide.insert(allMipHitsSide.end(), pixMipHitsSide.begin(), pixMipHitsSide.end());
442 
443  double medianTop = GetMedian(pixMipHitsTop);
444  double medianSide = GetMedian(pixMipHitsSide);
445  double medianTopSide = 0;
446 
447  if(medianTop + medianSide != 0) medianTopSide = (medianTop - medianSide) / (medianTop + medianSide);
448 
449  MIPTopSideSym[db-1]=medianTopSide;
450 
451  float maxAsym = MaxTopSideAsymActive;
452  float minAsym = MinTopSideAsymActive;
453 
454  if(db == 4){
455  maxAsym = MaxTopSideAsymMuonC;
456  minAsym = MinTopSideAsymMuonC;
457  }
458 
459  if(nGoodDCMInDB >= DCMsPerDB){
460  ngooddb++;
461  if(medianTopSide < maxAsym &&
462  medianTopSide > minAsym ){
463  ngoodmip++;
464  tempdb++;
465  }
466  else tempdb = 0;
467  }
468  else tempdb = 0;
469 
470  if(tempdb > ngoodcdb) ngoodcdb = tempdb;
471 
472  }
473 
474  double medianAllTop = GetMedian(allMipHitsTop);
475  double medianAllSide = GetMedian(allMipHitsSide);
476 
477  if(medianAllTop + medianAllSide != 0) mipasym = (medianAllTop - medianAllSide) / (medianAllTop + medianAllSide);
478  else mipasym = 0;
479 
480  if(nactivechannels!=0){
481  hitrate /= nactivechannels;
482  hitratepotscaled /= nactivechannels;
483  }
484  else{
485  hitrate = 0;
486  hitratepotscaled = 0;
487  }
488 
489  midhitrate = GetMedian(allHitRates);
490 
491  allMipHitsTop.insert(allMipHitsTop.end(), allMipHitsSide.begin(), allMipHitsSide.end());
492  midmiprate = GetMedian(allMipHitsTop);
493 
494  if(nevents != 0) midmiprate *= 1000./nevents;
495  else midmiprate = 0;
496 
497  midmipratepotscaled = midmiprate;
498 
499  if(pot > 0) midmipratepotscaled = midmiprate*(2.5e13*nevents/pot);
500  else midmipratepotscaled = -1;
501 
502  miprate = 0;
503 
504  for(int i=0;i<int(allMipHitsTop.size());i++){
505  miprate += allMipHitsTop[i];
506  }
507 
508  double chtime = allMipHitsTop.size() * nevents / 1000.;
509  if(chtime!=0) miprate /= chtime;
510  else miprate = 0;
511 
512  omFile.Close();
513 
514  // Start looking at reconstructed variables
515 
516  // Initialize variables. If no nearline file exists, subrun passes reco cut.
517  timingPeakStart = -9999;
518  timingPeakEnd = -9999;
519  pass_reco = 1;
520  pass_trk3D = 1;
521  pass_slc = 1;
522  corrupted = 0;
523  emptyspillfract = -1;
524  numslc = -1;
525  numslcpotscaled = -1;
526  trkfrac3D = -1;
527  goodbeam = 1.0;
528 
529  // Now Do Reco Processing
530  TDatime procDate;
531  procsec = procDate.Convert(kTRUE);
532 
533  if(hasReco){
534 
535  TFile recoFile(omFileNameReco,"read");
536 
537  // get goodbeam
538  TH1D *fgoodbeam = (TH1D*)recoFile.Get("nearlineana/fGoodBeam");
539  if(fgoodbeam) goodbeam = fgoodbeam->GetMean();
540 
541  //Get Timing Peak
542  //average time per slice
543 
544  TH1D *fTave = (TH1D*)recoFile.Get("nearlineana/fTave");
545  if(fTave){
546 
547  bool crossStart = false;
548  bool crossEnd = false;
549 
550  double bkgd = 0.06 * rectimesec;
551  double decay = 0.9e-16 * pot;
552  double effBkgd = TMath::Max(decay, bkgd);
553  double peak = 90e-16 * pot;
554 
555  double effThreshold = 0;
556  if(peak>0 && effBkgd>0) effThreshold = (peak - effBkgd) / log(peak / effBkgd);
557 
558  for(int iBin=1; iBin<=fTave->GetNbinsX() && effThreshold > 0; iBin++){
559 
560  int slices=fTave->GetBinContent(iBin);
561 
562  if(!crossStart && slices >= effThreshold){
563  timingPeakStart = fTave->GetBinLowEdge(iBin);
564  crossStart = true;
565  }
566 
567  if(crossStart && !crossEnd && slices < effThreshold){
568  timingPeakEnd = fTave->GetBinLowEdge(iBin);
569  crossEnd = true;
570  }
571 
572  if(crossStart && crossEnd && slices >= effThreshold){
573  timingPeakEnd = -9999;
574  crossEnd = false;
575  }
576 
577  }
578 
579  }
580 
581  // Get number of slices / spill
582  TH1D* fNumSlicesRaw = (TH1D*)recoFile.Get("nearlineana/fNumSlices");
583 
584  if(fNumSlicesRaw){
585  numslc = fNumSlicesRaw->GetMean();
586  if(pot > 0) numslcpotscaled = numslc*(2.5e13*nevents)/pot;
587  }
588  else corrupted = 1;
589 
590  // Get number of slices / spiill passing good beam cuts
591  TH1D* fNumSlices = (TH1D*)recoFile.Get("nearlineana/fNumSlicesGoodBeam");
592  if(fNumSlices->GetEntries()>0){
593  float numemptyspills = fNumSlices->GetBinContent(1);
594  if(nevents>0) emptyspillfract = numemptyspills/nevents;
595  }
596  else corrupted = 1;
597 
598  delete fNumSlices;
599 
600  TH1D* fTrackFractionAll3D = (TH1D*)recoFile.Get("nearlineana/fTrackFractionAll3D");
601  if(fTrackFractionAll3D){
602  trkfrac3D = fTrackFractionAll3D->GetMean();
603  }
604  else corrupted = 1;
605  delete fTrackFractionAll3D;
606 
607  recoFile.Close();
608 
609  }
610 
611  // Pass reco cuts if both metrics are good
612  if(apply_emptyspill) passemptyspill = ((emptyspillfract<emptyspillfract_up && emptyspillfract>=0) || !hasReco);
613  else passemptyspill=1;
614 
615  if(apply_timingpeak) pass_timingpeak = (abs(timingPeakStart-NomtimingPeakStart)<=TimingPeakWindow &&
616  abs(timingPeakEnd-NomtimingPeakEnd)<=TimingPeakWindow) || !hasReco;
617  else pass_timingpeak=1;
618 
619  if(apply_slc) pass_slc = ( (numslcpotscaled > numslc_down && numslcpotscaled < numslc_up) || !hasReco);
620  else pass_slc = 1;
621 
622  if(apply_trk) pass_trk3D = ( trkfrac3D > MinTrack3D || !hasReco);
623  else pass_trk3D = 1;
624 
625  if(apply_reco) pass_reco = ( pass_slc && pass_trk3D );
626  else pass_reco = 1;
627 
628  if (apply_hits) pass_hits = ( (midmipratepotscaled < hits_up && midmipratepotscaled > hits_down) || !hasReco);
629  else pass_hits = 1;
630 
631  if (apply_db) pass_db = ( ngoodcdb >= MinGoodDB );
632  else pass_db = 1;
633 
634  if(apply_duration) pass_duration = (nevents > minDuration);
635  else pass_duration=1;
636 
637  if (apply_time) pass_time = ( firstsec < lastsec && startDate.GetYear() >= 2013 && endDate.GetYear() >= 2013);
638  else pass_time = 1;
639 
640  if ( pass_hits && pass_db && pass_time && pass_reco && passemptyspill && pass_timingpeak && pass_duration){
641  pass_all = 1;
642  }
643  else pass_all = 0;
644 
645  if ( !pass_all ){
646 
647  cout << "Run " << run << "." << subrun << " has failed: ";
648  if(!pass_duration) cout << "Number of Triggers = " << nevents << "; ";
649  if(!passemptyspill) cout << "Empty Spill Fraction = " <<emptyspillfract << " " ;
650  if(!pass_hits) cout << "Median PoT scaled MIP/Hit Rate = " << midmipratepotscaled << "/" << midhitrate << "; ";
651  if(!pass_db) cout << "Good MIP/DiBlocks = " << ngoodmip << "/" << ngooddb << "; ";
652  if(!pass_slc) cout << "Num Slices = " << numslcpotscaled << "; ";
653  if(!pass_trk3D) cout << "Fraction of 3D tracks = " << trkfrac3D << "; ";
654  if (!pass_timingpeak) cout << "Timing peak locations " << timingPeakStart * 1e-3 << " us - " << 1e-3 * timingPeakEnd << " us; ";
655  if(!pass_time){
656  cout << "Bad timestamp => Diff = " << lastsec - firstsec << ";" << endl;
657  startDate.Print();
658  endDate.Print();
659  }
660  cout << endl;
661 
662  }
663 
664  mytree->Fill();
665 
666  procfile << TString::Format("%08d\t%03d",run,subrun) << endl;
667 
668  } //End of while loop
669 
670  fileListItr.close();
671 
672  procfile.close();
673 
674  outfile->cd();
675  mytree->Write("mytree");
676  outfile->Close();
677 
678 } //End of Good Runs
679 //////////////////////////////////////////////
681 {
682  //Hardcoded values for antineutrino running.
683  //Easy fix, but possible to switch automatically if horn current is read DIRECTLY from IFBeam_service
684  //database. Not a big deal since we don't suspect we will be switching between the two modes much.
685  //Run numbers for (anti)-neutrino running can be found on productions epoch/period page on the wiki
686  //https://cdcvs.fnal.gov/redmine/projects/novaart/wiki/Period_and_Epoch_Naming
687  //-Ryan Murphy August 2nd, 2016
688  if( (run >11631 && run < 11668) || (run > 11820 && run < 11927) || ( run > 12086 && run < 13131) )
689  return true;
690  else
691  return false;
692 }
693 //////////////////////////////////////////////
694 float GetMedian(vector<float> scores)
695 {
696  float median;
697  size_t size = scores.size();
698 
699  if(size==0) return 0;
700 
701  sort(scores.begin(), scores.end());
702 
703  if (size % 2 == 0)
704  {
705  median = (scores[size / 2 - 1] + scores[size / 2]) / 2;
706  }
707  else
708  {
709  median = scores[size / 2];
710  }
711 
712  return median;
713 }
714 ////////////////////////////////////////////////
715 int GetFileName(TString& filename, int run, int subrun, bool isOnMon){
716 
717  TString fDir;
718  vector<TString> fRel;
719  vector<TString> fTrigDir;
720  vector<TString> fForm;
721  vector<TString> fName;
722 
723  TString fSubDir = TString::Format("/%06d/%08d/",int(run)/100,run);
724 
725  if(isOnMon){
726 
727  fDir = "/nova/data/nearline-OnMon/NearDet/";
728 
729  fRel.push_back("S19-10-22");
730  fRel.push_back("S18-09-02");
731  fRel.push_back("S18-07-03");
732  fRel.push_back("S18-01-05");
733  fRel.push_back("S17-10-10");
734  fRel.push_back("S17-07-18");
735  fRel.push_back("S16-10-26");
736  fRel.push_back("S16-10-07");
737  fRel.push_back("S16-08-23");
738  fRel.push_back("S15-11-06");
739  fRel.push_back("S15-02-05");
740  fRel.push_back("FA14-09-23");
741 
742  fTrigDir.push_back("t00/");
743 
744  fForm.push_back( TString::Format("neardet_nearline_r%08d_s%03d.onmon.root",run,subrun) );
745 
746  fName.push_back(fDir + fRel[0] + fSubDir + fTrigDir[0] + fForm[0]);
747  fName.push_back(fDir + fRel[1] + fSubDir + fTrigDir[0] + fForm[0]);
748  fName.push_back(fDir + fRel[2] + fSubDir + fTrigDir[0] + fForm[0]);
749  fName.push_back(fDir + fRel[3] + fSubDir + fTrigDir[0] + fForm[0]);
750  fName.push_back(fDir + fRel[4] + fSubDir + fTrigDir[0] + fForm[0]);
751  fName.push_back(fDir + fRel[5] + fSubDir + fTrigDir[0] + fForm[0]);
752  fName.push_back(fDir + fRel[6] + fSubDir + fTrigDir[0] + fForm[0]);
753  fName.push_back(fDir + fRel[7] + fSubDir + fTrigDir[0] + fForm[0]);
754  fName.push_back(fDir + fRel[8] + fSubDir + fTrigDir[0] + fForm[0]);
755  fName.push_back(fDir + fRel[9] + fSubDir + fTrigDir[0] + fForm[0]);
756  fName.push_back(fDir + fRel[10] + fSubDir + fTrigDir[0] + fForm[0]);
757  fName.push_back(fDir + fRel[11] + fSubDir + fTrigDir[0] + fForm[0]);
758  }
759  else{
760 
761  fDir = "/nova/data/nearline-Ana/NearDet/";
762 
763  fRel.push_back("S19-10-22");
764  fRel.push_back("S18-09-02");
765  fRel.push_back("S18-07-03");
766  fRel.push_back("S18-01-05");
767  fRel.push_back("S17-10-10");
768  fRel.push_back("S17-07-18");
769  fRel.push_back("S16-10-26");
770  fRel.push_back("S16-10-07");
771  fRel.push_back("S16-08-23");
772  fRel.push_back("S15-11-06");
773  fRel.push_back("S15-02-05");
774  fRel.push_back("S14-12-02");
775  fRel.push_back("FA14-09-23");
776 
777  fTrigDir.push_back("t00/");
778 
779  fForm.push_back( TString::Format("neardet_r%08d_s%02d_reco_hist_t00.root",run,subrun) );
780 
781  fName.push_back(fDir + fRel[0] + fSubDir + fTrigDir[0] + fForm[0]);
782  fName.push_back(fDir + fRel[1] + fSubDir + fTrigDir[0] + fForm[0]);
783  fName.push_back(fDir + fRel[2] + fSubDir + fTrigDir[0] + fForm[0]);
784  fName.push_back(fDir + fRel[3] + fSubDir + fTrigDir[0] + fForm[0]);
785  fName.push_back(fDir + fRel[4] + fSubDir + fTrigDir[0] + fForm[0]);
786  fName.push_back(fDir + fRel[5] + fSubDir + fTrigDir[0] + fForm[0]);
787  fName.push_back(fDir + fRel[6] + fSubDir + fTrigDir[0] + fForm[0]);
788  fName.push_back(fDir + fRel[7] + fSubDir + fTrigDir[0] + fForm[0]);
789  fName.push_back(fDir + fRel[8] + fSubDir + fTrigDir[0] + fForm[0]);
790  fName.push_back(fDir + fRel[9] + fSubDir + fTrigDir[0] + fForm[0]);
791  fName.push_back(fDir + fRel[10] + fSubDir + fTrigDir[0] + fForm[0]);
792  fName.push_back(fDir + fRel[11] + fSubDir + fTrigDir[0] + fForm[0]);
793  fName.push_back(fDir + fRel[12] + fSubDir + fTrigDir[0] + fForm[0]);
794  }
795 
796  int ntries = 0;
797  int filetype = 0;
798 
799  while(ntries < int(fRel.size())){
800  filename = fName[ntries];
801  if(!gSystem->AccessPathName(filename)){
802  filetype = 1;
803  break;
804  }
805  ntries++;
806  }
807 
808  return filetype;
809 
810 
811 }
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)
void BuildMetricsTree_OnMon(TString filelist)
Int_t ngooddcm
Definition: SimpleIterate.C:30
Int_t pass_all
Definition: SimpleIterate.C:42
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
Double_t miprate
Definition: SimpleIterate.C:48
Double_t rectimesec
Definition: SimpleIterate.C:44
Int_t nevents
Definition: SimpleIterate.C:25
float abs(float number)
Definition: d0nt_math.hpp:39
::xsd::cxx::tree::duration< char, simple_type > duration
Definition: Database.h:188
Int_t nactivefeb
Definition: SimpleIterate.C:34
Double_t midhitrate
Definition: SimpleIterate.C:47
Int_t nactivedcm
Definition: SimpleIterate.C:35
#define pot
Near Detector in the NuMI cavern.
Int_t ngoodmip
Definition: SimpleIterate.C:32
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
Int_t firstsec
Definition: SimpleIterate.C:26
void FEBXY(unsigned int feb, unsigned int *ix, unsigned int *iy)
ifstream fileListItr("ManBadSubruns.txt")
Definition: decay.py:1
Int_t pass_db
Definition: SimpleIterate.C:40
ofstream procfile
Double_t hitrate
Definition: SimpleIterate.C:46
Int_t pass_slc
Definition: SimpleIterate.C:57
bool IsAntiNeutrinoRun(int run)
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 corrupted
Definition: SimpleIterate.C:59
FILE * outfile
Definition: dump_event.C:13
Int_t pass_time
Definition: SimpleIterate.C:41
void PixXY(unsigned int pix, unsigned int *ix, unsigned int *iy)
Int_t pass_hits
Definition: SimpleIterate.C:39