//  
//  deplVphiN.C
//  Script to plot the result of the depletion voltage measurement as a function of the fluence
//
//  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 <vector>
#include <map>
#include <stdlib.h>
#include <limits>
#include <string>
#include <time.h>


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


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


#include "lhcbstyle.C"


#include "deplTool.h"
#include "hamburgTool.h"

double vTop     = 280.0;
double vBottom  = 125.0;
double v1Bound  = 160.0;
double v2Bound  = 220.0;

if (det.EqualTo("IT")) {
  vTop    = 250.0;
  vBottom =  75.0;
  v1Bound = 150.0;
  v2Bound = 220.0;
}


// -d detector (TT/IT)
// -s set 0: all;
//        1: TT [125-160] V / IT [ 75-150] V;
//        1: TT [160-220] V / IT [150-220] V;
//        1: TT [220-280] V / IT [220-250] V;
//        Voltage values are the initial depletion voltage values measured after production


int deplVphiN(int argc, char* argv[]);


int main(int argc, char* argv[]){
  
  TStyle* style = lhcbStyle();
  
  int argc_copy = argc;
  
  TApplication* theApp = new TApplication("Analysis", &argc, argv);
  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();
    }
  }
  

  
  int exit = deplVphiN(argc,argv);
  
  
  
  Printf("End");
  if (!runBatch) {
    theApp->Run(!runBatch);
  }
  
  
  return exit;
}


//  Executable method
int deplVphiN(int argc, char* argv[]){
  
  TVectorD v_res(5);
  
  gStyle->SetPadLeftMargin(0.20);
  gStyle->SetTitleOffset(1.5,"Y");
  gStyle->SetTitleOffset(1.0,"X");

  const char* diskVar = std::getenv ("DISK");
  const char* homeVar = std::getenv ("HOME");
  

  TString det("TT");
  TString lay("TTaU");
  int set       = 0;
  
  
  for (int i=1; i<argc; i++) {
    if (strcmp(argv[i],"-d")==0 && i+1<argc) {
      det = argv[++i];
    }else if (strcmp(argv[i],"-s")==0 && i+1<argc) {
      set = std::atoi(argv[++i]);
    }
  }
  
  
  if (det.EqualTo("IT")) {
    lay = "T3X2";
  }else{
    lay = "TTaU";
  }
  
  // Directory and file names for the measured Vdepl values, the fluence values,
  // the Vdepl values after production and the luminosity (i.e. development of the flux over time)
  
  TString dn_data(Form("%s/data/ST/Aging",diskVar));
  TString fn_data(Form("%s/CCEScan.root",dn_data.Data()));
  
  TString dn_deplV(Form("%s/cmtuser/Vetra_v13r2/ST/STAging/data/db",homeVar));
  TString fn_deplV(Form("%s/vdepl.root",dn_deplV.Data()));
  
  TString dn_fill(Form("%s/cmtuser/Vetra_v13r2/ST/STAging/data/db",homeVar));
  TString fn_fill(Form("%s/lumi.root",dn_fill.Data()));
  
  TString dn_sector(Form("%s/cmtuser/Vetra_v13r2/ST/STAging/data/db",homeVar));
  TString fn_sector(Form("%s/fluence.root",dn_sector.Data()));
  
  TFile* f_data   = TFile::Open(fn_data.Data());
  TFile* f_deplV  = TFile::Open(fn_deplV.Data());
  TFile* f_fill   = TFile::Open(fn_fill.Data());
  TFile* f_sector = TFile::Open(fn_sector.Data());
  
  const Int_t nRTT = 4;
  
  
  // Define ranges of radii
  Double_t rTT[nRTT] = {-1.0,75.0,45.0,40.0};
  

  // Get information about the luminosity and the
  TVectorD* vp_fill  = (TVectorD*)f_fill->Get("v_fill");
  TVectorD* vp_lumi7 = (TVectorD*)f_fill->Get("v_lumi7");
  TVectorD* vp_lumi8 = (TVectorD*)f_fill->Get("v_lumi8");
  
  if (!vp_fill || !vp_lumi7 || !vp_lumi8) {
    Error("deplVphiN","Luminosity vectors are not available");
    return EXIT_FAILURE;
  }
  
  TVectorD v_fill    = *vp_fill;
  TVectorD v_lumi7   = *vp_lumi7;
  TVectorD v_lumi8   = *vp_lumi8;
  
  
  // Get numbers of available read-out sectors and corresponding initial depletion voltage values
  TVectorD* vp_sector  = (TVectorD*)f_deplV->Get(Form("v_sector_%s",lay.Data()));
  TVectorD* vp_vdepl   = (TVectorD*)f_deplV->Get(Form("v_vdepl_%s",lay.Data()));
  
  if (!vp_sector || !vp_vdepl) {
    Error("deplVphiN","Depletion Voltage vectors are not available");
    return EXIT_FAILURE;
  }
  
  TVectorD v_sector  = *vp_sector;
  TVectorD v_vdepl   = *vp_vdepl;
  
  Int_t effnR = nRTT;
  
  // No radius is taken into account in case of IT
  if (det.EqualTo("IT")) {
    effnR = 1;
  }
  
  // Define the starting depletion voltage for the predictions

  double startV   = 200.0;
  if (set==1) {
    startV = 160.0;
  }else if(set==3){
    startV = 240.0;
  }
  if (det.EqualTo("IT")) {
    startV = 100.0;
    if (set==1) {
      startV = 80.0;
    }else if(set==3){
      startV = 120.0;
    }
  }
  
  
  // Prediction calculation based on the stable damage part in the Hamburg model
  int nTheoPoints = 1000;
  double minTheo  = 1e+6;
  double maxTheo  = 1e+20;
  double binTheo  = TMath::Power(maxTheo/minTheo,1.0/nTheoPoints);
  
  std::vector<double> v_vdepl_v_theo;
  std::vector<double> v_vdepl_e_theo;
  std::vector<double> v_phiN_v_theo;
  std::vector<double> v_phiN_e_theo;
  
  Info("deplVphiN","Calculating theoretical prediction");
  for (int i=0; i<nTheoPoints; i++) {
    double phiN_point = minTheo*TMath::Power(binTheo,i);
    std::vector<double> v_func = STHamburg::getVeff(phiN_point,det.EqualTo("IT"));
    Printf("   phi = %20.2f",phiN_point);
    Printf("   V   = (%6.2f+/-%6.2f) V",TMath::Abs(startV-v_func[0]),v_func[1]);
    Printf("****");
    v_vdepl_v_theo.push_back(TMath::Abs(startV-v_func[0]));
    v_vdepl_e_theo.push_back(v_func[1]);
    v_phiN_v_theo.push_back(phiN_point);
    v_phiN_e_theo.push_back(0.0);
  }
  
  // Convert vectors related to predictions into root vectors
  
  TVectorD vr_vdepl_v_theo = getROOTVector<double>(v_vdepl_v_theo);
  TVectorD vr_vdepl_e_theo = getROOTVector<double>(v_vdepl_e_theo);
  TVectorD vr_phiN_v_theo  = getROOTVector<double>(v_phiN_v_theo);
  TVectorD vr_phiN_e_theo  = getROOTVector<double>(v_phiN_e_theo);
  
  
  
  TGraphErrors* ge_theo = new TGraphErrors(vr_phiN_v_theo,
                                           vr_vdepl_v_theo,
                                           vr_phiN_e_theo,
                                           vr_vdepl_e_theo);
  TGraph* ge_theo_noerr = new TGraph(vr_phiN_v_theo,
                                     vr_vdepl_v_theo);
  ge_theo->SetFillColor(kGray+2);
  
  TMultiGraph* mg_frame = new TMultiGraph("mg_frame","Frame multigraph");
  
  
  
  TCanvas* c_frame = new TCanvas("c_frame","Frame Canvas",800,600);
  c_frame->SetLogx();
  mg_frame->Draw("");
  
  // Vectors to store measurement points
  
  std::vector<double> v_vdepl_v_norm_low;
  std::vector<double> v_vdepl_v_norm_mid;
  std::vector<double> v_vdepl_v_norm_high;
  std::vector<double> v_vdepl_v_r_low;
  std::vector<double> v_vdepl_v_r_mid;
  std::vector<double> v_vdepl_v_r_high;
  std::vector<double> v_vdepl_e_norm_low;
  std::vector<double> v_vdepl_e_norm_mid;
  std::vector<double> v_vdepl_e_norm_high;
  std::vector<double> v_vdepl_e_r_low;
  std::vector<double> v_vdepl_e_r_mid;
  std::vector<double> v_vdepl_e_r_high;
  std::vector<double> v_phiN_v_norm_low;
  std::vector<double> v_phiN_v_norm_mid;
  std::vector<double> v_phiN_v_norm_high;
  std::vector<double> v_phiN_v_r_low;
  std::vector<double> v_phiN_v_r_mid;
  std::vector<double> v_phiN_v_r_high;
  std::vector<double> v_phiN_e_norm_low;
  std::vector<double> v_phiN_e_norm_mid;
  std::vector<double> v_phiN_e_norm_high;
  std::vector<double> v_phiN_e_r_low;
  std::vector<double> v_phiN_e_r_mid;
  std::vector<double> v_phiN_e_r_high;
  
  
  double vTop     = 280.0;
  double vBottom  = 125.0;
  double v1Bound  = 160.0;
  double v2Bound  = 220.0;
  
  if (det.EqualTo("IT")) {
    vTop    = 250.0;
    vBottom =  75.0;
    v1Bound = 150.0;
    v2Bound = 220.0;
  }
  
  double vMax     = vTop;
  double vMin     = vBottom;
  
  if (set==1) {
    vMax = v1Bound;
  }else if(set==2){
    vMax = v2Bound;
    vMin = v1Bound;
  }else if(set==3){
    vMin = v2Bound;
  }
  
  double mean_sys = 0.0;
  int nMeas       = 0;
  
  // Loop over all fills, sectors and radii
  
  for (int i=0; i<v_fill.GetNoElements(); i++) {
    for (int j=0; j<v_sector.GetNoElements(); j++) {
      for (int k=0; k<effnR; k++) {
        
        // Obtain corresponding vector names
        
        int sector = (int) v_sector(j);
        int fill   = (int) v_fill(i);
        TString s_rad("");
        TString vn_sector_flux(Form("v_sector_%s",lay.Data()));
        TString vn_flux7_val(Form("v_flux7_v_%s",lay.Data()));
        TString vn_flux8_val(Form("v_flux8_v_%s",lay.Data()));
        TString vn_flux7_err(Form("v_flux7_e_%s",lay.Data()));
        TString vn_flux8_err(Form("v_flux8_e_%s",lay.Data()));
        
        TString vn_data(Form("%s/%s/%d/%d/v_volt_val%d",
                             det.Data(),lay.Data(),sector,fill,val));
        
        // Add radius postfix
        if (rTT[k]>0.0) {
          s_rad = Form("%d",(int)rTT[k]);
          vn_sector_flux  += "_r"+s_rad;
          vn_flux7_val    += "_r"+s_rad;
          vn_flux8_val    += "_r"+s_rad;
          vn_flux7_err    += "_r"+s_rad;
          vn_flux8_err    += "_r"+s_rad;
          vn_data         += "_r"+s_rad;
        }
        
        
        // Extract them from the file
        TVectorD* vp_sector_flux = (TVectorD*)f_sector->Get(vn_sector_flux.Data());
        TVectorD* vp_flux7_val   = (TVectorD*)f_sector->Get(vn_flux7_val.Data());
        TVectorD* vp_flux8_val   = (TVectorD*)f_sector->Get(vn_flux8_val.Data());
        TVectorD* vp_flux7_err   = (TVectorD*)f_sector->Get(vn_flux7_err.Data());
        TVectorD* vp_flux8_err   = (TVectorD*)f_sector->Get(vn_flux8_err.Data());
        
        if (!vp_sector_flux ||
            !vp_flux7_val || !vp_flux7_err ||
            !vp_flux8_val || !vp_flux8_err) {
          Error("deplVphiN","Flux information for radius = %6.2f mm",rTT[k]>0.0?rTT[k]:2000.0);
          continue;
        }
        
        TVectorD* vp_data = (TVectorD*)f_data->Get(vn_data.Data());
        
        if (!vp_data){
          Warning("deplVphiN","No vector found for %s",vn_data.Data());
          continue;
        }
        
        TVectorD v_sector_flux = *vp_sector_flux;
        TVectorD v_flux7_val   = *vp_flux7_val;
        TVectorD v_flux7_err   = *vp_flux7_err;
        TVectorD v_flux8_val   = *vp_flux8_val;
        TVectorD v_flux8_err   = *vp_flux8_err;
        
        
        TVectorD v_data        = *vp_data;
        
        int i_flux = getClosestIndex(v_sector_flux,(double)sector);
        
        double flux7_val = v_flux7_val(i_flux);
        double flux7_err = v_flux8_err(i_flux);
        double flux8_val = v_flux7_val(i_flux);
        double flux8_err = v_flux8_err(i_flux);
        
        double lumi7     = v_lumi7(i);
        double lumi8     = v_lumi8(i);
        
        double vdepl_prod = v_vdepl(j);
        
        double vdepl_val = v_data(0);
        double vdepl_err = v_data(1);
        double vdepl_sys = v_data(2);
        
        double vdepl_tot = vdepl_err;
        
        double phiN_val = STTool::FLUKAConvFac;
        double phiN_err = STTool::FLUKAConvFac;
        
        phiN_val *= (flux7_val*lumi7+flux8_val*lumi8);
        phiN_err *= (flux7_err*lumi7+flux8_err*lumi8);
        
        
        // Define maximal and minimal Vdepl values to be drawn
        double plotMin = 40.0;
        double plotMax = 290.0;
        
        if (det.EqualTo("TT")) {
          plotMin = 80.0;
        }
        
        
        // Remove points with too high uncertainty
        if (vdepl_err>20.0 ){
          Warning("deplVphiN","Too high or low uncertainty on V_depl");
          continue;
        }
      
        
        
        if (vdepl_val>plotMax || vdepl_val<plotMin) {
          Warning("deplVphiN","Too high or low values for V_depl");
          continue;
        }
        
       
        if (vdepl_prod>vMax || vdepl_prod<=vMin) {
          continue;
        }
        
        // Add systematic uncertainty to get mean value
        mean_sys += vdepl_sys;
        nMeas++;
        
        // Remove excluded sectors
        if (STTool::IsExcluded(sector)) {
          Warning("deplVphiN","Sector excluded");
          continue;
        }
        
        Info("deplVphiN","Result for %s/%s sector %d, fill %d, radius %7.2f mm:",
             det.Data(),lay.Data(),sector,fill,rTT[k]>0.0?rTT[k]:2000.0);
        Printf("     phi_N        = (%10.2f+/-%7.2f)*10^12 cm^{-2}",
               phiN_val/1e12,phiN_err/1e12);
        Printf("     V_depl       = (%6.2f+/-%6.2f(stat)+/-%6.2f(syst)) V",
               vdepl_val,vdepl_err,vdepl_sys);
        Printf("     V_depl,prod  = %6.2f V",
               vdepl_prod);
        
        // Define data point in plot
        TMarker* m_res = new TMarker(phiN_val,vdepl_val,20);
        TLine* l_hor   = new TLine(phiN_val-phiN_err,vdepl_val,
                                   phiN_val+phiN_err,vdepl_val);
        TLine* l_ver   = new TLine(phiN_val,vdepl_val-vdepl_err,
                                   phiN_val,vdepl_val+vdepl_err);
        
        // Add them to the proper vector
        
        if (det.EqualTo("TT")) {
          if (vdepl_prod>v2Bound) {
            if (rTT[k]>0.0) {
              v_vdepl_v_r_high.push_back(vdepl_val);
              v_vdepl_e_r_high.push_back(vdepl_tot);
              v_phiN_v_r_high.push_back(phiN_val);
              v_phiN_e_r_high.push_back(phiN_err);
            }else{
              v_vdepl_v_norm_high.push_back(vdepl_val);
              v_vdepl_e_norm_high.push_back(vdepl_tot);
              v_phiN_v_norm_high.push_back(phiN_val);
              v_phiN_e_norm_high.push_back(phiN_err);
            }
          }else if(vdepl_prod>v1Bound) {
            if (rTT[k]>0.0) {
              v_vdepl_v_r_mid.push_back(vdepl_val);
              v_vdepl_e_r_mid.push_back(vdepl_tot);
              v_phiN_v_r_mid.push_back(phiN_val);
              v_phiN_e_r_mid.push_back(phiN_err);
            }else{
              v_vdepl_v_norm_mid.push_back(vdepl_val);
              v_vdepl_e_norm_mid.push_back(vdepl_tot);
              v_phiN_v_norm_mid.push_back(phiN_val);
              v_phiN_e_norm_mid.push_back(phiN_err);
            }
          }else{
            if (rTT[k]>0.0) {
              v_vdepl_v_r_low.push_back(vdepl_val);
              v_vdepl_e_r_low.push_back(vdepl_tot);
              v_phiN_v_r_low.push_back(phiN_val);
              v_phiN_e_r_low.push_back(phiN_err);
            }else{
              v_vdepl_v_norm_low.push_back(vdepl_val);
              v_vdepl_e_norm_low.push_back(vdepl_tot);
              v_phiN_v_norm_low.push_back(phiN_val);
              v_phiN_e_norm_low.push_back(phiN_err);
            }
          }
        }else{
          if (vdepl_prod>v2Bound) {
            v_vdepl_v_norm_high.push_back(vdepl_val);
            v_vdepl_e_norm_high.push_back(vdepl_tot);
            v_phiN_v_norm_high.push_back(phiN_val);
            v_phiN_e_norm_high.push_back(phiN_err);
          }else if(vdepl_prod>v1Bound) {
            v_vdepl_v_norm_mid.push_back(vdepl_val);
            v_vdepl_e_norm_mid.push_back(vdepl_tot);
            v_phiN_v_norm_mid.push_back(phiN_val);
            v_phiN_e_norm_mid.push_back(phiN_err);
          }else{
            v_vdepl_v_norm_low.push_back(vdepl_val);
            v_vdepl_e_norm_low.push_back(vdepl_tot);
            v_phiN_v_norm_low.push_back(phiN_val);
            v_phiN_e_norm_low.push_back(phiN_err);
          }
        }
      }
    }
  }
  
  // Calculate mean systematic uncertainty
  if (nMeas>0) {
    mean_sys /= nMeas;
  }
  
  // Add systmatic uncertainty from the extraction procedure of the depletion voltage
  if (det.EqualTo("IT")) {
    mean_sys = TMath::Sqrt(mean_sys*mean_sys+STTool::voltSyst_IT*STTool::voltSyst_IT);
  }
  if (det.EqualTo("TT")) {
    mean_sys = TMath::Sqrt(mean_sys*mean_sys+STTool::voltSyst_TT*STTool::voltSyst_TT);
  }
  
  
  vr_vdepl_v_theo+=mean_sys;
  
  TGraph* ge_theo_high   = new TGraph(vr_phiN_v_theo,
                                      vr_vdepl_v_theo);
  ge_theo_high->SetLineStyle(9);
  
  vr_vdepl_v_theo-=2*mean_sys;
  
  TGraph* ge_theo_low    = new TGraph(vr_phiN_v_theo,
                                      vr_vdepl_v_theo);
  ge_theo_low->SetLineStyle(9);
  
  
  // Transfrom STL vectors to ROOT vectors
  TVectorD vr_vdepl_v_norm_low   = getROOTVector<double>(v_vdepl_v_norm_low);
  TVectorD vr_vdepl_v_norm_mid   = getROOTVector<double>(v_vdepl_v_norm_mid);
  TVectorD vr_vdepl_v_norm_high  = getROOTVector<double>(v_vdepl_v_norm_high);
  TVectorD vr_vdepl_v_r_low      = getROOTVector<double>(v_vdepl_v_r_low);
  TVectorD vr_vdepl_v_r_mid      = getROOTVector<double>(v_vdepl_v_r_mid);
  TVectorD vr_vdepl_v_r_high     = getROOTVector<double>(v_vdepl_v_r_high);
  TVectorD vr_vdepl_e_norm_low   = getROOTVector<double>(v_vdepl_e_norm_low);
  TVectorD vr_vdepl_e_norm_mid   = getROOTVector<double>(v_vdepl_e_norm_mid);
  TVectorD vr_vdepl_e_norm_high  = getROOTVector<double>(v_vdepl_e_norm_high);
  TVectorD vr_vdepl_e_r_low      = getROOTVector<double>(v_vdepl_e_r_low);
  TVectorD vr_vdepl_e_r_mid      = getROOTVector<double>(v_vdepl_e_r_mid);
  TVectorD vr_vdepl_e_r_high     = getROOTVector<double>(v_vdepl_e_r_high);
  TVectorD vr_phiN_v_norm_low    = getROOTVector<double>(v_phiN_v_norm_low);
  TVectorD vr_phiN_v_norm_mid    = getROOTVector<double>(v_phiN_v_norm_mid);
  TVectorD vr_phiN_v_norm_high   = getROOTVector<double>(v_phiN_v_norm_high);
  TVectorD vr_phiN_v_r_low       = getROOTVector<double>(v_phiN_v_r_low);
  TVectorD vr_phiN_v_r_mid       = getROOTVector<double>(v_phiN_v_r_mid);
  TVectorD vr_phiN_v_r_high      = getROOTVector<double>(v_phiN_v_r_high);
  TVectorD vr_phiN_e_norm_low    = getROOTVector<double>(v_phiN_e_norm_low);
  TVectorD vr_phiN_e_norm_mid    = getROOTVector<double>(v_phiN_e_norm_mid);
  TVectorD vr_phiN_e_norm_high   = getROOTVector<double>(v_phiN_e_norm_high);
  TVectorD vr_phiN_e_r_low       = getROOTVector<double>(v_phiN_e_r_low);
  TVectorD vr_phiN_e_r_mid       = getROOTVector<double>(v_phiN_e_r_mid);
  TVectorD vr_phiN_e_r_high      = getROOTVector<double>(v_phiN_e_r_high);
  
  
  
  // Generate graphs
  double markerSize = 1.0;
  
  TGraphErrors* ge_norm_low      = new TGraphErrors(vr_phiN_v_norm_low,
                                                    vr_vdepl_v_norm_low,
                                                    vr_phiN_e_norm_low,
                                                    vr_vdepl_e_norm_low);
  ge_norm_low->SetMarkerSize(markerSize);
  ge_norm_low->SetMarkerStyle(23);
  ge_norm_low->SetMarkerColor(kGreen+3);
  ge_norm_low->SetLineColor(kGreen+3);
  
  TGraphErrors* ge_norm_mid      = new TGraphErrors(vr_phiN_v_norm_mid,
                                                    vr_vdepl_v_norm_mid,
                                                    vr_phiN_e_norm_mid,
                                                    vr_vdepl_e_norm_mid);
  ge_norm_mid->SetMarkerSize(markerSize);
  ge_norm_mid->SetMarkerStyle(22);
  ge_norm_mid->SetMarkerColor(kBlue);
  ge_norm_mid->SetLineColor(kBlue);
  
  TGraphErrors* ge_norm_high     = new TGraphErrors(vr_phiN_v_norm_high,
                                                    vr_vdepl_v_norm_high,
                                                    vr_phiN_e_norm_high,
                                                    vr_vdepl_e_norm_high);
  ge_norm_high->SetMarkerSize(markerSize);
  ge_norm_high->SetMarkerStyle(20);
  ge_norm_high->SetMarkerColor(kRed);
  ge_norm_high->SetLineColor(kRed);
  
  
  TGraphErrors* ge_r_low      = new TGraphErrors(vr_phiN_v_r_low,
                                                 vr_vdepl_v_r_low,
                                                 vr_phiN_e_r_low,
                                                 vr_vdepl_e_r_low);
  ge_r_low->SetMarkerSize(markerSize);
  ge_r_low->SetMarkerStyle(32);
  ge_r_low->SetMarkerColor(kGreen+3);
  ge_r_low->SetLineColor(kGreen+3);
  
  TGraphErrors* ge_r_mid      = new TGraphErrors(vr_phiN_v_r_mid,
                                                 vr_vdepl_v_r_mid,
                                                 vr_phiN_e_r_mid,
                                                 vr_vdepl_e_r_mid);
  ge_r_mid->SetMarkerSize(markerSize);
  ge_r_mid->SetMarkerStyle(26);
  ge_r_mid->SetMarkerColor(kBlue);
  ge_r_mid->SetLineColor(kBlue);
  
  TGraphErrors* ge_r_high     = new TGraphErrors(vr_phiN_v_r_high,
                                                 vr_vdepl_v_r_high,
                                                 vr_phiN_e_r_high,
                                                 vr_vdepl_e_r_high);
  ge_r_high->SetMarkerSize(markerSize);
  ge_r_high->SetMarkerStyle(24);
  ge_r_high->SetMarkerColor(kRed);
  ge_r_high->SetLineColor(kRed);
  
  mg_frame->Add(ge_norm_low,"pe");
  mg_frame->Add(ge_norm_mid,"pe");
  mg_frame->Add(ge_norm_high,"pe");
  mg_frame->Add(ge_r_low,"p");
  mg_frame->Add(ge_r_mid,"p");
  mg_frame->Add(ge_r_high,"p");
  mg_frame->Draw("a");
  
  // Plot graphs
  
  mg_frame->GetXaxis()->SetTitle("#phi_{1 MeV-n,eq} [cm^{-2}]");
  mg_frame->GetYaxis()->SetTitle("#it{V}_{depl} [V]");
  mg_frame->GetXaxis()->SetLimits(1e+8,1e+14);
  mg_frame->GetXaxis()->SetLabelOffset(-0.005);
  mg_frame->GetXaxis()->SetTitleOffset(1.2);
  mg_frame->SetMaximum(300.0);
  mg_frame->SetMinimum(0.0);
  
  
  
  // Draw legend
  
  double baseY = 0.30;
  if (det.EqualTo("IT")) {
    baseY = 0.90;
  }
  
  
  
  TMarker* m_leg = new TMarker(0.25,baseY,20);
  TLatex* t_leg  = new TLatex(0.27,baseY,"#it{V}_{depl} > 225 V");
  t_leg->SetTextSize(0.04);
  t_leg->SetTextAlign(12);
  t_leg->SetNDC();
  m_leg->SetNDC();
  
  if (det.EqualTo("TT") && (set==3 || set==0)) {
    m_leg->SetX(0.25);
    m_leg->SetY(baseY);
    m_leg->SetMarkerColor(kRed);
    m_leg->SetMarkerStyle(20);
    m_leg->DrawClone();
    m_leg->SetMarkerStyle(24);
    m_leg->SetX(0.55);
    if (det.EqualTo("TT")) {
      m_leg->DrawClone();
    }
    t_leg->DrawLatex(0.30,baseY,Form("#it{V}_{depl} #in [%3.0f,%3.0f] V",v2Bound,vTop));
    if (det.EqualTo("TT")) {
      t_leg->DrawLatex(0.60,baseY,"Innermost region");
    }
    baseY -= 0.05;
  }
  
  if (det.EqualTo("TT") && (set==2 || set==0)) {
    m_leg->SetX(0.25);
    m_leg->SetY(baseY);
    m_leg->SetMarkerColor(kBlue);
    m_leg->SetMarkerStyle(22);
    m_leg->DrawClone();
    m_leg->SetMarkerStyle(26);
    m_leg->SetX(0.55);
    if (det.EqualTo("TT")) {
      m_leg->DrawClone();
    }
    t_leg->DrawLatex(0.30,baseY,Form("#it{V}_{depl} #in [%3.0f,%3.0f] V",v1Bound,v2Bound));
    if (det.EqualTo("TT")) {
      t_leg->DrawLatex(0.60,baseY,"Innermost region");
    }
    baseY -= 0.05;
  }
  
  if (det.EqualTo("IT") && (set==1 || set==0)) {
    m_leg->SetX(0.25);
    m_leg->SetY(baseY);
    m_leg->SetMarkerColor(kGreen+3);
    m_leg->SetMarkerStyle(23);
    m_leg->DrawClone();
    m_leg->SetMarkerStyle(32);
    m_leg->SetX(0.55);
    if (det.EqualTo("TT")) {
      m_leg->DrawClone();
    }
    t_leg->DrawLatex(0.30,baseY,Form("#it{V}_{depl} #in [%3.0f,%3.0f] V",vBottom,v1Bound));
    if (det.EqualTo("TT")) {
      t_leg->DrawLatex(0.60,baseY,"Innermost region");
    }
    baseY -= 0.05;
  }
  
  
  
  TLine* l_line = new TLine(0.20,baseY,0.25,baseY);
  l_line->DrawLineNDC(0.23,baseY,0.27,baseY);
  t_leg->DrawLatex(0.30,baseY,Form("Hamburg Model (#it{V}_{depl} = %3.0f V)",startV));
  
  mg_frame->GetHistogram()->Draw("axissame");
  
  
  
  Printf("=========================================\n");

  return EXIT_SUCCESS;
}


