JFitMultiplicity.cc

Go to the documentation of this file.

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

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

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

#include "JROOT/JRootToolkit.hh"
#include "JTools/JRange.hh"
#include "JSupport/JMeta.hh"

#include "JCalibrate/JCalibrateK40.hh"
#include "JCalibrate/JFitK40.hh"

#include "JGizmo/JManager.hh"

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

const char* livetime_t = "LiveTime";

/**
 * \file
 *
 * Auxiliary program to fit multiplicity rates from JMonitorMultiplicity output.
 * \author mlincett
 */

int main(int argc, char **argv) {
  using namespace std;
  using namespace JPP;

  string    inputFile;
  string    outputFile;
  string    summaryFile;
  string    option;
  int       debug;
  bool      fitBackground;
  bool      nominalLiveTime;

  try { 
    JParser<> zap("Auxiliary program to fit multiplicity rates from JMonitorMultiplicity output.");
    
    zap['f'] = make_field(inputFile,         "input file");
    zap['o'] = make_field(outputFile,        "output file")                                           = "fitmultiplicity.root";
    zap['O'] = make_field(option,            "fitting option")                                        = "Q";
    zap['B'] = make_field(fitBackground,     "get signal by fit and subtraction from background from the integral");
    zap['L'] = make_field(nominalLiveTime,   "use nominal live time from LiveTime histogram, instead of individual DOMs");
    zap['d'] = make_field(debug) = 1;   

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

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

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

  
  //----------------------------------------------------------
  // general histograms
  //----------------------------------------------------------

  const int    nx   = 1 + NUMBER_OF_PMTS; 
  const double xmin = -0.5;
  const double xmax = nx - 0.5;
  
  typedef JManager<TString, TH1D> JManager_t;

  // inclusive
  JManager_t HTF(new TH1D("%F", "",  nx, xmin, xmax));
  // exclusive
  JManager_t HTFE(new TH1D("%FE", "", nx, xmin, xmax));

  HTF->Sumw2(kFALSE);

  //----------------------------------------------------------
  // get overall normalization constants
  //----------------------------------------------------------

  TH1D* hLiveTime = (TH1D*) in->Get(livetime_t);

  if (hLiveTime == NULL) {
    FATAL("File " << inputFile << " does not contain livetime histogram." << endl);
  }
  
  //----------------------------------------------------------
  // initialise output file
  //----------------------------------------------------------

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

  //----------------------------------------------------------
  // loop over bins of livetime histogram
  //----------------------------------------------------------

  TDirectory *fdir = out.mkdir("FIT");
  TDirectory *pdir = out.mkdir("MUL"); 

  // nominal live time
  double nT = hLiveTime->GetBinContent(1);

  for (int bin = 2; bin <= hLiveTime->GetNbinsX(); bin++) { 
 
    TString prefix = hLiveTime->GetXaxis()->GetBinLabel(bin);

    DEBUG("Processing data from " << prefix << endl);

    double T = nominalLiveTime ? nT : hLiveTime->GetBinContent(bin);

    TString h1D_label = prefix + TString("_P" );
    TString h2D_label = prefix + TString("_HT");

    if (T > 0) {
     
      TH2D* h1D = (TH2D*) in->Get(h1D_label);
      TH2D* h2D = (TH2D*) in->Get(h2D_label);

      if (h1D != NULL) {

        h1D->Scale(1.0 / T);

        pdir->cd();

        h1D->Write();
        
      } else {
        DEBUG(h1D_label << " not found" << endl);
      }

      if (h2D != NULL) {

        TDirectory* cdir = fdir->mkdir(h2D_label); cdir->cd();     

        for (int M = 2; M <= NUMBER_OF_PMTS; M++) {

          int Mbin = h2D->GetXaxis()->FindBin(M);

          TH1D* hI = h2D->ProjectionY(h2D_label + Form("_%d", M), Mbin, NUMBER_OF_PMTS);

          hI->Sumw2();

          if (hI->GetEntries() > 0) {
            
            const double W = (hI->GetXaxis()->GetXmax() - hI->GetXaxis()->GetXmin()); 
            const double N = hI->GetNbinsX();
            const double V = W / N;
            const double minRMS = 0.1;

            double inclusiveCountVal = 0, inclusiveCountErr = 0;
            
            if (hI->GetRMS() > minRMS) {

              hI->Scale(1.0 / V);

              TF1 f1("f1", "[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2]))/(TMath::Sqrt(2*TMath::Pi())*[2]) + [3]");

              f1.SetParameter(0, hI->GetMaximum());
              f1.SetParameter(1, hI->GetMean());
              f1.SetParameter(2, hI->GetRMS() * 0.25);
              f1.SetParameter(3, hI->GetMinimum());
              
              hI->Fit(&f1, option.c_str(), "same");

              double sg_v = f1.GetParameter(0);
              double sg_e = f1.GetParError( 0);
              
              double bg_v = f1.GetParameter(3);
              double bg_e = f1.GetParError( 3);

              // 2-sigma acceptance criteria for fitted values
              if (bg_v / bg_e <= 2.0) { bg_v = bg_e = 0; }
              if (sg_v / sg_e <= 2.0) { sg_v = sg_e = 0; }

              double integral = hI->GetSumOfWeights();

              if (fitBackground) {

                inclusiveCountVal = V * (integral - (bg_v * N));
                inclusiveCountErr = V * bg_e * N;

              } else {

                inclusiveCountErr = sg_e;
                inclusiveCountVal = sg_v;

              }               

              hI->Scale(1.0 / T);

              // re-fit after scaling
              hI->Fit(&f1, option.c_str(), "same");         
            
            }
                                       
            if (inclusiveCountVal <= 0) { inclusiveCountVal = 0; }

            // --------------------------------
            // write fitted histogram
            // --------------------------------

            hI->Write();

            // --------------------------------
            // build inclusive rates histograms
            // --------------------------------

            int bin = HTF[h2D_label]->FindBin(M);

            HTF[h2D_label]->SetBinContent(bin, inclusiveCountVal);

            double binError = sqrt(inclusiveCountVal + pow(inclusiveCountErr, 2));

            HTF[h2D_label]->SetBinError(bin, binError);
            
            // --------------------------------
            // build exclusive rates histograms
            // --------------------------------

            double exclusiveCountVal = 0;

            if (M == NUMBER_OF_PMTS) {
              exclusiveCountVal = inclusiveCountVal;
            } else {
              int prevBin = HTF[h2D_label]->GetXaxis()->FindBin(M+1);
              exclusiveCountVal = inclusiveCountVal - HTF[h2D_label]->GetBinContent(prevBin); 
            }
               
            HTFE[h2D_label]->SetBinContent(HTFE[h2D_label]->FindBin(M), exclusiveCountVal);
              
          
          }
          
        }

      } else {
        DEBUG(h1D_label << " not found" << endl);
      } 

      // HTF[hLabel]->Sumw2(kTRUE);
        
      HTFE[h2D_label]->Sumw2();

      HTF[h2D_label ]->Scale(1.0 / T);
      HTFE[h2D_label]->Scale(1.0 / T);
      
    }
  }

  fdir = out.mkdir("HTF" ); fdir->cd();

  HTF.Write(*fdir);

  fdir = out.mkdir("HTFE"); fdir->cd();

  HTFE.Write(*fdir);

  out.Close();

}