//  
//  flukaExtract.C
//  macro to extract the information from the FLUKA raw data
//  and store it in histograms
//
//  Created by Christian Elsasser on 31.05.13.
//  University of Zurich, elsasser@cern.ch
//  Copyright only for commercial use
//

//  General include
#include <iostream>
#include <fstream>
#include <string.h>
#include <vector>
#include <map>
#include <locale>


//  ROOT include
#include "../../include/incBasicROOT.h"
#include "../../include/incDrawROOT.h"
#include "../../include/incHistROOT.h"
#include "../../include/incIOROOT.h"




#include "../../include/basicROOTTools.h"
#include "../../include/basicRooFitTools.h"

#include "../../include/lhcbstyle.C"








int flukaExtract(int, char**);

int writeHisto(TString, TH2F*);
int setUncertainty(TH2F*);
int scale(TH2F*, TH2F*);
int merge(TH2F*, TH2F*);
int storeUncertainty(TH2F*);
int configProfile(TH1D*);

int main(int argc, char* argv[]){
  
  TStyle* style = lhcbStyle();
  
  int argc_copy = argc;
  
  TApplication theApp("Analysis", &argc, argv);
  argc = theApp.Argc();
  argv = theApp.Argv();
  
  style->SetPadLeftMargin(0.12);
  style->SetPadRightMargin(0.22);
  style->SetPadBottomMargin(0.15);
  style->SetPadTopMargin(0.05);
  
  int exit = flukaExtract(argc,argv);
  
  Printf("End");
  theApp.Run(kTRUE); 
  
  return exit;
}


//  Executable method
int flukaExtract(int argc, char* argv[]){
  
  //const char* homeVar = std::getenv ("HOME");
  const char* homeVar = "/afs/cern.ch/user/e/egraveri/cmtuser/STMonitoring/STAging";

  TString dn_lumi(Form("%s/data/fluka",homeVar));
  
  TString fn_output(Form("%s/fluka.root",
                         dn_lumi.Data()));
  
  
  TString fn_TT_old(Form("%s/TTSistrip_%s_Si1MeVne_cm-2_per_coll.txt",
                       dn_lumi.Data(),
                       "14TeV"));
  TString fn_T3_old(Form("%s/IT3_%s_Si1MeVne_cm-2_per_coll.txt",
                        dn_lumi.Data(),
                        "14TeV"));
  
  
  TString fn_TT_8(Form("%s/TTSistrip_%s_Si1MeVne_cm-2_per_coll.txt",
                       dn_lumi.Data(),
                       "8TeV"));
  TString fn_T3A_8(Form("%s/IT3As_%s_Si1MeVne_cm-2_per_coll.txt",
                       dn_lumi.Data(),
                       "8TeV"));
  TString fn_T3B_8(Form("%s/IT3bottom_%s_Si1MeVne_cm-2_per_coll.txt",
                       dn_lumi.Data(),
                       "8TeV"));
  TString fn_T3C_8(Form("%s/IT3Cs_%s_Si1MeVne_cm-2_per_coll.txt",
                       dn_lumi.Data(),
                       "8TeV"));
  TString fn_T3T_8(Form("%s/IT3top_%s_Si1MeVne_cm-2_per_coll.txt",
                       dn_lumi.Data(),
                       "8TeV"));
  
  TString fn_TT_7(Form("%s/TTSistrip_%s_Si1MeVne_cm-2_per_coll.txt",
                       dn_lumi.Data(),
                       "7TeV"));
  TString fn_T3A_7(Form("%s/IT3As_%s_Si1MeVne_cm-2_per_coll.txt",
                        dn_lumi.Data(),
                        "7TeV"));
  TString fn_T3B_7(Form("%s/IT3bottom_%s_Si1MeVne_cm-2_per_coll.txt",
                        dn_lumi.Data(),
                        "7TeV"));
  TString fn_T3C_7(Form("%s/IT3Cs_%s_Si1MeVne_cm-2_per_coll.txt",
                        dn_lumi.Data(),
                        "7TeV"));
  TString fn_T3T_7(Form("%s/IT3top_%s_Si1MeVne_cm-2_per_coll.txt",
                        dn_lumi.Data(),
                        "7TeV"));

  
  
  double minRange = -60.0;
  double maxRange =  60.0;
  double binRange =   1.0;
  int nBins = 120;
  
  TH2F* h_tt_old = new TH2F("h_tt_old","Fluka map",
                            nBins,minRange,maxRange,
                            nBins,minRange,maxRange);
  h_tt_old->SetXTitle("#it{x} position [cm]");
  h_tt_old->SetYTitle("#it{y} position [cm]");
  h_tt_old->SetZTitle("#phi_{1 MeV-n,eq} [cm^{-2} per coll]");
  h_tt_old->GetZaxis()->SetTitleOffset(1.30);
  h_tt_old->SetMinimum(1e-4);
  h_tt_old->SetMaximum(1e+0);
  
  TH2F* h_tt_7     = (TH2F*)h_tt_old->Clone("h_tt_7");
  TH2F* h_tt_8     = (TH2F*)h_tt_old->Clone("h_tt_8");
  TH2F* h_it_old   = (TH2F*)h_tt_old->Clone("h_it_old");
  TH2F* h_it_7     = (TH2F*)h_tt_7->Clone("h_it_7");
  TH2F* h_it_8     = (TH2F*)h_tt_8->Clone("h_it_8");
  
  writeHisto(fn_TT_old,h_tt_old);
  writeHisto(fn_T3_old,h_it_old);
  
  
  
  writeHisto(fn_TT_8, h_tt_8);
  writeHisto(fn_T3A_8,h_it_8);
  writeHisto(fn_T3B_8,h_it_8);
  writeHisto(fn_T3C_8,h_it_8);
  writeHisto(fn_T3T_8,h_it_8);
  
  writeHisto(fn_TT_7, h_tt_7);
  writeHisto(fn_T3A_7,h_it_7);
  writeHisto(fn_T3B_7,h_it_7);
  writeHisto(fn_T3C_7,h_it_7);
  writeHisto(fn_T3T_7,h_it_7);
  
  // Mirror old one
  int nBinsX = h_tt_old->GetNbinsX();
  int nBinsY = h_tt_old->GetNbinsY();
  for (int i=nBinsX/2+1; i<=nBinsX; i++) {
    for (int j=nBinsY/2+1; j<=nBinsY; j++) {
      double tt_val = h_tt_old->GetBinContent(i,j);
      double it_val = h_it_old->GetBinContent(i,j);
      h_tt_old->SetBinContent(nBinsX-i+1,j,tt_val);
      h_tt_old->SetBinContent(nBinsX-i+1,nBinsY-j+1,tt_val);
      h_tt_old->SetBinContent(i,nBinsY-j+1,tt_val);
      h_it_old->SetBinContent(nBinsX-i+1,j,it_val);
      h_it_old->SetBinContent(nBinsX-i+1,nBinsY-j+1,it_val);
      h_it_old->SetBinContent(i,nBinsY-j+1,it_val);
    }
  }
  
  setUncertainty(h_tt_8);
  setUncertainty(h_tt_7);
  setUncertainty(h_tt_old);
  setUncertainty(h_it_8);
  setUncertainty(h_it_7);
  setUncertainty(h_it_old);
  
  TH2F* h_tt_old_8 = (TH2F*)h_tt_old->Clone("h_tt_old_8");
  TH2F* h_tt_old_7 = (TH2F*)h_tt_old->Clone("h_tt_old_7");
  TH2F* h_it_old_8 = (TH2F*)h_it_old->Clone("h_it_old_8");
  TH2F* h_it_old_7 = (TH2F*)h_it_old->Clone("h_it_old_7");
  
  scale(h_tt_8,h_tt_old_8);
  scale(h_it_8,h_it_old_8);
  scale(h_tt_7,h_tt_old_7);
  scale(h_it_7,h_it_old_7);
  
  merge(h_tt_8,h_tt_old_8);
  merge(h_tt_7,h_tt_old_7);
  
  storeUncertainty(h_tt_8);
  storeUncertainty(h_tt_7);
  storeUncertainty(h_it_8);
  storeUncertainty(h_it_7);
  
  
  TH1D* h_xProf_7 = h_tt_7->ProjectionX("h_xProf_7",
                                         h_tt_7->GetNbinsY()/2,
                                         h_tt_7->GetNbinsY()/2);
  TH1D* h_yProf_7 = h_tt_7->ProjectionY("h_yProf_7",
                                         h_tt_7->GetNbinsX()/2,
                                         h_tt_7->GetNbinsX()/2);
  
  TH1D* h_xProf_8 = h_tt_8->ProjectionX("h_xProf_8",
                                         h_tt_8->GetNbinsY()/2,
                                         h_tt_8->GetNbinsY()/2);
  TH1D* h_yProf_8 = h_tt_8->ProjectionY("h_yProf_8",
                                         h_tt_8->GetNbinsX()/2,
                                         h_tt_8->GetNbinsX()/2);
  
  configProfile(h_xProf_7);
  configProfile(h_yProf_7);
  configProfile(h_xProf_8);
  configProfile(h_yProf_8);
  
  
  TCanvas* c_tt_8 = new TCanvas("c_tt_8","c_tt_8",1000,800);
  c_tt_8->SetLogz();
  h_tt_8->Draw("colz");
  
  TCanvas* c_tt_7 = new TCanvas("c_tt_7","c_tt_7",1000,800);
  c_tt_7->SetLogz();
  h_tt_7->Draw("colz");
  
  TCanvas* c_it_8 = new TCanvas("c_it_8","c_it_8",1000,800);
  c_it_8->SetLogz();
  h_it_8->Draw("colz");
  
  TCanvas* c_it_7 = new TCanvas("c_it_7","c_it_7",1000,800);
  c_it_7->SetLogz();
  h_it_7->Draw("colz");
  
  
  gStyle->SetPadLeftMargin(0.20);
  gStyle->SetPadRightMargin(0.10);
  gStyle->SetPadBottomMargin(0.15);
  gStyle->SetPadTopMargin(0.05);
  
  
  
  
  TCanvas* c_tt_pX7 = new TCanvas("c_tt_pX7","c_tt_pX7",1000,800);
  c_tt_pX7->SetLogy();
  h_yProf_7->SetLineColor(kBlue);
  h_yProf_7->SetMarkerColor(kBlue);
  h_yProf_7->SetMarkerStyle(25);
  h_xProf_7->SetXTitle("#it{x}/#it{y} position [cm]");
  h_xProf_7->Draw("H");
  h_yProf_7->Draw("H same");
  
  
  TCanvas* c_tt_pX8 = new TCanvas("c_tt_pX8","c_tt_pX8",1000,800);
  c_tt_pX8->SetLogy();
  h_yProf_8->SetLineColor(kBlue);
  h_yProf_8->SetMarkerColor(kBlue);
  h_yProf_8->SetMarkerStyle(25);
  h_xProf_8->SetXTitle("#it{x}/#it{y} position [cm]");
  h_xProf_8->Draw("H");
  h_yProf_8->Draw("H same");
  
  TCanvas* c_tt_pY7 = new TCanvas("c_tt_pY7","c_tt_pY7",1000,800);
  c_tt_pY7->SetLogy();
  h_yProf_7->Draw("H");
  
  TCanvas* c_tt_pY8 = new TCanvas("c_tt_pY8","c_tt_pY8",1000,800);
  c_tt_pY8->SetLogy();
  h_yProf_8->Draw("H");
  
  
  TFile* f_output = new TFile(fn_output.Data(),"RECREATE");
  h_tt_7->Write("h_tt_7");
  h_tt_8->Write("h_tt_8");
  
  h_it_7->Write("h_t3_7");
  h_it_8->Write("h_t3_8");
  
  f_output->Close();
  
  
  TLegend* leg = new TLegend(0.0,0.0,1.0,1.0);
  leg->SetFillStyle(0);
  leg->SetTextSize(0.3);
  leg->AddEntry(h_xProf_7,"#it{y} = 0 cm","p");
  leg->AddEntry(h_yProf_7,"#it{x} = 0 cm","p");
  
  TCanvas* c_leg = new TCanvas("c_leg","c_leg",150,150);
  c_leg->cd();
  leg->Draw();
  
  
  TH1F* h_func = new TH1F("h_func","h_func",100,0.0,100.0);
  h_func->SetXTitle("#it{r} [cm]");
  h_func->SetYTitle("#sigma_{#phi_{1 MeV-n,eq}}/#phi_{1 MeV-n,eq}");
  h_func->GetYaxis()->SetTitleOffset(1.3);
  h_func->SetMaximum(1.0e-0);
  h_func->SetMinimum(1.0e-3);
  
  double alpha = 0.918;
  double r1    = 0.7e-2;
  
  TF1* f_relunc = new TF1("f_relunc","[0]*TMath::Power(x,[1])",0.0,200.0);
  f_relunc->SetParameters(r1,alpha);
  f_relunc->SetLineColor(kRed);
  
  TCanvas* c_relunc = new TCanvas("c_relunc","c_relunc",1000,800);
  c_relunc->SetLogy();
  h_func->Draw("");
  f_relunc->Draw("Same");
  
  
  return 0;
                               
}







int writeHisto(TString f_name,TH2F* h_write){
  std::ifstream f_input(f_name.Data());
  
  std::vector<std::vector<std::string> > a_input;
  double binWidthX = 0.0;
  double binWidthY = 0.0;
  while (f_input) {
    
    
    std::string line;
    if (!getline(f_input,line)) {
      Error("flukaExtract","Problem in line! Leave read-in!");
      break;
    }
    
    
    std::string binSizeTag = "* Bin size: ";
    
    if (line[0]=='*') {
      std::size_t p = line.find(binSizeTag);
      if (p!=std::string::npos) {
        line.replace(line.find(binSizeTag),binSizeTag.length(),"");
        std::istringstream s_line(line);
        std::string bit;
        if (!getline(s_line,bit,',')) {
          Error("flukaExtract","Bin width X not available!");
          return 1;
        }
        binWidthX = atof(bit.c_str());
        if (!getline(s_line,bit,',')) {
          Error("flukaExtract","Bin width Y not available!");
          return 1;
        }
        binWidthY = atof(bit.c_str());
        
        Info("flukaExtract","Bin size (%4.2f,%4.2f)",binWidthX,binWidthY);
        continue;
      }
      
      Info("flukaExtract","Header line ignoring!");
      continue;
    }
    
    std::istringstream s_line(line);
    std::vector<std::string> v_line;
    double x,y,flux;
    for (int i=0; i<3 && s_line; i++) {
      std::string bit;
      if (!getline(s_line,bit,'\t')) {
        break;
      }
      if (bit.compare("")==0) {
        break;
      }
      if (i==0) {
        x = atof(bit.c_str());
      }else if(i==1){
        y = atof(bit.c_str());
      }else if(i==2){
        flux = atof(bit.c_str());
        for (int xcoord=0; 1.0*xcoord<binWidthX; xcoord++) {
          for (int ycoord=0; 1.0*ycoord<binWidthY; ycoord++) {
            Info("flukaExtract","(%4.2f,%4.2f): %4.4f",
                 x+0.5+1.0*xcoord,y+0.5+1.0*ycoord,flux);
            int bin = h_write->FindBin(x+0.5+1.0*xcoord,y+0.5+1.0*ycoord);
            h_write->SetBinContent(bin,flux);
          }
        }
      }
      
      
    }
    
  }
  return 0;
}

// Determines uncertainty on fluence
int setUncertainty(TH2F* h_write){
  // Determination of uncertainty
  double alpha = 0.918;
  double r1    = 0.7e-2;
  for (int i=1; i<=h_write->GetNbinsX(); i++) {
    for (int j=1; j<=h_write->GetNbinsY(); j++) {
      double x = h_write->GetXaxis()->GetBinCenter(i);
      double y = h_write->GetYaxis()->GetBinCenter(j);
      double r = TMath::Sqrt(x*x+y*y);
      double u = r1*TMath::Power(r,alpha);
      h_write->SetBinError(i,j,u);
    }
  }
  return 0;
}

// Scales old FLUKA map by chi^2 fit
int scale(TH2F* h_ref, TH2F* h_tar){
  if (h_tar->GetNbinsX()!=h_ref->GetNbinsX() ||
      h_tar->GetNbinsY()!=h_ref->GetNbinsY()) {
    Error("scale","Histograms do not have the same size!");
    return 1;
  }
  double xy = 0.0, xx = 0.0;
  for (int i=1; i<=h_tar->GetNbinsX(); i++) {
    for (int j=1; j<=h_tar->GetNbinsY(); j++) {
      if (TMath::Abs(h_tar->GetBinContent(i,j))<1e-9 ||
          TMath::Abs(h_ref->GetBinContent(i,j))<1e-9) {
        continue;
      }
      xy+=h_tar->GetBinContent(i,j)/(TMath::Power(h_ref->GetBinError(i,j),2)*h_ref->GetBinContent(i,j));
      xx+=TMath::Power(h_tar->GetBinContent(i,j),2)/(TMath::Power(h_ref->GetBinError(i,j)*h_ref->GetBinContent(i,j),2));
    }
  }
  double alpha = xy/xx;
  double chi2  = 0.0;
  int ndf = 0;
  for (int i=1; i<=h_tar->GetNbinsX(); i++) {
    for (int j=1; j<=h_tar->GetNbinsY(); j++) {
      h_tar->SetBinContent(i,j,alpha*h_tar->GetBinContent(i,j));
      if (TMath::Abs(h_tar->GetBinContent(i,j))<1e-9 ||
          TMath::Abs(h_ref->GetBinContent(i,j))<1e-9) {
        continue;
      }
      chi2 += TMath::Power((h_tar->GetBinContent(i,j)-h_ref->GetBinContent(i,j))/
                           (h_ref->GetBinContent(i,j)*h_ref->GetBinError(i,j)),2);
      ndf++;
    }
  }
  Info("scale","alpha = %4.2f",alpha);
  Info("scale","chi^2 = %4.2f/%d",chi2,(ndf-1));
  return 0;
}

// Adds values for empty bins from old FLUKA map
int merge(TH2F* h_tar, TH2F* h_sou){
  if (h_tar->GetNbinsX()!=h_sou->GetNbinsX() ||
      h_tar->GetNbinsY()!=h_sou->GetNbinsY()) {
    Error("merge","Histograms do not have the same size!");
    return 1;
  }
  for (int i=1; i<=h_tar->GetNbinsX(); i++) {
    for (int j=1; j<=h_tar->GetNbinsY(); j++) {
      if (TMath::Abs(h_tar->GetBinContent(i,j))>1e-9) {
        continue;
      }
      h_tar->SetBinContent(i,j,h_sou->GetBinContent(i,j));
      h_tar->SetBinError(i,j,h_sou->GetBinError(i,j));
    }
  }
  return 0;
}

// Computes absolute uncertainty from relative uncertainty
int storeUncertainty(TH2F* h){
  for (int i=1; i<=h->GetNbinsX(); i++) {
    for (int j=1; j<=h->GetNbinsY(); j++) {
      h->SetBinError(i,j,h->GetBinError(i,j)*h->GetBinContent(i,j));
    }
  }
  return 0;
}

// Configure Profiles
int configProfile(TH1D* h){
  h->SetYTitle("#phi_{1 MeV-n,eq} [cm^{-2} per coll]");
  h->SetLineColor(kRed);
  h->SetMarkerColor(kRed);
  h->GetYaxis()->SetTitleOffset(1.3);
  h->SetMaximum(1.e+0);
  h->SetMinimum(1.e-4);
  return 0;
}