{
#include <fstream.h>
gROOT->LoadMacro("AtlasUtils.C");
ATLAS->SetPadTickX(1);
ATLAS->SetPadTickY(0);
ATLAS->SetPadRightMargin(0.12);
ATLAS->SetPadLeftMargin(0.12);
c1 = new TCanvas("c1","A Simple Graph",600,10,700,500);
//c1->SetFillStyle(4100);
//TDatime db(2010,01,01,00,00,00);
TDatime db(2010,01,01,00,00,00);
int t1 = db.Convert();
//TDatime db(2011,04,01,00,00,00);
TDatime db(2021,12,31,00,00,00);
int t2 = db.Convert();
gStyle->SetOptStat(0);
//c1->SetGrid();
// create the arrays for the points
Double_t ymax= 200;
double Tscale=ymax*0.1/20., Toffset = ymax*0.8;
//double Tscale=0.01/10., Toffset = 0.02;
Int_t n = 3000;
double intL[n],Temp[4][n],Temp_plot[4][n], Tbias=3.6;
double time[n], Ttime[n];
double weekL[n], energy[n];
double leak[4][n],dleak[4][n],leakop[4][n],dleakop[4][n],
leak3[5][n],dleak3[5][n];
int week,year1,month1,day1,year2,month2,day2; 
//--------readin luminosity data-------------------
int iweek=1, iT=0;
ifstream fin("lhcall_7_temperature_plan1.txt");
while(fin>>year1>>month1>>day1>>year2>>month2>>day2>>energy[iweek]>>
weekL[iweek]>>intL[iweek]>>Temp[0][iweek]>>Temp[1][iweek]>>Temp[2][iweek]>>Temp[3][iweek] )
{
cout<<year1<<" "<<month1<<" "<<day1<<" "<<year2<<" "<<month2<<" "<<day2<<" "<<energy[iweek]<<" "<<
weekL[iweek]<<" "<<intL[iweek]<<" "<<Temp[0][iweek]<<" "<<Temp[1][iweek]<<" "<<Temp[2][iweek]<<" "<<Temp[3][iweek] <<endl;
  intL[iweek]= intL[iweek]/1000000.; 
  TDatime d1(year1,month1,day1,12,0,0);
  int time1=d1.Convert();
  TDatime d2(year2,month2,day2,12,0,0);
  time[iweek]=d2.Convert();
  Ttime[iT]=time1;
  for(int layer=0; layer<4; layer++){
    Temp[layer][iweek] = Temp[layer][iweek] - Tbias;
    Temp_plot[layer][iT] = Temp[layer][iweek]*Tscale+Toffset;
  }
  iT++;
  Ttime[iT]=time[iweek];
  for(int layer=0; layer<4; layer++) 
    Temp_plot[layer][iT] = Temp[layer][iweek]*Tscale+Toffset;
   cout<<"iweek, time="<<iweek<<" "<<time[iweek]<<endl;
  iT++;
  iweek++;
}
fin.close();
cout<<"number of time intervals= "<< iweek << endl;
cout<<"last week: from "<<year1<<" "<<month1<<" "<<day1<<" to "<<
year2<<" "<<month2<<" "<<day2<< endl;

//-------------------Calculate prediction -----------------
//

double c = 1.04e-13;
double gc = 0.017;
double ga = 0.018;
double taua = 2.290;
double Ea = 1.090;
double gy = 0.059;
double tauy = 480.;
double Ey = 1.330;
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;
double d = 0.0285, S = 150.7983;
double Vd0 = 70., dVd0 = 20.;
double Neff00 = Vd0 * 2.* epsilon_prime*epsilon0/e/pow(d,2); 
double dNeff0 = dVd0* 2.* epsilon_prime*epsilon0/e/pow(d,2);
double NC00 = 0.7 * Neff00;
cout<<"Neff00="<<Neff00<<endl;
double RHV = 11200., Vmax = 430.; 
double Fluence_7TeV[4] = {1.65e5, 1.29e5, 1.07e5, 9.0e4};
double Fluence_14TeV[4];
       Fluence_14TeV[0] = Fluence_7TeV[0] * 1.153;
       Fluence_14TeV[1] = Fluence_7TeV[1] * 1.14;
       Fluence_14TeV[2] = Fluence_7TeV[2] * 1.13;
       Fluence_14TeV[3] = Fluence_7TeV[3] * 1.114;
double TR=293.15, T0=273.15; 
double Neff[4][15][n], Vd[4][15][n];
double dNy[4][15][n], C[4][15][n];
double fluence0, fluence, total_fluence0[4][n], total_fluence[4][n];
double Tsensor, Tsensor0;
//------------initialization------------------------------
time[0] = time[1] - 24.* 60.*60;
for (int layer=0; layer<4; layer++) {
  total_fluence0[layer][0]= 0.;
  total_fluence[layer][0]= 0.;
  for (int ierror=0; ierror<15; ierror++) {
     C[layer][ierror][0] = 0.; 
     dNy[layer][ierror][0] = 0.;
  } 
}
//---------------------go over all weeks-----------------------------
for (int iw=1; iw<iweek; iw++) {
  double days = (time[iw] - time[iw-1])/(24.*60.*60);
  for (layer=0; layer<4; layer++) {
    if(energy[iw] < 10.) fluence0 = weekL[iw]*Fluence_7TeV[layer];
       else fluence0 = weekL[iw]*Fluence_14TeV[layer];
    total_fluence0[layer][iw] = total_fluence0[layer][iw-1] + fluence0;
    double Tsensor0 = Temp[layer][iw];
    cout<<iw<<" T="<<Tsensor0<<" E,weekL,fluence="<<energy[iw]
    <<" "<<weekL[iw]<<" "<<fluence0<<" total_fluence="
    <<total_fluence0[layer][iw]<<endl;
//---------start loop for error estimate (0 to 14)------------------------
    double isum[15], differenceError[15];
    double quadSum = 0.;
    for (int ierror=0; ierror<3; ierror++) {
      fluence = fluence0;
      Tsensor = Tsensor0;
      double Neff0 = Neff00;
      double NC0 = NC00;
      total_fluence[layer][iw] = total_fluence0[layer][iw];  
      if( ierror==1 ) Neff0 += dNeff0;
      if( ierror==2 ) Neff0 -= dNeff0;
//      if( ierror==12) fluence *= 1.045;
//      if( ierror==13 && Tsensor0 < 10.) Tsensor += 1.0;
//      if( ierror==14 && Tsensor0 >=10.) Tsensor += 2.0;
//------calculate each component of leak current, Hamburg model -----------------
      double ThetaTtau = exp(Ea/kB*(1/TR-1/(T0+Tsensor))) * days / taua;
      double A = Neff0 - NC0*(1-exp(-c*total_fluence[layer][iw]));
      double B = gc * total_fluence[layer][iw];
      C[layer][ierror][iw] = ga * fluence/ThetaTtau*(1-exp(-ThetaTtau))
              + C[layer][ierror][iw-1]*exp(-ThetaTtau);
      double Theta_y = exp( Ey / kB * (1./TR - 1./(T0 + Tsensor)));
      double gyF=gy * (total_fluence[layer][iw-1]
                 + total_fluence[layer][iw])/2.;
      dNy[layer][ierror][iw] = dNy[layer][ierror][iw-1] 
         + gyF * (1-1/(1+Theta_y * days / tauy ));
      Neff[layer][ierror][iw] = A-B-C[layer][ierror][iw]
                                - dNy[layer][ierror][iw];
      Vd[layer][ierror][iw] = e * d * d * fabs(Neff[layer][ierror][iw])
                      /(2*epsilon0*epsilon_prime);
      if(iw==1)Vd[layer][ierror][0] = Vd[layer][ierror][1];
      cout<<"iw,layer,ierror="<<iw<<" "<<layer<<" "<<ierror
      <<" days="<<days<<" T="<<Tsensor<<" dNy=" 
      <<dNy[layer][ierror][iw] <<" Neff="<<Neff[layer][ierror][iw]
      <<" Vd="<<Vd[layer][ierror][iw]<<endl;   
    }
  }
}

//-----------------------------------------------------
//
leg0 = new TLegend(0.20 ,0.55 ,0.55,0.75);
for (int layer=0; layer<4; layer++){
  hc = new TGraph(iweek,time,Vd[layer][0]);
  hu = new TGraph(iweek,time,Vd[layer][1]);
  hd = new TGraph(iweek,time,Vd[layer][2]);
  hu->SetMaximum(ymax);

  hu->GetYaxis()->SetTitle("V_{d} [V]");
  //hu->GetXaxis()->SetTitle("DATE");
  hu->GetXaxis()->SetTitleOffset(1.4);
  hu->GetXaxis()->SetLabelOffset(0.03);
  hu->GetXaxis()->SetLabelSize(0.04);
  hu->GetYaxis()->SetTitleSize(0.05);
  hu->GetYaxis()->SetTitleOffset(1.0);
  hu->GetYaxis()->SetLabelOffset(0.01);
//  hu->GetXaxis()->SetRangeUser(13,45);
  hu->GetXaxis()->SetTimeDisplay(1);
  hu->GetXaxis()->SetTimeFormat("#splitline{  %Y}{%b %d}");
  hu->GetXaxis()->SetLimits(t1,t2);
//  hu->GetXaxis()->SetNdivisions(206);
  hu->GetXaxis()->SetNdivisions(206);
  hc->SetLineColor(layer+2);
  hu->SetLineColor(layer+2);
  hd->SetLineColor(layer+2);

  if(layer==3){
  hc->SetLineColor(1);
  hu->SetLineColor(1);
  hd->SetLineColor(1);}
  if(layer==0)hu->Draw("AL");
  if(layer!=0)hu->Draw("SAME L");
  hc->Draw("SAME L");
  hd->Draw("SAME L");
  if(layer==0)leg0->AddEntry(hu,"   B3","l");
  if(layer==1)leg0->AddEntry(hu,"   B4","l");
  if(layer==2)leg0->AddEntry(hu,"   B5","l");
  if(layer==3)leg0->AddEntry(hu,"   B6","l");
}

gT3 = new TGraph(iT,Ttime,Temp_plot[0]);
gT4 = new TGraph(iT,Ttime,Temp_plot[1]);
gT5 = new TGraph(iT,Ttime,Temp_plot[2]);
gT6 = new TGraph(iT,Ttime,Temp_plot[3]);
gT3->SetLineStyle(2); gT3->SetLineColor(2);gT3->Draw("SAME L");
gT4->SetLineStyle(2); gT4->SetLineColor(3);gT4->Draw("SAME L");
gT5->SetLineStyle(2); gT5->SetLineColor(4);gT5->Draw("SAME L");
gT6->SetLineStyle(2); gT6->SetLineColor(1);gT6->Draw("SAME L");
//
gL = new TGraph(iweek,time,intL);
gL->SetLineStyle(3); gL->SetLineColor(1);gL->Draw("SAME L");
leg0->AddEntry(gL,"int-Lum (fb^{-1})","l");
//
leg0->SetFillColor(10); leg0->SetBorderSize(0); leg0->SetTextSize(0.04);
leg0->Draw();
//draw an axis on the right side
TGaxis *axis = new TGaxis(t2,140.,t2, ymax,-20.,40.,505,"+L");
axis->SetTitle("T_{sensor} [#circC]");
axis->SetLabelOffset(0.01);
axis->SetTitleOffset(1.0);
axis->SetTitleSize(0.05);
axis->SetLabelSize(0.05);
axis->SetLabelFont(42);
axis->SetTitleFont(42);
axis->Draw();
//-- write last day
TDatime db(2021,09,01,00,00,00);int t4 = db.Convert();
legA = new TText(t4 ,-10.4 ,"2022");
legA->SetTextSize(0.040);legA->Draw();
TDatime db(2021,08,01,00,00,00);int t4 = db.Convert();
legB = new TText(t4 ,-20.5 ,"Jan 01");
legB->SetTextSize(0.040);legB->Draw();
//-- write start day
TDatime db(2009,09,01,00,00,00);int t4 = db.Convert();
legA = new TText(t4 ,-10.4 ,"2010");
legA->SetTextSize(0.040);legA->Draw();
TDatime db(2009,08,01,00,00,00);int t4 = db.Convert();
legB = new TText(t4 ,-20.5 ,"Jan 01");
legB->SetTextSize(0.040);legB->Draw();
//
TDatime db(2015,01,01,00,00,00);int t4 = db.Convert();
legA = new TLatex(t4 ,185 ,"sensor temperature");
legA->SetTextSize(0.040);legA->Draw();
//
TDatime db(2015,01,01,00,00,00);
int t4 = db.Convert();
legA = new TText(t4 ,140 ,"Cooling plan 1");
legA->SetTextSize(0.050);
legA->Draw();

c1->Print("Vd_plan1.png");
}