diff --git a/macros/CCEScan/plotSector_paper.C b/macros/CCEScan/plotSector_paper.C
new file mode 100755
index 0000000..ed3e0dd
--- /dev/null
+++ b/macros/CCEScan/plotSector_paper.C
@@ -0,0 +1,411 @@
+//  
+//  plotSector.C
+//  Plot measured depletion voltage and Hamburg model for a certain read-out sector
+//
+//  Created by Christian Elsasser on 31.05.13.
+//  University of Zurich, elsasser@cern.ch
+//  Updated by Elena Graverini on 10.06.15
+//  University of Zurich, elena.graverini@cern.ch
+//
+//  compile with:
+//  g++ -o plotSector plotSector.C `root-config --glibs --cflags`
+//
+//  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 read-out sector number
+// -r consider only hit in a distance closer than r mm from the beam axis
+// -v number of summed up strips
+// -ns Do not save the value in the root file
+
+
+int plotSector(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 = plotSector(argc,argv);
+  
+  
+  
+  Printf("End");
+  if (!runBatch) {
+    theApp->Run(!runBatch);
+  }
+  
+  
+  return exit;
+}
+
+
+//  Executable method
+int plotSector(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 = "/afs/cern.ch/user/e/egraveri/cmtuser/STMonitoring/STAging";//std::getenv ("DISK");
+  const char* diskVar = std::getenv ("DISK");
+  //const char* homeVar = "/afs/cern.ch/user/e/egraveri/cmtuser/STMonitoring/STAging";//std::getenv ("HOME");
+  const char* homeVar = std::getenv ("CCEHOME");
+  
+
+  TString det("TT");
+  TString lay("TTaU");
+  int val        = 3;
+  int sector     = 2634;
+  double radius  = -1000;
+  bool notSave   = false;
+  bool use_height = false;
+  bool average_ratio = 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],"-v")==0 && i+1<argc) {
+      val = std::atoi(argv[++i]);
+    }else if (strcmp(argv[i],"-s")==0 && i+1<argc) {
+      sector = std::atoi(argv[++i]);
+    }else if (strcmp(argv[i],"-r")==0 && i+1<argc) {
+      radius = std::atof(argv[++i]);
+
+    }else if (strcmp(argv[i],"-use_height")==0) { // use height of pulse shape
+      use_height = true;
+    }else if (strcmp(argv[i],"-average_ratio")==0) { // use the average ratio found for each layer and val
+      average_ratio = true;
+
+    }else if (strcmp(argv[i],"-ns")==0) {
+      notSave = true;
+    }
+  }
+  
+  
+  if (det.EqualTo("IT")) {
+    lay = "T3X2";
+  }else{
+    lay = "TTaU";
+  }
+  
+  // Directory and file name for output graph
+  TString dn_graph(Form("%s/data/ST/graphics/%s/%s/%d",
+                        diskVar,det.Data(),lay.Data(),sector));
+  TString fn_graph(Form("%s/time_evolv_v%d",
+                        dn_graph.Data(),val));
+  
+  if (radius>0.0) {
+    fn_graph += Form("_r%d",(int)radius);
+  }
+  if(use_height){
+    fn_graph +=Form("_height");
+  }else{
+    fn_graph +=Form("_int");
+  }
+
+  if(average_ratio){
+    fn_graph +=Form("_avgratio");
+  }else{
+    fn_graph +=Form("_indratio");
+  }
+
+  fn_graph += ".pdf";
+  
+  
+  // Directory and file name for measured Vdepl, initial depletion voltage, luminosity and Hamburg model output
+  TString dn_data(Form("%s/data/ST/Aging",diskVar));
+  TString fn_data(Form("%s/CCEScan.root",dn_data.Data()));
+  //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_hamburg(Form("%s/data/ST/Aging",diskVar));
+  TString fn_hamburg(Form("%s/Hamburg.root",dn_hamburg.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_hamburg = TFile::Open(fn_hamburg.Data());
+  
+
+  TVectorD* vp_fill  = (TVectorD*)f_fill->Get("v_fill");//CHECK
+  TVectorD* vp_time  = (TVectorD*)f_fill->Get("v_time");//CHECK: needs to be updated
+
+  
+  if (!vp_fill || !vp_time) {
+    Error("plotSector","Luminosity vectors are not available");
+    return EXIT_FAILURE;
+  }
+  
+  TVectorD v_fill    = *vp_fill;
+  TVectorD v_time    = *vp_time;
+  
+  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("plotSector","Depletion Voltage vectors are not available");
+    return EXIT_FAILURE;
+  }
+  
+  TVectorD v_sector  = *vp_sector;
+  TVectorD v_vdepl   = *vp_vdepl;
+  
+  TString vn_input = TString::Format("v_Hamburg_%s_%d",
+                                      lay.Data(),sector);
+  
+  if (radius>0.0) {
+    vn_input += Form("_r%d",(int)radius);
+  }
+  
+  TVectorD* vp_deltaV  = (TVectorD*)f_hamburg->Get(Form("%s_deltaV",vn_input.Data()));
+  TVectorD* vp_time_H  = (TVectorD*)f_hamburg->Get(Form("%s_time",  vn_input.Data()));
+  TVectorD* vp_hig     = (TVectorD*)f_hamburg->Get(Form("%s_hig",   vn_input.Data()));
+  TVectorD* vp_low     = (TVectorD*)f_hamburg->Get(Form("%s_low",   vn_input.Data()));
+  
+  if (!vp_deltaV || !vp_time_H || !vp_hig || !vp_low) {
+    Error("plotSector","Hamburg model vectors are not available");
+    return EXIT_FAILURE;
+  }
+  
+  TVectorD v_deltaV   = *vp_deltaV;
+  TVectorD v_time_H   = *vp_time_H;
+  TVectorD v_hig      = *vp_hig;
+  TVectorD v_low      = *vp_low;
+  
+  double vdeplProdPrec = 2.5;
+  
+  TVectorD v_zero(v_deltaV.GetNoElements());
+  
+  v_deltaV *= -1.0;
+  v_deltaV += v_vdepl(getClosestIndex(v_sector,sector));
+  
+  TVectorD v_hig_Draw(v_deltaV);
+  v_hig_Draw += v_hig;
+  v_hig_Draw += vdeplProdPrec;
+  
+  TVectorD v_low_Draw(v_deltaV);
+  v_low_Draw -= v_low;
+  v_low_Draw -= vdeplProdPrec;
+  
+  
+  std::vector<double> v_time_val;
+  std::vector<double> v_time_err;
+  std::vector<double> v_vdepl_val;
+  std::vector<double> v_vdepl_sta;
+  std::vector<double> v_vdepl_tot;
+  
+  // Extract for each fill the measurement and the corresponding time
+  // and fill them into a STL vector
+  for (int i=0; i<v_fill.GetNoElements(); ++i) {
+    double time_val   = v_time(i);
+    double time_err   = 0.0;
+    TString s_rad("");
+    TString vn_data(Form("%s/%s/%d/%d/v_volt_val%d",
+                         det.Data(),lay.Data(),sector,(int)v_fill(i),val));
+    if(radius>0.0) {
+      s_rad = Form("%d",(int)radius);
+      vn_data         += "_r"+s_rad;
+    }
+    if(use_height){
+      vn_data  +=Form("_height");
+    }else{
+      vn_data  +=Form("_int");
+    }
+
+    if(average_ratio){
+      vn_data  +=Form("_avgratio");
+    }else{
+      vn_data  +=Form("_indratio");
+    }
+    
+    TVectorD* vp_data = (TVectorD*)f_data->Get(vn_data.Data());
+
+    if (!vp_data){
+      Warning("plotSector","No vector found for %s",vn_data.Data());
+      continue;
+    }
+
+
+    // Exclude fill 1616 -> used for Vdepl normalisation ------> This is not clear to me, what are we normalizing?
+    if((int)v_fill(i)==1616){
+      continue;
+    }
+    
+    TVectorD v_data   = *vp_data;
+    double vdepl_val  = v_data(0);
+    double vdepl_sta  = v_data(1);
+    double vdepl_sys  = v_data(2);
+
+    Printf("        V_depl( %d ) : %8.4f +/- %8.4f(stat) +/-%8.4f (sys)",(int)v_fill(i),vdepl_val, vdepl_sta,vdepl_sys);//MIchele    
+    if (det.EqualTo("IT")) {
+      vdepl_sys = TMath::Sqrt(vdepl_sys*vdepl_sys+STTool::voltSyst_IT*STTool::voltSyst_IT);
+    }
+    if (det.EqualTo("TT")) {
+      vdepl_sys = TMath::Sqrt(vdepl_sys*vdepl_sys+STTool::voltSyst_TT*STTool::voltSyst_TT);
+    }
+    
+    if (vdepl_sta>20.0) {
+      Info("plotSector","Statistical uncertainty (%5.2f) for %s too high",vdepl_sta,vn_data.Data());
+      //continue;
+    }
+    v_time_val.push_back(time_val);
+    v_time_err.push_back(time_err);
+    v_vdepl_val.push_back(vdepl_val);
+    v_vdepl_sta.push_back(vdepl_sta);
+    v_vdepl_tot.push_back(TMath::Sqrt(vdepl_sta*vdepl_sta+vdepl_sys*vdepl_sys));
+    
+  }
+  
+  // Convert STL vectors into ROOT vectors
+  TVectorD vr_time_val     = getROOTVector<double>(v_time_val);
+  TVectorD vr_time_err     = getROOTVector<double>(v_time_err);
+  TVectorD vr_vdepl_val    = getROOTVector<double>(v_vdepl_val);
+  TVectorD vr_vdepl_sta    = getROOTVector<double>(v_vdepl_sta);
+  TVectorD vr_vdepl_tot    = getROOTVector<double>(v_vdepl_tot);
+  
+  if (v_time_val.size()==0) {
+    Warning("plotSector","For sector %d %s there is no data point!",sector,radius>0.0?Form("and %d mm",(int)radius):"");
+    return EXIT_SUCCESS;
+  }
+  
+  // Define graphs representing Hamburg model
+  TGraph* g_hamburg = new TGraph(v_time_H,v_deltaV);
+  TGraphAsymmErrors* ge_hamburg = new TGraphAsymmErrors(v_time_H,v_deltaV,v_zero,v_zero,v_low,v_hig);
+  TGraph* g_hamburg_hig = new TGraph(v_time_H,v_hig_Draw);
+  TGraph* g_hamburg_low = new TGraph(v_time_H,v_low_Draw);
+  
+  ge_hamburg->SetFillColor(kGray+2);
+  g_hamburg_hig->SetLineStyle(9);
+  g_hamburg_low->SetLineStyle(9);
+  
+  
+  TGraphErrors* ge_vdepl_sta  = new TGraphErrors(vr_time_val,vr_vdepl_val,
+                                                 vr_time_err,vr_vdepl_sta);
+  ge_vdepl_sta->SetLineColor(kRed);
+  ge_vdepl_sta->SetMarkerColor(kRed);
+  ge_vdepl_sta->SetMarkerSize(1);
+  
+  TGraphErrors* ge_vdepl_tot  = new TGraphErrors(vr_time_val,vr_vdepl_val,
+                                                 vr_time_err,vr_vdepl_tot);
+  ge_vdepl_tot->SetLineColor(kRed);
+  ge_vdepl_tot->SetMarkerColor(kRed);
+  ge_vdepl_tot->SetMarkerSize(1);
+  
+  TMultiGraph* mg_frame = new TMultiGraph("mg_frame","Frame multigraph");
+  
+  // Plot all graphs
+  TCanvas* c_frame = new TCanvas("c_frame","Frame Canvas",800,600);
+  mg_frame->Draw("");
+  
+  mg_frame->Add(ge_hamburg,"E3");
+  mg_frame->Add(g_hamburg_hig,"L");
+  mg_frame->Add(g_hamburg_low,"L");
+  mg_frame->Add(g_hamburg,"L");
+  mg_frame->Add(ge_vdepl_sta,"P");
+  mg_frame->Add(ge_vdepl_tot,"[]");
+  
+  mg_frame->Draw("A");
+  
+  // mg_frame->GetXaxis()->SetLimits(1262304000,1372636800);
+  double jan2010 = 1262304000.0;
+  double dec2015 = 1450103679.0;
+  double feb2016 = 1455926400.0;
+  double may2016 = 1463040975.0;
+  double sep2016 = 1475243957.0;
+  double nov2016 = 1480427833.0;
+  double feb2018 = 1519776000.0;
+  double jan2019 = 1546300800.0;
+  mg_frame->GetXaxis()->SetLimits(jan2010, feb2018);
+
+  mg_frame->SetMinimum(100.0);
+  mg_frame->SetMaximum(300.0);
+
+  if (det.EqualTo("IT")) {
+    mg_frame->SetMinimum(  0.0);
+    mg_frame->SetMaximum(200.0);
+  }
+  mg_frame->GetXaxis()->SetTimeDisplay(1);
+  mg_frame->GetXaxis()->SetTimeOffset(0.0);
+  mg_frame->GetXaxis()->SetTimeFormat("#splitline{%d/%m}{%Y}");
+  mg_frame->GetXaxis()->SetTitle("Date");
+  mg_frame->GetYaxis()->SetTitle("#it{V}_{depl} [V]");
+  mg_frame->GetXaxis()->SetLabelOffset(0.02);
+  mg_frame->GetXaxis()->SetTimeOffset(0.0);
+  mg_frame->GetYaxis()->SetTitleOffset(1.3);
+  mg_frame->GetXaxis()->SetTitleOffset(1.4);
+  
+  mg_frame->Draw("A");
+  
+  mg_frame->GetHistogram()->Draw("axissame");
+  
+  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());
+  }
+  
+  
+  
+  Printf("=========================================\n");
+
+  
+}
+
+