Jpp 19.3.0-rc.1
the software that should make you happy
Loading...
Searching...
No Matches
Namespaces | Functions
JFitL1dt.cc File Reference

Auxiliary program to find the detector time offsets from FitL1dtSlices output. More...

#include "km3net-dataformat/online/JDAQ.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JMath/JMatrixNS.hh"
#include "Jeep/JParser.hh"
#include "TF2.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TMath.h"
#include "TMatrixDSym.h"
#include "TROOT.h"
#include "TVectorD.h"
#include <string>
#include <iostream>
#include <set>

Go to the source code of this file.

Namespaces

namespace  FITL1DT
 

Functions

JMatrixNS FITL1DT::removeRowColumn (const JMatrixNS &A, int C)
 
vector< double > FITL1DT::dotprod (const JMatrixNS &A, const vector< double > &x)
 
int main (int argc, char **argv)
 

Detailed Description

Auxiliary program to find the detector time offsets from FitL1dtSlices output.

Author
kmelis, lnauta

Definition in file JFitL1dt.cc.

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

Definition at line 71 of file JFitL1dt.cc.

72{
73 using namespace std;
74 using namespace JPP;
75 using namespace FITL1DT;
76
77 string inputFile;
78 string outputFile;
79 string detectorFile;
80 bool overwriteDetector;
81 bool noInterDU;
82 int maxFloors;
83 string option;
84 int idom_ref;
85 int debug;
86
87 try {
88
89 JParser<> zap("Program to calculate time offsets from FitL1dtSlices output");
90
91 zap['f'] = make_field(inputFile, "input file") = "fitslices.root";
92 zap['o'] = make_field(outputFile, "output file") = "fit.root";
93 zap['a'] = make_field(detectorFile, "detector file");
94 zap['A'] = make_field(overwriteDetector, "overwrite the input detector file");
95 zap['D'] = make_field(noInterDU, "do not use inter DU information in the calculation");
96 zap['r'] = make_field(idom_ref, "reference DOM set to t=0") = 9;
97 zap['F'] = make_field(maxFloors, "number of floors used in fit") = 9999;
98 zap['d'] = make_field(debug) = 0;
99
100 if (zap.read(argc, argv) != 0) {
101 return 1;
102 }
103 }
104 catch(const exception &error) {
105 FATAL(error.what() << endl);
106 }
107
108
109 JDetector detector;
110
111 try {
112 load(detectorFile, detector);
113 }
114 catch(const JException& error) {
115 FATAL(error);
116 }
117
118 set<int> DU_IDs;
119 for (JDetector::iterator moduleA = detector.begin(); moduleA != detector.end(); ++moduleA) {
120 DU_IDs.insert(moduleA->getString());
121 if (!noInterDU) {
122 break;
123 }
124 }
125
126 TFile* in = TFile::Open(inputFile.c_str(), "exist");
127
128 if (in == NULL || !in->IsOpen()) {
129 FATAL("File: " << inputFile << " not opened." << endl);
130 }
131
132 // Get fitting histograms
133 TH2D* h2A = (TH2D*)in->Get("Aij");
134 TH2D* h2B = (TH2D*)in->Get("Bij");
135
136 if (h2A == NULL || h2B == NULL) {
137 FATAL("File does not contain histogram with matrix. Please use FitL1dtSlices script first." << endl);
138 }
139
140 TFile out(outputFile.c_str(), "recreate");
141
142 h2A->Write();
143 h2B->Write();
144
145 vector<double> dt0(detector.size(), 0.0);
146 vector<double> vart0(detector.size(), 0.0);
147 vector<int> status(detector.size(), 0);
148
149 for (set<int>::const_iterator DU_ID = DU_IDs.begin(); DU_ID != DU_IDs.end(); ++DU_ID) {
150 JMatrixNS fitA(detector.size());
151 vector<double> fitB(detector.size());
152 vector<double> diag(detector.size());
153
154 // Fill the matrix and vector
155 for (JDetector::iterator moduleA = detector.begin(); moduleA != detector.end(); ++moduleA) {
156 const int idom = distance(detector.begin(), moduleA);
157 const int ibin = h2A->GetXaxis()->FindBin(TString(Form("%i", moduleA->getID())));
158
159 if (noInterDU && (moduleA->getString() != *DU_ID)) {
160 continue;
161 }
162
163 for (JDetector::iterator moduleB = detector.begin(); moduleB != detector.end(); ++moduleB) {
164 const int jdom = distance(detector.begin(), moduleB);
165 const int jbin = h2A->GetXaxis()->FindBin(TString(Form("%i", moduleB->getID())));
166
167 if (idom == jdom) {
168 continue;
169 }
170
171 if (noInterDU && (moduleB->getString() != *DU_ID)) {
172 continue;
173 }
174
175 if (fabs(moduleA->getFloor() - moduleB->getFloor()) > maxFloors) {
176 continue;
177 }
178
179 double Aij = h2A->GetBinContent(ibin, jbin);
180 double Bij = h2B->GetBinContent(ibin, jbin);
181
182 // If the parabola does not have a maximum: Aij > 0, this means the fit in FitL1dtSlices was unsuccesful.
183 if (Aij >= 0.0) {
184 continue;
185 }
186
187 NOTICE(" " << idom << " " << jdom << " " << moduleA->getID() << " " << moduleB->getID() << " Aij: " << Aij << " Bij: " << Bij << " bf: " << -0.5*Bij/Aij << endl);
188
189 fitA.at(idom, jdom) = -4.0*Aij;
190 fitA.at(idom, idom) += 4.0*Aij;
191 fitB[idom] += 2.0*Bij; // Here Bij = -2*Aij*Bij in the mathematical derivation due to the form of the TF1 in FitL1Slices: [0]x^2+[1]x+[2] instead of [0](x-[1])^2+[2]
192 diag[idom] += 4.0*Aij; // Wilks theorem: d^2L/dp^2 = 1/var(p), from this the diagonal contains the second derivative of logL.
193 }
194 }
195
196 // remove empty columns and rows (no data and thus no time info)
197 unsigned int icolumn = 0;
198 int idom = 0;
199 vector<bool> changet0(detector.size(), true);
200 while (icolumn < fitA.size()) {
201 if (fitA.at(icolumn, icolumn) == 0) {
202 fitA = removeRowColumn(fitA, icolumn);
203 fitB.erase(fitB.begin() + icolumn);
204
205 changet0[idom] = false;
206 } else {
207 icolumn++;
208 }
209 idom++;
210 }
211
212 // no data for this string
213 if (fitA.size() <= 1) {
214 continue;
215 }
216
217 // Use ROOT matrices from here to have more options
218 // runtime is not important for this program and the matrix functionality defined before main() is not trivial to port to ROOT.
219
220 TMatrixDSym fit_A(fitA.size());
221 fit_A.SetTol(1e-26); // Det is near 1E-40, so set the tolerance 10 orders lower for safety.
222 TVectorD fit_B(fitB.size());
223
224 for (unsigned int i = 0; i < fitA.size(); i++) {
225 for (unsigned int j = 0; j < fitA.size(); j++) {
226 fit_A(i,j) = fitA.at(i,j);
227 }
228 fit_B(i) = fitB[i];
229 }
230
231 // fix reference DOM time offset (the first DOM with data, i.e. with off-diagonal elements > 0.0)
232 for (int i = 0; i < fit_A.GetNrows(); ++i) {
233 fit_A(idom_ref-1, i) = 0.0;
234 fit_A(i, idom_ref-1) = 0.0;
235 }
236 fit_A(idom_ref-1, idom_ref-1) = 1.0;
237 fit_B(idom_ref-1) = 0.0;
238
239 if (debug >= 2) fit_A.Print();
240
241 // Fit best time offsets via matrix inversion
242 Double_t fit_A_det = -999; // dont initialize to 0. det == 0 means the matrix is non-invertable.
243 fit_A.Invert(&fit_A_det);
244 if (fit_A_det == 0) {
245 NOTICE("Matrix not invertible." << endl);
246 }
247
248 if (debug >= 2) fit_A.Print();
249
250 // dt = A^-1*B
251 TVectorD dt_i = fit_A * fit_B;
252
253 // This shape for writing out the dt0 is chose because row+columns are removed when they have a diagonal element of 0, i.e. make the matrix singular.
254 int k = 0;
255 for (JDetector::iterator moduleA = detector.begin(); moduleA != detector.end(); ++moduleA) {
256 const int idom = distance(detector.begin(), moduleA);
257 if (changet0[idom]) {
258 dt0[idom] = dt_i(k);
259 vart0[idom] = TMath::Sqrt(-1 / diag[idom]);
260 status[idom] = 2; // fitted
261 if (k == 0) {
262 status[idom] = 1;
263 } // fixed t0
264 k++;
265 }
266 }
267 }
268
269 TH1D h_dt0("dt0", "dt0", detector.size(), -0.5, detector.size()-0.5);
270
271 for (JDetector::iterator moduleA = detector.begin(); moduleA != detector.end(); ++moduleA) {
272 const int idom = distance(detector.begin(), moduleA);
273 NOTICE("Changing t0 of DOM " << moduleA->getID() << " (S"<<setw(2) << moduleA->getString() << "F" << setw(2) << moduleA->getFloor() << ") by " << setw(13) << dt0[idom] << " [ns] ");
274
275 if (status[idom] == 0) {
276 NOTICE(" ** unable to fit **" << endl);
277 }
278 else if (moduleA->getFloor() == idom_ref) {
279 NOTICE(" ** fixed **" << endl);
280 }
281 else {
282 NOTICE(endl);
283 }
284
285 if (overwriteDetector) {
286 moduleA->add(dt0[idom]); // minus sign difference with JHobbit.cc
287 }
288 h_dt0.GetXaxis()->SetBinLabel(idom+1, Form("%i", moduleA->getID()));
289 h_dt0.SetBinContent(idom+1, dt0[idom]);
290 h_dt0.SetBinError(idom+1, vart0[idom]);
291 }
292
293 h_dt0.Write();
294
295 out.Close();
296
297 if (overwriteDetector) {
298 NOTICE("Store calibration data on file " << detectorFile << endl);
299
300 store(detectorFile, detector);
301 }
302}
outputFile
string outputFile
Definition JDAQTimesliceSelector.cc:37
NOTICE
#define NOTICE(A)
Definition JMessage.hh:64
FATAL
#define FATAL(A)
Definition JMessage.hh:67
debug
int debug
debug level
Definition JSirene.cc:72
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
JPARSER::JParser
Utility class to parse command line options.
Definition JParser.hh:1698
FITL1DT::removeRowColumn
JMatrixNS removeRowColumn(const JMatrixNS &A, int C)
Definition JFitL1dt.cc:26
JDETECTOR::load
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Definition JDetectorToolkit.hh:508
JDETECTOR::store
void store(const std::string &file_name, const JDetector &detector)
Store detector to output file.
Definition JDetectorToolkit.hh:605
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
JMATH::JMatrixNS
N x N symmetric matrix.
Definition JMatrixNS.hh:30
Generated by doxygen 1.12.0