diff --git a/InnerDetector/InDetCalibAlgs/PixelCalibAlgs/RadDamage/HamburgModel/README.md b/InnerDetector/InDetCalibAlgs/PixelCalibAlgs/RadDamage/HamburgModel/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..6271ead6172e995fce03788e705e4d198b51ec39
--- /dev/null
+++ b/InnerDetector/InDetCalibAlgs/PixelCalibAlgs/RadDamage/HamburgModel/README.md
@@ -0,0 +1,33 @@
+Original author: Julien Beyer, 2019 (https://gitlab.cern.ch/jbeyer/RadDamSim)
+
+The code in data_processing.cpp is a c++ implementation of the Hamburg model of radiation damage. It calculates the effective doping concentration based on an introduction of stable/decaying acceptors and the removal of donors due to irradiation. Moreover, it handles the annealing processes depending in the temperature. According to the effective doping concentration, the depletion voltage is calculated.
+Besides the depletion voltage, the leakage current (inclusive alpha parameter) is calculated based on the Hamburg model. Again, leakage current increase due to irradiation and decrease due to annealing is handled. Two different options for the temperature averaging are available (use_CMS_paper_alpha0).
+
+In order to simulate all of these properties, a profile of irradiation dose and temperature is necessary. The profile has to have the following format:
+
+duration (in seconds!!!!!!!!!, integer) / temperature in K (float) / dose rate in neq/cm2s (integer !)
+
+Different examples are available in the temp_rad_data folder. Besides that, a program to create new profiles from two given temperature and irradiation profiles is available in this folder (help included in the cpp file).
+
+As a prerequisite make sure that Root (for plotting, tested with version 6.06) and Boost (for Date handling, tested with version 1.62) as well as a c++ 11 compatible compiler are installed. Now, the program can be compiled with:
+
+g++ data_processing.cpp -I /path/to/boost/boost_1_62_0 -Wall -std=c++11 -o output `root-config --cflags` `root-config --libs`
+
+After compilation the program can be called with:
+
+./output temp_rad_data/my_profile.txt
+
+By default, a variety of pdf files is generated and all corresponding root histograms are stored in a single root file. Furthermore, by default, the pdf files will be overwritten and the histograms will be added to the existing root file - this can be changed in the first lines of the code, a unique naming scheme based on the system clock is available!
+
+Some variables might be changed in the first part of the code, apart from that, no adaptions should be necessary.
+
+Attention:
+1.duration in the profile file HAS to be dividable by global timestep (parameter in the code, default=1) ->Integer result required for the ratio of duration and global time step
+
+2.profile file may not have additional (even/especially empty) lines after the last line or in the beginning
+
+3.the value for the global_layer_conversion has to be changed in order to fit to another layer (parameter transforms luminosity in neq) eg: B-Layer: 2.5e12, IBL 6e12
+
+4.due to some computational reasons the final luminosity might not be correct. The fluence (in neq/cm2) however, should be exact.
+
+5. to change from one detector type to another (e.g. for a different layer in the detector), change mainly 3 things: thickness (200 for IBL and 250 for PIXEL) and global_layer_conversion (global values) and Ndonor (different for PIXEL and IBL)
diff --git a/InnerDetector/InDetCalibAlgs/PixelCalibAlgs/RadDamage/HamburgModel/data_processing.cpp b/InnerDetector/InDetCalibAlgs/PixelCalibAlgs/RadDamage/HamburgModel/data_processing.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..80abb6405a800a7275323157c5ec8593e5bdfa8d
--- /dev/null
+++ b/InnerDetector/InDetCalibAlgs/PixelCalibAlgs/RadDamage/HamburgModel/data_processing.cpp
@@ -0,0 +1,770 @@
+/*
+Compile Code with: g++ data_processing.cpp -I /path/to/boost/boost_1_62_0 -Wall -std=c++11 -o output `root-config --cflags` `root-config --libs`
+
+Requirements: Root (for plotting) and Boost (for Date handling) have to be installed
+
+Call program with 1 Arguments: 1. name of the file with the temperature/radiation profile (you may want to uncomment the gA,Y,C to be read in again as additional parameters, around line 686)
+
+Radiation and temperature profile have to fulfil the following scheme: duration (in seconds!!!!!!!!!) / temperature in K / irradiation in neq/cm2s (ALL have to be integer !)
+
+Attention: timesteps in the profile file HAVE to be devideable by global timestep ->Integer result required
+ profile files may not have additional (even/especially empty) lines after the last line
+
+Attention 2: the value for the global_layer_conversion has to be changed in order to fit to another layer (parameter transforms luminosity in neq) eg: B-Layer: 2.5e12, IBL 6e12
+
+Attention 3: to change from one detector type to another, change mainly 3 things: thickness and global_layer_conversion (global values) and Ndonor
+*/
+
+
+#include <iostream> // Basic IO
+#include <fstream> // Read and write files
+#include <cstdlib> // convert arguments of function call into int, float,...
+#include <string>
+#include <vector>
+#include <math.h> // for fabs()
+#include <sstream> // to get string into stream
+//#include <TH2F.h> // root stuff for file reading
+#include <TFile.h> // more root stuff
+//#include <TH2D.h> // root stuff for file reading
+#include <TCanvas.h>
+#include <TROOT.h>
+#include <TGraphErrors.h>
+#include <TVector.h>
+#include <TStyle.h>
+#include <TAxis.h>
+#include <TF1.h>
+#include <TMath.h>
+#include <TLine.h>
+#include <TH1.h>
+#include <cmath>
+#include "TDatime.h"
+#include <time.h>
+#include <ctime>
+#include <boost/date_time/gregorian/gregorian.hpp>
+
+using namespace std;
+using namespace boost::gregorian;
+
+/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+bool use_CMS_paper_alpha0 = false; //if false: temperature averaging for leakage current calculation a la Moll (theta function for time scale change), if true: use CMS approach of temperature average (attention: constants can only be changed directly in the code in this case)
+
+bool debug=false; // additional debugging console outputs
+bool plot_pdf=true; // plots are by default saved as root file, but can also be directly exported as pdf
+bool overwrite_files=true; // false = append internal clock identifier to file name, false = overwrite pdfs, append root files // files are by default created with a unique name depending on the internal clock and you can switch it of here such that pdf files are overwritten and the root file is extended (old plots will stay in the root file as well)
+double timestep=1; // step size is 1 second, do not change this!
+double donorremovalfraction=0.99; // fraction of donors which can be removed by donor removal
+double userTref=273.15; // set a reference temperature for the volume corrected leakage current plot (it will only effect this one plot!) Now: implemented!
+double bandGap=1.21; // eV used for scaling temperatures
+
+string output_file_name = "simulation_results"; // set unique file name for each simulation
+
+date d(2015,May,1); // IBL //set a date for the plots to begin (to be correct it has to be equal to the beginning of your (!) temp/irr profile)
+//date d(2011,Feb,11); // PIXEL
+
+
+const double Ndonor_0 = 1.7e12; // IBL // initial donor concentration (right now the code only works for originally n-type bulk material!)
+//const double Ndonor_0 = 1.4e12; // B-Layer Layer1/2 Disks
+
+double thickness=200; // IBL // sensor thickness
+//double thickness=250; // B-Layer Layer1/2 Disks
+
+//double global_layer_conversion=0.92e12; // Layer 1 //conversion factor from luminosity to neq (1fb-1=2.3e12neq/cm2) - is used only for computation of total luminosity, will be wrong if there are multiple conversion factor in the original profiles as for example due to different center of mass energies
+//double global_layer_conversion=1.1e12; // L1 average
+//double global_layer_conversion=0.571e12; // Layer 2
+//double global_layer_conversion=0.7e12; // Layer 2 average
+//double global_layer_conversion=0.582e12; // Disks
+double global_layer_conversion=6.262e12; // IBL
+//double global_layer_conversion=2.8e12;//2.929e12; // B-Layer
+
+int limit_file_size_input = 0; // how many lines should be read from the profile, 0 for everything
+
+float DoseRateScaling=1.; // this parameter is multiplied to all fluence rates from the input profile
+
+/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////Do not change things below this line without thinking carefully!
+
+//Definition of Hamburg model depletion voltage functions
+
+TF1 *Naccept_rev_TF1 = new TF1("Naccept_rev_TF1","[0]*[1]*(1-TMath::Exp(-[2]*x))/[2] + [3]*TMath::Exp(-[2]*x)",0,10000000);
+TF1 *Naccept_const_TF1 = new TF1("Naccept_const_TF1","[0] * [1]*x",0,10000000);
+TF1 *Nneutrals_rev_TF1 = new TF1("Ndonor_rev_TF1","[0]*[1]*(1-TMath::Exp(-[2]*x))/[2] + [3]*TMath::Exp(-[2]*x)",0,10000000);
+TF1 *Ndonor_const_TF1 = new TF1("Ndonor_const_TF1","-[0]*(1-TMath::Exp(-[1]*[2]*x))",0,10000000);
+
+TF1 *Naccept_rev_TF1_approx = new TF1("Naccept_rev_TF1","[0]*[1]*x + [3]*TMath::Exp(-[2]*x)",0,10000000);
+TF1 *Ndonor_neutrals_rev_TF1_approx = new TF1("Ndonor_rev_TF1","[0]*[1]*x + [3]*TMath::Exp(-[2]*x)",0,10000000);
+
+TF1 *Ndonor_TF1 = new TF1("Ndonor_TF1","[0]*[1]/[2] * ([2]*x+TMath::Exp(-[2]*x)-1) +[3]*(1-TMath::Exp(-[2]*x)) ",0,10000000);
+
+/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+struct DataElement
+{
+ int duration;
+ float temperature;
+ long int dose_rate;
+};
+
+
+/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+//calculate depletion voltage from given effective doping concentration
+
+double NtoV_function(double doping_conc)
+{
+ return 1.6021766208e-13/(2*11.68*8.854187817)*fabs(doping_conc)*thickness*thickness; // q/2*epsilon*epsilon_0 * Neff*D^2
+}
+
+/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+
+class Annealing_constants { // Class to store a set of annealing specific constants and compute temperature dependent values
+private:
+ double gA;
+ double gY;
+ double gC;
+ double ka;
+ double ky1;
+ double Ea;
+ double Ey;
+ double cc;
+
+public:
+ Annealing_constants(double,double,double,double,double,double,double,double);
+ double get_gA(int) const;
+ double get_gY(int) const;
+ double get_gC(int) const;
+ double get_ka(int) const;
+ double get_ky1(int) const;
+ double get_Ea() const;
+ double get_Ey() const;
+ double get_cc() const;
+};
+
+Annealing_constants::Annealing_constants(double a ,double b ,double c ,double d ,double e ,double f ,double g , double h)
+{
+ gA=a ;
+ gY=b ;
+ gC=c ;
+ ka=d ;
+ ky1=e ;
+ Ea=f ;
+ Ey=g ;
+ cc=h ;
+}
+
+double Annealing_constants::get_gA(int T) const
+{
+ return gA;
+}
+
+double Annealing_constants::get_gY(int T) const
+{
+ return gY;
+}
+
+double Annealing_constants::get_gC(int T) const
+{
+ return gC;
+}
+
+double Annealing_constants::get_ka(int T) const
+{
+ return ka*exp(-Ea/(8.6173303e-5*(double)T));
+ //return ka;
+}
+
+double Annealing_constants::get_ky1(int T) const
+{
+ return ky1*exp(-Ey/(8.6173303e-5*(double)T));
+ //return ky1;
+}
+
+double Annealing_constants::get_Ea() const
+{
+ return Ea;
+}
+
+double Annealing_constants::get_Ey() const
+{
+ return Ey;
+}
+
+double Annealing_constants::get_cc() const
+{
+ return cc;
+}
+
+
+/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+// leakage current constants will be stored in this struct
+
+struct leakage_current_consts
+{
+ double alpha_1;
+ double alpha_0_star;
+ double beta;
+ double k01;
+ double E1;
+ double E1_star;
+};
+
+/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+// all atributes of a sensor (like irradiation, doping concentration, leakage currents, etc) will be stored in this class
+
+class Sensor //Class to handle Sensors
+{
+private: //Content of a sensor
+ double Nacceptor;
+ double Ndonor;
+ double Ndonor_const;
+ double Nacceptor_reversible;
+ double Nneutral_reversible;
+ double Nacceptor_stable_constdamage;
+ double Ndonor_stable_donorremoval;
+ double Nacceptor_stable_reverseannealing;
+
+ int Temperature;
+ vector<vector<double>> leakage_current_alpha;
+ vector<vector<double>> time_history;
+ vector<double> G_i;
+ vector<double> leakage_current;
+ double volume;
+ vector<double> alpha_vec;
+ vector<double> powerconsumption;
+ vector<double> inverse_weighted_temp;
+
+
+public:
+ Sensor(double, double, int, double, double); //Constructor
+ Sensor();
+ double get_Neff() const;
+ void set_Nacceptor(double);
+ void set_Ndonor(double);
+ double get_Nacceptor() const;
+ double get_Ndonor() const;
+ void set_temp(int);
+ double get_temp() const;
+ void irradiate(leakage_current_consts,Annealing_constants, long int, float, long double);
+ vector<double> get_G_i() const;
+ vector<double> get_leakage_current() const;
+ vector<double> get_alpha_vec() const;
+ vector<double> get_powerconsumption() const;
+};
+
+Sensor::Sensor( double a, double b, int c, double d, double e) //Constructordefinition
+{
+ Nacceptor=a;
+ Ndonor=b;
+ Temperature=c;
+ Nacceptor_reversible=d;
+ Nneutral_reversible=e;
+}
+
+Sensor::Sensor() //Constructordefinition
+{
+ Nacceptor=0;
+ Ndonor=Ndonor_0;
+ Temperature=0;
+ Nneutral_reversible=0;
+ Nacceptor_reversible=0;
+ Nneutral_reversible=0;
+ volume= 0.135;
+ Nacceptor_stable_constdamage=0;
+ Ndonor_stable_donorremoval=donorremovalfraction*fabs(Nacceptor-Ndonor);
+ Ndonor_const=Ndonor-Ndonor_stable_donorremoval;
+ Nacceptor_stable_reverseannealing=0;
+}
+
+double Sensor::get_Neff() const
+{
+ return Nacceptor-Ndonor;
+}
+
+double Sensor::get_temp() const
+{
+ return Temperature;
+}
+
+void Sensor::set_temp(int value)
+{
+ Temperature=value;
+ return;
+}
+
+double Sensor::get_Nacceptor() const
+{
+ return Nacceptor;
+}
+
+double Sensor::get_Ndonor() const
+{
+ return Ndonor;
+}
+
+void Sensor::set_Nacceptor(double new_value)
+{
+ Nacceptor=new_value;
+ return;
+}
+
+void Sensor::set_Ndonor(double new_value)
+{
+ Ndonor=new_value;
+ return;
+}
+
+vector<double> Sensor::get_G_i() const
+{
+ return G_i;
+}
+
+vector<double> Sensor::get_leakage_current() const
+{
+ return leakage_current;
+}
+
+vector<double> Sensor::get_alpha_vec() const
+{
+ return alpha_vec;
+}
+
+vector<double> Sensor::get_powerconsumption() const
+{
+ return powerconsumption;
+}
+
+void Sensor::irradiate(leakage_current_consts leconsts, Annealing_constants constants,long int phi, float t, long double totalDose)
+{
+ //t=t*3600; //conversion from hours to seconds
+ double a=1e-30;
+//calculating the effective doping concentration
+
+if(debug) cout << t<<" "<< constants.get_gA(Temperature)<<" "<<phi<<" "<<constants.get_ka(Temperature)<<" "<<Nacceptor<<" "<< constants.get_gC(Temperature)<<" "<< constants.get_gY(Temperature)<<" "<< constants.get_ky1(Temperature) <<" "<<Ndonor <<endl;
+
+ Naccept_rev_TF1->SetParameter(0,constants.get_gA(Temperature));
+ Naccept_rev_TF1->SetParameter(1,phi);
+ Naccept_rev_TF1->SetParameter(2,constants.get_ka(Temperature));
+ Naccept_rev_TF1->SetParameter(3,Nacceptor_reversible);
+
+ Naccept_rev_TF1_approx->SetParameter(0,constants.get_gA(Temperature));
+ Naccept_rev_TF1_approx->SetParameter(1,phi);
+ Naccept_rev_TF1_approx->SetParameter(2,constants.get_ka(Temperature));
+ Naccept_rev_TF1_approx->SetParameter(3,Nacceptor_reversible);
+
+ Naccept_const_TF1->SetParameter(0,constants.get_gC(Temperature));
+ Naccept_const_TF1->SetParameter(1,phi);
+
+ Nneutrals_rev_TF1->SetParameter(0,constants.get_gY(Temperature));
+ Nneutrals_rev_TF1->SetParameter(1,phi);
+ Nneutrals_rev_TF1->SetParameter(2,constants.get_ky1(Temperature));
+ Nneutrals_rev_TF1->SetParameter(3,Nneutral_reversible);
+
+ Ndonor_neutrals_rev_TF1_approx->SetParameter(0,constants.get_gY(Temperature));
+ Ndonor_neutrals_rev_TF1_approx->SetParameter(1,phi);
+ Ndonor_neutrals_rev_TF1_approx->SetParameter(2,constants.get_ky1(Temperature));
+ Ndonor_neutrals_rev_TF1_approx->SetParameter(3,Nneutral_reversible);
+
+ Ndonor_const_TF1->SetParameter(0,Ndonor_stable_donorremoval);
+ Ndonor_const_TF1->SetParameter(1,constants.get_cc());
+ Ndonor_const_TF1->SetParameter(2,phi);
+
+ Ndonor_TF1->SetParameter(0,constants.get_gY(Temperature));
+ Ndonor_TF1->SetParameter(1,phi);
+ Ndonor_TF1->SetParameter(2,constants.get_ky1(Temperature));
+ Ndonor_TF1->SetParameter(3,Nneutral_reversible);
+
+ if(constants.get_ka(Temperature)>a && constants.get_ky1(Temperature)>a)
+ {
+ Nacceptor_reversible = Naccept_rev_TF1->Eval(t);
+ Nacceptor_stable_constdamage += Naccept_const_TF1->Eval(t);
+ Nneutral_reversible = Nneutrals_rev_TF1->Eval(t);
+ Ndonor_stable_donorremoval += Ndonor_const_TF1->Eval(t);
+ Nacceptor_stable_reverseannealing += Ndonor_TF1->Eval(t);
+ }
+ else
+ {
+ cout << "Potential numerical problem due to ultra low temp and therby caused very small ky1 and ka values. Using approach in order to perform calculation. In general, no problem!"<<endl;
+
+ Nacceptor_reversible = Naccept_rev_TF1_approx->Eval(t);
+ Nacceptor_stable_constdamage += Naccept_const_TF1->Eval(t);
+ Nneutral_reversible = Ndonor_neutrals_rev_TF1_approx->Eval(t);
+ Ndonor_stable_donorremoval += Ndonor_const_TF1->Eval(t);
+ Nacceptor_stable_reverseannealing += Ndonor_TF1->Eval(t);
+ }
+
+ Nacceptor = Nacceptor_reversible + Nacceptor_stable_constdamage + Nacceptor_stable_reverseannealing;
+ Ndonor = Ndonor_const + Ndonor_stable_donorremoval;
+
+//calculating the leackage current in the following part
+
+ vector<double> tmp(3);
+ vector<double> tmp2(2);
+ vector<double> tmp3(2);
+
+ tmp3.at(0)=0;
+ tmp3.at(1)=0;
+
+//////////////////////////////
+
+ double G_i_tmp=0; // G_i_tmp needs to be outside of the if clause to be existend also afterwards when it is pushed back in the G_i vector, can be put back to pos1 after finished with the if clauses
+
+ if(use_CMS_paper_alpha0)
+ {
+ double CMS_alpha0_c1=-8.9e-17;
+ double CMS_alpha0_c2=4.6e-14;
+ double CMS_beta=2.9e-18;
+
+ double inverse_weighted_temp_tmp=0;
+ inverse_weighted_temp_tmp = t/(double)Temperature;
+ inverse_weighted_temp.push_back(inverse_weighted_temp_tmp);
+
+ tmp.at(0)=phi*t*leconsts.alpha_1;
+ tmp.at(1)=phi*t; //needs to be updated at every itaration in the while loop since the averaged temperature needs to be recalculated
+ tmp.at(2)=phi*t*CMS_beta;
+
+ leakage_current_alpha.push_back(tmp);
+
+ tmp2.at(0)=t*leconsts.k01*exp(-leconsts.E1/(8.6173303e-5*(double)Temperature)); // time comes in seconds and k01 is seconds-1 so units are fine here
+ tmp2.at(1)=t / 60.0; // t/60 to convert time from seconds to minutes since this is the proposed unit here
+
+ if(time_history.size() > 0.1) tmp3=time_history.at(time_history.size()-1);
+
+ tmp2.at(0)+=tmp3.at(0);
+ tmp2.at(1)+=tmp3.at(1);
+
+ time_history.push_back(tmp2);
+
+
+ int i=0;
+ double temperature_average_i=0;
+
+ while(i<leakage_current_alpha.size())
+ {
+ temperature_average_i=0;
+
+ for(int j=i; j<leakage_current_alpha.size(); j++) //calculate the average inverse temperature weighted with the individual time from the moment of irradiation (i) to now (leakage_current_alpha.size())
+ {
+ temperature_average_i += inverse_weighted_temp.at(j);
+ }
+
+ if(i>=1) //for i=0 the time is not the difference but just the total time
+ {
+ temperature_average_i /= (double)(time_history.back().at(1)*60.0-time_history.at(i-1).at(1)*60.0);
+ G_i_tmp += leakage_current_alpha.at(i).at(0) * exp(-(time_history.back().at(0)-time_history.at(i-1).at(0))) + leakage_current_alpha.at(i).at(1)*(CMS_alpha0_c1+CMS_alpha0_c2*temperature_average_i) - leakage_current_alpha.at(i).at(2) * log(time_history.back().at(1)-time_history.at(i-1).at(1) ) ;
+ }
+ else
+ {
+ temperature_average_i /= (double)(time_history.back().at(1)*60.0);
+ G_i_tmp += leakage_current_alpha.at(i).at(0) * exp(-time_history.back().at(0)) + leakage_current_alpha.at(i).at(1)*(CMS_alpha0_c1+CMS_alpha0_c2*temperature_average_i) - leakage_current_alpha.at(i).at(2) * log(time_history.back().at(1)) ;
+ }
+
+ i++;
+ }
+ }
+
+ else
+ {
+ tmp.at(0)=phi*t*leconsts.alpha_1;
+ tmp.at(1)=phi*t*leconsts.alpha_0_star; //temperature dependence of alpha0 is in the theta function of the time calculation - high temp = longer times and vice versa
+ tmp.at(2)=phi*t*leconsts.beta;
+
+ leakage_current_alpha.push_back(tmp);
+
+ tmp2.at(0)=t*leconsts.k01*exp(-leconsts.E1/(8.6173303e-5*(double)Temperature)); // time comes in seconds and k01 is seconds-1 so units are fine here
+ tmp2.at(1)=t / 60. * exp(-leconsts.E1_star*(1.0/(double)Temperature-1.0/294.15)/(8.6173303e-5)); //T-ref is hardcoded and fixed, it cannot be changed in a trivial way since the other parameters, especially alpha 0 star were evaluated at this reference temperature!!! // t/60 to convert time from seconds to minutes since this is the proposed unit here
+
+ if(time_history.size() > 0.1) tmp3=time_history.at(time_history.size()-1);
+
+ tmp2.at(0)+=tmp3.at(0);
+ tmp2.at(1)+=tmp3.at(1);
+
+ time_history.push_back(tmp2);
+
+ //pos1
+ int i=0;
+
+ while(i<leakage_current_alpha.size())
+ {
+ if(i>=1) G_i_tmp += leakage_current_alpha.at(i).at(0) * exp(-(time_history.back().at(0)-time_history.at(i-1).at(0))) + leakage_current_alpha.at(i).at(1) - leakage_current_alpha.at(i).at(2) * log(time_history.back().at(1)-time_history.at(i-1).at(1)) ;
+ else G_i_tmp += leakage_current_alpha.at(i).at(0) * exp(-time_history.back().at(0)) + leakage_current_alpha.at(i).at(1) - leakage_current_alpha.at(i).at(2) * log(time_history.back().at(1)) ;
+ i++;
+ }
+ }
+
+//////////////////////////////
+
+
+ G_i.push_back(G_i_tmp);
+ if(G_i_tmp/(totalDose*1.0e6+(double)phi+0.001)>5e-16)alpha_vec.push_back(1e-17); //totalDose was devided by 1e6 to fit into a long double
+ else alpha_vec.push_back(G_i_tmp/(totalDose*1.0e6+(double)phi+0.001));
+ leakage_current.push_back(G_i_tmp*1000.0*thickness*1.e-4); //insert volume here for current per module
+ powerconsumption.push_back(leakage_current.back()*NtoV_function(get_Neff()));
+
+ return;
+}
+
+
+/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+//function to read in the temp/rad profile
+
+vector<DataElement> get_Profile(string filename) // reads in a file of format (timestep/temperature/radiation_dose_rate) and gives back a vector of DataElements containing each the timestep/temperature/dose_rate information of one line of the input file
+{
+ ifstream input_file(filename.c_str()); // open stream to data file
+ vector<DataElement> Temperature_Profile_vector;
+ DataElement tmp_data;
+
+ if (input_file.is_open()) // check whether file opening was succesfull
+ {
+ if(debug) cout<< "Reading temperature/radiation file: " << filename << endl;
+
+ int timestep;
+ float temp1;
+ long int rad1;
+
+ while(true)
+ {
+ input_file >> timestep >> temp1 >> rad1;
+
+ //if(temp1<294 && temp1>290) //TODO: remove!
+ //{
+ // cout << "modifying temperature from: " << temp1 << endl;
+ // //temp1=temp1+5;
+ //}
+ //else //temp1=temp1+5;
+ //if (rad1!=0)temp1=temp1+5;
+
+ tmp_data.duration=timestep;
+ tmp_data.temperature=temp1;
+ tmp_data.dose_rate=rad1*DoseRateScaling;
+ Temperature_Profile_vector.push_back(tmp_data);
+
+ if(debug) cout << "Line read: " << timestep << " " << temp1 << " " << rad1 << endl;
+
+ if(input_file.eof())break; // if end of file is reached: brake the loop
+
+ limit_file_size_input--;
+ if(limit_file_size_input==0)break;
+ }
+ }
+
+ else
+ {
+ if(debug) cout << "Error opening the temperature/radiation file!"<<endl;
+ return Temperature_Profile_vector;
+ }
+
+ input_file.close();
+
+ return Temperature_Profile_vector;
+}
+
+/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+//function to plot different informations in root and pdf files
+
+void plot_vectors(vector<double> vecx, vector<double> vecy, string yName, string xName, long double totalDose, string plotname, int mode)
+{
+ string output_file_name_pdf=output_file_name;
+ output_file_name_pdf.insert(output_file_name_pdf.size(),"_");
+ output_file_name_pdf.insert(output_file_name_pdf.size(),plotname);
+ output_file_name_pdf.insert(output_file_name_pdf.size(),".pdf");
+
+ string output_file_name_root=output_file_name;
+ output_file_name_root.insert(output_file_name_root.size(),".root");
+
+
+ TFile *file = new TFile(output_file_name_root.c_str(),"UPDATE");
+ //file->cd();
+
+ if(mode==1 || mode==2 || mode==66 || mode==8) //in case of plot as a function of time, convert the days from simulation into a more handy date
+ {
+ for(int k=0; k<vecx.size();k++)
+ {
+ int day=vecx.at(k);
+
+ date d2 = d + days(day);
+
+ TDatime da (d2.year(), d2.month(), d2.day(),0,0,0);
+
+ vecx.at(k)=da.Convert();
+ }
+ }
+ else //otherwise the axis is fluence and we need to reconvert it since for safty reasons it was stored devided by 1e6 in order not to overflow the long double
+ {
+ //for(int k=0; k<vecx.size();k++)
+ //{
+ // vecx.at(k)*=1.0e6;
+ //}
+ }
+
+ if(mode==8 || mode==9) //in case of 8, change from unit area to unit surface for the leakage current (IMPLEMENT USER TREF HERE)
+ {
+ for(int k=0; k<vecy.size();k++)
+ {
+ vecy.at(k)/=thickness*1.0e-4; //convert from surface to volume normalisation
+ //vecy.at(k)*=(userTref*userTref/(294.15*294.15)*exp(-6500*(1/userTref - 1/294.15))); //scale to user defined temperature
+ vecy.at(k)*=(userTref*userTref/(294.15*294.15)*exp(-(bandGap/(2*8.6173303e-5))*(1/userTref - 1/294.15))); //scale to user defined temperature
+ }
+ }
+
+ const TVectorD t_timevector(vecx.size(),&vecx[0]);
+ const TVectorD t_vdepvector(vecy.size(),&vecy[0]);
+ const TVectorD t_time_error(vecx.size());
+ const TVectorD t_vdep_error(vecy.size());
+
+ TCanvas * c1 = new TCanvas("TCanvas_name","Title",0,0,1024,668);
+ gStyle->SetOptTitle(0);
+
+ TGraphErrors *gr = new TGraphErrors(t_timevector,t_vdepvector,t_time_error,t_vdep_error);
+
+ if(mode==1 || mode==2 || mode==66 || mode==8) //in case of plot as a function of time, convert the days from simulation into a more handy date
+ {
+ gr->GetXaxis()->SetTimeDisplay(1);
+ gr->GetXaxis()->SetNdivisions(6,2,0);
+ gr->GetXaxis()->SetTimeFormat("%d/%m/%Y");
+ gr->GetXaxis()->SetTimeOffset(0,"gmt");
+ }
+
+ gr->GetXaxis()->SetTitle(xName.c_str());
+ gr->GetYaxis()->SetTitle(yName.c_str());
+
+ if(mode!=66)
+ {
+ gr->SetName(plotname.c_str());
+ gr->Draw("AP");
+ if(plot_pdf) c1->Print(output_file_name_pdf.c_str());
+ gr->Write();
+ }
+
+ if(mode==1)
+ {
+ cout << "Final " << plotname << " is: " << vecy.at(vecy.size()-1)<< yName <<endl;
+ cout << "Total collected fluence is: " << totalDose*1.0e6/global_layer_conversion << "fb-1 (only one conversion factor supported - if more than one was used for the profile, this value will be wrong)" << endl;
+ cout << "Total collected dose is: " << totalDose*1.0e6 << "neq/cm2" << endl;
+ }
+ else
+ {
+ cout << "Final " << plotname << " is: " << vecy.at(vecy.size()-1)<<endl;
+ }
+
+ file->Close();
+ delete c1;
+}
+
+
+/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+int main(int argc, char *argv[]) {
+
+ if(argc!=2) // Check for correct number of arguments
+ {
+ cout << "Wrong number of arguments! There are: " << argc << " arguments." << endl;
+ return -1;
+ }
+
+ if(!overwrite_files)
+ {
+ time_t current_time = time(nullptr);
+ output_file_name.insert(output_file_name.size(),to_string(current_time));
+ }
+
+ string profilinput_file_name = argv[1];
+
+ leakage_current_consts LeCo = { 1.23e-17, //alpha_1
+ 7.07e-17, //alpha_0_star
+ 3.3e-18, //beta
+ 1.2e13, //k01
+ 1.11, //E1
+ 1.3}; //E1_star
+
+ Annealing_constants AnCo( 1.0e-2, //atof(argv[2]), //gA
+ 1.6e-2, //atof(argv[3]), //gY
+ 1.0e-2, //atof(argv[4]), //gC
+ 2.4e13, //ka
+ 7.4e14, //ky
+ 1.09, //Ea
+ 1.325, //Ey
+ 6.4118e-14); //cc
+
+ //DoseRateScaling = atof(argv[5]) ;
+
+
+ vector<DataElement> temp_rad_profile = get_Profile(profilinput_file_name); //read in the input data profile, a vector containing the individual data elements (duration, temperature, dose rate) is returned
+
+ int max_steps=temp_rad_profile.size();
+
+ vector<double> V_dep_vector;
+ vector<double> Neff_vector;
+ vector<double> Ndonor_vector;
+ vector<double> Nacceptor_vector;
+ vector<double> Temperature_vector;
+
+ vector<double> time_vector;
+ vector<double> fluence_vector;
+
+ Sensor *sensor = new Sensor();
+
+ long double time=0.0;
+ long double totalDose=0;
+
+ // main loop where all irradiation steps happen
+
+ cout << "Profile succesfully read. Length of the profile is: " << max_steps << endl;
+
+ for (int t=0; t<max_steps; t++) // iterate through the profile and irradiate the sensor at each step
+ {
+ sensor->set_temp(temp_rad_profile.at(t).temperature);
+ sensor->irradiate(LeCo,AnCo,temp_rad_profile.at(t).dose_rate,temp_rad_profile.at(t).duration,totalDose);
+
+ Temperature_vector.push_back(temp_rad_profile.at(t).temperature);
+
+ Neff_vector.push_back(sensor->get_Neff());
+ V_dep_vector.push_back(NtoV_function(sensor->get_Neff()));
+ Ndonor_vector.push_back(sensor->get_Ndonor());
+ Nacceptor_vector.push_back(sensor->get_Nacceptor());
+
+ totalDose+=temp_rad_profile.at(t).dose_rate/1.0e3*temp_rad_profile.at(t).duration/1.0e3; // time (seconds) * dose rate (neq/cm2/s)
+
+ time+=temp_rad_profile.at(t).duration;
+ time_vector.push_back(time/(24.0*3600.0)); // time vector in days
+ fluence_vector.push_back(totalDose*1.0e6);
+
+ if(t%(int)(max_steps/100.)==0) cout << (int)((int)t*100./max_steps) << " percent done..." << endl;
+ }
+
+ // data output, plotting and visualisation
+
+ cout << "Processing finished, writing data..." << endl;
+
+ // plots as function of time
+ plot_vectors(time_vector,Neff_vector,"N_{eff} [1/cm^{3}]","Date [days]",totalDose,"Neff",1);
+ plot_vectors(time_vector,Ndonor_vector,"N_{donor} [1/cm^{3}]","Date [days]",totalDose,"Ndonors",2);
+ plot_vectors(time_vector,Nacceptor_vector,"N_{acceptor} [1/cm^{3}]","Date [days]",totalDose,"Nacceptors",2);
+ plot_vectors(time_vector,V_dep_vector,"U_{dep} [V]","Date [days]",totalDose,"U_dep",66);
+ plot_vectors(time_vector,sensor->get_alpha_vec(),"#alpha [A/cm]","Date [days]",totalDose,"alpha",2);
+ plot_vectors(time_vector,sensor->get_leakage_current(),"I_{leak} (@21C) [mA/cm^{2}]","Date [days]",totalDose,"I_leak",2);
+ plot_vectors(time_vector,sensor->get_leakage_current(),"I_{leak} (@userTref) [mA/cm^{3}]","Date [days]",totalDose,"I_leak_volume",8);
+ plot_vectors(time_vector,sensor->get_G_i(),"G_{i} [A/cm^{3}]","Date [days]",totalDose,"G_i",2);
+ plot_vectors(time_vector,sensor->get_powerconsumption(),"P [mW/cm^{2}]","Date [days]",totalDose,"powerconsumption",2);
+ plot_vectors(time_vector,Temperature_vector,"Temperature [K]","Date [days]",totalDose,"temperature",2);
+ plot_vectors(time_vector,fluence_vector,"Fluence [n_{eq}/cm^{2}]","Date [days]",totalDose,"fluence",2);
+
+ // plots as function of dose
+ plot_vectors(fluence_vector,Neff_vector,"N_{eff} [1/cm^{3}]","Fluence [n_{eq}/cm^{2}]",totalDose,"Neff_vs_fluence",3);
+ plot_vectors(fluence_vector,Ndonor_vector,"N_{donor} [1/cm^{3}]","Fluence [n_{eq}/cm^{2}]",totalDose,"Ndonors_vs_fluence",4);
+ plot_vectors(fluence_vector,Nacceptor_vector,"N_{acceptor} [1/cm^{3}]","Fluence [n_{eq}/cm^{2}]",totalDose,"Nacceptors_vs_fluence",4);
+ plot_vectors(fluence_vector,V_dep_vector,"U_{dep} [V]","Fluence [n_{eq}/cm^{2}]",totalDose,"U_dep_vs_fluence",4);
+ plot_vectors(fluence_vector,sensor->get_alpha_vec(),"#alpha [A/cm]","Fluence [n_{eq}/cm^{2}]",totalDose,"alpha_vs_fluence",4);
+ plot_vectors(fluence_vector,sensor->get_leakage_current(),"I_{leak} [mA/cm^{2}]","Fluence [n_{eq}/cm^{2}]",totalDose,"I_leak_vs_fluence",4);
+ plot_vectors(fluence_vector,sensor->get_leakage_current(),"I_{leak} (@userTref) [mA/cm^{3}]","Fluence [n_{eq}/cm^{2}]",totalDose,"I_leak_volume_vs_fluence",9);
+ plot_vectors(fluence_vector,sensor->get_G_i(),"G_{i} [A/cm^{3}]","Fluence [n_{eq}/cm^{2}]",totalDose,"G_i_vs_fluence",4);
+ plot_vectors(fluence_vector,sensor->get_powerconsumption(),"P [mW/cm^{2}]","Fluence [n_{eq}/cm^{2}]",totalDose,"powerconsumption_vs_fluence",4);
+ plot_vectors(fluence_vector,Temperature_vector,"Temperature [K]","Fluence [n_{eq}/cm^{2}",totalDose,"temperature_vs_fluence",4);
+
+ return 0;
+}
diff --git a/InnerDetector/InDetCalibAlgs/PixelCalibAlgs/RadDamage/HamburgModel/output b/InnerDetector/InDetCalibAlgs/PixelCalibAlgs/RadDamage/HamburgModel/output
new file mode 100755
index 0000000000000000000000000000000000000000..748d625712c48f0ff9381005ac7784701d7010d7
Binary files /dev/null and b/InnerDetector/InDetCalibAlgs/PixelCalibAlgs/RadDamage/HamburgModel/output differ
diff --git a/InnerDetector/InDetCalibAlgs/PixelCalibAlgs/RadDamage/HamburgModel/temp_rad_data/converter b/InnerDetector/InDetCalibAlgs/PixelCalibAlgs/RadDamage/HamburgModel/temp_rad_data/converter
new file mode 100755
index 0000000000000000000000000000000000000000..8ff05b1700dc24f730ab6ebb9564f47dad1185a1
Binary files /dev/null and b/InnerDetector/InDetCalibAlgs/PixelCalibAlgs/RadDamage/HamburgModel/temp_rad_data/converter differ
diff --git a/InnerDetector/InDetCalibAlgs/PixelCalibAlgs/RadDamage/HamburgModel/temp_rad_data/profile_combiner.cpp b/InnerDetector/InDetCalibAlgs/PixelCalibAlgs/RadDamage/HamburgModel/temp_rad_data/profile_combiner.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..411eacd59027a1f465171751c0c35db125f5ea72
--- /dev/null
+++ b/InnerDetector/InDetCalibAlgs/PixelCalibAlgs/RadDamage/HamburgModel/temp_rad_data/profile_combiner.cpp
@@ -0,0 +1,380 @@
+/*
+Compile Code with: g++ profile_combiner.cpp -Wall -std=c++11 -o converter
+
+Program reads two files (see pattern below) in and orders them and prints out a file with at least 'time_spacing' hours spacing between the data points in the format time temperature(in K) delta_radiation(in neq/cm2/sec)
+
+temperature file:
+21-11-2011 15:45:40:838 +01.99
+21-11-2011 19:45:40:838 +02.99
+
+radiation file:
+2011-11-21 15:45:40.538 00001.234
+2011-11-21 22:48:52.330 00010.142
+
+with luminosity in pb-1
+
+both files have to start at a date later than 1.1.2011 01:01:01.001 and should start at the same date!!!
+
+what needs to be changed before starting:
+ 1. change conversion factor (from pb-1 to neq/cm2) before usage (two conversion factors are supported for 7/8 and 13 TeV)
+ 2. change output name
+ 3. call program after compiling with ./converter temperature_file_name irradiation_file_name
+*/
+
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+#include <iostream> // Basic IO
+#include <fstream> // Read and write files
+#include <cstdlib> // convert arguments of function call into int, float,...
+#include <string>
+#include <sstream> // to get string into stream
+#include <vector>
+#include <cmath>
+
+using namespace std;
+
+///////////////////////////////////////////////////////
+
+string define_output_file_name = "profile_output.txt";
+long int conversion_factor1 = 0.582e9; //conversion factor from pb-1 to neq/cm2 for 7TeV
+long int conversion_factor2 = 0.839e9; //for 13TeV
+long int changeinenergy = 2014; //year< since which there is 13TeV instead of 7TeV
+int time_spacing = 4; //minimal required spacing between two data points in hours
+int replacement88 = 2200; //replace no temperature reading from DCS (eg during shutdown) (=-88 C) with this temperature
+
+bool debug = false; //switch debugging information on/off
+
+///////////////////////////////////////////////////////
+
+struct newData
+{
+ unsigned long long int timestamp;
+ int temperature;
+ int luminosity;
+};
+
+///////////////////////////////////////////////////////
+
+struct data
+{
+ int year;
+ int month;
+ int day;
+ int hour;
+ int minute;
+ int second;
+ int milisecond;
+ unsigned long long int timestamp;
+ int payload;
+};
+
+///////////////////////////////////////////////////////
+
+unsigned long long int diffintime(int time1, int time2, int offset, bool uberold, unsigned long long int conversion)
+{
+ long long int diff=0;
+ bool uber = false;
+
+ diff = time1 - time2;
+
+ if(diff<0)
+ {
+ diff+=offset;
+ uber = true;
+ }
+
+ unsigned long long diff_long=abs(diff);
+
+ if(uberold)diff_long-=1;
+
+ diff_long*=conversion;
+
+ return diff_long;
+}
+
+///////////////////////////////////////////////////////
+
+unsigned long long int BuildTimeStamp(data newpoint, data oldpoint)
+{
+ bool next_uber=false;
+ long long int difference_milisecond = diffintime(newpoint.milisecond, oldpoint.milisecond, 1000, false, 1);
+ if(newpoint.milisecond-oldpoint.milisecond<0) next_uber = true;
+
+ unsigned long long int difference_second = diffintime(newpoint.second, oldpoint.second, 60, next_uber, 1000);
+ next_uber=false;
+ if(newpoint.second-oldpoint.second<0) next_uber = true;
+
+ unsigned long long int difference_minute = diffintime(newpoint.minute, oldpoint.minute, 60, next_uber, 60*1000);
+ next_uber=false;
+ if(newpoint.minute-oldpoint.minute<0) next_uber = true;
+
+ unsigned long long int difference_hour = diffintime(newpoint.hour, oldpoint.hour, 24, next_uber, 60*60*1000);
+ next_uber=false;
+ if(newpoint.hour-oldpoint.hour<0) next_uber = true;
+
+ unsigned long long int dayoffset=0;
+ if(oldpoint.month == 0) dayoffset=0;
+ else if(oldpoint.month == 2 )
+ {
+ dayoffset=28;
+ if(oldpoint.year==2008 || oldpoint.year==2012 || oldpoint.year==2016 || oldpoint.year==2020) dayoffset=29;
+ }
+ else if(oldpoint.month == 4 || oldpoint.month == 6 || oldpoint.month == 9 || oldpoint.month == 11)dayoffset=30;
+ else if(oldpoint.month == 1 || oldpoint.month == 3 || oldpoint.month == 5 || oldpoint.month == 7 || oldpoint.month == 8 || oldpoint.month == 10 || oldpoint.month == 12)dayoffset=31;
+ else
+ {
+ cout<<"Strange month, aborting..."<<oldpoint.month<<endl;
+ return -1;
+ }
+
+
+ unsigned long long int difference_day = diffintime(newpoint.day, oldpoint.day, dayoffset, next_uber, 24*60*60*1000);
+ //cout<<"difference in days is: "<<difference_day<< " "<<next_uber<<endl;
+ next_uber=false;
+ if(newpoint.day-oldpoint.day<0) next_uber = true;
+ //cout<<difference_day<< "ddd"<<endl;
+
+ unsigned long long int difference_month = diffintime(newpoint.month, oldpoint.month, 12, next_uber, dayoffset*24*60*60*1000);
+ next_uber=false;
+ if(newpoint.month-oldpoint.month<0) next_uber = true;
+
+ unsigned long long int difference_year = diffintime(newpoint.year, oldpoint.year, 100000, next_uber, 365.25*24*60*60*1000);
+
+
+ unsigned long long int newtimestamp = oldpoint.timestamp + difference_milisecond + difference_second + difference_minute + difference_hour + difference_day + difference_month + difference_year;
+
+ if(debug)cout<<"Timestamp calculation: "<<oldpoint.timestamp<<" "<<newtimestamp<<" "<<difference_milisecond<<" "<<difference_second<<" "<<difference_minute<<" "<<difference_hour<<" "<<difference_day<<" "<<difference_month<<" "<<difference_year<<endl;
+
+ return newtimestamp;
+}
+
+///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+int main(int argc, char *argv[]) {
+
+ string temp_file_name = argv[1];
+ string irr_file_name = argv[2];
+
+ ofstream output_file;
+ output_file.open(define_output_file_name.c_str());
+ ifstream temp_file(temp_file_name.c_str());
+ ifstream irr_file(irr_file_name.c_str());
+
+
+ if (!temp_file.is_open())
+ {
+ cout<<"Temp file not found"<<endl;
+ return -1;
+ }
+ if (!irr_file.is_open())
+ {
+ cout<<"Irr file not found"<<endl;
+ return -1;
+ }
+
+
+ cout<<"Temperature File: " << temp_file_name << endl;
+ cout<<"Irradiation File: " << irr_file_name << endl;
+
+
+ char column_dump_temp[100]; // max 100 charachters per line!!!
+ char column_dump_irr[100]; // max 100 charachters per line!!!
+
+ vector<data> temp_data_vec;
+ vector<data> irr_data_vec;
+
+ data temp_data_element;
+ data irr_data_element;
+
+ int year_irr, year_temp;
+ int month_irr, month_temp;
+ int day_irr, day_temp;
+ int hour_irr, hour_temp;
+ int minute_irr, minute_temp;
+ int second_irr, second_temp;
+ int miliseconds_irr, miliseconds_temp;
+ unsigned long long int timestamp_irr=0;
+ unsigned long long int timestamp_temp=0;
+ int irr_i, temp_i;
+
+ irr_data_element={2011,1,1,1,1,1,1,0,0};
+ temp_data_element={2011,1,1,1,1,1,1,0,0};
+
+ irr_data_vec.push_back(irr_data_element);
+ temp_data_vec.push_back(temp_data_element);
+
+
+//read in the two files (the two while loops) and put them into a vector each, calculate a unique timestamp for each point (in ms)
+
+ cout<<"Start to read in irr file..."<<endl;
+
+ int counter=0;
+
+ while(true)
+ {
+ irr_file.getline(column_dump_irr,100);
+
+ if(debug) cout << "Originally read line is: " << column_dump_irr << endl;
+
+
+ year_irr = (static_cast<int>(column_dump_irr[0])-48)*1000+(static_cast<int>(column_dump_irr[1])-48)*100+(static_cast<int>(column_dump_irr[2])-48)*10+static_cast<int>(column_dump_irr[3])-48;
+ month_irr = (static_cast<int>(column_dump_irr[5])-48)*10+static_cast<int>(column_dump_irr[6])-48;
+ day_irr = (static_cast<int>(column_dump_irr[8])-48)*10+static_cast<int>(column_dump_irr[9])-48;
+ hour_irr = (static_cast<int>(column_dump_irr[11])-48)*10+static_cast<int>(column_dump_irr[12])-48;
+ minute_irr = (static_cast<int>(column_dump_irr[14])-48)*10+static_cast<int>(column_dump_irr[15])-48;
+ second_irr = (static_cast<int>(column_dump_irr[17])-48)*10+static_cast<int>(column_dump_irr[18])-48;
+ miliseconds_irr = (static_cast<int>(column_dump_irr[20])-48)*100+(static_cast<int>(column_dump_irr[21])-48)*10+(static_cast<int>(column_dump_irr[22])-48);
+ irr_i = (static_cast<int>(column_dump_irr[24])-48)*10000000+(static_cast<int>(column_dump_irr[25])-48)*1000000+(static_cast<int>(column_dump_irr[26])-48)*100000+(static_cast<int>(column_dump_irr[27])-48)*10000+(static_cast<int>(column_dump_irr[28])-48)*1000+(static_cast<int>(column_dump_irr[30])-48)*100+(static_cast<int>(column_dump_irr[31])-48)*10+(static_cast<int>(column_dump_irr[32])-48);
+
+ if(debug) cout<<"Datapoint irr read: "<<irr_i<<" "<<miliseconds_irr<<" "<<second_irr<<" "<<minute_irr<<" "<<hour_irr<<" "<<day_irr<<" "<<month_irr<<" "<<year_irr<<endl;
+
+ if(year_irr>changeinenergy)irr_i=(int)round((double)irr_i*(double)conversion_factor2/(double)conversion_factor1);
+
+ irr_data_element={year_irr,month_irr,day_irr,hour_irr,minute_irr,second_irr,miliseconds_irr,timestamp_irr,irr_i};
+ timestamp_irr = BuildTimeStamp(irr_data_element, irr_data_vec.back());
+ irr_data_element={year_irr,month_irr,day_irr,hour_irr,minute_irr,second_irr,miliseconds_irr,timestamp_irr,irr_i};
+
+ if(timestamp_irr == -1) return -1;
+
+ irr_data_vec.push_back(irr_data_element);
+
+ if(irr_file.eof()) break;
+
+ counter++;
+ if(counter%200 == 0) cout << counter << " lines read. Year: " << irr_data_element.year << " Month:" << irr_data_element.month << " Day:" << irr_data_element.day << endl;
+ if(debug && (counter == 1 || counter ==2))system("read");
+ }
+
+ cout<<"Irr file completed! Starting with temperature file..."<<endl;
+
+ counter=0;
+ system("read");
+
+ int numberofreplacements=0;
+
+ while(true)
+ {
+ temp_file.getline(column_dump_temp,100);
+
+ if(debug) cout << "Originally read line is: " << column_dump_temp << endl;
+
+
+ year_temp = (static_cast<int>(column_dump_temp[6])-48)*1000+(static_cast<int>(column_dump_temp[7])-48)*100+(static_cast<int>(column_dump_temp[8])-48)*10+static_cast<int>(column_dump_temp[9])-48;
+ month_temp = (static_cast<int>(column_dump_temp[3])-48)*10+static_cast<int>(column_dump_temp[4])-48;
+ day_temp = (static_cast<int>(column_dump_temp[0])-48)*10+static_cast<int>(column_dump_temp[1])-48;
+ hour_temp = (static_cast<int>(column_dump_temp[11])-48)*10+static_cast<int>(column_dump_temp[12])-48;
+ minute_temp = (static_cast<int>(column_dump_temp[14])-48)*10+static_cast<int>(column_dump_temp[15])-48;
+ second_temp = (static_cast<int>(column_dump_temp[17])-48)*10+static_cast<int>(column_dump_temp[18])-48;
+ miliseconds_temp = (static_cast<int>(column_dump_temp[20])-48)*100+(static_cast<int>(column_dump_temp[21])-48)*10+(static_cast<int>(column_dump_irr[22])-48);
+ temp_i =(static_cast<int>(column_dump_temp[25])-48)*1000+(static_cast<int>(column_dump_temp[26])-48)*100+(static_cast<int>(column_dump_temp[28])-48)*10+(static_cast<int>(column_dump_temp[29])-48);
+
+ if(column_dump_temp[24]=='-')
+ {
+ temp_i*=-1;
+ if(debug)cout<<"negative temp"<<temp_i<<endl;
+ }
+
+ if(temp_i+8800<1)
+ {
+ numberofreplacements++;
+ temp_i=replacement88;
+ if(debug)cout << "temperature replaced! " << temp_i<< " "<< numberofreplacements<<"/"<<counter<<endl;
+ }
+ if(debug) cout<<"Datapoint temp read: "<<temp_i<<" "<<miliseconds_temp<<" "<<second_temp<<" "<<minute_temp<<" "<<hour_temp<<" "<<day_temp<<" "<<month_temp<<" "<<year_temp<<endl;
+
+ temp_data_element={year_temp,month_temp,day_temp,hour_temp,minute_temp,second_temp,miliseconds_temp,timestamp_temp,temp_i};
+ timestamp_temp = BuildTimeStamp(temp_data_element, temp_data_vec.back());
+ temp_data_element={year_temp,month_temp,day_temp,hour_temp,minute_temp,second_temp,miliseconds_temp,timestamp_temp,temp_i};
+
+ if(timestamp_temp == -1) return -1;
+
+ temp_data_vec.push_back(temp_data_element);
+
+ if(temp_file.eof()) break;
+ counter++;
+ if(counter%1000 == 0) cout << counter << " lines read. Year: " << temp_data_element.year << " Month:" << temp_data_element.month << " Day:" << temp_data_element.day << endl;
+ }
+
+ cout<<"Both files succesfully read! Beginning with sorting the data..."<<endl;
+
+//create a new vector with the ordered values in the format: timestamp Temperature luminosity
+
+ irr_data_vec.erase(irr_data_vec.begin()); // erase first element each because this was only a space holder
+ temp_data_vec.erase(temp_data_vec.begin());
+
+ int temp_iterator = 1;
+ int irr_iterator = 1;
+
+ newData element;
+ element={0,0,0};
+
+ vector<newData> output_vec;
+
+ while(true)
+ {
+ if(irr_data_vec.at(irr_iterator).timestamp>temp_data_vec.at(temp_iterator).timestamp)
+ {
+ element = {temp_data_vec.at(temp_iterator).timestamp, temp_data_vec.at(temp_iterator).payload, irr_data_vec.at(irr_iterator-1).payload};
+ temp_iterator++;
+ }
+
+ else
+ {
+ element = {irr_data_vec.at(irr_iterator).timestamp, temp_data_vec.at(temp_iterator-1).payload, irr_data_vec.at(irr_iterator).payload};
+ irr_iterator++;
+ }
+
+ output_vec.push_back(element);
+
+ if(temp_data_vec.size()<=temp_iterator || irr_data_vec.size()<=irr_iterator) break;
+ }
+
+
+ cout<<"Data is sorted now! Writing file with the data with fixed time spacing..."<<endl;
+
+
+ if(debug)
+ {
+ for(int k=0;k<output_vec.size();k++)
+ {
+ cout<<output_vec.at(k).timestamp<<" "<<output_vec.at(k).temperature<<" "<<output_vec.at(k).luminosity<<endl;
+ }
+ }
+
+
+//write out a file with a constant spacing of time_spacing hours, temperature and the converted irradiation is rounded to int
+
+ output_file<<1*60*60<<" "<<(int)round(273.0+(output_vec.at(0).temperature/100.0))<<" "<<(long int)round((conversion_factor1*(output_vec.at(0).luminosity)/1000.0/60.0/60.0))<<endl;
+
+ int old_i = 0;
+ int i = 1;
+ int time_difference=0;
+ int radiation_difference=0;
+
+ while(i<output_vec.size()-2)
+ {
+ while(i<output_vec.size()-1 && output_vec.at(i).timestamp-output_vec.at(old_i).timestamp<time_spacing*60*60*1000) //wait until difference between to points is at least time_spacing hour
+ {
+ i++;
+
+ if(debug) cout<<"Position in output file:"<<i << " and file size is: " << output_vec.size() <<endl;
+ }
+
+ time_difference=(int)round((output_vec.at(i).timestamp-output_vec.at(old_i).timestamp)/1000.0);
+ radiation_difference=(long int)round((conversion_factor1*(output_vec.at(i).luminosity-output_vec.at(old_i).luminosity)/1000.0/(float)time_difference));
+
+ if(output_vec.at(i).luminosity-(int)round((double)conversion_factor2/(double)conversion_factor1*output_vec.at(old_i).luminosity)==0)
+ {
+ radiation_difference=0;
+ cout << "corrected for change in conversion factor, should not happen more than once!"<<endl;
+ }
+ output_file<<time_difference<<" "<<(int)round(273.0+(output_vec.at(i).temperature/100.0))<<" "<<radiation_difference<<endl;
+
+ old_i=i;
+ i++;
+ }
+
+ output_file.close();
+
+ return 0;
+}
diff --git a/InnerDetector/InDetCalibAlgs/PixelCalibAlgs/RadDamage/HamburgModel/temp_rad_data/temp_IBL_2015-jan2018.txt b/InnerDetector/InDetCalibAlgs/PixelCalibAlgs/RadDamage/HamburgModel/temp_rad_data/temp_IBL_2015-jan2018.txt
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391