Jpp  17.3.0
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JDepthDependence.cc
Go to the documentation of this file.
1 #include <iomanip>
2 #include <iostream>
3 #include <sstream>
4 #include <string>
5 
6 #include "TFile.h"
7 #include "TH1D.h"
8 #include "TMath.h"
9 #include "TROOT.h"
10 
11 #include "TApplication.h"
12 #include "TLegend.h"
13 #include "TGraphErrors.h"
14 
16 #include "JDetector/JDetector.hh"
19 #include "Jeep/JParser.hh"
20 #include "Jeep/JMessage.hh"
21 #include "JLang/JString.hh"
23 #include "JSupport/JSupport.hh"
24 
25 
26 namespace {
27 
28  using namespace JPP;
29  using namespace std;
30 
31  double getCount (TH1D* hptr, int muon_threshold) {
32  double count = 0;
33  for (int i = hptr->GetXaxis()->FindBin(muon_threshold); i <= hptr->GetNbinsX(); i++) {
34  count += hptr->GetBinContent(i);
35  }
36  return count;
37  }
38 
39  double getLiveTime (TH1D* hptr, const JModule& module) {
40  string label = getLabel(module);
41  double livetime = hptr->GetBinContent(hptr->GetXaxis()->FindBin(label.c_str()));
42  return livetime;
43  }
44 
45 }
46 
47 /**
48  * \file
49  *
50  * Auxiliary program to find depth dependence of multiplicity rates
51  * \author mlincetto
52  */
53 int main(int argc, char **argv)
54 {
55  using namespace std;
56  using namespace JPP;
57 
58  //-----------------------------------------------------------
59  // parameter interface
60  //-----------------------------------------------------------
61 
62 
63  string inputFile;
64  string outputFile;
65  string detectorFile;
66  int line;
67  int debug;
68  int minMultiplicity;
69  int minLiveTime_s;
70  // string summaryFile;
71 
72  try {
73 
74  JParser<> zap("Auxiliary program to find depth dependence of multiplicity rates");
75 
76  zap['f'] = make_field(inputFile, "JMM input file");
77  zap['o'] = make_field(outputFile, "output file") = "depthdependence.root";
78  // zap['s'] = make_field(summaryFile, "summary file") = "depthdependence.txt";
79  zap['a'] = make_field(detectorFile);
80  zap['d'] = make_field(debug) = 1;
81  zap['l'] = make_field(line) = 1;
82  zap['M'] = make_field(minMultiplicity, "Minimum multiplicity") = 8;
83  zap['T'] = make_field(minLiveTime_s, "Minimum DOM livetime [s] to be eligible for plotting") = 3600;
84 
85  zap(argc, argv);
86  }
87 
88  catch(const exception &error) {
89  FATAL(error.what() << endl);
90  }
91 
92 
93  //----------------------
94  // loading detector file
95  //----------------------
96 
98 
99  try {
100  load(detectorFile, detector);
101  }
102  catch(const JException& error) {
103  FATAL(error);
104  }
105 
106  if (detector.empty())
107  FATAL("Empty detector." << endl);
108 
109  const JModuleRouter router(detector);
110 
111  double utm_z = detector.getUTMZ();
112 
113  NOTICE("Detector base UTM z [m]: " << utm_z << endl);
114 
115  //-----------------------------------------------------------
116  // load input files
117  //-----------------------------------------------------------
118 
119 
120  DEBUG("Loading input file " << inputFile << endl);
121 
122  TFile* in = TFile::Open(inputFile.c_str(), "exist");
123 
124  if (in == NULL || !in->IsOpen()) {
125  FATAL("File: " << inputFile << " not opened." << endl);
126  }
127 
128  //-----------------------------------------------------------
129  // recover rates from histograms
130  //-----------------------------------------------------------
131 
132 
133  DEBUG("Loading livetime histogram from " << inputFile << endl);
134 
135  TH1D* liveTime = (TH1D*)in->Get("LiveTime");
136 
137  if (liveTime == NULL) {
138  FATAL("Missing live time histogram.");
139  }
140 
142  vector<double> rate_val;
143  vector<double> rate_err;
144 
145  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
146 
147  if (module->getString() != line) { continue; }
148 
149  STATUS("Loading data histogram for from input file: " << getLabel(*module) << "\t ID = " << module->getID() << endl);
150 
151  TH1D* data_histogram = (TH1D*)in->Get(TString(getLabel(*module)) + "_P");
152 
153  if (data_histogram != NULL) {
154 
155  double data_count = getCount(data_histogram, minMultiplicity);
156  double data_livetime = getLiveTime(liveTime, *module);
157 
158  double module_depth = - utm_z - module->getZ();
159 
160  if (data_livetime > minLiveTime_s) {
161 
162  double val = data_count / data_livetime;
163  double err = sqrt(data_count) / data_livetime;
164 
165  NOTICE(module_depth << " " << val << endl);
166 
167  rate_val.push_back(val);
168  rate_err.push_back(err);
169  depth.push_back( module_depth);
170 
171  }
172 
173  }
174 
175  }
176 
177 
178  TGraph* gr_data = new TGraphErrors(depth.size(), &depth[0], &rate_val[0], 0, &rate_err[0]);
179 
180  gr_data->SetTitle(TString("KM3NeT Preliminary; Depth [m]; Inclusive ") + Form("%d", minMultiplicity) + TString("-fold coincidence rate [Hz]"));
181 
182  TFile out(outputFile.c_str(), "recreate");
183 
184  gr_data->Write();
185 
186  out.Close();
187 
188  return 0;
189 
190 }
Utility class to parse command line options.
Definition: JParser.hh:1517
General exception.
Definition: JException.hh:23
int main(int argc, char *argv[])
Definition: Main.cc:15
ROOT TTree parameter settings of various packages.
Data structure for a composite optical module.
Definition: JModule.hh:68
Depth.
Definition: JHead.hh:967
std::string getLabel(const JLocation &location)
Get module label for monitoring and other applications.
Definition: JLocation.hh:246
#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.
string outputFile
Data structure for detector geometry and calibration.
Detector file.
Definition: JHead.hh:226
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1993
int getID() const
Get identifier.
Definition: JObjectID.hh:50
#define NOTICE(A)
Definition: JMessage.hh:64
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.
int getString() const
Get string number.
Definition: JLocation.hh:134
Normalisation of MUPAGE events.
Definition: JHead.hh:822
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Utility class to parse command line options.
int getCount(const T &hit)
Get hit count.
do set_variable DETECTOR_TXT $WORKDIR detector
KM3NeT DAQ constants, bit handling, etc.
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable INPUT_FILE $argv[2] eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY JAcoustics sh $DETECTOR_ID source JAcousticsToolkit sh CHECK_EXIT_CODE typeset A EMITTERS get_tripods $WORKDIR tripod txt EMITTERS get_transmitters $WORKDIR transmitter txt EMITTERS for EMITTER in
Definition: JCanberra.sh:46
double getZ() const
Get z position.
Definition: JVector3D.hh:115
int debug
debug level
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62