BuildMetricsTree_OnMon.C

Go to the documentation of this file.

#include <iostream>
#include <fstream>
#include <iomanip>
#include <sstream>
#include <string>
#include <streambuf>
#include <cstring>
#include <algorithm>

#include <TStyle.h>
#include <TSystem.h>
#include <TFile.h>
#include <TString.h>
#include <TTree.h>
#include <TH2F.h>
#include <TH1F.h>
#include <TMath.h>

#include <TSQLServer.h>
#include <TSQLResult.h>
#include <TSQLStatement.h>

using namespace std;

struct RunInfo{

  int run;
  bool isLast;
  bool isIncomplete;
  vector< double > livetime;
  vector< int > subrun;
  vector< TString > fnames;
  vector< TString > rnames;
  vector< bool > isNewRelom;
  vector< bool > isNewRelana;

};

double GetMedian(TH1D* hist);

float GetMedian(vector<float> scores);

bool GetDBPartition(int run, int &par);

int GetFileName(TString& filename, int run, int subrun, bool isOnMon = true);

RunInfo GetRunInfo(int run);

vector< RunInfo > GetSubrunSets(RunInfo wholerun, double run_len);

void BuildMetricsTree_OnMon(TString filelist){

  gStyle->SetTimeOffset(0);

  Int_t     apply_runlen = 1;     //1 if want to apply live time cut;      0 if do not 
  Int_t     apply_empty  = 1;     //1 if want to apply empty cut;          0 if do not 
  Int_t     apply_time   = 1;     //1 if want to apply negative time cut;  0 if do not (start time before end time) 
  Int_t     apply_db     = 1;     //1 if want to apply diblock cut;        0 if do not
  Int_t     apply_hits   = 1;     //1 if want to apply hit rate cut;       0 if do not
  Int_t     apply_reco   = 1;     //1 if want to apply reco cut;           0 if do not

  //Cut Parameter Values

  Double_t run_len          = 1.00;  //Minimum live time of run in seconds. 
  Double_t hits_up          = 23.0;  //Median MIP rate smaller than this parameter
  Double_t hits_down        = 13.0;  //Median MIP rate larger than this parameter
  Double_t PixelsPerFEB     = 26;    //Minimum # of good pixels per FEB
  Double_t FEBsPerDCM       = 56;    //Minimum # of good FEBs per DCM
  Double_t DCMsPerDB        = 12;    //Minimum # of good DCMs per DiBlock
  Double_t MaxTopSideAsym   = 0.5;   //Maximum Top/Side MIP hit rate Asymmetry
  Double_t MinTopSideAsym   = 0.1;   //Minimum Top/Side MIP hit rate Asymmetry
  Double_t MinGoodDB        = 2;     //Minimum # of good consecutive DiBlocks
  Double_t empty_cut        = 100.0; //The percentage of empty events less than this value. This only cuts out files with ONLY empty events.
  Double_t MaxTrkFrac2D     = 0.15;  //Maximum 2D track fraction
  Double_t MinNumSlc        = 1.2;   //Minimum # of slices / trigger / 10^4 channels
  Double_t MaxNumSlc        = 3.2;   //Maximum # of slices / trigger / 10^4 channels

  // MIP rate cuts. Entries are for DCMs {01 - 06 (Top), 07, 08, 09, 10, 11, 12}

  float lowpixmip[7]  = {13,  5,  4,  4,  2,  2, 1.5}; //Minimum MIP rate for pixels in different DCMs
  float highpixmip[7] = {31, 45, 36, 30, 26, 23, 20.}; //Maximum MIP rate for pixels in different DCMs

  float lowfebmip[7]  = {400,  300,  200, 200, 100, 100, 100};  //Minimum MIP rate for FEBs in different DCMs
  float highfebmip[7] = {1000, 1200, 1000, 850, 700, 600, 550}; //Maximum MIP rate for FEBs in different DCMs
  

  // **************END OF TYPICAL USER'S NEED TO ALTER THIS SCRIPT*******************************************

  // Open file to keep track of processed subruns
  ofstream procfile;
  procfile.open("listProcessedTemp.txt");

  // Variables to fill Good Run metrics
  Int_t run, subrun, par, nevents, pass_runlen, pass_hits, pass_all, pass_db, pass_time;
  Int_t pass_empty, firstsec, lastsec, ngoodpix, ngoodfeb, ngooddcm, ngooddb, ngoodmip;
  Int_t nactivefeb, nactivedcm, nactivedb, ngoodcdb, dbencoded, dbaencoded, setsize;

  Double_t hitrate, nactivechannels, rectimesec, emptypercentage, setlivetime, 
           mipratio, midhitrate, midmiprate, mipasym, miprate;

  Int_t pass_reco, pass_slc, pass_trk, corrupted, procsec;
  Double_t numslc, trkfrac2D;

  Int_t preliminary;

  // Create root file to store tree
  TFile *outfile = new TFile("MetricsNew.root","recreate");

  // Create tree to store Good Run metrics
  TTree *mytree = new TTree("evts","");
  mytree->Branch("run",&run,"run/I");
  mytree->Branch("subrun",&subrun,"subrun/I");
  mytree->Branch("par",&par,"par/I");
  mytree->Branch("nevents",&nevents,"nevents/I");
  mytree->Branch("firstsec",&firstsec,"firstsec/I");
  mytree->Branch("lastsec",&lastsec,"lastsec/I");
  mytree->Branch("ngoodpix",&ngoodpix,"ngoodpix/I");
  mytree->Branch("ngoodfeb",&ngoodfeb,"ngoodfeb/I");
  mytree->Branch("ngooddcm",&ngooddcm,"ngooddcm/I");
  mytree->Branch("ngooddb",&ngooddb,"ngooddb/I");
  mytree->Branch("ngoodmip",&ngoodmip,"ngoodmip/I");
  mytree->Branch("ngoodcdb",&ngoodcdb,"ngoodcdb/I");
  mytree->Branch("nactivefeb",&nactivefeb,"nactivefeb/I");
  mytree->Branch("nactivedcm",&nactivedcm,"nactivedcm/I");
  mytree->Branch("nactivedb",&nactivedb,"nactivedb/I");
  mytree->Branch("pass_runlen",&pass_runlen,"pass_runlen/I");
  mytree->Branch("pass_empty",&pass_empty,"pass_empty/I");
  mytree->Branch("pass_hits",&pass_hits,"pass_hits/I");
  mytree->Branch("pass_db",&pass_db,"pass_db/I");
  mytree->Branch("pass_time",&pass_time,"pass_time/I");
  mytree->Branch("pass_all",&pass_all,"pass_all/I");
  mytree->Branch("setsize",&setsize,"setsize/I");
  mytree->Branch("rectimesec",&rectimesec,"rectimesec/D");
  mytree->Branch("setlivetime",&setlivetime,"setlivetime/D");
  mytree->Branch("hitrate",&hitrate,"hitrate/D");
  mytree->Branch("midhitrate",&midhitrate,"midhitrate/D");
  mytree->Branch("miprate",&miprate,"miprate/D");
  mytree->Branch("midmiprate",&midmiprate,"midmiprate/D");
  mytree->Branch("mipratio",&mipratio,"mipratio/D");
  mytree->Branch("mipasym",&mipasym,"mipasym/D");
  mytree->Branch("nactivechannels",&nactivechannels,"nactivechannels/D");
  mytree->Branch("emptypercentage",&emptypercentage,"emptypercentage/D");

  mytree->Branch("dbencoded",&dbencoded,"dbencoded/I");
  mytree->Branch("dbaencoded",&dbaencoded,"dbaencoded/I");

  mytree->Branch("pass_reco",&pass_reco,"pass_reco/I");
  mytree->Branch("pass_slc",&pass_slc,"pass_slc/I");
  mytree->Branch("pass_trk",&pass_trk,"pass_trk/I");
  mytree->Branch("corrupted",&corrupted,"corrupted/I");
  mytree->Branch("procsec",&procsec,"procsec/I");
  mytree->Branch("numslc",&numslc,"numslc/D");
  mytree->Branch("trkfrac2D",&trkfrac2D,"trkfrac2D/D");

  mytree->Branch("preliminary",&preliminary,"preliminary/I");

  // Open list of subruns to process
  ifstream fileListItr(filelist.Data());

  // Loop over subruns
  while(fileListItr >> run){

    // Get a vector of subrun sets
    vector< RunInfo > setvec = GetSubrunSets( GetRunInfo(run), run_len );

    // If vector is empty, fill tree with empty set
    if(setvec.size()==0){

      subrun = 0; par = 4; nevents = 0; pass_runlen = 0; pass_hits = 0;
      pass_all = 0; pass_db = 0; pass_time = 0; pass_empty = 0;
      firstsec = 0; lastsec = 0; ngoodpix = 0; ngoodfeb = 0; 
      ngooddcm = 0; ngooddb = 0; ngoodmip = 0; nactivefeb = 0;
      nactivedcm = 0; nactivedb = 0; ngoodcdb = 0; dbencoded = 0;
      dbaencoded = 0; setsize = 0;
      
      hitrate = 0; nactivechannels = 0; rectimesec = 0; emptypercentage = 0;
      setlivetime = 0; mipratio = 0; midhitrate = 0; midmiprate = 0; 
      mipasym = 0; miprate =0;

      pass_reco = 0; pass_slc = 0; pass_trk = 0; corrupted = 0; procsec = 0;
      numslc = 0; trkfrac2D = 0;

      preliminary = 0;

      mytree->Fill();

    }

    // Loop over subrun sets
    for(int seti=0; seti<int(setvec.size()); seti++){

      // 2D vector of Pixel MIP hits maps (DB x DCM)
      vector< vector< TH2F* > > TotPixMapMip(14, vector< TH2F* >(12,0));

      // Reset live time of set
      setlivetime = 0;

      // First loop over subruns in set to add Mip hits 
      for(int subruni=0; subruni<int(setvec[seti].subrun.size()); subruni++){

        // Open OnMon file
        TFile omFile(setvec[seti].fnames[subruni],"read");

        // Add MIP hits for whole run
        for(int db=1;db<=14;db++){
        for(int dcm=1;dcm<=12;dcm++){

          // Temporary histogram to store maps
          TH2F *h2temp = 0;

          // Get map from OnMon file
          h2temp = (TH2F*)omFile.Get(TString::Format("MipADCPixelsDCM_%02d_%02d",db,dcm));
          
          // Check that map exists in file, otherwise continue
          if(!h2temp){
           delete h2temp;
           continue;
          }
          
          // If it's the first map, clone it, else add histograms
          if(!TotPixMapMip[db-1][dcm-1]){
            TotPixMapMip[db-1][dcm-1] = (TH2F*)h2temp->Clone("");
          }
          else{
            TotPixMapMip[db-1][dcm-1]->Add(h2temp);
          }

          // Delete temporary histogram
          delete h2temp;

        }}

        // Close OnMon File
        omFile.Clear();
        omFile.Close();
        
        // Add the live time for this subrun
        setlivetime += setvec[seti].livetime[subruni];

      }// End of first subruns loop

      // Fill run and preliminary status
      run = setvec[seti].run;
      preliminary = setvec[seti].isIncomplete;

      // Now loop over subruns again to apply cuts
      for(int subruni=0; subruni<int(setvec[seti].subrun.size()); subruni++){

        // Fill subrun number and set size
        subrun = setvec[seti].subrun[subruni];
        setsize = setvec[seti].subrun.size();

        // Fill single subrun live time
        rectimesec = setvec[seti].livetime[subruni]; 

        bool isNewRelom = setvec[seti].isNewRelom[subruni];

        // Get Header information
        double startHour, endHour;
        Int_t Nevents;
        UInt_t Partition;
        UInt_t startYear, endYear;
        UShort_t startMonth, endMonth;
        UShort_t startDay, endDay;
  
        TFile omFile(setvec[seti].fnames[subruni],"read");
  
        TTree *header = (TTree*)omFile.Get("Header");
        header->SetBranchAddress("Nevents",&Nevents);
        header->SetBranchAddress("StartHour",&startHour);
        header->SetBranchAddress("StartYear",&startYear);
        header->SetBranchAddress("StartMonth",&startMonth);
        header->SetBranchAddress("StartDay",&startDay);
        header->SetBranchAddress("EndHour",&endHour);
        header->SetBranchAddress("EndYear",&endYear);
        header->SetBranchAddress("EndMonth",&endMonth);
        header->SetBranchAddress("EndDay",&endDay);
        // Partition info only exists in newer releases
        if(isNewRelom) header->SetBranchAddress("Partition",&Partition);
  
        header->GetEntry(0);
  
        // Number of triggers
        nevents = Nevents;
  
        // Get times in seconds. Need to think about time zones.
        double hour = startHour;
        int hourInt = int(hour);
        int minutes = int((hour - hourInt)*60);
        int seconds = int(((hour-hourInt)*60 - minutes)*60);
        TDatime startDate;
        if(isNewRelom) startDate.Set(startYear,startMonth,startDay,hourInt,minutes,seconds);
        else         startDate.Set(startYear+1900,startMonth+1,startDay,hourInt,minutes,seconds);
  
        hour = endHour;
        hourInt = int(hour);
        minutes = int((hour - hourInt)*60);
        seconds = int(((hour-hourInt)*60 - minutes)*60);
        TDatime endDate;
        if(isNewRelom) endDate.Set(endYear,endMonth,endDay,hourInt,minutes,seconds);
        else         endDate.Set(endYear+1900,endMonth+1,endDay,hourInt,minutes,seconds);
  
        firstsec = startDate.Convert();
        lastsec = endDate.Convert();
  
        delete header;
  
        bool connected = false;
  
        // Get partition number
        if(isNewRelom){
          par = Partition;
          connected = true;
        }
        // Old releases need to use the database
        else{
          connected = GetDBPartition(run,par);
        }

        // Skip if database connection failed  
        if(!connected) continue;
  
        // Get number of empty events
        TH2F *NhitVsHour = (TH2F*)omFile.Get("NhitVsHour");
  
        int nemptyevents = NhitVsHour->Integral(61,1500,1,1); 
        int nallevents = NhitVsHour->Integral(61,1500); 
        if(nallevents!=0) emptypercentage  = 100.0*nemptyevents/nallevents;
        else              emptypercentage  = 0;
  
        delete NhitVsHour;
      
        // Initialize variables
        nactivechannels = 0;
        hitrate = 0;
  
        nactivefeb = 0; nactivedb = 0; nactivedcm = 0; ngoodcdb = 0;
        ngoodpix = 0; ngoodfeb = 0; ngooddb = 0; ngooddcm = 0; ngoodmip = 0;
        dbencoded = 0; dbaencoded = 0;
  
        // Vectors to store all pixel hits
        vector<float> allMipHitsTop;
        vector<float> allMipHitsSide;
        vector<float> allHitRates;
  
        // Nested loops to determine good diblock cut 

        // This is to keep track of good consecutive diblocks
        int tempdb = 0;
  
        // Loop over diblocks
        for(int db=1;db<=14;db++){
  
          // Count good dcms in this db
          int nGoodDCMInDB = 0;
  
          // Vectors to store pixel hits for this db
          vector<float> pixMipHitsTop;
          vector<float> pixMipHitsSide;

          // Keep track of active dcms  
          int tempdcm = nactivedcm;
          
          // Loop over DCMs
          for(int dcm=1;dcm<=12;dcm++){
          
            // Check if MIP hit map exists for this DCM
            if(!TotPixMapMip[db-1][dcm-1]) continue;
  
            // Get overall hit rate map for this DCM
            TH2F *PixelRate = (TH2F*)omFile.Get(TString::Format("PixelsRateDCM_%02d_%02d",db,dcm));
            if(!PixelRate) continue; // Skip if it doesn't exist
  
            // Count good FEBs in this DCM
            int nGoodFEBInDCM = 0;
  
            // Keep track of active FEBs
            int tempfeb = nactivefeb;
  
            // Loop over FEBs
            for(int feb=0;feb<64;feb++){
  
              // Count good Pixels in this FEB
              int nGoodPixInFEB = 0;
  
              // Keep track of active Pixels
              int temppix = nactivechannels;
              
              // Sum of MIP hits in this FEB
              double febmiphits = 0;
  
              // Loop over Pixels
              for(int pix=0;pix<32;pix++){
  
                // Figure out bins x and y for this pixel
                int ix = (feb%16)*4 + pix%4 + 1;
                int iy = int(feb/16)*8 + int(pix/4) + 1;
  
                // Check this pixel has MIP hits, i.e. it is active
                if( TotPixMapMip[db-1][dcm-1]->GetBinContent(ix,iy)==0 ) continue;
  
                // Increment active channels
                nactivechannels++;
  
                // Store overall hit rate and MIP hits in this Pixel
                double pixhitrate = PixelRate->GetBinContent(ix,iy);
                double pixmiphits = TotPixMapMip[db-1][dcm-1]->GetBinContent(ix,iy);
                
                // Increment MIP hits in this FEB with a max threshold of 4x the pixel high MIP rate cut
                // The threshold mitigates the effect of few pixels with very large rates 
                febmiphits += TMath::Min(pixmiphits, highpixmip[TMath::Max(dcm-6,0)] * setlivetime * 4);
  
                // Fill total hit rate variable and vectors
                hitrate += pixhitrate;
                allHitRates.push_back(pixhitrate);
                if(dcm<7) pixMipHitsTop.push_back(pixmiphits);
                else      pixMipHitsSide.push_back(pixmiphits);
  
                // Determine if pixel is good, i.e. neither cold nor hot
                if(pixhitrate<pow(10,3.5) && pixhitrate>pow(10,0.5) &&
                   TMath::Gamma(lowpixmip[TMath::Max(dcm-6,0)] * setlivetime , pixmiphits) > 0.1 &&
                   TMath::Gamma(highpixmip[TMath::Max(dcm-6,0)] * setlivetime , pixmiphits) < 0.9 ){

                  // Increment good pixels in total and in this FEB
                  ngoodpix++;
                  nGoodPixInFEB++;

                }

              }// End of Pixel loop
  
              // If active channels has increased, this FEB is active.
              // Increment active FEBs 
              if(temppix < nactivechannels) nactivefeb++;
  
              // Determine if FEB is good, i.e. good MIP rate and enough good pixels
              if(nGoodPixInFEB >= PixelsPerFEB &&
                 TMath::Gamma(lowfebmip[TMath::Max(dcm-6,0)] * setlivetime , febmiphits) > 0.1 &&
                 TMath::Gamma(highfebmip[TMath::Max(dcm-6,0)] * setlivetime , febmiphits) < 0.9 ){

                // Increment good FEBs in total and in this DCM
                ngoodfeb++;
                nGoodFEBInDCM++;

              }

            }// End of FEB loop
  
            // If active FEBs has increased, this DCM is active.
            // Increment active DCMs 
            if(tempfeb < nactivefeb) nactivedcm++;
  
            // Determine if DCM is good, i.e. enough good FEBs
            if(nGoodFEBInDCM >= FEBsPerDCM){
              // Increment good DCMs in total and in this diblock
              ngooddcm++;
              nGoodDCMInDB++;
            }
  
            // Delete pixel rate map
            delete PixelRate;
  
          }// End of DCM loop
  
          // If active DCMs has increased, this diblock is active.
          // Increment active diblocks 
          if(tempdcm < nactivedcm){
            nactivedb++;
            // Also, fill active diblock mask
            dbaencoded += pow(2,db-1);
          }
  
          // Fill vectors for whole detector with hits from this diblock
          allMipHitsTop.insert(allMipHitsTop.end(), pixMipHitsTop.begin(), pixMipHitsTop.end());
          allMipHitsSide.insert(allMipHitsSide.end(), pixMipHitsSide.begin(), pixMipHitsSide.end());
  
          // Find median of MIP hits for top and side of this diblock
          double medianTop  = GetMedian(pixMipHitsTop);
          double medianSide = GetMedian(pixMipHitsSide);
  
          // Calculate top-side asymmetry if there are MIP hits
          if(medianTop + medianSide != 0) medianTop = (medianTop - medianSide) / (medianTop + medianSide);
          else                            medianTop = 0;
  
          // Determine if diblock is good
          if(nGoodDCMInDB >= DCMsPerDB){
            // Increment good diblocks (DCMs statuses only)
            ngooddb++;
            // Check for bad top-side asymmetry
            if(medianTop < MaxTopSideAsym && 
               medianTop > MinTopSideAsym ){
              // Increment good diblocks (including asymmetry)
              ngoodmip++;
              // Counting consecutive good diblocks
              tempdb++;
              // Fill good diblock mask
              dbencoded += pow(2,db-1);
            } 
            // Else, reset consecutive diblocks 
            else tempdb = 0;
          }
          // Else, reset consecutive diblocks 
          else tempdb = 0;
  
          // Fill largest number of consecutive diblocks
          if(tempdb > ngoodcdb) ngoodcdb = tempdb;
  
        }// End of diblock loop
  
  
        // Compute median of all MIP hits for top and side
        double medianAllTop  = GetMedian(allMipHitsTop);
        double medianAllSide = GetMedian(allMipHitsSide);
  
        // Calculate overall top-side asymmetry
        if(medianAllTop + medianAllSide != 0) mipasym = (medianAllTop - medianAllSide) / (medianAllTop + medianAllSide);
        else                                  mipasym = 0;
  
        // Calculate top-side ratio
        if(medianAllSide!=0) mipratio = medianAllTop / medianAllSide;
        else                 mipratio = -1;
  
        // Normalise hit rate per active channel (mean)
        if(nactivechannels!=0) hitrate /= nactivechannels;
        else                   hitrate = 0;
  
        // Get overall median of hit rates
        midhitrate = GetMedian(allHitRates);
  
        // Add top and side MIP hits and compute median
        allMipHitsTop.insert(allMipHitsTop.end(), allMipHitsSide.begin(), allMipHitsSide.end());
        midmiprate = GetMedian(allMipHitsTop);
        if(setlivetime != 0) midmiprate /= setlivetime;
        else                midmiprate = 0; 
  
        // Reset mean MIP rate
        miprate = 0;
  
        // Add all MIP hits
        for(int i=0;i<int(allMipHitsTop.size());i++){
          miprate += allMipHitsTop[i];
        }
  
        // Compute mean MIP rate
        double chtime = allMipHitsTop.size() * setlivetime;
        if(chtime!=0) miprate /= chtime;
        else          miprate = 0;
  
        // Close OnMon file
        omFile.Clear();
        omFile.Close();
  

        // Start looking at reconstructed variables
  
        // Initialize variables. If no nearline file exists, subrun passes reco cut. 
        pass_reco = 1;
        pass_trk = 1;
        pass_slc = 1;
        corrupted = 0;
        numslc = -1;
        trkfrac2D = -1;
        
        // Get file name and release type
        TString anaFileName = setvec[seti].rnames[subruni];
        bool isNewRelana = setvec[seti].isNewRelana[subruni];
        
        // Check if file exists
        bool hasReco = true;
        if(anaFileName == "") hasReco = false;
  
        // Fill timestamp in which this tree was filled
        TDatime procDate;
        procsec = procDate.Convert(kTRUE);
  
        // If file exists, get reco info
        if(hasReco){
  
          // Open Ana file
          TFile recoFile(anaFileName,"read");
  
          // Get number of slices / trigger
          TH1D* fNumSlices = (TH1D*)recoFile.Get("nearlineana/fNumSlices");
          // Check if histogram exists
          if(fNumSlices){
            // Make sure there are active channels
            if(nactivechannels!=0){
              // Fill slicing rate normalised to 10k active channels
              numslc = 1e4 * fNumSlices->GetMean() / nactivechannels;
              // Fill pass slicing rate cut
              pass_slc = ((MinNumSlc < numslc) && (numslc < MaxNumSlc));
            }
          }
          // If histogram failed, file is corrupted
          else corrupted = 1;
          //Delete histogram
          delete fNumSlices;
          
          // If this is a new release, also get fraction of tracks that are 2D
          if(isNewRelana){
            // Get 2D track fraction histogram
            TH1D* fTrackFractionAll2D = (TH1D*)recoFile.Get("nearlineana/fTrackFractionAll2D");
            // Check histogram exists
            if(fTrackFractionAll2D){
              // Overall fraction is just the mean
              trkfrac2D = fTrackFractionAll2D->GetMean();
              // Fill pass tracking rate cut
              pass_trk = trkfrac2D < MaxTrkFrac2D;
            }
            // If histogram failed, file is corrupted
            else corrupted = 1;
            //Delete histogram
            delete fTrackFractionAll2D;
          }
          
          // Close Ana file
          recoFile.Clear();
          recoFile.Close();
  
        }// End of reco if
  
  
        // Apply all cuts
  
        // Pass reco cuts if both metrics are good
        if(apply_reco) pass_reco = ( pass_trk && pass_slc );
        else pass_reco = 1;
  
        if (apply_runlen) pass_runlen = ( setlivetime > run_len );
        else pass_runlen = 1;
  
        if (apply_empty) pass_empty = ( emptypercentage < empty_cut );
        else pass_empty = 1;
  
        if (apply_hits) pass_hits = ( midmiprate < hits_up && midmiprate > hits_down );
        else pass_hits = 1;
  
        // Partition 2 requires only 1 good diblock
        if (apply_db) pass_db = ( (ngoodcdb >= MinGoodDB) || (par == 2 && ngoodcdb > 0) );
        else pass_db = 1;
  
        if (apply_time) pass_time = ( firstsec < lastsec && startDate.GetYear() >= 2013 && endDate.GetYear() >= 2013);
        else pass_time = 1;
  
        // Fill variable for passing all cuts
        if ( pass_runlen && pass_hits && pass_db && pass_time && pass_empty && pass_reco ){ 
  
          pass_all = 1;
  
        }
        else pass_all = 0;
  
        // If subrun is bad, cout more info
        if ( !pass_all ){ 
  
          cout << "Run " << run << "." << subrun << " has failed: ";
          if(!pass_runlen) cout << "Run Live Time = " << setlivetime << " seconds; ";
          if(!pass_empty)  cout << "Empty Percentage = " << emptypercentage << "%; ";
          if(!pass_hits)   cout << "Median MIP/Hit Rate = " << midmiprate << "/" << midhitrate << "; ";
          if(!pass_db)     cout << "Good MIP/DiBlocks = " << ngoodmip << "/" << ngooddb << "; ";
          if(!pass_slc)    cout << "Num Slices / 1e4 Pixels = " << numslc << "; ";
          if(!pass_trk)    cout << "Fraction of 2D tracks = " << trkfrac2D << "; ";
          if(!pass_time){   
            cout << "Bad timing => Diff = " << lastsec - firstsec << ";" << endl; 
            startDate.Print();
            endDate.Print();
          }
          cout << endl;
  
        }
  
        // Fill tree
        mytree->Fill();
  
      } //End of subrun loop
  
    } //End of set loop
  
    // Write this run to processed list
    procfile << TString::Format("%08d",run) << endl;
  
  } //End of while loop over run list

  // Close run list file
  fileListItr.close();
  
  // Close processed runs file
  procfile.close();
  
  // Write tree to root file
  outfile->cd();
  mytree->Write("mytree");
  outfile->Close();
  

} //End of Good Runs

///////////////////////////////////////////////////////////////

int GetFileName(TString& filename, int run, int subrun, bool isOnMon){

  TString fDir;
  vector<TString> fRel;
  vector<TString> fTrigDir;
  vector<TString> fForm;
  vector<TString> fName;
  vector<int> fRelType;
  
  TString fSubDir = TString::Format("/%06d/%08d/",int(run)/100,run);
           
  if(isOnMon){

    fDir = "/nova/data/nearline-OnMon/FarDet/";

    fRel.push_back("S15-02-05");
    fRel.push_back("FA14-09-23");
    fRel.push_back("S14-08-01");
    fRel.push_back("S14-02-05");
    fRel.push_back("S13-09-17");

    fTrigDir.push_back("");
    fTrigDir.push_back("t02-t00/"); 

    fForm.push_back( TString::Format("FDNL-onmon-endsubrun-%08d-%03d.root",run,subrun) );
    fForm.push_back( TString::Format("fardet_nearline_r%08d_s%03d.onmon.root",run,subrun) );

    fName.push_back(fDir + fRel[0] + fSubDir + fTrigDir[1] + fForm[1]);
    fName.push_back(fDir + fRel[1] + fSubDir + fTrigDir[1] + fForm[1]);
    fName.push_back(fDir + fRel[2] + fSubDir + fTrigDir[0] + fForm[1]);
    fName.push_back(fDir + fRel[3] + fSubDir + fTrigDir[0] + fForm[1]);
    fName.push_back(fDir + fRel[4] + fSubDir + fTrigDir[0] + fForm[0]);
    
    fRelType.push_back(2);
    fRelType.push_back(2);
    fRelType.push_back(2);
    fRelType.push_back(2);
    fRelType.push_back(1);
  
  }
  else{
  
    fDir = "/nova/data/nearline-Ana/FarDet/";

    fRel.push_back("S15-02-05");
    fRel.push_back("FA14-09-23");
    fRel.push_back("S14-08-01");
    fRel.push_back("S14-02-24");
    fRel.push_back("S14-01-20");

    fTrigDir.push_back("");
    fTrigDir.push_back("t02/"); 

    fForm.push_back( TString::Format("fardet_r%08d_s%02d_reco_hist_t02.root",run,subrun) );

    fName.push_back(fDir + fRel[0] + fSubDir + fTrigDir[1] + fForm[0]);
    fName.push_back(fDir + fRel[1] + fSubDir + fTrigDir[1] + fForm[0]);
    fName.push_back(fDir + fRel[2] + fSubDir + fTrigDir[0] + fForm[0]);
    fName.push_back(fDir + fRel[3] + fSubDir + fTrigDir[0] + fForm[0]);
    fName.push_back(fDir + fRel[4] + fSubDir + fTrigDir[0] + fForm[0]);

    fRelType.push_back(2);
    fRelType.push_back(2);
    fRelType.push_back(2);
    fRelType.push_back(2);
    fRelType.push_back(1);

  }

  int ntries = 0;
  int filetype = 0;
  
  while(!filetype && ntries < int(fRel.size())){

    filename = fName[ntries];
    
    if(!gSystem->AccessPathName(filename)){ 
    
      filetype = fRelType[ntries];
    
    }
    
    ntries++;
    
  }
  
  return filetype;

}

///////////////////////////////////////////////////////////////

RunInfo GetRunInfo(int run){

    RunInfo myrun;
    
    myrun.run = run;

    bool incomplete = false;
    
    for(int subrun=0; subrun<64; subrun++){

      // Find OnMon file for subrun. Use older release if needed

      TString omFileName;

      int omfiletype = GetFileName(omFileName, run, subrun, true);

      bool isNewRelom = false;
      
      if(omfiletype == 0) continue;
      else if(omfiletype > 1) isNewRelom = true;

      ////////////////////////////////////////////////////////////////
      // HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK//
      // HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK//
      // HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK//
      /* HACK */  if((run == 17647 && subrun > 25)           // HACK//
      /* HACK */  || (run == 18531 && subrun == 13))         // HACK//
      /* HACK */  continue;                                  // HACK//
      // HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK//
      // HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK//
      // HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK HACK//
      ////////////////////////////////////////////////////////////////

      TFile omFile(omFileName,"read");

      TDatime tnow;
      TDatime tfile = omFile.GetCreationDate();
    
      if(tnow.Convert() < tfile.Convert() + 4500){
        incomplete = true;
      }

      // Get live time
      TH1F *RecordedTime = (TH1F*)omFile.Get("RecordedTime");

      if(RecordedTime) myrun.livetime.push_back(RecordedTime->Integral()*500e-9);
      else{
        incomplete = true;
        delete RecordedTime;
        continue;
      }                   
                  
      delete RecordedTime;

      omFile.Clear();
      omFile.Close();

      TString anaFileName;

      int anafiletype = GetFileName(anaFileName, run, subrun, false);

      bool isNewRelana = false;
      
      if(anafiletype == 0){
        anaFileName = "";
        if(tnow.Convert() < tfile.Convert() + 3600*24){
          incomplete = true;
        }
      }
      else if(anafiletype > 1) isNewRelana = true;
      
      myrun.subrun.push_back(subrun);
      myrun.fnames.push_back(omFileName);
      myrun.rnames.push_back(anaFileName);
      myrun.isNewRelom.push_back(isNewRelom);
      myrun.isNewRelana.push_back(isNewRelana);

    }// End of subrun loop

    if(incomplete){
      cout << "Run " << run << " is incomplete. Marking as preliminary..." << endl;
    }

    myrun.isIncomplete = incomplete;
    
    return myrun;

}

///////////////////////////////////////////////////////////////

vector< RunInfo > GetSubrunSets(RunInfo wholerun, double run_len){

    int nsubruns = wholerun.subrun.size();
    double remlt = 0;
    for(int subruni=0; subruni<nsubruns; subruni++){
      remlt += wholerun.livetime[subruni];
    }

    double setlt = 0;
    int nsets = 0;

    vector< RunInfo > setvec;

    RunInfo emptystruct;
    RunInfo thisset;

    for(int subruni=0; subruni<nsubruns; subruni++){

      thisset.subrun.push_back(wholerun.subrun[subruni]);
      thisset.livetime.push_back(wholerun.livetime[subruni]);
      thisset.fnames.push_back(wholerun.fnames[subruni]);
      thisset.rnames.push_back(wholerun.rnames[subruni]);
      thisset.isNewRelom.push_back(wholerun.isNewRelom[subruni]);
      thisset.isNewRelana.push_back(wholerun.isNewRelana[subruni]);

      setlt += thisset.livetime.back();
      remlt -= thisset.livetime.back();

      if((setlt > run_len && remlt > run_len) || subruni == nsubruns-1){
        if(subruni<nsubruns-1) thisset.isLast = false;
        else                   thisset.isLast = true;
        thisset.run = wholerun.run;
        thisset.isIncomplete = wholerun.isIncomplete;
        setvec.push_back(thisset);
        thisset = emptystruct;
        setlt = 0;
        nsets++;
      }

    }

    cout << TString::Format("Run %d; %d subruns; %d sets",wholerun.run,nsubruns,nsets) << endl;
    
    return setvec;

}

///////////////////////////////////////////////////////////////

float GetMedian(vector<float> scores)
{
  float median;
  size_t size = scores.size();

  if(size==0) return 0;

  sort(scores.begin(), scores.end());

  if (size  % 2 == 0)
  {
      median = (scores[size / 2 - 1] + scores[size / 2]) / 2;
  }
  else 
  {
      median = scores[size / 2];
  }

  return median;
}

double GetMedian(TH1D* hist){

  int nbins = hist->GetNbinsX();

  double *x = new double[nbins];
  hist->GetXaxis()->GetCenter(x);

  double *y = hist->GetArray();

  double out = TMath::Median(nbins,x,y);

  delete x;

  return out;

}

///////////////////////////////////////////////////////////////

bool GetDBPartition(int run, int &par){

  par = 4;

  std::ifstream t;
  t.open(gSystem->Getenv("NOVADBPWDFILE"));      // open input file
  std::string buffer;
  getline(t,buffer);
  //cout<<buffer<<endl;
  char * dbpass;
  dbpass = new char [8];
  strcpy(dbpass, buffer.c_str());
  //cout<<dbpass<<endl;

  std::stringstream dbstream;
  std::string dbstring;
  dbstream << "select run, partition from fardet.runs where run = "
           << run
           << " order by run asc";
  dbstring = dbstream.str();
  // cout<<dbstring<<endl;
  char * dbcstr;
  dbcstr = new char [dbstring.size()+1];
  strcpy(dbcstr, dbstring.c_str());
  TSQLServer *dbserv = TSQLServer::Connect("pgsql://ifdbprod.fnal.gov:5433/nova_prod", "nova_reader", dbpass);

  int attempts = 1;

  while ( (dbserv==0 || !dbserv->IsConnected()) && attempts<100 ) {

    cout << "Try to connect again..." << endl;

    dbserv = TSQLServer::Connect("pgsql://ifdbprod.fnal.gov:5433/nova_prod", "nova_reader", dbpass);

    attempts++;

  }

  bool connected = false;

  if ( dbserv!=0 && dbserv->IsConnected() ) {

    connected = true;

    TSQLStatement *stmt = dbserv->Statement(dbcstr,1000);

    if(stmt->Process()) {
      // store result of statement in buffer
      stmt->StoreResult();
      while (stmt->NextResultRow()) {
        // Extract rows one after another
        par = stmt->GetUInt(1);
      }
    }
    // std::cout << "\n" << run << " " << subrun << " " << par;

    delete stmt;
  }

  delete dbserv;
  delete dbcstr;
  delete dbpass;

  if(attempts>1 && connected) cout << "Connected after " << attempts << " attempts." << endl;
  if(!connected)              cout << "Failed after " << attempts << " attempts." << endl;

  return connected;

}