{
#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,08,01,00,00,00);
int t1 = db.Convert();
//TDatime db(2011,04,01,00,00,00);
TDatime db(2011,05,22,00,00,00);
int t2 = db.Convert();
gStyle->SetOptStat(0);
//c1->SetGrid();
// create the arrays for the points
Double_t ymax= 1.0;
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 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]/1000.; 
  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 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---------------------------
    double Tplot = -10.;
    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();
data[L][m]=tempdata;
cout<<L<<" "<<m<<" "<<dtime[L][m]<<" "<<data[L][m]<<endl;
idata[L]++;
}
fin.close();
//-----------------------------------------------------
//  add one more null data for color-area plotting
iweek++;
TDatime dc(2011,07,01,00,00,00);
time[iweek-1]=dc.Convert();
for (int layer=0; layer<4; layer++) leak[layer][iweek-1] = 0.;
//
legd = new TLegend(0.18 ,0.69 ,0.65,0.91);
leg0 = new TLegend(0.20 ,0.70 ,0.25,0.90,"Layer");
leg2 = new TLegend(0.28 ,0.70 ,0.50,0.90,"Data");
leg3 = new TLegend(0.35 ,0.70 ,0.50,0.90,"T_{sensor}");
leg1 = new TLegend(0.45 ,0.70 ,0.60,0.90,"  Prediction");
leg4 = new TLegend(0.55 ,0.70 ,0.65,0.86);
for (int layer=0; layer<4; layer++){
  for (int iw=0; iw<iweek; iw++) leak[layer][iw] /= 1000.;
  hc = new TGraph(iweek,time,leak[layer]);
  for (int iw=0; iw<iweek; iw++) leak[layer][iw] *= (1.+dleak[layer][iw]);
  hu = new TGraph(iweek,time,leak[layer]);
  for (int iw=0; iw<iweek; iw++) leak[layer][iw] *= (1.-dleak[layer][iw])
                                                    /(1.+dleak[layer][iw]);
  hd = new TGraph(iweek,time,leak[layer]);
  hu->SetMaximum(ymax);
  hu->GetYaxis()->SetTitle("I _{leak} [#muA] @-10#circC");
  //hu->GetXaxis()->SetTitle("DATE");
  for (int j=0; j<idata[layer]; j++) data[layer][j] /= 1000.;
  g2 = new TGraph(idata[layer],dtime[layer],data[layer]);
  hu->GetXaxis()->SetTitleOffset(1.4);
  hu->GetXaxis()->SetLabelOffset(0.05);
  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);
  if(layer==0){
    hc->SetLineColor(kRed+2);hu->SetFillColor(kRed-9);
    g2->SetMarkerColor(kRed+2);hd->SetLineColor(kRed+2);}
  if(layer==1){
    hc->SetLineColor(kGreen+2);hu->SetFillColor(kGreen-9);
    g2->SetMarkerColor(kGreen+2);hd->SetLineColor(kGreen+2);}
  if(layer==2){
    hc->SetLineColor(kBlue+2);hu->SetFillColor(kBlue-9);
    g2->SetMarkerColor(kBlue+2);hd->SetLineColor(kBlue+2);}
  if(layer==3){
    hc->SetLineColor(kMagenta+2);hu->SetFillColor(kMagenta-9);
    g2->SetMarkerColor(kMagenta+2);hd->SetLineColor(kMagenta+2);}
  hu->SetFillStyle(3001);hu->SetLineColor(0);
  if(layer==0)hu->Draw("ALF");
  if(layer!=0)hu->Draw("SAME LF");
  hc->Draw("SAME L");
  hd->SetFillColor(0);
  hd->Draw("SAME LF");
  g2->SetMarkerSize(1.0);
  g2->Draw("SAME P");
  leg1->AddEntry(hd,"  #pm1#sigma","L");
  leg4->AddEntry(hu," ","F");
  if(layer==0)leg0->AddEntry(hu,"B3","");
  if(layer==1)leg0->AddEntry(hu,"B4","");
  if(layer==2)leg0->AddEntry(hu,"B5","");
  if(layer==3)leg0->AddEntry(hu,"B6","");
  leg2->AddEntry(g2,"   ","P");
}
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]);
leg3->AddEntry(gT3," ","L");
leg3->AddEntry(gT4," ","L");
leg3->AddEntry(gT5," ","L");
leg3->AddEntry(gT6," ","L");
gT3->SetLineStyle(2); gT3->SetLineColor(kRed+2);gT3->Draw("SAME L");
gT4->SetLineStyle(2); gT4->SetLineColor(kGreen+2);gT4->Draw("SAME L");
gT5->SetLineStyle(2); gT5->SetLineColor(kBlue+2);gT5->Draw("SAME L");
gT6->SetLineStyle(2); gT6->SetLineColor(kMagenta+2);gT6->Draw("SAME L");
//
p = new TPave(0.2,225,7.2,345); 
//
legd->SetFillColor(10); legd->SetBorderSize(1); 
leg0->SetFillColor(10); leg0->SetBorderSize(0); leg0->SetTextSize(0.04);
leg1->SetFillColor(10); leg1->SetBorderSize(0); leg1->SetTextSize(0.04);
leg2->SetFillColor(10); leg2->SetBorderSize(0); leg2->SetTextSize(0.04);
leg3->SetFillColor(10); leg3->SetBorderSize(0); leg3->SetTextSize(0.04);
leg4->SetFillColor(10); leg4->SetBorderSize(0); leg4->SetTextSize(0.01);
legd->Draw();
leg0->Draw();
leg2->Draw();
leg3->Draw();
leg1->Draw();
leg4->Draw();
//draw an axis on the right side
TGaxis *axis = new TGaxis(t2,0.,t2, ymax,-Toffset/Tscale,(ymax-Toffset)/Tscale,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();
//draw the tick marks again on X axis  
//TGaxis ax:
//ax = hu.Xaxis();
//axis->SetLabelSize(0);
//axis->Draw();
//
//--- write tick mark on bottom X axis
TDatime db(2010,10,01,00,00,00); int t4 = db.Convert();
tick = new TLine(t4, 0, t4, 0.01); tick->Draw();
TDatime db(2010,11,01,00,00,00); int t4 = db.Convert();
tick = new TLine(t4, 0, t4, 0.02); tick->Draw();
TDatime db(2010,12,01,00,00,00); int t4 = db.Convert();
tick = new TLine(t4, 0, t4, 0.01); tick->Draw();
TDatime db(2011,01,01,00,00,00); int t4 = db.Convert();
tick = new TLine(t4, 0, t4, 0.02); tick->Draw();
TDatime db(2011,02,01,00,00,00); int t4 = db.Convert();
tick = new TLine(t4, 0, t4, 0.01); tick->Draw();
TDatime db(2011,03,03,00,00,00); int t4 = db.Convert();
tick = new TLine(t4, 0, t4, 0.02); tick->Draw();
TDatime db(2011,05,03,00,00,00); int t4 = db.Convert();
tick = new TLine(t4, 0, t4, 0.02); tick->Draw();
TDatime db(2011,04,03,00,00,00); int t4 = db.Convert();
tick = new TLine(t4, 0, t4, 0.01); tick->Draw();
TDatime db(2011,05,03,00,00,00); int t4 = db.Convert();
tick = new TLine(t4,0.02,t4,0.34); tick->SetLineStyle(3);tick->Draw();
TDatime db(2011,03,03,00,00,00); int t4 = db.Convert();
tick = new TLine(t4,0.02,t4,0.04); tick->SetLineStyle(3);tick->Draw();
TDatime db(2011,01,01,00,00,00); int t4 = db.Convert();
tick = new TLine(t4,0.02,t4,0.04); tick->SetLineStyle(3);tick->Draw();
TDatime db(2010,11,01,00,00,00); int t4 = db.Convert();
tick = new TLine(t4,0.02,t4,0.05); tick->SetLineStyle(3);tick->Draw();
TDatime db(2011,04,18,00,00,00); int t4 = db.Convert();
tick = new TLine(t4, 0.1, t2, 0.1);tick->SetLineStyle(3); tick->Draw();
TDatime db(2011,04,25,00,00,00); int t4 = db.Convert();
tick = new TLine(t4, 0.2, t2, 0.2);tick->SetLineStyle(3); tick->Draw();
TDatime db(2011,05,01,00,00,00); int t4 = db.Convert();
tick = new TLine(t4, 0.3, t2, 0.3);tick->SetLineStyle(3); tick->Draw();

//

c1->Print("leak_current.png");
c1->Print("leak_current.eps");
c1->Print("leak_current.pdf");
}


: