//   HiLumSimDataC.C
//   24 May 2013
//   Plot HiLum drift data and its C-History
//
//   Run this as follows:  root -l -L HiLumSimDataC.C++
//   
#include "Riostream.h"
#include <string>
#include <fstream>
#include <iomanip>
#include <cmath>
#include <TCanvas.h>
#include <TStyle.h>
#include <TH1D.h>
#include <TH2D.h>
#include <TGraph.h>
#include <TLatex.h>
#include <TFrame.h>
#include <TLine.h>
#include <TTimeStamp.h>
#include <TImage.h>
#include <TGaxis.h>
//#include "/home/rutherfo/main/lhc/atlas/lhcupgrade/Intas/protvino/Mar2013Data/HiLumUtils"
#include "HiLumUtils"

void HiLumSimDataC() {
  //int SimNo=39; // Simulation data number for run 1141
  //int SimNo=29; // Simulation data number for run 1092
  int SimNo=47; // Simulation data from other
  int SplNo= 1; // Spill number (For Run 1092 use SplNo=1, for Run 1141 use SplNo=2
  double Area=60.*60.*6; // Area of EMEC electrodes
  int ListNo;
  int RunNo;
  char asc[200];
  char DateTime[200];
  char parms[200];
  char fn[200]; // Simulation file name to read in
  char fnS[]="drifter-SimRuns.logbook";
  ifstream inS;
  inS.open(fnS); // Read in list of simulation file names
    for(int i=0;i<8;i++) {
      inS.getline(asc,200); // Skip over header lines
    }
    do {
      inS >> ListNo >> DateTime >> RunNo; inS.getline(parms,200);
    } while(SimNo!=ListNo);
  inS.close();
  sprintf(fn,"drifter-%s.dat",DateTime);
  printf("\n\nSimulated data file %s\n",fn);
  char GateCFileName[100];
  Cerenkov(SplNo,GateCFileName); // Read in Cerenkov data
  double xlo,xhi,ylo,yhi;
  xlo=0; xhi=39*165.; ylo=0; yhi=1000;
  //xlo=0; xhi=44*165.; ylo=0; yhi=1000;
  //xlo=0.0; xhi=4.0; ylo=0; yhi=1000;
  double eq=1.602E-10; // Charge on electron in nC
  double Jc,V0=1200.,a=2.0,kappa=1.51,eps=8.854E-15,muplus=0.1;
  Jc=4.*V0*V0*kappa*eps*muplus/pow(a,3)*1.E9; // Critical current in EMEC in nA/mm^2
  printf(" Jc= %10.3f nA/mm^2\n",Jc);

  int EvtNo;
  printf("For Run %4d enter Event Number: ",RunNo);
  //cout << "Enter Event Number: ";
  cin >> EvtNo;
  if(EvtNo<0) {cout << "Enter a positive integer Event Number: "<<endl; return;}
  if(EvtNo>=kevents) {cout << "This spill has only "<<kevents<<" events:  Try again!"
     <<endl; return;}

  const int nGADC=39;
  double tGADC[nGADC],cGADC[nGADC],cGADCmax=-1.E9,cGADCmin=+1.E9;
  //printf("  i    tGADC    cGADC\n");
  for(int i=0;i<nGADC;i++) {
    tGADC[i]=165.*double(i);
    cGADC[i]=0;
    if(i>0) cGADC[i]=double(Cbuf[EvtPnt[EvtNo]-1+i])-AvePedBeg;
    if(cGADC[i]>cGADCmax) cGADCmax=cGADC[i];
    if(cGADC[i]<cGADCmin) cGADCmin=cGADC[i];
    //printf("%3d %8.1f %8.1f\n",i,tGADC[i],cGADC[i]);
  }
  printf(" Max/Min Cherenkov: %10.2f %10.2f\n",cGADCmax,cGADCmin);
  cGADCmax*=1.1;
  cGADCmin*=(cGADCmin>=0.0) ? 0.9 : 1.1;  

  double tm,tm0,jh1,jh2,jh3,t0,t1;
  long it,it0;
  const int nIter=15000;
  double time1[nIter],Jh1[nIter],Jh2[nIter],Jh3[nIter];
  long itime1[nIter];
  int nSamp1=0;
  char dum1[30],dum2[30],dum3[30],dum4[30],dum5[30],dum6[30],dum7[30];
  double emm3ADC; // Ionization density per ADC count in units of e/mm^3/ADC count
  double ppADC;   // Protons per ADC count
  ifstream in;
  in.open(fn);
    for(int j=0;j<2;j++) in.getline(asc,200); // Skip first two header lines
    in >> dum1 >> dum2 >> dum3 >> dum4 >> dum5 >> ppADC >> dum6 >> dum7 >> emm3ADC;
          in.getline(asc,200);
    printf(" ppADC: %15.5f  emm3ADC: %15.5f\n",ppADC,emm3ADC);
    for(int j=0;j<3;j++) in.getline(asc,200); // Skip last two header lines
    for(Int_t j=0;j<=EvtNo;j++) {
      int i=0;
      in >> it0 >> tm0 >> jh1 >> jh2 >> jh3; in.getline(asc,200);
      t0=1E9*tm0; t1=t0-165.0;
      itime1[i]=it0; time1[i]=t0-t1; Jh1[i]=eq*jh1; Jh2[i]=eq*jh2; Jh3[i]=eq*jh3;
      i++; // printf("Here: j,i=%4d %4d\n",j,i);
      while(1) {
        in >> it >> tm >> jh1 >> jh2 >> jh3; in.getline(asc,200);
	if(in.eof()) {
	  if(i<2) {printf("***** No such event!\n"); return;}
	  break;
	}
	if((it!=(itime1[i-1]+1))||in.eof()) {break;}
	itime1[i]=it; time1[i]=1E9*tm-t1;
        Jh1[i]=eq*Area*jh1/1000.; Jh2[i]=eq*Area*jh2/1000.; Jh3[i]=eq*Area*jh3/1000.;
        // if(i<20) printf("%12ld %14.10f %9.4f %14.7f\n",itime1[i],tm,time1[i],Jh1[i]);
	i++;  if(i>=nIter) printf("Warning! i exceeds bounds\n");
      }
      nSamp1=i-1;
      // cout << endl;
    }
  double DT=tm/double(it)*1.E9;
  printf(" DT: %10.7f\n",DT);
  in.close();
  // Iteraction steps to back up from peak to get current before the "step"
  int ITback=int(30./DT);

  printf("Number of samples %5d\n",nSamp1);

  TCanvas* c1=new TCanvas("c1","HV Plateau",350,20,1000,470);
  //TCanvas* c1=new TCanvas("c1","HV Plateau",700,0,1000,700);
  c1->SetFillColor(10);   //white
  c1->Divide(1,2);
  c1->cd(1)->SetGrid();          //put grid on graph
  c1->cd(2)->SetGrid();

  // Find maximum current in event
  int Im=0; // RFbunch in this event at first isolated filled bunch
  double Stepm;
  yhi=-1000.; Stepm=0.0;
  for(int i=0;i<nSamp1;i++) {
    if(Jh1[i]>yhi) {
      yhi=Jh1[i];
      if(i>=4) {
        Stepm=yhi-Jh1[i-ITback];
        Im=i;
      }
    }
  }
  yhi*=1.1;
  for(int i=0;i<nGADC;i++) {cGADC[i]*=yhi/cGADCmax;}
  ylo=cGADCmin*yhi/cGADCmax;

  TH2D *hpx1=new TH2D("hpx1","",10,xlo,xhi,10,ylo,yhi);
  hpx1->SetStats(kFALSE);
  hpx1->GetXaxis()->SetTitle("                 Time, ns");
  hpx1->GetXaxis()->SetTitleSize(0.07);
  hpx1->GetXaxis()->SetTitleOffset(0.53);
  hpx1->GetXaxis()->CenterTitle();
  hpx1->GetXaxis()->SetLabelSize(0.07);
  hpx1->GetYaxis()->SetLabelSize(0.07);
  hpx1->GetYaxis()->SetNdivisions(1005);
  hpx1->GetYaxis()->SetTitle("Signal Current, #muA");
  hpx1->GetYaxis()->SetTitleSize(0.07);
  hpx1->GetYaxis()->SetTitleOffset(0.60);
  hpx1->GetYaxis()->CenterTitle();
  hpx1->SetFillColor(10);  //white

  double tlo=-20.0, thi=1.0; // Time in ms
  double RFcycleT=165.0E-9; // RF period
  double Fwrd=+thi/1000./RFcycleT; // Number of RF bunches after Event
  double Back=-tlo/1000./RFcycleT; // Number of RF bunches before Event
  int nBack=int(Back);
  int nFwrd=int(Fwrd);
  // No. of deca-turns (bin of 10 turns) in history
  int ntBins=int((Fwrd+Back)/double(nBnTrn)/10.);
  //TStyle *tst=new TStyle();
  //tst->SetHistFillStyle(4000);
  TH1D *hpx2=new TH1D("hpx2","",ntBins,tlo,thi);
  hpx2->SetStats(kFALSE);
  hpx2->GetXaxis()->SetTitle("Time, ms                                                         ");
  hpx2->GetXaxis()->SetTitleSize(0.07);
  hpx2->GetXaxis()->SetTitleOffset(0.53);
  hpx2->GetXaxis()->CenterTitle();
  hpx2->GetXaxis()->SetLabelSize(0.07);
  hpx2->GetYaxis()->SetLabelSize(0.07);
  hpx2->GetYaxis()->SetNdivisions(1005);
  hpx2->GetYaxis()->SetTitle("Filled bunches - 10 turn Average");
  hpx2->GetYaxis()->SetTitleSize(0.07);
  hpx2->GetYaxis()->SetTitleOffset(0.60);
  hpx2->GetYaxis()->CenterTitle();
  hpx2->SetFillColor(10);  //white
  hpx2->SetLineColor(2);   //red

  double aNorm=10.*double(mModFlg); // Number of filled RF bunches per bin of 10 turns
  double CAve=0.0; // Average Cherenkov pulse height
  double f1,f2;
  int ij;
  for(int i=0;i<nFwrd+nBack;i++) {
    ij=JEnd[EvtNo]+(nFwrd-i);
    f1=float(nFwrd-i)*1000.*RFcycleT; // Time in ms relative to event
    f2=0.0;
    if(ij>=0) f2=(double(Cbuf[ij])-AvePedBeg)/aNorm; 
    hpx2->Fill(f1,f2);
    CAve+=f2;
  }
  // Average sum of filled bunches in a deca-turn, averaged over deca-turns in history
  CAve/=double(ntBins);
  printf("CAve=%10.3f\n",CAve);
  double sCAve=CAve*double(mModFlg); // GADC sum in a turn, averaged over turns in history
  double r=sCAve*emm3ADC*a/(RFcycleT*double(nBnTrn))*eq/Jc;

  double tEvt1,tEvt2; // Time of events (relative to trigger) at and just before trigger
  tEvt1=(EvtPnt[EvtNo]-JEnd[EvtNo])*165.E-9*1000.;
  tEvt2=(EvtNo>0) ? (EvtPnt[EvtNo-1]-JEnd[EvtNo])*165.E-9*1000.: -40.;
  double Hbin, Hmax=-1.E9, Hmin=+1.E9;
  for(int i=1;i<=ntBins;i++) {
    Hbin=hpx2->GetBinContent(i);
    if(Hbin>Hmax) Hmax=Hbin;
    if(Hbin<Hmin) Hmin=Hbin;
  }
  if(Hmin>0.0) Hmin=0.0; // Don't suppress zero
  double Hdiff=Hmax-Hmin;
  Hmax+=0.02*Hdiff; Hmin-=0.02*Hdiff;
  hpx2->SetMaximum(Hmax);
  hpx2->SetMinimum(Hmin);

  TGraph *gr1=new TGraph(nSamp1,time1,Jh1);
  gr1->SetLineColor(2);
  gr1->SetLineWidth(2);  // Linewidth of error bars
  gr1->SetLineStyle(1);  // Linestyle solid

  TGraph *gr2=new TGraph(nSamp1,time1,Jh2);
  gr2->SetLineColor(3);
  gr2->SetLineWidth(2);  // Linewidth of error bars
  gr2->SetLineStyle(1);  // Linestyle solid

  TGraph *gr3=new TGraph(nSamp1,time1,Jh3);
  gr3->SetLineColor(4);
  gr3->SetLineWidth(2);  // Linewidth of error bars
  gr3->SetLineStyle(1);  // Linestyle solid

  c1->cd(1); hpx1->Draw(); gr1->Draw("L"); gr2->Draw("L"); gr3->Draw("L");
  c1->cd(2); hpx2->Draw();
  TLine *LAve = new TLine(tlo,CAve,thi,CAve);
  LAve->Draw();

// Draw a histogram
  c1->cd(1);
  double tblo,tbhi,yblo=0.0,ybhi;
  TLine *LVert = new TLine();
  TLine *LHorz = new TLine();
  for(int i=0;i<nGADC;i++) {
    tblo=tGADC[i]; tbhi=tblo+165.;
    ybhi=cGADC[i];
    LVert->DrawLine(tblo,yblo,tblo,ybhi);
    LHorz->DrawLine(tblo,ybhi,tbhi,ybhi);
    yblo=ybhi;
  }

// y-axis at right of plot 
  c1->cd(1);
  TGaxis *Haxis = new TGaxis(xhi,ylo,xhi,yhi,cGADCmin,cGADCmax,1005,"+"); 
  Haxis->SetLabelSize(0.07); 
  Haxis->SetLabelOffset(0.04);
  Haxis->SetLabelFont(42);
  Haxis->SetTitle("GADC counts above pedestal");
  Haxis->SetTitleSize(0.07);
  Haxis->SetTitleOffset(0.50);
  Haxis->SetTitleFont(42);
  Haxis->CenterTitle();
  Haxis->Draw(); 

  int itChi=TnPntr[EvtNo]+iModFlg[0];
  double tChi=165.0*(double(itChi)+0.5)-t1;
  if(tChi>30.0*165.0) {tChi-=30.0*165.0; itChi-=30;}
  //printf(" itChi: %3d\n",itChi);
  double yChi=double(Cbuf[itChi]-AvePedBeg);
  //printf("ITIME =%10d, yChi =%8.1f\n",TnPntr[EvtNo]+iModFlg[0],yChi);
  // Find maximum current near the first isolated filled bunch
  int In=0;
  double ynhi=-1000., Stepn=0.0;
  for(int i=0;i<nIter;i++) {if(fabs(time1[i]-tChi)<DT) Im=i;}
  int iDT=int(165.0/DT);
  int itSt=Im-iDT,itEn=Im+iDT;
  double tSt=time1[itSt],tEn=time1[itEn];
  printf(" tSt,tEn %15.10f %15.10f\n",time1[itSt],time1[itEn]);
  for(int i=itSt;i<itEn;i++) {
    if(Jh1[i]>ynhi) {
      ynhi=Jh1[i];
      if(i>=4) {
        Stepn=ynhi-Jh1[i-ITback];
        In=i;
      }
    }
  }
  printf(" itStart,itEnd,In: %10d %10d %10d\n",itSt,itEn,In);
  TLine *fIso = new TLine();
  fIso->SetLineWidth(2);
  fIso->SetLineColor(3);
  fIso->DrawLine(tChi,0.0,tChi,0.1*yhi);
  fIso->DrawLine(tSt,ylo,tEn,ylo);
  double phC=Stepn/yChi; // Pulse height per GADC
  double php=phC/ppADC; // Pulse height per proton
  char Start1[200],Start2[200];
  sprintf(Start1,"Iter: %12ld, Time %13.4f #mus   %s",it0,1E6*tm0,fn);
  sprintf(Start2,
    "EvntNo:%4d,  Peak:%9.3f,  Stepn:%9.3f  Ceren:%7.1f  Ratio:%9.4f  php3:%9.4f  Hstry1:%8.1f",
    EvtNo,ynhi,Stepn,yChi,phC,php*1.E3,CHist1[EvtNo]);
  double Dy=yhi-ylo;
  double yl1=yhi+0.07*Dy;
  double yl2=yhi+0.01*Dy;
  TLatex *TI = new TLatex();
  TI->SetTextAlign(11); 
  TI->SetTextFont(42); 
  TI->SetTextSize(0.055); 
  TI->SetLineWidth(1); 
  c1->cd(1);
  TI->DrawLatex(xlo,yl1,Start1); 
  TI->DrawLatex(xlo,yl2,Start2);
  c1->cd(2);
  TLine *Levt = new TLine();
  Levt->SetLineWidth(2);
  Levt->SetLineColor(3);
  Levt->DrawLine(tEvt1,0.0,tEvt1,0.1*Hmax);
  Levt->DrawLine(tEvt2,0.0,tEvt2,0.1*Hmax);
  char rFactor[100];
  sprintf(rFactor,"r = %7.2f for EMEC;   SimNo = %3d;   parms =  %s",r,SimNo,parms);
  printf("Max/Min GADC History: %8.1f %8.1f\n",Hmax,Hmin);
  yl1=Hmax+0.06*(Hmax-Hmin);
  xlo=tlo;
  TI->SetTextAlign(13);
  TI->DrawLatex(xlo,yl1,rFactor); 
  TI->DrawLatex(xlo+15,yl1,"");
  TI->DrawLatex(xlo+25,yl1,"");
  char c1Name[100];
  sprintf(c1Name,"SimR%04dS%03dE%03d.pdf",RunNo,SplNo,EvtNo);
  //c1->Print(c1Name,"pdf");
  sprintf(c1Name,"SimR%04dS%03dE%03d.gif",RunNo,SplNo,EvtNo);
  c1->Print(c1Name,"gif");
  // Record data at end of table
  char fileout[]="HiLumSimDataC.dat";
  FILE *fout=fopen(fileout,"a");
  fprintf(fout,"%10ld%13.4f%6d%7.1f%7.1f%8.3f%10.1f\n",
      it0,1E6*tm0,EvtNo,Stepm,yChi,phC,CHist1[EvtNo]);
  fclose(fout);

} // End of Root macro