Jpp  test_elongated_shower_pde
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
examples/JSirene/JDomino.cc
Go to the documentation of this file.
1 
2 #include <string>
3 #include <iostream>
4 #include <iomanip>
5 #include <vector>
6 
7 #include "TROOT.h"
8 #include "TFile.h"
9 #include "TH1D.h"
10 
13 
14 #include "JAAnet/JHead.hh"
15 #include "JAAnet/JHeadToolkit.hh"
16 #include "JAAnet/JAAnetToolkit.hh"
17 
18 #include "JDetector/JDetector.hh"
20 #include "JDetector/JPMTRouter.hh"
21 
24 #include "JSupport/JSupport.hh"
25 
26 #include "JLang/JPredicate.hh"
27 #include "JTools/JRange.hh"
28 
29 #include "JROOT/JManager.hh"
30 #include "JROOT/JRootToolkit.hh"
31 #include "JGizmo/JGizmoToolkit.hh"
32 
33 #include "Jeep/JPrint.hh"
34 #include "Jeep/JParser.hh"
35 #include "Jeep/JMessage.hh"
36 
37 
38 /**
39  * \file
40  *
41  * Example program to verify Monte Carlo data.
42  * \author mdejong
43  */
44 int main(int argc, char **argv)
45 {
46  using namespace std;
47  using namespace JPP;
48 
49  typedef JRange<size_t> JRange_t;
50 
51  JMultipleFileScanner<Evt> inputFile;
52  JLimit_t& numberOfEvents = inputFile.getLimit();
53  string outputFile;
54  string detectorFile;
55  JRange_t M;
56  int debug;
57 
58  try {
59 
60  JParser<> zap("Example program to verify Monte Carlo data.");
61 
62  zap['f'] = make_field(inputFile);
63  zap['n'] = make_field(numberOfEvents) = JLimit::max();
64  zap['o'] = make_field(outputFile) = "domino.root";
65  zap['a'] = make_field(detectorFile);
66  zap['M'] = make_field(M) = JRange_t();
67  zap['d'] = make_field(debug) = 2;
68 
69  zap(argc, argv);
70  }
71  catch(const exception &error) {
72  FATAL(error.what() << endl);
73  }
74 
75  cout.tie(&cerr);
76 
78 
79  if (detectorFile != "") {
80 
81  try {
82  load(detectorFile, detector);
83  }
84  catch(const JException& error) {
85  FATAL(error);
86  }
87  }
88 
89  const JPMTRouter router(detector);
90 
91  const JHead header = getHeader(inputFile);
92 
93  JPosition3D center = get<JPosition3D>(header);
94 
95  NOTICE("Apply detector offset " << center << endl);
96 
97  detector -= center;
98 
100 
101  double x = -100.0;
102 
103  for ( ; x < -10.0; x += 5.0) { X.push_back(x); }
104  for ( ; x < +20.0; x += 1.0) { X.push_back(x); }
105  for ( ; x < +50.0; x += 2.0) { X.push_back(x); }
106  for ( ; x < +100.0; x += 5.0) { X.push_back(x); }
107  for ( ; x < +250.0; x += 10.0) { X.push_back(x); }
108  for ( ; x < +900.0; x += 50.0) { X.push_back(x); }
109 
110  JManager<int, TH1D> H1(new TH1D("h[%]", NULL, X.size() - 1, X.data()));
111 
112  size_t number_of_events = 0;
113 
114  while (inputFile.hasNext()) {
115 
116  STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
117 
118  const Evt* event = inputFile.next();
119 
120  if (!event->mc_hits.empty()) {
121  ++number_of_events;
122  }
123 
124  size_t number_of_muons = 0;
125 
126  for (vector<Trk>::const_iterator track = event->mc_trks.begin(); track != event->mc_trks.end(); ++track) {
127  if (is_muon(*track)) {
128  number_of_muons += 1;
129  }
130  }
131 
132  if (!M(number_of_muons)) {
133  continue;
134  }
135 
136  for (vector<Hit>::const_iterator hit = event->mc_hits.begin(); hit != event->mc_hits.end(); ++hit) {
137 
138  const double t1 = getTime(*hit);
139 
140  vector<Trk>::const_iterator track = find_if(event->mc_trks.begin(),
141  event->mc_trks.end(),
142  make_predicate(&Trk::id, hit->origin));
143 
144  if (track != event->mc_trks.end() && is_muon(*track)) {
145 
146  const JTrack3E muon = getTrack(*track);
147 
148  if (router.hasPMT(hit->pmt_id)) {
149 
150  const JPMT& pmt = router.getPMT(hit->pmt_id);
151  const double t0 = muon.getT(pmt);
152 
153  H1[pmt.getID()]->Fill(t1 - t0, hit->a);
154 
155  H1->Fill(t1 - t0, hit->a);
156  }
157  }
158  }
159  }
160  STATUS(endl);
161 
162  if (number_of_events != 0) {
163 
164  const double W = 1.0 / (double) number_of_events;
165 
166  convertToPDF(*H1, "WE", W);
167 
168  for (auto& i : H1) {
169  convertToPDF(*i.second, "WE", W);
170  }
171  }
172 
173  TFile out(outputFile.c_str(), "recreate");
174 
175  out << *H1 << H1;
176 
177  out.Write();
178  out.Close();
179 }
JDetector detector
Definition: JRunAnalyzer.hh:23
Router for direct addressing of PMT data in detector data structure.
Definition: JPMTRouter.hh:35
Utility class to parse command line options.
Definition: JParser.hh:1500
General exception.
Definition: JException.hh:23
JPredicate< JResult_t T::*, JComparison::eq > make_predicate(JResult_t T::*member, const JResult_t value)
Helper method to create predicate for data member.
Definition: JPredicate.hh:128
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
JTrack3E getTrack(const Trk &track)
Get track.
#define STATUS(A)
Definition: JMessage.hh:63
Detector data structure.
Definition: JDetector.hh:89
bool is_muon(const Trk &track)
Test whether given track is a (anti-)muon.
bool hasPMT(const JObjectID &id) const
Has PMT.
Definition: JPMTRouter.hh:116
const JPMT & getPMT(const JPMTAddress &address) const
Get PMT.
Definition: JPMTRouter.hh:92
Dynamic ROOT object management.
double getTime(const Hit &hit)
Get true time of hit.
string outputFile
3D track with energy.
Definition: JTrack3E.hh:30
Data structure for detector geometry and calibration.
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
Type definition of range.
Definition: JHead.hh:39
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:224
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
int getID() const
Get identifier.
Definition: JObjectID.hh:50
void convertToPDF(TH1 &h1, const std::string &option="NW", const double factor=1.0)
Convert 1D histogram to PDF.
Data structure for PMT geometry, calibration and status.
Definition: JPMT.hh:43
Definition of hit and track types and auxiliary methods for handling Monte Carlo data.
#define NOTICE(A)
Definition: JMessage.hh:64
then break fi done getCenter read X Y Z let X
void Write(TDirectory &out, const bool wm=false)
Write objects to file.
Definition: JManager.hh:295
Direct access to PMT in detector data structure.
int debug
debug level
Definition: JSirene.cc:68
double getT(const JVector3D &pos) const
Get arrival time of Cherenkov light at given position.
Definition: JTrack3D.hh:126
Range of values.
Definition: JRange.hh:38
General purpose messaging.
Monte Carlo run header.
Definition: JHead.hh:1164
#define FATAL(A)
Definition: JMessage.hh:67
Scanning of objects from multiple files according a format that follows from the extension of each fi...
int id
track identifier
Definition: Trk.hh:16
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
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:73
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition: Evt.hh:19