{
//-----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 
//    2012.12.19 change for 2012 and beam data
//    2013.1.16  use fitted values from Temperature/Days/LCfit_selected.txt
//-----------------------------------------------
//-----------------------------------------------
#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,1400,1000);
//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(2013,1,01,23,59,59);
int t2 = db.Convert();
gStyle->SetOptStat(0);
gPad->SetLogy();
int isGrid = 1;
int isOldincluded = 0;
//c1->SetGrid();
//------------------create the arrays for the points-------------------------
char cid[40];
double Tplot = 0.; 
Double_t ymax= 1000.0;
Double_t ymin= 0.0001;
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;
//--------readin luminosity data-------------------
int iweek=1, iT=0;
ifstream fin("lhcall_7_temperature_hour_2012.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]*10.; 
   weekL[iweek] = weekL[iweek]*1000.;
  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 factor8TeV=1.05, factor14TeV=1.24;  // see 2012 Nov 8 talk at SCT Performance meeting
//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 d = 0.0285, width = 6.144, length = 6.136, nsensor = 4.0;
double Ssensor = width*length*nsensor;
double Vsensor = Ssensor * d;
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++) {
    fluence0 = weekL[iw]*Fluence_7TeV[layer];
    if(energy[iw]>7. && energy[iw] < 9.) fluence0 *= factor8TeV;
    if(energy[iw]>13.) fluence0 *= factor14TeV;
    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;
}

//---------readin leak current data-------------------
//double pedestal[4]={0.0230429, 0.0250301, 0.0265945, 0.0316359};
double pedestal[4]={0.0, 0.0, 0.0, 0.0};
int j=0,idata[2][4]={{0,0,0,0},{0,0,0,0}};
int bec, layer, FSIon, Beam, nfit[4];
double temptime,dtime[2][4][1000],data[2][4][1000], error, fitv[4];
//ifstream fin("Temperature/Days/LCfit.txt");
ifstream fin("LCfit.txt");
char line[255]; 
//fin.getline(line,sizeof(line));  // skip 1st line
while(fin>>year2>>month2>>day2>>FSIon>>nfit[0]>>fitv[0]>>nfit[1]>>fitv[1]>>nfit[2]>>fitv[2]>>nfit[3]>>fitv[3]){
   TDatime da(year2,month2,day2,12.,0.,0);
//-------just plot all
   if(FSIon==0 || FSIon!=0) {    
       for (int L=0; L<4; L++) {
            int m = idata[0][L];
            dtime[0][L][m]=da.Convert();
            data[0][L][m]=(fitv[L]-pedestal[L])/Vsensor;
            idata[0][L]++;
       }
    }
   else {
      for (int L=0; L<4; L++) {
            int m = idata[1][L];
            dtime[1][L][m]=da.Convert();
            data[1][L][m]=(fitv[L]-pedestal[L])/Vsensor;
            idata[1][L]++;
       }
   }
}
fin.close();
cout<<"data points(ALL):"<<idata[0][0]<<","<<idata[0][1]<<","<<idata[0][2]<<","<<idata[0][3]<<endl;
cout<<"data points (FSIOFF):"<<idata[1][0]<<","<<idata[1][1]<<","<<idata[1][2]<<","<<idata[1][3]<<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);
  hu->GetXaxis()->SetTitle("Date");
  hu->GetXaxis()->SetTitleOffset(2.0);
  hu->GetXaxis()->SetTitleSize(0.03);
  hu->GetXaxis()->SetLabelOffset(0.05);
  hu->GetXaxis()->SetLabelSize(0.0);
  sprintf(cid,"I_{vol} @ %d^{o}C (#muA/cm^{3})\0",(int)Tplot); 
  hu->GetYaxis()->SetTitle(cid);
  hu->GetYaxis()->SetTitleSize(0.035);
  hu->GetYaxis()->SetLabelSize(0.035);
  hu->GetYaxis()->SetTitleOffset(1.3);
  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 = kCyan+2;
  if(layer==3) kColor = kBlue;
  hu->SetFillColor(kColor-10);
  //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[0][layer],dtime[0][layer],data[0][layer]);
  //g3 = new TGraph(idata[1][layer],dtime[1][layer],data[1][layer]);
  int kColor=kRed;
  if(layer==1) kColor = kGreen;
  if(layer==2) kColor = kCyan+2;
  if(layer==3) kColor = kBlue;
  hc->SetLineColor(kColor+2);
  g2->SetMarkerColor(kColor+2);
  if(layer==2)   g2->SetMarkerColor(kColor);
  hc->Draw("SAME L");
  g2->SetMarkerSize(0.8);
  g2->Draw("SAME P");
  //g3->SetMarkerColor(1);
  //g3->SetMarkerSize(1.5);   // special for FSI OFF points
  //g3->SetMarkerStyle(20);
  //if(isOldincluded == 1) g3->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);gT3->Draw("SAME L");
gT4->SetLineStyle(1); gT4->SetLineColor(kGreen);gT4->Draw("SAME L");
gT5->SetLineStyle(1); gT5->SetLineColor(kCyan+2);gT5->Draw("SAME L");
gT6->SetLineStyle(1); gT6->SetLineColor(kBlue);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,10,1,0,0,0);
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=10;
double ylabel1= ymin*0.5, ylabel2=ymin*0.3;
for (int iy=2010; iy<2014; iy++) {
  for (int m=1; m<13; m+=3){
    TDatime db(iy,m-1,day1,00,00,00);int t4 = db.Convert();
    if(t4<t1-3000000 || t4> t2+3000000) continue;
    sprintf(cid,"%d/01\0",m); legA = new TText(t4 ,ylabel1 ,cid);
    legA->SetTextSize(0.030); legA->Draw();
    sprintf(cid,"%d\0",iy); legA = new TText(t4 ,ylabel2 ,cid);
    legA->SetTextSize(0.030);legA->Draw();
  }
//--- tick marks------
  for (int i=1;i<13;i++) {
     TDatime db(iy,i,01,00,00,00); int t4 = db.Convert();
     if(t4<=t1 || t4>=t2) continue;
     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();
     if(isGrid != 0 && i%3 == 1) {
        tick = new TLine(t4, 0, t4, ymax); tick->SetLineStyle(3);  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.03);
leg1->SetFillColor(10); leg1->SetBorderSize(0); leg1->SetTextSize(0.03);
leg2->SetFillColor(10); leg2->SetBorderSize(0); leg2->SetTextSize(0.03);
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.03);
axis->SetLabelSize(0.03);
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.03);
axis->SetLabelSize(0.03);
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.035);
lum->Draw();
//-----------------------OLD point description------------------
if( isOldincluded!=0) {
    TDatime db(2012,01,15,00,00,00);
    int t5 = db.Convert();
    t=new TLatex(t5,0.1 ,"old method with FSI OFF");
    t->SetTextSize(0.035); t->Draw();
    TDatime db(2012,01,0,00,00,00);
    t5 = db.Convert();
    TMarker *mm=new TMarker(t5,0.13,20);
    mm->SetMarkerSize(1.5); mm->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=new TLatex(t5,ymax*0.05,"preliminary");
t->SetTextSize(0.04);
t->Draw();
//-----------------------print -------------------------
if(isOldincluded==1) c1->Print("Barrel_log_2012end_OLDincluded_fit.png");
else {
  c1->Print("Barrel_log_2012end_fit.png");
  //c1->Print("Barrel_log_2012end_fit.pdf");
  //c1->Print("Barrel_log_2012end_fit.eps");
  } 
}