//
// 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 "TH2.h"
#include "TF1.h"
#include "TROOT.h"
#include "TStyle.h"
#include "TMath.h"
#include "TCanvas.h"

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


  using namespace std;

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


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

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


  extern void readN();               // read input file

  extern Double_t myFunc(Double_t ee);   // my function for N-res


  extern void histos();

  extern void getConv(Int_t ih,  Double_t temp);

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

  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

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


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

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


   //histos();                             // define histos

   

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


   histos();                             // define 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;

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

   cout << " " << endl;
   cout << "    .......myFunc(4.15) =  " << myFunc(4.15) << endl;
   cout << " " << endl;
 

   if(ifl == 0) {  // full region
      mygausn = new TF1("mygausn", "gausn", 1, 1e+09);
      //voigt = new TF1("voigt", "CONV(bw, mygausn)", 1, 1e+09);
      func = new TF1("func",myFunc,1,1e+09);

      //voigt = new TF1("voigt", "CONV(myFunc, mygausn)", 1, 1e+09);
      voigt = new TF1("voigt", "CONV(func, mygausn)", 1, 1e+09);
   }

   if(ifl == 1) {  // E = 1 - 10 eV
      mygausn = new TF1("mygausn", "gausn", 1, 10);
      //voigt = new TF1("voigt", "CONV(bw, mygausn)", 1, 10);
      func = new TF1("func",myFunc,1,10);

      //voigt = new TF1("voigt", "CONV(myFunc, mygausn)", 1, 10);
      voigt = new TF1("voigt", "CONV(func, mygausn)", 1, 10);
   }

   if(ifl == 0) {  // E = 1 - 60 eV
      mygausn = new TF1("mygausn", "gausn", 1, 60);
      //voigt = new TF1("voigt", "CONV(bw, mygausn)", 1, 60);
      func = new TF1("func",myFunc,1,60);

      //voigt = new TF1("voigt", "CONV(myFunc, mygausn)", 1, 60);
      voigt = new TF1("voigt", "CONV(func, mygausn)", 1, 60);
   }


   ih = ifl;                   // ifl = 1, E= 1- 10
  //ih = 12;                            // ifl = 2, E= 1-60 eV
  //ih = 10;                              // ifl = 0, E= all for W182

  if(ifl != 1)   
       getConv(ih, temp);             //  do convolution and fill histos


   //==========================
   // fill histos ...
   //

   if(ifl == 1 ) { 

     ih = 101;              // for ifl  =  1
     temp = 300.0;
     getConv(ih, temp);             //  do convolution and fill histos

     ih = 102;              
     temp = 500.0;
     getConv(ih, temp);             //  do convolution and fill histos

     ih = 103;            
     temp = 1000.0;
     getConv(ih, temp);             //  do convolution and fill histos

     //----------------
     //      ordinary scale
     // 1) create initial hist = 5
     //

 
     for (int i=0; i < nbin-1; i++) {
       xx= hh[5]->GetBinLowEdge(i+1);
       vv= myFunc(xx);
       hh[5]-> SetBinContent(i+1,vv);
     }

     // 2) call convolution
     //
   
     ih = 104;              // for ifl  =  1
     temp = 300.0;
     getConv(ih, temp);             //  do convolution and fill histos

     ih = 105;              
     temp = 500.0;
     getConv(ih, temp);             //  do convolution and fill histos

     ih = 106;            
     temp = 1000.0;
     getConv(ih, temp);             //  do convolution and fill histos


   } // ifl == 1



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

   //hh[1]->SetLineColor(kRed);   // kYellow, kBlue, kGreen, kMagenta, kBlack
   //hh[11]->SetLineColor(kBlue);   

   //hh[0]->SetLineColor(kRed);
   //hh[1]->SetLineColor(kRed);
   //hh[2]->SetLineColor(kRed);
   //
   //hh[10]->SetLineColor(kBlue);
   //hh[11]->SetLineColor(kBlue);
   //hh[12]->SetLineColor(kBlue);

   
   //----------------
   // now draw
   //

   //new TCanvas("c1","c1",800,1000);
   //TCanvas *c1 = new TCanvas("c1");
   //-> TCanvas *c1 = new TCanvas("c1","Page 1",800,600);
   //TCanvas *c1 = new TCanvas("c1","1-D ",600,800);
   //c1->Divide(2,1);
   //c1->SetLogx();
   //c1->SetLogy();

   //TCanvas *c1 = new TCanvas("c1","Function draw",10,10,600,450);
   //c1->Show();

   //gPad->SetGridx(); 
   //gPad->SetGridy(); 
   //gPad->SetLogx();

      // Global style settings
   //gStyle->SetOptStat(1111);
   //gStyle->SetOptFit(111);
   //gStyle->SetLabelSize(0.03,"x");
   //gStyle->SetLabelSize(0.03,"y");
 
   //hh->GetXaxis()->SetRange(1,10);
   //----------------


   //hh[1]->Draw();
   //hh[11]->Draw("same");   // hsame, lsame,

   //-> hh[ifl]->Draw();
   //-> hh[ifl+10]->Draw("same");   // hsame, lsame,

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

   //bw->Draw();                                           // red
   //mygausn->SetLineColor(1);      mygausn->Draw("same"); // black
   //voigt->SetLineColor(4);        voigt->Draw("same");    // blue

   //=========== 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 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
   getConv(Int_t ih,  Double_t temp) {
   //--------------------------------
   //  for given ih and T fill histos with conv. result
   //

  double xa[1];
  double xx, vv, vres, vt, conv, sigD;

  int nbin;
  //--------------
   //----------
   // axis->GetNbins();                  h->GetNbinsX();
   // axis->GetBinLowEdge(bin_number);   h->GetBinLowEdge(1)
   // axis->GetBinCenter(bin_number);    h->GetBinCenter(i)
   // axis->GetBinUpEdge(bin_number);    h->GetBinWidth(i)
   //  h->GetYaxis()->GetBinCenter(bin);   h->GetBinContent(i)
   //  h->GetYaxis()->GetBinCenter(biny);
   // bin_number = axis->FindBin(x_value);
   // axis->SetRange(firstbin,lastbin);
   //
   // h->SetBinContent(ibin,0);
   // b->SetBinContent(i+1, h->GetBinContent(i+1));
  //-------------------

   // 
   // sigma(Doppler) = sqrt(4* E* K_B*  T / A )
   //
   // use T(Debaya) = 3/8 *TET * coth[3/8 * TET/ T)
   // then
   //      sigma(Doppler) = 0.0186 * sqrt( E* Teff / (A+1) ),  eV
   //   TET = 315 for W
   //            4* 8.61733e-05 * 3./8. *315. /182 = 2.24e-04
   //
   // then  sigma(Doppler) = 0.0150 * sqrt( E / tanh(118.125 / T) ),  eV
   //       =====            =====================================    
   //
   // use tanh()   or  TMath::TanH() ...


  //--------------------
  // sigma Gauus with T :

  xa[0]= 4.15;

  sigD = 0.0150 * sqrt( xa[0] / tanh(118.125 / temp) );
  
  //-------------------------------------


   cout << " " << endl;
   cout << " ===== getConv:   hist=  " << ih << " T = " << temp << endl;
   if( ih > 103) cout << " =====   Ordinary scale "  << endl;
   cout << " ===========================================" << endl;
   cout << " " << endl;
   //--------------------------

   
   cout << " " << endl;
   cout << " myFunc(4.15)=  " << myFunc(4.15) << endl;
   cout << "   func(4.15)=  " << func->Eval(4.15) << endl;
   cout << " " << endl;

   mygausn->SetParameters(1, xa[0], sigD);
   //cout << " myGaus(4.15,4.15,temp)=  " << myGaus(xa,4.15,temp) << endl;
   cout << " myGaus(4.15,4.15,temp)=  " << mygausn->Eval(4.15) << endl;
   cout << " " << endl;

   cout << " voigt(4.15)=  " << voigt->Eval(4.15) << endl;

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


   
   //nbin = hh[1]->GetXaxis()->GetNbins(); 
   //nbin = hh[1]->GetNbinsX();   
   nbin = hh[ifl]->GetNbinsX();   

   cout << " ..nbin= " << nbin << endl;

   
   //=============================
   // fill also histos for GAUSS
   //
   if((ifl == 1) && (temp < 500.0) && (ih < 104)) {
                                   // when we have 1 resonance only

     
     cout << " " << endl;
     cout << "  ... for E= 4.15 eV : " << endl;
     cout << " " << endl;

     vt= 300.0;
     sigD = 0.0150 * sqrt( xa[0] / tanh(118.125 / vt) );
     mygausn->SetParameters(1, xa[0], sigD);
     cout << " -- T= 300, sigD = " << sigD << endl;
     //cout << "            myGauss(4.15) = " << myGaus(xa,4.15,vt)  << endl;
     cout << "            myGauss(4.15) = " <<  mygausn->Eval(4.15)  << endl;

     vt= 500.0;
     sigD = 0.0150 * sqrt( xa[0] / tanh(118.125 / vt) );
     mygausn->SetParameters(1, xa[0], sigD);
     cout << " -- T= 500, sigD = " << sigD << endl;
     //cout << "            myGauss(4.15) = " << myGaus(xa,4.15,vt)  << endl;
     cout << "            myGauss(4.15) = " << mygausn->Eval(4.15)  << endl;

     vt= 1000.0;
     sigD = 0.0150 * sqrt( xa[0] / tanh(118.125 / vt) );
     mygausn->SetParameters(1, xa[0], sigD);
     cout << " -- T= 1000, sigD = " << sigD << endl;
     //cout << "            myGauss(4.15) = " << myGaus(xa,4.15,vt)  << endl;
     cout << "            myGauss(4.15) = " << mygausn->Eval(4.15)  << endl;
     cout << " " << endl;
     
     for (int i=0; i < nbin-1; i++) {
       xx= hh[ifl]->GetBinLowEdge(i+1);
       xa[0]= xx;
       vres= 4.15;
       //-----
       
       vt=   300.0;
       sigD = 0.0150 * sqrt( xa[0] / tanh(118.125 / vt) );
       //vv= myGaus(xa,vres,vt);
       mygausn->SetParameters(1, vres, sigD);
       vv= mygausn->Eval(xx);
       hh[21]-> SetBinContent(i+1,vv);
       
       vt=   500.0;
       sigD = 0.0150 * sqrt( xa[0] / tanh(118.125 / vt) );
       //vv= myGaus(xa,vres,vt);
       mygausn->SetParameters(1, vres, sigD);
       vv= mygausn->Eval(xx);
       hh[22]-> SetBinContent(i+1,vv);
       
       vt=   1000.0;
       sigD = 0.0150 * sqrt( xa[0] / tanh(118.125 / vt) );
       //vv= myGaus(xa,vres,vt);
       mygausn->SetParameters(1, vres, sigD);
       vv= mygausn->Eval(xx);
       hh[23]-> SetBinContent(i+1,vv);
     }

     //-----------------
     // and ordinary scale


     cout << " ....... ordinary scale ..." << endl;
     
     for (int i=0; i < nbin-1; i++) {
       xx= hh[24]->GetBinLowEdge(i+1);
       xa[0]= xx;
       vres= 4.15;
       //----
       
       vt=   300.0;
       sigD = 0.0150 * sqrt( xa[0] / tanh(118.125 / vt) );
       //vv= myGaus(xa,vres,vt);
       mygausn->SetParameters(1, vres, sigD);
       vv= mygausn->Eval(xx);
       hh[24]-> SetBinContent(i+1,vv);
       
       vt=   500.0;
       sigD = 0.0150 * sqrt( xa[0] / tanh(118.125 / vt) );
       //vv= myGaus(xa,vres,vt);
       mygausn->SetParameters(1, vres, sigD);
       vv= mygausn->Eval(xx);
       hh[25]-> SetBinContent(i+1,vv);
       
       vt=   1000.0;
       sigD = 0.0150 * sqrt( xa[0] / tanh(118.125 / vt) );
       //vv= myGaus(xa,vres,vt);
       mygausn->SetParameters(1, vres, sigD);
       vv= mygausn->Eval(xx);
       hh[26]-> SetBinContent(i+1,vv);
     }
     
   }                         //  ifl = 1


   //
   //=============================
   // and main procedure
   //

   cout << " ....... main procedure ..." << endl;


   if( (ih >= 104 ) && (ih <= 106) ) { // ordinary scale

     nbin = hh[5]->GetNbinsX();
     for (int i=0; i < nbin-1; i++) {
       xx= hh[5]->GetBinLowEdge(i+1);
       vv= hh[5]->GetBinContent(i+1);

       //conv = myConvFun(xx,temp);
       xa[0]= xx;
       vt=   300.0;
       sigD = 0.0150 * sqrt( xa[0] / tanh(118.125 / vt) );
       mygausn->SetParameters(1, xx, sigD);
       //conv = myConvFun(xa,temp);
       conv = voigt->Eval(xx);

       //hh[11]-> SetBinContent(i+1,conv);
       //if(temp < 500.0) hh[ifl+10]-> SetBinContent(i+1,conv);
       hh[ih]-> SetBinContent(i+1,conv);
     }

   } else {            // input scale

     for (int i=0; i < nbin-1; i++) {
       //xx= hh[1]->GetBinLowEdge(i+1);
       //vv= hh[1]->GetBinContent(i+1);
       xx= hh[ifl]->GetBinLowEdge(i+1);
       vv= hh[ifl]->GetBinContent(i+1);

       //conv = myConvFun(xx,temp);
       xa[0]= xx;
       vt=   300.0;
       sigD = 0.0150 * sqrt( xa[0] / tanh(118.125 / vt) );
       mygausn->SetParameters(1, xx, sigD);
       //conv = myConvFun(xa,temp);
       conv = voigt->Eval(xx);

       //hh[11]-> SetBinContent(i+1,conv);
       if(temp < 500.0) hh[ifl+10]-> SetBinContent(i+1,conv);
       hh[ih]-> SetBinContent(i+1,conv);

       if( i%500 == 0) cout << "... i= " << i << endl;
     }

   } // if...


   cout << " " << endl;
   cout << " " << endl;
   cout << "========= END conv " << endl;
   cout << " " << endl;

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

   return;

 } // getConv   




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


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

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