Jpp  17.0.0
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JRunAnalyzer.hh
Go to the documentation of this file.
1 #ifndef __JRUNANALYZER__
2 #define __JRUNANALYZER__
3 
4 #include "JSupport/JSupport.hh"
9 #include "JROOT/JManager.hh"
10 #include "JTools/JQuantile.hh"
11 #include "JTools/JRange.hh"
12 #include "JSupport/JTreeScanner.hh"
14 #include "TFile.h"
15 #include "TH1D.h"
16 #include "JRunHistograms.hh"
17 
18 /**
19  * \author rgruiz, adomi
20  */
21 
23 
24 /**
25  * Class dedicated to the analysis of %KM3NeT runs.
26  *
27  */
28 class JRunAnalyzer {
29 
38  // JDETECTOR::JDetector detector;
39 
40 public :
41 
42  /**
43  * Constructor.
44  *
45  * \param file file name(s)
46  * \param detector detector
47  * \param nTimeslices number of timeslices to be read
48  * \param nSummaryslices number of summary slices to be read
49  * \param nEvents number of events to be read
50  * \param pmtanalysis option for analysis of PMT data
51  */
52  JRunAnalyzer (const JMultipleFileScanner<>& file, const JDetector& detector , JLimit_t nTimeslices , JLimit_t nSummaryslices , JLimit_t nEvents, bool pmtanalysis):
53  inputFile (file),
54  router (detector),
55  numberOfTimeslices (nTimeslices),
56  numberOfSummaryslices (nSummaryslices),
57  numberOfEvents (nEvents),
58  pmt_analysis (pmtanalysis)
59  {
60  using namespace JPP;
61  using namespace std;
62 
63  histograms = JRA_Histograms(detector);
65  }
66 
67  /**
68  * Destructor.
69  */
71 
72  /*
73  * Function template used to read events from a tree.
74  *
75  * \param scanner A JTreeScanner
76  * \param
77  */
79 
80  while (scanner.hasNext()) {
81 
82  JDAQEvent event = *(scanner.next());
83 
84  histograms.h_trigger.h_Snapshot_hits -> Fill((Double_t) event.size<JDAQSnapshotHit> ());
85  histograms.h_trigger.h_Triggered_hits -> Fill((Double_t) event.size<JDAQTriggeredHit> ());
86  histograms.h_trigger.h_Triggered_over_Snapshot_hits -> Fill((Double_t) event.size<JDAQTriggeredHit> () / event.size<JDAQSnapshotHit > ());
87  histograms.h_trigger.h_Number_of_overlays -> Fill(event.getOverlays());
88 
89  for (unsigned int i = 0; i != NUMBER_OF_TRIGGER_BITS; ++i) {
90 
91  if (event.hasTriggerBit(i)) {
92 
93  histograms.h_trigger.h_Trigger_bit_event -> Fill((Double_t) i);
94  }
95  }
96 
97  int counter_3dmuon = 0;
98  typedef JRange<double> JRange_t;
99  JRange_t range(JRange_t::DEFAULT_RANGE);
100 
101  for (JDAQEvent::const_iterator<JDAQTriggeredHit> hit = event.begin<JDAQTriggeredHit>(); hit != event.end<JDAQTriggeredHit>(); ++hit) {
102 
103  const JModule& module = router.getModule(hit->getModuleID());
104  const double t1 = getTime(hit->getT(), module.getPMT(hit->getPMT()));
105 
106  range.include(t1);
107 
108  JDAQTriggerMask trigger_mask(event.getTriggerMask(*hit));
109 
110  if(trigger_mask.hasTriggerBit(getTriggerBit<JTrigger3DMuon>())) counter_3dmuon++;
111 
112  if (router.hasModule(hit->getModuleID())) {
113 
114  const JModule& module = router.getModule (hit->getModuleID());
115 
118 
119  int String = module.getString();
120  int Floor = module.getFloor();
121 
122  histograms.h_trigger.h_Triggered_hits_per_module -> Fill(String , Floor);
123 
124  if(trigger_mask.hasTriggerBit(getTriggerBit<JTrigger3DMuon>())){
126  }
127 
128  for (unsigned int i = 0; i != NUMBER_OF_TRIGGER_BITS; ++i) {
129 
130  if (hit -> hasTriggerBit(i)) {
131 
132  histograms.h_trigger.h_Trigger_bit_hit->Fill((Double_t) i);
133  }
134  }
135  }else{
136  FATAL("JModuleRouter trying to access non existing identifier: "<< hit->getModuleID());
137  }
138  }
139 
141 
142  histograms.h_trigger.h_pmt_distribution_triggered_hits->Scale(1. / (Double_t) event.size<JDAQTriggeredHit> ());
143 
144  if(counter_3dmuon != 0) histograms.h_trigger.h_Triggered_hits_3dmuon->Fill(counter_3dmuon);
145 
146  for (JDAQEvent::const_iterator<JDAQSnapshotHit> hit = event.begin<JDAQSnapshotHit>() ; hit != event.end<JDAQSnapshotHit>() ; ++hit){
147 
148  if (router.hasModule(hit->getModuleID())) {
149 
150  const JModule& module = router.getModule (hit->getModuleID());
151 
152  int String = module.getString();
153  int Floor = module.getFloor();
154  int pmt = hit-> getPMT();
155 
156  if(pmt_analysis == true){
157 
158  (*histograms.h_trigger.m_Snapshot_hits_per_pmt)[MAKE_STRING("Detector/DU" + to_string(String))] -> Fill(pmt, Floor);
159 
160  }
161 
162  histograms.h_trigger.h_pmt_distribution_snapshot_hits -> Fill(hit->getPMT());
164  histograms.h_trigger.h_Snapshot_hits_per_module -> Fill(String, Floor);
165 
166 
167  }else{
168  FATAL("JModuleRouter trying to access non existing identifier: "<< hit->getModuleID());
169  }
170  }
171  }
172  }
173 
174 
175 
176  struct dom_type :
177  public std::vector< JQuantile>
178  {
181  {}
182  };
183 
184  /*
185  * Function template to read and process time ordered summary slices.
186  *
187  * \param scanner A JTreeScanner from which time ordered summary slices are accessed.
188  */
190 
191  using namespace std;
192  using namespace JPP;
193 
194  map<int, map<int, dom_type> > PMT_rate_quantiles;
195  map<int, map<int, JQuantile> > DOM_rate_quantiles;
196 
197  while (scanner.hasNext()){
198 
199  JDAQSummaryslice slice = *(scanner.next());
200 
201  for (JDAQSummaryslice::const_iterator frame = slice.begin() ; frame != slice.end() ; ++ frame) {
202 
203  if (router.hasModule(frame->getModuleID())) {
204 
205  const JModule& module = router.getModule (frame->getModuleID());
206 
207  int string = module.getString();
208  int floor = module.getFloor ();
209 
210  JDAQFrameStatus status = frame -> getDAQFrameStatus();
211  int nFIFOcount = status.countFIFOStatus();
212  int nHRVcount = status.countHighRateVeto();
213 
214  if (nHRVcount > 0) {
215  histograms.h_summary.h_hrv_per_dom->Fill(string , floor, nHRVcount);
216  }
217 
218  if (status.testDAQStatus() == false) {
219  histograms.h_summary.h_daq_status_per_dom->Fill(string , floor, (1./distance(scanner.begin(), scanner.end()))*100);
220  }
221 
222  histograms.h_summary.h_fifo_per_dom->Fill(string , floor , nFIFOcount);
223 
224  if(pmt_analysis == true){
225 
226  TH2D* h2 = (*histograms.h_summary.m_summary_rate_distribution)[MAKE_STRING("Detector/DU" + to_string(string) + "/F" + to_string(floor))];
227 
228  double rate = 0;
229  const double factor = 1.0e-3; // [kHz]
230 
231  for (int i = 0 ; i < NUMBER_OF_PMTS ; i++){
232 
233  rate += frame->getRate(i, factor);
234 
235  h2->Fill(i , frame->getRate(i, factor));
236 
237  PMT_rate_quantiles[string][floor][i].put(frame->getRate(i, factor));
238 
239  histograms.h_summary.h_pmt_rate_distribution->Fill(frame->getRate(i, factor), frame->getWeight(i, factor));
240 
241  }
242 
244 
245  DOM_rate_quantiles[string][floor].put(rate);
246 
247  } else {
248 
249  double rate = 0;
250  const double factor = 1.0e-3; // [kHz]
251 
252  for (int i = 0 ; i < NUMBER_OF_PMTS ; i++){
253 
254  rate += frame->getRate(i, factor);
255 
256  histograms.h_summary.h_pmt_rate_distribution->Fill(frame->getRate(i, factor), frame->getWeight(i, factor));
257  }
258 
260 
261  DOM_rate_quantiles[string][floor].put(rate);
262 
263  }
264 
265  }else{
266  FATAL("JModuleRouter trying to access non existing identifier: "<< frame->getModuleID());
267  }
268  }
269  }
270 
271  for(auto& i : DOM_rate_quantiles){
272 
273  for (int j=1; j <= highest_floor; j++) {
274 
275  if (i.second[j].getCount() > 0) histograms.h_summary.h_rate_summary->Fill(i.first, j, i.second[j].getMean());
276  }
277  }
278 
279  if(pmt_analysis == true){
280 
281  for(auto& i : PMT_rate_quantiles){
282 
283  for (int j = 1; j <= highest_floor; j++){
284 
285  for (int k=0 ; k < NUMBER_OF_PMTS ; k++){
286 
287  if (i.second[j][k].getCount() > 0){
288 
289  (*histograms.h_summary.m_mean_summary_rate) [MAKE_STRING("Detector/DU" + to_string(i.first))]->Fill(k, j, i.second[j][k].getMean());
290  (*histograms.h_summary.m_mean_summary_rate_distribution)[MAKE_STRING("Detector/DU" + to_string(i.first))]->Fill( i.second[j][k].getMean());
291  }
292  }
293  }
294  }
295  }
296  }
297 
298  /*
299  * Function template to read and process time ordered timeslices.
300  *
301  * \param scanner A JTreeScanner from which time ordered timeslices are accessed.
302  */
303  template <class T>
305 
307 
308  std::map < int , std::map <int , JQuantile> > DOM_rate_quantiles;
309 
310  double inverseFrameTimeSec = 1. / (1.0e-9 * getFrameTime());
311 
312  while (scanner.hasNext()){
313 
314  T slice = *(scanner.next());
315 
316  for(auto & frame : slice) {
317  if (router.hasModule(frame.getModuleID())) {
318 
319  const JModule& module = router.getModule (frame.getModuleID());
320  double rate = frame.numberOfHits * inverseFrameTimeSec;
321  int string = module.getString();
322  int floor = module.getFloor ();
323 
324  DOM_rate_quantiles[string][floor].put(rate);
325 
326  vector <int> pmt_hit_count (NUMBER_OF_PMTS , 0) ;
327 
328  if(pmt_analysis == true){
329 
330  TH2D* h2 = (*histograms.h_timeslice.m_pmt_tot_distributions[ts_type])[MAKE_STRING("Detector/DU" + to_string(string) + "/F" + to_string(floor) + "/" + to_string(frame.getModuleID()))];
331 
332  for (JDAQSuperFrame::const_iterator hit = frame.begin() ; hit != frame.end() ; ++hit){
333 
334  h2 -> Fill(hit->getPMT() , hit->getToT()) ;
335 
336  pmt_hit_count[hit->getPMT()]++;
337 
338  }
339 
340  for (int pmt = 0 ; pmt != NUMBER_OF_PMTS ; ++pmt) {
341 
342  (*histograms.h_timeslice.m_pmt_rate_distributions[ts_type])[MAKE_STRING("Detector/DU" + to_string(string) + "/F" + to_string(floor))] -> Fill(pmt , 1e-3 * pmt_hit_count[pmt] * inverseFrameTimeSec);
343 
344  }
345 
346  }
347 
348  }else{
349  FATAL("JModuleRouter trying to access non existing identifier: "<< frame.getModuleID());
350  }
351  }
352  }
353 
354 
355  for (std::map< int , std::map<int,JQuantile>>::const_iterator i = DOM_rate_quantiles.begin() ; i != DOM_rate_quantiles.end() ; i++){
356 
357  for (std::map<int , JQuantile>::const_iterator j = i->second.begin() ; j != i->second.end() ; j++){
358 
359  if (j->second.getCount() > 0) histograms.h_timeslice.h_dom_mean_rates[ts_type] -> Fill (i->first , j->first , j->second.getMean() ) ;
360  }
361  }
362  }
363 
364 
365  /*
366  * Checks through the use of a JTreeScanner, if summary slices are present in the run file.
367  * In case they are, the corresponding tree is iterated.
368  *
369  */
371 
373  if (scanner.hasNext()) {
375  scanner.rewind();
376  iterateSummarysliceTree(scanner);
377  }
378  }
379 
380 
381  /*
382  * Function template that checks through the use of a JTreeScanner, if objects of different timeslice classes are present in the run file.
383  * In case they are, the corresponding tree is iterated.
384  *
385  */
386  template <class T>
388 
390  if(scanner.hasNext()) {
391 
393  scanner.rewind();
394  iterateTimesliceTree(scanner);
395  }
396  }
397 
398 
399  /*
400  * Function template that checks through the use of a JTreeScanner, if JDAQEvent objects are present in the run file.
401  * In case they are, the corresponding tree is read.
402  *
403  */
404  void readEvents() {
405 
407 
408  if(scanner.hasNext()) {
409 
411  scanner.rewind();
412  iterateEventTree(scanner);
413  }
414  }
415 
416  /*
417  * Returns the histograms produced from the data in the file.
418  */
420 
421  return histograms ;
422  }
423 
424  /*
425  * Writes the histograms to a root file.
426  * \param f The root file.
427  */
428  void writeToFile(TFile & f){
429 
430  f.mkdir("Detector");
431  f .cd("Detector");
432 
439 
441 
442  for (typename vector < JManager < string , TH2D >* >::const_iterator i = histograms.h_timeslice.m_pmt_tot_distributions.begin();
444 
445  if ((*i)){
446 
447  for (typename JManager < string , TH2D >::const_iterator j = (*i) -> begin() ; j != (*i) -> end() ; ++j){
448 
449  j->second->Scale(1., "width");
450  }
451  }
452  }
453 
456 
458 
461 
463 
466 
467  f.mkdir ( MAKE_STRING ("JDAQEvent").c_str());
468  f.cd ("JDAQEvent");
469 
486  }
487 };
488 
489 #endif
TH1D * h_pmt_distribution_snapshot_hits
TH2D * h_Snapshot_hits_per_module
DAQ key hit.
Definition: JDAQKeyHit.hh:19
then fatal No hydrophone data file $HYDROPHONE_TXT fi sort gr k
void iterateTimesliceTree(JTreeScanner< T > &scanner)
ROOT TTree parameter settings of various packages.
static const unsigned int NUMBER_OF_TRIGGER_BITS
Number of trigger bits.
int countFIFOStatus() const
Count FIFO status.
TH1D * h_pmt_distribution_triggered_hits
vector< JManager< string, TH2D > * > m_pmt_tot_distributions
JManager< string, TH1D > * m_mean_summary_rate_distribution
int getFloor() const
Get floor number.
Definition: JLocation.hh:145
const JModule & getModule(const JObjectID &id) const
Get module parameters.
Data structure for a composite optical module.
Definition: JModule.hh:68
Auxiliary class for trigger mask.
range_type & include(argument_type x)
Include given value to range.
Definition: JRange.hh:397
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
vector< TH2D * > h_dom_mean_rates
JRA_Histograms histograms
Definition: JRunAnalyzer.hh:32
TH2D * h_Triggered_hits_per_module
Auxiliary data structure for running average, standard deviation and quantiles.
Definition: JQuantile.hh:43
Detector data structure.
Definition: JDetector.hh:89
Template const_iterator.
Definition: JDAQEvent.hh:68
Router for direct addressing of module data in detector data structure.
then set_variable singlesRate set_variable doublesRate set_variable numberOfSlices echo Generating random background echo Singles rate
o $QUALITY_ROOT d $DEBUG!JPlot1D f
Definition: JDataQuality.sh:66
JLimit_t numberOfTimeslices
Definition: JRunAnalyzer.hh:34
Dynamic ROOT object management.
void initialize_summary_histograms()
double getTime(const Hit &hit)
Get true time of hit.
void iterateSummarysliceTree(JTreeScanner< JDAQSummaryslice > &scanner)
JManager< string, TH2D > * m_mean_summary_rate
Class dedicated to the analysis of KM3NeT runs.
Definition: JRunAnalyzer.hh:28
floor_range getRangeOfFloors(const JDetector &detector)
Get range of floors.
Template definition for direct access of elements in ROOT TChain.
TH1D * h_dom_rate_distribution
#define MAKE_STRING(A)
Make string.
Definition: JPrint.hh:142
void writeToFile(TFile &f)
void readTimesliceData()
T getLength() const
Get length (difference between upper and lower limit).
Definition: JRange.hh:289
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
Type definition of range.
Definition: JHead.hh:39
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys...
Definition: JManager.hh:43
TimesliceHistograms h_timeslice
~JRunAnalyzer()
Destructor.
Definition: JRunAnalyzer.hh:70
JModuleRouter router
Definition: JRunAnalyzer.hh:31
Detector file.
Definition: JHead.hh:224
TH1D * h_tot_distribution_snapshot_hits
Hit data structure.
Definition: JDAQHit.hh:34
TriggerHistograms h_trigger
double getFrameTime()
Get frame time duration.
Definition: JDAQClock.hh:162
do set_variable OUTPUT_DIRECTORY $WORKDIR T
void Write_manager_table_in_key_dir(TFile &f, vector< JManager< T, V > * > table)
T getUpperLimit() const
Get upper limit.
Definition: JRange.hh:213
Support methods.
TH1D * h_pmt_rate_distribution
JMultipleFileScanner inputFile
Definition: JRunAnalyzer.hh:30
void Write_histogram_table_to_file(TFile &f, string dirname, vector< vector< T * > > table)
void readSummaryData()
const JPMT & getPMT(const int index) const
Get PMT.
Definition: JModule.hh:173
void iterateEventTree(JTreeScanner< JDAQEvent > &scanner)
Definition: JRunAnalyzer.hh:78
vector< JManager< string, TH1D > * > m_mean_ToT_distribution
#define FATAL(A)
Definition: JMessage.hh:67
Scanning of objects from multiple files according a format that follows from the extension of each fi...
SummaryHistograms h_summary
Direct access to module in detector data structure.
z range($ZMAX-$ZMIN)< $MINIMAL_DZ." fi fi typeset -Z 4 STRING typeset -Z 2 FLOOR JPlot1D -f $
int getString() const
Get string number.
Definition: JLocation.hh:134
Auxiliary class to define a range between two values.
int countHighRateVeto() const
Count high-rate veto status.
General purpose class for object reading from a list of file names.
JManager< string, TH2D > * m_summary_rate_distribution
JLimit_t numberOfEvents
Definition: JRunAnalyzer.hh:36
void Write_manager_in_key_dir(TFile &f, JManager< T, V > *manager)
General purpose string class.
Definition: JHead.hh:150
std::string to_string(const T &value)
Convert value to string.
vector< JManager< string, TH2D > * > m_mean_ToT
vector< JManager< string, TH2D > * > m_pmt_rate_distributions
bool hasModule(const JObjectID &id) const
Has module.
JRunAnalyzer(const JMultipleFileScanner<> &file, const JDetector &detector, JLimit_t nTimeslices, JLimit_t nSummaryslices, JLimit_t nEvents, bool pmtanalysis)
Constructor.
Definition: JRunAnalyzer.hh:52
TH1D * h_tot_distribution_triggered_hits
TH1D * h_Triggered_hits_3dmuon
JLimit_t numberOfSummaryslices
Definition: JRunAnalyzer.hh:35
int j
Definition: JPolint.hh:682
TH2D * h_Triggered_hits_3dmuon_per_module
Indexing of data type in type list.
Definition: JTypeList.hh:310
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
Definition: JDAQ.hh:26
bool testDAQStatus() const
Test DAQ status of packets.
JRA_Histograms getHistograms()
void readEvents()
void initialize_trigger_histograms()
JManager< string, TH2D > * m_Snapshot_hits_per_pmt
void initialize_timeslice_histograms()
TH1D * h_Triggered_over_Snapshot_hits