#include <stdio.h>

/**************************************
compare true and reconstructed spectra
Arguments:
1.)rs_nom = the center of mass energy
**************************************/

void lumi_compare(TString fileName, TString opt="" , TString filtOpt="theta" ){

  //parse the options string to make a case switch on int
  //this is faster than doing it in the loop.

  Int_t fOpt[]={1,0,0,0,0,0,0};
  if (filtOpt.Contains("none")){
    fOpt[0]=0;
  }
  if (filtOpt.Contains("phi")){
    fOpt[1]=1;
  }
  if (filtOpt.Contains("0phot")){
    fOpt[2]=1;
  }
  if (filtOpt.Contains("1phot")){
    fOpt[3]=1;
  }
  if (filtOpt.Contains("2phot")){
    fOpt[4]=1;
  }
  if (filtOpt.Contains("bs")){
    fOpt[5]=1;
  }

  //define pi = pi
  const Double_t pi=TMath::Pi();
  //load the data
  TNtuple *lnt = makeBhwNtuple(fileName.Data());

  //make histos
  TH1D *rs_true = new TH1D("rs_true","rs_true",100,0.9,1.0);
  TH1D *rs_true1 = new TH1D("rs_true1","rs_true from e1 e2",100,0.9,1.0);
  TH1D *rs_rec0 = new TH1D("rs_rec0","rs_rec old",100,0,1.3);
  TH1D *rs_rec = new TH1D("rs_rec","rs_rec moenig",100,0.9,1.0);

  TH2D *rs_rec_rs_true = new TH2D("rs_rec_true","rs_rec (y axis) vs rs_true (x axis)",200,0.9,1.0,200,0.9,1.0);
  TH2D *rs_rec_th1 = new TH2D("rs_rec_th1","rs_rec (y axis) vs theta1 (x axis)",200,0.1,0.301,200,0.9,1.0);
  TH2D *rs_rec_th2 = new TH2D("rs_rec_th2","rs_rec (y axis) vs theta2 (x axis)",200,0.1,0.301,200,0.9,1.0);
  TH2D *rs_true_th1 = new TH2D("rs_true_th1","rs_true (y axis) vs theta1 (x axis)",200,0.1,0.301,200,0.9,1.0);
  TH2D *rs_true_th2 = new TH2D("rs_true_th2","rs_true (y axis) vs theta2 (x axis)",200,0.1,0.301,200,0.9,1.0);

  TH1D *beamstr1 =  new TH1D("beamstr1","beamsstrahlung 1",100,0.6,1.01);
  TH1D *beamstr2 =  new TH1D("beamstr2","beamsstrahlung 2",100,0.6,1.01);
  TH1D *beamspr1 =  new TH1D("beamspr1","beamspread 1",100,0.999,1.001);
  TH1D *beamspr2 =  new TH1D("beamspr2","beamspread 2",100,0.999,1.001);

  TH1D *theta1 = new TH1D("theta1","theta1",100,-pi,pi);
  TH1D *theta2 = new TH1D("theta2","theta2",100,-pi,pi);
  TH1D *theta_a = new TH1D("theta_a","acollinearity",100,-pi,pi);
  TH1D *theta_as= new TH1D("theta_as","acollinearity zoom",101,-0.1,0.1);
  TH1D *theta_as1= new TH1D("theta_as1","acollinearity zoom 1",101,-0.1,0.1);

  TH2D *theta12 = new TH2D("theta12","theta2 (y axis) vs theta1 (x axis)",200,0,pi,200,0,pi);

  TH1D *rs_recdjm = new TH1D("rs_recdjm","rs_rec miller",100,0.9,1.0);

  TH1D *phi1 = new TH1D("phi1","phi1",100,-pi,pi);
  TH1D *phi2 = new TH1D("phi2","phi2",100,-pi,pi);

  //the filestructure
  //"e:se:sp:be:bp:isr:tot:px1:py1:pz1:e1:px2:py2:pz2:e2:ngzpos:ngzneg"
  Float_t e1,e2,px1,py1,pz1,px2,py2,pz2,be,bp,se,sp,tot;
  Float_t nzpos,nzneg;
  lnt->SetBranchAddress("tot",&tot);
  lnt->SetBranchAddress("be",&be);
  lnt->SetBranchAddress("bp",&bp);
  lnt->SetBranchAddress("se",&se);
  lnt->SetBranchAddress("sp",&sp);
  lnt->SetBranchAddress("e1",&e1);
  lnt->SetBranchAddress("e2",&e2);
  lnt->SetBranchAddress("px1",&px1);
  lnt->SetBranchAddress("py1",&py1);
  lnt->SetBranchAddress("pz1",&pz1);
  lnt->SetBranchAddress("px2",&px2);
  lnt->SetBranchAddress("py2",&py2);
  lnt->SetBranchAddress("pz2",&pz2);
  lnt->SetBranchAddress("nzneg",&nzneg);
  lnt->SetBranchAddress("nzpos",&nzpos);


  //main loop
  cout << "filling histograms" << endl;
  Int_t j=0;
  Int_t filtered=0;
  for(Int_t i=0;i<lnt->GetEntries();i++) {

    lnt->GetEntry(i);

    TLorentzVector p1 = TLorentzVector(px1,py1,pz1,e1);
    TLorentzVector p2 = TLorentzVector(px2,py2,pz2,e2);


    //Moenig reconstructed rs
    Double_t rsrec = 0;
    if(p1->Theta()+p2->Theta() >= pi)
      rsrec = sqrt((1/tan(p1->Theta()/2))*(1/tan(p2->Theta()/2)));
    else
      rsrec = sqrt((1/tan((pi-p1->Theta())/2))*(1/tan((pi-p2->Theta())/2)));

    //Miller reconstructed
    Double_t th1, th2, thav, theta_aco,ph1,ph2;
    th1=p1->Theta();
    th2=p2->Theta();
    thav=pi- th2 +th1;
    theta_aco=pi - th1 - th2;
    ph1=p1->Phi();
    ph2=p2->Phi();
    Double_t rsrecdjm=1-(theta_aco)/(2*sin(thav/2));


    //true rs, not sqrt(tot)!
    Double_t rstrue = tot;
    //without angles
    Double_t rsrec0 = 1/sqrt(tan(p1->Theta()/2 )*tan(p2->Theta()/2));
    //the approximate ecms value
    Double_t rstrue1 = 2*sqrt(e1*e2)/350.0;

    bool go=1;
    //filter particles

    if (fOpt[0]){
      if (angleRange(th1,(pi-th2),0.099,0.301)){
   go=0;
      }
    }
    if (fOpt[1]){
      if (phiRange(ph1,ph2,0.00)){
   go=0;
      }
    }
    if (fOpt[2]){
      if (photFilter0(nzpos,nzneg)){
   go=0;
      }
    }
    if (fOpt[3]){
      if (photFilter1(nzpos,nzneg)){
   go=0;
      }
    }
    if (fOpt[4]){
      if (photFilter2(nzpos,nzneg)){
   go=0;
      }
    }
    if (fOpt[5]){
      if (beamsstrahlFilter(be,bp)){
   go=0;
      }
    }

    if (go==0) filtered++;

    else{
      //fill histograms
      rs_rec->Fill(rsrec);
      rs_true->Fill(rstrue);

      rs_rec0->Fill(rsrec0);
      rs_true1->Fill(rstrue1);
      rs_rec_rs_true->Fill(rstrue,rsrec);
      rs_recdjm->Fill(rsrecdjm);

      phi1->Fill(ph1);
      phi2->Fill(ph2);

      theta1->Fill(th1);
      theta2->Fill(th2);
      theta_a->Fill(theta_aco);
      theta_as->Fill(theta_aco);
      theta_as1->Fill(-theta_aco);

      theta12->Fill(th1,pi-th2);

      rs_rec_th1->Fill(th1,rsrec);
      rs_rec_th2->Fill(pi-th2,rsrec);
      rs_true_th1->Fill(th1,rstrue);
      rs_true_th2->Fill(pi-th2,rstrue);

      beamstr1->Fill(be);
      beamstr2->Fill(bp);
      beamspr1->Fill(se);
      beamspr2->Fill(sp);

      //just increment counter as well
      j++;
    }
  }

  cout <<"\nmaking plots "<< endl;
  cout <<"number of data points: " << j<< endl;
  cout <<"number of filtered points: " << filtered << endl;
  //normalise
  rs_true->Sumw2();
  rs_rec->Sumw2();
  rs_rec0->Sumw2();
  rs_true1->Sumw2();
  rs_recdjm->Sumw2();
  cout << "integral true, rec: "<<rs_true->Integral()<<" "<<rs_rec->Integral()<<endl;

  rs_true->Scale(1./rs_true->Integral());
  rs_rec->Scale(1./rs_rec->Integral());
  rs_rec0->Scale(1./rs_rec0->Integral());
  rs_true1->Scale(1./rs_true1->Integral());
  rs_recdjm->Scale(1./rs_recdjm->Integral());

  beamspr1->Scale(1./beamspr1->Integral());
  beamspr2->Scale(1./beamspr2->Integral());
  beamstr1->Scale(1./beamstr1->Integral());
  beamstr2->Scale(1./beamstr2->Integral());

  if(opt.Contains("fit")){
    cout << "\nfitting==============================================" <<endl;

    //fit to our user-defined function
    Double_t d_rstrue=lumiFit(rs_true,"fit",3);
    Double_t d_rsrec=lumiFit(rs_rec,"fit1",2);
    Double_t d_rsrecdjm=lumiFit(rs_recdjm,"fit2",4);

    cout << "========================================finito fitting" <<endl;


  }


  //make a nice filename to save the ps

  TString fName=fileName.ReplaceAll("/home/fz/gp/runucl/","").ReplaceAll("/fort.51","").ReplaceAll("-","").ReplaceAll("/","_");

  if(filtered!=0)fName=fName.Append("_filtered_").Append(filtOpt).Append("_").Append(numFormat(j-filtered));
  cout << "name: " << fName << endl;


  /////////////////////////////////
  // plot things

  //old things
  if (opt.Contains("old")){

    rs_true->SetLineColor(3);
    rs_true1->SetLineColor(5);
    rs_rec->SetLineColor(2);

    TCanvas *can3 = new TCanvas("c3","c3",600,400);
    //  can->Divide(1,1);
    can3->cd(0);
    gPad->SetLogy();
    rs_rec0.Draw("hist");
    rs_true.Draw("histsame");
    rs_rec.Draw("histsame");

  }
  //2-d plot
  if (opt.Contains("2D")){
    gStyle->SetOptStat(0);

    TCanvas *can2 = new TCanvas("c2","c2",500,450);
    gPad->SetLogz();
    rs_rec_rs_true->GetXaxis()->SetTitle("true");
    rs_rec_rs_true->GetYaxis()->SetTitle("reconstructed");
    rs_rec_rs_true.Draw("colorz");//colorz or surf3
    gStyle->SetPalette(1);
  }
  if (opt.Contains("angreco")){
    gStyle->SetOptStat(0);

    TCanvas *can7 = new TCanvas("c7","c7",700,550);
    can7->Divide(2,2);
    can7->cd(1);
    gPad->SetLogz();
    rs_rec_th1->GetXaxis()->SetTitle("#theta_{1}");
    rs_rec_th1->GetYaxis()->SetTitle("rs rec");
    rs_rec_th1.Draw("colorz");//colorz or surf3
    gStyle->SetPalette(1);
    can7->cd(2);
    gPad->SetLogz();
    rs_rec_th2->GetXaxis()->SetTitle("#theta_{2}");
    rs_rec_th2->GetYaxis()->SetTitle("rs rec");
    rs_rec_th2.Draw("colorz");//colorz or surf3
    gStyle->SetPalette(1);
    can7->cd(3);
    gPad->SetLogz();
    rs_true_th1->GetXaxis()->SetTitle("#theta_{1}");
    rs_true_th1->GetYaxis()->SetTitle("rs true");
    rs_true_th1.Draw("colorz");//colorz or surf3
    gStyle->SetPalette(1);
    can7->cd(4);
    gPad->SetLogz();
    rs_true_th2->GetXaxis()->SetTitle("#theta_{2}");
    rs_true_th2->GetYaxis()->SetTitle("rs true");
    rs_true_th1.Draw("colorz");//colorz or surf3
    gStyle->SetPalette(1);

    can7->Print("../ps/"+fName+"angles-reco.eps");

  }
  //look at the angles
  if (opt.Contains("angles")){
    //scale acollinearity histogram
    theta_as->Scale(1./theta_as->Integral());
    theta_as1->Scale(1./theta_as1->Integral());
    //make pave for acollinearity distribution
    TPaveText *pa = new TPaveText(0.05,0.4,0.4,0.6,"br");

    TString vals1[10];
    vals1[0]=" mean: ";
    vals1[0]+=numFormat(histMoment(theta_as,1));
    vals1[1]="true rms: ";
    vals1[1]+=numFormat(theta_as.GetRMS());
    vals1[2]="true skew: ";
    vals1[2]+=numFormat(histMoment(theta_as,3));
    vals1[3]="true kurt: ";
    vals1[3]+=numFormat(histMoment(theta_as,4));
    //fill pave
    pa->AddText("Acollinearity stats");
    for (Int_t i=0;i<4;i++){
      pa->AddText(vals1[i]);
    }

    gStyle->SetOptStat(0);
    gStyle->SetOptTitle(1);
    //make canvas
    TCanvas *can6 = new TCanvas("c6","c6",800,650);
    can6.Divide(2,2);

    can6.cd(1);
    theta12->GetXaxis()->SetTitle("#theta_{1}[radians]");
    theta12->GetYaxis()->SetTitle("#theta_{2}[radians]");
    theta12->Draw("");

    can6.cd(2);
    theta_as->GetXaxis()->SetTitle("#theta_{aco}[radians]");
    gPad->SetLogy();
    theta_as->Draw("hist");
    theta_as1->SetLineColor(3);
    theta_as1->Draw("histsame");

    can6.cd(3);
    pa.Draw();

    can6.cd(4);
    //gPad->SetLogy();
    TLegend *phileg = new TLegend(0,-1.2,3,-0.5,"#phi_{1} and #phi_{2} distributions","");
    phileg.AddEntry(phi1,"#phi_{1}");
    phileg.AddEntry(phi2,"#phi_{2}");
    phileg.Draw();
    phi1.Draw();
    phi2.SetLineStyle(2);
    phi2.SetLineColor(2);
    phi2.Draw("same");

    //print as eps
    can6->Print("../ps/"+fName+"angles.eps");
  }

  if(opt.Contains("spreads")){
    TCanvas *can10 = new TCanvas("c10","c10",500,650);
    can10->Divide(1,2);
    beamstr1->SetXTitle("beamsstrahlung ");
    beamspr1->SetXTitle("beamspread ");
    beamspr2->SetLineColor(3);
    beamstr2->SetLineColor(3);
    beamspr2->SetLineStyle(3);
    beamstr2->SetLineStyle(3);
    can10->cd(1);
    gPad->SetLogy();
    beamstr1->Draw();
    beamstr2->Draw("same");
    can10->cd(2);
    beamspr1->Draw();
    beamspr2->Draw("same");
  }


  //plot the miller recon on its own
  if (opt.Contains("djm")){
    TCanvas *can5 = new TCanvas("c5","c5",500,650);
    gStyle->SetOptStat(111111);
    can5->Divide(1,2);
    can5->cd(1);
    gPad->SetLogy();
    theta1->SetLineColor(2);
    theta2->SetLineColor(4);
    theta_a->Draw("hist");
    theta1->Draw("samehist");
    theta2->Draw("samehist");
    can5->cd(2);
    gPad->SetLogy();
    theta_as->Draw();
  }


  //output stats
  cout <<"mean: "<< TString(numFormat(1e-8)) << " " << histMoment(rs_true,4)<<" " << rs_true.GetMean()<<endl;

  //text and stats output
  //xleft,bottom,xright,ytop
  TPaveText *p = new TPaveText(0.1,0,0.9,1,"br");
  TPaveText *pt = new TPaveText(0.05,0.1,0.25,0.7,"br");
  TPaveText *pttl = new TPaveText(0.05,0.8,0.6,0.95,"br");
  //p->AddText(fileName.Data());

  //pttl->AddText("Using fitting function #sqrt{s}/#sqrt{s_{nom}} = ax^{b}(1-x)^{c} + d#delta(1-x)");
  pttl->AddText("Using fitting function #sqrt{s}/#sqrt{s_{nom}} = ax^{b}(1-x)^{c}  :: last bin from integral");

  TString vals[10];
  vals[0]="true mean: ";
  vals[0]+=numFormat(histMoment(rs_true,1));
  vals[1]="rec mean: ";
  vals[1]+=numFormat(histMoment(rs_rec,1));
  vals[2]="true rms: ";
  vals[2]+=numFormat(rs_true.GetRMS());
  vals[3]="rec rms: ";
  vals[3]+=numFormat(rs_rec.GetRMS());
  vals[4]="true skew: ";
  vals[4]+=numFormat(histMoment(rs_true,3));
  vals[5]="rec skew: ";
  vals[5]+=numFormat(histMoment(rs_rec,3));
  vals[6]="true kurt: ";
  vals[6]+=numFormat(histMoment(rs_true,4));
  vals[7]="rec kurt: ";
  vals[7]+=numFormat(histMoment(rs_rec,4));
  vals[8]="kept ";
  vals[8]+=numFormat(j);
  vals[8]+=" points";

  for (Int_t i=0;i<10;i++){
    //if (i>1) cout << vals[i] << endl;
    p->AddText(vals[i]);
    pt->AddText(vals[i]);
  }

  if (opt.Contains("fit")){
    //the three paves with fit data
    TPaveText *pst1 = new TPaveText(0.3,0.4,0.5,0.7,"br");
    TPaveText *pst2 = new TPaveText(0.5,0.4,0.7,0.7,"br");
    TPaveText *pst3 = new TPaveText(0.7,0.4,0.9,0.7,"br");

    pst1.SetTextColor(3);
    pst2.SetTextColor(2);
    pst3.SetTextColor(4);

    TString pst[5];
    pst[0]="a = ";
    pst[1]="b = ";
    pst[2]="c = ";
    pst[3]="last bin = ";

    pst1->AddText("True Lum Spectrum");
    pst2->AddText("Rec Moenig");
    pst3->AddText("Rec Miller");

    pst1->AddText("fit: #chi^{2} = "+numFormat(fit->GetChisquare()));
    pst2->AddText("fit: #chi^{2} = "+numFormat(fit1->GetChisquare()));
    pst3->AddText("fit: #chi^{2} = "+numFormat(fit2->GetChisquare()));

    for(Int_t j=0;j<4;j++){
      TString s1=pst[j];
      TString s2=pst[j];
      TString s3=pst[j];

      //this will not work with an if statement
      //weird cint bug!!
      switch (j){
      case 3:
   cout << "b1" << endl;
   s1+=numFormat(100*d_rstrue);
   s1+=" % of area";
   s2+=numFormat(100*d_rsrec);
   s2+=" % of area";
   s3+=numFormat(100*d_rsrecdjm);
   s3+=" % of area";
   break;
      default:
   s1+=numFormat( fit->GetParameter(j) );
   s2+=numFormat( fit1->GetParameter(j) );
   s3+=numFormat( fit2->GetParameter(j) );
   s1+=" #pm ";
   s2+=" #pm ";
   s3+=" #pm ";
   s1+=numFormat(fit->GetParError(j));
   s2+=numFormat(fit1->GetParError(j));
   s3+=numFormat(fit2->GetParError(j));
      }

      //cout << s1.Data()<<s2.Data()<<s3.Data() << endl;
      pst1.AddText(s1);
      pst2.AddText(s2);
      pst3.AddText(s3);

    }

  }

  if (!opt.Contains("noplot")){
    TString fName1=fName;
    TPaveLabel *lab = new TPaveLabel(0.89,-0.2,1,1.5,
                 fName1->ReplaceAll("_"," ").Data(),"br");
    p->SetFillColor(10);

    //formatting
    rs_true->SetLineColor(3);
    rs_rec->SetLineColor(2);
    rs_true1->SetLineColor(5);
    rs_recdjm->SetLineColor(4);

    //legend
    TLegend *leg = new TLegend(0.91,-1.2,0.95,-0.5,"true vs reconstructed","");
    leg.AddEntry(rs_true,"true spectrum");
    leg.AddEntry(rs_rec,"reconstructed (Moenig)");
    leg.AddEntry(rs_recdjm,"reconstructed (Miller)");

    //stats
    gStyle->SetOptStat(0);
    gStyle->SetOptTitle(0);

    //make canvas and draw
    TCanvas *can = new TCanvas("c0","c0",500,700);
    can->Divide(1,2);

    can->cd(1);
    gPad->SetLogy();
    rs_true.Draw("");
    rs_rec.Draw("samee");
    rs_recdjm.Draw("samee");
    //rs_true1.Draw("smoothsame");
    leg->Draw();
    //p->Draw();
    lab->Draw();

    can->cd(2);
    pt->Draw();
    pttl->Draw();
    if (opt.Contains("fit")){
      pst1.Draw();
      pst2.Draw();
      pst3.Draw();
    }
    can->Print("../ps/"+fName+".eps");
  }

  //invoke the top graph function
  if (opt.Contains("top")){
    cout << "\ndrawing top cross section" << endl;

    //invoke the function for all different histograms needed.
    const Int_t graphSize=150;
    Double_t xtrue[graphSize];
    Double_t ytrue[graphSize];
    Double_t y1true[graphSize];
    topGraphs(rs_true,xtrue,ytrue,y1true,graphSize);

    //Double_t xrec[graphSize];
    Double_t yrec[graphSize];
    Double_t y1rec[graphSize];
    topGraphs(rs_rec,xtrue,yrec,y1rec,graphSize);

    //Double_t xrecdjm[graphSize];
    Double_t yrecdjm[graphSize];
    Double_t y1recdjm[graphSize];
    topGraphs(rs_recdjm,xtrue,yrecdjm,y1recdjm,graphSize);

    TGraph *cleanSec=new TGraph(graphSize,xtrue,ytrue);
    TGraph *smearSec=new TGraph(graphSize,xtrue,y1true);
    TGraph *smearSec_rec=new TGraph(graphSize,xtrue,y1rec);
    TGraph *smearSec_recdjm=new TGraph(graphSize,xtrue,y1recdjm);

    //make canv BUSINESSas
    TCanvas *topcan=new TCanvas("topcan",
            "top true and observed x-section",600,400);

    //draw all the different graphs
    TPaveLabel *toplab = new TPaveLabel(345,1.5,352,1.7,
               fName1->ReplaceAll("_"," ").Data(),"br");


    TLegend *leg1=new TLegend(349.6,0.15,352.4,0.5,"","");
    cleanSec->SetLineColor(1);
    cleanSec->SetLineStyle(4);
    smearSec->SetLineColor(3);
    smearSec_rec->SetLineColor(2);
    smearSec_recdjm->SetLineColor(4);
    leg1->AddEntry(cleanSec,"\"real\" top cross section");
    leg1->AddEntry(smearSec, "degraded with guinea-pig Lumi Spectrum");
    leg1->AddEntry(smearSec_rec,"degraded with Moenig reconstructed Lumi Spectrum");
    leg1->AddEntry(smearSec_recdjm,
         "degraded with Frary-Miller reconstructed Lumi Spectrum");
    //draw 
    cleanSec->Draw("AC");
    smearSec->Draw("C");
    smearSec_rec->Draw("C");
    smearSec_recdjm->Draw("C");
    toplab->Draw();
    cleanSec->GetHistogram()->SetXTitle("#sqrt{s}");
    cleanSec->GetHistogram()->SetYTitle("#sigma");
    leg1->Draw();


    //print a ps file
    topcan->Print("../ps/"+fName+"topsection.eps");

  }

  return;
}