Jpp  19.1.0
the software that should make you happy
Namespaces | Functions
JFitL1dtSlices.cc File Reference

Auxiliary program to fit L1dt output data to an L1dt model, both created with JMonitorL1dt. More...

#include "km3net-dataformat/online/JDAQ.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "Jeep/JParser.hh"
#include "TF2.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TMath.h"
#include "JROOT/JMinimizer.hh"
#include <iostream>
#include <string>

Go to the source code of this file.

Namespaces

 FITL1DTSLICES
 

Functions

double FITL1DTSLICES::logPoison (double n, double nhat, double logP_min=-999999.9)
 
double FITL1DTSLICES::getLogQuality (const TH1D *data, const TH1D *model, int di, double data_bkg=0.0001, double model_bkg=0.0001)
 
int main (int argc, char **argv)
 

Detailed Description

Auxiliary program to fit L1dt output data to an L1dt model, both created with JMonitorL1dt.

Author
kmelis, lnauta

Definition in file JFitL1dtSlices.cc.

Function Documentation

◆ main()

int main ( int  argc,
char **  argv 
)

Definition at line 81 of file JFitL1dtSlices.cc.

82 {
83  using namespace std;
84  using namespace FITL1DTSLICES;
85  using namespace JPP;
86 
87  string inputFile;
88  string modelFile;
89  string outputFile;
90  string detectorFile;
91  double dt_fitrange;
92  int rebin;
93  double TMax_ns;
94  bool normaliseMC;
95  int debug;
96 
97  try {
98 
99  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");
100 
101  zap['f'] = make_field(inputFile, "input file") = "monitor.root";
102  zap['m'] = make_field(modelFile, "model file");
103  zap['o'] = make_field(outputFile, "output file") = "fitslices.root";
104  zap['a'] = make_field(detectorFile, "detector file");
105  zap['T'] = make_field(TMax_ns, "scan range around 0 in data-model comparison [ns]") = 50.0;
106  zap['t'] = make_field(dt_fitrange, "fitrange of polynomial to quality [ns]") = 5.0;
107  zap['r'] = make_field(rebin, "rebin the histograms") = 1;
108  zap['N'] = make_field(normaliseMC, "normalize the MC histogram");
109  zap['d'] = make_field(debug) = 0;
110 
111  if (zap.read(argc, argv) != 0) {
112  return 1;
113  }
114  }
115  catch(const exception &error) {
116  FATAL(error.what() << endl);
117  }
118 
119 
120  const double TSignal_ns = 750.0; // The width of the signal for integrating out the background
121 
122  JDetector detector;
123 
124  try {
125  load(detectorFile, detector);
126  }
127  catch(const JException& error) {
128  FATAL(error);
129  }
130 
131 
132  TFile* in = TFile::Open(inputFile.c_str(), "exist");
133  TFile* in_model = TFile::Open(modelFile.c_str(), "exist");
134 
135  if (in == NULL || !in->IsOpen()) {
136  FATAL("File: " << inputFile << " not opened." << endl);
137  }
138  if (in_model == NULL || !in_model->IsOpen()) {
139  FATAL("File: " << modelFile << " not opened." << endl);
140  }
141 
142  TFile out(outputFile.c_str(), "recreate");
143 
144  TH2D h2A("Aij", "Aij", detector.size(), -0.5, detector.size()-0.5, detector.size(), -0.5, detector.size()-0.5);
145  TH2D h2B("Bij", "Bij", detector.size(), -0.5, detector.size()-0.5, detector.size(), -0.5, detector.size()-0.5);
146  TH2D h2C("Cij", "Cij", detector.size(), -0.5, detector.size()-0.5, detector.size(), -0.5, detector.size()-0.5);
147 
148  // The bare histogram can be used for debugging, it contains the actual time offsets.
149  // Histogram Bij contains Bij/2Aij due to a different representation used in the code compared to the mathematical formulation.
150  TH2D h2Bbare("Bij_bare", "Bij_bare", detector.size(), -0.5, detector.size()-0.5, detector.size(), -0.5, detector.size()-0.5);
151  for (JDetector::iterator moduleA = detector.begin(); moduleA != detector.end(); ++moduleA) {
152  const int idom = distance(detector.begin(), moduleA);
153  TString label = Form("%i", moduleA->getID());
154  h2A.GetXaxis()->SetBinLabel(idom+1, label);
155  h2A.GetYaxis()->SetBinLabel(idom+1, label);
156  h2B.GetXaxis()->SetBinLabel(idom+1, label);
157  h2B.GetYaxis()->SetBinLabel(idom+1, label);
158  h2C.GetXaxis()->SetBinLabel(idom+1, label);
159  h2C.GetYaxis()->SetBinLabel(idom+1, label);
160 
161  h2Bbare.GetXaxis()->SetBinLabel(idom+1, label);
162  h2Bbare.GetYaxis()->SetBinLabel(idom+1, label);
163  }
164 
165  for (JDetector::iterator moduleA = detector.begin(); moduleA != detector.end(); ++moduleA) {
166 
167  DEBUG("Module " << moduleA->getID() << endl);
168  const JString title = JString(moduleA->getID());
169  const int idom = distance(detector.begin(), moduleA);
170 
171  // data
172  TH2D* d2s = (TH2D*) in->Get((title + ".2S").c_str());
173  TH1D* d1l = (TH1D*) in->Get((title + ".1L").c_str());
174 
175  if (d2s == NULL || d1l == NULL) {
176  WARNING("No data in data histogram " << title << endl);
177  continue;
178  }
179 
180  // model
181  TH2D* m2s = (TH2D*) in_model->Get((title + ".2S").c_str());
182  TH1D* m1l = (TH1D*) in_model->Get((title + ".1L").c_str());
183 
184  if (m2s == NULL || m1l == NULL) {
185  WARNING("No mata in model histogram " << title << endl);
186  continue;
187  }
188 
189 
190  for (JDetector::iterator moduleB = moduleA; moduleB != detector.end(); ++moduleB) {
191  const int jdom = distance(detector.begin(), moduleB);
192 
193  if (moduleB->getID() == moduleA->getID()) {
194  continue;
195  }
196  if (moduleB->getString() != moduleA->getString()) {
197  continue;
198  }
199 
200  TH1D* d1s = d2s->ProjectionY((title + Form(".%i.d1S", moduleB->getID())).c_str(),
201  d2s->GetXaxis()->FindBin(TString(Form("%i", moduleB->getID()))),
202  d2s->GetXaxis()->FindBin(TString(Form("%i", moduleB->getID()))), "e");
203  TH1D* m1s = m2s->ProjectionY((title + Form(".%i.m1S", moduleB->getID())).c_str(),
204  m2s->GetXaxis()->FindBin(TString(Form("%i", moduleB->getID()))),
205  m2s->GetXaxis()->FindBin(TString(Form("%i", moduleB->getID()))), "e");
206 
207  if (d1s->Integral() <= 0.0 || m1s->Integral() <= 0.0) {
208  continue;
209  }
210 
211  // background fit: data
212  double binCenterMaxBin = d1s->GetXaxis()->GetBinCenter(d1s->GetMaximumBin());
213  double backgroundrate_s = 0.0;
214  int nbins = 0;
215  for (int j = 1; j <= d1s->GetXaxis()->GetNbins(); ++j) {
216  if (fabs(d1s->GetXaxis()->GetBinCenter(j) - binCenterMaxBin) > TSignal_ns ) {
217  backgroundrate_s += d1s->GetBinContent(j);
218  nbins++;
219  }
220  }
221  backgroundrate_s /= nbins;
222 
223  // background fit: model
224  binCenterMaxBin = m1s->GetXaxis()->GetBinCenter(m1s->GetMaximumBin());
225  double backgroundrate_m = 0.0;
226  nbins = 0;
227  for (int j = 1; j <= m1s->GetXaxis()->GetNbins(); ++j) {
228  if (fabs(m1s->GetXaxis()->GetBinCenter(j) - binCenterMaxBin) > TSignal_ns ) {
229  backgroundrate_m += m1s->GetBinContent(j);
230  nbins++;
231  }
232  }
233  backgroundrate_m /= nbins;
234 
235  // scale livetimes and remove background
236  const double Ld = d1l->GetBinContent(d1l->GetXaxis()->FindBin(TString(Form("%i", moduleB->getID()))));
237  const double Lm = m1l->GetBinContent(m1l->GetXaxis()->FindBin(TString(Form("%i", moduleB->getID()))));
238 
239  for (int j = 1; j <= d1s->GetXaxis()->GetNbins(); ++j) {
240  const double y = d1s->GetBinContent(j);
241  const double dy = d1s->GetBinError(j);
242 
243  d1s->SetBinContent(j, y);
244  d1s->SetBinError(j, dy);
245  }
246 
247  // scale the livetime of the model
248  for (int j = 1; j <= m1s->GetXaxis()->GetNbins(); ++j) {
249  const double y = m1s->GetBinContent(j);
250  const double dy = m1s->GetBinError(j);
251 
252  m1s->SetBinContent(j, y * Ld/Lm);
253  m1s->SetBinError(j, dy*Ld/Lm);
254  }
255 
256 
257  if (normaliseMC) {
258  const double scalefactor = d1s->Integral()/m1s->Integral();
259  m1s->Scale(scalefactor);
260  }
261 
262  // rebin
263  d1s->Rebin(rebin);
264  m1s->Rebin(rebin);
265 
266  // write background subtracted, normalised slices
267  d1s->Write();
268  m1s->Write();
269 
270  // get quality
271  int ddt_min = -0.5*TMax_ns; // maximum time shift [ns]
272  int ddt_max = 0.5*TMax_ns; // maximum time shift [ns]
273 
274  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);
275 
276  for (int di = 1; di <= q1.GetNbinsX(); di++) {
277  q1.SetBinContent(di, getLogQuality(d1s, m1s, (int)(q1.GetBinCenter(di)/rebin), 0.0, 0.0));
278  }
279 
280  // Fit peak with parabola
281  const double dt_fitmid = q1.GetBinCenter(q1.GetMaximumBin());
282 
283  TF1 q1f((title + Form(".%i.q1f", moduleB->getID())).c_str(), "[0]*x*x + [1]*x + [2]", dt_fitmid - dt_fitrange, dt_fitmid + dt_fitrange);
284  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
285 
286  string fitoptions = "R";
287  if (debug < JEEP::debug_t) {
288  fitoptions += " Q";
289  }
290  q1.Fit(&q1f, fitoptions.c_str());
291 
292  // write quality and fit
293  q1.Write();
294 
295  // Fit results to global fit matrix
296  const double Aij = q1f.GetParameter(0);
297  const double Bij = q1f.GetParameter(1);
298  const double Cij = q1f.GetParameter(2);
299 
300  double bareBij = Bij / Aij;
301  if (Aij == 0) bareBij = 0;
302 
303  if (fabs(-0.5*Bij/Aij - dt_fitmid) > 3 || fabs(-0.5*Bij/Aij) > dt_fitrange) {
304  WARNING(" Please check fit of pair " << idom << ", " << jdom << " : " << moduleA->getID() << ", " << moduleB->getID() << endl);
305  WARNING(" DOMpair " << idom << "; " << jdom << " A=" << Aij << " B=" << Bij << " , max at: " << dt_fitmid << " [ns], best fit at " << -0.5*Bij/Aij << " ns " << endl);
306  }
307  DEBUG("DOMpair " << idom << "; " << jdom << " A=" << Aij << " B=" << Bij << " , max at: " << dt_fitmid << " [ns], best fit at " << -0.5*Bij/Aij << " ns " << endl);
308 
309  h2A.SetBinContent(idom+1, jdom+1, Aij);
310  h2A.SetBinContent(jdom+1, idom+1, Aij);
311  h2B.SetBinContent(idom+1, jdom+1, Bij);
312  h2B.SetBinContent(jdom+1, idom+1,-Bij);
313  h2C.SetBinContent(idom+1, jdom+1, Cij);
314  h2C.SetBinContent(jdom+1, idom+1, Cij);
315 
316  h2Bbare.SetBinContent(idom+1, jdom+1, bareBij);
317  h2Bbare.SetBinContent(jdom+1, idom+1, 0);
318  }
319  }
320 
321  out.cd();
322 
323  h2A.Write();
324  h2B.Write();
325  h2C.Write();
326 
327  h2Bbare.Write();
328 
329  out.Close();
330 }
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
DEBUG
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
debug
int debug
debug level
Definition: JSirene.cc:69
WARNING
#define WARNING(A)
Definition: JMessage.hh:65
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:2142
distance
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Definition: PhysicsEvent.hh:434
JDETECTOR::JDetector
Detector data structure.
Definition: JDetector.hh:96
JLANG::JException
General exception.
Definition: JException.hh:24
JLANG::JString
Wrapper class around STL string class.
Definition: JString.hh:29
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1698
FITL1DTSLICES::getLogQuality
double getLogQuality(const TH1D *data, const TH1D *model, int di, double data_bkg=0.0001, double model_bkg=0.0001)
Definition: JFitL1dtSlices.cc:47
JDETECTOR::load
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Definition: JDetectorToolkit.hh:505
JEEP::debug_t
@ debug_t
debug
Definition: JMessage.hh:29
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JAAnetToolkit.hh:43
JTOOLS::j
int j
Definition: JPolint.hh:792
std
Definition: JSTDTypes.hh:14
JAANET::detector
Detector file.
Definition: JHead.hh:227
Generated by doxygen 1.9.1