{
//===========================================================================
//   Sample ROOT programme for Vd evolution, Taka Kondo (KEK, ATLAS) 2011.10.24
//===========================================================================
gROOT->LoadMacro("AtlasUtils.C");
c1 = new TCanvas("c1","A Simple Graph",500,30,700,500);
gPad->SetLogy();
//------NIEL = 1 MeV n-eq fluence [/cm^2] at the sensor for 1 nb^-1 pp collisions ------------
Double_t NIEL= 1.65e5, dNIEL= 0.2 * NIEL; // e.g. ATLAS B3, error=20% is a pure guess
//Double_t NIEL= 21.50e5, dNIEL= 0.2 * NIEL; // e.g. ATLAS Pixel B0, error=20% is a pure guess
//
//---------set start and end days , yminPlot and ymaxPlot for plot --------------- 
TDatime db(2010,07,02,00,00,00);
Int_t t1 = db.Convert();
TDatime db(2011,08,14,23,00,00);
Int_t t2 = db.Convert();
Double_t yminPlot = 0.0001, ymaxPlot= 10.0; // min max of plot
//
//===========================================================================
//
//----------------------------------------------------------------
Int_t n = 3000;
Double_t time[n], weekL[n], dweekL[n], intL[n], Temp[n], dTemp[n];
Double_t leak[n], dleak[n], leakop[n], dleakop[n], leakP[n];
//
//--------readin days_luminosity_temperature.txt------------------
//
Int_t year1,month1,day1,hour1,min1,year2,month2,day2,hour2,min2;
Int_t iweek=1;
ifstream fin("days_lumi_temp_Plan1.txt");
char line[255];
fin.getline(line,sizeof(line));   //------read 1st line---------
while(fin>>year1>>month1>>day1>>hour1>>min1>>year2>>month2>>day2>>hour2>>min2>>weekL[iweek]>>dweekL[iweek]>>intL[iweek]>>Temp[iweek]>>dTemp[iweek])
{
  if(iweek<10) cout<<year1<<"-"<<month1<<"-"<<day1<<" "<<hour1<<":"<<min1<<" "<<year2<<"-"
  <<month2<<"-"<<day2<<" "<<hour2<<":"<<min2<<"  "<<weekL[iweek]<<" "
  <<dweekL[iweek]<<" "<<intL[iweek]<<" "<<Temp[iweek]<<" "<<dTemp[iweek]<<endl;
  if(iweek==1) {
    TDatime d1(year1,month1,day1,hour1,min1,0.);
    time[0]=d1.Convert();
  }
  TDatime d2(year2,month2,day2,hour2,min2,60.);
  time[iweek]=d2.Convert();
  intL[iweek] /= 1000000.;
  iweek++;
}
fin.close();
//
//---------Calculate leak current prediction -----------------
//-------- based on R. Harper, Thesis of University of Sheffield, Oct. 2001, page 109-----
//
Double_t A0[5] = {0.42, 0.10, 0.23, 0.21, 0.04}, A[5];
Double_t dA[5] = {0.11, 0.01, 0.02, 0.02, 0.03};
Double_t tau0[5] = {833.33, 28.47, 2.57, 0.086, 0.0056}, tau[5];
Double_t dtau[5] = {138.89, 4.17, 0.21, 0.017, 0.0035};
Double_t Eg = 1.21, EI=1.09, kB=8.617343E-5, TR=293.15, T0=273.15; 
Double_t alpha0 = 6.90e-18, dalpha0 = 0.20e-18, Talpha= -7.0; // value at -7 degreeC 
Double_t ileak[n][5][15];
for (Int_t k=0; k<5; k++) {
  for (Int_t ierror=0; ierror<15; ierror++) ileak[0][k][ierror] = 0.;
}
//
//---------------------go over all weeks-----------------------------
//
for (Int_t iw=1; iw<iweek; iw++) {
  Double_t days = (time[iw] - time[iw-1])/(24.*60.*60);
  Double_t luminosity0 = weekL[iw];
  Double_t dluminosity0 = dweekL[iw];
  Double_t Tsensor0 = Temp[iw]; 
  Double_t dTsensor0 = dTemp[iw];
//
//---------start loop for error estimate (erro item = 0 to 14)-----------------
//  (ierror==0) central value
//  (ierror==1) alpha = aplha0 + dalpha0
//  (ierror==2) A[0] = A0[0] + dA[0]
//  (ierror==3) A[1] = A0[1] + dA[1]
//  (ierror==4) A[2] = A0[2] + dA[2]
//  (ierror==5) A[3] = A0[3] + dA[3]
//  (ierror==6) A[4] = A0[4] + dA[4]
//  (ierror==7) tau[0] = tau[0] + dtau[0]
//  (ierror==8) tau[1] = tau[1] + dtau[1]
//  (ierror==9) tau[2] = tau[2] + dtau[2]
//  (ierror==10) tau[3] = tau[3] + dtau[3]
//  (ierror==11) tau[4] = tau[4] + dtau[4]
//  (ierror==12) luminosity = luminosity0 + dluminosity0
//  (ierror==13) Tsensor = Tsensor + dTsensor0
//  (ierror==14) fluence = NIEL + dNIEL;
//
    Double_t isum[15];
    Double_t quadSum = 0.;
    for (Int_t ierror=0; ierror<15; ierror++) {
//----set back to central values for each time----
      Double_t C;
      Double_t alpha = alpha0;
      Double_t luminosity = luminosity0;
      Double_t fluence = NIEL;
      Double_t Tsensor = Tsensor0; 
      for(Int_t k=0;k<5;k++) A[k] = A0[k];
      for(Int_t k=0;k<5;k++) tau[k] = tau0[k];
//----change corresponding item-----------------
      switch(ierror){  
            case 1: alpha += dalpha0; break;
            case 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]; break;}
            case 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]; break;}
            case 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]; break;}
            case 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]; break;}
            case 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]; break;}
            case 7: { tau[0] += dtau[0]; break;}
            case 8: { tau[1] += dtau[1]; break;}
            case 9: { tau[2] += dtau[2]; break;}
            case 10: { tau[3] += dtau[3]; break;}
            case 11: { tau[4] += dtau[4]; break;}
            case 12: { luminosity += dluminosity0; break;}
            case 13: { Tsensor += dTsensor0; break;}
            case 14: { fluence += dNIEL; break;}
            default: { break ;}
      }
      Double_t Arrhenius = exp(EI/kB*(1/TR-1/(Tsensor+T0)));
      isum[ierror] = 0.;
//------calculate each component of leak current-----------------
      for (Int_t k=0; k<5; k++) {
        Double_t leak_week = alpha * A[k] * tau[k] / Arrhenius / days
                  *(1-exp(-Arrhenius*days/tau[k])) * luminosity * fluence;
        ileak[iw][k][ierror] = ileak[iw-1][k][ierror] * 
             exp(-Arrhenius*days/tau[k]) + leak_week;
        isum[ierror] += ileak[iw][k][ierror];  
      }
//--------------------------------------------------------------
      quadSum += pow(isum[ierror]-isum[0],2);
      if(iw==94) 
        cout<<"iw, ierror="<<iw<<" "<<ierror<<" : "
        <<isum[ierror]<<" "<<(isum[ierror]-isum[0])/isum[0]<<endl;
    }
//---------end of error estimate loop---------------------------
    Double_t Tplot = 0.;
    Double_t W=pow((Tplot+T0)/(T0+Talpha),2)*exp(Eg/2./kB*(1./(T0+Talpha)-1./(Tplot+T0)));
    leak[iw] = isum[0] * 1.e9 * W ;
    dleak[iw] = sqrt(quadSum) * 1.e9 * W; 
    if(iw==94) cout<<"iw, ="<<iw<<" total error="<<dleak[iw]/leak[iw]<<endl; 
}
//
//---------readin leak current data-------------------
//
Int_t idata=0;
cout<<idata<<" "<<year2<<endl;
Double_t tempdata,temptime,dtime[500],data[500],ddata[500];
ifstream fin("leak_data_sample.txt");
char line[255];
fin.getline(line,sizeof(line));   //------readr 1st line---------
while(fin>>year2>>month2>>day2>>data[idata]>>ddata[idata])
{
  TDatime da(year2,month2,day2,12, 0,0);
  dtime[idata]=da.Convert();
//  cout<<idata<<" "<<dtime[idata]<<" "<<data[idata]<<endl;
  idata++;
}
fin.close();
//
//-----------------add one more null data for color-area plotting--------------
//
iweek++;
TDatime dc(2011,12,01,00,00,00);
time[iweek-1]=dc.Convert();
leak[iweek-1] = 0.;
//
//-----------------now make plots-------------------------;
//
  for (Int_t iw=0; iw<iweek; iw++) leakP[iw] = (leak[iw] + dleak[iw])/1000.;
  hu = new TGraph(iweek,time,leakP);  // curve for +1 sigma
  for (Int_t iw=0; iw<iweek; iw++) leakP[iw] = (leak[iw] - dleak[iw])/1000.;
  hd = new TGraph(iweek,time,leakP);  // curve for -1 sigma
  hu->SetMaximum(ymaxPlot);
  hu->SetMinimum(yminPlot);
  hu->SetFillColor(kRed-10);
  hu->GetXaxis()->SetTitle("DATE");
  hu->GetYaxis()->SetTitle("I _{leak}@0#circC (#muA/cm^{3}) ");
  hu->GetXaxis()->SetLabelSize(0.04);   hu->GetYaxis()->SetLabelSize(0.05);
  hu->GetXaxis()->SetLabelOffset(0.03); hu->GetYaxis()->SetLabelOffset(0.01);
  hu->GetXaxis()->SetTitleSize(0.045);  hu->GetYaxis()->SetTitleSize(0.05);
  hu->GetXaxis()->SetTitleOffset(1.5);  hu->GetYaxis()->SetTitleOffset(1.2);
  hu->GetXaxis()->SetTimeDisplay(1);
  hu->GetXaxis()->SetTimeFormat("#splitline{  %Y}{%b %d}");
  hu->GetXaxis()->SetLimits(t1,t2);
  hu->Draw("ALF");
  hd->SetFillColor(0);
  hd->Draw("SAME LF");
//---------plot center values-------------
  for (Int_t iw=0; iw<iweek; iw++) leakP[iw] = leak[iw]/1000.;  
  hc = new TGraph(iweek,time,leakP);  // curve for center values
  hc->SetLineColor(kRed+2);
  hc->Draw("SAME L");
//-------- overlay integrated luminosity--------------
gL = new TGraph(iweek,time,intL);
gL->SetLineStyle(3); gL->SetLineColor(1);gL->Draw("SAME L");
//-------draw legend---------------
leg = new TLegend(0.20 ,0.80 ,0.5,0.90);
leg->AddEntry(gL,"Integrated Luminosity (fb^{-1})","l");
leg->AddEntry(hc,"  #pm1#sigma","L");
leg->SetFillColor(10); leg->SetTextSize(0.04);leg->SetBorderSize(0);
leg->Draw(); 
leg = new TLegend(0.33 ,0.80 ,0.45,0.85);
leg->AddEntry(hu,"  prediction","F");
leg->SetFillColor(10); leg->SetTextSize(0.04);leg->SetBorderSize(0);
leg->Draw();
//----------draw data points------------------------
  for (Int_t j=0; j<idata; j++) {data[j] /= 1000., ddata[j] /=1000.;}
  g2 = new TGraphErrors(idata,dtime,data,0,ddata);
  g2->SetMarkerColor(kRed+2);hd->SetLineColor(kRed+2);
  g2->SetMarkerSize(1.0);
  g2->Draw("SAME P");
//--- write tick mark on bottom X axis-------------------
for (Int_t month=8; month<13; month++){
  TDatime db(2010,month,01,00,00,00); Int_t t4 = db.Convert();
  tick = new TLine(t4, yminPlot, t4, yminPlot*1.4);tick->Draw();
}
for (Int_t month=1; month<9; month++){
  TDatime db(2011,month,2,00,00,00); Int_t t4 = db.Convert();
  tick = new TLine(t4, yminPlot, t4, yminPlot*1.4);tick->Draw();
}
//-----------x bottom axis, draw once morer------------------------
tick = new TLine(t1, yminPlot, t2, yminPlot); tick->SetLineStyle(1);tick->Draw();
//-----------draw an axis on the right side---------------------
TGaxis *axis = new TGaxis(t2,yminPlot,t2, ymaxPlot,yminPlot,ymaxPlot,505,"+LG");
axis->SetLabelSize(0.0);
axis->Draw();
//-----------------------print -------------------------
c1->Print("leak_current_sample.png");
c1->Print("leak_current_sample.pdf");
}


: