#define CHECK_DIST
#include <iomanip>

using namespace std;

const Int_t kNlayer = 28;

void PlotGMResolAllLayer(Int_t module = 2)
{
//TFile *fIn = new TFile("ntuple.root");

  gROOT->LoadMacro("SetGlobalStyle.C");
  SetGlobalStyle();

  gStyle->SetHistMinimumZero(kTRUE);
  gStyle->SetStatFontSize(0.05);
  gStyle->SetHistTopMargin(0.40);

  /////

  Double_t xdata[2][kNlayer], ydata[2][kNlayer], dxdata[2][kNlayer], dydata[2][kNlayer];
  Double_t x, y, dx, dy;
  Int_t n = 0;

  const Char_t *kInOut[2] = { "in", "ot" };

  for (Int_t inout = 0; inout < 2; inout++) {
    n = 0;
    for (Int_t ilayer = 0; ilayer < kNlayer; ilayer++) {
      Process(module, ilayer, x, y, dx, dy, kInOut[inout]);
      xdata[inout][n] = x; ydata[inout][n] = y; dxdata[inout][n] = dx; dydata[inout][n] = dy; n++; 
    }
  }

  Double_t xpoint[kNlayer], ypoint[kNlayer], dxpoint[kNlayer], dypoint[kNlayer];
  for (Int_t ilayer = 0; ilayer < kNlayer; ilayer++) {
    xpoint [ilayer] = xdata[0][ilayer];
    dxpoint[ilayer] = dxdata[0][ilayer];
#if 1
    cerr << "x [layer" << kNlayer * (module - 1) + ilayer + 1 << "] = " << xpoint [ilayer] << endl;
    cerr << "dx[layer" << kNlayer * (module - 1) + ilayer + 1 << "] = " << dxpoint[ilayer] << endl;
#endif
    // Calculate geometric mean of the point resolution
    // ydata[0][ilayer]: sigma for dx_in
    // ydata[1][ilayer]: sigma for dx_out
    ypoint [ilayer] = TMath::Sqrt(ydata[0][ilayer] * ydata[1][ilayer]);
    dypoint[ilayer] = (dydata[0][ilayer] + dydata[1][ilayer]) / 2.;
#if 1
    cerr << "y [layer" << kNlayer * (module - 1) + ilayer + 1 << "] = " << ypoint [ilayer] << endl;
    cerr << "yx[layer" << kNlayer * (module - 1) + ilayer + 1 << "] = " << dypoint[ilayer] << endl;
#endif
  }

  TCanvas *c1 = new TCanvas("c1", "Canvas 1",   0,   0, 1200, 400);
  c1->cd();

  TGraphErrors *gr = new TGraphErrors(n, xpoint, ypoint, dxpoint, dypoint);
  gr->GetHistogram()->GetXaxis()->SetLimits(0., 29.);
//gr->GetHistogram()->GetYaxis()->SetLimits(0., 0.15);
  gr->GetHistogram()->SetMinimum(0.);
  gr->GetHistogram()->SetMaximum(0.15);
  gr->GetHistogram()->GetXaxis()->SetTitle("# layer");
  gr->GetHistogram()->GetYaxis()->SetTitle("#sigma_{x} [mm]");
  gr->Draw("ap");


  TCanvas *c2 = new TCanvas("c2", "Canvas 2",   0,   0, 1200, 400);
  TCanvas *c3 = new TCanvas("c3", "Canvas 3",   0,   0, 1200, 400);

  c2->Divide(7, 4);
  c3->Divide(7, 4);


#if 0 // doesn't work
  for (Int_t ilayer = 0; ilayer < kNlayer; ilayer++) {
    c2->cd(ilayer + 1);
    stringstream hnamein;
    hnamein << "hdxin" << setw(2) << setfill('0') << kNlayer * (module - 1) + ilayer << ends;
    TH1F *hinPtr = static_cast<TH1F *>(gROOT->FindObject(hnamein.str().data()));
    hinPtr->Draw(ep);

    c3->cd(ilayer + 1);
    stringstream hnameot;
    hnamein << "hdxot" << setw(2) << setfill('0') << kNlayer * (module - 1) + ilayer << ends;
    TH1F *hotPtr = static_cast<TH1F *>(gROOT->FindObject(hnameot.str().data()));
    hotPtr->Draw(ep);
  }
#endif
}

void Process(Int_t module, Int_t layer, Double_t &x, Double_t &y, Double_t &dx, Double_t &dy, const Char_t *inout)
{
  TFile *fIn = new TFile("ntuple.root");

  stringstream item;
  item << "dx" << inout << setw(2) << setfill('0') << kNlayer * (module - 1) + layer << ends;

  stringstream cut;

  ///// Cuts ////////////////////////
  const Double_t kResXmin   = -0.3;
  const Double_t kResXmax   = +0.3;
  const Int_t    kNdfCut    =  40;
  const Double_t kChi2Cut   = 10000.;
  const Double_t kCpaMinCut = -0.8;
  const Double_t kCpaMaxCut =  0.0;
  ///////////////////////////////////

  cut << item.str() << ">" << kResXmin << "&&"
      << item.str() << "<" << kResXmax << "&&"
      << "ndf>"  << kNdfCut    << "&&"
      << "chi2<" << kChi2Cut   << ends;

  cerr << item.str().data() << endl;
  cerr <<  cut.str().data() << endl;

  TDirectory *last = gDirectory;
  stringstream hstr;
  hstr << "c" << "_" << layer << "_" << inout << ends;
  TCanvas *cp = new TCanvas(hstr.str().data(), "", 400, 400);
  cp->cd();

  hResXin->Draw(item.str().data(), cut.str().data());
  htemp->Rebin(1);
  htemp->Fit("gaus", "", "", -0.3, 0.3);
  Double_t mu = htemp->GetFunction("gaus")->GetParameter(1);
  Double_t sg = htemp->GetFunction("gaus")->GetParameter(2);
  htemp->Fit("gaus", "", "", mu - 2.5*sg, mu + 2.5*sg);

  stringstream hname;

  hname << "h" << item.str().data() << ends;

  TH1F *htemp2 = new TH1F();
  htemp->Copy(*htemp2);
  htemp2->SetName(hname.str().data());

  x  = layer + 1;
  dx = 0.;
  y  = htemp->GetFunction("gaus")->GetParameter(2);
  dy = htemp->GetFunction("gaus")->GetParError(2);
  y  *= 10.; // cm -> mm
  dy *= 10.; // cm -> mm
#ifdef CHECK_DIST
  htemp->SetMinimum(-0.1);
  htemp->SetMarkerStyle(4);
  htemp->SetMarkerSize(0.7);
  htemp->SetMarkerColor(2);
  htemp->SetLineColor(2);
  htemp->Draw("pe");
#endif
}