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();
}