Jpp
JPulsar.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 #include <iomanip>
4 #include <map>
5 
6 #include "TROOT.h"
7 #include "TFile.h"
8 #include "TH1D.h"
9 #include "TF1.h"
10 #include "TFitResult.h"
11 
12 #include "JDAQ/JDAQEvent.hh"
13 #include "JDAQ/JDAQTimeslice.hh"
14 
15 #include "JDetector/JDetector.hh"
16 #include "JDetector/JDetectorToolkit.hh"
17 #include "JDetector/JModuleRouter.hh"
18 #include "JDetector/JPMTIdentifier.hh"
19 
20 #include "JTrigger/JHit.hh"
21 #include "JTrigger/JFrame.hh"
22 #include "JTrigger/JTimeslice.hh"
23 #include "JTrigger/JHitL0.hh"
24 #include "JTrigger/JBuildL0.hh"
25 
26 #include "JSupport/JMultipleFileScanner.hh"
27 #include "JSupport/JSupport.hh"
28 #include "JSupport/JMeta.hh"
29 
30 #include "JLang/JObjectMultiplexer.hh"
31 #include "JROOT/JROOTClassSelector.hh"
32 
33 #include "Jeep/JParser.hh"
34 #include "Jeep/JMessage.hh"
35 
36 
37 /**
38  * \file
39  *
40  * Application for dark room time calibration.
41  * \author mbouwhuis
42  */
43 int main(int argc, char **argv)
44 {
45  using namespace std;
46  using namespace JPP;
47  using namespace KM3NETDAQ;
48 
49  typedef multimap<int, int> JTDC_t;
50  typedef std::pair<JTDC_t::const_iterator, JTDC_t::const_iterator> range_type;
51 
52  JMultipleFileScanner<JDAQTimesliceTypes_t> inputFile;
53  JLimit_t& numberOfEvents = inputFile.getLimit();
54  string outputFile;
55  string detectorFile;
56  JTDC_t TDC;
57  double laserFrequency_Hz;
58  bool overwriteDetector;
59  JROOTClassSelector selector;
60  string option;
61  int debug;
62 
63  try {
64 
65  JParser<> zap("Application for dark room time calibration.");
66 
67  zap['f'] = make_field(inputFile, "input file (output from JCalibrateK40).");
68  zap['o'] = make_field(outputFile, "output file.") = "pulsar.root";
69  zap['a'] = make_field(detectorFile, "detector file.");
70  zap['!'] = make_field(TDC,
71  "Set reference PMTs, e.g."
72  "\n-! \"808969848 0 808982077 23\" sets PMT 0 of module 808969848 and PMT 23 of module 808982077 as references.");
73  zap['n'] = make_field(numberOfEvents) = JLimit::max();
74  zap['l'] = make_field(laserFrequency_Hz, "laser frequency [Hz]") = 10000;
75  zap['A'] = make_field(overwriteDetector, "overwrite detector file provided through '-a' with correct time offsets.");
76  zap['O'] = make_field(option, "ROOT fit option, see TH1::Fit.") = "L";
77  zap['C'] = make_field(selector, "timeslice selector, e.g. JDAQTimesliceL1.") = getROOTClassSelection<JDAQTimesliceTypes_t>();
78  zap['d'] = make_field(debug, "debug.") = 1;
79 
80  zap(argc, argv);
81  }
82  catch(const exception& error) {
83  FATAL(error.what() << endl);
84  }
85 
86 
87  for (JTDC_t::const_iterator tdc = TDC.begin(); tdc != TDC.end(); ++tdc) {
88  if (tdc->second < 0 || tdc->second >= NUMBER_OF_PMTS) {
89  FATAL("Illegal TDC (= readout channel) identifier: " << tdc->second << " {0, .., " << NUMBER_OF_PMTS - 1 << "}" << endl);
90  }
91  }
92 
93  if (laserFrequency_Hz < 0.0) {
94  FATAL("Invalid laser frequency " << laserFrequency_Hz << endl);
95  }
96 
97 
98  const double laserPeriod_ns = 1.0e9 / laserFrequency_Hz;
99 
100  cout.tie(&cerr);
101 
102  JDetector detector;
103 
104  try {
105  load(detectorFile, detector);
106  }
107  catch(const JException& error) {
108  FATAL(error);
109  }
110 
111  const JModuleRouter moduleRouter(detector);
112 
113 
114  // ROOT I/O
115 
116  TFile out(outputFile.c_str(), "recreate");
117 
118  putObject(&out, JMeta(argc, argv));
119 
120  typedef map<JPMTIdentifier, TH1D*> map_type;
121 
122  map_type zmap;
123 
124  const int nx = 2 * (int) laserPeriod_ns;
125  const double xmin = -0.5 * laserPeriod_ns;
126  const double xmax = 0.5 * laserPeriod_ns;
127 
128 
129  for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
130 
131  const range_type range = TDC.equal_range(module->getID());
132 
133  for (JTDC_t::const_iterator i = range.first; i != range.second; ++i) {
134 
135  NOTICE("Reference PMT at"
136  << " string " << setw(3) << module->getString()
137  << " floor " << setw(2) << module->getFloor()
138  << " module " << setw(8) << module->getID()
139  << " channel " << setw(2) << i->second << endl);
140 
141  ostringstream os;
142 
143  os << "s" << module->getString()
144  << "_f" << module->getFloor()
145  << "_m" << module->getID()
146  << "_c" << setfill('0') << setw(2) << i->second;
147 
148  zmap.insert(make_pair(JPMTIdentifier(module->getID(),i->second), new TH1D(os.str().c_str(), NULL, nx, xmin, xmax)));
149  }
150  }
151 
152 
153  typedef vector<JHitL0> JDataL0_t;
154 
155  const JBuildL0<JHitL0> buildL0;
156 
157 
158  JObjectMultiplexer<JDAQTimesliceTypes_t, JDAQTimeslice> in(inputFile, selector);
159 
160  for (counter_type counter = 0; in.hasNext() && counter != inputFile.getLimit(); ++counter) {
161 
162  STATUS("event: " << setw(10) << counter << '\r'); DEBUG(endl);
163 
164  const JDAQTimeslice* timeslice = in.next();
165 
166  for (JDAQTimeslice::const_iterator frame = timeslice->begin(); frame != timeslice->end(); ++frame) {
167 
168  const int id = frame->getModuleID();
169 
170  DEBUG("Module id =" << id << endl);
171 
172  double t0 = getTimeOfFrame(frame->getFrameIndex());
173 
174  while (t0 > xmax) {
175  t0 -= laserPeriod_ns;
176  }
177 
178  JDataL0_t dataL0;
179 
180  buildL0(*frame, moduleRouter.getModule(id), back_inserter(dataL0));
181 
182  for (JDataL0_t::const_iterator hit = dataL0.begin(); hit != dataL0.end(); ++hit) {
183 
184  map_type::const_iterator p = zmap.find(JPMTIdentifier(id,hit->getPMTAddress()));
185 
186  if (p != zmap.end()) {
187 
188  double t1 = t0 + hit->getT();
189 
190  while (t1 > xmax) {
191  t1 -= laserPeriod_ns;
192  }
193  while (t1 < xmin) {
194  t1 += laserPeriod_ns;
195  }
196 
197  p->second->Fill(t1);
198  }
199  }
200  }
201  }
202 
203  STATUS(endl);
204 
205 
206  map<JPMTIdentifier, double> t0; // fitted time offsets
207 
208 
209  TF1 f1("f1", "[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2])) + [3]");
210 
211  for (map_type::iterator it = zmap.begin(); it != zmap.end(); ++it) {
212 
213  JPMTIdentifier pmt = it->first;
214  TH1D* h1 = it->second;
215 
216  if (h1->GetEntries() == 0) {
217  WARNING("Histogram " << h1->GetName() << " empty" << endl);
218  continue;
219  }
220 
221  STATUS("--- PMT = " << pmt << "; histogram " << h1->GetName() << endl);
222 
223  // start values
224 
225  Double_t ymax = 0.0;
226  Double_t x0 = 0.0; // [ns]
227  Double_t sigma = 2.0; // [ns]
228  Double_t ymin = 0.0;
229 
230  for (int i = 1; i != h1->GetNbinsX(); ++i) {
231 
232  const Double_t x = h1->GetBinCenter (i);
233  const Double_t y = h1->GetBinContent(i);
234 
235  if (y > ymax) {
236  ymax = y;
237  x0 = x;
238  }
239  }
240 
241  f1.SetParameter(0, ymax);
242  f1.SetParameter(1, x0);
243  f1.SetParameter(2, sigma);
244  f1.SetParameter(3, ymin);
245 
246 
247  // fit
248 
249  TFitResultPtr result = h1->Fit(&f1, option.c_str(), "same", x0 - 3 * sigma, x0 + 3 * sigma);
250 
251 
252  if (debug >= notice_t || !result->IsValid()) {
253  cout << "Histogram " << h1->GetName() << " fit " << (result->IsValid() ? "" : "failed") << endl;
254  cout << "\tx0 = " << FIXED(12,3) << x0 << endl;
255  cout << "\tt0 = " << FIXED(12,3) << f1.GetParameter(1) << endl;
256  }
257 
258  t0[pmt] = f1.GetParameter(1);
259  }
260 
261 
262  if (overwriteDetector) {
263 
264  NOTICE("Store calibration data on file " << detectorFile << endl);
265 
266  detector.comment.add(JMeta(argc, argv));
267 
268  for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
269 
270  const range_type range = TDC.equal_range(module->getID());
271 
272  if (range.first != range.second) {
273 
274  const double t1 = t0[JPMTIdentifier(range.first->first,
275  range.first->second)]; // t0 of the first reference PMT in this optical module
276 
277  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
278 
279  map<JPMTIdentifier, double>::const_iterator p = t0.find(JPMTIdentifier(module->getID(),pmt));
280 
281  if (p != t0.end())
282  module->getPMT(pmt).subT0(p->second); // offset with t0 of this reference PMT
283  else
284  module->getPMT(pmt).subT0(t1); // offset with t0 of first reference PMT in same optical module
285  }
286  }
287  }
288 
289  store(detectorFile, detector);
290  }
291 
292  out.Write();
293  out.Close();
294 }
JMeta.hh
FIXED
Auxiliary data structure for floating point format specification.
Definition: JPrint.hh:481
JMessage.hh
JEEP::notice_t
notice
Definition: JMessage.hh:32
JTimeslice.hh
JDETECTOR::load
void load(const JString &file_name, JDetector &detector)
Load detector from input file.
Definition: JDetectorToolkit.hh:456
std::vector< JHitL0 >
JROOTClassSelector.hh
JSUPPORT::JLimit_t
JLimit JLimit_t
Type definition of limit.
Definition: JLimit.hh:215
KM3NETDAQ::NUMBER_OF_PMTS
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
Definition: JDAQ.hh:26
KM3NETDAQ::JDAQTimeslice
Data time slice.
Definition: JDAQTimeslice.hh:36
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1493
NOTICE
#define NOTICE(A)
Definition: JMessage.hh:64
JDAQTimeslice.hh
JROOT::putObject
bool putObject(TDirectory *dir, const T &object)
Write object to ROOT directory.
Definition: JRootFileWriter.hh:38
JSupport.hh
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JAAnetToolkit.hh:37
JBuildL0.hh
WARNING
#define WARNING(A)
Definition: JMessage.hh:65
debug
int debug
debug level
Definition: JSirene.cc:59
JTOOLS::result
return result
Definition: JPolint.hh:695
JPMTIdentifier.hh
JSUPPORT::JLimit::getLimit
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:73
JHit.hh
JModuleRouter.hh
JObjectMultiplexer.hh
JMultipleFileScanner.hh
std::map
Definition: JSTDTypes.hh:16
std::pair
Definition: JSTDTypes.hh:15
STATUS
#define STATUS(A)
Definition: JMessage.hh:63
JParser.hh
JDetectorToolkit.hh
JDETECTOR::store
void store(const JString &file_name, const JDetector &detector)
Store detector to output file.
Definition: JDetectorToolkit.hh:532
std::multimap
Definition: JSTDTypes.hh:17
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1954
JSUPPORT::JMeta
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:71
DEBUG
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
std
Definition: jaanetDictionary.h:36
KM3NETDAQ
KM3NeT DAQ data structures and auxiliaries.
Definition: DataQueue.cc:39
JDAQEvent.hh
JDetector.hh
JHitL0.hh
main
int main(int argc, char **argv)
Definition: JPulsar.cc:43
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
JROOT::counter_type
Long64_t counter_type
Type definition for counter.
Definition: JCounter.hh:24
KM3NETDAQ::getTimeOfFrame
double getTimeOfFrame(const int frame_index)
Get start time of frame in ns since start of run for a given frame index.
Definition: JDAQClock.hh:185
JFrame.hh
Generated by   doxygen 1.8.16