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
JHalibut.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 #include <fstream>
4 #include <iomanip>
5 #include <vector>
6 #include <map>
7 
8 #include "TROOT.h"
9 #include "TFile.h"
10 #include "TH1D.h"
11 #include "TF1.h"
12 
13 #include "JDAQ/JDAQTimesliceIO.hh"
15 #include "JDAQ/JDAQEvaluator.hh"
16 
17 #include "JDetector/JDetector.hh"
20 
21 #include "JTrigger/JHitR0.hh"
22 #include "JTrigger/JMatchL0.hh"
23 #include "JTrigger/JBuildL0.hh"
26 
27 #include "JROOT/JManager.hh"
29 
31 #include "JTools/JQuantile.hh"
32 #include "JTools/JRange.hh"
33 
35 #include "JSupport/JTreeScanner.hh"
37 #include "JSupport/JSupport.hh"
38 
40 #include "JMath/JMathToolkit.hh"
41 
42 #include "Jeep/JPrint.hh"
43 #include "Jeep/JParser.hh"
44 #include "Jeep/JMessage.hh"
45 
46 
47 /**
48  * \file
49  * Example program to determine N-fold coincidence rates.
50  * \author mdejong
51  */
52 int main(int argc, char **argv)
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
int main(int argc, char *argv[])
Definition: Main.cc:15
ROOT TTree parameter settings of various packages.
do $JPP JMEstimator M
Definition: JMEstimator.sh:37
Auxiliary methods for geometrical methods.
Basic data structure for L0 hit.
const JModule & getModule(const JObjectID &id) const
Get module parameters.
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
const int getIndex(const JObjectID &id) const
Get index of module.
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.
Dynamic ROOT object management.
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
Data structure for detector geometry and calibration.
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
I/O formatting auxiliaries.
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
Match operator for consecutive hits.
General purpose messaging.
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
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.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Auxiliary class to define a range between two values.
Utility class to parse command line options.
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