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 "JDAQ/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;

  try { 

    JParser<> zap("Auxiliary program to fit multiplicity rates from JMonitorMultiplicity output.");
    
    zap['f'] = make_field(inputFile,         "input file (output from JMonitorMultiplicity).");
    zap['o'] = make_field(outputFile,        "output file.")                                           = "fitmultiplicity.root";
    zap['O'] = make_field(option,            "fitting option")                                         = "Q";
    zap['B'] = make_field(fitBackground,     "fit and subtract background instead of fitting signal");
    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;

  JManager_t HTF(new TH1D("%F", "",  nx, xmin, xmax));
  
  HTF->Sumw2(kFALSE);
  
  JManager_t HTFE(new TH1D("%FE", "", nx, xmin, xmax));

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

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

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

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

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

  //       bin = 1 -> nominal livetime
  for (int bin = 2; bin <= hLiveTime->GetNbinsX(); bin++) { 
 
    TString hLabel = hLiveTime->GetXaxis()->GetBinLabel(bin) + TString("_HT");

    DEBUG("Reading livetime for " << hLabel << endl);

    double T = hLiveTime->GetBinContent(bin);

    if (T > 0) {
     
      TDirectory *cdir = out.mkdir(hLabel);
      
      cdir->cd();
     
      TH2D* hDOM = (TH2D*) in->Get(hLabel);

      if (hDOM == NULL) {
        // inconsistent file, should be FATAL?
        ERROR("Histogram " << hLabel << " not found, skipping DOM." << endl);
        continue;
      }

      //----------------------------------------------------------
      // loop over multiplicities
      //----------------------------------------------------------
        
      TH1D *hI, *hIF;

      for (int M = NUMBER_OF_PMTS; M >= 2; M--) {

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

        TH1D* hM = hDOM->ProjectionY(hLabel + Form("_%d", M), Mbin, Mbin);

        if (M == NUMBER_OF_PMTS) {
          hI = (TH1D*) hM->Clone(hLabel + "I");
        } else {
          hI->Add(hM);
        }

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

          double inclusiveCountVal = 0, inclusiveCountErr = 0;

          hIF = (TH1D*) hI->Clone(hLabel + "IF" + Form("_%d", M)); hIF->Sumw2();

          if (hIF->GetRMS() > 0.1) {

            hIF->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, hIF->GetMaximum());
            f1.SetParameter(1, hIF->GetMean());
            f1.SetParameter(2, hIF->GetRMS() * 0.25);
            f1.SetParameter(3, hIF->GetMinimum());
              
            hIF->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 = hIF->GetSumOfWeights();

            if (fitBackground) {

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

            } else {

              inclusiveCountErr = sg_e;
              inclusiveCountVal = sg_v;

            }         

            hIF->Scale(1.0 / T);
            hIF->Fit(&f1, option.c_str(), "same");
            
            
          }
                                       
          if (inclusiveCountVal <= 0) { inclusiveCountVal = 0; }

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

          hIF->Write();

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

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

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

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

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

          double exclusiveCountVal = 0;

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

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

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

  HTF.Write(out);

  HTFE.Write(out);

  out.Close();

}