Jpp  master_rocky
the software that should make you happy
JOffset_extract.cc
Go to the documentation of this file.
1 \
2 #include <iostream>
3 #include <string>
4 #include <map>
5 
6 #include "TFile.h"
7 #include "TH1D.h"
8 #include "TH2D.h"
9 #include "TF1.h"
10 #include "TMath.h"
11 
12 #include "JROOT/JMinimizer.hh"
13 
14 #include "JDetector/JDetector.hh"
17 
18 #include "Jeep/JPrint.hh"
19 #include "Jeep/JParser.hh"
20 
21 
22 /**
23  * \file
24  *
25  * Auxiliary program to extract string-string correlation information from the output created with JMonitorL1dt
26  *
27  * \author dsamtleben
28  */
29 int main(int argc, char **argv)
30 {
31  using namespace std;
32  using namespace JPP;
33 
34  string inputFile;
35  string outputFile;
36  string detectorFile;
37  int neighbour, setmax;
38  int debug;
39 
40  try {
41 
42  JParser<> zap("Program to extract time offsets of DOM-DOM correlations");
43 
44  zap['f'] = make_field(inputFile, "input file") = "monitor.root";
45  zap['o'] = make_field(outputFile, "output file") = "corr_max.root";
46  zap['a'] = make_field(detectorFile, "detector file");
47  zap['n'] = make_field(neighbour, "neighbour level") = 5;
48  zap['m'] = make_field(setmax, "minimal entries") = 5;
49  zap['d'] = make_field(debug) = 0;
50 
51  if (zap.read(argc, argv) != 0) {
52  return 1;
53  }
54  }
55  catch(const exception &error) {
56  FATAL(error.what() << endl);
57  }
58 
59 
61 
62  try {
63  load(detectorFile, detector);
64  }
65  catch(const JException& error) {
66  FATAL(error);
67  }
68 
69 
70  const JStringRouter string(detector);
71 
73 
74  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
75  zmap[module->getString()][module->getFloor()] = module->getID();
76  }
77 
78  struct h1_t {
79  h1_t() :
80  p(NULL) // pointer must be initialised
81  {}
82 
83  TH1D* p;
84  };
85 
86  TF1 f1("f1", "[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2])) + [3]");
87 
88  const int number_of_strings = getNumberOfStrings(detector);
89 
90 // cout << "number of strings " << number_of_strings << endl;
91 
92  double maxarr[number_of_strings][number_of_strings][3];
93  int duarr[number_of_strings];
94 
96 
97  TFile* in = TFile::Open(inputFile.c_str(), "exist");
98 
99  if (in == NULL || !in->IsOpen()) {
100  FATAL("File: " << inputFile << " not opened." << endl);
101  }
102 
103 
104  for (const auto& string_1 : zmap) {
105  duarr[string.getIndex(string_1.first)]=string_1.first;
106 
107  for (const auto& floor_1 : string_1.second) {
108 
109  const int module_1 = floor_1.second; // parent module
110 
111  TH2D* h2 = (TH2D*) in->Get(MAKE_CSTRING(module_1 << ".2S"));
112 
113  if (h2 == NULL) {
114  continue;
115  }
116 
117  for (const auto& string_2 : zmap) {
118  duarr[string.getIndex(string_2.first)]=string_2.first;
119 
120  if (string_1.first != string_2.first) {
121 
122  TH1D* h1 = H1[string_1.first][string_2.first].p;
123 
124  if (h1 == NULL) {
125 
126  h1 = new TH1D(MAKE_CSTRING(string_1.first << "_" << string_2.first << "_" << neighbour << ".2T"), NULL,
127  h2->GetYaxis()->GetNbins(), h2->GetYaxis()->GetXmin(), h2->GetYaxis()->GetXmax());
128 
129  H1[string_1.first][string_2.first].p = h1; // book keeping
130  }
131 
132  for (const auto& floor_2 : string_2.second) {
133  if (floor_1.first > floor_2.first && (floor_1.first - floor_2.first) == neighbour) { // floor condition
134 
135  const int module_2 = floor_2.second; // daughter module
136 
137  TH1D* py = h2->ProjectionY("__py",
138  h2->GetXaxis()->FindBin(TString(Form("%i", module_2))),
139  h2->GetXaxis()->FindBin(TString(Form("%i", module_2))), "e");
140 
141  h1->Add(py); // add data
142 
143  delete py;
144  }
145  }
146 //------------
147  if (floor_1.first == 18){
148 
149  double mm=h1->GetXaxis()->GetBinCenter(h1->GetMaximumBin());
150 
151  f1.SetParameter(0, double(h1->GetMaximum()));
152  f1.SetParameter(1, double(mm));
153  f1.SetParameter(2, 50.);
154 
155  for (Int_t i = 0; i != f1.GetNpar(); ++i) {
156  f1.SetParError(i, 0.0);
157  }
158 
159  // fit
160 
161  if (h1->GetMaximum()>0){
162 // cout << "we are fitting " << string_1.first << " " << string_2.first << endl;
163  h1->Fit(&f1,"LQ", "same");
164 
165 // maxarr[ string.getIndex(string_1.first) ][ string.getIndex(string_2.first) ][0] = h1->GetXaxis()->GetBinCenter(h1->GetMaximumBin());
166 // maxarr[ string.getIndex(string_1.first) ][ string.getIndex(string_2.first) ][1] = h1->GetMaximum();
167  maxarr[ string.getIndex(string_1.first) ][ string.getIndex(string_2.first) ][0] = h1->GetFunction("f1")->GetParameter(1);
168  maxarr[ string.getIndex(string_1.first) ][ string.getIndex(string_2.first) ][1] = h1->GetFunction("f1")->GetParameter(0);
169  maxarr[ string.getIndex(string_1.first) ][ string.getIndex(string_2.first) ][2] = h1->GetFunction("f1")->GetParameter(2);
170  }
171  }
172  } // if string 1 != string 2
173  } // loop over string 2
174  } // loop over modules of string 1
175  } // loop over string 1
176 
177 
178  for (int i=0;i<number_of_strings;i++){
179  for (int j=i+1;j<number_of_strings;j++){
180  if (maxarr[i][j][1]>setmax && maxarr[j][i][1]>setmax && maxarr[i][j][2]>20. && maxarr[j][i][2]>20.){
181  cout << i << " " << j << " " << maxarr[i][j][0] << " " << maxarr[j][i][0] << " " << (maxarr[i][j][0]-maxarr[j][i][0])/2. << " " << maxarr[i][j][1] << " " << maxarr[j][i][1] << " " << " " << duarr[i] << " " << duarr[j] << " " << neighbour << endl;
182  }
183  }
184  }
185 
186 
187 
188  TFile out(outputFile.c_str(), "recreate");
189 
190  out.cd();
191 
192  for (auto& h1 : H1) {
193  for (auto& h2 : h1.second) {
194  h2.second.p->Write();
195  }
196  }
197 
198  out.Close();
199 }
string outputFile
Data structure for detector geometry and calibration.
#define FATAL(A)
Definition: JMessage.hh:67
int debug
debug level
Definition: JSirene.cc:69
int main(int argc, char **argv)
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:2142
I/O formatting auxiliaries.
#define MAKE_CSTRING(A)
Make C-string.
Definition: JPrint.hh:72
Direct access to string in detector data structure.
Detector data structure.
Definition: JDetector.hh:96
General exception.
Definition: JException.hh:24
Utility class to parse command line options.
Definition: JParser.hh:1698
int read(const int argc, const char *const argv[])
Parse the program's command line options.
Definition: JParser.hh:1992
const JPolynome f1(1.0, 2.0, 3.0)
Function.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
static const JStringCounter getNumberOfStrings
Function object to count unique strings.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
int j
Definition: JPolint.hh:792
Definition: JSTDTypes.hh:14
Detector file.
Definition: JHead.hh:227
Router for mapping of string identifier to index.