Jpp
 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 
15 #include "JDAQ/JDAQ.hh"
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"
25 
26 using namespace JPP;
27 
28 double getCount (TH1D* hptr, int muon_threshold) {
29  double count = 0;
30  for (int i = hptr->GetXaxis()->FindBin(muon_threshold); i <= hptr->GetNbinsX(); i++) {
31  count += hptr->GetBinContent(i);
32  }
33  return count;
34 }
35 
36 double getLiveTime (TH1D* hptr, const JModule& module) {
37  string label = getModuleLabel(module);
38  double livetime = hptr->GetBinContent(hptr->GetXaxis()->FindBin(label.c_str()));
39  // trick for retrocompatibility with old JMM
40  if (livetime == 0) {
41  std::replace(label.begin(), label.end() - 1, '0', ' ');
42  livetime = hptr->GetBinContent(hptr->GetXaxis()->FindBin(label.c_str()));
43  }
44  return livetime;
45 }
46 
47 /**
48  * \file
49  *
50  * Auxiliary program to analyse coincidence multiplicity;
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  const int nFiles = 2;
63  vector<string> inputFile(nFiles,"");
64  string outputFile;
65  string detectorFile;
66  int line;
67  int debug;
68  int muonThreshold;
69  int minLiveTime_s;
70  bool dataLiveTime;
71  bool interactive;
72 
73  try {
74 
75  JParser<> zap("Example program to plot the muon depth-dependence relationship from JMonitorMultiplicity output");
76 
77  zap['f'] = make_field(inputFile[0], "JMM data input file");
78  zap['g'] = make_field(inputFile[1], "JMM MC input file") = "";
79  zap['o'] = make_field(outputFile, "Name for the .pdf and .root output files") = "depthdependence";
80  zap['a'] = make_field(detectorFile);
81  zap['d'] = make_field(debug) = 1;
82  zap['l'] = make_field(line) = 1;
83  zap['M'] = make_field(muonThreshold, "Minimum multiplicity") = 8;
84  zap['T'] = make_field(minLiveTime_s, "Minimum DOM livetime [s] to be eligible for plotting") = 3600;
85  zap['I'] = make_field(interactive, "Launch TApplication for interactive display.");
86  zap['L'] = make_field(dataLiveTime);
87 
88  zap(argc, argv);
89  }
90 
91  catch(const exception &error) {
92  FATAL(error.what() << endl);
93  }
94 
95  cout.tie(&cerr);
96 
97 
98  //----------------------
99  // loading detector file
100  //----------------------
101 
103 
104  try {
105  load(detectorFile, detector);
106  }
107  catch(const JException& error) {
108  FATAL(error);
109  }
110 
111  if (detector.empty())
112  FATAL("Empty detector." << endl);
113 
114  const JModuleRouter router(detector);
115 
116  double utm_z = detector.getUTMZ();
117 
118  NOTICE("Detector base UTM z [m]: " << utm_z << endl);
119 
120  //-----------------------------------------------------------
121  // load input files
122  //-----------------------------------------------------------
123 
124  vector<TFile*> in;
125  TFile* fptr;
126 
127  for (unsigned i = 0; i < inputFile.size() && inputFile[i] != ""; i++) {
128 
129  DEBUG("Loading input file " << inputFile[i] << endl);
130 
131  fptr = TFile::Open(inputFile[i].c_str(), "exist");
132 
133  if (fptr == NULL || !fptr->IsOpen()) {
134  FATAL("File: " << inputFile[i] << " not opened." << endl);
135  }
136 
137  in.push_back(fptr);
138  }
139 
140  //-----------------------------------------------------------
141  // recover rates from histograms
142  //-----------------------------------------------------------
143 
144  vector<TH1D*> LiveTime;
145 
146  for (unsigned i = 0; i < in.size(); i++) {
147  DEBUG("Loading livetime histogram from " << inputFile[i] << endl);
148  TH1D* hptr = (TH1D*)in[i]->Get("LiveTime");
149 
150  if (hptr == NULL) {
151  FATAL("Missing live time histogram.");
152  }
153 
154  LiveTime.push_back(hptr);
155  }
156 
157  vector<double> depth;
158  vector<double> data_rate;
159  vector<double> data_err;
160  vector<double> mc_rate;
161  vector<double> mc_err;
162 
163  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
164 
165  if (module->getString() != line) { continue; }
166 
167  STATUS("Loading data histogram for from input file: " << getModuleLabel(*module) << "\t ID = " << module->getID() << endl);
168 
169  TH1D* data_histogram = (TH1D*)in[0]->Get(TString(getModuleLabel(*module)) + "_P");
170 
171  if (data_histogram != NULL) {
172 
173  double data_count = getCount(data_histogram, muonThreshold);
174  double data_livetime = getLiveTime(LiveTime[0], *module);
175 
176  if (data_livetime > minLiveTime_s) {
177 
178  double rate = data_count / data_livetime;
179  double err = sqrt(data_count) / data_livetime;
180 
181  NOTICE(getModuleLabel(*module) << ": z_rel [m] = " << module->getZ() << ", count [#] = " << data_count << ", livetime [s] = " << data_livetime << ", rate [Hz] = " << rate << endl);
182 
183  data_rate.push_back(rate);
184  data_err.push_back(err);
185  depth.push_back((-utm_z) - module->getZ());
186 
187  if (in.size() > 1) {
188 
189  LiveTime[1] = (TH1D*)in[1]->Get("LiveTime");
190  TH1D* mc_histogram = (TH1D*)in[1]->Get(TString(getModuleLabel(*module)) + "_P");
191  double mc_count = getCount(mc_histogram, muonThreshold);
192  double mc_livetime = (dataLiveTime ? data_livetime : getLiveTime(LiveTime[1], *module));
193 
194  double rate = mc_count / mc_livetime;
195  double err = sqrt(mc_count) / mc_livetime;
196 
197  NOTICE("\t MC livetime [s] = " << mc_livetime << ", MC rate [Hz] = " << rate << endl);
198 
199  mc_rate.push_back(rate);
200  mc_err.push_back(err);
201  }
202  }
203  }
204 
205  }
206 
207  //-----------------------------------------------------------
208  // time to graph!
209  //-----------------------------------------------------------
210 
211  TApplication* theApp = NULL;
212 
213  if (interactive) {
214  theApp = new TApplication("Plotter", &argc, argv);
215  }
216 
217  TCanvas *c1 = new TCanvas("c1","KM3NeT preliminary");
218 
219  TGraph* gr_data = new TGraphErrors(data_rate.size(), &depth[0], &data_rate[0], 0, &data_err[0]);
220  gr_data->SetMarkerColor(kRed);
221  gr_data->SetLineColor(kRed);
222  TGraph* gr_mc;
223  TMultiGraph* mgr = new TMultiGraph();
224  mgr->Add(gr_data);
225 
226  TLegend* leg = new TLegend(0.7, 0.7, 0.9, 0.9);
227  leg->AddEntry(gr_data, "Data", "lp");
228 
229 
230  if (in.size() > 1) {
231  DEBUG("Loading MC plots" << endl);
232  gr_mc = new TGraphErrors(mc_rate.size(), &depth[0], &mc_rate[0], 0, &mc_err[0]);
233  mgr->Add(gr_mc);
234  DEBUG("MC plot loaded" << endl);
235  leg->AddEntry(gr_mc, "MC", "lp");
236 
237  }
238 
239  mgr->SetTitle(TString("KM3NeT Preliminary; Depth [m]; Inclusive ") + Form("%d", muonThreshold) + TString("-fold coincidence rate [Hz]"));
240 
241  mgr->Draw("AL");
242  mgr->GetYaxis()->SetTitleOffset(1.5);
243 
244  leg->SetFillStyle(0);
245  leg->Draw();
246 
247  c1->Print(TString(outputFile.c_str()) + ".pdf");
248 
249  if (interactive) {
250  theApp->Run();
251  }
252 
253  //-----------------------------------------------------------
254  // writing output
255  //-----------------------------------------------------------
256 
257  // writing of ROOT file should be included
258 
259  return 0;
260 
261 }
Utility class to parse command line options.
Definition: JParser.hh:1410
General exception.
Definition: JException.hh:40
Data structure for a composite optical module.
Definition: JModule.hh:47
std::string replace(const std::string &input, const std::string &target, const std::string &replacement)
Replace tokens in string.
double getCount(TH1D *hptr, int muon_threshold)
#define STATUS(A)
Definition: JMessage.hh:61
Detector data structure.
Definition: JDetector.hh:77
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:126
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1836
#define NOTICE(A)
Definition: JMessage.hh:62
void load(const JString &file_name, JDetector &detector)
Load detector from input file.
int debug
debug level
Definition: JSirene.cc:59
General purpose messaging.
#define FATAL(A)
Definition: JMessage.hh:65
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.
std::string getModuleLabel(const JModuleLocation &location)
Get module label (DU-floor) for JMonitor applications.
TCanvas * c1
Global variables to handle mouse events.
Normalisation of MUPAGE events.
Definition: JHead.hh:537
Utility class to parse command line options.
ROOT TTree parameter settings.
KM3NeT DAQ constants, bit handling, etc.
double getLiveTime(TH1D *hptr, const JModule &module)
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:60
int main(int argc, char *argv[])
Definition: Main.cpp:15