JMergeCalibrateK40.cc

Go to the documentation of this file.

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

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TKey.h"
#include "TString.h"

#include "km3net-dataformat/online/JDAQHeader.hh"

#include "JROOT/JRootFileReader.hh"
#include "JROOT/JRootFileWriter.hh"
#include "JSupport/JMeta.hh"

#include "JCalibrate/JCalibrateK40.hh"

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

namespace {

  double LIVETIME_S       =  10.0;        //!< minimal live time [s]
  double ERROR            =   1.0;        //!< default error [Hz/ns]
}


/**
 * \file
 * Auxiliary program to merge K40 data.
 *
 * \author mkarel, mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;

  vector<string>     inputFile;
  string             outputFile;
  int                debug;

  try { 

    JProperties properties;

    properties.insert(gmake_property(LIVETIME_S));
    properties.insert(gmake_property(ERROR));

    JParser<> zap("Auxiliary program to merge K40 data.");
    
    zap['@'] = make_field(properties)        = JPARSER::initialised();
    zap['f'] = make_field(inputFile,         "input file (output from JCalibrateK40).");
    zap['o'] = make_field(outputFile,        "output file (input to JFitK40).")   = "mergek40.root";
    zap['d'] = make_field(debug,             "debug.")                            = 1;

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


  gErrorIgnoreLevel = kError;

  TH1* weights_hist = NULL;

  typedef map<TString, TH1*>  map_type;

  map_type zmap;


  for (vector<string>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {

    DEBUG("Processing " << *i << endl) ;

    TFile in(i->c_str(), "read");
    
    // check if file contains weights vector

    if (!in.GetListOfKeys()->Contains(weights_hist_t)) {
      FATAL("Histogram does not contain weights histogram." << endl);
    }

    TH1* weights_hist_i = (TH1*) in.Get(weights_hist_t);

    if (weights_hist == NULL) {

      weights_hist = (TH1*) weights_hist_i->Clone();

    } else {

      // check normalization

      TAxis* x_axis = weights_hist->GetXaxis();

      if (weights_hist_i->GetBinContent(x_axis->FindBin(WS_t)) != 
          weights_hist  ->GetBinContent(x_axis->FindBin(WS_t))) {
        FATAL("Number of y-bins does not match for this file " << *i << endl);
      }

      if (weights_hist_i->GetBinContent(x_axis->FindBin(WB_t)) != 
          weights_hist  ->GetBinContent(x_axis->FindBin(WB_t))) {
        FATAL("TMax_ns is not the same for this file " << *i << endl);
      }

      // add normalization

      weights_hist->SetBinContent(x_axis->FindBin(W1_overall_t), 
                                  weights_hist  ->GetBinContent(x_axis->FindBin(W1_overall_t)) + 
                                  weights_hist_i->GetBinContent(x_axis->FindBin(W1_overall_t)));
    }


    TIter iter(in.GetListOfKeys());

    for (TKey* key; (key = (TKey*) iter.Next()) != NULL; ) {

      if (TString(key->GetName()).EndsWith(_2S) ||
          TString(key->GetName()).EndsWith(_1B) ||
          TString(key->GetName()).EndsWith(_1L)) {

        TH1* h1 = dynamic_cast<TH1*>(key->ReadObj());

        map_type::iterator p = zmap.find(h1->GetName());

        if (p == zmap.end()) {

          DEBUG("Clone " << h1->GetName() << endl);

          p = zmap.insert(make_pair(h1->GetName(), (TH1*) h1->Clone())).first;

        } else {

          DEBUG("Add   " << h1->GetName() << endl);

          p->second->Add(h1);
        }
      }
    }

    weights_hist->SetDirectory(0);
    
    for (map_type::iterator i = zmap.begin(); i != zmap.end(); ++i) {
      i->second->SetDirectory(0);
    }
    
    in.Close();
  }


  if (weights_hist == NULL) {
    FATAL("Missing histogram " << weights_hist_t << endl);
  }


  // write file

  TFile out(outputFile.c_str(), "recreate");

  const Double_t WS = weights_hist->GetBinContent(weights_hist->GetXaxis()->FindBin(WS_t)) ;  // [ns]
  const Double_t WB = weights_hist->GetBinContent(weights_hist->GetXaxis()->FindBin(WB_t)) ;  // [ns]

  for (map_type::iterator i = zmap.begin(); i != zmap.end(); ++i) {

    if (i->first.EndsWith(_2S)) {

      TH2D* h2s = (TH2D*) i->second;

      if (h2s->GetEntries() != 0) {

        TH1D* h1b = (TH1D*) zmap.find(TString(i->first).ReplaceAll(_2S, _1B))->second;
        TH1D* h1L = (TH1D*) zmap.find(TString(i->first).ReplaceAll(_2S, _1L))->second;

        //----------------------------------------------------------
        // normalisation and background estimation
        //----------------------------------------------------------

        for (int ix = 1; ix <= h2s->GetXaxis()->GetNbins(); ++ix) {

          // normalised background estimation
          // if livetime is compatible with zero then set bin contents
          // and erros to zero and default value, respectively

          if (h1L->GetBinContent(ix) <= LIVETIME_S) {

            h1b->SetBinContent(ix, 0.0);
            h1b->SetBinError  (ix, 0.0);

            for (int iy = 1; iy <= h2s->GetYaxis()->GetNbins(); ++iy) {
              h2s->SetBinContent(ix, iy, 0.0);
              h2s->SetBinError  (ix, iy, ERROR);
            }

          } else {

            const double Y = h1b->GetBinContent(ix) / h1L->GetBinContent(ix) / WB;
            const double W = h1b->GetBinError  (ix) / h1L->GetBinContent(ix) / WB;

            h1b->SetBinContent(ix, Y);
            h1b->SetBinError  (ix, W);

            for (int iy = 1; iy <= h2s->GetYaxis()->GetNbins(); ++iy) {

              double y = h2s->GetBinContent(ix, iy) / h1L->GetBinContent(ix) / WS;
              double w = h2s->GetBinError  (ix, iy) / h1L->GetBinContent(ix) / WS;

              if (w == 0.0) {

                // for chi2 minimisation; accordingly set error of empty bins

                w = 1.0 / h1L->GetBinContent(ix) / WS;
              }

              h2s->SetBinContent(ix, iy, y - Y);
              h2s->SetBinError  (ix, iy, sqrt(w*w + W*W));
            }
          } 
        }

        // reset under and overflows.

        for (int ix = 0; ix <= h2s->GetXaxis()->GetNbins() + 1; ++ix) {
          h2s->SetBinContent(ix, 0, 0.0);  h2s->SetBinContent(ix, h2s->GetYaxis()->GetNbins() + 1, 0.0);
          h2s->SetBinError  (ix, 0, 0.0);  h2s->SetBinError  (ix, h2s->GetYaxis()->GetNbins() + 1, 0.0);
        }

        for (int iy = 0; iy <= h2s->GetYaxis()->GetNbins() + 1; ++iy) {
          h2s->SetBinContent(0, iy, 0.0);  h2s->SetBinContent(h2s->GetXaxis()->GetNbins() + 1, iy, 0.0);
          h2s->SetBinError  (0, iy, 0.0);  h2s->SetBinError  (h2s->GetXaxis()->GetNbins() + 1, iy, 0.0);
        }

        NOTICE("Biased coincidence rate [Hz] " << h2s->GetName() << ' ' << h2s->GetSumOfWeights() * WS << endl);

        h2s->SetName(TString(i->first).ReplaceAll(_2S, _2R));

        h2s->Write();
      }
    } 
  }

  putObject(out, JMeta(argc, argv));

  for (vector<string>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {

    JMeta::copy(i->c_str(), out);

    for (JRootFileReader<JDAQHeader> in(i->c_str()); in.hasNext(); ) {
      putObject(out, *in.next());
    }
  }

  out.Write();
  out.Close();
}