JFitSinglesRates.cc

Go to the documentation of this file.

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

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TF1.h"
#include "TF2.h"
#include "TMath.h"

#include "JPhysics/JConstants.hh"
#include "km3net-dataformat/online/JDAQ.hh"
#include "JDAQ/JDAQSummarysliceIO.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"

#include "JTools/JRange.hh"

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

/**
 * \file
 *
 * Auxiliary program to fit singles rate distributions. The fitted parameters are the mean singles rate and the spread.
 * \author mkarel
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  string             inputFile;
  string             outputFile;
  string             detectorFile;
  bool               setDefaultLimits;
  double             peakFraction;
  string             option;
  bool               writeFits;
  int                debug;

  //Parser module handles input parameters
  try { 

    JParser<> zap("Auxiliary program to fit singles rate distributions.");
    
    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)          = "fit.root";
    zap['a'] = make_field(detectorFile);
    zap['L'] = make_field(setDefaultLimits);
    zap['p'] = make_field(peakFraction)        = 0.1;
    zap['O'] = make_field(option)              = "R W M";
    zap['w'] = make_field(writeFits);
    zap['d'] = make_field(debug)               = 1;

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


  using namespace KM3NETDAQ;

  JDetector detector;

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

  TFile* in = TFile::Open(inputFile.c_str(), "exist");

  if (in == NULL || !in->IsOpen()) {
    FATAL("File: " << inputFile << " not opened." << endl);
  }

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

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

    DEBUG("Module " << module->getID() << endl);
    const JString title = JString(module->getID());

    TH2D* h2s = (TH2D*) in->Get((title+".2S").c_str());

    if (h2s == NULL) {

      WARNING("No data in histogram " << title << endl);

      continue;
    }

    const double factor = 1.0/1000 ;

    TH1D mean((title+".1mean").c_str(), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS-0.5) ;
    TH1D sigma((title+".1sigma").c_str(), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS-0.5) ;

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

      TH1D slice((title+Form(".%i.1S", i-1)).c_str(), NULL, h2s->GetYaxis()->GetNbins(), JDAQRate::getData(factor)) ;

      slice.Sumw2() ;

      for (int j = 1 ; j <= h2s->GetNbinsY() ; ++j) {

        slice.SetBinContent(j, h2s->GetBinContent(i, j)) ;
        slice.SetBinError  (j, sqrt(h2s->GetBinContent(i, j))) ;

      }

      if (slice.Integral() <= 0.0) {

        WARNING("No data in PMT " << i-1 << " of module " << title << ", skipping fit" << endl) ;

        continue ;
      }

      slice.Scale(1.0/slice.Integral()) ;

      // determine fit range
      int    binmax = slice.GetMaximumBin() ;
      double ratemax = slice.GetBinContent(binmax) ;

      int bin_l = binmax ;
      int bin_r = binmax ;
      while ((slice.GetBinContent(bin_l)==0 || slice.GetBinContent(bin_l) >= peakFraction*ratemax) && bin_l > 0)                 { bin_l-- ; } ; bin_l++ ;
      while ((slice.GetBinContent(bin_r)==0 || slice.GetBinContent(bin_r) >= peakFraction*ratemax) && bin_r < slice.GetNbinsX()) { bin_r++ ; } ; bin_r-- ;

      JRange<double> fitrange(slice.GetXaxis()->GetBinCenter(bin_l), slice.GetXaxis()->GetBinCenter(bin_r)) ;

      TF1 f1("f1", "[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2]))/(sqrt(2*pi)*[2])");

      // default parameters for start of fit
      f1.SetParameter(0, 0.2) ;
      f1.SetParameter(1, 6.0) ; // [kHz]
      f1.SetParameter(2, 0.3) ; // [kHz]

      if (setDefaultLimits) {
        f1.SetParLimits(0,     0.0,  10.0);
        f1.SetParLimits(1,     0.0,  20.0); // [kHz]
        f1.SetParLimits(2,     0.0,  1.0);  // [kHz]
      }

      if (debug < JEEP::debug_t) {
        option += " Q";
      }

      DEBUG("Fit histogram " << slice.GetName() << " in range " << fitrange << " kHz " << endl ) ;

      slice.Fit(&f1, option.c_str(), "same", max(0.0, fitrange.first), fitrange.second) ;

      DEBUG( f1.GetParameter(0)<<" "<<f1.GetParameter(1)<<" "<<f1.GetParameter(2)<<endl ) ;

      mean.SetBinContent(i, f1.GetParameter(1)) ;
      sigma.SetBinContent(i, f1.GetParameter(2)) ;

      if (writeFits) {
        slice.Write() ;
      }

    }

    mean.Write() ;
    sigma.Write() ;

  }

  out.Close() ;

}