//======
// ROOT
//======
#include <TString.h>
#include <TMath.h>
#include <TFile.h>
#include <TTree.h>
#include <TCanvas.h>
#include <TH1.h>
#include <TVector3.h>
#include <TStyle.h>
#include <TF1.h>
#include <TGaxis.h>
#include <TFrame.h>


//=========
// STD C++
//=========
#include <iostream>

using namespace std;



void recoCherenkovAngle( TString beamtestFilename = "",  TString kBarFilename = "", Int_t bining = 100 )
{
//   TString beamtestFilename = "";
//   TString kBarFilename     = "";

  if( beamtestFilename == "" || kBarFilename == "" )
  {
    cout << "Usage: recoCherenkovAngleTrue( beamtest-filename, kBar-filename, bining )" << endl;
    cout << "  screenPix.cc        => beamtest-File" << endl;
    cout << "  run_KBarAnalysis.cc => kBar-File    " << endl;
    return;
  }

  cout << "Check yourself if beamtest-file and kBar-file were generated by the same parameter!" << endl;


  Double_t rad2degree = 180./TMath::Pi();
  Double_t degree2rad = TMath::Pi()/180.;
  Double_t deg2mrad = degree2rad*1000;
  Double_t mrad2deg = rad2degree/1000;

  // plot range
  Double_t xmin_c = 43. *deg2mrad; // number in degree
  Double_t xmax_c = 51. *deg2mrad;



//==============================================================================
// Access to the input ROOT-files
//==============================================================================
  TFile *beamtestFile = new TFile( beamtestFilename );
  TTree *pixel_beamtest = (TTree*) beamtestFile->Get("pixel");
  TTree *info = (TTree*) beamtestFile->Get("info");


  Double_t parDirX, parDirY, parDirZ;
  Int_t parDir_size;
  Double_t thetaC;
  Int_t x_bins, y_bins;

  info->SetBranchAddress( "parDirX", &parDirX );
  info->SetBranchAddress( "parDirY", &parDirY );
  info->SetBranchAddress( "parDirZ", &parDirZ );
  info->SetBranchAddress( "parDir_size", &parDir_size );
  info->SetBranchAddress( "thetaC" , &thetaC );
  info->SetBranchAddress( "x_bins" , &x_bins );
  info->SetBranchAddress( "y_bins" , &y_bins );
  info->GetEntry( 0 );


  Bool_t opt_angleAcceptance = false;

  if( parDirX == -666 )
  {
    opt_angleAcceptance = true;
    cout << "*** included particles with different incidence angles ***" << endl;
  }


  Int_t hitsPerPixel_helper;
  if( opt_angleAcceptance )
    hitsPerPixel_helper = parDir_size;
  else
    hitsPerPixel_helper = 1;

  const Int_t hitsPerPixel = hitsPerPixel_helper;


  TString parDir_size_str;
  parDir_size_str += hitsPerPixel;
  parDir_size_str.Remove( TString::kLeading, ' ' );
  TString parDirX_str = "parDirX[" + parDir_size_str + "]";
  TString parDirY_str = "parDirY[" + parDir_size_str + "]";
  TString parDirZ_str = "parDirZ[" + parDir_size_str + "]";

  Double_t phot_parDirX[hitsPerPixel];
  Double_t phot_parDirY[hitsPerPixel];
  Double_t phot_parDirZ[hitsPerPixel];

  for( int i = 0; i < hitsPerPixel; i++)
  {
    phot_parDirX[i] = 0;
    phot_parDirY[i] = 0;
    phot_parDirZ[i] = 0;
  }


  Int_t freq_beamtest;

  pixel_beamtest->SetBranchAddress( "freq", &freq_beamtest );
  if( opt_angleAcceptance )
  {
    pixel_beamtest->SetBranchAddress( parDirX_str, phot_parDirX );
    pixel_beamtest->SetBranchAddress( parDirY_str, phot_parDirY );
    pixel_beamtest->SetBranchAddress( parDirZ_str, phot_parDirZ );
  }
  pixel_beamtest->GetEntry(0);

  Int_t nPixel_beamtest = pixel_beamtest->GetEntries();


  TFile *kBarFile = new TFile( kBarFilename );
  TTree *pixel_kBar = (TTree*) kBarFile->Get("pixel");

  Double_t kBarX[6], kBarY[6], kBarZ[6];
  Double_t kBarXerr[6], kBarYerr[6], kBarZerr[6];
  Int_t freq[6];

  pixel_kBar->SetBranchAddress( "kBarX[6]"   , kBarX );
  pixel_kBar->SetBranchAddress( "kBarY[6]"   , kBarY );
  pixel_kBar->SetBranchAddress( "kBarZ[6]"   , kBarZ );
  pixel_kBar->SetBranchAddress( "kBarXerr[6]", kBarXerr );
  pixel_kBar->SetBranchAddress( "kBarYerr[6]", kBarYerr );
  pixel_kBar->SetBranchAddress( "kBarZerr[6]", kBarZerr );
  pixel_kBar->SetBranchAddress( "freq[6]"    , freq );

  Int_t nPixel_kBar = pixel_kBar->GetEntries();


  if( nPixel_kBar > 1000 && opt_angleAcceptance )
  {
    cout << "*******************************************" << endl;
    cout << "It's recommended to compile this macro.    " << endl;
    cout << "Compile procedure:                         " << endl;
    cout << "root[] .L recoCherenkovAngle.cc++          " << endl;
    cout << "*******************************************" << endl;
  }



//==============================================================================
// Event loop
//==============================================================================
  TCanvas *canvas = new TCanvas( "canvas", "", 1 );
  canvas->Draw();
  canvas->SetTopMargin( 0.14 ); // for 2 x-axis

  TH1F *cherenkov = new TH1F( "cherenkov_angle", "", bining, xmin_c, xmax_c);

  cherenkov->GetXaxis()->SetTitle( "#Theta_{c} [mrad]" );
  cherenkov->GetXaxis()->CenterTitle();
  cherenkov->SetMinimum(0);


  if( nPixel_beamtest != nPixel_kBar )
  {
    cout << "Different number of pixel in beamtest data and kBar tree! (" << nPixel_beamtest << ", " << nPixel_kBar << ")" << endl;
    return;
  }


  Int_t pxY = 0; // row
  Int_t pxX = 0; // col


  for( int i = 0; i < nPixel_kBar; i++ )
  {
    pixel_kBar->GetEntry( i );
    pixel_beamtest->GetEntry( i );

    if( i%1000 == 0 && opt_angleAcceptance )
      cout << "Pixel: " << i+1 << " / " << nPixel_kBar << endl;


    if( i%y_bins == 0 && i != 0 ) // next column
    {
      pxX++;
      pxY = 0;
    }

//     cout << "Pixel: " << i+1 << " / " << nPixel_kBar << "  col: " << pxX << "  row: " << pxY << endl;


    TVector3 kBar[6];

    Int_t loopForDiffParctiles = hitsPerPixel;

    for( int l =0; l < loopForDiffParctiles; l++)
    {
      TVector3 parDir;

      if( opt_angleAcceptance )
      {
        parDir.SetX( phot_parDirX[l] );
        parDir.SetY( phot_parDirY[l] );
        parDir.SetZ( phot_parDirZ[l] );
      }
      else
      {
        parDir.SetX(parDirX);
        parDir.SetY(parDirY);
        parDir.SetZ(parDirZ);
      }

//       cout << "Particle: " << i+1 << " / " << loopForDiffParctiles
//           << "  parDir: (" + parDir.X() << ", " << parDir.Y() << ", " << parDir.Z() << ")" << endl;


      for( int j = 1; j < 6; j++ ) // 0: all, 1: direct, 2: left, 3: right, 4: up, 5: down
      {
        Double_t angleC[8]; // up, down, left, right, forward, backward; total 2*2*2 possibilities in bar

        kBar[j].SetX( kBarX[j] );
        kBar[j].SetY( kBarY[j] );
        kBar[j].SetZ( kBarZ[j] );

        Double_t mag = parDir.Mag() * kBar[j].Mag();

        if( mag == 0 && !opt_angleAcceptance )
        {
          cout << "Zero vector!" << endl;
          return;
        }


        angleC[0] = TMath::ACos( kBar[j].Dot(parDir) / mag );

        kBar[j].SetX( -kBarX[j] );
        kBar[j].SetY( kBarY[j] );
        kBar[j].SetZ( kBarZ[j] );
        angleC[1] = TMath::ACos( kBar[j].Dot(parDir) / mag );

        kBar[j].SetX( kBarX[j] );
        kBar[j].SetY( -kBarY[j] );
        kBar[j].SetZ( kBarZ[j] );
        angleC[2] = TMath::ACos( kBar[j].Dot(parDir) / mag );

        kBar[j].SetX( kBarX[j] );
        kBar[j].SetY( kBarY[j] );
        kBar[j].SetZ( -kBarZ[j] );
        angleC[3] = TMath::ACos( kBar[j].Dot(parDir) / mag );

        kBar[j].SetX( -kBarX[j] );
        kBar[j].SetY( -kBarY[j] );
        kBar[j].SetZ( kBarZ[j] );
        angleC[4] = TMath::ACos( kBar[j].Dot(parDir) / mag );

        kBar[j].SetX( -kBarX[j] );
        kBar[j].SetY( kBarY[j] );
        kBar[j].SetZ( -kBarZ[j] );
        angleC[5] = TMath::ACos( kBar[j].Dot(parDir) / mag );

        kBar[j].SetX( kBarX[j] );
        kBar[j].SetY( -kBarY[j] );
        kBar[j].SetZ( -kBarZ[j] );
        angleC[6] = TMath::ACos( kBar[j].Dot(parDir) / mag );

        kBar[j].SetX( -kBarX[j] );
        kBar[j].SetY( -kBarY[j] );
        kBar[j].SetZ( -kBarZ[j] );
        angleC[7] = TMath::ACos( kBar[j].Dot(parDir) / mag );


//         angleC[2] = 0; // phi is zero; symmetry case
//         angleC[4] = 0;
//         angleC[6] = 0;
//         angleC[7] = 0;
        //
//         angleC[3] = 0; // no mirrored photons
//         angleC[5] = 0;
//         angleC[6] = 0;
//         angleC[7] = 0;


        for( int k = 0; k < 8; k++ )
        {
          if( TMath::Abs( kBar[j].X() ) == 666 || TMath::Abs( kBar[j].Y() ) == 666 || TMath::Abs( kBar[j].Z() ) == 666 )
            angleC[k] = 0;
        }


        for( int k = 0; k < 8; k++ )
        {
          if( j > 0 && angleC[k] != 0 )
          {
            if( loopForDiffParctiles > 1 )
              cherenkov->Fill( angleC[k]*1000, 1 );
            else
              cherenkov->Fill( angleC[k]*1000, freq_beamtest );
          }
        } // solution loop
      } // add. ambigiuity loop
    } // particle loop

    pxY++;
  } // pixel loop



//==============================================================================
// Fit
//==============================================================================
  gStyle->SetStatX(0.97);
  gStyle->SetStatY(0.81);
  gStyle->SetOptStat("mr");
  gStyle->SetOptFit(111);
  gStyle->SetStatFontSize(0.05);

  cherenkov->Draw();

  cout << "RMS:  " << cherenkov->GetRMS()  << " mrad" << endl;
  cout << "Mean: " << cherenkov->GetMean() << " mrad" << endl;


  // Fit
  TF1 *fit;

  if( opt_angleAcceptance )
  {
    fit = new TF1( "fit", "gaus(0)+pol0(3)" );
    fit->SetParNames("A","#mu","#sigma","c");
    fit->SetParameters( 1000, cherenkov->GetMean(), 10, 1 );
  }
  else
  {
    fit = new TF1( "fit", "gaus(0)" );
    fit->SetParNames("A","#mu","#sigma");
    fit->SetParameters( 1000, cherenkov->GetMean(), 10 );
  }

  fit->SetParLimits(2, 0, xmax_c - xmin_c ); // to avoid getting negative sigma
  fit->SetLineColor(2);

  cherenkov->Fit( "fit", "B", "", xmin_c, xmax_c );


  cout << (fit->GetParameter(2)) << " +- " << (fit->GetParError(2)) << " mrad" << endl;


  canvas->Update();
  Double_t ymax = canvas->GetFrame()->GetY2();

  TF1 *degreefunc = new TF1( "degreefunc", "x", xmin_c *mrad2deg, xmax_c *mrad2deg);
  TGaxis *degreeAxis = new TGaxis( xmin_c, ymax, xmax_c, ymax, "degreefunc", 510, "-");
  degreeAxis->SetLabelFont(132);
  degreeAxis->SetLabelSize(0.05);
  degreeAxis->SetTitle("#Theta_{c} [#circ]");
  degreeAxis->SetTitleFont(132);
  degreeAxis->SetTitleSize(0.06);
  degreeAxis->CenterTitle(true);
  degreeAxis->Draw();
}