JProcessKexing2D.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 "JPhysics/JConstants.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"
 
#include <string>
 
#include "JSupernova/JKexing2D.hh"
 
#include "JROOT/JRootToolkit.hh"
#include "JMath/JMathToolkit.hh"
 
/**
 * \file
 *
 * Auxiliary program to build supernova background from JKexing2D output
 * \author mlincett
 */
int main(int argc, char **argv) {
  using namespace std;
  using namespace JPP;
 
  string      inputFile;
  string      outputFile;
  string      summaryFile;
  int         windowSize;
  double      fractionThreshold;
  int         debug;
  JRange<int> multiplicity;
 
  try { 
 
    JParser<> zap("Auxiliary program to build supernova background from JKexing2D output");
    
    zap['f'] = make_field(inputFile,   "input file (JKexing2D).");
    zap['o'] = make_field(outputFile,  "output file.") = "sn_background.root";
    zap['w'] = make_field(windowSize,  "size of the sliding window to test") = 5;
    zap['M'] = make_field(multiplicity, "final multiplicity range") = JRange<int>(7,11);
    zap['F'] = make_field(fractionThreshold, "minimum fraction of active channels to compute distribution") = 0.99;
    zap['d'] = make_field(debug) = 1; 
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
  typedef JManager<int,    TH1D>  JManager_i1D_t;
  typedef JManager<string, TH1D>  JManager_s1D_t;
  typedef JManager<string, TH2D>  JManager2D_t;
 
  // input
  
  TFile in(inputFile.c_str(), "exist");
 
  TParameter<int>* runNumber;
 
  in.GetObject("RUNNR", runNumber);
 
  JManager2D_t   MT =   JManager2D_t::Read(in, mul_p   , '%');
  JManager_s1D_t ST = JManager_s1D_t::Read(in, status_p, '%');
 
  in.Close();
 
  // output
  
  const int    factoryLimit_peak = 250;
  const int    factoryLimit_runs = 10000;
 
  // distributions of individual multiplicities
 
  vector<string> filters = { "RAW", "TAV", "TA1" };
 
  map<string, JManager_i1D_t> mmap;
  
  for (vector<string>::const_iterator f = filters.begin(); f != filters.end(); ++f) {
    string title = "MD_" + (*f) + "_%";
    mmap[*f] = JManager_i1D_t(new TH1D(title.c_str(), NULL, factoryLimit_peak, 0, factoryLimit_peak)); 
  }
 
  const int    nx   = 1 + NUMBER_OF_PMTS; 
  const double xmin = -0.5;
  const double xmax = nx - 0.5;
 
  // distributions of the trigger level
 
  JManager_s1D_t TD(new TH1D(Form("SNT_[%d,%d]_", multiplicity.first, multiplicity.second) + TString("%"), NULL, factoryLimit_peak, -0.5, -0.5 + factoryLimit_peak));
 
  double epsilon = 1e-4;
 
  // bioluminescence benchmarking
  JManager2D_t  BL(new TH2D("BL_%", NULL, 100, epsilon, 1 + epsilon, factoryLimit_peak, -0.5, -0.5 + factoryLimit_peak));
  JManager2D_t  MUL_EFF(new TH2D("MUL_EFF_%", NULL, 100, epsilon, 1 + epsilon, nx, xmin, xmax));
 
  // history (value vs run number)
  JManager_s1D_t H(new TH1D("H_%", NULL, factoryLimit_runs, 0, factoryLimit_runs));
 
  // livetime vs fraction of active channels
 
  TH1D* LT = new TH1D("LIVETIME_ACF", NULL, 100, epsilon, 1 + epsilon);
 
  // fraction of active PMT
 
  TGraph* activeChannelFraction = histogramToGraph(*ST["PMT"]);
 
  const int       nb        = activeChannelFraction->GetN();
  const Double_t* arr_acf_y = activeChannelFraction->GetY();
 
  const vector<Double_t> vec_acf_y(arr_acf_y, arr_acf_y + nb);
 
  LT->FillN(nb, arr_acf_y, NULL);
 
  for (int M = 0; M <= NUMBER_OF_PMTS; M++) {
 
    map<string, TGraph*> count_vs_frame;
 
    // individual multiplicities
 
    for (vector<string>::const_iterator f = filters.begin(); f != filters.end(); ++f) {
 
      // number of events with multiplicity M as a function of the frame index
 
      const TH1* px = projectHistogram(*MT[*f], M, M, 'x');
 
      count_vs_frame[*f] = histogramToGraph(*px);
 
      delete px;
 
      // distribution of the number of events with multiplicity M
      // selection criterion based on the minimum fraction of active PMTs
 
      vector<double> select;
      vector<double> threshold(nb, fractionThreshold);
 
      // stores in 'select' the result of the element-wise comparison between 'vec_acf_y' and 'threshold'
 
      transform(vec_acf_y.begin(),
                vec_acf_y.end(),
                threshold.begin(),
                back_inserter(select),
                greater_equal<double> {});
      
      mmap[*f][M]->FillN(nb, count_vs_frame[*f]->GetY(), &select[0]);
 
      // multiplicity spectrum as a function of number of active PMTs
 
      const vector<double> nb_M(nb, M);
 
      MUL_EFF[*f]->FillN(nb,
                         arr_acf_y, 
                         &nb_M[0],
                         count_vs_frame[*f]->GetY());
 
    }
  }
 
  int RUNNR = runNumber->GetVal();
 
  for (vector<string>::const_iterator f = filters.begin(); f != filters.end(); ++f) {
 
    // distribution of number of events per frame in the selected multiplicity range (A)
 
    const TH1* px = projectHistogram(*MT[*f], multiplicity.getLowerLimit(), multiplicity.getUpperLimit(), 'x');
 
    TGraph* events_per_frame = histogramToGraph(*px);
 
    delete px;
 
    TD[*f]->FillN(events_per_frame->GetN(), events_per_frame->GetY(), NULL);
 
    // distribution of number of events per frame in the selected multiplicity range (A)
    // sliding window applied
 
    vector<double> vec(events_per_frame->GetY(), events_per_frame->GetY() + nb);
    vector<double> ker(windowSize, 1.0);
 
    vector<double> events_per_frame_sliding = convolve(vec, ker);
 
    TD[(*f) + "_nTS"]->FillN(events_per_frame_sliding.size(),
                             &events_per_frame_sliding[0],
                             NULL);
    
    // distribution of number of events per frame in the selected multiplicity range (A)
    // as a function of the fraction of active channels
    BL[*f]->FillN(nb,
                  arr_acf_y,
                  events_per_frame->GetY(),
                  NULL, 1);
 
    // peak as a function of the run number
 
    px = projectHistogram(*MT[*f], multiplicity.getLowerLimit(), multiplicity.getUpperLimit(), 'x');
 
    int peak = px->GetMaximum();
 
    delete px;
 
    H["PK_" + (*f)]->Fill(RUNNR, peak);
 
  }
 
  // history of bioluminescence
 
  double BI = (ST["PMT"]->GetSumOfWeights() / ST["PMT"]->GetEntries());
 
  H["BIOLUM"]->Fill(RUNNR, BI);
 
  // write output file
 
  TFile out(outputFile.c_str(), "recreate");
 
  TD.Write(out);
  
  TDirectory* bm = out.mkdir("BIOLUM");
  TDirectory* hs = out.mkdir("HISTORY");  
 
  H.Write(*hs);
  BL.Write(*bm);
  
  MUL_EFF.Write(*bm);
 
  for (vector<string>::const_iterator f = filters.begin(); f != filters.end(); ++f) {
    string dir_name = "MUL_" + (*f); 
    TDirectory* dir = out.mkdir(dir_name.c_str());
    mmap[*f].Write(*dir);
  }
 
  LT->Write();
 
  out.Close();
 
}