Jpp  15.0.4
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
RunAnalyzer.hh
Go to the documentation of this file.
1 #ifndef __RUNANALYZER__
2 #define __RUNANALYZER__
3 
4 #include "JSupport/JSupport.hh"
9 #include "JROOT/JManager.hh"
10 #include "JTools/JQuantile.hh"
11 #include "TFile.h"
12 #include "TH1D.h"
13 #include "JRunHistograms.hh"
14 
15 using namespace std ;
16 using namespace KM3NETDAQ ;
17 using namespace JPP ;
18 
19 /**
20  * \author rgruiz
21  */
22 
23 /**
24  * Class dedicated to the analysis of %KM3NeT runs.
25  *
26  */
27 class RunAnalyzer {
28 
29  string inputFile;
37 
38 public :
39 
40  /**
41  * Default constructor.
42  */
44 
45  /**
46  * Constructor.
47  *
48  * \param file Path to the file to be analyzed
49  * \param detectorFile path to a detector file
50  * \param nTimeslices Number of timeslices to be read
51  * \param nSummaryslices Number of summary slices to be read
52  * \param nEvents Number of frames to be read
53  */
54  RunAnalyzer (string file , string detectorFile , JLimit_t nTimeslices , JLimit_t nSummaryslices , JLimit_t nEvents):
55  inputFile (file),
56  router (NULL),
57  numberOfTimeslices (nTimeslices),
58  numberOfSummaryslices (nSummaryslices),
59  numberOfEvents (nEvents)
60  {
61  try {
62  load (detectorFile, detector);
63  }
64  catch(const JException & error) {
65  cerr << "FATAL ERROR. Could not open detector file '" << detectorFile << "'." << endl;
66  exit(1);
67  }
68 
69  histograms = JRA_Histograms(detector);
70  router = new JModuleRouter(detector);
72  }
73 
74  /**
75  * Destructor.
76  */
78 
79  /*
80  * Function template used to read time ordered events from a tree.
81  *
82  * \param scanner A JTreeScanner
83  * \param
84  */
86 
87  while (scanner.hasNext()) {
88 
89  JDAQEvent event = *(scanner.next());
90 
91  histograms.h_trigger.h_Snapshot_hits -> Fill((Double_t) event.size<JDAQSnapshotHit > ());
92  histograms.h_trigger.h_Triggered_hits -> Fill((Double_t) event.size<JDAQTriggeredHit> ());
93  histograms.h_trigger.h_Number_of_overlays -> Fill(event.getOverlays());
94 
95  for (unsigned int i = 0; i != NUMBER_OF_TRIGGER_BITS; ++i) {
96 
97  if (event.hasTriggerBit(i)) {
98 
99  histograms.h_trigger.h_Trigger_bit_event -> Fill((Double_t) i);
100  }
101  }
102 
103  histograms.h_trigger.h_Trigger_time->Fill(event.begin<JDAQTriggeredHit>()->getT() * 1e-6);
104 
105  for (JDAQEvent::const_iterator<JDAQTriggeredHit> hit = event.begin<JDAQTriggeredHit>(); hit != event.end<JDAQTriggeredHit>(); ++hit) {
106 
107  if (router->hasModule(hit->getModuleID())) {
108 
109  const JModule& module = router->getModule (hit->getModuleID());
110 
111  histograms.h_trigger.h_pmt_distribution_triggered_hits->Fill(hit->getPMT());
112  histograms.h_trigger.h_tot_distribution_triggered_hits->Fill(hit->getToT());
113 
114  int String = module.getString();
115  int Floor = module.getFloor();
116 
117  histograms.h_trigger.h_Triggered_hits_per_module -> Fill(String , Floor);
118 
119  for (unsigned int i = 0; i != NUMBER_OF_TRIGGER_BITS; ++i) {
120 
121  if (hit -> hasTriggerBit(i)) {
122 
123  histograms.h_trigger.h_Trigger_bit_hit->Fill((Double_t) i);
124  }
125  }
126  }else{
127  FATAL("JModuleRouter trying to access non existing identifier: "<< hit->getModuleID());
128  }
129  }
130 
131  for (JDAQEvent::const_iterator<JDAQSnapshotHit> hit = event.begin<JDAQSnapshotHit>() ; hit != event.end<JDAQSnapshotHit>() ; ++hit){
132 
133  if (router->hasModule(hit->getModuleID())) {
134 
135  const JModule& module = router->getModule (hit->getModuleID());
136 
137  int String = module.getString();
138  int Floor = module.getFloor();
139  int pmt = hit-> getPMT();
140 
141  (*histograms.h_trigger.m_Snapshot_hits_per_pmt)[MAKE_STRING("Detector/DU" + to_string(String))] -> Fill(pmt, Floor);
142  histograms.h_trigger.h_Snapshot_hits_per_module -> Fill(String, Floor);
143  histograms.h_trigger.h_pmt_distribution_snapshot_hits -> Fill(hit->getPMT());
144  histograms.h_trigger.h_tot_distribution_snapshot_hits -> Fill(hit->getToT());
145 
146  }else{
147  FATAL("JModuleRouter trying to access non existing identifier: "<< hit->getModuleID());
148  }
149  }
150  }
151  }
152 
153 
154  /*
155  * Function template to read and process time ordered summary slices.
156  *
157  * \param scanner A JTreeScanner from which time ordered summary slices are accessed.
158  */
160 
161  std::map <int,vector<vector <JQuantile>>> PMT_rate_quantiles;
162  std::map <int,vector<JQuantile>> DOM_rate_quantiles;
163 
164  while (scanner.hasNext()){
165 
166  JDAQSummaryslice slice = *(scanner.next());
167 
168  for (JDAQSummaryslice::const_iterator frame = slice.begin() ; frame != slice.end() ; ++ frame) {
169 
170  if (router->hasModule(frame->getModuleID())) {
171 
172  const JModule& module = router->getModule (frame->getModuleID());
173 
174  int string = module.getString();
175  int floor = module.getFloor ();
176 
177  PMT_rate_quantiles[string].resize(modulesPerString , vector <JQuantile>(NUMBER_OF_PMTS));
178 
179  JDAQFrameStatus status = frame -> getDAQFrameStatus();
180  int nFIFOcount = status.countFIFOStatus();
181  int nHRVcount = status.countHighRateVeto();
182 
183  if (nHRVcount > 0) {
184  histograms.h_summary.h_hrv_per_dom->Fill(string , floor);
185  }
186 
187  if (status.testDAQStatus() == false) {
188  histograms.h_summary.h_daq_status_per_dom->Fill(string , floor);
189  }
190 
191  histograms.h_summary.h_fifo_per_dom->Fill(string , floor , nFIFOcount);
192 
193  TH2D* h2 = (*histograms.h_summary.m_summary_rate_distribution)[MAKE_STRING("Detector/DU" + to_string(string) + "/F" + to_string(floor))];
194  double rate = 0;
195 
196  for (int i = 0 ; i < NUMBER_OF_PMTS ; i++){
197 
198  h2->Fill(i , frame->getRate(i) * 1e-3);
199 
200  rate += frame->getRate(i);
201  histograms.h_summary.h_pmt_rate_distribution->Fill(frame->getRate(i)*1e-3);
202  PMT_rate_quantiles[string][floor-1][i].put(frame->getRate(i)*1e-3);
203  }
204 
205  histograms.h_summary.h_dom_rate_distribution->Fill(rate*1e-3);
206 
207  DOM_rate_quantiles[string].resize(modulesPerString);
208 
209  DOM_rate_quantiles[string][floor-1].put(rate * 1e-3);
210 
211  }else{
212  FATAL("JModuleRouter trying to access non existing identifier: "<< frame->getModuleID());
213  }
214  }
215  }
216  for (std::map<int , vector<JQuantile> >::const_iterator i = DOM_rate_quantiles.begin() ; i!= DOM_rate_quantiles.end() ; ++i) {
217 
218  for (int j=1 ; j < modulesPerString + 1 ; j++) {
219 
220  if (i->second[j-1].getCount() > 0) histograms.h_summary.h_rate_summary->Fill(i->first, j, i->second[j-1].getMean());
221  }
222  }
223 
224  for (std::map<int , vector < vector<JQuantile> > >::const_iterator i = PMT_rate_quantiles.begin() ; i!= PMT_rate_quantiles.end() ; ++i) {
225 
226  for (int j=0 ; j < modulesPerString ; j++){
227 
228  for (int k=0 ; k < NUMBER_OF_PMTS ; k++){
229 
230  if (i -> second[j][k].getCount() > 0){
231 
232  (*histograms.h_summary.m_mean_summary_rate) [MAKE_STRING("Detector/DU" + to_string(i->first))]->Fill(k, j+1, i->second[j][k].getMean());
233  (*histograms.h_summary.m_mean_summary_rate_distribution)[MAKE_STRING("Detector/DU" + to_string(i->first))]->Fill( i->second[j][k].getMean());
234  }
235  }
236  }
237  }
238  }
239 
240  /*
241  * Function template to read and process time ordered timeslices.
242  *
243  * \param scanner A JTreeScanner from which time ordered timeslices are accessed.
244  * \param frame_index_range the range of the indices of the frames to be accessed.
245  */
246  template <class T>
248 
250 
251  std::map < int , std::map <int , JQuantile> > DOM_rate_quantiles;
252 
253  double inverseFrameTimeSec = 1. / (1.0e-9 * getFrameTime());
254 
255  while (scanner.hasNext()){
256 
257  T slice = *(scanner.next());
258 
259  for(auto & frame : slice) {
260  if (router->hasModule(frame.getModuleID())) {
261 
262  const JModule& module = router->getModule (frame.getModuleID());
263  double rate = frame.numberOfHits * inverseFrameTimeSec;
264  int string = module.getString();
265  int floor = module.getFloor ();
266 
267  DOM_rate_quantiles[string][floor].put(rate);
268 
269  vector <int> pmt_hit_count (NUMBER_OF_PMTS , 0) ;
270 
271  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()))];
272 
273  for (JDAQSuperFrame::const_iterator hit = frame.begin() ; hit != frame.end() ; ++hit){
274 
275  h2 -> Fill(hit->getPMT() , hit->getToT()) ;
276 
277  pmt_hit_count[hit->getPMT()]++;
278 
279  }
280 
281  for (int pmt = 0 ; pmt != NUMBER_OF_PMTS ; ++pmt) {
282 
283  (*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);
284 
285  }
286 
287  }else{
288  FATAL("JModuleRouter trying to access non existing identifier: "<< frame.getModuleID());
289  }
290  }
291  }
292 
293  for (std::map< int , std::map<int,JQuantile>>::const_iterator i = DOM_rate_quantiles.begin() ; i != DOM_rate_quantiles.end() ; i++){
294 
295  for (std::map<int , JQuantile>::const_iterator j = i->second.begin() ; j != i->second.end() ; j++){
296 
297  if (j->second.getCount() > 0) histograms.h_timeslice.h_dom_mean_rates[ts_type] -> Fill (i->first , j->first , j->second.getMean() ) ;
298  }
299  }
300  }
301 
302 
303  /*
304  * Checks through the use of a JTreeScanner, if summary slices are present in the run file.
305  * In case they are, the corresponding tree is iterated.
306  *
307  */
309 
310  JTreeScanner <JDAQSummaryslice> scanner(inputFile, numberOfSummaryslices);
311  if (scanner.hasNext()) {
312  histograms.initialize_summary_histograms();
313  scanner.rewind();
314  iterateSummarysliceTree(scanner);
315  }
316  }
317 
318 
319  /*
320  * Function template that checks through the use of a JTreeScanner, if objects of different timeslice classes are present in the run file.
321  * In case they are, the corresponding tree is iterated.
322  *
323  */
324  template <class T>
326 
327  JTreeScanner <T> scanner(inputFile, numberOfTimeslices);
328  if(scanner.hasNext()) {
329 
330  histograms.initialize_timeslice_histograms <T>();
331  scanner.rewind();
332  iterateTimesliceTree(scanner);
333  }
334  }
335 
336 
337  /*
338  * Function template that checks through the use of a JTreeScanner, if JDAQEvent objects are present in the run file.
339  * In case they are, the corresponding tree is read.
340  *
341  */
342  void readEvents() {
343 
344  JTreeScanner <JDAQEvent> scanner(inputFile, numberOfEvents);
345 
346  if(scanner.hasNext()) {
347 
348  histograms.initialize_trigger_histograms();
349  scanner.rewind();
350  iterateEventTree(scanner);
351  }
352  }
353 
354  /*
355  * Returns the histograms produced from the data in the file.
356  */
358 
359  return histograms ;
360  }
361 };
362 
363 #endif
JLimit_t numberOfSummaryslices
Definition: RunAnalyzer.hh:35
JRA_Histograms histograms
Definition: RunAnalyzer.hh:32
General exception.
Definition: JException.hh:23
void readEvents()
Definition: RunAnalyzer.hh:342
DAQ key hit.
Definition: JDAQKeyHit.hh:19
then fatal No hydrophone data file $HYDROPHONE_TXT fi sort gr k
void iterateEventTree(JTreeScanner< JDAQEvent > &scanner)
Definition: RunAnalyzer.hh:85
ROOT TTree parameter settings of various packages.
static const unsigned int NUMBER_OF_TRIGGER_BITS
Number of trigger bits.
void iterateSummarysliceTree(JTreeScanner< JDAQSummaryslice > &scanner)
Definition: RunAnalyzer.hh:159
int countFIFOStatus() const
Count FIFO status.
JLimit_t numberOfTimeslices
Definition: RunAnalyzer.hh:34
int getFloor() const
Get floor number.
Definition: JLocation.hh:145
Data structure for a composite optical module.
Definition: JModule.hh:68
RunAnalyzer(string file, string detectorFile, JLimit_t nTimeslices, JLimit_t nSummaryslices, JLimit_t nEvents)
Constructor.
Definition: RunAnalyzer.hh:54
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.
JLimit_t numberOfEvents
Definition: RunAnalyzer.hh:36
then set_variable singlesRate set_variable doublesRate set_variable numberOfSlices echo Generating random background echo Singles rate
exit
Definition: JPizza.sh:36
Dynamic ROOT object management.
void iterateTimesliceTree(JTreeScanner< T > &scanner)
Definition: RunAnalyzer.hh:247
JTDC_t getT() const
Get time.
Definition: JDAQHit.hh:86
void readTimesliceData()
Definition: RunAnalyzer.hh:325
Template definition for direct access of elements in ROOT TChain.
Definition: JTreeScanner.hh:91
#define MAKE_STRING(A)
Make string.
Definition: JPrint.hh:142
string inputFile
Definition: RunAnalyzer.hh:29
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
Detector file.
Definition: JHead.hh:196
~RunAnalyzer()
Destructor.
Definition: RunAnalyzer.hh:77
Class dedicated to the analysis of KM3NeT runs.
Definition: RunAnalyzer.hh:27
Hit data structure.
Definition: JDAQHit.hh:34
double getFrameTime()
Get frame time duration.
Definition: JDAQClock.hh:162
JRA_Histograms getHistograms()
Definition: RunAnalyzer.hh:357
do set_variable OUTPUT_DIRECTORY $WORKDIR T
int modulesPerString
Definition: RunAnalyzer.hh:33
Support methods.
static const JStringCounter getNumberOfStrings
Function object to count unique strings.
#define FATAL(A)
Definition: JMessage.hh:67
Scanning of objects from multiple files according a format that follows from the extension of each fi...
Direct access to module in detector data structure.
int getString() const
Get string number.
Definition: JLocation.hh:134
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
int countHighRateVeto() const
Count high-rate veto status.
void readSummaryData()
Definition: RunAnalyzer.hh:308
General purpose string class.
Definition: JHead.hh:122
std::string to_string(const T &value)
Convert value to string.
int getCount(const T &hit)
Get hit count.
JModuleRouter * router
Definition: RunAnalyzer.hh:31
int j
Definition: JPolint.hh:666
JDetector detector
Definition: RunAnalyzer.hh:30
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.
int getNumberOfModules(const JDetector &detector)
Get number of modules.
RunAnalyzer()
Default constructor.
Definition: RunAnalyzer.hh:43
then usage $script< string identifier >< detectorfile > event file(toashort file)+" "\nNote that the event files and toashort files should be one-to-one related." fi if (( $