//
// Hamburg.C
// Script to calculate development
//
// 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/incStatROOT.h"
#include "../../include/basicROOTTools.h"
#include "../../include/basicRooFitTools.h"
#include "../../include/lhcbstyle.C"
#include "deplTool.h"
#include "hamburgTool.h"
// -d detector (TT/IT)
// -l layer (TTaX,TTaU,TTbV,TTbX, T[1-3]{X1,U,V,X2})
// -s read-out sector number
// -r consider only hit in a distance closer than r mm from the beam axis
int Hamburg(int argc, char* argv[]);
double rungeKutta4(HMap, double, double, double, double,double,bool);
double GetNc(HMap, double);
double GetNa(HMap, double, double, double, double, double);
double GetNy(HMap, double, double, double, double, double);
double conv_Neff_DeltaV(HMap);
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 = Hamburg(argc,argv);
Printf("End");
if (!runBatch) {
theApp->Run(!runBatch);
}
return exit;
}
// Executable method
int Hamburg(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");
// const char* homeVar = "/afs/cern.ch/user/e/egraveri/cmtuser/STMonitoring/STAging";
TString det("TT");
TString lay("TTaU");
double radius = -1000.0;
int sector = 2634;
for (int i=1; i<argc; i++) {
if (strcmp(argv[i],"-d")==0 && i+1<argc) {
det = argv[++i];
}else if (strcmp(argv[i],"-l")==0 && i+1<argc) {
lay = argv[++i];
}else if (strcmp(argv[i],"-r")==0 && i+1<argc) {
radius = std::atof(argv[++i]);
}else if (strcmp(argv[i],"-s")==0 && i+1<argc) {
sector = std::atoi(argv[++i]);
}
}
TString vn_output = TString::Format("v_Hamburg_%s_%d",
lay.Data(),sector);
if (radius>0.0) {
vn_output += Form("_r%d",(int)radius);
}
// Directory and file name for the output
TString dn_out(Form("%s/data/ST/Aging",
diskVar));
TString fn_out(Form("%s/Hamburg.root",
dn_out.Data()));
// Directory and file name for the temperature
TString dn_temp(Form("%s/data/temperature",homeVar));
TString fn_temp(Form("%s/temperature.root",dn_temp.Data()));
// Directory and file name for the luminosity
TString dn_lumi(Form("%s/data/lumi",homeVar));
//TString dn_lumi(Form("%s/data/db",homeVar));
TString fn_lumi(Form("%s/lumi.root",dn_lumi.Data()));
// // Directory and file name for the fluence
TString dn_sector(Form("%s/data/db",homeVar));
TString fn_sector(Form("%s/fluence.root",dn_sector.Data()));
TFile* f_temp = TFile::Open(fn_temp.Data());
TFile* f_lumi = TFile::Open(fn_lumi.Data());
TFile* f_sector = TFile::Open(fn_sector.Data());
if (!f_temp || !f_lumi || !f_sector) {
Error("Hamburg","Some files are not available");
return EXIT_FAILURE;
}
// Extract vectors to get beam periods
TVectorF* vp_fillTime = (TVectorF*)f_lumi->Get("v_fillTime_tot");
TVectorF* vp_beamEner = (TVectorF*)f_lumi->Get("v_beamEner_tot");
TVectorF* vp_stopTime = (TVectorF*)f_lumi->Get("v_stopTime_tot");
TVectorF* vp_peakLumi = (TVectorF*)f_lumi->Get("v_peakLumi_tot");
if (!vp_fillTime || !vp_beamEner || !vp_stopTime || !vp_peakLumi) {
Error("Hamburg","Luminosity vectors are not available");
return EXIT_FAILURE;
}
TVectorD v_fillTime = *vp_fillTime;
TVectorD v_beamEner = *vp_beamEner;
TVectorD v_stopTime = *vp_stopTime;
TVectorD v_peakLumi = *vp_peakLumi;
// Extract temperature measurements
TVectorD* vp_temp1 = (TVectorD*)f_temp->Get("v_TT_CT6_temp");
TVectorD* vp_temp2 = (TVectorD*)f_temp->Get("v_TT_AT6_temp");
TVectorD* vp_time1 = (TVectorD*)f_temp->Get("v_TT_CT6_time");
TVectorD* vp_time2 = (TVectorD*)f_temp->Get("v_TT_AT6_time");
if (det.EqualTo("IT")) {
vp_temp1 = (TVectorD*)f_temp->Get("v_TT_T3_temp");
vp_temp2 = (TVectorD*)f_temp->Get("v_TT_B3_temp");
vp_time1 = (TVectorD*)f_temp->Get("v_IT_T3_time");
vp_time2 = (TVectorD*)f_temp->Get("v_TT_B3_time");
}
if (!vp_temp1 || !vp_temp2 || !vp_time1 || !vp_time2) {
Error("Hamburg","Temperature vectors are not available");
return EXIT_FAILURE;
}
TVectorD v_temp1 = *vp_temp1;
TVectorD v_temp2 = *vp_temp2;
TVectorD v_time1 = *vp_time1;
TVectorD v_time2 = *vp_time2;
// Extract information about the flux in the corresponding layer
TString vn_sector("v_sector_");
TString vn_flux7_v("v_flux7_v_");
TString vn_flux7_e("v_flux7_e_");
TString vn_flux8_v("v_flux8_v_");
TString vn_flux8_e("v_flux8_e_");
TString vn_flux13_v("v_flux13_v_");
TString vn_flux13_e("v_flux13_e_");
vn_sector += lay.Data();
vn_flux7_v += lay.Data();
vn_flux7_e += lay.Data();
vn_flux8_v += lay.Data();
vn_flux8_e += lay.Data();
vn_flux13_v += lay.Data();
vn_flux13_e += lay.Data();
// Add radius postfix
if (det.EqualTo("IT")) {
}else{
if (radius>0.0) {
vn_sector += Form("_r%d",(int)radius);
vn_flux7_v += Form("_r%d",(int)radius);
vn_flux7_e += Form("_r%d",(int)radius);
vn_flux8_v += Form("_r%d",(int)radius);
vn_flux8_e += Form("_r%d",(int)radius);
vn_flux13_v += Form("_r%d",(int)radius);
vn_flux13_e += Form("_r%d",(int)radius);
}
}
TVectorD* vp_sector = (TVectorD*)f_sector->Get(vn_sector.Data());
TVectorD* vp_flux7_v = (TVectorD*)f_sector->Get(vn_flux7_v.Data());
TVectorD* vp_flux7_e = (TVectorD*)f_sector->Get(vn_flux7_e.Data());
TVectorD* vp_flux8_v = (TVectorD*)f_sector->Get(vn_flux8_v.Data());
TVectorD* vp_flux8_e = (TVectorD*)f_sector->Get(vn_flux8_e.Data());
TVectorD* vp_flux13_v = (TVectorD*)f_sector->Get(vn_flux13_v.Data());
TVectorD* vp_flux13_e = (TVectorD*)f_sector->Get(vn_flux13_e.Data());
if (!vp_sector || !vp_flux7_v || !vp_flux7_e || !vp_flux8_v || !vp_flux8_e || !vp_flux13_v || !vp_flux13_e) {
Error("Hamburg","Flux vectors are not available");
return EXIT_FAILURE;
}
TVectorD v_sector = *vp_sector;
TVectorD v_flux7_v = *vp_flux7_v;
TVectorD v_flux7_e = *vp_flux7_e;
TVectorD v_flux8_v = *vp_flux8_v;
TVectorD v_flux8_e = *vp_flux8_e;
TVectorD v_flux13_v = *vp_flux13_v;
TVectorD v_flux13_e = *vp_flux13_e;
// Define number of calculated pathes
int nErr = 1000;
// Parameters for Hamburg model
HMap para_base = STHamburg::hamburgPara;
HMap unce_base = STHamburg::hamburgUnce;
if (det.EqualTo("IT")) {
para_base["d"] = 410.0e-4; // So far only 410 µm sensors considered
}
std::vector<HMap> v_para;
std::vector<double> v_tempUnce;
std::vector<double> v_fluxUnce7;
std::vector<double> v_fluxUnce8;
std::vector<double> v_fluxUnce13;
std::vector<double> v_deltaV;
std::vector<TVectorD> v_N;
std::vector<double> v_time;
std::vector<double> v_val;
std::vector<double> v_low;
std::vector<double> v_hig;
std::vector<double> v_flux;
int i_flux = getClosestIndex(v_sector,(double)sector);
double flux7_v = v_flux7_v(i_flux);
double flux7_e = v_flux7_e(i_flux);
double flux8_v = v_flux8_v(i_flux);
double flux8_e = v_flux8_e(i_flux);
double flux13_v = v_flux13_v(i_flux);
double flux13_e = v_flux13_e(i_flux);
// Temperature uncertainty of 2 degree
double temp_e = 2.0;
Info("Hamburg","Filling parameter vectors");
for (int i=0; i<nErr; i++) {
HMap temp_para = para_base;
// Define parameter for each path
if (i==0) {
v_para.push_back(para_base);
v_fluxUnce7.push_back(0.0);
v_fluxUnce8.push_back(0.0);
v_fluxUnce13.push_back(0.0);
v_tempUnce.push_back(0.0);
}else{
// Make gaussian variation of parameters
std::map<std::string, double>::iterator it = temp_para.begin();
for(;it!=temp_para.end(); ++it){
it->second += gRandom->Gaus()*unce_base[it->first];
}
v_para.push_back(temp_para);
v_fluxUnce7.push_back(flux7_e*gRandom->Gaus());
v_fluxUnce8.push_back(flux8_e*gRandom->Gaus());
v_fluxUnce13.push_back(flux13_e*gRandom->Gaus());
v_tempUnce.push_back(temp_e*(2.0*gRandom->Uniform()-1.0));
}
v_deltaV.push_back(0.0);
v_flux.push_back(0.0);
v_N.push_back(TVectorD(3));
}
double flux = 0.0;
double dFlux = 0.0;
int i_temp1 = 0;
int i_temp2 = 0;
int i_lumi = 0;
double tStart = 1262325600; // 1-1-2010
// double tStop = 1372654800; // 1-7-2013
double tStop = 1448928000; // 1-12-2015
double tRun = tStart;
double tFillEnd = tStart;
double tStep = 60000.0;
time_t raw_tStart = (long)tStart;
time_t raw_tStop = (long)tStop;
Info("Hamburg","Starting on: %20s",ctime(&raw_tStart));
Info("Hamburg","Stopping on: %20s",ctime(&raw_tStop ));
Info("Hamburg","Steps are %6.2f sec",tStep);
Info("Hamburg","Start running");
// Run between the start and the stop time and check if the index for temperature and luminosity needs to be increased
while (tRun<tStop) {
double temp = 273.15;
while (i_temp1<v_time1.GetNoElements() && tRun>v_time1(i_temp1)) {
++i_temp1;
}
while (i_temp2<v_time2.GetNoElements() && tRun>v_time2(i_temp2)) {
++i_temp2;
}
temp = (v_temp1(i_temp1)+v_temp2(i_temp2))/2.0+273.15;
if(i_lumi<v_fillTime.GetNoElements() && tRun>v_fillTime(i_lumi)) {
tFillEnd = v_stopTime(i_lumi);
}
while (i_lumi<v_fillTime.GetNoElements() && tRun>v_stopTime(i_lumi)){
++i_lumi;
}
// Calculate for each path the evolution of the depletion voltage
for (int i=0; i<v_deltaV.size(); i++) {
temp += v_tempUnce[i];
dFlux = 0.0;
if (i_lumi<v_fillTime.GetNoElements()) {
if (tRun<v_fillTime(i_lumi)) {
dFlux = 0.0;
}else{
if (v_beamEner(i_lumi)==7.0) {
dFlux = (flux7_v+v_fluxUnce7[i])*v_peakLumi(i_lumi)*STTool::FLUKAConvFac;
}else if(v_beamEner(i_lumi)==8.0) {
dFlux = (flux8_v+v_fluxUnce8[i])*v_peakLumi(i_lumi)*STTool::FLUKAConvFac;
}else if(v_beamEner(i_lumi)==13.0) {
dFlux = (flux13_v+v_fluxUnce13[i])*v_peakLumi(i_lumi)*STTool::FLUKAConvFac;
}else{
dFlux = 0.0;
}
}
}
// Calculate change in the effective doping concentration
v_N[i](0) = GetNc(v_para[i],v_flux[i]);
v_N[i](1) = GetNa(v_para[i],v_N[i](1),v_flux[i],dFlux,tStep,temp);
v_N[i](2) = GetNy(v_para[i],v_N[i](2),v_flux[i],dFlux,tStep,temp);
// Convert it to a change in the depletion voltage
v_deltaV[i] = v_N[i].Sum()*conv_Neff_DeltaV(v_para[i]);
v_flux[i] += dFlux*tStep;
// Reset temperature
temp -= v_tempUnce[i];
}
v_time.push_back(tRun);
v_val.push_back(v_deltaV[0]);
std::vector<double> v_temp = v_deltaV;
std::sort(v_temp.begin(),v_temp.end());
// Extract 68% confidence interval
v_low.push_back(TMath::Max(0.0,v_deltaV[0]-v_temp[(int)(0.158*v_temp.size())]));
v_hig.push_back(TMath::Max(0.0,v_temp[(int)(0.842*v_temp.size())]-v_deltaV[0]));
Printf("Delta V = (%6.2f+%6.2f/-%6.2f) V",
v_deltaV[0],
TMath::Max(0.0,v_temp[(int)(0.842*v_temp.size())]-v_deltaV[0]),
TMath::Max(0.0,v_deltaV[0]-v_temp[(int)(0.158*v_temp.size())]));
tRun+=tStep;
time_t raw_tRun = (long)tRun;
Info("Hamburg","Date/Time: %s",ctime(&raw_tRun));
}
// Reduce the number of data points to be plotted
Info("Hamburg","Reduce number of data points");
std::vector<double> v_val_red;
std::vector<double> v_time_red;
std::vector<double> v_hig_red;
std::vector<double> v_low_red;
int redFac = (int)(86400/tStep);
Printf("%d",v_val.size());
for (int i=0; i<v_val.size(); i++) {
Printf("%d",i);
if(i%redFac==0){
v_val_red.push_back(v_val[i]);
v_time_red.push_back(v_time[i]);
v_hig_red.push_back(v_hig[i]);
v_low_red.push_back(v_low[i]);
}
}
TVectorD vr_val = getROOTVector<double>(v_val_red);
TVectorD vr_time = getROOTVector<double>(v_time_red);
TVectorD vr_hig = getROOTVector<double>(v_hig_red);
TVectorD vr_low = getROOTVector<double>(v_low_red);
// Store result
TFile* f_out = new TFile(fn_out.Data(),"UPDATE");
f_out->Delete(Form("%s_deltaV;*",vn_output.Data()));
f_out->Delete(Form("%s_time;*" ,vn_output.Data()));
f_out->Delete(Form("%s_hig;*" ,vn_output.Data()));
f_out->Delete(Form("%s_low;*" ,vn_output.Data()));
vr_val.Write(Form("%s_deltaV",vn_output.Data()));
vr_time.Write(Form("%s_time",vn_output.Data()));
vr_hig.Write(Form("%s_hig",vn_output.Data()));
vr_low.Write(Form("%s_low",vn_output.Data()));
f_out->Close();
Info("Hamburg","Data written to %s",fn_out.Data());
Printf("=========================================\n");
return EXIT_SUCCESS;
}
// Runge-Kutta 4th order algorithm for Hamburg model (not used)
double rungeKutta4(HMap para, double flux, double dFlux, double temp, double dT,double T_noBeam,bool isBeam){
double dNc = 0.0;
double dNa = 0.0;
double dNy = 0.0;
double tau_a = 1.0/(para["k0a"]*TMath::Exp(-para["Eaa"]/(para["kB"]*temp)));
double tau_y = 1.0/(para["k0y"]*TMath::Exp(-para["Ey" ]/(para["kB"]*temp)));
Printf("tau_a: %f,tau_y: %f",tau_a,tau_y);
Printf("flux: %f, dFlux: %f, t: %f",flux,dFlux,T_noBeam);
dNc += (para["gc"]+para["Nc"]*para["c"]*TMath::Exp(-para["c"]*flux))*dFlux*dT/2.0;
dNa += (flux*para["ga"]*-1.0/tau_a+para["ga"]*dFlux)*TMath::Exp(-T_noBeam/tau_a)*dT/2.0;
dNy += (dFlux*para["gy"]*(1.0-1.0/(1.0+T_noBeam/tau_y))+
flux*para["gy"]/TMath::Power(1.0+T_noBeam/tau_y,2)/tau_y)*dT/2.0;
flux += dFlux*dT/2.0;
dNc += (para["gc"]+para["Nc"]*para["c"]*TMath::Exp(-para["c"]*flux))*dFlux*dT/2.0;
dNa += (flux*para["ga"]*-1.0/tau_a+para["ga"]*dFlux)*TMath::Exp(-T_noBeam/tau_a)*dT/2.0;
dNy += (dFlux*para["gy"]*(1.0-1.0/(1.0+T_noBeam/tau_y))+
flux*para["gy"]/TMath::Power(1.0+T_noBeam/tau_y,2)/tau_y)*dT/2.0;
Printf("dNc: %f,dNa: %f,dNy: %f",dNc,dNa,dNy);
return dNc+dNa+dNy;
}
// Calculate stable damage part
double GetNc(HMap para, double flux){
return para["Nc"]*(1.0-TMath::Exp(-para["c"]*flux))+para["gc"]*flux;
}
// Calculate change in the annealing part
double GetNa(HMap para, double Na, double flux, double dFlux, double dT, double temp){
double tau_a = 1.0/(para["k0a"]*TMath::Exp(-para["Eaa"]/(para["kB"]*temp)));
Na += (para["ga"]*dFlux-Na/tau_a)*dT/2.0;
flux += dFlux*dT/2.0;
Na += (para["ga"]*dFlux-Na/tau_a)*dT/2.0;
return Na;
}
// Calculate change in the reverse annealing part
double GetNy(HMap para, double Ny, double flux, double dFlux, double dT, double temp){
double tau_y = 1.0/(para["k0y"]*TMath::Exp(-para["Ey" ]/(para["kB"]*temp)));
Ny += (TMath::Power((flux==0.0?0.0:Ny/(flux*para["gy"]))-1,2)*flux*para["gy"]/tau_y)*dT/2.0;
flux += dFlux*dT/2.0;
Ny += (TMath::Power((flux==0.0?0.0:Ny/(flux*para["gy"]))-1,2)*flux*para["gy"]/tau_y)*dT/2.0;
return Ny;
}
// Convert change in the effective doping concentration to change in the depletion voltage
double conv_Neff_DeltaV(HMap para){
return para["q"]*para["d"]*para["d"]/(2.0*para["eps"]*para["eps0"]);
}
elena