{
//-----plot_Dortmund_with_noError_1.C------------------------
//    2012.1.26  change log term correctly
//    2012.1.30 bug fixed
//    2012.1.31 add hour and minutes in input
//    2012.2.22 added both polarity in error, 4.5%->3.7% in fluence error
//    2012.3.1  data and model with 1.21
//    2012.3.3. Atlas Perliminary
//    2012.3.15 int. lum deliv. -> right axis
//    2012.3.17 Y axis title: leak-> vol
//    2012.3.22 change to sigma band and cleanup 
//-----------------------------------------------
//-----------------------------------------------
#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,30,700,500);
//c1->SetGrid();
//TDatime db(2010,06,01,00,00,00);
TDatime db(2010,07,01,00,00,00);
int t1 = db.Convert();
//TDatime db(2011,04,01,00,00,00);
TDatime db(2012,1,01,23,59,59);
int t2 = db.Convert();
gStyle->SetOptStat(0);
gPad->SetLogy();
int isGrid = 1;
//c1->SetGrid();
//------------------create the arrays for the points-------------------------
char cid[40];
double Tplot = 0.; 
Double_t ymax= 1000.0;
Double_t ymin= 0.001;
double Tscale=ymax*0.1/10., Toffset = ymax*0.2;
//double Tscale=0.01/10., Toffset = 0.02;
Int_t n = 3000; // number of weeks
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];
for (int i=0; i<n; i++) {
  for (int k=0; k<4; k++) {leak[k][i]=0.; dleak[k][i]= 0.;}
}
double leak_exp[4][n],dleak_exp[4][n];
Int_t year1,month1,day1,hour1,min1,year2,month2,day2,hour2,min2;
double d = 0.0285, width = 6.144, length = 6.136, nsensor = 4.0;
double Ssensor = width*length*nsensor;
double Vsensor = Ssensor * d;
//---------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_1.21eV.txt");
while(fin>>bec>>layer>>year2>>month2>>day2>>hour2>>min2>>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,hour2, min2,0);
dtime[L][m]=da.Convert();
data[L][m]=tempdata/Vsensor;
cout<<L<<" "<<m<<" "<<dtime[L][m]<<" LC="<<data[L][m]<<endl;
idata[L]++;
}
fin.close();
//-----------------------------------------------------
//--------readin luminosity data-------------------
int iweek=1, iT=0;
ifstream fin("lhcall_7_temperature_hour.txt");
while(fin>>year1>>month1>>day1>>hour1>>min1>>year2>>month2>>day2>>hour2>>min2>>energy[iweek]>>
weekL[iweek]>>intL[iweek]>>Temp[0][iweek]>>Temp[1][iweek]>>Temp[2][iweek]>>Temp[3][iweek] )
{
if(iweek<10) 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]/100.; 
  TDatime d1(year1,month1,day1,hour1,min1,0);
  int time1=d1.Convert();
  if(time1 > t2) continue;
  if(iweek==1) time[0] = time1;
  TDatime d2(year2,month2,day2,hour2,min2,60);
  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;
    Temp_plot[layer][iT] = pow(10., 2.+(Temp[layer][iweek]+10.)/40.);
  }
  iT++;
  Ttime[iT]=time[iweek];
  for(int layer=0; layer<4; layer++) 
    //Temp_plot[layer][iT] = Temp[layer][iweek]*Tscale+Toffset;
    Temp_plot[layer][iT] = pow(10., 2.+(Temp[layer][iweek]+10)/40.);
  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 fluence0, dfluence=0.037;
double Tsensor, Tsensor0;
double EI0=1.09, dEI=0.14, kB=8.617343E-5, T0=273.15;
double Eg0 = 1.21, dEg_plus=0.04, dEg_minus=0.09;
//-----------Dortmund model parameters-------------
double alphaI0 = 1.23e-17, dalphaI = 0.06e-17;
double k0I0=1.2E13*60.*60.*24., dk0I_plus=5.3E13*60.*60.*24., 
       dk0I_minus=1.0E13*60.*60.*24.;
double EI0=1.11, dEI=0.05;
double alpha0star0 = 7.07E-17, dalpha0star = alpha0star0*dalphaI/alphaI0;
double beta0 = 3.29E-18, dbeta = 0.3E-18; 
double EIstar0 = 1.30, dEIstar = 0.14;
double Tref = T0 + 21.;
double t0 = 1./(24.*60.);  // 1 min. for log term
//--------------------------------------------------
int nError=10,  n2Error = nError*2; 
double ileak_exp[4][n][n2Error],ileak_log[4][n][n2Error] ;
double days_eff[4][n][n2Error], fluences[4][n][n2Error];
for (int i=0; i<n; i++) {
  for(int j=0; j<n2Error; j++) {
     for (int layer=0; layer<4; layer++) {
        ileak_exp[layer][i][j] = 0.;
        ileak_log[layer][i][j] = 0.;
        days_eff[layer][i][j] = 0.;
        fluences[layer][i][j] = 0.;
     }
  }
}
time[0] = time[1] - 24.* 60.*60;
//---------------------go over all weeks-----------------------------
for (int iw=1; iw<iweek; iw++) {
  double days = (time[iw] - time[iw-1])/(24.*60.*60);  // day in current week
  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 nError-1)------------------------
    double itotal=0., quadSum = 0.;
    for (double pol=-1.; pol<2.; pol+=2.) {
      for (int jerror=0; jerror<nError; jerror++) {
        int ierror = jerror;
        if( pol > 0.) ierror = jerror + nError;
//---------------------------------------------------
        double alphaI = alphaI0;
        double k0I = k0I0;
        double EI = EI0;
        double alpha0star = alpha0star0;
        double beta = beta0;
        double EIstar = EIstar0; 
        Tsensor = Tsensor0; 
        double fluence_thisweek = fluence0;
        double Eg = Eg0;
//
        if( jerror==1 ) alphaI += pol*dalphaI;
        if( jerror==2 && pol>0.) k0I += dk0I_plus;
        if( jerror==2 && pol<0.) k0I -= dk0I_minus;
        if( jerror==3 ) EI += pol*dEI;
        if( jerror==4 ) alpha0star += pol*dalpha0star;
        if( jerror==5 ) beta += pol * dbeta;
        if( jerror==6 ) EIstar += pol * dEIstar;
        if( jerror==7 && pol<0.) Eg -= dEg_minus;
        if( jerror==7 && pol>0.) Eg += dEg_plus;
        if( jerror==8 && Tsensor0 < 10.) Tsensor += pol * 1.0;
        if( jerror==8 && Tsensor0 >=10.) Tsensor += pol * 2.0;
        if( jerror==9 ) fluence_thisweek *= (1. + pol*dfluence) ;
        fluences[layer][iw][ierror] = fluence_thisweek;
//-----leak_week = new leakcurrent due to new flux during this week -----------
        double Arrhenius = exp(-EIstar/kB*(1/(Tsensor+T0)-1/Tref));
        double W=pow((Tplot+T0)/Tref,2)*exp(-Eg/2./kB*(1./(Tplot+T0)-1./Tref));
        double tauI = 1./ (k0I * exp(-EI/kB/(Tsensor+T0)));
        double leak_week_exp = W * fluence_thisweek * alphaI 
                             *exp(-0.5*days/tauI);
        double leak_week_log = W * fluence_thisweek * 
                             (alpha0star - beta * log(0.368*days/t0));
        ileak_exp[layer][iw][ierror] = leak_week_exp 
                      + ileak_exp[layer][iw-1][ierror] * exp(-days/tauI);
        days_eff[layer][iw][ierror] = days * Arrhenius;
        double sum_fluence = 0.;
        double beta_term = 0.;
        if( iw != 1) {
           for (int jw = 1; jw<iw; jw++) {
             sum_fluence += fluences[layer][jw][ierror];
             double daysij = 0.;
             for (int kw=jw+1; kw<iw+1; kw++) daysij += days_eff[layer][kw][ierror];
             beta_term += fluences[layer][jw][ierror]*log(daysij/t0); 
           }
        }
        ileak_log[layer][iw][ierror] = leak_week_log + 
              W * (alpha0star*sum_fluence - beta*beta_term); 
//--------------------------------------------------------------
        double diff = 0.;
        itotal = ileak_exp[layer][iw][0] + ileak_log[layer][iw][0];
        if(iw==100 && layer==0) {
         cout<<"iw,layer,ierror,pol="<<iw<<" "<<layer<<" "<<" "<<ierror<<" "<<pol<<
          " (exp,log)=("<<leak_week_exp*1.e9<<","<<leak_week_log*1.e9<<
          ")-->"<<ileak_exp[layer][iw][ierror]*1.e9<<", "
          <<ileak_log[layer][iw][ierror]*1.e9<<" ="<<itotal*1.e9<<endl; 
         }   
        if(itotal != 0.) diff += sqrt(
         (pow(ileak_exp[layer][iw][ierror] + ileak_log[layer][iw][ierror]-itotal,2) 
        +pow(ileak_exp[layer][iw][ierror] + ileak_log[layer][iw][ierror]-itotal,2))         *0.5)/itotal;
//      ------factor 0.5 because +1 and -1 sigma cases----
         quadSum += pow(diff,2);
//         cout<<"iw,layer,ierror,%error="<<iw<<" "<<layer<<" "<<ierror
//             <<" "<<diff*100.<<"%"<<endl;
      }  //---------end of error estimate loop----------------------
    }  // loop end for polarity loop

//---------end of error estimate loop---------------------------
    leak[layer][iw] = itotal * 1.e9  ;
    dleak[layer][iw] = itotal* sqrt(quadSum) * 1.e9  ;
  }
//  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;
}
//-----------------now make plots-------------------------
legd = new TLegend(0.55 ,0.18 ,0.86,0.40);
leg0 = new TLegend(0.56 ,0.19 ,0.68,0.39,"Layer");
leg2 = new TLegend(0.62 ,0.19 ,0.83,0.39,"Data");
leg1 = new TLegend(0.68 ,0.19 ,0.85,0.39,"model prediction");
leg4 = new TLegend(0.80 ,0.19 ,0.84,0.35);
for (int layer=0; layer<4; layer++){
  double xsigma[2*iweek], ysigma[2*iweek];
  for (int iw=0; iw<iweek; iw++) {
      xsigma[iw] = time[iw]; xsigma[2*iweek-1-iw] = time[iw];
      ysigma[iw] = (leak[layer][iw]+dleak[layer][iw])/1000.;
      ysigma[2*iweek-1-iw] =  (leak[layer][iw]-dleak[layer][iw])/1000.;
  }
  hu = new TGraph(2*iweek,xsigma,ysigma);
  hu->SetMaximum(ymax);
  hu->SetMinimum(ymin);
  sprintf(cid,"I_{vol}@%d^{o}C (#muA/cm^{3})\0",(int)Tplot); 
  hu->GetYaxis()->SetTitle(cid);
  hu->GetXaxis()->SetTitle("Date");
  hu->GetXaxis()->SetTitleOffset(1.5);
  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()->SetTimeDisplay(1);
  hu->GetXaxis()->SetTimeFormat("#splitline{  %Y}{%b %d}");
  hu->GetXaxis()->SetLimits(t1,t2);
  hu->GetXaxis()->SetNdivisions(207);
  int kColor=kRed;
  if(layer==1) kColor = kGreen;
  if(layer==2) kColor = kBlue;
  if(layer==3) kColor = kMagenta;
  hu->SetFillColor(kColor-9);
  hu->SetFillStyle(3001);
  hu->GetXaxis()->SetNdivisions(0);
  if(layer==0)hu->Draw("ALF");
  if(layer!=0)hu->Draw("SAME LF");
  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","");
}
//----------draw central values and data points -------------------------
for (int layer=0; layer<4; layer++){
  double yplot[iweek];
  for (int iw=0; iw<iweek; iw++) yplot[iw] = leak[layer][iw]/1000.;
  hc = new TGraph(iweek,time,yplot);
  g2 = new TGraph(idata[layer],dtime[layer],data[layer]);
  int kColor=kRed;
  if(layer==1) kColor = kGreen;
  if(layer==2) kColor = kBlue;
  if(layer==3) kColor = kMagenta;
  hc->SetLineColor(kColor+2);
  g2->SetMarkerColor(kColor+2);
  hc->Draw("SAME L");
  g2->SetMarkerSize(0.8);
  g2->Draw("SAME P");
  leg1->AddEntry(hc,"  #pm1#sigma","L");
  leg2->AddEntry(g2,"   ","P");
}
//-------------------------plot temperature--------------------------
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(1); gT3->SetLineColor(kRed+2);gT3->Draw("SAME L");
gT4->SetLineStyle(1); gT4->SetLineColor(kGreen+2);gT4->Draw("SAME L");
gT5->SetLineStyle(1); gT5->SetLineColor(kBlue+2);gT5->Draw("SAME L");
gT6->SetLineStyle(1); gT6->SetLineColor(kMagenta+2);gT6->Draw("SAME L");
//-------------------plot integrated luminosity---------------
for (int i=0; i<iweek; i++) intL[i]/=1000.;
hintL = new TGraph(iweek,time,intL);
hintL->Draw("SAME L");
TDatime db(2010,11,15,00,00,00);
int t4 = db.Convert();
lum = new TLatex(t4 ,0.60 ,"int. lum. deliv.");
lum->SetTextSize(0.035);
lum->Draw();
//-------------------- write years and days on X label--------------
int day1=16;
double ylabel1= ymin*0.4, ylabel2=ymin*0.2;
for (int m=1; m<13; m+=3){
  TDatime db(2010,m-1,day1,00,00,00);int t4 = db.Convert();
  if(t4>t1-3000000) {
    sprintf(cid,"%d/01\0",m); legA = new TText(t4 ,ylabel1 ,cid);
    legA->SetTextSize(0.040); legA->Draw();
    legA = new TText(t4 ,ylabel2 ,"2010");
    legA->SetTextSize(0.040);legA->Draw();
  }
  TDatime db(2011,m-1,day1,00,00,00);int t4 = db.Convert();
  sprintf(cid,"%d/01\0",m); legA = new TText(t4 ,ylabel1 ,cid);
  legA->SetTextSize(0.040);legA->Draw();
  legA = new TText(t4 ,ylabel2 ,"2011");
  legA->SetTextSize(0.040);legA->Draw();
  TDatime db(2012,m-1,day1,00,00,00);int t4 = db.Convert();
  if(t4<t2) {
     sprintf(cid,"%d/01\0",m); legA = new TText(t4 ,ylabel1 ,cid);
     legA->SetTextSize(0.040); if( t4<t2) legA->Draw();
     legA = new TText(t4 ,ylabel2 ,"2012");
     legA->SetTextSize(0.040);legA->Draw();
  }
//----grid-----
  if(isGrid !=0) {
     TDatime db(2010,m,01,00,00,00); t4 = db.Convert();
     tick = new TLine(t4, 0, t4, ymax); tick->SetLineStyle(3);
     if(t4>t1)tick->Draw();
     TDatime db(2011,m,01,00,00,00); t4 = db.Convert();
     tick = new TLine(t4, 0, t4, ymax); tick->SetLineStyle(3);tick->Draw();
   }
}
//--- tick marks------
for (int i=1;i<13;i++) {
 TDatime db(2010,i,01,00,00,00); int t4 = db.Convert();
 if(t4>t1){
 tick = new TLine(t4, ymin, t4, ymin*1.3); tick->SetLineStyle(1);tick->Draw();
 tick = new TLine(t4, ymax*0.8, t4, ymax); tick->SetLineStyle(1);tick->Draw();
 }
 TDatime db(2011,i,01,00,00,00); int t4 = db.Convert();
 tick = new TLine(t4, ymin, t4, ymin*1.3); tick->SetLineStyle(1);tick->Draw();
 tick = new TLine(t4, ymax*0.8, t4, ymax); tick->SetLineStyle(1);tick->Draw();
}
//-----------horizontal bar------------------------
if(isGrid !=0) {
   for (double y=ymin; y<ymax; y*=10.) {
      tick = new TLine(t1, y, t2, y); tick->SetLineStyle(3);tick->Draw();
   }
}
//----------Legend-----------------
legd->SetFillColor(10); legd->SetBorderSize(1); 
leg0->SetFillColor(10); leg0->SetBorderSize(0); leg0->SetTextSize(0.035);
leg1->SetFillColor(10); leg1->SetBorderSize(0); leg1->SetTextSize(0.035);
leg2->SetFillColor(10); leg2->SetBorderSize(0); leg2->SetTextSize(0.035);
leg4->SetFillColor(10); leg4->SetBorderSize(0); leg4->SetTextSize(0.01);
legd->Draw();
leg0->Draw();
leg2->Draw();
leg1->Draw();
leg4->Draw();
//---------------draw Temperature axis on the right side-----------------
TGaxis *axis = new TGaxis(t2,ymax/10.*pow(10.,0.25),t2, ymax,0.,30.,404,"+L");
axis->SetTitle("T (#circC)");
axis->SetLabelOffset(0.01);
axis->SetTitleOffset(1.0);
axis->SetTitleSize(0.035);
axis->SetLabelSize(0.035);
axis->SetLabelFont(42);
axis->SetTitleFont(42);
axis->Draw();
//-----------------------draw an axis on the right side----------------
TGaxis *axis = new TGaxis(t2,ymin,t2, ymax/10.,ymin/10.,ymax/100.,505,"+LG");
axis->SetTitle("Int. Lum. Deliv. (fb^{-1})     ");
axis->SetLabelOffset(0.0);
axis->SetTitleOffset(1.0);
axis->SetTitleSize(0.035);
axis->SetLabelSize(0.035);
axis->SetLabelFont(42);
axis->SetTitleFont(42);
axis->Draw();
//-------------------write sensor temperature ---------------
TDatime db(2011,2,10,00,00,00); int t4 = db.Convert();
TDatime db(2011,7,25,00,00,00); int t5 = db.Convert(); 
p = new TPave(t4, ymax*0.6,t5, ymax*0.95); 
p->SetFillColor(10); p->SetBorderSize(0), p->Draw();
TDatime db(2011,2,15,00,00,00);
int t4 = db.Convert();
lum = new TLatex(t4 ,600 ,"sensor temperature");
lum->SetTextSize(0.040);
lum->Draw();
//----------------- write ATLAS preliminary--------------------------
TDatime db(2010,07,15,00,00,00);
int t5 = db.Convert();
t=new TLatex(t5,ymax*0.05,"ATLAS preliminary");
t->SetTextSize(0.05);
t->Draw();
//-----------------------print -------------------------
c1->Print("Barrel_log_2011end_2.png");
//c1->Print("Barrel_log_2011end_2.pdf");
//c1->Print("Barrel_log_2011end_2.eps");
}