//
// 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 "../../include/incBasicROOT.h"
#include "../../include/incDrawROOT.h"
#include "../../include/incGraphROOT.h"
#include "../../include/incHistROOT.h"
#include "../../include/incIOROOT.h"
#include "../../include/incRooFit.h"
#include "../../include/basicROOTTools.h"
#include "../../include/basicRooFitTools.h"
#include "../../include/lhcbstyle.C"
#include "deplTool.h"
#include "hamburgTool.h"
// -d detector (TT/IT)
// -s set 0: all;
// 1: TT [125-160] V / IT [ 75-150] V;
// 2: TT [160-220] V / IT [150-220] V;
// 3: 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 ("CCEHOME");
TString det("TT");
TString lay("TTaU");
int set = 0;
int val = 3;
bool notSave = false;
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]);
}else if (strcmp(argv[i],"-v")==0 && i+1<argc) {
val = std::atoi(argv[++i]);
}else if (strcmp(argv[i],"-ns")==0) {
notSave = true;
}
}
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/data/db",homeVar));
TString fn_deplV(Form("%s/vdepl.root",dn_deplV.Data()));
TString dn_fill(Form("%s/data/db",homeVar));
TString fn_fill(Form("%s/lumi.root",dn_fill.Data()));
TString dn_sector(Form("%s/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());
// Define ranges of radii
const Int_t nRTT = 4;
Double_t rTT[nRTT] = {-1.0,75.0,45.0,40.0};
// Define ranges of radii
const Int_t nRIT = 2;
Double_t rIT[nRIT] = {-1.0,175.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");
TVectorD* vp_lumi13 = (TVectorD*)f_fill->Get("v_lumi13");
if (!vp_fill || !vp_lumi7 || !vp_lumi8 || !vp_lumi13) {
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;
TVectorD v_lumi13 = *vp_lumi13;
// 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;
/* --- Comment this bit out (Elena)
// No radius is taken into account in case of IT
--- */
if (det.EqualTo("IT")) {
// effnR = 1;
effnR = nRIT;
}
// We plot only one prediction... not one per sector!
// 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("****");
// Normalize plot to the starting depletion voltage
v_vdepl_v_theo.push_back(TMath::Abs(startV-v_func[0]) / startV);
v_vdepl_e_theo.push_back(v_func[1] / startV);
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);
// Make graphs out of them
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;
// Boundaries for sectors sets
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,
for (int i=0; i<v_fill.GetNoElements(); i++) {
// Loop over all sectors
for (int j=0; j<v_sector.GetNoElements(); j++) {
// Loop over all radii
for (int k=0; k<effnR; k++) {
// Obtain corresponding fluence 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_flux13_val(Form("v_flux13_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_flux13_err(Form("v_flux13_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 fluence vectors 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_flux13_val = (TVectorD*)f_sector->Get(vn_flux13_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());
TVectorD* vp_flux13_err = (TVectorD*)f_sector->Get(vn_flux13_err.Data());
if (!vp_sector_flux ||
!vp_flux7_val || !vp_flux7_err ||
!vp_flux13_val || !vp_flux13_err ||
!vp_flux8_val || !vp_flux8_err) {
Error("deplVphiN","Missing flux information for radius = %6.2f mm",rTT[k]>0.0?rTT[k]:2000.0);
continue;
}
// Read data from CCEScan.root
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_flux13_val = *vp_flux13_val;
TVectorD v_flux13_err = *vp_flux13_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_flux7_err(i_flux);
double flux8_val = v_flux8_val(i_flux);
double flux8_err = v_flux8_err(i_flux);
double flux13_val = v_flux13_val(i_flux);
double flux13_err = v_flux13_err(i_flux);
double lumi7 = v_lumi7(i);
double lumi8 = v_lumi8(i);
double lumi13 = v_lumi13(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 + flux13_val*lumi13);
phiN_err *= (flux7_err*lumi7 + flux8_err*lumi8 + flux13_err*lumi13);
// 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;
}
// Remove excluded sectors
if (STTool::IsExcluded(sector)) {
Warning("deplVphiN","Sector excluded");
continue;
}
// Add systematic uncertainty to get mean value of all uncertainties. Normalize
mean_sys += vdepl_sys/vdepl_prod;
nMeas++;
// Print data
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 (NOT NEEDED)
//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 data to the proper vectors and normalize
if (det.EqualTo("TT")) {
if (vdepl_prod>v2Bound) {
if (rTT[k]>0.0) {
v_vdepl_v_r_high.push_back(vdepl_val/vdepl_prod);
v_vdepl_e_r_high.push_back(vdepl_tot/vdepl_prod);
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/vdepl_prod);
v_vdepl_e_norm_high.push_back(vdepl_tot/vdepl_prod);
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/vdepl_prod);
v_vdepl_e_r_mid.push_back(vdepl_tot/vdepl_prod);
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/vdepl_prod);
v_vdepl_e_norm_mid.push_back(vdepl_tot/vdepl_prod);
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/vdepl_prod);
v_vdepl_e_r_low.push_back(vdepl_tot/vdepl_prod);
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/vdepl_prod);
v_vdepl_e_norm_low.push_back(vdepl_tot/vdepl_prod);
v_phiN_v_norm_low.push_back(phiN_val);
v_phiN_e_norm_low.push_back(phiN_err);
}
}
}else{
if (vdepl_prod>v2Bound) { // Sectors with higher Vdepl after production
if (rIT[k]>0.0) { // Inner circles
v_vdepl_v_r_high.push_back(vdepl_val/vdepl_prod);
v_vdepl_e_r_high.push_back(vdepl_tot/vdepl_prod);
v_phiN_v_r_high.push_back(phiN_val);
v_phiN_e_r_high.push_back(phiN_err);
}else{ // Whole sectors
v_vdepl_v_norm_high.push_back(vdepl_val/vdepl_prod);
v_vdepl_e_norm_high.push_back(vdepl_tot/vdepl_prod);
v_phiN_v_norm_high.push_back(phiN_val);
v_phiN_e_norm_high.push_back(phiN_err);
}
}else if(vdepl_prod>v1Bound) { // Sectors with lower Vdepl after production
if (rIT[k]>0.0) { // Inner circles
v_vdepl_v_r_mid.push_back(vdepl_val/vdepl_prod);
v_vdepl_e_r_mid.push_back(vdepl_tot/vdepl_prod);
v_phiN_v_r_mid.push_back(phiN_val);
v_phiN_e_r_mid.push_back(phiN_err);
}else{ // Whole sectors
v_vdepl_v_norm_mid.push_back(vdepl_val/vdepl_prod);
v_vdepl_e_norm_mid.push_back(vdepl_tot/vdepl_prod);
v_phiN_v_norm_mid.push_back(phiN_val);
v_phiN_e_norm_mid.push_back(phiN_err);
}
}else{ // Sectors with lowest Vdepl after production
if (rIT[k]>0.0) { // Inner circles
v_vdepl_v_r_low.push_back(vdepl_val/vdepl_prod);
v_vdepl_e_r_low.push_back(vdepl_tot/vdepl_prod);
v_phiN_v_r_low.push_back(phiN_val);
v_phiN_e_r_low.push_back(phiN_err);
}else{ // Whole sectors
v_vdepl_v_norm_low.push_back(vdepl_val/vdepl_prod);
v_vdepl_e_norm_low.push_back(vdepl_tot/vdepl_prod);
v_phiN_v_norm_low.push_back(phiN_val);
v_phiN_e_norm_low.push_back(phiN_err);
}
}
} // end if (IT or TT)
} // end loop on radii
} // end loop on sectors
} // end loop on fills
// Calculate mean systematic uncertainty
if (nMeas>0) {
mean_sys /= nMeas;
}
// Wait, what constant is this one?!
// Add systmatic uncertainty from the extraction procedure of the depletion voltage. Normalize
if (det.EqualTo("IT")) {
mean_sys = TMath::Sqrt(mean_sys*mean_sys + (STTool::voltSyst_IT*STTool::voltSyst_IT / (startV*startV)));
}
if (det.EqualTo("TT")) {
mean_sys = TMath::Sqrt(mean_sys*mean_sys + (STTool::voltSyst_TT*STTool::voltSyst_TT / (startV*startV)));
}
// We are summing systematics from the data to systematics from the predictions!!!!
// Drawing prediction limits
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_theo,"E3");
mg_frame->Add(ge_theo_noerr,"L");
mg_frame->Add(ge_theo_high,"L");
mg_frame->Add(ge_theo_low,"L");
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");
// Plot graphs
mg_frame->Draw("a");
mg_frame->GetXaxis()->SetTitle("#phi_{1 MeV-n,eq} [cm^{-2}]");
// Normalize
if (det.EqualTo("TT")) {
mg_frame->GetYaxis()->SetTitle("TT sectors #it{V}_{depl} / #it{V}_{0}");
} else {
mg_frame->GetYaxis()->SetTitle("IT sectors #it{V}_{depl} / #it{V}_{0}");
}
mg_frame->GetXaxis()->SetLimits(1e+8,1e+14);
mg_frame->GetXaxis()->SetLabelOffset(-0.005);
mg_frame->GetXaxis()->SetTitleOffset(1.2);
// Normalize
mg_frame->SetMaximum(300.0 / startV);
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");
// Directory and file name for output graph
TString dn_graph(Form("%s/data/ST/graphics/%s/%s",
diskVar, det.Data(), lay.Data()));
TString fn_graph(Form("%s/deplVphiN_normalized_v%d",
dn_graph.Data(), val));
fn_graph += ".pdf";
if (!notSave) {
// Save picture
gSystem->Exec(Form("mkdir -p %s", dn_graph.Data()));
gSystem->Exec(Form("rm %s &> /dev/null", fn_graph.Data()));
c_frame->SaveAs(fn_graph.Data());
}
return EXIT_SUCCESS;
}
Elena Graverini