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 "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 (int ind_neighbour=1; ind_neighbour<=neighbour ; ind_neighbour++){


  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 << ".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) == ind_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.);
          
         // 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] << " " << ind_neighbour << endl;
           }
    }
  }

  } // loop over neighbours

  TFile out(outputFile.c_str(), "recreate");

  out.cd();

  for (auto& h1 : H1) {
    for (auto& h2 : h1.second) {
      h2.second.p->Write();
    }
  }
  
  out.Close();
}