JGlow.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 "TMath.h"
#include "TRandom3.h"
 
#include "JTools/JConstants.hh"
#include "JTools/JQuantile.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"
 
#include "JLightCurveBackgroundGenerator.hh"
 
using namespace JSUPERNOVA;
 
/**
 * \file
 *
 * Auxiliary program to test simulation of L0 background
 *
 * \author mlincett
 */
 
int main(int argc, char **argv) {
  using namespace std;
  using namespace JPP;
 
  string    inputFile;
  string    outputFile;
  int       debug;
  int       duration_ms;
  int       outNumberOfLines;
 
  try { 
 
    JParser<> zap("Auxiliary program to test generation of L0 background.");
    
    zap['f'] = make_field(inputFile,         "input file (output from JRipple).");
    zap['o'] = make_field(outputFile,        "output file.")                                       = "glow.root";
    zap['T'] = make_field(duration_ms,          "duration in ms of the background sample.")           = 500;
    zap['N'] = make_field(outNumberOfLines,  "number of lines of the simulated (output) detector") = 115;
    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);
  }
 
  NOTICE("Initialising generator." << endl);
 
  JLightCurveBackgroundGenerator simbad(in, true);
 
  in->Close();
 
  NOTICE("File loaded." << endl);
 
  int T_ms = duration_ms;
 
  simbad.configureTimeWindow(T_ms);
  simbad.setSeed();
 
  int inNumberOfLines = 1;
 
  simbad.configureRatio(inNumberOfLines, outNumberOfLines);
 
  JManager<int, TH1D> RT_SEQ(new TH1D("RT_SEQ_%", "Hit count"             , T_ms, 0, T_ms));
  JManager<int, TH1D> NC_SEQ(new TH1D("NC_SEQ_%", "Active channels"       , T_ms, 0, T_ms));  
  
  JManager<int, TH1D> RT_SHF(new TH1D("RT_SHF_%", "Hit count"             , T_ms, 0, T_ms));
  JManager<int, TH1D> NC_SHF(new TH1D("NC_SHF_%", "Active channels"       , T_ms, 0, T_ms));  
 
  JManager<string, TH1D> RT_FIT(new TH1D("RT_FIT_%", "Hit count"             , T_ms, 0, T_ms));
 
  NOTICE("Generation SEQ" << endl);
 
  for (int i = 0; i < 10; i++) {
    bg_type seq = simbad.generate();
 
    for (int j = 0; j < T_ms; j++) {
      RT_SEQ[i]->Fill(j, seq[0][j]);
      NC_SEQ[i]->Fill(j, seq[1][j]);
 
    }
 
  }
 
  NOTICE("Generation SHF" << endl);
 
  for (int i = 0; i < 10; i++) {
    bg_type shf = simbad.generate_shuffled();
      
    for (int j = 0; j < T_ms; j++) {
 
      RT_SHF[i]->Fill(j, shf[0][j]);
      NC_SHF[i]->Fill(j, shf[1][j]);
 
    }
  }
    
  NOTICE("Generation H2D + FIT" << endl);
 
  h2d_t* test1 = simbad.generate_H2D();
 
  bg_type fit1 = simbad.generate_fitted();
 
  for (int j = 0; j < T_ms; j++) {
    RT_FIT["SG_1"]->Fill(j, fit1[0][j]);
  }
    
  NOTICE("Generation H2D + FIT - rebin = 5" << endl);
 
  h2d_t* test5 = simbad.generate_H2D(5);
 
 
  bg_type fit5 = simbad.generate_fitted(5);
 
  for (size_t j = 0; j < fit5[0].size(); j++) {
    cout << fit5[0][j] << endl;
    RT_FIT["SG_5"]->Fill(j, fit5[0][j]);
  }
 
 
  /*
    NOTICE("Generation FIT" << endl);
 
    for (int i = 0; i < 10; i++) {
      
    // RT_FIT["BG"]->Fill(j, fit[1][j]);
    // RT_FIT["ER"]->Fill(j, fit[2][j]);
    }
 
    }*/
 
 
  TFile out(outputFile.c_str(), "recreate");
 
  TDirectory* dseq = out.mkdir("SEQ");
  TDirectory* dshf = out.mkdir("SHF");
  TDirectory* dfit = out.mkdir("FIT");
 
  RT_SEQ.Write(*dseq);
  NC_SEQ.Write(*dseq);
 
  RT_SHF.Write(*dshf);
  NC_SHF.Write(*dshf);
 
  RT_FIT.Write(*dfit);
 
  dfit->cd();
  test1->Write();
  test5->Write();
 
  out.Close();
  
}