Jpp  15.0.1-rc.2-highQE
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  cout.tie(&cerr);
74 
75  using namespace KM3NETDAQ;
76 
77  if (!multiplicity.is_valid()) { FATAL("Invalid multiplicity " << multiplicity << endl); }
78  if ( multiplicity.getLowerLimit() < 2) { FATAL("Invalid multiplicity " << multiplicity << endl); }
79  if (!totVeto_ns .is_valid()) { FATAL("Invalid time over threshold " << totVeto_ns << endl); }
80 
81  JHit::setSlewing(slewing);
82 
84 
85  try {
86  load(detectorFile, detector);
87  }
88  catch(const JException& error) {
89  FATAL(error);
90  }
91 
92  if (detector.empty()) {
93  FATAL("Empty detector." << endl);
94  }
95 
96  const JModuleRouter router(detector);
97 
98 
99  TFile out(outputFile.c_str(), "recreate");
100 
101  TH1D h0("h0", NULL, 41, -TMax_ns, +TMax_ns);
102  TH1D h1("h1", NULL, 50, -0.5, 49.5);
103  TProfile2D hx("hx", NULL, 50, -0.5, 49.5, 50, -0.5, 49.5, -2*TMax_ns, +2*TMax_ns);
104  TProfile h2("h2", NULL, 100, 0.5, 100.5, -2*TMax_ns, +2*TMax_ns);
105 
106  h0.Sumw2();
107  h1.Sumw2();
108 
109 
110  vector<JHitR0> buffer;
111 
112  while (inputFile.hasNext()) {
113 
114  STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
115 
116  const JDAQTimeslice* timeslice = inputFile.next();
117 
118  for (JDAQTimeslice::const_iterator super_frame = timeslice->begin(); super_frame != timeslice->end(); ++super_frame) {
119 
120  if (router.hasModule(super_frame->getModuleID()) && !super_frame->empty()) {
121 
122  const JModule& module = router.getModule(super_frame->getModuleID());
123 
124  // process data
125 
126  buffer.resize(super_frame->size());
127 
128  vector<JHitR0>::iterator out = buffer.begin();
129 
130  for (JDAQSuperFrame::const_iterator i = super_frame->begin(); i != super_frame->end(); ++i) {
131  *out = JHitR0(i->getPMT(), JHitToolkit<JHit>::getHit(*i, module.getPMT(i->getPMT())));
132  ++out;
133  }
134 
135  sort(buffer.begin(), buffer.end());
136 
137  for (vector<JHitR0>::iterator p = buffer.begin(); p != buffer.end(); ) {
138 
140 
141  while (++q != buffer.end() && q->getT() - p->getT() <= TMax_ns ) {}
142 
143  if (multiplicity(distance(p,q))) {
144 
145  random_shuffle(p,q);
146 
147  for (vector<JHitR0>::const_iterator __p = p; __p != q; ++__p) {
148  for (vector<JHitR0>::const_iterator __q = __p; ++__q != q; ) {
149 
150  const double dot = JMATH::getDot(module.getPMT(__p->getPMT()),
151  module.getPMT(__q->getPMT()));
152 
153  if (ct(dot)) {
154 
155  h0.Fill(__q->getT() - __p->getT());
156  h1.Fill(__q->getToT());
157  hx.Fill(__p->getToT(), __q->getToT(), __q->getT() - __p->getT());
158 
159  if (totVeto_ns(__p->getToT())) {
160  h2.Fill(__q->getToT(), __q->getT() - __p->getT());
161  }
162  }
163  }
164  }
165 
166  sort(p,q);
167  }
168 
169  p = q;
170  }
171  }
172  }
173  }
174  STATUS(endl);
175 
176 
177  // Fit function
178 
179  TF1 f1("f1", "[0]*exp([1]*sqrt(x) + [2]*x) + [3]");
180 
181  f1.SetParameter(0, h2.GetMaximum());
182  f1.SetParameter(1, -0.01);
183  f1.SetParameter(2, -0.05);
184  f1.SetParameter(3, h2.GetMinimum());
185 
186  h2.ProjectionX()->Fit(&f1, "", "same");
187 
188  for (int i = 0; i != f1.GetNpar(); ++i) {
189  cout << "\tstatic double p" << i << "() { return " << setw(9) << setprecision(5) << f1.GetParameter(i) << "; }" << endl;
190  }
191 
192 
193  out.Write();
194  out.Close();
195 }
Utility class to parse command line options.
Definition: JParser.hh:1500
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
Detector data structure.
Definition: JDetector.hh:89
Router for direct addressing of module data in detector data structure.
then for HISTOGRAM in h0 h1
Definition: JMatrixNZ.sh:71
string outputFile
Data structure for detector geometry and calibration.
Tools for handling different hit types.
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
Detector file.
Definition: JHead.hh:196
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:1961
bool is_valid(const json &js)
Check validity of JSon data.
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
int debug
debug level
Definition: JSirene.cc:63
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.
then usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
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:73
do set_variable DETECTOR_TXT $WORKDIR detector
KM3NeT DAQ constants, bit handling, etc.