//
// it is from BOOK:  TFormula  or RooFit ...
//
// Listing 3 presents an example of the use of this operator to perform 
//        the convolution of a Breit-Wigner function and a Gaussian function. 
//        The convolution can be performed as an fast Fourier transform 
//       convolution (by default) or by applying numerical integration.
//
//Listing 3: Code necessary for building a TF1 convolution of 
//           the Breit-Winger and Gaussian functions.
//---------------------------------------
//
// it works at the cluster
//
//
//{
//  TF1 *bw = new TF1("bw", "breitwigner", -15, 15);
//  bw->SetParameters(1, 0, 1);
//
//  TF1 *mygausn = new TF1("mygausn", "gausn", -15, 15);
//  mygausn->SetParameters(1, 0, 1);
//
//  TF1 *voigt = new TF1("voigt", "CONV(bw, mygausn)", -15, 15);
//
//  bw->Draw();                          // red
//  mygausn->SetLineColor(1);      mygausn->Draw("same"); // black
//  voigt->SetLineColor(4);        voigt->Draw("same");    // blue
//}

//===============================


#include "TH1.h"
#include "TH1F.h"
#include "TH2.h"
#include "TF1.h"
#include <TH1F.h>
#include <TGraph.h>
#include "TROOT.h"
#include "TStyle.h"
#include "TMath.h"
#include "TCanvas.h"

#include <TObject.h>
#include <TRandom.h>
#include <math.h>
#include <TF1Convolution.h>

#include <fstream>
#include <iostream>
#include "TFile.h"
#include "TTree.h"

#include <vector>
#include "TVectorD.h"


using namespace RooFit;
 
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooDataHist.h"
#include "RooGaussian.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "TTree.h"
#include "TH1D.h"




  using namespace std;

  //----------------------


  ifstream in1, in2;        // C++ - example

  //--------------------------


  extern void readN();               // read input file

  extern void histos();
  //-----------------------------

  TF1 *mygausn;
  //TF1 *mygausn = new TF1("mygausn", "gausn", 1, 1e+09);


  TF1 *bw;
  TF1 *voigt;

  //TF1* func = new TF1("Multiple",myFunc,xlo,xhi,nMipsMax+2);
  TF1* func;
  

  //----------------------------

  double W74xx[500000];                // for input data
  double W74vv[500000];

  double w1_10x[500];       // values for E = 1 - 10 eV
  double w1_10v[500];
  
  double w1_60x[1000];       // values for E = 1 - 60 eV
  double w1_60v[1000];

  int n10, n60;
  int    nW74;
 
  int ifl;                           // define input hist-array

  //-------------------------------


  TH1F  *hh[200];                        // --- general histos
  TH2F *hh2[200];                        // --- general histos

  //=========================
//////////////////////////////

  // TGraph to TF1
  TGraph *gg;

  double myfunc(double *xx, double *)
  //double myfunc(double *xx)
  {
   return gg->Eval(xx[0]);
  };



//===========================
////////////////////////////////////

// Main function

 void 
  Convolution()
 {//-----------------------------
  //  the start point of the program
  //

   int nbin, ih;

  //  float told, tnew;
  float rtime, ctime;
  
  //double sigD;
  double con1, con2, temp;
  double xlo,  xhi, xx, vv, conv;
  double xa[1], xp[1];

 //=====> initialize ROOT

 //   TROOT neut("Neutrons", "Read neutron data");
 //  gROOT->SetBatch();    // must be set before TApplication is created
 //  TApplication theApp("App", 0, 0);

   TFile* getN= 
          new TFile("exBwGausConv.root","RECREATE","read and convolutes");

   //=========================

            // ----  get to the file-scope with histograms !
            //-----------------------------------------

   getN->cd();   
   //getN.cd();   
 
              //--------------------------
              //  read input ascii data
              //  ... fill histos...
              //--------------------------


   readN();                   // read data here
                              // and fill histos ...




              //--------------------------
              //  now fill hist for my input hist
              //    convoluted with myGaus ...
              //--------------------------

   //                 check interval :   1 - 10 eV first:
   // T= 300 degrees K
   // N-res = 4.15 eV
   //

   temp= 300.0;
   //xx= 4.15;
   xa[0]= 4.15;
   xp[0]= 0;

   //ifl = 0;                   // use data total range, hist = 0
   ifl = 1;                   // use data for 1-10 eV, hist = 1
   //ifl = 2;                   // use data for 1-60 eV, hist = 2

   
   cout << " " << endl;
   cout << " flag =  " << ifl << endl;
   cout << " ---------------- " << endl;


   //------------------------------------

  


   histos();                             // define histos

   //---------------------

 
   cout << " " << endl;
   //cout << "    .......myFunc(4.15) =  " << myFunc(4.15) << endl;
   //cout << "    .......myfunc(4.15) =  " << myfunc(xa) << endl;
   cout << "    .......myfunc(4.15) =  " << myfunc(xa,xp) << endl;
   //cout << "    .......myFunc(4.15) =  " << myFunc(xa) << endl;
   cout << " " << endl;
 
     mygausn = new TF1("mygausn", "gausn", 1, 10);
     mygausn->SetParameters(1, 4.15, 0.10);         // initial values

     //Double_t emin = 1.0;
     //Double_t emax = 10.0;

     Double_t emin = 1.07151;
     Double_t emax = 9.81108;

     TF1 *func = new TF1("func","myfunc",emin,emax,0);

   cout << " " << endl;
   cout << " ...func->Eval(4.15) " <<  func->Eval(4.15) << endl;
   cout << " " << endl;
   cout << " ...initial mygausn(4.15) " << mygausn->Eval(4.15)  << endl;
   cout << " " << endl;


   //---------------------
   //=====================
   // S e t u p   c o m p o n e n t   p d f s
   // ---------------------------------------

   cout << " " << endl;
   cout << " Now RooFit "  << endl;
   cout << " =============" << endl;
   cout << " " << endl;
   cout << " " << endl;

   // Construct observable
   //RooRealVar t("t", "t", 1, 10);
   //RooRealVar t1("t1", "t1", 1, 60);

   RooRealVar t("t", "t", 1.07151, 9.81108);
   RooRealVar t1("t1", "t1", 1.07151, 59.7634);


  // Create an observable 
  //RooRealVar x3("x3","x3",0.01,20) ;

  // Create binding of TF1 object to above observable
  RooAbsReal* rfunc = bindFunction(func,t) ;

  // Print rfa1 definition
  cout << " " << endl;
  cout << " rfunc->Print() "  << endl;
  rfunc->Print() ;
  cout << " " << endl;

  //------------------------

   cout << " " << endl;
   cout << " Now PLOT "  << endl;
   cout << " =============" << endl;
   cout << " " << endl;

   //------- ROOT method
   //
   //  TCanvas* c1 = new TCanvas("c1, ROOT example","my c1",800,600) ;
   //
   //  gPad->SetLeftMargin(0.15) ; 
   //  hh[1]->GetYaxis()->SetTitleOffset(1.6) ; 
   //  hh[1]->Draw() ;                                // plot 1-10
   //  ...


   //------- RooFit method
   // 1) first plot it into frame...
   // 2) then canvas and plot...
   //

   //RooPlot* frame = t.frame(Title("TF1 bound as RooFit function")) ;
   //rfunc->plotOn(frame) ;


   // Create a binned dataset that imports contents of TH1 and 
   //          associates its contents to observable 'x'
   //   RooDataHist dh1("dh1","dh1",t,Import(*hh[1])) ;    //  1-10
   RooDataHist dh1("dh1","dh1",t,Import(*hh[1],kFALSE)) ;    //  1-10

   RooDataHist dh2("dh2","dh2",t1,Import(*hh[2])) ;    //  1-60

   // ...now frames: 1-3
   //
  // If histogram has custom error (i.e. its contents is does not originate 
  //  from a Poisson process
  // but e.g. is a sum of weighted events) you can data with symmetric
   //  'sum-of-weights' error instead
  // (same error bars as shown by ROOT)

   RooPlot* frame = t.frame(Title("TF1 bound as RooFit function")) ;
   rfunc->plotOn(frame) ;

   RooPlot* frame1 = t.frame(Title(" 1-10 range")) ;
   dh1.plotOn(frame1,DataError(RooAbsData::SumW2)) ; 
   dh1.plotOn(frame1) ;
   //dh1->plotOn(frame1) ;

   RooPlot* frame2 = t1.frame(Title(" 1-60 range")) ;
   dh2.plotOn(frame2,DataError(RooAbsData::SumW2)) ; 
   dh2.plotOn(frame2) ;
   //dh2->plotOn(frame2) ;

   // ...now canvas: 1-3
   //
   //TCanvas* c = new TCanvas("root-example","rf105_funcbinding",1200,400) ;
  TCanvas* c = new TCanvas("root-example","rf105_funcbinding",800,800) ;
  //c->Divide(3) ;
  //c->Divide(2,3) ;
  c->Divide(2,2) ;

 
  c->cd(1) ; 
  gPad->SetLeftMargin(0.15) ; 
  frame->GetYaxis()->SetTitleOffset(1.6) ;           // plot rfunc !
  frame->Draw() ;
 

  c->cd(2) ; 
  gPad->SetLeftMargin(0.15) ; 
  frame1->GetYaxis()->SetTitleOffset(1.6) ;            // plot 1-10
  frame1->Draw() ;

  c->cd(3) ; 
  gPad->SetLeftMargin(0.15) ; 
  frame2->GetYaxis()->SetTitleOffset(1.6) ;            // plot 1-60
  frame2->Draw() ;


  c->cd(4) ; 
  gPad->SetLeftMargin(0.15) ; 
  hh[1]->GetYaxis()->SetTitleOffset(1.6) ;            // and naitive hist
  hh[1]->Draw() ;





   //  // Make plot frame in observable, plot TF1 binding function
   //  TCanvas* c1 = new TCanvas("c1","my c1",800,600) ;

  //  RooPlot* frame = t.frame(Title("TF1 bound as RooFit function")) ;
  // rfunc->plotOn(frame) ;
  //
  //  RooPlot* frame1 = t.frame(Title(" for my hist")) ;
  //  hh[1]->plotOn(frame1);
  //
  //  RooPlot* frame2 = t.frame(Title(" for my hist")) ;
  //  hh[2]->plotOn(frame2);

  //----------------

  //  TCanvas* c = new TCanvas("root-example","rf105_funcbinding",1200,400) ;
  //  c->Divide(3) ;
  //c->Divide(2,3) ;
 
  //  c->cd(1) ; 
 // gPad->SetLeftMargin(0.15) ; 
  //  frame->GetYaxis()->SetTitleOffset(1.6) ;           // plot rfunc !
  //  frame->Draw() ;
  //hh[0]->Draw();  

  //  c->cd(2) ; 
  //  gPad->SetLeftMargin(0.15) ; 
  //  frame1->GetYaxis()->SetTitleOffset(1.6) ; 
  //  frame1->Draw() ;
  //hh[1]->Draw();                                   // plot 1-10

  //  c->cd(3) ; 
  //  gPad->SetLeftMargin(0.15) ; 
  //  frame2->GetYaxis()->SetTitleOffset(1.6) ; 
  //  frame2->Draw() ;
  //hh[2]->Draw();                                   // plot 1-60

















   //=========== close calculations and save results...

  //   getN->Write();
  //   getN->Close();            // do not close it if You want look at.. them

   //getN.Write();
   //getN.Close();            // do not close it if You want look at.. them
   
   //=============

   return;

 }; // Convolution()



//========================================
//////////////////////////////////////////

void 
readN()
{ //-----------------------------
  //-----------------------------
  //       Input description :   format
  //
  //          E(eV) - V (b) :  (crossection)
  //-------------------------------------

   //--------------------
   // read data file: ../data/W182-in-W-R-full.txt
   // ... examples at Resonance/G4-W/
   //----------------------------------

  int ntot;

  double xx, vv;

  TString aaa[20];

  //double w1_10x[10000];       // values for E = 1 - 10 eV
  //double w1_10v[10000];
  //
  //double w1_60x[10000];       // values for E = 1 - 60 eV
  //double w1_60v[10000];
  //
  //int n10, n60;

  //----------------------------


  //aaa[0]= "../data/W182-in-W-R-full.txt";      //  W182 
  aaa[0]= "data/W182-in-W-R-full.txt";      //  W182 

  in1.open(aaa[0], ios::in);       // open input file
  //-----------------------------

      if(!in1.is_open()) {
        cout << "  "  << endl;
        cout << " !!!! kak tak !!! input is  NOT open !!!  "  << endl;
        cout << "  "  << endl;

        return;
	//=========   
      } else {
        cout << "  "  << endl;
        cout << " ... input FILE OPEN : "  << endl;
      }

  cout << " - file : " << aaa[0] << endl;
  cout << "  "  << endl;

  //------------------


  ntot = 0;

  n10 = 0;
  n60 = 0;

  while (1) {                           // ... LOOP on input lines...
  //-----------------
    ntot += 1;

    //if( ntot < 2) {
    //   getline( in1 , line ) ;
    //   cout << " n= "<<  endl;
    //   cout << "    file: " << ntot << line << endl;
    //
    //   continue;                       // skip this line
    //   //--------
    // }
    //-----------------

    in1 >> xx >> vv ;      //  read line ...

    if( ntot < 5) {
      printf(" %14.8e   %14.8e  \n", xx,vv);
    }

    if (!in1.good()) break;                // end of LOOP
    //------------------------------

    if(vv > 0.0) {      
      //  W74[ntot-1][0]= xx;          
      //  W74[ntot-1][1]= vv;          
      W74xx[ntot-1]= xx; 
      W74vv[ntot-1]= vv; 
    }

    if( (xx >= 1.0) &&  (xx <= 10.0) ) {
      n10 += 1;
      w1_10x[n10-1] = xx;
      w1_10v[n10-1] = vv;
    }

    if( (xx >= 1.0) &&  (xx <= 60.0) ) {
      n60 += 1;
      w1_60x[n60-1] = xx;
      w1_60v[n60-1] = vv;
    }

    //-----------------

  } // while

  cout << "  "  << endl;
  cout << " N-total : " << ntot-1 << endl;
  cout << "  "  << endl;
  cout << " N-1-10 : " << n10-1 << endl;
  cout << "  "  << endl;
  cout << " N-1-60 : " << n60-1 << endl;
  cout << "  "  << endl;

  nW74= ntot - 1;

  //--------------------


  in1.close();

  //--------------


  //---------------------
  // fill graph
  //

  TVectorD ene(n10), sigma(n10);

  for (int i=0; i < n10; i++) {
    ene(i)   =  w1_10x[i];
    sigma(i) =  w1_10v[i];
  }

  gg = new TGraph(ene,sigma);

  // check
 
  cout << " " << endl;
  cout << " ... from readN ... " << endl;
  cout << " " << endl;
  cout << " - graph gg : gg(1.0)=  " << gg->Eval(1.0) << endl;
  cout << " - graph gg : gg(4.15)=  " << gg->Eval(4.15) << endl;
  cout << " - graph gg : gg(10.0)=  " << gg->Eval(10.0) << endl;
  cout << " " << endl;


  //--------------------
  // fill hist with these data
  //


  hh[0]= new TH1F("hh_00", "G4, 74W-182, full, n vs E, eV", nW74-1, W74xx);
  hh[1]= new TH1F("hh_01", "G4, 74W-182, 1-10, n vs E, eV", n10-1, w1_10x);
  hh[2]= new TH1F("hh_02", "G4, 74W-182, 1-60, n vs E, eV", n60-1, w1_60x);



  //for (int i=0; i < 10; i++) {
  //  cout << " x= " << w1_10x[i] << " V =  " << w1_10v[i] << endl;
  //}
  //cout << " " << endl;


  for (int i=0; i < nW74; i++) {
    hh[0]->SetBinContent(i+1, W74vv[i]);
  }

  for (int i=0; i < n10; i++) {
    hh[1]->SetBinContent(i+1, w1_10v[i]);
  }

  for (int i=0; i < n60; i++) {
    hh[2]->SetBinContent(i+1, w1_60v[i]);
  }


  //for (int i=0; i < 10; i++) {
  //  cout << " x= " << hh[1]->GetBinLowEdge(i+1) << " V =  " << 
  //                    hh[1]->GetBinContent(i+1) << endl;
  //}


  //---------------------
   cout << "     ----------------  " << endl;
   cout << "       " << endl;

  //         ============= end of reading ==============
  //         ============================================

  return;

} // end of readN();



//========================================
//////////////////////////////////////////


void 
 histos()
{ //---------------
  // create the tree of directories:
  //  .... sector -> SM ....
  //-------------------------------------

  //   char dirname[50];
  //   char  htitle[80];
  //


   printf( "  --- HISTOS -----   \n ");
   printf( "     \n ");


//---------------------------------------------------


   //hh[0]= new TH1F("hh_00", "G4, 74W-182, full, n vs E, eV", nW74-1, W74xx);
   //hh[1]= new TH1F("hh_01", "G4, 74W-182, 1-10, n vs E, eV", n10-1, w1_10x);
   //hh[2]= new TH1F("hh_02", "G4, 74W-182, 1-60, n vs E, eV", n60-1, w1_60x);



  hh[10]= new TH1F("hh_10", "G4, 74W-182, full, convol, T=300", 
                                                            nW74-1, W74xx);
  hh[11]= new TH1F("hh_11", "G4, 74W-182, 1-10, convol. T= 300", 
                                                            n10-1, w1_10x);
  hh[12]= new TH1F("hh_12", "G4, 74W-182, 1-60, convol. T= 300", 
                                                            n60-1, w1_60x);

  //---------------

  if( ifl == 1) {

    hh[5]= new TH1F("hh_05", " simple Res= 4.15, n vs E, eV", 200, 3.6,4.6);


    hh[21]= new TH1F("hh_21", " Gauss-4.15, T= 300, n vs E, eV", 
                                                        n10-1, w1_10x);
    hh[22]= new TH1F("hh_22", " Gauss-4.15, T= 500, n vs E, eV", 
                                                        n10-1, w1_10x);
    hh[23]= new TH1F("hh_23", " Gauss-4.15, T= 1000, n vs E, eV", 
                                                        n10-1, w1_10x);

     
    hh[24]= new TH1F("hh_24", " Gauss-simple-4.15, T= 300, n vs E, eV", 
                                                         200, 3.6, 4.6);
    hh[25]= new TH1F("hh_25", " Gauss-simple-4.15, T= 500, n vs E, eV",  
                                                         200, 3.6, 4.6);
    hh[26]= new TH1F("hh_26", " Gauss-simple-4.15, T= 1000, n vs E, eV", 
                                                         200, 3.6, 4.6);



    hh[101]= new TH1F("hh_101", " ifl= 1, T= 300, n vs E, eV", n10-1, w1_10x);
    hh[102]= new TH1F("hh_102", " ifl= 1, T= 500, n vs E, eV", n10-1, w1_10x);
    hh[103]= new TH1F("hh_103", " ifl= 1, T= 1000, n vs E, eV", n10-1, w1_10x);


    hh[104]= new TH1F("hh_104", " Scale-Ordinary-4.15, T= 300, n vs E, eV", 
                                                      200, 3.6,4.6);
    hh[105]= new TH1F("hh_105", " Scale-Ordinary-4.15, T= 500, n vs E, eV",  
                                                      200, 3.6,4.6);
    hh[106]= new TH1F("hh_106", " Scale-Ordinary-4.15, T= 1000, n vs E, eV", 
                                                      200, 3.6,4.6);
  }

   //-----------------------------------------
   //
   // check bins of histos
   //

   //   TAxis* xAxisE = hh[41]->GetXaxis(); 
   //
   //   cout << "   2 bins on E-hist: " << endl;      
   //   cout << "\tLower edge: " << xAxisE->GetBinLowEdge(1) << endl;
   //   cout << "          \tUpper edge: " << xAxisE->GetBinUpEdge(1) << endl;
   //   cout << "\tLower edge: " << xAxisE->GetBinLowEdge(2) << endl;
   //   cout << "          \tUpper edge: " << xAxisE->GetBinUpEdge(2) << endl;
   //   cout << " " << endl;
   //
   //   TAxis* xAxisT = hh[43]->GetXaxis();
   //
   //   cout << "  2 bins on T-hist: " << endl;      
   //   cout << "\tLower edge: " << xAxisT->GetBinLowEdge(1) << endl;
   //   cout << "          \tUpper edge: " << xAxisT->GetBinUpEdge(1) << endl;
   //   cout << "\tLower edge: " << xAxisT->GetBinLowEdge(2) << endl;
   //   cout << "          \tUpper edge: " << xAxisT->GetBinUpEdge(2) << endl;
   //   cout << " " << endl;

   //--------------------------------------------

  //-----------------------------------------

   //   hTop->cd();                              // the h-tree head

 
   //--------------------------------

 
   cout << "   ... histogrames are booked..." << endl;
   cout << " " << endl;


  return;

} // histos()

//===============================
//////////////////////////////////////////