{
#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",500,10,700,500);
c1->SetGrid();
//TDatime db(2010,06,01,00,00,00);
TDatime db(2010,01,01,00,00,00);
int t1 = db.Convert();
//TDatime db(2021,12,31,00,00,00);
TDatime db(2021,12,31,00,00,00);
int t2 = db.Convert();
gStyle->SetOptStat(0);
gPad->SetLogy();
//c1->SetGrid();
// create the arrays for the points
Double_t ymax= 10000., ymin=0.1;
double Tscale=ymax*0.1/10., Toffset = ymax*0.2;
//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 Tplot0[4]={1.35,1.12,1.83,8.75};
//double Tplot0[4]={-7.4,-7.4,-7.4,-7.4};
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] = pow(10., 3+(Temp[layer][iweek]+10)/40.);
  }
  iT++;
  Ttime[iT]=time[iweek];
  for(int layer=0; layer<4; layer++) 
    Temp_plot[layer][iT] =  pow(10., 3+(Temp[layer][iweek]+10)/40.);
   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 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 A0[5] = {0.42, 0.10, 0.23, 0.21, 0.04}, A[5];
double dA[5]= {0.11, 0.01, 0.02, 0.02, 0.03};
double tau0[5] = {833.33, 28.47, 2.57, 0.086, 0.0056}, tau[5];
double dtau[5]= {138.89, 4.17, 0.21, 0.017, 0.0035};
double Eg = 1.12, EI=1.09, kB=8.617343E-5, TR=293.15, T0=273.15; 
double alpha0 = 6.90e-18, dalpha = 0.20e-18;
double d = 0.0285, width = 6.144, length = 6.136, nsensor = 4.0;
double Ssensor = width*length*nsensor;
double Vsensor = Ssensor * d;
double ileak[4][n][5][15], leak[4][n], dleak[4][n];
double fluence0, fluence;
double Tsensor, Tsensor0;
time[0] = time[1] - 24.* 60.*60;
for (int layer=0; layer<4; layer++) {
  for (int k=0; k<5; k++) {
    for (int ierror=0; ierror<15; ierror++) ileak[layer][0][k][ierror] = 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];
    double Tsensor0 = Temp[layer][iw];
//    cout<<"Temp="<<Temp[layer][iw]<<" energy,weekL,fluence="<<energy[iw]
//    <<" "<<weekL[iw]<<" "<<fluence<<endl;
//---------start loop for error estimate (0 to 14)------------------------
    double isum[15], differenceError[15];
    double quadSum = 0.;
    for (int ierror=0; ierror<15; ierror++) {
      double C;
      double alpha = alpha0;
      for(int k=0;k<5;k++) A[k] = A0[k];
      for(int k=0;k<5;k++) tau[k] = tau0[k];
      fluence = fluence0;
      Tsensor = Tsensor0; 
      if( ierror==1 ) alpha += dalpha;
      if( ierror==2 ){ C=dA[0]/(1.-A0[0]); A[0]+=  dA[0]; A[1]-=C*A0[1];
                       A[2]-=C*A0[2];      A[3]-=C*A0[3]; A[4]-=C*A0[4];}
      if( ierror==3 ){ C=dA[1]/(1.-A0[1]); A[0]-=C*A0[0]; A[1]+=  dA[1];
                       A[2]-=C*A0[2];      A[3]-=C*A0[3]; A[4]-=C*A0[4];}
      if( ierror==4 ){ C=dA[2]/(1.-A0[2]); A[0]-=C*A0[0]; A[1]-=C*A0[1];
                       A[2]+=  dA[2];      A[3]-=C*A0[3]; A[4]-=C*A0[4];}
      if( ierror==5 ){ C=dA[3]/(1.-A0[3]); A[0]-=C*A0[0]; A[1]-=C*A0[1];
                       A[2]-=C*A0[2];      A[3]+=  dA[3]; A[4]-=C*A0[4];}
      if( ierror==6 ){ C=dA[4]/(1.-A0[4]); A[0]-=C*A0[0]; A[1]-=C*A0[1];
                       A[2]-=C*A0[2];      A[3]-=C*A0[3]; A[4]+=  dA[4];}
      if( ierror==7 ) tau[0] += dtau[0];
      if( ierror==8 ) tau[1] += dtau[1];
      if( ierror==9 ) tau[2] += dtau[2];
      if( ierror==10) tau[3] += dtau[3]; 
      if( ierror==11) tau[4] += dtau[4];
      if( ierror==12) fluence *= 1.045;
      if( ierror==13 && Tsensor0 < 10.) Tsensor += 1.0;
      if( ierror==14 && Tsensor0 >=10.) Tsensor += 2.0;
      double Arrhenius = exp(EI/kB*(1/TR-1/(Tsensor+T0)));
      isum[ierror] = 0.;
//------calculate each component of leak current-----------------
      for (int k=0; k<5; k++) {
        double leak_week = alpha*Vsensor*A[k]*tau[k]/Arrhenius/days
                  *(1-exp(-Arrhenius*days/tau[k]))*fluence;
        ileak[layer][iw][k][ierror] = ileak[layer][iw-1][k][ierror] * 
             exp(-Arrhenius*days/tau[k]) + leak_week;
        isum[ierror] += ileak[layer][iw][k][ierror];
//        cout<<"iw,layer,k,ierror="<<iw<<" "<<layer<<" "<<k<<" "
//        <<Arrhenius<<" days="<<days<<" "<<tau[k]<<" LC="
//        <<leak_week*1e9<<" "<<ileak[layer][iw][k][ierror]*1e9<<endl;   
      }
//--------------------------------------------------------------
      differenceError[ierror] = 0.;
      if(isum[0] != 0.) differenceError[ierror] = (isum[ierror]-isum[0])/isum[0];
      quadSum += pow(differenceError[ierror],2);
//      cout<<"iw, layer, ierror, %error="<<iw<<" "<<layer<<" "<<ierror
//          <<" "<<differenceError[ierror]*100.<<"%"<<endl;
    }
//---------end of error estimate loop---------------------------
//---------shift to operation temperature-----------------------
    double Tplot = Tplot0[layer] - Tbias;
    double W=pow((Tplot+T0)/(T0-7),2)*exp(Eg/2./kB*(1./(T0-7.)-1./(Tplot+T0)));
    leak[layer][iw] = isum[0] * 1.e9 * W ;
    dleak[layer][iw] = sqrt(quadSum);
  }
  cout<<"iweek, weekL="<<iw<<" "<<weekL[iw]<<" : leak (=" <<Tplot<<"C) "
  <<leak[0][iw]<<" "<<leak[1][iw] <<" "<<leak[2][iw]
  <<" "<<leak[3][iw]<<endl;
  cout<<"             "<<iw<<" "<<weekL[iw]<<" : error="
  <<dleak[0][iw]<<" "<<dleak[1][iw] <<" "<<dleak[2][iw]
  <<" "<<dleak[3][iw]<<endl;
}

//---------readin leak current data-------------------
int j=0,idata[4]={0,0,0,0};
int bec, layer;
double tempdata,temptime,dtime[4][100],data[4][100],error;
ifstream fin("leak_data.txt");
while(fin>>bec>>layer>>year2>>month2>>day2>>tempdata>>error)
{
//cout<<bec<<" "<<layer<<" "<<year2<<" "<<month2<<" "<<day2<<" "<<tempdata<<endl;
if(bec!=1)continue;
int L=layer-3;
if(L<0 || L>3) continue;
int m = idata[L];
TDatime da(year2,month2,day2,12, 0,0);
dtime[L][m]=da.Convert();
//   shift data(-10C) to Tplot----
double Tplot = Tplot0[L] - Tbias;
double W=pow((Tplot+T0)/(T0-10.),2)*exp(Eg/2./kB*(1./(T0-10.)-1./
(Tplot+T0)));
data[L][m]=tempdata * W;
cout<<L<<" "<<m<<" "<<dtime[L][m]<<" "<<data[L][m]<<endl;
idata[L]++;
}
fin.close();
//-----------------------------------------------------
//
leg0 = new TLegend(0.65 ,0.20 ,0.88,0.40);
for (int layer=3; layer>-1; layer--){
  for (int iw=0; iw<iweek; iw++) leak[layer][iw] /= 1000.;
  hc = new TGraph(iweek,time,leak[layer]);
  hc->SetMaximum(ymax);
  hc->SetMinimum(ymin);
//  hc->GetYaxis()->SetTitle("I _{leak} [#muA] @-10#circC");
  hc->GetYaxis()->SetTitle("I _{leak} (#muA) @ operation");
  //hc->GetXaxis()->SetTitle("DATE");
  for (int j=0; j<idata[layer]; j++) data[layer][j] /= 1000.;
  g2 = new TGraph(idata[layer],dtime[layer],data[layer]);
  hc->GetXaxis()->SetTitleOffset(1.4);
  hc->GetXaxis()->SetLabelOffset(0.05);
  hc->GetXaxis()->SetLabelSize(0.04);
  hc->GetYaxis()->SetTitleSize(0.05);
  hc->GetYaxis()->SetTitleOffset(1.0);
  hc->GetYaxis()->SetLabelOffset(0.01);
//  hc->GetXaxis()->SetRangeUser(13,45);
  hc->GetXaxis()->SetTimeDisplay(1);
  hc->GetXaxis()->SetTimeFormat("#splitline{  %Y}{%b %d}");
  hc->GetXaxis()->SetLimits(t1,t2);
//  hc->GetXaxis()->SetNdivisions(206);
  hc->GetXaxis()->SetNdivisions(206);
  g2->SetMarkerSize(0.9); hc->SetMarkerSize(0.9);
    hc->SetLineColor(layer+2); hc->SetMarkerColor(layer+2);g2->SetMarkerColor(layer+2);
if(layer==3) {
    hc->SetLineColor(1); hc->SetMarkerColor(1); g2->SetMarkerColor(1);}

//  hu->SetFillStyle(3001);
  if(layer==3)hc->Draw("AL");
  if(layer!=3)hc->Draw("SAME L");
//  hd->SetFillColor(0);
//  hd->Draw("SAME LF");
  g2->Draw("SAME P");
  if(layer==0)leg0->AddEntry(hc,"    B3","LP");
  if(layer==1)leg0->AddEntry(hc,"    B4","LP");
  if(layer==2)leg0->AddEntry(hc,"    B5","LP");
  if(layer==3)leg0->AddEntry(hc,"    B6","LP");
}
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]);

gT6->SetLineStyle(2); gT6->SetLineColor(1);gT6->Draw("SAME L");
gT5->SetLineStyle(2); gT5->SetLineColor(4);gT5->Draw("SAME L");
gT4->SetLineStyle(2); gT4->SetLineColor(3);gT4->Draw("SAME L");
gT3->SetLineStyle(2); gT3->SetLineColor(2);gT3->Draw("SAME L");
//
p = new TPave(0.2,225,7.2,345); 
//------------------integrated luminosity--------------------
gL = new TGraph(iweek,time,intL);
gL->SetLineStyle(2); gL->SetLineColor(46);gL->Draw("SAME L");
leg0->AddEntry(gL,"int-Lumi","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,1000.,t2, ymax,-10.,30.,404,"+L");
axis->SetTitle("T_{sensor} (#circC)");
axis->SetLabelOffset(0.01);
axis->SetTitleOffset(1.2);
axis->SetTitleSize(0.040);
axis->SetLabelSize(0.04);
axis->SetLabelFont(42);
axis->SetTitleFont(42);
axis->Draw();
//
//draw an axis on the right side
TGaxis *axis = new TGaxis(t2,0.1,t2, 1000.,0.1,1000.,505,"+LG");
axis->SetLabelOffset(0.01);
axis->SetTitleOffset(1.0);
axis->SetTitleSize(0.040);
axis->SetLabelSize(0.0);
axis->SetLabelFont(42);
axis->SetTitleFont(42);
axis->Draw();

//-- write last and first days
TDatime db(2021,09,01,00,00,00);int t4 = db.Convert();
legA = new TText(t4 ,0.04 ,"2022");
legA->SetTextSize(0.040);legA->Draw();
TDatime db(2021,08,01,00,00,00);int t4 = db.Convert();
legB = new TText(t4 ,0.0225 ,"Jan 01");
legB->SetTextSize(0.040);legB->Draw();
//
TDatime db(2009,09,01,00,00,00);int t4 = db.Convert();
legA = new TText(t4 ,0.04 ,"2010");
legA->SetTextSize(0.040);legA->Draw();
TDatime db(2009,08,01,00,00,00);int t4 = db.Convert();
legB = new TText(t4 ,0.0225 ,"Jan 01");
legB->SetTextSize(0.040);legB->Draw();
//
TDatime db(2015,01,01,00,00,00);int t4 = db.Convert();
legA = new TLatex(t4 ,6000 ,"sensor temperature");
legA->SetTextSize(0.040);legA->Draw();
//
TDatime db(2018,06,01,00,00,00);int t4 = db.Convert();
legA = new TLatex(t4 ,5 ,"Cooling plan 1");
legA->SetTextSize(0.050);legA->Draw();
//
c1->Print("leak_current_10years_plan1.png");
//c1->Print("leak_current_new.eps");
//c1->Print("leak_current_new.pdf");
}


: