Jpp  18.0.0-rc.3
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JSlewingK40.cc
Go to the documentation of this file.
1 
2 #include <string>
3 #include <iostream>
4 #include <iomanip>
5 
6 #include "TROOT.h"
7 #include "TFile.h"
8 #include "TH1D.h"
9 #include "TH2D.h"
10 #include "TF1.h"
11 #include "TProfile.h"
12 #include "TProfile2D.h"
13 
15 #include "JDAQ/JDAQTimesliceIO.hh"
16 #include "JDetector/JDetector.hh"
19 #include "JMath/JMathToolkit.hh"
20 #include "JTrigger/JHitR0.hh"
21 #include "JTrigger/JHitToolkit.hh"
22 #include "JTools/JRange.hh"
24 #include "JSupport/JSupport.hh"
25 #include "Jeep/JParser.hh"
26 #include "Jeep/JMessage.hh"
27 
28 
29 /**
30  * \file
31  *
32  * Auxiliary program to determine time slewing from K40 data.
33  * \author mdejong
34  */
35 int main(int argc, char **argv)
36 {
37  using namespace std;
38  using namespace JPP;
39  using namespace KM3NETDAQ;
40 
42  JLimit_t& numberOfEvents = inputFile.getLimit();
43  string outputFile;
44  string detectorFile;
45  Double_t TMax_ns;
46  JRange<double> totVeto_ns;
47  JRange<int> multiplicity;
48  JRange<double> ct;
49  bool slewing;
50  int debug;
51 
52  try {
53 
54  JParser<> zap("Auxiliary program to determine time slewing from K40 data.");
55 
56  zap['f'] = make_field(inputFile);
57  zap['o'] = make_field(outputFile) = "monitor.root";
58  zap['a'] = make_field(detectorFile);
59  zap['n'] = make_field(numberOfEvents) = JLimit::max();
60  zap['T'] = make_field(TMax_ns) = 20.5; // [ns]
61  zap['t'] = make_field(totVeto_ns) = JRange<double>(0, 1e4); // time over threshold of reference hits [ns]
62  zap['M'] = make_field(multiplicity) = JRange<int>(2, 2);
63  zap['C'] = make_field(ct) = JRange<double>(0.0, 1.0); // cosine space angle between PMTs
64  zap['S'] = make_field(slewing);
65  zap['d'] = make_field(debug) = 1;
66 
67  zap(argc, argv);
68  }
69  catch(const exception &error) {
70  FATAL(error.what() << endl);
71  }
72 
73 
74  using namespace KM3NETDAQ;
75 
76  if (!multiplicity.is_valid()) { FATAL("Invalid multiplicity " << multiplicity << endl); }
77  if ( multiplicity.getLowerLimit() < 2) { FATAL("Invalid multiplicity " << multiplicity << endl); }
78  if (!totVeto_ns .is_valid()) { FATAL("Invalid time over threshold " << totVeto_ns << endl); }
79 
80  JHit::setSlewing(slewing);
81 
83 
84  try {
85  load(detectorFile, detector);
86  }
87  catch(const JException& error) {
88  FATAL(error);
89  }
90 
91  if (detector.empty()) {
92  FATAL("Empty detector." << endl);
93  }
94 
95  const JModuleRouter router(detector);
96 
97 
98  TFile out(outputFile.c_str(), "recreate");
99 
100  TH1D h0("h0", NULL, 41, -TMax_ns, +TMax_ns);
101  TH1D h1("h1", NULL, 50, -0.5, 49.5);
102  TProfile2D hx("hx", NULL, 50, -0.5, 49.5, 50, -0.5, 49.5, -2*TMax_ns, +2*TMax_ns);
103  TProfile h2("h2", NULL, 100, 0.5, 100.5, -2*TMax_ns, +2*TMax_ns);
104 
105  h0.Sumw2();
106  h1.Sumw2();
107 
108 
109  vector<JHitR0> buffer;
110 
111  while (inputFile.hasNext()) {
112 
113  STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
114 
115  const JDAQTimeslice* timeslice = inputFile.next();
116 
117  for (JDAQTimeslice::const_iterator super_frame = timeslice->begin(); super_frame != timeslice->end(); ++super_frame) {
118 
119  if (router.hasModule(super_frame->getModuleID()) && !super_frame->empty()) {
120 
121  const JModule& module = router.getModule(super_frame->getModuleID());
122 
123  // process data
124 
125  buffer.resize(super_frame->size());
126 
127  vector<JHitR0>::iterator out = buffer.begin();
128 
129  for (JDAQSuperFrame::const_iterator i = super_frame->begin(); i != super_frame->end(); ++i) {
130  *out = JHitR0(i->getPMT(), JHitToolkit<JHit>::getHit(*i, module.getPMT(i->getPMT())));
131  ++out;
132  }
133 
134  sort(buffer.begin(), buffer.end());
135 
136  for (vector<JHitR0>::iterator p = buffer.begin(); p != buffer.end(); ) {
137 
139 
140  while (++q != buffer.end() && q->getT() - p->getT() <= TMax_ns ) {}
141 
142  if (multiplicity(distance(p,q))) {
143 
144  random_shuffle(p,q);
145 
146  for (vector<JHitR0>::const_iterator __p = p; __p != q; ++__p) {
147  for (vector<JHitR0>::const_iterator __q = __p; ++__q != q; ) {
148 
149  const double dot = JMATH::getDot(module.getPMT(__p->getPMT()),
150  module.getPMT(__q->getPMT()));
151 
152  if (ct(dot)) {
153 
154  h0.Fill(__q->getT() - __p->getT());
155  h1.Fill(__q->getToT());
156  hx.Fill(__p->getToT(), __q->getToT(), __q->getT() - __p->getT());
157 
158  if (totVeto_ns(__p->getToT())) {
159  h2.Fill(__q->getToT(), __q->getT() - __p->getT());
160  }
161  }
162  }
163  }
164 
165  sort(p,q);
166  }
167 
168  p = q;
169  }
170  }
171  }
172  }
173  STATUS(endl);
174 
175 
176  // Fit function
177 
178  TF1 f1("f1", "[0]*exp([1]*sqrt(x) + [2]*x) + [3]");
179 
180  f1.SetParameter(0, h2.GetMaximum());
181  f1.SetParameter(1, -0.01);
182  f1.SetParameter(2, -0.05);
183  f1.SetParameter(3, h2.GetMinimum());
184 
185  h2.ProjectionX()->Fit(&f1, "", "same");
186 
187  for (int i = 0; i != f1.GetNpar(); ++i) {
188  cout << "\tstatic double p" << i << "() { return " << setw(9) << setprecision(5) << f1.GetParameter(i) << "; }" << endl;
189  }
190 
191 
192  out.Write();
193  out.Close();
194 }
Utility class to parse command line options.
Definition: JParser.hh:1514
General exception.
Definition: JException.hh:23
int main(int argc, char *argv[])
Definition: Main.cc:15
ROOT TTree parameter settings of various packages.
Auxiliary methods for geometrical methods.
Basic data structure for L0 hit.
const JModule & getModule(const JObjectID &id) const
Get module parameters.
Data structure for a composite optical module.
Definition: JModule.hh:68
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
#define STATUS(A)
Definition: JMessage.hh:63
bool is_valid(const json &js)
Check validity of JSon data.
Detector data structure.
Definition: JDetector.hh:89
Router for direct addressing of module data in detector data structure.
string outputFile
Data structure for detector geometry and calibration.
Tools for handling different hit types.
const JPolynome f1(1.0, 2.0, 3.0)
Function.
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
Detector file.
Definition: JHead.hh:226
Reduced data structure for L0 hit.
Definition: JHitR0.hh:25
Hit data structure.
Definition: JDAQHit.hh:34
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1989
Data time slice.
double getDot(const JFirst_t &first, const JSecond_t &second)
Get dot product of objects.
then $JPP_DIR examples JDetector JSlewing P slewing
Definition: JSlewing.sh:43
const JPMT & getPMT(const int index) const
Get PMT.
Definition: JModule.hh:173
General purpose messaging.
#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.
General purpose class for object reading from a list of file names.
Utility class to parse command line options.
bool hasModule(const JObjectID &id) const
Has module.
Template definition of hit toolkit.
Definition: JHitToolkit.hh:60
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:84
do set_variable DETECTOR_TXT $WORKDIR detector
KM3NeT DAQ constants, bit handling, etc.
int debug
debug level
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62