JMatrix_sol.cc

Go to the documentation of this file.

#include "km3net-dataformat/online/JDAQ.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "Jeep/JParser.hh"
#include "JTools/JCombinatorics.hh"
 
#include <TVectorD.h>
#include <TMatrixD.h>
 
#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 <string>
#include <iostream>
#include <set>
 
 
/**
 * \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<number_of_strings;ii++){
    cout << duarr[ii] << " " << offsets[ii] << endl;
  }
 
  double sumres=0.;
  for (int ii=0;ii<npairs;ii++){
    sumres+=fabs(residues[ii]);
  }
  cout << -1 << " " << sumres/double(npairs) << endl;
 
  return 1;
 
 
}