#include #include #include #include #include #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 RT_SEQ(new TH1D("RT_SEQ_%", "Hit count" , T_ms, 0, T_ms)); JManager NC_SEQ(new TH1D("NC_SEQ_%", "Active channels" , T_ms, 0, T_ms)); JManager RT_SHF(new TH1D("RT_SHF_%", "Hit count" , T_ms, 0, T_ms)); JManager NC_SHF(new TH1D("NC_SHF_%", "Active channels" , T_ms, 0, T_ms)); JManager 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(); }