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 "JROOT/JMinimizer.hh"
 
#include "JPhysics/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 "JROOT/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();
 
}