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 "JGizmo/JManager.hh"

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

#include "JL0BackgroundSimulator.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);

  JL0BackgroundSimulator 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 = JL0BackgroundSimulator::fit_H2D(test5);

  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();
  
}