//
// 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
//}

//===============================
//
//  ... TRY to do convolution : use numeric Conv
//
//=================================

#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 "RooGenericPdf.h"

#include "RooTFnBinding.h"
//#include "RooCFunction1Binding.h"
#include "RooFFTConvPdf.h"

#include "TTree.h"
#include "TH1D.h"




  using namespace std;

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


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

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


  extern void readN();               // read input file

  extern void histos();


  extern Double_t myFunc(Double_t ee);

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

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


  TF1 *bw;
  TF1 *voigt;

  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]);
  };

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


//  // TGraph to TF1
//  TGraph *gg;
//
//  double myfunc(double xx)
//  {
//    return gg->Eval(xx);
//  };
//
//--------------


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


  double myf1(double xx)
  {                                 // check this variant
    double a = 1.0;
    return a;
  };



  double myf2(double xx)
  {                                 // check this variant
    double a = 1.0;
    return a*xx;
  };


  double myf3(double xx)
  {                                 // check this variant
  // Lorentzian Peak function

    double p1= 0.5;
    //double p2= 0.1;
    double p2= 0.05;
    double p3= 4.15;
    
    //    return   (0.5*p1*p2/TMath::Pi()) / TMath::Max(1.e-10,
    //      (xx-p3)*(xx-p3)+ .25*p2*p2);
    
    return   (0.5*p1*p2/TMath::Pi()) / ((xx-p3)*(xx-p3) + .25*p2*p2);
  };




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

// 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("test.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;            // resonance
   xp[0]= 0;               // parameter

   //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;

   //   xx= 4.15;
   cout << "    .......myfunc(4.15) =  " << myfunc(xa,xp) << endl;


   cout << " " << endl;
 

     mygausn = new TF1("mygausn", "gausn", 1, 10);
     mygausn->SetParameters(1, 4.15, 0.10);         // initial values

   xx= 4.15;
   //   cout << " ...initial mygausn(4.15) " << mygausn->Eval(4.15)  << endl;
   cout << " ...initial mygausn(4.15) " << mygausn->Eval(xx)  << endl;
   cout << " " << endl;


   auto c1 = new TCanvas("c1", "mygausn", 700,500);
   mygausn->Draw();

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

     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);
   //   TF1 *func = new TF1("func","myfunc(x)",emin,emax,0);
   //   ... this variant looks for LIBRARY ! whet it is not present ..

   //   TF1 *func = new TF1("func","myFunc(x)",emin,emax,0);
   TF1 *func = new TF1("func",myfunc,emin,emax,0);
   // general definition of any function :
   // (*x, *par) ! -> so we need follow it !
   //

   //TF1 *func = new TF1("func",myfunc,emin,emax,0);
   // this variant is my definition : no LIBRARY

   cout << " " << endl;
   //   cout << " ...func->Eval(4.15) " <<  func->Eval(4.15) << endl;
   cout << " ...func->Eval(4.15) = " <<  func->Eval(xx) << endl;
   cout << " " << endl;


   auto c2 = new TCanvas("c2", "func", 700,500);
   func->Draw();


   //---------------------
   //=====================
   // 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", emin, emax);
   t.Print();

   cout << " " << endl;

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

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

   //  RooAbsReal* rfunc = bindFunction(myFunc,t) ;
   // ... it does not work... 

  cout << " " << endl;
  cout << " rfunc->Print() "  << endl;
  rfunc->Print() ;
  cout << " " << endl;


  t.setVal(1.4);
  cout << " rfunc.getVal(), t= 1.4 " <<  rfunc->getVal()  << endl;
  t.setVal(1);
  cout << " rfunc.getVal(), t= 1   " <<  rfunc->getVal()  << endl;
  t.setVal(4.15);
  cout << " rfunc.getVal(), t= 4.15 " <<  rfunc->getVal()  << endl;
  t.setVal(4.5);
  cout << " rfunc.getVal(), t= 4.5  " <<  rfunc->getVal()  << endl;
  t.setVal(4.6);
  cout << " rfunc.getVal(), t= 4.6  " <<  rfunc->getVal()  << endl;
  t.setVal(4.7);
  cout << " rfunc.getVal(), t= 4.7  " <<  rfunc->getVal()  << endl;
  t.setVal(10);
  cout << " rfunc.getVal(), t= 10 " <<  rfunc->getVal()  << endl;
  cout << " " << endl;



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

  // Define physics my rfunc p.d.f 
  //RooGenericPdf mypdf("mypdf",rfunc,t);
  //RooGenericPdf mypdf("mypdf","rfunc",t);
  //RooGenericPdf mypdf("mypdf","mypdf",rfunc,t);
  //RooGenericPdf mypdf("mypdf","mypdf","rfunc",t);
  //RooGenericPdf mypdf("mypdf","mypdf",rfunc,RooArgSet(t));
  //RooGenericPdf mypdf("mypdf","mypdf",rfunc,RooArgList(t));
  //  ---> error
  //RooGenericPdf mypdf("mypdf","mypdf",rfunc,RooArgSet(t));
  //RooGenericPdf mypdf("mypdf","mypdf",rfunc,RooArgList(t));
  //RooGenericPdf mypdf("mypdf",rfunc,RooArgList(t));
  //RooGenericPdf mypdf("mypdf","rfunc",RooArgList(t));
  //  ---> error
  //RooGenericPdf mypdf("mypdf",func,t);
  //  ---> error
  //  RooGenericPdf mypdf("mypdf",myfunc,t);

  //  RooGenericPdf mypdf("mypdf","myFunc(t)",t);
  //RooGenericPdf mypdf("mypdf","myfunc(t)",t);
  //RooGenericPdf mypdf("mypdf","func(t)",t);
  //RooGenericPdf mypdf("mypdf","rfunc(t)",t);
  // ...ura !...



  //  RooGenericPdf mypdf("mypdf","myf1(t)",t);

  //  RooGenericPdf mypdf("mypdf","myf2(t)",t);

  RooGenericPdf mypdf("mypdf","myf3(t)",t);


  cout << " ....myf3{4.15) = " <<  myf3(4.15)  << endl;
  cout << " ....myf3{4.3)  = " <<  myf3(4.3)  << endl;
  cout << " ....myf3{4.4)  = " <<  myf3(4.4)  << endl;
  cout << " ....myf3{4.5)  = " <<  myf3(4.5)  << endl;
  cout << " ....myf3{4.6)  = " <<  myf3(4.6)  << endl;
  cout << " ....myf3{4.7)  = " <<  myf3(4.7)  << endl;
  cout << " " << endl;


  t.setVal(4.15);
  cout << " ....mypdf.getVal()( t= 4.15) = " <<  mypdf.getVal()  << endl;
  t.setVal(4.3);
  cout << " ....mypdf.getVal()( t= 4.3)  = " <<  mypdf.getVal()  << endl;
  t.setVal(4.4);
  cout << " ....mypdf.getVal()( t= 4.4)  = " <<  mypdf.getVal()  << endl;
  t.setVal(4.5);
  cout << " ....mypdf.getVal()( t= 4.5) = " <<  mypdf.getVal()  << endl;
  t.setVal(4.6);
  cout << " ....mypdf.getVal()( t= 4.6) = " <<  mypdf.getVal()  << endl;
  t.setVal(4.7);
  cout << " ....mypdf.getVal()( t= 4.7) = " <<  mypdf.getVal()  << endl;
  cout << " " << endl;
 

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

   // Construct gauss(t,mg,sg)
   RooRealVar mg("mg", "mg", 0);
   //  RooRealVar mg("mg", "mg", 4.15);
   //   RooRealVar sg("sg", "sg", 0.01, 0.0001, 1);
   //   RooRealVar sg("sg", "sg", 0.1, 0.0001, 1);
   RooRealVar sg("sg", "sg", 0.1, 0.0001, 10.0);
   RooGaussian gauss("gauss", "gauss", t, mg, sg);
   
  t.setVal(4.15);
  cout << " ....gauss.getVal()( t= 4.15) = " <<  gauss.getVal()  << endl;
  t.setVal(4.5);
  cout << " ....gauss.getVal()( t= 4.5) = " <<  gauss.getVal()  << endl;
  t.setVal(4.6);
  cout << " ....gauss.getVal()( t= 4.6) = " <<  gauss.getVal()  << endl;
  t.setVal(4.7);
  cout << " ....gauss.getVal()( t= 4.7) = " <<  gauss.getVal()  << endl;
  cout << " " << endl;
 

   auto c3 = new TCanvas("c3", "gauss for conv.", 700,500);
   gauss->Draw();

   
  //--------------
  // now convolution :
  //

  // Define sampling frequency ==> 100 !  -> break at X= 4.6 !
  // t.setBins(10000) ;                 //  !  the same !
                                      // number of point is VERY large,
                                      // BUT break is here !
                                    //  !  the same !
   t.setBins(500) ;       
  //t.setBins(10000,"cache") ; 


  // Construct convolution
  //RooFFTConvPdf fxg("fxg","myfunc (X) mygausn",t,func,mygausn) ;
  //RooFFTConvPdf fxg("fxg","rfunc (X) gaus",t,rfunc,gauss) ;

  // RooFFTConvPdf fxg("fxg","mypdf (X) gaus",t,mypdf,gauss) ;
  //------------
  //  numeric...
  //
  //RooNumConvPdf fxg("fxg","fxg",t,myfunc,gauss);
  // --> error
  //RooNumConvPdf fxg("fxg","fxg",t,"myfunc",gauss);
  // --> error
  //RooNumConvPdf fxg("fxg","fxg",t,rfunc,gauss);
  // --> error
  //RooNumConvPdf fxg("fxg","fxg",t,"rfunc",gauss);
  // --> error

  RooNumConvPdf fxg("fxg","fxg",t,mypdf,gauss);
    // mypdf error...

  //  RooNumConvPdf fxg("fxg","fxg",t,rfunc,gauss);



 // Set the buffer fraction to zero to obtain a true cyclical convolution
  //  Mpsi.setBufferFraction(0) ;

  cout << " .... ConvPDF is ready ...."   << endl;
  cout << " " << endl;

  //  return;


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

  // check convolution : create hist for it !
  //  RooDataSet* data = fxg.generate(t,10000) ;
  RooDataSet* data = fxg.generate(t,50000) ;
  // RooDataSet* data = fxg.generate(t,1000) ;


  cout << " .... generate data for Conv ...."   << endl;
  cout << " " << endl;



  //  // Fit gxlx to data
  //  fxg.fitTo(*data) ;

  // Plot data, fitted p.d.f
  //  RooPlot* frame = t.frame(Title("myfunc (x) mygaussn convolution")) ;
  RooPlot* frame = t.frame(Title("data: myf3->mpdf (x) gauss conv.")) ;
  data->plotOn(frame) ;
  //fxg.plotOn(frame) ;

  cout << " .... plot data for Conv ...."   << endl;
  cout << " " << endl;

  //-------

  //RooPlot* frame1 = t.frame() ;
   //RooPlot* frame = t.frame(Title("TF1 bound as RooFit function"),2,9) ;
  //   RooPlot* frame1 = t.frame(Title("rfunc function")) ;
   RooPlot* frame1 = t.frame(Title("rfunc: func-> myfunc : graf == gg")) ;
   rfunc->plotOn(frame1) ;

   //----
   
   //RooPlot* frame2 = t.frame() ;
   //   RooPlot* frame2 = t.frame(Title(" mypdf function")) ;
   RooPlot* frame2 = t.frame(Title(" mypdf: myf3 and conv")) ;
   //mypdf->plotOn(frame2) ;
   mypdf.plotOn(frame2) ;
   //   fxg.plotOn(frame2) ;
   fxg.plotOn(frame2, LineColor(kRed), LineStyle(kDashed)) ;

   
   //---------------
   // ...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) ;   
  frame->Draw() ;
 

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

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


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





   //=========== 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;

      cout << xx << "  " << vv << endl;
    }

    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();



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


//double
// myFunc(double ee)
Double_t
 myFunc(Double_t ee)
{ //---------------
  // get  Sig for given E = ee
  // no parameters here
  //

  int nn;
  Double_t xx, xx1, xx2, f1, f2, ff, aa;

  //-------------------
  
  //nn= nW74;             // n Ei points
  //
  //xx1 = W74xx[0];
  //f1=   W74vv[0];
  //xx2 = W74xx[nn-1];
  //f2=   W74vv[nn-1];

  if(ifl == 0) nn= nW74;
  if(ifl == 1) nn= n10;
  if(ifl == 2) nn= n60;

  xx1= hh[ifl]->GetBinLowEdge(1);
  f1 = hh[ifl]->GetBinContent(1);
  xx2= hh[ifl]->GetBinLowEdge(nn-1);
  f2 = hh[ifl]->GetBinContent(nn-1);

  //cout << " -myFunc- , ee = " << ee << " nn= " << nn << endl; 
  //cout << "    x1= " << xx1 << " x2= " << xx2 << endl;

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

  //if( (ee < xx1) || (ee > xx2) ) return 0;
  if( (ee <= xx1) || (ee >= xx2) ) return 0;
  //                             ---------
  //----------------------------------------
   
  //for (int i=0; i < nn-1; i++) {    // find E interval
  for (int i=0; i < nn-2; i++) {    // find E interval
    //xx1= W74xx[i];
    //f1=  W74vv[i];
    //xx2= W74xx[i+1];
    //f2=  W74vv[i+1];

    xx1= hh[ifl]->GetBinLowEdge(i+1);
    f1 = hh[ifl]->GetBinContent(i+1);
    
    xx2= hh[ifl]->GetBinLowEdge(i+2);
    f2 = hh[ifl]->GetBinContent(i+2);
     
    
    if( (ee >= xx1) && (ee <= xx2) )  break;
    //                              ==========
  }

  aa = (ee - xx1) / (xx2 - xx1);
  ff= f1 + aa * (f2 - f1);              // linear interpolation

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

  return ff;

} // myFunc







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


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()

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