#include "km3net-dataformat/online/JDAQ.hh" #include "JDetector/JDetector.hh" #include "JDetector/JDetectorToolkit.hh" #include "Jeep/JParser.hh" #include "JTools/JCombinatorics.hh" #include #include #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 "TGraph.h" #include #include #include /** * \file * * Auxiliary program to determine best inter-DU time offsets from hit time correlations * * \author dsamtleben */ int main(int argc, char **argv) { using namespace std; using namespace JPP; using namespace JMATH; string inputFile; string detectorFile; string option; double min_cont; int idu_ref; int debug; try { JParser<> zap("Program to calculate time offsets from FitL1dtSlices output"); zap['f'] = make_field(inputFile, "input file") = "nall.txt"; zap['a'] = make_field(detectorFile, "detector file"); zap['r'] = make_field(idu_ref, "reference DU set to t=0") = 0; zap['m'] = make_field(min_cont, "minimal content" ) = 5; zap['d'] = make_field(debug) = 0; if (zap.read(argc, argv) != 0) { return 1; } } catch(const exception &error) { FATAL(error.what() << endl); } JDetector detector; try { load(detectorFile, detector); } catch(const JException& error) { FATAL(error); } const int number_of_strings = getNumberOfStrings(detector); TVectorD means(1); TMatrixD variance(1,1); TMatrixD combimatrix(1,number_of_strings); int du1,du2,du1_name,du2_name,nei; int duarr[number_of_strings]; double off1,off2,offset,max1,max2; means[0]=0.; combimatrix(0,idu_ref)=1; int npairs = 1; std::ifstream file(inputFile); // --- read line with DU pair du1 & du2 and offset while (file >> du1 >> du2 >> off1 >> off2 >> offset >> max1 >> max2 >> du1_name >> du2_name >> nei) { if (max1 > min_cont and max2 > min_cont and max1 < 50. and max2 < 50.){ combimatrix.ResizeTo(npairs+1,number_of_strings); combimatrix(npairs,du1) = +1; combimatrix(npairs,du2) = -1; means.ResizeTo(npairs+1); means[npairs]=offset; duarr[du1]=du1_name; duarr[du2]=du2_name; npairs+=1; } } // Solve matrix equation: // // combimatrix: A // means: b // // A * string_time_offsets = b // string_time_offset = (A^T * A)^(-1) * (A^T * b) TMatrixD combimatrix_T(number_of_strings,npairs); combimatrix_T.Transpose(combimatrix); // A^T TMatrixD pseudo = combimatrix_T * combimatrix; // A^T * A TMatrixD pseudoinverse = pseudo.Invert(); // (A^T * A)^(-1) TVectorD offsets = pseudoinverse * combimatrix_T * means; // (A^T * A)^(-1) * (A^T * b) TVectorD residues = (combimatrix * offsets) - means; std::cout.precision(3); for (int ii=0;ii