Jpp  pmt_effective_area_update_2
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
JHalibut.cc File Reference

Example program to determine N-fold coincidence rates. More...

#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <vector>
#include <map>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TF1.h"
#include "JDAQ/JDAQTimesliceIO.hh"
#include "km3net-dataformat/online/JDAQClock.hh"
#include "JDAQ/JDAQEvaluator.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.hh"
#include "JTrigger/JHitR0.hh"
#include "JTrigger/JMatchL0.hh"
#include "JTrigger/JBuildL0.hh"
#include "JTrigger/JSuperFrame1D.hh"
#include "JTrigger/JSuperFrame2D.hh"
#include "JROOT/JManager.hh"
#include "JCalibrate/JCalibrateK40.hh"
#include "JROOT/JROOTClassSelector.hh"
#include "JTools/JQuantile.hh"
#include "JTools/JRange.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JTreeScanner.hh"
#include "JSupport/JAutoTreeScanner.hh"
#include "JSupport/JSupport.hh"
#include "JLang/JObjectMultiplexer.hh"
#include "JMath/JMathToolkit.hh"
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Detailed Description

Example program to determine N-fold coincidence rates.

Author
mdejong

Definition in file JHalibut.cc.

Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 52 of file JHalibut.cc.

53 {
54  using namespace std;
55  using namespace JPP;
56  using namespace KM3NETDAQ;
57 
58  typedef JRange<int> JRange_t;
59 
60  JMultipleFileScanner<> inputFile;
61  JLimit_t& numberOfEvents = inputFile.getLimit();
62  string outputFile;
63  string detectorFile;
64  double TMax_ns;
65  JROOTClassSelector selector;
66  JRange_t M;
67  string summaryFile;
68  string option;
69  int debug;
70 
71  try {
72 
73  JParser<> zap("Example program to determine N-fold coincidence rates.");
74 
75  zap['f'] = make_field(inputFile);
76  zap['o'] = make_field(outputFile) = "halibut.root";
77  zap['n'] = make_field(numberOfEvents) = JLimit::max();
78  zap['a'] = make_field(detectorFile);
79  zap['T'] = make_field(TMax_ns) = 20.0;
80  zap['C'] = make_field(selector) = getROOTClassSelection<JDAQTimesliceTypes_t>();
81  zap['M'] = make_field(M) = JRange_t(2, 31);
82  zap['s'] = make_field(summaryFile) = "halibut.txt";
83  zap['O'] = make_field(option) = "";
84  zap['d'] = make_field(debug) = 1;
85 
86  zap(argc, argv);
87  }
88  catch(const exception &error) {
89  FATAL(error.what() << endl);
90  }
91 
92  cout.tie(&cerr);
93 
94 
96 
97  try {
98  load(detectorFile, detector);
99  }
100  catch(const JException& error) {
101  FATAL(error);
102  }
103 
104  const JModuleRouter router(detector);
105 
106  map<int, double> TMax_s; // histogram upper limit
107 
108  TMax_s[2] = 5.0e-2;
109  TMax_s[3] = 0.2e+0;
110  TMax_s[4] = 1.0e+0;
111  TMax_s[5] = 1.0e+1;
112  TMax_s[6] = 1.0e+2;
113 
114  typedef JSuperFrame2D<JHitR0> JSuperFrame2D_t;
115  typedef JSuperFrame1D<JHitR0> JSuperFrame1D_t;
116 
117  const JMatchL0<JHitR0> match(TMax_ns); // time window self-coincidences [ns]
118 
119  typedef JManager<int, TH1D> JManager_t;
120 
121  JManager_t H1(new TH1D("H1[%]", NULL, 100, -TMax_ns, +TMax_ns));
122  JManager_t T1(new TH1D("T1[%]", NULL, 100, 0.0, TMax_s[M.getLowerLimit()]));
123 
124  H1->Sumw2();
125  T1->Sumw2();
126 
128 
130 
131  JTreeScannerInterface<JDAQTimeslice>* ps = zmap[selector];
132 
133  counter_type counter = 0;
134 
135  // process file-by-file to speed up JTreeScanner::configure();
136 
137  for (JMultipleFileScanner<>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {
138 
139  STATUS("Processing " << *i << endl);
140 
141  ps->configure(*i);
142 
143  vector<double> t0(detector.size(), 0.0);
144 
145  for ( ; ps->hasNext() && counter != inputFile.getLimit(); ++counter) {
146 
147  STATUS("event: " << setw(10) << counter << '\r'); DEBUG(endl);
148 
149  const JDAQTimeslice* timeslice = ps->next();
150 
151  for (JDAQTimeslice::const_iterator frame = timeslice->begin(); frame != timeslice->end(); ++frame) {
152 
153  JSuperFrame2D_t& buffer = JSuperFrame2D_t::demultiplex(*frame, router.getModule(frame->getModuleID()));
154 
155  for (JSuperFrame2D_t::iterator i = buffer.begin(); i != buffer.end(); ++i) {
156  i->join(match);
157  }
158 
159  JSuperFrame1D_t& data = JSuperFrame1D_t::multiplex(buffer);
160 
161  data.pop_back(); // remove end marker
162 
163  if (data.size() > 1) {
164 
165  TH1D* h1 = H1[frame->getModuleID()];
166  TH1D* t1 = T1[frame->getModuleID()];
167 
168  for (vector<JHitR0>::const_iterator p = data.begin(); p != data.end(); ) {
169 
171 
172  while (++q != data.end() && q->getT() - p->getT() <= TMax_ns ) {}
173 
174  const int N = distance(p,q);
175 
176  if (M(N)) {
177 
178  const int i = router.getIndex(frame->getModuleID());
179  const double ts = getTimeOfRTS(frame->getFrameIndex()) + p->getT();
180 
181  t1->Fill((ts - t0[i]) * 1.0e-9); // [s]
182 
183  t0[i] = ts;
184 
185  const double W = factorial(N, 2);
186 
187  for (vector<JHitR0>::const_iterator __p = p; __p != q; ++__p) {
188  for (vector<JHitR0>::const_iterator __q = __p; ++__q != q; ) {
189  h1->Fill(JCombinatorics::getSign(__p->getPMT(),__q->getPMT()) * (__q->getT() - __p->getT()), 1.0/W);
190  }
191  }
192  }
193 
194  p = q;
195  }
196  }
197  }
198  }
199  STATUS(endl);
200  }
201 
202  if (counter != 0) {
203 
204  const double V = (H1->GetXaxis()->GetXmax() - H1->GetXaxis()->GetXmin()) / (double) H1->GetXaxis()->GetNbins(); // [ns]
205  const double W = counter * getFrameTime() * 1.0e-9; // [s]
206 
207  for (JManager_t::iterator i = H1.begin(); i != H1.end(); ++i) {
208  i->second->Scale(1.0/(V*W));
209  }
210 
211  for (JManager_t::iterator i = T1.begin(); i != T1.end(); ++i) {
212  i->second->Scale(1.0/i->second->GetMaximum());
213  }
214  }
215 
216  if (summaryFile != "") {
217 
218  const double V = (H1->GetXaxis()->GetXmax() - H1->GetXaxis()->GetXmin()) / (double) H1->GetXaxis()->GetNbins(); // [ns]
219 
220  const int number_of_strings = getNumberOfStrings(detector);
221  const int number_of_floors = getNumberOfFloors (detector);
222  const int PRECISION = (M.getLowerLimit() > 2 ? 4 : 3);
223 
224  ofstream out(summaryFile.c_str());
225 
226  out << "Multiplicity " << M << endl;
227  out << "-------------------------------------------------------" << endl;
228  out << " location | Gauss | S - B | Total | slope " << endl;
229  out << " | [Hz] | [Hz] | [Hz] | [Hz] " << endl;
230  out << "-------------------------------------------------------" << endl;
231 
232  JQuantile Q[4];
233 
234  for (int string = 1; string <= number_of_strings; ++string) {
235  for (int floor = number_of_floors; floor >= 1; --floor) {
236 
237  const int id = detector.getModule(JLocation(string,floor)).getID();
238 
239  out << " " << setw(3) << string << ' ' << setw(2) << floor << " ";
240 
241  TH1D* h1 = (H1.find(id) != H1.end() ? H1[id] : NULL);
242  TH1D* t1 = (T1.find(id) != T1.end() ? T1[id] : NULL);
243 
244  if (h1 != NULL) {
245 
246  TF1 f1("f1", "[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2]))/(TMath::Sqrt(2*TMath::Pi())*[2]) + [3]");
247 
248  f1.SetParameter(0, h1->GetMaximum());
249  f1.SetParameter(1, 0.0);
250  f1.SetParameter(2, h1->GetRMS() * 0.25);
251  f1.SetParameter(3, h1->GetMinimum());
252 
253  h1->Fit(&f1, option.c_str(), "same");
254 
255  out << " | " << FIXED(8,PRECISION) << f1.GetParameter(0);
256  out << " | " << FIXED(8,PRECISION) << (h1->GetSumOfWeights() - f1.GetParameter(3) * h1->GetNbinsX()) * V;
257  out << " | " << FIXED(8,PRECISION) << h1->GetSumOfWeights() * V;
258 
259  Q[0].put( f1.GetParameter(0));
260  Q[1].put((h1->GetSumOfWeights() - f1.GetParameter(3) * h1->GetNbinsX()) * V);
261  Q[2].put( h1->GetSumOfWeights() * V);
262  }
263 
264  if (t1 != NULL) {
265 
266  TF1 f1("f1", "[0]*exp(-[1]*x)");
267 
268  f1.SetParameter(0, t1->GetMaximum());
269  f1.SetParameter(1, 1.0 / t1->GetRMS());
270 
271  t1->Fit(&f1, option.c_str(), "same");
272 
273  out << " | " << FIXED(8,PRECISION) << f1.GetParameter(1);
274 
275  Q[3].put(f1.GetParameter(1));
276  }
277 
278  out << endl;
279  }
280 
281  out << endl;
282  }
283 
284  if (Q[0].getCount() != 0) {
285 
286  out << "-------------------------------------------------------" << endl;
287  out << setw(10) << left << " average";
288 
289  for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
290  out << " | " << FIXED(8,PRECISION) << Q[i].getMean();
291  }
292 
293  out << endl;
294  }
295 
296  out.close();
297  }
298 
299  if (outputFile != "") {
300 
301  TFile out(outputFile.c_str(), "RECREATE");
302 
303  H1.Write(out);
304  T1.Write(out);
305 
306  out.Close();
307  }
308 }
Utility class to parse command line options.
Definition: JParser.hh:1500
General exception.
Definition: JException.hh:23
Q(UTCMax_s-UTCMin_s)-livetime_s
do $JPP JMEstimator M
Definition: JMEstimator.sh:37
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
L0 match criterion.
Definition: JMatchL0.hh:27
Auxiliary data structure for running average, standard deviation and quantiles.
Definition: JQuantile.hh:43
do rm f tmp H1
#define STATUS(A)
Definition: JMessage.hh:63
Detector data structure.
Definition: JDetector.hh:81
Auxiliary class to select ROOT class based on class name.
Router for direct addressing of module data in detector data structure.
then for HISTOGRAM in h0 h1
Definition: JMatrixNZ.sh:71
Long64_t counter_type
Type definition for counter.
Auxiliary class for a type holder.
Definition: JType.hh:19
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
double getTimeOfRTS(const JDAQChronometer &chronometer)
Get time of last RTS in ns since start of run for a given chronometer.
string outputFile
Template definition for direct access of elements in ROOT TChain.
Definition: JTreeScanner.hh:91
long long int factorial(const long long int n)
Determine factorial.
Definition: JMathToolkit.hh:42
Auxiliary interface for direct access of elements in ROOT TChain.
double getMean() const
Get mean value.
Definition: JQuantile.hh:252
1-dimensional frame with time calibrated data from one optical module.
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys...
Definition: JManager.hh:43
Detector file.
Definition: JHead.hh:196
Logical location of module.
Definition: JLocation.hh:37
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
double getFrameTime()
Get frame time duration.
Definition: JDAQClock.hh:162
int getNumberOfFloors(const JDetector &detector)
Get number of floors.
void put(const double x, const double w=1.0)
Put value.
Definition: JQuantile.hh:133
Data time slice.
int debug
debug level
Definition: JSirene.cc:63
static const JStringCounter getNumberOfStrings
Function object to count unique strings.
Auxiliary class to select JTreeScanner based on ROOT class name.
#define FATAL(A)
Definition: JMessage.hh:67
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
then usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
int getCount(const T &hit)
Get hit count.
2-dimensional frame with time calibrated data from one optical module.
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:73
do set_variable DETECTOR_TXT $WORKDIR detector
then usage $script[input file[working directory[option]]] nWhere option can be N
Definition: JMuonPostfit.sh:36