{
//===========================================================================
//   Sample ROOT programme for Vd evolution, Taka Kondo (KEK, ATLAS) 2011.10.24
//===========================================================================
//------------------set up for plot --------------------------------------
gROOT->LoadMacro("AtlasUtils.C");
ATLAS->SetPadTickY(0);
ATLAS->SetPadRightMargin(0.12);
ATLAS->SetPadLeftMargin(0.12);
c1 = new TCanvas("c1","A Simple Graph",500,30,700,500);
gStyle->SetOptStat(0);
//--------1 MeV n-eq fluence [/cm^2] at 1 nb^-1 7 TeV pp collision------------
//
Double_t NIEL0 = 1.65e5, dNIEL0 = 0.20 * NIEL0; // e.g. ATLAS B3, error=20% is a pure guess
//Double_t NIEL0 = 21.50e5, dNIEL0 = 0.20 * NIEL0; //e.g. ATLAS Pixel B0, error=20% pure guess
Double_t d = 0.0285; // sensor thickness [cm]
Double_t Vd0 = 70., dVd0 = 0.;   // initial voltage and its spread [V]
//---------set start and end days for plot --------------- 
TDatime db(2010,01,01,00,00,00);
Int_t t1 = db.Convert();
TDatime db(2022,01,02,00,00,00);
Int_t t2 = db.Convert();
Double_t yminPlot = 0., ymaxPlot= 200;  // min max of plot
Double_t Tscale=ymaxPlot*0.025/10., Toffset = ymaxPlot*0.9; // just for temperature plot
//===========================================================================
Int_t n = 3000;
Double_t intL[n],Temp[n],dTemp[n], Temp_plot[n];
Double_t time[n], Ttime[n];
Double_t weekL[n], dweekL[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,23,59,59);
  time[iweek]=d2.Convert();
  intL[iweek] /= 1000000.;
  Temp_plot[iweek] = Temp[iweek]*Tscale+Toffset;
  iweek++;
}
fin.close();
cout<<"Finished reading days_luminosity_temperature.txt, no of week="<<iweek<<endl;
//
//-------------------Calculate depletion voltage prediction -----------------
//
Double_t kB = 8.617343e-5;
Double_t TR=293.15, T0=273.15; 
Double_t epsilon0 = 8.854e-14, epsilon_prime=11.6;
Double_t e = 1.6020e-19;
Double_t TR=293.15, T0=273.15; 
//-------parameters and their errors------------------------
Double_t Neff00 =  Vd0 * 2.* epsilon_prime*epsilon0/e/pow(d,2); 
Double_t dNeff00 = dVd0* 2.* epsilon_prime*epsilon0/e/pow(d,2);
cout<<"Neff00="<<Neff00<<endl;
Double_t C0coeff0 = 0.7, dC0coeff0 = 0.05 ;// just for trial
//Double_t c0 = 1.04e-13,  dc0 = 0.;
Double_t ccoeff0 = 0.075, dccoeff0 = 0.005; // just for trial
Double_t ga0 = 0.018,    dga0 = 0.0014;
Double_t taua0 = 2.290,  dtaua0 = 0.2;  // just for trial
Double_t Ea0 = 1.090,    dEa0 = 0.03;
Double_t gc0 = 0.017,    dgc0 = 0.0005;
Double_t gy0 = 0.059,    dgy0 = 0.001;
Double_t tauy0 = 480.,   dtauy0 = 10.;  // just for trial
Double_t Ey0 = 1.330,    dEy0 = 0.03;
//
//------------initialization------------------------------
//
Double_t Neff[14][n], dNy[14][n], C[14][n], total_fluence[14][n];
Double_t Vdc[n], Vdu[n], Vdd[n], SigmaNeff[n];
for (Int_t ierror=0; ierror<14; ierror++) {
  Neff[ierror][0] = Neff00;
  dNy[ierror][0] = 0.;
  C[ierror][0] = 0.; 
  total_fluence[ierror][0] = 0.;
}
// 
//---------------------go over all weeks-----------------------------
//
for (Int_t iw=1; iw<iweek; iw++) {
    Double_t days = (time[iw] - time[iw-1])/(24.*60.*60);
//
//---------start loop for error estimate (erro item = 0 to 13)-----------------
    Double_t isum[14], differenceError[14];
    Double_t quadSum = 0.;
    for (Int_t ierror=0; ierror<14; ierror++) {
//----set back to central values for each time----
      Double_t Neff0 = Neff00;
      Double_t C0coeff = C0coeff0;
      Double_t ccoeff = ccoeff0;
      Double_t ga = ga0;
      Double_t taua = taua0;
      Double_t Ea = Ea0;
      Double_t gc = gc0;
      Double_t gy = gy0;
      Double_t tauy = tauy0;
      Double_t Ey = Ey0;
      Double_t Tsensor = Temp[iw];
      Double_t luminosity = weekL[iw];
      Double_t NIEL = NIEL0; 
//----change corresponding item-----------------
      switch(ierror){  
            case 1:  Neff0      += dNeff00;     break;
            case 2:  C0coeff    += dC0coeff0;   break;
            case 3:  ccoeff     += dccoeff0;    break;
            case 4:  ga         += dga0;        break;
            case 5:  taua       += dtaua0;      break;
            case 6:  Ea         += dEa0;        break;
            case 7:  gc         += dgc0;        break;
            case 8:  gy         += dgy0;        break;
            case 9:  tauy       += dtauy0;      break;
            case 10: Ey         += dEy0;        break;
            case 11: Tsensor    += dTemp[iw];   break;
            case 12: luminosity += dweekL[iw];  break;
            case 13: NIEL       += dNIEL0;       break;
            default: break; // central value
      }
//------calculate each component of the Hamburg model ---------------
      Double_t NC0 = C0coeff * Neff0;
      total_fluence[ierror][iw] = total_fluence[ierror][iw-1] + luminosity * NIEL;  
      Double_t ThetaTtau = exp(Ea/kB*(1/TR-1/(T0+Tsensor))) * days / taua;
      Double_t c = ccoeff / NC0;
      Double_t A = Neff0 - NC0*(1-exp(-c*total_fluence[ierror][iw]));
      Double_t B = gc * total_fluence[ierror][iw];
      C[ierror][iw] = ga * luminosity * NIEL / ThetaTtau*(1-exp(-ThetaTtau))
               + C[ierror][iw-1]*exp(-ThetaTtau);
      Double_t Theta_y = exp( Ey / kB * (1./TR - 1./(T0 + Tsensor)));
      Double_t gyF=gy * (total_fluence[ierror][iw-1]
               + total_fluence[ierror][iw])/2.;
      dNy[ierror][iw] = dNy[ierror][iw-1] 
               + gyF * (1-1/(1+Theta_y * days / tauy ));
      Neff[ierror][iw] = A - B - C[ierror][iw] - dNy[ierror][iw];
      if ( iw==210 || iw==314) {
        Double_t deltaNeff = Neff00 - Neff[ierror][iw];
        cout<<"iw,ierror="<<iw<<" "<<ierror<<" days="<<days<<" T="<<Tsensor<<" Neff=" 
        <<Neff[ierror][iw] <<" dN="<<deltaNeff<<"("<<
        (Neff[ierror][iw]-Neff[0][iw])/(Neff00-Neff[0][iw])*100.<<" %)"<<endl;
      } 
//------calculate quadratice error for each week ----------------------- 
      if(isum[0] != 0.) quadSum += pow(Neff[ierror][iw]-Neff[0][iw],2);
    } 
//--------------------------------------
    SigmaNeff[iw] = sqrt(quadSum); 
    if(iw==210 || iw==314) cout<<"result at"<<iw<<" "
     <<Neff[0][iw]<<"+-"<<SigmaNeff[iw]<<"("<<SigmaNeff[iw]/(Neff00-Neff[0][iw])*100.
     <<" %)"<<endl;  
    Vdc[iw] = e * d * d * fabs(Neff[0][iw])/(2*epsilon0*epsilon_prime);
    Double_t V1 = e * d * d * fabs(Neff[0][iw]+SigmaNeff[iw])/(2*epsilon0*epsilon_prime);     
    Double_t V2 = e * d * d * fabs(Neff[0][iw]-SigmaNeff[iw])/(2*epsilon0*epsilon_prime); 
    if( V1 > V2) { Vdu[iw] = V1; Vdd[iw] = V2; }
    else { Vdu[iw] = V2; Vdd[iw] = V1; }
    if((Neff[0][iw]+SigmaNeff[iw])* (Neff[0][iw]-SigmaNeff[iw])<0) Vdd[iw]=0.;
}

//
//-----------------add one more null data for color-area plotting--------------
//
iweek++;
TDatime dextra(2022,12,31,00,00,00);
time[iweek-1]=dextra.Convert();
Vdu[iweek-1] = 0.;
Vdd[iweek-1] = 0.;
//
Vdc[0] = Vdc[1]; Vdu[0] = 0.; Vdd[0] = 0.;

//-----------------------------------------------------
//
  hc = new TGraph(iweek,time,Vdc);
  hu = new TGraph(iweek,time,Vdu);
  hd = new TGraph(iweek,time,Vdd);
  hu->SetMaximum(ymaxPlot);
  hu->SetMinimum(yminPlot);
  hu->SetFillColor(kRed-10);
  hu->GetYaxis()->SetTitle("V_{d} [V]");
  hu->GetXaxis()->SetTitle("DATE");
  hu->GetXaxis()->SetTitleOffset(1.4);
  hu->GetXaxis()->SetLabelOffset(0.03);
  hu->GetXaxis()->SetLabelSize(0.04);
  hu->GetYaxis()->SetTitleSize(0.05);
  hu->GetYaxis()->SetTitleOffset(1.2);
  hu->GetYaxis()->SetLabelOffset(0.01);
  hu->GetXaxis()->SetTimeDisplay(1);
  hu->GetXaxis()->SetTimeFormat("#splitline{  %Y}{%b %d}");
  hu->GetXaxis()->SetLimits(t1,t2);
  hc->SetLineColor(2);
  hu->SetLineColor(1);
  hd->SetLineColor(4);
  hu->GetXaxis()->SetNdivisions(207);
  hu->GetYaxis()->SetNdivisions(505);
  hu->Draw("ALF");
  hd->SetFillColor(0);
  hd->Draw("SAME LF"); 
  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");
//-------------------------plot temperature--------------------------
Temp_plot[0] = Temp_plot[1];
gT = new TGraph(iweek-1,time,Temp_plot);
gT->SetLineColor(1);gT->Draw("SAME L");
//-------------draw an axis on the right side for temperature--------
TGaxis *axis = new TGaxis(t2,ymaxPlot*0.8,t2, ymaxPlot,-40.,40.,202,"+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();
//--- write tick mark on bottom X axis-------------------
for (Int_t year=2011; year<2022; year++){
  TDatime db(year,01,01,00,00,00); Int_t t4 = db.Convert();
  tick = new TLine(t4, 0, t4, ymaxPlot*0.03);tick->Draw();
}
//draw an axis on the left side
TGaxis *axis = new TGaxis(t1,yminPlot,t1, ymaxPlot,yminPlot,ymaxPlot,505,"L");
axis->SetLabelSize(0.0);
axis->Draw();
//-------draw legend---------------
leg = new TLegend(0.20 ,0.80 ,0.5,0.92);
leg->AddEntry(gL,"Integrated Luminosity (fb^{-1})","l");
leg->SetFillColor(10); leg->SetTextSize(0.04);leg->SetBorderSize(0);
//leg->Draw();  
//
c1->Print("Vd_sample_Plan1.pdf");
c1->Print("Vd_sample_Plan1.png");
}