JCheckK40.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <limits>
#include <vector>
#include <map>

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TF1.h"
#include "TFitResult.h"

#include "km3net-dataformat/online/JDAQ.hh"
#include "JDAQ/JDAQHeaderIO.hh"
#include "JROOT/JRootFileReader.hh"

#include "JDB/JDB.hh"
#include "JDB/JSelector.hh"
#include "JDB/JSelectorSupportkit.hh"
#include "JDB/JDBToolkit.hh"
#include "JDB/JPMTHVRunSettings.hh"
#include "JDB/JLocation_t.hh"
#include "JDB/JSupport.hh"

#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"

#include "JCalibrate/JCalibrateK40.hh"
#include "JROOT/JManager.hh"
#include "JTools/JRange.hh"
#include "JTools/JQuantile.hh"

#include "Jeep/JProperties.hh"
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"


namespace {

  using namespace std;
  using namespace JPP;

  bool     usePMTID;                    // option for old data for which correction for PMT cable swaps does not apply

  /**
   * Setup corresponding to an output file from JCalibrateK40.
   */
  struct JSetup {
    
    /**
     * Constructor.
     *
     * \param  file_name        file name
     */
    JSetup(const string& file_name) :
      file_name(file_name)
    {
      fp = TFile::Open(file_name.c_str(), "read");
      
      int counter = 0;
      
      for (JRootFileReader<JDAQHeader> in(file_name.c_str()); in.hasNext(); ) {

        const JDAQHeader* p = in.next();

        if (counter == 0)
          header = *p;
        else
          THROW(JException, "Multiple headers in file " << file_name);
      }

      ResultSet& rs  = getResultSet(getTable   <JPMTHVRunSettings>(),
                                    getSelector<JPMTHVRunSettings>(getDetector(header.getDetectorID()), header.getRunNumber()));
    
      for (JPMTHVRunSettings parameters; rs >> parameters; ) {

        JLocation_t location;

        if (usePMTID)
          location = JLocation_t(parameters.DUID, parameters.FLOORID, parameters.PMTINTID);
        else
          location = JLocation_t(parameters.DUID, parameters.FLOORID, parameters.CABLEPOS);
        
        HV[location] = parameters.HV_VALUE;
      }
      
      rs.Close();
    }

    /**
     * Get histogram from file.
     *
     * \param  key              histogram name
     * \return                  pointer to histogram
     */
    TH2D* get(const string& key) const
    {
      return (TH2D*) fp->Get(key.c_str());
    }
    
    string                   file_name;
    JDAQHeader               header;
    map<JLocation_t, double> HV;

  private:
    TFile*                   fp;
  };
}


/**
 * \file
 *
 * Auxiliary program to check t0's.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;

  typedef JRange<size_t>  JRange_t;
  
  JServer               server;
  string                usr;
  string                pwd;
  string                cookie;
  string                detectorFile;
  pair<string, string>  inputFile;
  string                outputFile;
  double                precision;
  double                Wmin         = 100.0;
  JRange_t              X;
  string                option;
  int                   debug;

  try {

    JProperties properties;

    properties.insert(gmake_property(Wmin));
    
    JParser<> zap("Auxiliary program to check t0's.");
    
    zap['s'] = make_field(server)        = getServernames();
    zap['u'] = make_field(usr)           = "";
    zap['!'] = make_field(pwd)           = "";
    zap['C'] = make_field(cookie)        = "";
    zap['a'] = make_field(detectorFile);
    zap['f'] = make_field(inputFile,         "pair of input files (output of JCalibrateK40)");
    zap['o'] = make_field(outputFile)    = "k40.root";
    zap['e'] = make_field(precision,         "precision for HV comparison")             = 0.5;
    zap['@'] = make_field(properties)    = JPARSER::initialised();
    zap['x'] = make_field(X,                 "ROOT fit range (PMT pairs).")             = JRange_t(300, numeric_limits<size_t>::max());
    zap['O'] = make_field(option,            "ROOT fit option, see TH1::Fit.")          = "";
    zap['U'] = make_field(usePMTID);
    zap['d'] = make_field(debug)         = 2;

    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }

  
  try {
    JDB::reset(usr, pwd, cookie);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }

  JSetup setups[] = {
    JSetup(inputFile.first), 
    JSetup(inputFile.second)
  };

  for (int i = 0; i != 2; ++i) {
    DEBUG(setw(32) << setups[i].file_name << ' ' << setups[i].header.getDetectorID() << ' ' << setups[i].header.getRunNumber() << endl);
  }
    
  JDetector detector;

  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }

  if (option.find('S') == string::npos) { option += 'S'; }
  if (option.find('Q') == string::npos && debug < JEEP::debug_t) { option += 'Q'; }


  JManager<int, TH1D>  H1(new TH1D("[%]", NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5));


  TF1 f1("f1", "[0]*TMath::Gaus(x,[1],[2]) + [3]");

  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {

    TH2D* h2[] = {
      setups[0].get(MAKE_CSTRING(module->getID() << _2S)),
      setups[1].get(MAKE_CSTRING(module->getID() << _2S))
    };

    DEBUG("Module " << setw(10) << module->getID() << ' ' << (h2[0] != NULL) << (h2[0] != NULL) << endl);
    
    if (h2[0] == NULL ||
        h2[1] == NULL) {
      continue;
    }

    JCombinatorics combinatorics;

    combinatorics.configure(module->size());

    combinatorics.sort(JPairwiseComparator(*module));

    vector<JQuantile> Q(module->size());
    
    for (size_t ip = max(X.getLowerLimit(), size_t(0)); ip != min(X.getUpperLimit(), combinatorics.getNumberOfPairs()); ++ip) {

      const JCombinatorics::pair_type pair = combinatorics.getPair(ip);

      const JLocation_t location_1(module->getString(), module->getFloor(), pair.first);
      const JLocation_t location_2(module->getString(), module->getFloor(), pair.second);

      const bool hv_1 = (fabs(setups[0].HV[location_1] - setups[1].HV[location_1]) < precision);
      const bool hv_2 = (fabs(setups[0].HV[location_2] - setups[1].HV[location_2]) < precision);

      double t1[] = {
        numeric_limits<double>::max(),
        numeric_limits<double>::max()
      };
      
      const Int_t ix = ip + 1;

      for (int i = 0; i != 2; ++i) {
      
        TH1D h1("__py", NULL, h2[i]->GetYaxis()->GetNbins(), h2[i]->GetYaxis()->GetXmin(), h2[i]->GetYaxis()->GetXmax());

        // start values

        Double_t ymin   =  numeric_limits<double>::max();
        Double_t ymax   =  numeric_limits<double>::lowest();
        Double_t mean   =  0.0;
        Double_t sigma  =  4.5;
        Double_t W      =  0.0;

        for (int iy = 1; iy <= h1.GetNbinsX(); ++iy) {

          const Double_t x = h1.GetBinCenter(iy);
          const Double_t y = h2[i]->GetBinContent(ix,iy);
        
          h1.SetBinContent(iy, y);
          h1.SetBinError  (iy, sqrt(y));

          if (y > ymax) { 
            mean = x;
            ymax = y;
          }

          if (y < ymin) {
            ymin = y;
          }

          W += y;
        }

        if (W >= Wmin) {
        
          f1.SetParameter(0, ymax);
          f1.SetParameter(1, mean);
          f1.SetParameter(2, sigma);
          f1.SetParameter(3, ymin);
        
          TFitResultPtr result = h1.Fit(&f1, option.c_str(), "same");

          if (result.Get() == NULL) {
            FATAL("Invalid TFitResultPtr " << h1.GetName() << endl);
          }

          if (debug >= debug_t || !result->IsValid()) {
            cout << "Histogram slice: "
                 << setw(3)    << ix                    << ' '
                 << FIXED(7,3) << f1.GetParameter(1)    << " +/- "
                 << FIXED(7,3) << f1.GetParError(1)     << ' '
                 << FIXED(7,3) << result->Chi2()        << '/'
                 << result->Ndf()                       << ' '
                 << (result->IsValid() ? "" : "failed") << endl;
          }

          t1[i] = f1.GetParameter(1);
        }
      }

      if (t1[0] != numeric_limits<double>::max() &&
          t1[1] != numeric_limits<double>::max()) {

        if (hv_1 != hv_2) {

          JCombinatorics::pair_type p2;
          
          if (hv_1) {
            p2 = JCombinatorics::pair_type(pair.second, pair.first);
          }

          if (hv_2) {
            p2 = JCombinatorics::pair_type(pair.first, pair.second);
          }

          if (debug >= debug_t) {
            cout << setw(10) << module->getID() << "." << FILL(2,'0') << p2.first  << FILL() << ' ';
            cout << "(" << FILL(2,'0') << p2.second << FILL() << ")" << ' ';
            cout << FIXED(6,2) << (combinatorics.getSign(p2) * (t1[1] - t1[0])) << endl;
          }

          Q[p2.first].put(combinatorics.getSign(p2) * (t1[1] - t1[0]));
        }
      }
    }

    for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
      if (Q[i].getCount() >= 2) {
        H1[module->getID()]->SetBinContent(i+1, Q[i].getMean());
        H1[module->getID()]->SetBinError  (i+1, Q[i].getSTDev());
      }
    }
  }

  if (outputFile != "") {
    H1.Write(outputFile.c_str());
  }
}