Jpp  debug
the software that should make you happy
Functions
JMergeRunAnalyzer.cc File Reference
#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <map>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TKey.h"
#include "TString.h"
#include "TMath.h"
#include "TTimeStamp.h"
#include "km3net-dataformat/online/JDAQHeader.hh"
#include "JROOT/JRootFileReader.hh"
#include "JROOT/JRootFileWriter.hh"
#include "JSupport/JFileRecorder.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JSingleFileScanner.hh"
#include "JSupport/JTreeScanner.hh"
#include "JSupport/JMeta.hh"
#include "JTrigger/JTriggerInterface.hh"
#include "JTrigger/JTriggerBits.hh"
#include "JGizmo/JGizmoToolkit.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
#include "JRunAnalyzer/JRunInfo.hh"

Go to the source code of this file.

Functions

void fill2DHistogram (int run, TH1D *h1, TH2D *h2, bool w_mc, double daq_livetime=1, double mc_livetime=1)
 
int main (int argc, char **argv)
 

Function Documentation

◆ fill2DHistogram()

void fill2DHistogram ( int  run,
TH1D *  h1,
TH2D *  h2,
bool  w_mc,
double  daq_livetime = 1,
double  mc_livetime = 1 
)
inline

Definition at line 45 of file JMergeRunAnalyzer.cc.

45  {
46 
47  if(h1 != NULL && h2 != NULL){
48 
49  double w = 1;
50 
51  for(int i = 1; i <= h1->GetNbinsX(); ++i){
52 
53  if(w_mc){
54  w = daq_livetime / mc_livetime;
55  }
56 
57  h2->Fill(run, h1->GetXaxis()->GetBinCenter(i), h1->GetBinContent(i) * w);
58  }
59  }
60 }
JTOOLS::w
data_type w[N+1][M+1]
Definition: JPolint.hh:867

◆ main()

int main ( int  argc,
char **  argv 
)

Definition at line 63 of file JMergeRunAnalyzer.cc.

64 {
65  using namespace std;
66  using namespace JPP;
67  using namespace KM3NETDAQ;
68 
69  vector<string> inputFile;
70  JFileRecorder<JTYPELIST<run_info, JMeta, TH1D, TH2D>::typelist> outputFile;
71  bool weight_mc;
72  int debug;
73 
74  try {
75 
76  JParser<> zap("Auxiliary program to merge JRunAnalyzer histograms.");
77 
78  zap['f'] = make_field(inputFile, "input file (output from JRunAnalyzer).");
79  zap['o'] = make_field(outputFile, "output file.") = "merge-jra.root";
80  zap['w'] = make_field(weight_mc, "weight mc events");
81  zap['d'] = make_field(debug, "debug.") = 1;
82 
83  zap(argc, argv);
84  }
85  catch(const exception &error) {
86  FATAL(error.what() << endl);
87  }
88 
89  int min_run = numeric_limits<int>::max(), max_run = 0;
90 
91  size_t min_date = numeric_limits<int>::max(), max_date = 0;
92 
93  run_info runinfo;
94 
95  double total_mc_livetime = 0;
96  int run = 0;
97 
98  for (vector<string>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {
99 
100  {
101  JSingleFileScanner<run_info> in(*i);
102  if(in.hasNext()){
103  runinfo = *in.next();
104  if(in.hasNext()){
105  FATAL("MULTIPLE RUN INFO");
106  }
107 
108  int run_id = runinfo.run_id;
109 
110  if(run_id != run){
111 
112  run = run_id;
113 
114  total_mc_livetime = runinfo.mc_livetime;
115 
116  } else {
117 
118  total_mc_livetime += runinfo.mc_livetime;
119  }
120 
121  } else {
122  FATAL("NO RUN INFO");
123  }
124  }
125 
126  for (JRootFileReader<JDAQHeader> in(i->c_str()); in.hasNext(); ) {
127 
128  const JDAQHeader* p = in.next();
129 
130  if(p->getTimesliceStart().getUTCseconds() < min_date) min_date = p->getTimesliceStart().getUTCseconds();
131  if(p->getTimesliceStart().getUTCseconds() > max_date) max_date = p->getTimesliceStart().getUTCseconds();
132  if(p->getRunNumber() < min_run) min_run = p->getRunNumber();
133  if(p->getRunNumber() > max_run) max_run = p->getRunNumber();
134  }
135  }
136 
137  time_t tmin = min_date, tmax = max_date;
138 
139  TTimeStamp date_min(tmin, 0);
140  TTimeStamp date_max(tmax, 0);
141 
142  cout << "Minimum RUN: " << min_run << ", "; date_min.Print("s");
143  cout << "Maximum RUN: " << max_run << ", "; date_max.Print("s");
144 
145  int xbins = (max_run - min_run) + 1;
146  string title = string("Data taken from: ") + date_min.AsString("s") + " - " + date_max.AsString("s");
147 
148  TH2D* h_pmt_rate_distribution = new TH2D("h_pmt_rate_distribution", string(title + "; RUN; PMT Rate Distribution [kHz]").c_str(), xbins, (double) min_run-0.5, (double) max_run+0.5, JDAQRate::getN(), JDAQRate::getData(1.0e-4));
149  TH2D* h_Trigger_bit_hit = new TH2D("h_Trigger_bit_hit", string(title + "; RUN; Hits trigger bit").c_str(), xbins, (double) min_run-0.5, (double) max_run+0.5, NUMBER_OF_TRIGGER_BITS, -0.5, NUMBER_OF_TRIGGER_BITS - 0.5);
150  TH2D* h_Triggered_hits = new TH2D("h_Triggered_hits", string(title + "; RUN; Number of triggered hits").c_str(), xbins, (double) min_run-0.5, (double) max_run+0.5, 50, 0, 4);
151  setLogarithmicY(h_Triggered_hits);
152  TH2D* h_Triggered_hits_3dmuon = new TH2D("h_Triggered_hits_3dmuon", string(title + "; RUN; Number of triggered hits - JTRIGGER3DMUON").c_str(), xbins, (double) min_run-0.5, (double) max_run+0.5, 50, 0, 3);
153  setLogarithmicY(h_Triggered_hits_3dmuon);
154  TH2D* h_Snapshot_hits = new TH2D("h_Snapshot_hits", string(title + "; RUN; Number of snapshot hits").c_str(), xbins, (double) min_run-0.5, (double) max_run+0.5, 50, 0, 4);
155  setLogarithmicY(h_Snapshot_hits);
156  TH2D* h_Triggered_over_Snapshot_hits = new TH2D("h_Triggered_over_Snapshot_hits", string(title + "; RUN; Triggered/Snapshot hits").c_str(), xbins, (double) min_run-0.5, (double) max_run+0.5, 100, 0, 0.5);
157  TH2D* h_Number_of_overlays = new TH2D("h_Number_of_overlays", string(title + ";RUN;Number of Overlays").c_str(), xbins, (double) min_run-0.5, (double) max_run+0.5, 1000, -0.5, 1000 - 0.5);
158  TH2D* h_pmt_distribution_triggered_hits = new TH2D("h_pmt_distribution_triggered_hits", string(title + ";RUN;Normalised PMT Distrib triggered hits (upper PMTs=[0-11])").c_str(), xbins, (double) min_run-0.5, (double) max_run+0.5, NUMBER_OF_PMTS, -0.5 , NUMBER_OF_PMTS - 0.5);
159  TH2D* h_pmt_distribution_snapshot_hits = new TH2D("h_pmt_distribution_snapshot_hits", string(title + ";RUN; Normalised PMT Distrib snapshot hits (upper PMTs=[0-11])").c_str(), xbins, (double) min_run-0.5, (double) max_run+0.5, NUMBER_OF_PMTS, -0.5 , NUMBER_OF_PMTS - 0.5);
160  TH2D* h_event_duration = new TH2D("h_event_duration", string(title + ";RUN;Event Duration [ns]").c_str(), xbins, (double) min_run-0.5, (double) max_run+0.5, 60, 1, 6);
161  setLogarithmicY(h_event_duration);
162  TH2D* h_tres = new TH2D("h_tres", string(title + ";RUN;Time residuals [ns]").c_str(), xbins, (double) min_run-0.5, (double) max_run+0.5, 100, -50, 150);
163  TH2D* h_likelihood = new TH2D ("h_likelihood", string(title + ";RUN;Likelihood").c_str(), xbins, (double) min_run-0.5, (double) max_run+0.5, 100, -1000, 1000);
164  TH2D* h_beta0 = new TH2D("h_beta0", string(title + ";RUN;beta0").c_str(), xbins, (double) min_run-0.5, (double) max_run+0.5, 20, 0, 1);
165  TH2D* h_energy = new TH2D ("h_energy", string(title + ";RUN;Energy [GeV];").c_str(), xbins, (double) min_run-0.5, (double) max_run+0.5, 65, 0, 9);
166  setLogarithmicY(h_energy);
167  TH2D* h_zenith = new TH2D("h_zenith", string(title + ";RUN;cosZenith").c_str(), xbins, (double) min_run-0.5, (double) max_run+0.5, 20, -1, 1);
168  TH2D* h_azimuth = new TH2D("h_azimuth", string(title + ";RUN;cosAzimuth").c_str(), xbins, (double) min_run-0.5, (double) max_run+0.5, 20, -1, 1);
169  TH2D* h_radial_position = new TH2D ("h_radial_position", string(title + ";RUN;Radial Position [m]").c_str(), xbins, (double) min_run-0.5, (double) max_run+0.5, 60, 0, 4.2);
170  setLogarithmicY(h_radial_position);
171  TH2D* h_z_position = new TH2D ("h_z_position", string(title + ";RUN;Z Position [m]").c_str(), xbins, (double) min_run-0.5, (double) max_run+0.5, 50, 0, 3.2);
172  setLogarithmicY(h_z_position);
173 
174  for (vector<string>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {
175 
176  DEBUG("Processing " << *i << endl) ;
177 
178  {
179  JSingleFileScanner<run_info> in(*i);
180  if(in.hasNext()){
181  runinfo = *in.next();
182  if(in.hasNext()){
183  FATAL("MULTIPLE RUN INFO");
184  }
185  } else {
186  FATAL("NO RUN INFO");
187  }
188  }
189 
190  TFile in(i->c_str(), "read");
191 
192  cout << i->c_str() << endl;
193 
194  int run = runinfo.run_id;
195 
196  TH1D* h1d_prd = new TH1D();
197  TH1D* h1d_thb = new TH1D();
198  TH1D* h1d_th = new TH1D();
199  TH1D* h1d_th3D = new TH1D();
200  TH1D* h1d_sh = new TH1D();
201  TH1D* h1d_tosh = new TH1D();
202  TH1D* h1d_noo = new TH1D();
203  TH1D* h1d_pdth = new TH1D();
204  TH1D* h1d_pdsh = new TH1D();
205  TH1D* h1d_ed = new TH1D();
206  TH1D* h1d_tr = new TH1D();
207  TH1D* h1d_lk = new TH1D();
208  TH1D* h1d_b0 = new TH1D();
209  TH1D* h1d_en = new TH1D();
210  TH1D* h1d_zn = new TH1D();
211  TH1D* h1d_az = new TH1D();
212  TH1D* h1d_rp = new TH1D();
213  TH1D* h1d_zp = new TH1D();
214 
215  in.GetDirectory("Detector")->GetObject("h_pmt_rate_distribution", h1d_prd);
216  fill2DHistogram(run, h1d_prd, h_pmt_rate_distribution, weight_mc, runinfo.daq_livetime, total_mc_livetime);
217 
218  TDirectory* dir = in.GetDirectory("JDAQEvent");
219 
220  dir->GetObject("h_Trigger_bit_hit", h1d_thb);
221  fill2DHistogram(run, h1d_thb, h_Trigger_bit_hit, weight_mc, runinfo.daq_livetime, total_mc_livetime);
222 
223  dir->GetObject("h_Triggered_hits", h1d_th);
224  fill2DHistogram(run, h1d_th, h_Triggered_hits, weight_mc, runinfo.daq_livetime, total_mc_livetime);
225 
226  dir->GetObject("h_Triggered_hits_3dmuon", h1d_th3D);
227  fill2DHistogram(run, h1d_th3D, h_Triggered_hits_3dmuon, weight_mc, runinfo.daq_livetime, total_mc_livetime);
228 
229  dir->GetObject("h_Snapshot_hits", h1d_sh);
230  fill2DHistogram(run, h1d_sh, h_Snapshot_hits, weight_mc, runinfo.daq_livetime, total_mc_livetime);
231 
232  dir->GetObject("h_Triggered_over_Snapshot_hits", h1d_tosh);
233  fill2DHistogram(run, h1d_tosh, h_Triggered_over_Snapshot_hits, weight_mc, runinfo.daq_livetime, total_mc_livetime);
234 
235  dir->GetObject("h_Number_of_overlays", h1d_noo);
236  fill2DHistogram(run, h1d_noo, h_Number_of_overlays, weight_mc, runinfo.daq_livetime, total_mc_livetime);
237 
238  dir->GetObject("h_pmt_distribution_triggered_hits", h1d_pdth);
239  fill2DHistogram(run, h1d_pdth, h_pmt_distribution_triggered_hits, weight_mc, runinfo.daq_livetime, total_mc_livetime);
240 
241  dir->GetObject("h_pmt_distribution_snapshot_hits", h1d_pdsh);
242  fill2DHistogram(run, h1d_pdsh, h_pmt_distribution_snapshot_hits, weight_mc, runinfo.daq_livetime, total_mc_livetime);
243 
244  dir->GetObject("h_event_duration", h1d_ed);
245  fill2DHistogram(run, h1d_ed, h_event_duration, weight_mc, runinfo.daq_livetime, total_mc_livetime);
246 
247  TDirectory* dir_r = in.GetDirectory("Reco");
248 
249  dir_r->GetObject("h_tres", h1d_tr);
250  fill2DHistogram(run, h1d_tr, h_tres, weight_mc, runinfo.daq_livetime, total_mc_livetime);
251 
252  dir_r->GetObject("h_likelihood", h1d_lk);
253  fill2DHistogram(run, h1d_lk, h_likelihood, weight_mc, runinfo.daq_livetime, total_mc_livetime);
254 
255  dir_r->GetObject("h_beta0", h1d_b0);
256  fill2DHistogram(run, h1d_b0, h_beta0, weight_mc, runinfo.daq_livetime, total_mc_livetime);
257 
258  dir_r->GetObject("h_energy", h1d_en);
259  fill2DHistogram(run, h1d_en, h_energy, weight_mc, runinfo.daq_livetime, total_mc_livetime);
260 
261  dir_r->GetObject("h_zenith", h1d_zn);
262  fill2DHistogram(run, h1d_zn, h_zenith, weight_mc, runinfo.daq_livetime, total_mc_livetime);
263 
264  dir_r->GetObject("h_azimuth", h1d_az);
265  fill2DHistogram(run, h1d_az, h_azimuth, weight_mc, runinfo.daq_livetime, total_mc_livetime);
266 
267  dir_r->GetObject("h_radial_position", h1d_rp);
268  fill2DHistogram(run, h1d_rp, h_radial_position, weight_mc, runinfo.daq_livetime, total_mc_livetime);
269 
270  dir_r->GetObject("h_z_position", h1d_zp);
271  fill2DHistogram(run, h1d_zp, h_z_position, weight_mc, runinfo.daq_livetime, total_mc_livetime);
272  }
273 
274  outputFile.open();
275 
276  outputFile.put(*h_event_duration);
277  outputFile.put(*h_pmt_distribution_snapshot_hits);
278  outputFile.put(*h_pmt_distribution_triggered_hits);
279  outputFile.put(*h_Number_of_overlays);
280  outputFile.put(*h_Triggered_over_Snapshot_hits);
281  outputFile.put(*h_Snapshot_hits);
282  outputFile.put(*h_Triggered_hits_3dmuon);
283  outputFile.put(*h_Triggered_hits);
284  outputFile.put(*h_Trigger_bit_hit);
285  outputFile.put(*h_tres);
286  outputFile.put(*h_likelihood);
287  outputFile.put(*h_beta0);
288  outputFile.put(*h_energy);
289  outputFile.put(*h_zenith);
290  outputFile.put(*h_azimuth);
291  outputFile.put(*h_radial_position);
292  outputFile.put(*h_z_position);
293 
294  JMultipleFileScanner<run_info> io(inputFile);
295 
296  io >> outputFile;
297 
298  outputFile.close();
299 }
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
fill2DHistogram
void fill2DHistogram(int run, TH1D *h1, TH2D *h2, bool w_mc, double daq_livetime=1, double mc_livetime=1)
Definition: JMergeRunAnalyzer.cc:45
DEBUG
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
debug
int debug
debug level
Definition: JSirene.cc:69
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:2158
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1714
JROOT::JRootFileReader
ROOT file reader.
Definition: JRootFileReader.hh:54
JSUPPORT::JFileRecorder
Object writing to file.
Definition: JFileRecorder.hh:44
JSUPPORT::JMultipleFileScanner
General purpose class for object reading from a list of file names.
Definition: JMultipleFileScanner.hh:306
JSUPPORT::JSingleFileScanner
Object reading from a list of files.
Definition: JSingleFileScanner.hh:172
KM3NETDAQ::JDAQChronometer::getRunNumber
int getRunNumber() const
Get run number.
Definition: JDAQChronometer.hh:121
KM3NETDAQ::JDAQChronometer::getTimesliceStart
JDAQUTCExtended getTimesliceStart() const
Get start of timeslice.
Definition: JDAQChronometer.hh:144
KM3NETDAQ::JDAQUTCExtended::getUTCseconds
uint32_t getUTCseconds() const
Get major time.
Definition: JDAQUTCExtended.hh:87
JGIZMO::setLogarithmicY
void setLogarithmicY(TList *list)
Make y-axis of objects in list logarithmic (e.g. after using log10()).
Definition: JGizmoToolkit.hh:991
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JAAnetToolkit.hh:43
JTOOLS::getN
double getN(const JRange< T > &range, const double R)
Get expected number of occurrences due to given rate within specified interval.
Definition: JRange.hh:704
KM3NETDAQ
KM3NeT DAQ data structures and auxiliaries.
Definition: DataQueue.cc:39
KM3NETDAQ::NUMBER_OF_TRIGGER_BITS
static const unsigned int NUMBER_OF_TRIGGER_BITS
Number of trigger bits.
Definition: JDAQTriggerMask.hh:27
KM3NETDAQ::NUMBER_OF_PMTS
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
Definition: JDAQ.hh:26
std
Definition: JSTDTypes.hh:14
run_info
Definition: JRunInfo.hh:10
run_info::run_id
int run_id
Definition: JRunInfo.hh:31
run_info::mc_livetime
double mc_livetime
Definition: JRunInfo.hh:33
run_info::daq_livetime
double daq_livetime
Definition: JRunInfo.hh:32
Generated by doxygen 1.9.1