JOffset_extract.cc

Go to the documentation of this file.

\
#include <iostream>
#include <string>
#include <map>
 
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TF1.h"
#include "TMath.h"
 
#include "JROOT/JMinimizer.hh"
 
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JStringRouter.hh"
 
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
 
 
         /**
          * \file
          *
          * Auxiliary program to extract string-string correlation information from the output created with JMonitorL1dt
          *
          * \author dsamtleben
          */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  string   inputFile;
  string   outputFile;
  string   detectorFile;
  int      neighbour, setmax;
  int      debug;
 
  try { 
 
    JParser<> zap("Program to extract time offsets of DOM-DOM correlations");
    
    zap['f'] = make_field(inputFile,    "input file")                                        = "monitor.root";
    zap['o'] = make_field(outputFile,   "output file")                                       = "corr_max.root";
    zap['a'] = make_field(detectorFile, "detector file");
    zap['n'] = make_field(neighbour,     "neighbour level")                              = 5;
    zap['m'] = make_field(setmax,        "minimal entries")                                 = 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 JStringRouter string(detector);
 
  map<int, map<int, int> >  zmap;
  
  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
    zmap[module->getString()][module->getFloor()] = module->getID();
  }
 
  struct h1_t {
    h1_t() :
      p(NULL)  // pointer must be initialised
    {}
    
    TH1D* p;
  };
  
  TF1 f1("f1", "[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2])) + [3]");
 
  const int number_of_strings = getNumberOfStrings(detector);
 
  //  cout << "number of strings " << number_of_strings << endl;
 
  double maxarr[number_of_strings][number_of_strings][3];
  int duarr[number_of_strings];
 
  map<int, map<int, h1_t > > H1;
  
  TFile* in = TFile::Open(inputFile.c_str(), "exist");
 
  if (in == NULL || !in->IsOpen()) {
    FATAL("File: " << inputFile << " not opened." << endl);
  }
 
 
  for (const auto& string_1 : zmap) {
    duarr[string.getIndex(string_1.first)]=string_1.first;
 
    for (const auto& floor_1 : string_1.second) {
 
      const int module_1 = floor_1.second;                  // parent module
      
      TH2D* h2 = (TH2D*) in->Get(MAKE_CSTRING(module_1 << ".2S"));
 
      if (h2 == NULL) {
        continue;
      }
      
      for (const auto& string_2 : zmap) {
        duarr[string.getIndex(string_2.first)]=string_2.first;
 
        if (string_1.first != string_2.first) {
 
          TH1D* h1 = H1[string_1.first][string_2.first].p;
          
          if (h1 == NULL) {
 
            h1 = new TH1D(MAKE_CSTRING(string_1.first << "_" << string_2.first << "_" << neighbour << ".2T"), NULL,
                          h2->GetYaxis()->GetNbins(), h2->GetYaxis()->GetXmin(), h2->GetYaxis()->GetXmax());
 
            H1[string_1.first][string_2.first].p = h1;      // book keeping
          }
 
          for (const auto& floor_2 : string_2.second) {
            if (floor_1.first > floor_2.first && (floor_1.first - floor_2.first) == neighbour) {            // floor condition
 
              const int module_2 = floor_2.second;          // daughter module
 
              TH1D* py = h2->ProjectionY("__py",
                                         h2->GetXaxis()->FindBin(TString(Form("%i", module_2))), 
                                         h2->GetXaxis()->FindBin(TString(Form("%i", module_2))), "e");
              
              h1->Add(py);                                  // add data
 
              delete py;
            }
          }
          //------------
          if (floor_1.first == 18){
 
            double mm=h1->GetXaxis()->GetBinCenter(h1->GetMaximumBin());
         
            f1.SetParameter(0, double(h1->GetMaximum()));
            f1.SetParameter(1, double(mm));
            f1.SetParameter(2, 50.);
 
            for (Int_t i = 0; i != f1.GetNpar(); ++i) {
              f1.SetParError(i, 0.0);
            }
          
            // fit
    
            if (h1->GetMaximum()>0){
              //                cout << "we are fitting " << string_1.first << " " << string_2.first << endl;
              h1->Fit(&f1,"LQ", "same");
 
              //            maxarr[ string.getIndex(string_1.first) ][ string.getIndex(string_2.first) ][0] = h1->GetXaxis()->GetBinCenter(h1->GetMaximumBin());
              //            maxarr[ string.getIndex(string_1.first) ][ string.getIndex(string_2.first) ][1] = h1->GetMaximum();
              maxarr[ string.getIndex(string_1.first) ][ string.getIndex(string_2.first) ][0] = h1->GetFunction("f1")->GetParameter(1);
              maxarr[ string.getIndex(string_1.first) ][ string.getIndex(string_2.first) ][1] = h1->GetFunction("f1")->GetParameter(0);
              maxarr[ string.getIndex(string_1.first) ][ string.getIndex(string_2.first) ][2] = h1->GetFunction("f1")->GetParameter(2);
            }
          }
        } // if string 1 != string 2
      } // loop over string 2
    } // loop over modules of string 1
  } // loop over string 1
 
 
  for (int i=0;i<number_of_strings;i++){
    for (int j=i+1;j<number_of_strings;j++){
      if (maxarr[i][j][1]>setmax && maxarr[j][i][1]>setmax && maxarr[i][j][2]>20. && maxarr[j][i][2]>20.){
        cout << i << " " << j << " " << maxarr[i][j][0] << " " << maxarr[j][i][0] << " " << (maxarr[i][j][0]-maxarr[j][i][0])/2. << " " << maxarr[i][j][1] << " " << maxarr[j][i][1] << " " << " " << duarr[i] << " " << duarr[j] << " " << neighbour << endl;
      }
    }
  }
 
 
 
  TFile out(outputFile.c_str(), "recreate");
 
  out.cd();
 
  for (auto& h1 : H1) {
    for (auto& h2 : h1.second) {
      h2.second.p->Write();
    }
  }
  
  out.Close();
}