Jpp  15.0.3
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 "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, adomi
21  */
22 
23 /**
24  * Class dedicated to the analysis of %KM3NeT runs.
25  *
26  */
27 class JRunAnalyzer {
28 
29  string inputFile;
38 
39 public :
40 
41  /**
42  * Default constructor.
43  */
45 
46  /**
47  * Constructor.
48  *
49  * \param file Path to the file to be analyzed
50  * \param detectorFile path to a detector file
51  * \param nTimeslices Number of timeslices to be read
52  * \param nSummaryslices Number of summary slices to be read
53  * \param nEvents Number of frames to be read
54  */
55  JRunAnalyzer (string file , string detectorFile , JLimit_t nTimeslices , JLimit_t nSummaryslices , JLimit_t nEvents, bool pmtanalysis):
56  inputFile (file),
57  router (NULL),
58  numberOfTimeslices (nTimeslices),
59  numberOfSummaryslices (nSummaryslices),
60  numberOfEvents (nEvents),
61  pmt_analysis (pmtanalysis)
62  {
63  try {
64  load (detectorFile, detector);
65  }
66  catch(const JException & error) {
67  cerr << "FATAL ERROR. Could not open detector file '" << detectorFile << "'." << endl;
68  exit(1);
69  }
70 
71  histograms = JRA_Histograms(detector);
72  router = new JModuleRouter(detector);
74  }
75 
76  /**
77  * Destructor.
78  */
80 
81  /*
82  * Function template used to read time ordered events from a tree.
83  *
84  * \param scanner A JTreeScanner
85  * \param
86  */
88 
89  while (scanner.hasNext()) {
90 
91  JDAQEvent event = *(scanner.next());
92 
93  histograms.h_trigger.h_Snapshot_hits -> Fill((Double_t) event.size<JDAQSnapshotHit > ());
94  histograms.h_trigger.h_Triggered_hits -> Fill((Double_t) event.size<JDAQTriggeredHit> ());
95  histograms.h_trigger.h_Triggered_over_Snapshot_hits -> Fill((Double_t) event.size<JDAQTriggeredHit> () / event.size<JDAQSnapshotHit > ());
96  histograms.h_trigger.h_Number_of_overlays -> Fill(event.getOverlays());
97 
98  for (unsigned int i = 0; i != NUMBER_OF_TRIGGER_BITS; ++i) {
99 
100  if (event.hasTriggerBit(i)) {
101 
102  histograms.h_trigger.h_Trigger_bit_event -> Fill((Double_t) i);
103  }
104  }
105 
106  int counter_3dmuon = 0;
107 
108  for (JDAQEvent::const_iterator<JDAQTriggeredHit> hit = event.begin<JDAQTriggeredHit>(); hit != event.end<JDAQTriggeredHit>(); ++hit) {
109 
110  int trigger_mask = event.getTriggerMask(*hit);
111 
112  if(trigger_mask == 16 || trigger_mask == 18 ||
113  trigger_mask == 20 || trigger_mask == 22 || trigger_mask == 2.147e+09) counter_3dmuon++;
114 
115  if (router->hasModule(hit->getModuleID())) {
116 
117  const JModule& module = router->getModule (hit->getModuleID());
118 
119  histograms.h_trigger.h_pmt_distribution_triggered_hits->Fill(hit->getPMT());
120  histograms.h_trigger.h_tot_distribution_triggered_hits->Fill(hit->getToT());
121 
122  int String = module.getString();
123  int Floor = module.getFloor();
124 
125  histograms.h_trigger.h_Triggered_hits_per_module -> Fill(String , Floor);
126 
127  if(trigger_mask == 16 || trigger_mask == 18 ||
128  trigger_mask == 20 || trigger_mask == 22 || trigger_mask == 2.147e+09){
129 
130  histograms.h_trigger.h_Triggered_hits_3dmuon_per_module -> Fill(String , Floor);
131  }
132 
133  for (unsigned int i = 0; i != NUMBER_OF_TRIGGER_BITS; ++i) {
134 
135  if (hit -> hasTriggerBit(i)) {
136 
137  histograms.h_trigger.h_Trigger_bit_hit->Fill((Double_t) i);
138  }
139  }
140  }else{
141  FATAL("JModuleRouter trying to access non existing identifier: "<< hit->getModuleID());
142  }
143  }
144 
145  if(counter_3dmuon != 0) histograms.h_trigger.h_Triggered_hits_3dmuon->Fill(counter_3dmuon);
146 
147  for (JDAQEvent::const_iterator<JDAQSnapshotHit> hit = event.begin<JDAQSnapshotHit>() ; hit != event.end<JDAQSnapshotHit>() ; ++hit){
148 
149  if (router->hasModule(hit->getModuleID())) {
150 
151  const JModule& module = router->getModule (hit->getModuleID());
152 
153  int String = module.getString();
154  int Floor = module.getFloor();
155  int pmt = hit-> getPMT();
156 
157  if(pmt_analysis == true){
158 
159  (*histograms.h_trigger.m_Snapshot_hits_per_pmt)[MAKE_STRING("Detector/DU" + to_string(String))] -> Fill(pmt, Floor);
160 
161  }
162 
163  histograms.h_trigger.h_pmt_distribution_snapshot_hits -> Fill(hit->getPMT());
164  histograms.h_trigger.h_tot_distribution_snapshot_hits -> Fill(hit->getToT());
165  histograms.h_trigger.h_Snapshot_hits_per_module -> Fill(String, Floor);
166 
167 
168  }else{
169  FATAL("JModuleRouter trying to access non existing identifier: "<< hit->getModuleID());
170  }
171  }
172  }
173  }
174 
175 
176  /*
177  * Function template to read and process time ordered summary slices.
178  *
179  * \param scanner A JTreeScanner from which time ordered summary slices are accessed.
180  */
182 
183  std::map <int,vector<vector <JQuantile>>> PMT_rate_quantiles;
184  std::map <int,vector<JQuantile>> DOM_rate_quantiles;
185 
186  while (scanner.hasNext()){
187 
188  JDAQSummaryslice slice = *(scanner.next());
189 
190  for (JDAQSummaryslice::const_iterator frame = slice.begin() ; frame != slice.end() ; ++ frame) {
191 
192  if (router->hasModule(frame->getModuleID())) {
193 
194  const JModule& module = router->getModule (frame->getModuleID());
195 
196  int string = module.getString();
197  int floor = module.getFloor ();
198 
199  PMT_rate_quantiles[string].resize(modulesPerString , vector <JQuantile>(NUMBER_OF_PMTS));
200 
201  JDAQFrameStatus status = frame -> getDAQFrameStatus();
202  int nFIFOcount = status.countFIFOStatus();
203  int nHRVcount = status.countHighRateVeto();
204 
205  if (nHRVcount > 0) {
206  histograms.h_summary.h_hrv_per_dom->Fill(string , floor, nHRVcount);
207  }
208 
209  if (status.testDAQStatus() == false) {
210  histograms.h_summary.h_daq_status_per_dom->Fill(string , floor);
211  }
212 
213  histograms.h_summary.h_fifo_per_dom->Fill(string , floor , nFIFOcount);
214 
215  if(pmt_analysis == true){
216 
217  TH2D* h2 = (*histograms.h_summary.m_summary_rate_distribution)[MAKE_STRING("Detector/DU" + to_string(string) + "/F" + to_string(floor))];
218 
219  double rate = 0;
220 
221  for (int i = 0 ; i < NUMBER_OF_PMTS ; i++){
222 
223  rate += frame->getRate(i);
224 
225  h2->Fill(i , frame->getRate(i) * 1e-3);
226 
227  PMT_rate_quantiles[string][floor-1][i].put(frame->getRate(i)*1e-3);
228 
229  }
230 
231  DOM_rate_quantiles[string].resize(modulesPerString);
232 
233  DOM_rate_quantiles[string][floor-1].put(rate * 1e-3);
234 
235  } else {
236 
237  double rate = 0;
238 
239  for (int i = 0 ; i < NUMBER_OF_PMTS ; i++){
240  rate += frame->getRate(i);
241  histograms.h_summary.h_pmt_rate_distribution->Fill(frame->getRate(i)*1e-3);
242  }
243 
244  histograms.h_summary.h_dom_rate_distribution->Fill(rate*1e-3);
245 
246  DOM_rate_quantiles[string].resize(modulesPerString);
247 
248  DOM_rate_quantiles[string][floor-1].put(rate * 1e-3);
249 
250  }
251 
252  }else{
253  FATAL("JModuleRouter trying to access non existing identifier: "<< frame->getModuleID());
254  }
255  }
256  }
257 
258 
259  for (std::map<int , vector<JQuantile> >::const_iterator i = DOM_rate_quantiles.begin() ; i!= DOM_rate_quantiles.end() ; ++i) {
260 
261  for (int j=1 ; j < modulesPerString + 1 ; j++) {
262 
263  if (i->second[j-1].getCount() > 0) histograms.h_summary.h_rate_summary->Fill(i->first, j, i->second[j-1].getMean());
264  }
265  }
266 
267 
268  if(pmt_analysis == true){
269 
270  for (std::map<int , vector < vector<JQuantile> > >::const_iterator i = PMT_rate_quantiles.begin() ; i!= PMT_rate_quantiles.end() ; ++i) {
271 
272  for (int j=0 ; j < modulesPerString ; j++){
273 
274  for (int k=0 ; k < NUMBER_OF_PMTS ; k++){
275 
276  if (i -> second[j][k].getCount() > 0){
277 
278  (*histograms.h_summary.m_mean_summary_rate) [MAKE_STRING("Detector/DU" + to_string(i->first))]->Fill(k, j+1, i->second[j][k].getMean());
279  (*histograms.h_summary.m_mean_summary_rate_distribution)[MAKE_STRING("Detector/DU" + to_string(i->first))]->Fill( i->second[j][k].getMean());
280  }
281  }
282  }
283  }
284  }
285  }
286 
287  /*
288  * Function template to read and process time ordered timeslices.
289  *
290  * \param scanner A JTreeScanner from which time ordered timeslices are accessed.
291  * \param frame_index_range the range of the indices of the frames to be accessed.
292  */
293  template <class T>
295 
297 
298  std::map < int , std::map <int , JQuantile> > DOM_rate_quantiles;
299 
300  double inverseFrameTimeSec = 1. / (1.0e-9 * getFrameTime());
301 
302  while (scanner.hasNext()){
303 
304  T slice = *(scanner.next());
305 
306  for(auto & frame : slice) {
307  if (router->hasModule(frame.getModuleID())) {
308 
309  const JModule& module = router->getModule (frame.getModuleID());
310  double rate = frame.numberOfHits * inverseFrameTimeSec;
311  int string = module.getString();
312  int floor = module.getFloor ();
313 
314  DOM_rate_quantiles[string][floor].put(rate);
315 
316  vector <int> pmt_hit_count (NUMBER_OF_PMTS , 0) ;
317 
318  if(pmt_analysis == true){
319 
320  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()))];
321 
322  for (JDAQSuperFrame::const_iterator hit = frame.begin() ; hit != frame.end() ; ++hit){
323 
324  h2 -> Fill(hit->getPMT() , hit->getToT()) ;
325 
326  pmt_hit_count[hit->getPMT()]++;
327 
328  }
329 
330  for (int pmt = 0 ; pmt != NUMBER_OF_PMTS ; ++pmt) {
331 
332  (*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);
333 
334  }
335 
336  }
337 
338  }else{
339  FATAL("JModuleRouter trying to access non existing identifier: "<< frame.getModuleID());
340  }
341  }
342  }
343 
344 
345  for (std::map< int , std::map<int,JQuantile>>::const_iterator i = DOM_rate_quantiles.begin() ; i != DOM_rate_quantiles.end() ; i++){
346 
347  for (std::map<int , JQuantile>::const_iterator j = i->second.begin() ; j != i->second.end() ; j++){
348 
349  if (j->second.getCount() > 0) histograms.h_timeslice.h_dom_mean_rates[ts_type] -> Fill (i->first , j->first , j->second.getMean() ) ;
350  }
351  }
352  }
353 
354 
355  /*
356  * Checks through the use of a JTreeScanner, if summary slices are present in the run file.
357  * In case they are, the corresponding tree is iterated.
358  *
359  */
361 
362  JTreeScanner <JDAQSummaryslice> scanner(inputFile, numberOfSummaryslices);
363  if (scanner.hasNext()) {
364  histograms.initialize_summary_histograms();
365  scanner.rewind();
366  iterateSummarysliceTree(scanner);
367  }
368  }
369 
370 
371  /*
372  * Function template that checks through the use of a JTreeScanner, if objects of different timeslice classes are present in the run file.
373  * In case they are, the corresponding tree is iterated.
374  *
375  */
376  template <class T>
378 
379  JTreeScanner <T> scanner(inputFile, numberOfTimeslices);
380  if(scanner.hasNext()) {
381 
382  histograms.initialize_timeslice_histograms <T>();
383  scanner.rewind();
384  iterateTimesliceTree(scanner);
385  }
386  }
387 
388 
389  /*
390  * Function template that checks through the use of a JTreeScanner, if JDAQEvent objects are present in the run file.
391  * In case they are, the corresponding tree is read.
392  *
393  */
394  void readEvents() {
395 
396  JTreeScanner <JDAQEvent> scanner(inputFile, numberOfEvents);
397 
398  if(scanner.hasNext()) {
399 
400  histograms.initialize_trigger_histograms();
401  scanner.rewind();
402  iterateEventTree(scanner);
403  }
404  }
405 
406  /*
407  * Returns the histograms produced from the data in the file.
408  */
410 
411  return histograms ;
412  }
413 };
414 
415 #endif
General exception.
Definition: JException.hh:23
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.
JRunAnalyzer(string file, string detectorFile, JLimit_t nTimeslices, JLimit_t nSummaryslices, JLimit_t nEvents, bool pmtanalysis)
Constructor.
Definition: JRunAnalyzer.hh:55
int getFloor() const
Get floor number.
Definition: JLocation.hh:145
Data structure for a composite optical module.
Definition: JModule.hh:68
int modulesPerString
Definition: JRunAnalyzer.hh:33
JRA_Histograms histograms
Definition: JRunAnalyzer.hh:32
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
exit
Definition: JPizza.sh:36
JDetector detector
Definition: JRunAnalyzer.hh:30
string inputFile
Definition: JRunAnalyzer.hh:29
JLimit_t numberOfTimeslices
Definition: JRunAnalyzer.hh:34
Dynamic ROOT object management.
JModuleRouter * router
Definition: JRunAnalyzer.hh:31
void iterateSummarysliceTree(JTreeScanner< JDAQSummaryslice > &scanner)
Class dedicated to the analysis of KM3NeT runs.
Definition: JRunAnalyzer.hh:27
Template definition for direct access of elements in ROOT TChain.
Definition: JTreeScanner.hh:91
#define MAKE_STRING(A)
Make string.
Definition: JPrint.hh:142
void readTimesliceData()
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
JRunAnalyzer()
Default constructor.
Definition: JRunAnalyzer.hh:44
~JRunAnalyzer()
Destructor.
Definition: JRunAnalyzer.hh:79
Detector file.
Definition: JHead.hh:224
Hit data structure.
Definition: JDAQHit.hh:34
double getFrameTime()
Get frame time duration.
Definition: JDAQClock.hh:162
do set_variable OUTPUT_DIRECTORY $WORKDIR T
Support methods.
void readSummaryData()
void iterateEventTree(JTreeScanner< JDAQEvent > &scanner)
Definition: JRunAnalyzer.hh:87
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.
JLimit_t numberOfEvents
Definition: JRunAnalyzer.hh:36
General purpose string class.
Definition: JHead.hh:150
std::string to_string(const T &value)
Convert value to string.
int getCount(const T &hit)
Get hit count.
JLimit_t numberOfSummaryslices
Definition: JRunAnalyzer.hh:35
int j
Definition: JPolint.hh:666
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()
int getNumberOfModules(const JDetector &detector)
Get number of modules.