//  
//  fluenceExtract.C
//  macro to determine the fluence in all the read-out sectors from the FLUKA map and the sensor geometries
//
//  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>
#include <stdlib.h>

//  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"


// Remember to execute the method with all the possible options:
// ./fluenceExtract -d TT -b && ./fluenceExtract -d TT -r 40 -b && ./fluenceExtract -d TT -r 45 -b && ./fluenceExtract -d TT -r 75 -b
// ./fluenceExtract -d IT -b && ./fluenceExtract -d IT -r 40 -b && ./fluenceExtract -d IT -r 175 -b


// -d detector (TT/IT)
// -r radius


int fluenceExtract(int, char**);

int main(int argc, char* argv[]){
  
  TStyle* style = lhcbStyle();
  
  int argc_copy = argc;
  
  TApplication* theApp = new TApplication("Analysis", &argc, argv,0,-1);
  argc = theApp->Argc();
  argv = theApp->Argv();

  bool runBatch = false;
  
  for (int i = 1; i<argc; i++) {
    if (strcmp(argv[i],"-b")==0) {
      runBatch = true;
      gROOT->SetBatch();
    }
  }
  
  style->SetPadLeftMargin(0.12);
  style->SetPadRightMargin(0.22);
  style->SetPadBottomMargin(0.15);
  style->SetPadTopMargin(0.05);
  
  int exit = fluenceExtract(argc,argv);
  
  Printf("End");
  if (!runBatch) {
    theApp->Run(!runBatch);
  }
  
  return exit;
}

//  Executable method
int fluenceExtract(int argc, char* argv[]){
  
  double radius  = -10.0;
  double radmin = 0.0;
  double radmax = 2000.0;
  bool useRadius = false;
  
  TString det("TT");
  for (int i=1; i<argc; i++) {
    // Switch to IT if set from the command line
    if (strcmp(argv[i],"-d")==0 && i+1<argc) {
      det = argv[++i];
    }else if (strcmp(argv[i],"-r")==0 && i+1<argc) {
      useRadius = true;
      radius = std::atof(argv[++i]);
    }
  }

  // Set the radius to an irrelevant value
  if (radius<0.0) {
    radius = 2000.0;
  } else {
    // Elena - implement circular crowns
    if (radius == 45) {
        radmax = 45.0;
    } else if (radius == 75) {
        radmin = 45.0;
        radmax = 75.0;
    } else if (radius == 175) {
        // This setting is the only one for the IT
        radmin = 0.0;
        radmax = 175.0;
    } else if (radius == 40) {
        // Placeholder for "rest of the sector"
        // r < 40 mm has too little statistics anyways!
        radmin = 75.0;
        if (strcmp(det.Data(),"IT")==0) {
            // The only radius for the IT is 175 mm
            radmin = 175.0;
        }
        radmax = 2000.0;
    }
  }

  
  TString lay;
  if (strcmp(det.Data(),"TT")==0) {
    lay = "TTaU";
  }else{
    lay = "T3X2";
  }
  
  
  Info("fluenceExtract","Extract fluence for %s / %s",lay.Data(),det.Data());
  if (useRadius) {
    Info("fluenceExtract","Radius restriction to %4.2f mm < r < %4.2f mm", radmin, radmax);
  }
  
  
  const char* homeVar = std::getenv ("CCEHOME");
  const char* diskVar = std::getenv ("DISK");
  TString dn_fluka(Form("%s/data/fluka",homeVar));
  TString fn_fluka(Form("%s/fluka.root",
                        dn_fluka.Data()));
  
  TString dn_hit(Form("%s/data/ST/Aging",diskVar));
  TString fn_hit;
  
  if (strcmp(det.Data(),"TT")==0) {
    fn_hit = Form("%s/TTHits.root",
                  dn_hit.Data());
  }else{
    fn_hit = Form("%s/T3Hits.root",
                  dn_hit.Data());
  }
  
  TString tn_hit(Form("STHitIntegrator/%s",lay.Data()));

  TString dn_lumi(Form("%s/data/lumi",homeVar));
  TString fn_lumi(Form("%s/lumi.root",
                       dn_lumi.Data()));
  
  TString dn_out(Form("%s/data/db",homeVar));
  TString fn_out(Form("%s/fluence.root",
                      dn_out.Data()));
  
  TString fn_db(Form("%s/%s_db.txt",
                      dn_out.Data(),
                      det.Data()));
  
  Printf("%s",fn_db.Data());
  
  std::ifstream f_db(fn_db.Data());
  
  std::vector<double> v_sector;
  std::vector<double> v_flux7;
  std::vector<double> v_flux8;
  std::vector<double> v_flux13;
  std::vector<double> v_flux7_e;
  std::vector<double> v_flux8_e;
  std::vector<double> v_flux13_e;
  
  while (f_db) {
    std::string line;
    if (!getline(f_db,line)) {
      break;
    }
    std::istringstream s_line(line);
    std::vector<std::string> v_line;
    while (s_line) {
      std::string bit;
      if (!getline(s_line,bit,'\t')) {
        break;
      }
      v_line.push_back(bit);
    }
    if (v_line.size()>4 && strcmp(v_line[4].c_str(),lay.Data())==0) {
      Info("fluenceExtract","Using sector %s",v_line[0].c_str());
      v_sector.push_back(atof(v_line[0].c_str()));
    }
  }
  
  
  TFile* f_lumi = TFile::Open(fn_fluka.Data());
  
  TH2F* h_map7;
  TH2F* h_map8;
  TH2F* h_map13;
  
  if (strcmp(det.Data(),"TT")==0) {
    h_map7 = (TH2F*)f_lumi->Get("h_tt_7");
    h_map8 = (TH2F*)f_lumi->Get("h_tt_8");
    h_map13 = (TH2F*)f_lumi->Get("h_tt_13");
  }else{
    h_map7 = (TH2F*)f_lumi->Get("h_t3_7");
    h_map8 = (TH2F*)f_lumi->Get("h_t3_8");
    h_map13 = (TH2F*)f_lumi->Get("h_t3_13");
  }
  
  TFile* f_hit = TFile::Open(fn_hit.Data());
  TTree* t_hit = (TTree*)f_hit->Get(tn_hit.Data());
  
  
  
  for (int i=0; i<v_sector.size(); i++) {
    int sec = (double)v_sector[i];
    
    TH2F* h_hit7 = (TH2F*)h_map7->Clone("h_hit7");
    h_hit7->Reset("ICES");
    
    TCut c_hit(Form("(Sec==%d)",sec));
    
    if (useRadius) {
      // Elena - implement circular crowns
      c_hit=c_hit&&Form("(TMath::Sqrt(x*x+y*y)>%f/10 && (TMath::Sqrt(x*x+y*y)<%f/10))", radmin, radmax);
    }
    
    t_hit->Draw("y:x>>h_hit7",c_hit,"goff");
    h_hit7->Scale(1.0/h_hit7->Integral());
    
  
    TH2F* h_hit8 = (TH2F*)h_hit7->Clone("h_hit8");
    TH2F* h_hit13 = (TH2F*)h_hit7->Clone("h_hit13");
    h_hit7->Multiply(h_map7);
    h_hit8->Multiply(h_map8);
    h_hit13->Multiply(h_map13);
    
    if (h_hit7->Integral()<1e-9 && h_hit8->Integral()<1e-9 && h_hit13->Integral()<1e-9) {
      Info("fluenceExtract","Negligible flux for sector %d",sec);
    }
   
    
    
    double flux7_v, flux8_v, flux13_v;
    double flux7_e, flux8_e, flux13_e;
    
    
    flux7_v = h_hit7->IntegralAndError(1,h_hit7->GetNbinsX(),
                                       1,h_hit7->GetNbinsY(),
                                       flux7_e);
    flux8_v = h_hit8->IntegralAndError(1,h_hit8->GetNbinsX(),
                                       1,h_hit8->GetNbinsY(),
                                       flux8_e);
    flux13_v = h_hit13->IntegralAndError(1,h_hit13->GetNbinsX(),
                                         1,h_hit13->GetNbinsY(),
                                         flux13_e);
    
    if (TMath::IsNaN(flux7_v)) {
      flux7_v = 0.0;
      flux7_e = 0.0;
    }
    
    if (TMath::IsNaN(flux8_v)) {
      flux8_v = 0.0;
      flux8_e = 0.0;
    }
    
    if (TMath::IsNaN(flux13_v)) {
      flux13_v = 0.0;
      flux13_e = 0.0;
    }
    
    Info("fluenceExtract","Sector %d:",sec);
    Printf("             Average flux: %10.7f+/-%10.7f 1-MeV n / cm^{2} collision for 7 TeV",
         sec,flux7_v,flux7_e);
    Printf("             Average flux: %10.7f+/-%10.7f 1-MeV n / cm^{2} collision for 8 TeV",
           sec,flux8_v,flux8_e);
    Printf("             Average flux: %10.7f+/-%10.7f 1-MeV n / cm^{2} collision for 13 TeV",
           sec,flux13_v,flux13_e);
    
    v_flux7.push_back(flux7_v);
    v_flux7_e.push_back(flux7_e);
    v_flux8.push_back(flux8_v);
    v_flux8_e.push_back(flux8_e);
    v_flux13.push_back(flux13_v);
    v_flux13_e.push_back(flux13_e);
    delete h_hit7;
    delete h_hit8;
    delete h_hit13;
  }
  
  TFile* f_out = new TFile(fn_out.Data(),"UPDATE");
  
  TVectorD vr_sector  = getROOTVector<double>(v_sector);
  TVectorD vr_flux7_v = getROOTVector<double>(v_flux7);
  TVectorD vr_flux7_e = getROOTVector<double>(v_flux7_e);
  TVectorD vr_flux8_v = getROOTVector<double>(v_flux8);
  TVectorD vr_flux8_e = getROOTVector<double>(v_flux8_e);
  TVectorD vr_flux13_v = getROOTVector<double>(v_flux13);
  TVectorD vr_flux13_e = getROOTVector<double>(v_flux13_e);
  
  if (useRadius) {
    lay+=Form("_r%d",(int)radius);
  }
  
  f_out->Delete(Form("v_sector_%s;*", lay.Data()));
  f_out->Delete(Form("v_flux7_v_%s;*",lay.Data()));
  f_out->Delete(Form("v_flux7_e_%s;*",lay.Data()));
  f_out->Delete(Form("v_flux8_v_%s;*",lay.Data()));
  f_out->Delete(Form("v_flux8_e_%s;*",lay.Data()));
  f_out->Delete(Form("v_flux13_v_%s;*",lay.Data()));
  f_out->Delete(Form("v_flux13_e_%s;*",lay.Data()));
  
  vr_sector.Write(Form("v_sector_%s", lay.Data()));
  vr_flux7_v.Write(Form("v_flux7_v_%s",lay.Data()));
  vr_flux7_e.Write(Form("v_flux7_e_%s",lay.Data()));
  vr_flux8_v.Write(Form("v_flux8_v_%s",lay.Data()));
  vr_flux8_e.Write(Form("v_flux8_e_%s",lay.Data()));
  vr_flux13_v.Write(Form("v_flux13_v_%s",lay.Data()));
  vr_flux13_e.Write(Form("v_flux13_e_%s",lay.Data()));
  
  f_out->Close();
  
  Info("fluenceExtract","Results written to %s",fn_out.Data());
  
  return EXIT_SUCCESS;
}