examples/JReconstruction/JMuonCompass.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <map>
 
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TF1.h"
#include "TFitResult.h"
#include "TError.h"
 
#include "JROOT/JMinimizer.hh"
 
#include "JSupport/JMeta.hh"
 
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.hh"
#include "JDetector/JStringRouter.hh"
 
#include "JReconstruction/JMuonCompass.hh"
 
#include "Jeep/JProperties.hh"
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
/**
 * \file
 * Program to fit histogram in output of JMuonCompass.cc.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  struct {
    double          XMIN   = -0.2*PI;   //!< minimal abscissa with respect to minimum [rad]
    double          XMAX   = +0.2*PI;   //!< maximal abscissa with respect to minimum [rad]
    double          FACTOR =  1.0;      //!< scaling factor for error
    double          YMIN   = 10.0;      //!< minimal difference in likelihood
  } parameters;
    
  vector<string>  inputFile;
  string          outputFile;
  string          detectorFile;
  JCounter        overwriteDetector;
  string          option;
  int             debug;
 
  try {
 
    JProperties properties;
 
    properties.insert(gmake_property(parameters.XMIN));
    properties.insert(gmake_property(parameters.XMAX));
    properties.insert(gmake_property(parameters.FACTOR));
    properties.insert(gmake_property(parameters.YMIN));
 
    JParser<> zap("Program to fit histogram in output of JMuonCompass.cc.");
 
    zap['@'] = make_field(properties)        = JPARSER::initialised();
    zap['f'] = make_field(inputFile,         "input files (output from JMuonCompass).");
    zap['o'] = make_field(outputFile,        "output file.")                                                   = "hobbit.root";
    zap['a'] = make_field(detectorFile,      "detector file.");
    zap['A'] = make_field(overwriteDetector, "overwrite detector file provided through '-a' with corrected orientations;"\
                          "\nUse option -AA to also correct compass calibrations.");
    zap['O'] = make_field(option,            "ROOT fit option, see TH1::Fit.")                                 = "";
    zap['d'] = make_field(debug)          = 1;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
 
  gErrorIgnoreLevel = kFatal;
 
 
  JDetector detector;
 
  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }
 
  const JModuleRouter router(detector);
 
  if (option.find('S') == string::npos) { option += 'S'; }
  if (option.find('Q') == string::npos && debug < JEEP::debug_t) { option += 'Q'; }
 
  TH2D* h2 = NULL;
 
  for (vector<string>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {
 
    NOTICE("Processing " << *i << endl);
 
    TFile in(i->c_str(), "exist");
 
    if (!in.IsOpen()) {
      FATAL("File " << *i << " not opened." << endl);
    }
 
    TH2D* p = dynamic_cast<TH2D*>(in.Get(h2_t));
      
    if (p == NULL) {
      FATAL("File " << *i << " has no histogram <" << h2_t << ">." << endl);
    }
 
    if (h2 == NULL)
      h2 = (TH2D*) p->Clone();
    else
      h2->Add(p);
 
    h2->SetDirectory(0);
 
    in.Close();
  }
 
  if (h2 == NULL) {
    FATAL("No histogram <" << h2_t << ">." << endl);
  }
 
  // auxiliary data structure for fit result
 
  struct result_type {
 
    result_type() :
      value(0.0),
      error(0.0)
    {}
 
    result_type(double value,
                double error) :
      value(value),
      error(error)
    {}
 
    double value;
    double error;
  };
 
  typedef map<int, result_type> map_type;    // index -> fit result
 
  map_type zmap;
 
 
  // histogram fit results
 
  const TAxis* x_axis = h2->GetXaxis();
  const TAxis* y_axis = h2->GetYaxis();
 
  TH1D h0("h0", NULL, x_axis->GetNbins(), -0.5, x_axis->GetNbins() - 0.5);
 
  const floor_range   floors = getRangeOfFloors(detector);
  const JStringRouter strings(detector);
 
  TH2D H2d("detector", NULL,
           strings.size(), -0.5, strings.size() - 0.5,
           floors.getUpperLimit(), 1 - 0.5, floors.getUpperLimit() + 0.5);
 
  for (Int_t ix = 1; ix <= H2d.GetXaxis()->GetNbins(); ++ix) {
    H2d.GetXaxis()->SetBinLabel(ix, MAKE_CSTRING(strings.at(ix-1)));
  }
  for (Int_t iy = 1; iy <= H2d.GetYaxis()->GetNbins(); ++iy) {
    H2d.GetYaxis()->SetBinLabel(iy, MAKE_CSTRING(iy));
  }
 
  TH2D* H2c = (TH2D*) H2d.Clone("chi2");
  TH2D* H2q = (TH2D*) H2d.Clone("status");
  
 
  TFile out(outputFile.c_str(), "recreate");
 
  for (Int_t ix = 1; ix <= x_axis->GetNbins(); ++ix) {
 
    const int      index  = ix - 1; 
    const JModule& module = detector[index];
 
    if (module.getFloor() == 0) {
      continue;
    }
 
    DEBUG("Bin " << setw(4) << ix << " -> " << setw(8) << module.getID() << endl);
 
    TH1D h1(MAKE_CSTRING(module.getID() << ".1D"), NULL, y_axis->GetNbins(), y_axis->GetXmin(), y_axis->GetXmax());
 
    // start values
    
    Double_t ymin   =  numeric_limits<double>::max();
    Double_t ymax   =  numeric_limits<double>::lowest();
    Double_t x0     =  0.0;       // [rad]
 
    for (Int_t iy = 1; iy <= y_axis->GetNbins(); ++iy) {
 
      h1.SetBinContent(iy, h2->GetBinContent(ix,iy));
      h1.SetBinError  (iy, h2->GetEntries() * 1.0e-7 * parameters.FACTOR);
 
      const Double_t x = h1.GetBinCenter (iy);
      const Double_t y = h1.GetBinContent(iy);
 
      if (y < ymin) {
        x0 = x;
      }
 
      if (y > ymax) { ymax = y; }
      if (y < ymin) { ymin = y; }
    }
 
    if (ymax - ymin >= parameters.YMIN) {
    
      Double_t xmin = x0 + parameters.XMIN;
      Double_t xmax = x0 + parameters.XMAX;
 
      Double_t Rx   = 0.0;
 
      if (xmin < y_axis->GetXmin()) { Rx = +PI; }
      if (xmax > y_axis->GetXmax()) { Rx = -PI; }
 
      // apply offset
 
      if (Rx != 0.0) {
 
        NOTICE("Module " << setw(8) << module.getID() << " offset " << Rx/PI << " [pi]" << endl);
 
        x0   += Rx;
        xmin += Rx;
        xmax += Rx;
 
        h1.Reset();
 
        const Double_t dx = (y_axis->GetXmax() - y_axis->GetXmin()) / y_axis->GetNbins();
 
        for (Int_t iy = 1; iy <= y_axis->GetNbins(); ++iy) {
 
          Int_t i = iy + (Int_t) (Rx / dx);
 
          if (i < 1)                  { i += (Int_t) (2.0*PI / dx); }
          if (i > y_axis->GetNbins()) { i -= (Int_t) (2.0*PI / dx); }
 
          h1.SetBinContent(i, h2->GetBinContent(ix,iy));
          h1.SetBinError  (i, h2->GetEntries() * 1.0e-7 * parameters.FACTOR);
        }
      }
 
      if (xmin < y_axis->GetXmin()) { xmin = y_axis->GetXmin();  xmax = xmin + 0.8*PI; }
      if (xmax > y_axis->GetXmax()) { xmax = y_axis->GetXmax();  xmin = xmax - 0.8*PI; }
 
      TF1* f1 = new TF1("f1", "[0] - [1]*cos(x - [2])");
 
      f1->SetParameter(0, 0.5 * (ymin + ymax));
      f1->SetParameter(1, 0.5 * (ymax - ymin));
      f1->SetParameter(2, x0);
 
      f1->SetParError(0, 0.01*(ymin+ymax));
      f1->SetParError(1, 0.01*(ymin+ymax));
      f1->SetParError(2, 0.05);
    
      DEBUG("Start values:" << endl);
      
      for (int i = 0; i != f1->GetNpar(); ++i) {
        DEBUG(left << setw(12) << f1->GetParName  (i) << ' ' << 
              SCIENTIFIC(12,5) << f1->GetParameter(i) << endl);
      }
 
      // fit
      
      TFitResultPtr result = h1.Fit(f1, option.c_str(), "same", xmin, xmax);
 
      const bool status = result->IsValid();
 
      // revert offset
      
      if (Rx != 0.0) {
        f1->SetParameter(2, f1->GetParameter(2) - Rx); 
      }
 
      if (debug >= notice_t || !result->IsValid()) {
        cout << "Slice: " 
             << setw(4)    << ix                    << ' '
             << setw(16)   << h1.GetName()          << ' '
             << FIXED(7,3) << f1->GetParameter(2)   << " +/- " 
             << FIXED(7,3) << f1->GetParError(2)    << " [rad] "
             << FIXED(9,2) << result->Chi2()        << '/'
             << result->Ndf()                       << ' '
             << result->IsValid()
             << (status ? "" : "failed")            << endl;
      }
 
      if (status) {
        zmap[index] = result_type(f1->GetParameter(2), f1->GetParError(2)); 
      }
 
      struct {
        const int ix;
        const int iy;
      } H2 = {
        strings.getIndex(module.getString()) + 1,
        module .getFloor()
      };
 
      if (status) {
        H2d .SetBinContent(H2.ix, H2.iy, f1->GetParameter(2));
      }
      if (result->Ndf() > 0) {
        H2c->SetBinContent(H2.ix, H2.iy, result->Chi2() / result->Ndf());
      }
      {
        H2q->SetBinContent(H2.ix, H2.iy, status ? +1.0 : -1.0);
      }
 
      h1.Write();
 
      delete f1;
    }
  }
 
  if (zmap.empty()) {
    NOTICE("No calibration results." << endl);
  }
 
  for (map_type::const_iterator i = zmap.begin(); i != zmap.end(); ++i) {
    h0.SetBinContent(i->first, i->second.value);
    h0.SetBinError  (i->first, i->second.error);
  }
 
  if (overwriteDetector) {
 
    NOTICE("Store calibration data on file " << detectorFile << endl);
 
    detector.comment.add(JMeta(argc, argv));
 
    for (map_type::const_iterator i = zmap.begin(); i != zmap.end(); ++i) {
 
      JModule&        module = detector[i->first];
      const JVector3D center = module.getPosition();
 
      module.sub(center);
 
      module.rotate(JRotation3Z(i->second.value));
 
      module.add(center);
 
      if (overwriteDetector > 1) {
        module.setQuaternion(module.getQuaternion() * JQuaternion3Z(i->second.value));
      }
    }
      
    store(detectorFile, detector);
  }
 
  h0  .Write();
  h2 ->Write();
  H2d .Write();  
  H2c->Write();
  H2q->Write();
 
  out.Close();
}