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;


}