Jpp 19.3.0-rc.1
the software that should make you happy
Loading...
Searching...
No Matches
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

namespace  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:72
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:508
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:801
JAANET::detector
Detector file.
Definition JHead.hh:227
Generated by doxygen 1.12.0