{
//===========================================================================
//  Sample ROOT programme for Vd evolution, Taka Kondo (KEK, ATLAS) 2011.10.24
//  2012.11.8 updated
//===========================================================================
//------------------set up for plot --------------------------------------
ATLAS->SetPadTickY(0);
ATLAS->SetPadRightMargin(0.12);
ATLAS->SetPadLeftMargin(0.12);
c1 = new TCanvas("c1","plot_Vd",500,30,700,500);
c1->SetGrid();
//--------1 MeV n-eq fluence [/cm^2] at 1 pb^-1 7 TeV pp collision------------//https://twiki.cern.ch/twiki/bin/viewauth/Atlas/BenchmarkingAtTheLHC
double NIEL0[8] = {21.5e8, 8.91e8, 5.57e8, 1.66e8, 1.30e8,1.07e8, 0.90e8, 21.5e8};
double d[8]     = {0.0250, 0.0250,0.0250,0.0285,0.0285,0.0285,0.0285, 0.0200}; 
double Vd0[8]   = {65.,65.,65.,70.,70.,70.,70.,65.}; // initial Vd values
double Factor7_14= 1.24,  Factor7_8=1.05;
//---------set start and end days for plot ---------------
TDatime TX(2010,01,01,00,00,00);
int X0 = TX.Convert();
gStyle->SetTimeOffset(X0); 
TDatime db(2010,01,01,00,00,00);
int t1 = db.Convert()-X0;
TDatime db(2022,01,01,00,00,00);
int t2 = db.Convert()-X0;
double yminPlot = 0., ymaxPlot= 800;                      // min max of plot
double Tscale=ymaxPlot*0.025/10., Toffset = ymaxPlot*0.9; // for temperature plot
//=============================== cooling plan ==============================
int cooling_plan=1 ;
//===========================================================================
int n = 1000;
double Lmax = 8;          // L=0(B0)......6(B6),7(B0 with d=200um)
double Temp[Lmax][n],Temp_plot[Lmax][n];
double time[n], Ttime[n], weekL[n]];
double Neff[n], dNy[n], C[n], total_fluence[n];
double Vdc[Lmax][n]; 
//--------readin days_luminosity_temperature.txt------------------
int    year1,month1,day1,hour1,min1,year2,month2,day2,hour2,min2;
double T1, T2, T3, T4, T5, T6;
int    iweek=1, iweekT=1;
char   line[255];
ifstream fin("intL_B0_B3.txt");
fin.getline(line,sizeof(line));   //------skip the 1st dummy line---------
while(fin>>year1>>month1>>day1>>hour1>>min1>>year2>>
      month2>>day2>>hour2>>min2>>weekL[iweek]>>T1>>T2>>T3>>T4>>T5>>T6) {
  if(cooling_plan==1) {Temp[0][iweek]=T1; Temp[3][iweek]=T2;}
  if(cooling_plan==2) {Temp[0][iweek]=T3; Temp[3][iweek]=T4;}
  if(cooling_plan==3) {Temp[0][iweek]=T5; Temp[3][iweek]=T6;}
//---- check if time is contiguos--------------
  TDatime d1(year1,month1,day1,hour1,min1,0.);
  double time1 =d1.Convert()-X0;
  if(iweek > 1 && time1 != time[iweek-1]) {
     cout<<"discontinuoe time found;"<<year1<<" "<<month1<<" "<<day1
       <<" "<<hour1<<" "<<min1<<endl;
     stop;
  }
  if(iweek==1) {
    time[0] = time1; 
    Ttime[0]= time1;
    Temp[0][0] = Temp[0][iweek]; Temp[1][0] = Temp[0][iweek]; 
    Temp[2][0] = Temp[0][iweek]; Temp[3][0] = Temp[3][iweek];
    Temp[4][0] = Temp[3][iweek]; Temp[5][0] = Temp[3][iweek];
    Temp[6][0] = Temp[3][iweek]; Temp[7][0] = Temp[0][iweek];
  }
  TDatime d2(year2,month2,day2,hour2,min2,0.);
  time[iweek]=d2.Convert()-X0;
  Temp[1][iweek] = Temp[0][iweek]; Temp[2][iweek] = Temp[0][iweek];
  Temp[4][iweek] = Temp[3][iweek]; Temp[5][iweek] = Temp[3][iweek];
  Temp[6][iweek] = Temp[3][iweek]; Temp[7][iweek] = Temp[0][iweek];
  if(Temp[3][iweek] < 10.)  Temp[6][iweek] = Temp[3][iweek]+8.0;
//-----adjust fluence due to colllison energy ----------
  if(year1==2012) weekL[iweek] *=  Factor7_8 ; 
  if(year1> 2012) weekL[iweek] *=  Factor7_14; 
//-----for temperature plot----------------------
  Ttime[iweekT] = time[iweek];
  for (int i=0; i<8; i++){
    if(iweekT==1) Temp_plot[i][0]=Temp[i][1]*Tscale+Toffset;
    Temp_plot[i][iweekT] = Temp[i][iweek-1]*Tscale+Toffset;
  }
  iweekT++;
  Ttime[iweekT] = time[iweek];
  for (int i=0; i<7; i++) Temp_plot[i][iweekT] = Temp[i][iweek]*Tscale+Toffset;
  iweekT++;
  iweek++;
}
fin.close();
cout<<"Finished reading input file, No of weeks = "<<iweek<<endl;
//
//---------------- Calculate full depletion voltage -----------------
//
double kB = 8.617343e-5;
double TR=293.15, T0=273.15; 
double epsilon0 = 8.854e-14, epsilon_prime=11.6;
double e = 1.6020e-19;
//------------initialization------------------------------
//
for (int L=0; L<Lmax; L++){
  if(L>2 && L<7) {
     //-------parameters and their errors  standard Si----------------------
     double C0coeff0 = 0.7;
     double ccoeff0 = 0.075; 
     double ga0 = 0.018;
     double taua0 = 2.290; 
     double Ea0 = 1.090;
     double gc0 = 0.017;
     double gy0 = 0.059;
     double tauy0 = 480.;  // just for trial
     double Ey0 = 1.330;
  }
  else {
     //-------parameters for Oxygenated Silicon------------------------
     double C0coeff0 = 0.45;
     double ccoeff0 = 0.075;
     double ga0 = 0.014;
     double taua0 = 2.917;
     double Ea0 = 1.090;
     double gc0 = 0.020;
     double gy0 = 0.048;
     double tauy0 = 800.;
     double Ey0 = 1.330;  
  }
  double Neff00 =  Vd0[L] * 2.* epsilon_prime*epsilon0/e/pow(d[L],2); 
  cout<<"Layer="<<L<<", initial Neff0="<<Neff00<<endl;
  Neff[0] = Neff00;
  dNy[0] = 0.;
  C[0] = 0.; 
  total_fluence[0] = 0.;
// 
//---------------------go over all weeks-----------------------------
//
  for (int iw=1; iw<iweek; iw++) {
      double day = float(time[iw] - time[iw-1])/(24.*60.*60.);
      double Neff0 = Neff00;
      double C0coeff = C0coeff0;
      double ccoeff = ccoeff0;
      double ga = ga0;
      double taua = taua0;
      double Ea = Ea0;
      double gc = gc0;
      double gy = gy0;
      double tauy = tauy0;
      double Ey = Ey0;
      double Tsensor = Temp[L][iw];
      double luminosity = weekL[iw];
      double NIEL = NIEL0[L]; 
//------calculate each component of the Hamburg model ---------------
      double NC0 = C0coeff * Neff0;
      total_fluence[iw] = total_fluence[iw-1] + luminosity * NIEL;  
      double ThetaTtau = exp(Ea/kB*(1/TR-1/(T0+Tsensor))) * day / taua;
      double c = ccoeff / NC0;
      double A = Neff0 - NC0*(1-exp(-c*total_fluence[iw])); //Donor removal
      double B = gc * total_fluence[iw];               // stable acceptors
      C[iw] = ga * luminosity * NIEL / ThetaTtau*(1-exp(-ThetaTtau))
               + C[iw-1]*exp(-ThetaTtau);   // unstable acceptors
      double Theta_y = exp( Ey / kB * (1./TR - 1./(T0 + Tsensor)));
      double gyF=gy * (total_fluence[iw-1]  + total_fluence[iw])/2.;
      dNy[iw] = dNy[iw-1] + gyF * (1-1/(1+Theta_y * day/tauy));// reverse annealing
      Neff[iw] = A - B - C[iw] - dNy[iw];
//---------------calculate full depletion voltage-------------------------  
      Vdc[L][iw] = e * d[L] * d[L] * fabs(Neff[iw])/(2*epsilon0*epsilon_prime);
//---------print out last values------------------------
      if (iw==iweek-1)   
         cout<<"iw="<<iw<<" time="<<time[iw]<<" L="<<L<<", days="<<day<<", T="
             <<Tsensor<<" Neff=" <<Neff[iw] <<" Vd="<<Vdc[L][iw]<< endl;
  }
}  // end of loop of L 
//
for (int L=0; L<Lmax; L++) Vdc[L][0] = Vdc[L][1]; 
//
//-----------------------------------------------------
//----------draw the 1st plot---------------
h1 = new TGraph(iweek,time,Vdc[0]);
h1->SetMaximum(ymaxPlot);
h1->SetMinimum(yminPlot);
h1->SetMarkerColor(2);   
h1->SetLineColor(1);   
h1->GetYaxis()->SetTitle("V_{d} [V]");
h1->GetXaxis()->SetLabelOffset(0.03);
h1->GetXaxis()->SetLabelSize(0.00);   
h1->GetYaxis()->SetLabelSize(0.04);
h1->GetXaxis()->SetTitleSize(0.04);
h1->GetYaxis()->SetTitleSize(0.04);
h1->GetXaxis()->SetTitleOffset(1.7);
h1->GetYaxis()->SetTitleOffset(1.2);
h1->GetYaxis()->SetLabelOffset(0.01);
h1->GetXaxis()->SetTimeDisplay(1);
h1->GetXaxis()->SetTimeFormat("#splitline{  %Y}{%b %d}");
h1->GetXaxis()->SetLimits(t1,t2);
h1->SetLineWidth(2);
h1->Draw("AL");  
//-----------------------------------------------------------
sprintf(line,"Cooling Plan %d\0",cooling_plan);
leg = new TLegend(0.15,0.50,0.43,0.80,line);
leg->AddEntry(h1,"pixel B0 Oxygenated Si","l");
for (int i=1; i<Lmax; i++) {
   h2 = new TGraph(iweek,time,Vdc[i]);
   h2->SetLineColor(i+1);
   if(i>=4) h2->SetLineColor(i+2);
   h2->SetLineWidth(2);
   if(i==7) {h2->SetLineColor(1); h2->SetLineWidth(1);}
   h2->Draw("SAME L"); 
   if(i<3) sprintf(line,"pixel B%d Oxygenated Si\0",i);
   else sprintf(line,"SCT B%d Standard Si\0",i);
   if(i==7)sprintf(line,"B0 Oxygen. 200#mum thick\0",i);
   leg->AddEntry(h2,line,"l");
  }
//-------------------------plot temperature--------------------------
for (int i=0; i<7; i+=3){
  gT = new TGraph(iweekT,Ttime,Temp_plot[i]);
  gT->SetLineColor(i+1);gT->SetLineWidth(2);gT->Draw("SAME L"); 
}
//-- write years----
for (int iy=2011; iy<2022; iy+=2){
  TDatime db(iy-1,08,01,00,00,00);int t4 = db.Convert()-X0;
  sprintf(line,"%d\0",iy);
  legB = new TText(t4 ,-40 ,line);legB->SetTextSize(0.040);legB->Draw();
  TDatime db(iy-1,06,01,00,00,00);int t4 = db.Convert()-X0;
  legC = new TText(t4 ,-70 ,"Jan 01");legC->SetTextSize(0.040);legC->Draw();
}
//-- LS1 LS2 box----
TDatime db(2013,02,11,00,00,00);int t4 = db.Convert()-X0;
TDatime db(2015,02,08,00,00,00);int t5 = db.Convert()-X0;
z = new TLine(t4,0.,t4,ymaxPlot); z->Draw();
z = new TLine(t5,0.,t5,ymaxPlot); z->Draw();
t=new TText(t4,ymaxPlot*0.95,"     LS1"); t->SetTextSize(0.04); t->Draw(); 
TDatime db(2017,12,11,00,00,00);int t4 = db.Convert()-X0;
TDatime db(2019,02,17,00,00,00);int t5 = db.Convert()-X0;
z = new TLine(t4,0.,t4,ymaxPlot); z->Draw();
z = new TLine(t5,0.,t5,ymaxPlot); z->Draw();
t=new TText(t4,ymaxPlot*0.95,"  LS2"); t->SetTextSize(0.04); t->Draw();
//-------------draw an axis on the right side for temperature--------
TGaxis *axis = new TGaxis(t2,0.,t2,ymaxPlot*0.8,0.,ymaxPlot*0.8,510,"+L");
axis->SetLabelSize(0.00); axis->Draw();
TGaxis *axis = new TGaxis(t2,ymaxPlot*0.8,t2, ymaxPlot,-40.,40.,202,"+L");
axis->SetTitle("T_{sensor} [#circC]");
axis->SetLabelOffset(0.01);
axis->SetTitleOffset(1.0);
axis->SetTitleSize(0.04);
axis->SetLabelSize(0.04);
axis->SetLabelFont(42);
axis->SetTitleFont(42);
axis->Draw();
z = new TLine(t2,0.,t2,ymaxPlot); z->Draw();
//---------write legend -----------------------
leg->SetFillColor(10); leg->SetBorderSize(2); leg->SetTextSize(0.03);
leg->Draw();  
//
sprintf(line,"Vd_plan-%d.png\0",cooling_plan);
c1->Print(line);
}