Jpp  15.0.4
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JFitL1dtSlices.cc
Go to the documentation of this file.
2 #include "JDetector/JDetector.hh"
4 #include "Jeep/JParser.hh"
5 
6 #include "TF2.h"
7 #include "TFile.h"
8 #include "TH1D.h"
9 #include "TH2D.h"
10 #include "TMath.h"
11 
12 #include <iostream>
13 #include <string>
14 
15 
16 namespace FITL1DTSLICES {
17 
18  /* Calculate the log-poisson
19  *
20  * This function has some assumptions: since the input data from JMonitorL1dt use SN datastream
21  * most background is strongly suppressed. This means many bins on both data and model will
22  * have no entries: n = 0 or nhat = 0. This will result in infinities when calculating the
23  * log-poisson. To avoid the infities the options below are used.
24  * This has been crosschecked with L1 datastream which gives a significant background and the
25  * shapes of the likelihood distributions are similar.
26  */
27  double logPoison(double n, double nhat, double logP_min = -999999.9) {
28  if (nhat < 0.0 || n < 0.0) {
29  FATAL("logPoisson: n (" << n <<") or nhat (" << nhat << ") is < 0.0" << std::endl);
30  }
31  if (n == 0.0 && nhat == 0.0) {
32  return 0.0; // ok.
33  }
34  if (n == 0.0 && nhat > 0.0) {
35  return -1*nhat; // log( lab^0 * e^-lab / 0! ) = log( e^-lab ) = -lab, ok.
36  }
37  if (n >= 0.0 && nhat == 0.0) {
38  return 0.0; // should be -inf, we only consider bins where we expect non-zero entries: no expected entries we throw away.
39  }
40  Double_t poisson = TMath::Poisson(n, nhat);
41  if (poisson == 0) return 0.0; // only when model and data are close enough do we care about their values, so if P=0, do nothing, since we are not in the region we care about.
42  else return TMath::Log(poisson);
43  }
44 
45  double getLogQuality(const TH1D* data, const TH1D* model, int di, double data_bkg = 0.0001, double model_bkg = 0.0001) {
46 
47  double q = 0.0;
48 
49  // If you calculate over the whole histogram you lose bins at edges: only when dt=0 are all N bins considered, everywhere else
50  // either data or model bins get lost at the edge. So integrate over a smaller range:
51  int i_low = data->GetNbinsX() / 4; // Start at a quarter of the hist
52  int i_hgh = data->GetNbinsX() / 4 * 3; // End at three quarters
53  for (int i = i_low; i <= i_hgh; ++i) {
54  double n = data_bkg;
55  double nhat = model_bkg;
56  if (i >= 1 && i <= data->GetNbinsX()) {
57  // Count in bin i for the data
58  n = data->GetBinContent(i);
59  }
60  if (i+di >= 1 && i+di <= model->GetNbinsX()) {
61  // Count in bin i+di for the model
62  nhat = model->GetBinContent(i+di);
63  }
64  double logP = logPoison(n, nhat);
65  q += logP; // Added 0 for when the background is zero in model and data, but then P(n|0) gives inf! So removed overwritten logPmin
66  }
67  return q;
68  }
69 }
70 
71 
72 /**
73  * \file
74  *
75  * Auxiliary program to fit L1dt output data to an L1dt model, both created with JMonitorL1dt
76  *
77  * \author kmelis, lnauta
78  */
79 int main(int argc, char **argv)
80 {
81  using namespace std;
82  using namespace FITL1DTSLICES;
83  using namespace JPP;
84 
85  string inputFile;
86  string modelFile;
87  string outputFile;
88  string detectorFile;
89  double dt_fitrange;
90  int rebin;
91  double TMax_ns;
92  bool normaliseMC;
93  int debug;
94 
95  try {
96 
97  JParser<> zap("Program to calculate log-likelihood distributions between DOM pairs and fit a poly2 to find the shape, used in FitL1dt to find time offsets");
98 
99  zap['f'] = make_field(inputFile, "input file") = "monitor.root";
100  zap['m'] = make_field(modelFile, "model file");
101  zap['o'] = make_field(outputFile, "output file") = "fitslices.root";
102  zap['a'] = make_field(detectorFile, "detector file");
103  zap['T'] = make_field(TMax_ns, "scan range around 0 in data-model comparison [ns]") = 50.0;
104  zap['t'] = make_field(dt_fitrange, "fitrange of polynomial to quality [ns]") = 5.0;
105  zap['r'] = make_field(rebin, "rebin the histograms") = 1;
106  zap['N'] = make_field(normaliseMC, "normalize the MC histogram");
107  zap['d'] = make_field(debug) = 0;
108 
109  if (zap.read(argc, argv) != 0) {
110  return 1;
111  }
112  }
113  catch(const exception &error) {
114  FATAL(error.what() << endl);
115  }
116 
117 
118 
119  const double TSignal_ns = 750.0; // The width of the signal for integrating out the background
120 
122 
123  try {
124  load(detectorFile, detector);
125  }
126  catch(const JException& error) {
127  FATAL(error);
128  }
129 
130 
131  TFile* in = TFile::Open(inputFile.c_str(), "exist");
132  TFile* in_model = TFile::Open(modelFile.c_str(), "exist");
133 
134  if (in == NULL || !in->IsOpen()) {
135  FATAL("File: " << inputFile << " not opened." << endl);
136  }
137  if (in_model == NULL || !in_model->IsOpen()) {
138  FATAL("File: " << modelFile << " not opened." << endl);
139  }
140 
141  TFile out(outputFile.c_str(), "recreate");
142 
143  TH2D h2A("Aij", "Aij", detector.size(), -0.5, detector.size()-0.5, detector.size(), -0.5, detector.size()-0.5);
144  TH2D h2B("Bij", "Bij", detector.size(), -0.5, detector.size()-0.5, detector.size(), -0.5, detector.size()-0.5);
145  TH2D h2C("Cij", "Cij", detector.size(), -0.5, detector.size()-0.5, detector.size(), -0.5, detector.size()-0.5);
146 
147  // The bare histogram can be used for debugging, it contains the actual time offsets.
148  // Histogram Bij contains Bij/2Aij due to a different representation used in the code compared to the mathematical formulation.
149  TH2D h2Bbare("Bij_bare", "Bij_bare", detector.size(), -0.5, detector.size()-0.5, detector.size(), -0.5, detector.size()-0.5);
150  for (JDetector::iterator moduleA = detector.begin(); moduleA != detector.end(); ++moduleA) {
151  const int idom = distance(detector.begin(), moduleA);
152  TString label = Form("%i", moduleA->getID());
153  h2A.GetXaxis()->SetBinLabel(idom+1, label);
154  h2A.GetYaxis()->SetBinLabel(idom+1, label);
155  h2B.GetXaxis()->SetBinLabel(idom+1, label);
156  h2B.GetYaxis()->SetBinLabel(idom+1, label);
157  h2C.GetXaxis()->SetBinLabel(idom+1, label);
158  h2C.GetYaxis()->SetBinLabel(idom+1, label);
159 
160  h2Bbare.GetXaxis()->SetBinLabel(idom+1, label);
161  h2Bbare.GetYaxis()->SetBinLabel(idom+1, label);
162  }
163 
164  for (JDetector::iterator moduleA = detector.begin(); moduleA != detector.end(); ++moduleA) {
165 
166  DEBUG("Module " << moduleA->getID() << endl);
167  const JString title = JString(moduleA->getID());
168  const int idom = distance(detector.begin(), moduleA);
169 
170  // data
171  TH2D* d2s = (TH2D*) in->Get((title + ".2S").c_str());
172  TH1D* d1l = (TH1D*) in->Get((title + ".1L").c_str());
173 
174  if (d2s == NULL || d1l == NULL) {
175  WARNING("No data in data histogram " << title << endl);
176  continue;
177  }
178 
179  // model
180  TH2D* m2s = (TH2D*) in_model->Get((title + ".2S").c_str());
181  TH1D* m1l = (TH1D*) in_model->Get((title + ".1L").c_str());
182 
183  if (m2s == NULL || m1l == NULL) {
184  WARNING("No mata in model histogram " << title << endl);
185  continue;
186  }
187 
188 
189  for (JDetector::iterator moduleB = moduleA; moduleB != detector.end(); ++moduleB) {
190  const int jdom = distance(detector.begin(), moduleB);
191 
192  if (moduleB->getID() == moduleA->getID()) {
193  continue;
194  }
195  if (moduleB->getString() != moduleA->getString()) {
196  continue;
197  }
198 
199  TH1D* d1s = d2s->ProjectionY((title + Form(".%i.d1S", moduleB->getID())).c_str(),
200  d2s->GetXaxis()->FindBin(TString(Form("%i", moduleB->getID()))),
201  d2s->GetXaxis()->FindBin(TString(Form("%i", moduleB->getID()))), "e");
202  TH1D* m1s = m2s->ProjectionY((title + Form(".%i.m1S", moduleB->getID())).c_str(),
203  m2s->GetXaxis()->FindBin(TString(Form("%i", moduleB->getID()))),
204  m2s->GetXaxis()->FindBin(TString(Form("%i", moduleB->getID()))), "e");
205 
206  if (d1s->Integral() <= 0.0 || m1s->Integral() <= 0.0) {
207  continue;
208  }
209 
210  // background fit: data
211  double binCenterMaxBin = d1s->GetXaxis()->GetBinCenter(d1s->GetMaximumBin());
212  double backgroundrate_s = 0.0;
213  int nbins = 0;
214  for (int j = 1; j <= d1s->GetXaxis()->GetNbins(); ++j) {
215  if (fabs(d1s->GetXaxis()->GetBinCenter(j) - binCenterMaxBin) > TSignal_ns ) {
216  backgroundrate_s += d1s->GetBinContent(j);
217  nbins++;
218  }
219  }
220  backgroundrate_s /= nbins;
221 
222  // background fit: model
223  binCenterMaxBin = m1s->GetXaxis()->GetBinCenter(m1s->GetMaximumBin());
224  double backgroundrate_m = 0.0;
225  nbins = 0;
226  for (int j = 1; j <= m1s->GetXaxis()->GetNbins(); ++j) {
227  if (fabs(m1s->GetXaxis()->GetBinCenter(j) - binCenterMaxBin) > TSignal_ns ) {
228  backgroundrate_m += m1s->GetBinContent(j);
229  nbins++;
230  }
231  }
232  backgroundrate_m /= nbins;
233 
234  // scale livetimes and remove background
235  const double Ld = d1l->GetBinContent(d1l->GetXaxis()->FindBin(TString(Form("%i", moduleB->getID()))));
236  const double Lm = m1l->GetBinContent(m1l->GetXaxis()->FindBin(TString(Form("%i", moduleB->getID()))));
237 
238  for (int j = 1; j <= d1s->GetXaxis()->GetNbins(); ++j) {
239  const double y = d1s->GetBinContent(j);
240  const double dy = d1s->GetBinError(j);
241 
242  d1s->SetBinContent(j, y);
243  d1s->SetBinError(j, dy);
244  }
245 
246  // scale the livetime of the model
247  for (int j = 1; j <= m1s->GetXaxis()->GetNbins(); ++j) {
248  const double y = m1s->GetBinContent(j);
249  const double dy = m1s->GetBinError(j);
250 
251  m1s->SetBinContent(j, y * Ld/Lm);
252  m1s->SetBinError(j, dy*Ld/Lm);
253  }
254 
255 
256  if (normaliseMC) {
257  const double scalefactor = d1s->Integral()/m1s->Integral();
258  m1s->Scale(scalefactor);
259  }
260 
261  // rebin
262  d1s->Rebin(rebin);
263  m1s->Rebin(rebin);
264 
265  // write background subtracted, normalised slices
266  d1s->Write();
267  m1s->Write();
268 
269  // get quality
270  int ddt_min = -0.5*TMax_ns; // maximum time shift [ns]
271  int ddt_max = 0.5*TMax_ns; // maximum time shift [ns]
272 
273  TH1D q1((title + Form(".%i.q1", moduleB->getID())).c_str(), (title + Form(".%i.q1", moduleB->getID())).c_str(), (ddt_max-ddt_min)/rebin+1, ddt_min-0.5*rebin, ddt_max+0.5*rebin);
274 
275  for (int di = 1; di <= q1.GetNbinsX(); di++) {
276  q1.SetBinContent(di, getLogQuality(d1s, m1s, (int)(q1.GetBinCenter(di)/rebin), 0.0, 0.0));
277  }
278 
279  // Fit peak with parabola
280  const double dt_fitmid = q1.GetBinCenter(q1.GetMaximumBin());
281 
282  TF1 q1f((title + Form(".%i.q1f", moduleB->getID())).c_str(), "[0]*x*x + [1]*x + [2]", dt_fitmid - dt_fitrange, dt_fitmid + dt_fitrange);
283  q1f.SetParLimits(0, -1e5, 0.0); // Added as to improve fits. Aij should always be l.t. 0 to fit a maximum to the LLH
284 
285  string fitoptions = "R";
286  if (debug < JEEP::debug_t) {
287  fitoptions += " Q";
288  }
289  q1.Fit(&q1f, fitoptions.c_str());
290 
291  // write quality and fit
292  q1.Write();
293 
294  // Fit results to global fit matrix
295  const double Aij = q1f.GetParameter(0);
296  const double Bij = q1f.GetParameter(1);
297  const double Cij = q1f.GetParameter(2);
298 
299  double bareBij = Bij / Aij;
300  if (Aij == 0) bareBij = 0;
301 
302  if (fabs(-0.5*Bij/Aij - dt_fitmid) > 3 || fabs(-0.5*Bij/Aij) > dt_fitrange) {
303  WARNING(" Please check fit of pair " << idom << ", " << jdom << " : " << moduleA->getID() << ", " << moduleB->getID() << endl);
304  WARNING(" DOMpair " << idom << "; " << jdom << " A=" << Aij << " B=" << Bij << " , max at: " << dt_fitmid << " [ns], best fit at " << -0.5*Bij/Aij << " ns " << endl);
305  }
306  DEBUG("DOMpair " << idom << "; " << jdom << " A=" << Aij << " B=" << Bij << " , max at: " << dt_fitmid << " [ns], best fit at " << -0.5*Bij/Aij << " ns " << endl);
307 
308  h2A.SetBinContent(idom+1, jdom+1, Aij);
309  h2A.SetBinContent(jdom+1, idom+1, Aij);
310  h2B.SetBinContent(idom+1, jdom+1, Bij);
311  h2B.SetBinContent(jdom+1, idom+1,-Bij);
312  h2C.SetBinContent(idom+1, jdom+1, Cij);
313  h2C.SetBinContent(jdom+1, idom+1, Cij);
314 
315  h2Bbare.SetBinContent(idom+1, jdom+1, bareBij);
316  h2Bbare.SetBinContent(jdom+1, idom+1, 0);
317  }
318  }
319 
320  out.cd();
321 
322  h2A.Write();
323  h2B.Write();
324  h2C.Write();
325 
326  h2Bbare.Write();
327 
328  out.Close();
329 }
Utility class to parse command line options.
Definition: JParser.hh:1500
General exception.
Definition: JException.hh:23
#define WARNING(A)
Definition: JMessage.hh:65
debug
Definition: JMessage.hh:29
int main(int argc, char *argv[])
Definition: Main.cc:15
Wrapper class around STL string class.
Definition: JString.hh:27
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Detector data structure.
Definition: JDetector.hh:89
string outputFile
Data structure for detector geometry and calibration.
const int n
Definition: JPolint.hh:660
Detector file.
Definition: JHead.hh:196
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
int debug
debug level
Definition: JSirene.cc:63
#define FATAL(A)
Definition: JMessage.hh:67
int read(const int argc, const char *const argv[])
Parse the program&#39;s command line options.
Definition: JParser.hh:1798
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
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
int j
Definition: JPolint.hh:666
do set_variable DETECTOR_TXT $WORKDIR detector
KM3NeT DAQ constants, bit handling, etc.
double logPoison(double n, double nhat, double logP_min=-999999.9)
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 source JAcoustics sh $DETECTOR_ID CHECK_EXIT_CODE typeset A TRIPODS get_tripods $WORKDIR tripod txt TRIPODS for EMITTER in
Definition: JCanberra.sh:41
double getLogQuality(const TH1D *data, const TH1D *model, int di, double data_bkg=0.0001, double model_bkg=0.0001)