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

Auxiliary program to determine time slewing from K40 data. More...

#include <string>
#include <iostream>
#include <iomanip>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TF1.h"
#include "TProfile.h"
#include "TProfile2D.h"
#include "km3net-dataformat/online/JDAQ.hh"
#include "JDAQ/JDAQTimesliceIO.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.hh"
#include "JMath/JMathToolkit.hh"
#include "JTrigger/JHitR0.hh"
#include "JTrigger/JHitToolkit.hh"
#include "JTools/JRange.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JSupport.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

Auxiliary program to determine time slewing from K40 data.

Author
mdejong

Definition in file JSlewingK40.cc.

Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 35 of file JSlewingK40.cc.

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 }
JDetector detector
Definition: JRunAnalyzer.hh:23
Utility class to parse command line options.
Definition: JParser.hh:1500
General exception.
Definition: JException.hh:23
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
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
Detector file.
Definition: JHead.hh:224
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:68
const JPMT & getPMT(const int index) const
Get PMT.
Definition: JModule.hh:173
#define FATAL(A)
Definition: JMessage.hh:67
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
General purpose class for object reading from a list of file names.
then usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
Template definition of hit toolkit.
Definition: JHitToolkit.hh:60
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:73