JFitSinglesRates.cc File Reference

Auxiliary program to fit singles rate distributions. More...

#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 "JROOT/JMinimizer.hh"
#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"

Go to the source code of this file.

Functions
int	main (int argc, char **argv)

Detailed Description

Auxiliary program to fit singles rate distributions.

The fitted parameters are the mean singles rate and the spread.

Author

mkarel

Definition in file JFitSinglesRates.cc.

Function Documentation

main()

int main	(	int	argc,
		char **	argv )

Definition at line 33 of file JFitSinglesRates.cc.

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