mcr/calib/tofino/MakePeakTable.C

/*
======================================================================
 FillTofFromHdt(const TString &runtype,const Int_t &run_num_start,const Int_t &run_num_stop,const TString &root_filename,const Int_t &nevs, Int_t B)
   After reading from the specified range of hdt files.
   * fill the ADC histograms for TOFINO
   * find the pedestal, the Bhabha(=LowTh*^HighTh) peak, the Kaon(=HighTh) peak
   The histogram is saved in root_filename.
   
 FillTofFromHdt(const TString &runtype,const Int_t &run_num_start,const Int_t &run_num_stop,const TString &root_filename,const Int_t &nevs)
   If the magnetic field is not given,
   1 Tesla (FINU) or 0 Tesla (COSM) is given by default.

 MakePeakTable(const TString &runtype,const Int_t &run_num_start,const Int_t &run_num_stop, const TString &directory)
   "Peak Table" (which includes the information on pedetal, Bhabha and Kaon
   peak for each channel) is produced run-by-run.

 TofinoStability(const TString &runtype,const Int_t &run_num_start,const Int_t &run_num_stop,const TString &root_filename)
   By executing FillTofFromHdt for each run, the trends of pedestal,
   Bhabha and Kaon peak are saved in histograms.
   It is recommended to use StabilityFromPeakTable instead of TofinoStability,
   if you have Peak Tables already.

void TofinoStability(const TString &runtype,const Int_t &run_num_start,const Int_t &run_num_stop, const TString &calib_filename, const TString &root_filename)
   If you have calibration table, the trends of equalized Bhabha and Kaon peaks
   are also saved. The spectum of sqrt(ADCe*ADCp) (ADC is equalized) is also
   created (summing up all the events in the period), which is useful in
   determining the threshold for kaon selection.
======================================================================
*/


#include "TH2.h"
#include "TCanvas.h"
#include "TFile.h"
#include "TApplication.h"
#include "TTree.h"
#include "TF1.h"
#include "TNamed.h"
#include "TGraph.h"
#include "TProfile.h"

#include "TFndHdt.h"
//#include "FROOT.h"

using namespace FROOT;

Double_t langaufun(Double_t *x, Double_t *par) {

   //Fit parameters:
   //par[0]=Width (scale) parameter of Landau density
   //par[1]=Most Probable (MP, location) parameter of Landau density
   //par[2]=Total area (integral -inf to inf, normalization constant)
   //par[3]=Width (sigma) of convoluted Gaussian function
   //
   //In the Landau distribution (represented by the CERNLIB approximation),
   //the maximum is located at x=-0.22278298 with the location parameter=0.
   //This shift is corrected within this function, so that the actual
   //maximum is identical to the MP parameter.

      // Numeric constants
      Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
      Double_t mpshift  = -0.22278298;       // Landau maximum location

      // Control constants
      Double_t np = 100.0;      // number of convolution steps
      Double_t sc =   5.0;      // convolution extends to +-sc Gaussian sigmas

      // Variables
      Double_t xx;
      Double_t mpc;
      Double_t fland;
      Double_t sum = 0.0;
      Double_t xlow,xupp;
      Double_t step;
      Double_t i;


      // MP shift correction
      mpc = par[1] - mpshift * par[0];

      // Range of convolution integral
      xlow = x[0] - sc * par[3];
      xupp = x[0] + sc * par[3];

      step = (xupp-xlow) / np;

      // Convolution integral of Landau and Gaussian by sum
      for(i=1.0; i<=np/2; i++) {
         xx = xlow + (i-.5) * step;
         fland = TMath::Landau(xx,mpc,par[0]) / par[0];
         sum += fland * TMath::Gaus(x[0],xx,par[3]);

         xx = xupp - (i-.5) * step;
         fland = TMath::Landau(xx,mpc,par[0]) / par[0];
         sum += fland * TMath::Gaus(x[0],xx,par[3]);
      }

      return (par[2] * step * sum * invsq2pi / par[3]);
}



TF1 *langaufit(TH1F *his, Double_t *fitrange, Double_t *startvalues, Double_t *parlimitslo, Double_t *parlimitshi, Double_t *fitparams, Double_t *fiterrors, Double_t *ChiSqr, Int_t *NDF)
{
   // Once again, here are the Landau * Gaussian parameters:
   //   par[0]=Width (scale) parameter of Landau density
   //   par[1]=Most Probable (MP, location) parameter of Landau density
   //   par[2]=Total area (integral -inf to inf, normalization constant)
   //   par[3]=Width (sigma) of convoluted Gaussian function
   //
   // Variables for langaufit call:
   //   his             histogram to fit
   //   fitrange[2]     lo and hi boundaries of fit range
   //   startvalues[4]  reasonable start values for the fit
   //   parlimitslo[4]  lower parameter limits
   //   parlimitshi[4]  upper parameter limits
   //   fitparams[4]    returns the final fit parameters
   //   fiterrors[4]    returns the final fit errors
   //   ChiSqr          returns the chi square
   //   NDF             returns ndf

   Int_t i;
   Char_t FunName[100];

   sprintf(FunName,"Fitfcn_%s",his->GetName());

   TF1 *ffitold = (TF1*)gROOT->GetListOfFunctions()->FindObject(FunName);
   if (ffitold) delete ffitold;

   TF1 *ffit = new TF1(FunName,langaufun,fitrange[0],fitrange[1],4);
   ffit->SetParameters(startvalues);
   ffit->SetParNames("Width","MP","Area","GSigma");

   for (i=0; i<4; i++) {
      ffit->SetParLimits(i, parlimitslo[i], parlimitshi[i]);
   }

   his->Fit(FunName,"RB");   // fit within specified range, use ParLimits, do not plot

   ffit->GetParameters(fitparams);    // obtain fit parameters
   for (i=0; i<4; i++) {
      fiterrors[i] = ffit->GetParError(i);     // obtain fit parameter errors
   }
   ChiSqr[0] = ffit->GetChisquare();  // obtain chi^2
   NDF[0] = ffit->GetNDF();           // obtain ndf

   return (ffit);              // return fit function

}

TTree *fCalibTree;

TH1F *h_LowBhaMult;
TH1F *h_LowHypMult;
TH1F *h_HighBhaMult;
TH1F *h_HighHypMult;

TH1F *h_adce[K_N_TOFINO_SLABS];
TH1F *h_adce_ped[K_N_TOFINO_SLABS];
TH1F *h_adce_bha[K_N_TOFINO_SLABS];
TH1F *h_adce_kaon[K_N_TOFINO_SLABS];
TH1F *h_adce_bha_nosel[K_N_TOFINO_SLABS];
TH1F *h_adce_kaon_nosel[K_N_TOFINO_SLABS];
TH1F *h_adce_bhakaon[K_N_TOFINO_SLABS];
TH1F *h_adce_subtr[K_N_TOFINO_SLABS];
TH1F *h_adce_ped_subtr[K_N_TOFINO_SLABS];
TH1F *h_adce_bha_subtr[K_N_TOFINO_SLABS];
TH1F *h_adce_kaon_subtr[K_N_TOFINO_SLABS];
TH1F *h_adce_bha_subtr_nosel[K_N_TOFINO_SLABS];
TH1F *h_adce_kaon_subtr_nosel[K_N_TOFINO_SLABS];
TH1F *h_adce_bhakaon_subtr[K_N_TOFINO_SLABS];
TH1F *h_adce_equ[K_N_TOFINO_SLABS];
TH1F *h_adce_ped_equ[K_N_TOFINO_SLABS];
TH1F *h_adce_bha_equ[K_N_TOFINO_SLABS];
TH1F *h_adce_kaon_equ[K_N_TOFINO_SLABS];
TH1F *h_adce_bha_equ_nosel[K_N_TOFINO_SLABS];
TH1F *h_adce_kaon_equ_nosel[K_N_TOFINO_SLABS];
TH1F *h_adce_bhakaon_equ[K_N_TOFINO_SLABS];

TH1F *h_adcp[K_N_TOFINO_SLABS];
TH1F *h_adcp_ped[K_N_TOFINO_SLABS];
TH1F *h_adcp_bha[K_N_TOFINO_SLABS];
TH1F *h_adcp_kaon[K_N_TOFINO_SLABS];
TH1F *h_adcp_bha_nosel[K_N_TOFINO_SLABS];
TH1F *h_adcp_kaon_nosel[K_N_TOFINO_SLABS];
TH1F *h_adcp_bhakaon[K_N_TOFINO_SLABS];
TH1F *h_adcp_subtr[K_N_TOFINO_SLABS];
TH1F *h_adcp_ped_subtr[K_N_TOFINO_SLABS];
TH1F *h_adcp_bha_subtr[K_N_TOFINO_SLABS];
TH1F *h_adcp_kaon_subtr[K_N_TOFINO_SLABS];
TH1F *h_adcp_bha_subtr_nosel[K_N_TOFINO_SLABS];
TH1F *h_adcp_kaon_subtr_nosel[K_N_TOFINO_SLABS];
TH1F *h_adcp_bhakaon_subtr[K_N_TOFINO_SLABS];
TH1F *h_adcp_equ[K_N_TOFINO_SLABS];
TH1F *h_adcp_ped_equ[K_N_TOFINO_SLABS];
TH1F *h_adcp_bha_equ[K_N_TOFINO_SLABS];
TH1F *h_adcp_kaon_equ[K_N_TOFINO_SLABS];
TH1F *h_adcp_bha_equ_nosel[K_N_TOFINO_SLABS];
TH1F *h_adcp_kaon_equ_nosel[K_N_TOFINO_SLABS];
TH1F *h_adcp_bhakaon_equ[K_N_TOFINO_SLABS];

TH1F *h_lowtdce[K_N_TOFINO_SLABS];
TH1F *h_lowtdcp[K_N_TOFINO_SLABS];
TH1F *h_hightdce[K_N_TOFINO_SLABS];
TH1F *h_hightdcp[K_N_TOFINO_SLABS];

TH1F *h_adc_bha_sqr[K_N_TOFINO_SLABS];
TH1F *h_adc_kaon_sqr[K_N_TOFINO_SLABS];
TH1F *h_adc_bhakaon_sqr[K_N_TOFINO_SLABS];
TH1F *h_adc_bha_rawsqr[K_N_TOFINO_SLABS];
TH1F *h_adc_kaon_rawsqr[K_N_TOFINO_SLABS];
TH1F *h_adc_bhakaon_rawsqr[K_N_TOFINO_SLABS];

TH1F *h_adce_bha_equ_all[K_N_TOFINO_SLABS];
TH1F *h_adcp_bha_equ_all[K_N_TOFINO_SLABS];
TH1F *h_adc_bha_sqr_all[K_N_TOFINO_SLABS];
TH1F *h_adce_kaon_equ_all[K_N_TOFINO_SLABS];
TH1F *h_adcp_kaon_equ_all[K_N_TOFINO_SLABS];
TH1F *h_adc_kaon_sqr_all[K_N_TOFINO_SLABS];

TH1F *p_ped_e[12];
TH1F *p_ped_p[12];
TH1F *p_bha_e[12];
TH1F *p_bha_p[12];
TH1F *p_kaon_e[12];
TH1F *p_kaon_p[12];
TH1F *p_eq_e[12];
TH1F *p_eq_p[12];

TH1F *p_EQbha_e[12];
TH1F *p_EQbha_p[12];
//TH1F *p_EQbha_sqr[12];
TH1F *p_EQkaon_e[12];
TH1F *p_EQkaon_p[12];
//TH1F *p_EQkaon_sqr[12];

TH1F *p_EQ_e[12];
TH1F *p_EQ_p[12];
TH1F *p_OverKaonThreshold[12];
Int_t CurrentRunNumber;

Int_t start_unixtime, stop_unixtime;
Float_t ped_e[K_N_TOFINO_SLABS];
Float_t ped_p[K_N_TOFINO_SLABS];
Float_t bhapeak_e[K_N_TOFINO_SLABS];
Float_t bhapeak_p[K_N_TOFINO_SLABS];
Float_t kaonpeak_e[K_N_TOFINO_SLABS];
Float_t kaonpeak_p[K_N_TOFINO_SLABS];
Float_t eq_e[K_N_TOFINO_SLABS];
Float_t eq_p[K_N_TOFINO_SLABS];

Float_t ped_table_e[K_N_TOFINO_SLABS];
Float_t ped_table_p[K_N_TOFINO_SLABS];
Float_t eq_table_e[K_N_TOFINO_SLABS];
Float_t eq_table_p[K_N_TOFINO_SLABS];

Int_t KaonThreshold[K_N_TOFINO_SLABS];
Int_t NevOverThreshold[K_N_TOFINO_SLABS];

TFile *fHdtFile = 0;

TFndHdt *fCurHdtEv = 0;
TTree *fHdtTree = 0;

void DeclareHistos(){

  TString hname;
  TString htitle;
  
  hname.Form("LowBhaMult");
  htitle.Form("Multiplicity (Low threshold & Bhabha trigger)");
  h_LowBhaMult = new TH1F(hname,htitle,12,0,12);
  Printf("Histo name: \"%s\"",hname.Data());

  hname.Form("LowHypMult");
  htitle.Form("Multiplicity (Low threshold & Hyp trigger)");
  h_LowHypMult = new TH1F(hname,htitle,12,0,12);
  Printf("Histo name: \"%s\"",hname.Data());

  hname.Form("HighBhaMult");
  htitle.Form("Multiplicity (High threshold & Bhabha trigger)");
  h_HighBhaMult = new TH1F(hname,htitle,12,0,12);
  Printf("Histo name: \"%s\"",hname.Data());
  
  hname.Form("HighHypMult");
  htitle.Form("Multiplicity (High threshold Hyp trigger)");
  h_HighHypMult = new TH1F(hname,htitle,12,0,12);
  Printf("Histo name: \"%s\"",hname.Data());

  for(Int_t sl=0;sl<K_N_TOFINO_SLABS;sl++ ){

    hname.Form("TofinoAdcE_%d",sl+1);
    htitle.Form("Tofino slab: %d e",sl+1);
    h_adce[sl] = new TH1F(hname,htitle,1000,0,4000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcEPed_%d",sl+1);
    htitle.Form("Tofino Pedestal slab: %d e",sl+1);
    h_adce_ped[sl] = new TH1F(hname,htitle,1000,0,4000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcEBha_%d",sl+1);
    htitle.Form("Tofino Bhabha slab: %d e",sl+1);
    h_adce_bha[sl] = new TH1F(hname,htitle,1000,0,4000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcEKaon_%d",sl+1);
    htitle.Form("Tofino Kaon slab: %d e",sl+1);
    h_adce_kaon[sl] = new TH1F(hname,htitle,1000,0,4000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcEBhaNosel_%d",sl+1);
    htitle.Form("Tofino Bhabha slab: %d e",sl+1);
    h_adce_bha_nosel[sl] = new TH1F(hname,htitle,1000,0,4000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcEKaonNosel_%d",sl+1);
    htitle.Form("Tofino Kaon slab: %d e",sl+1);
    h_adce_kaon_nosel[sl] = new TH1F(hname,htitle,1000,0,4000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcEBhaKaon_%d",sl+1);
    htitle.Form("Tofino Bhabha+Kaon slab: %d e",sl+1);
    h_adce_bhakaon[sl] = new TH1F(hname,htitle,1000,0,4000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcESubtr_%d",sl+1);
    htitle.Form("Tofino subtracted slab: %d e",sl+1);
    h_adce_subtr[sl] = new TH1F(hname,htitle,1050,-200,4000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcEPedSubtr_%d",sl+1);
    htitle.Form("Tofino Pedestal subtracted slab: %d e",sl+1);
    h_adce_ped_subtr[sl] = new TH1F(hname,htitle,1050,-200,4000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcEBhaSubtr_%d",sl+1);
    htitle.Form("Tofino Bhabha subtracted slab: %d e",sl+1);
    h_adce_bha_subtr[sl] = new TH1F(hname,htitle,1050,-200,4000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcEKaonSubtr_%d",sl+1);
    htitle.Form("Tofino Kaon subtracted slab: %d e",sl+1);
    h_adce_kaon_subtr[sl] = new TH1F(hname,htitle,1050,-200,4000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcEBhaSubtrNosel_%d",sl+1);
    htitle.Form("Tofino Bhabha subtracted slab: %d e",sl+1);
    h_adce_bha_subtr_nosel[sl] = new TH1F(hname,htitle,1050,-200,4000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcEKaonSubtrNosel_%d",sl+1);
    htitle.Form("Tofino Kaon subtracted slab: %d e",sl+1);
    h_adce_kaon_subtr_nosel[sl] = new TH1F(hname,htitle,1050,-200,4000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcEBhaKaonSubtr_%d",sl+1);
    htitle.Form("Tofino Bhabha+Kaon subtracted slab: %d e",sl+1);
    h_adce_bhakaon_subtr[sl] = new TH1F(hname,htitle,1050,-200,4000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcEEqu_%d",sl+1);
    htitle.Form("Tofino equalized slab: %d e",sl+1);
    h_adce_equ[sl] = new TH1F(hname,htitle,300,-200,1000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcEPedEqu_%d",sl+1);
    htitle.Form("Tofino Pedestal equalized slab: %d e",sl+1);
    h_adce_ped_equ[sl] = new TH1F(hname,htitle,300,-200,1000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcEBhaEqu_%d",sl+1);
    htitle.Form("Tofino Bhabha equalized slab: %d e",sl+1);
    h_adce_bha_equ[sl] = new TH1F(hname,htitle,300,-200,1000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcEKaonEqu_%d",sl+1);
    htitle.Form("Tofino Kaon equalized slab: %d e",sl+1);
    h_adce_kaon_equ[sl] = new TH1F(hname,htitle,300,-200,1000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcEBhaEquNosel_%d",sl+1);
    htitle.Form("Tofino Bhabha equalized slab: %d e",sl+1);
    h_adce_bha_equ_nosel[sl] = new TH1F(hname,htitle,300,-200,1000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcEKaonEquNosel_%d",sl+1);
    htitle.Form("Tofino Kaon equalized slab: %d e",sl+1);
    h_adce_kaon_equ_nosel[sl] = new TH1F(hname,htitle,300,-200,1000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcEBhaKaonEqu_%d",sl+1);
    htitle.Form("Tofino Bhabha+Kaon equalized slab: %d e",sl+1);
    h_adce_bhakaon_equ[sl] = new TH1F(hname,htitle,300,-200,1000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcP_%d",sl+1);
    htitle.Form("Tofino slab: %d p",sl+1);
    h_adcp[sl] = new TH1F(hname,htitle,1000,0,4000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcPPed_%d",sl+1);
    htitle.Form("Tofino Pedestal slab: %d p",sl+1);
    h_adcp_ped[sl] = new TH1F(hname,htitle,1000,0,4000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcPBha_%d",sl+1);
    htitle.Form("Tofino Bhabha slab: %d p",sl+1);
    h_adcp_bha[sl] = new TH1F(hname,htitle,1000,0,4000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcPKaon_%d",sl+1);
    htitle.Form("Tofino Kaon slab: %d p",sl+1);
    h_adcp_kaon[sl] = new TH1F(hname,htitle,1000,0,4000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcPBhaNosel_%d",sl+1);
    htitle.Form("Tofino Bhabha slab: %d p",sl+1);
    h_adcp_bha_nosel[sl] = new TH1F(hname,htitle,1000,0,4000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcPKaonNosel_%d",sl+1);
    htitle.Form("Tofino Kaon slab: %d p",sl+1);
    h_adcp_kaon_nosel[sl] = new TH1F(hname,htitle,1000,0,4000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcPBhaKaon_%d",sl+1);
    htitle.Form("Tofino Bhabha+Kaon slab: %d p",sl+1);
    h_adcp_bhakaon[sl] = new TH1F(hname,htitle,1000,0,4000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcPSubtr_%d",sl+1);
    htitle.Form("Tofino subtracted slab: %d p",sl+1);
    h_adcp_subtr[sl] = new TH1F(hname,htitle,1050,-200,4000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcPPedSubtr_%d",sl+1);
    htitle.Form("Tofino Pedestal subtracted slab: %d p",sl+1);
    h_adcp_ped_subtr[sl] = new TH1F(hname,htitle,1050,-200,4000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcPBhaSubtr_%d",sl+1);
    htitle.Form("Tofino Bhabha subtracted slab: %d p",sl+1);
    h_adcp_bha_subtr[sl] = new TH1F(hname,htitle,1050,-200,4000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcPKaonSubtr_%d",sl+1);
    htitle.Form("Tofino Kaon subtracted slab: %d p",sl+1);
    h_adcp_kaon_subtr[sl] = new TH1F(hname,htitle,1050,-200,4000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcPBhaSubtrNosel_%d",sl+1);
    htitle.Form("Tofino Bhabha subtracted slab: %d p",sl+1);
    h_adcp_bha_subtr_nosel[sl] = new TH1F(hname,htitle,1050,-200,4000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcPKaonSubtrNosel_%d",sl+1);
    htitle.Form("Tofino Kaon subtracted slab: %d p",sl+1);
    h_adcp_kaon_subtr_nosel[sl] = new TH1F(hname,htitle,1050,-200,4000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcPBhaKaonSubtr_%d",sl+1);
    htitle.Form("Tofino Bhabha+Kaon subtracted slab: %d p",sl+1);
    h_adcp_bhakaon_subtr[sl] = new TH1F(hname,htitle,1050,-200,4000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcPEqu_%d",sl+1);
    htitle.Form("Tofino equalized slab: %d p",sl+1);
    h_adcp_equ[sl] = new TH1F(hname,htitle,300,-200,1000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcPPedEqu_%d",sl+1);
    htitle.Form("Tofino Pedestal equalized slab: %d p",sl+1);
    h_adcp_ped_equ[sl] = new TH1F(hname,htitle,300,-200,1000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcPBhaEqu_%d",sl+1);
    htitle.Form("Tofino Bhabha equalized slab: %d p",sl+1);
    h_adcp_bha_equ[sl] = new TH1F(hname,htitle,300,-200,1000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcPKaonEqu_%d",sl+1);
    htitle.Form("Tofino Kaon equalized slab: %d p",sl+1);
    h_adcp_kaon_equ[sl] = new TH1F(hname,htitle,300,-200,1000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcPBhaEquNosel_%d",sl+1);
    htitle.Form("Tofino Bhabha equalized slab: %d p",sl+1);
    h_adcp_bha_equ_nosel[sl] = new TH1F(hname,htitle,300,-200,1000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcPKaonEquNosel_%d",sl+1);
    htitle.Form("Tofino Kaon equalized slab: %d p",sl+1);
    h_adcp_kaon_equ_nosel[sl] = new TH1F(hname,htitle,300,-200,1000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcPBhaKaonEqu_%d",sl+1);
    htitle.Form("Tofino Bhabha+Kaon equalized slab: %d p",sl+1);
    h_adcp_bhakaon_equ[sl] = new TH1F(hname,htitle,300,-200,1000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcBhaSqr_%d",sl+1);
    htitle.Form("Tofino Bhabha equalized sqrt(e*p) slab: %d",sl+1);
    h_adc_bha_sqr[sl] = new TH1F(hname,htitle,300,-200,1000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcKaonSqr_%d",sl+1);
    htitle.Form("Tofino Kaon equalized sqrt(e*p) slab: %d",sl+1);
    h_adc_kaon_sqr[sl] = new TH1F(hname,htitle,300,-200,1000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcBhaKaonSqr_%d",sl+1);
    htitle.Form("Tofino Bhabha+Kaon equalized sqrt(e*p) slab: %d",sl+1);
    h_adc_bhakaon_sqr[sl] = new TH1F(hname,htitle,300,-200,1000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcBhaRawSqr_%d",sl+1);
    htitle.Form("Tofino Bhabha equalized RAW sqrt(e*p) slab: %d",sl+1);
    h_adc_bha_rawsqr[sl] = new TH1F(hname,htitle,300,0,3000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcKaonRawSqr_%d",sl+1);
    htitle.Form("Tofino Kaon equalized RAW sqrt(e*p) slab: %d",sl+1);
    h_adc_kaon_rawsqr[sl] = new TH1F(hname,htitle,300,0,3000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcBhaKaonRawSqr_%d",sl+1);
    htitle.Form("Tofino Bhabha+Kaon equalized RAW sqrt(e*p) slab: %d",sl+1);
    h_adc_bhakaon_rawsqr[sl] = new TH1F(hname,htitle,300,0,3000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoLowTdcE_%d",sl+1);
    htitle.Form("Low TDC Tofino slab: %d e",sl+1);    h_lowtdce[sl] = new TH1F(hname,htitle,1000,0,4000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoLowTdcP_%d",sl+1);
    htitle.Form("Low TDC Tofino slab: %d p",sl+1);
    h_lowtdcp[sl] = new TH1F(hname,htitle,1000,0,4000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoHighTdcE_%d",sl+1);
    htitle.Form("High TDC Tofino slab: %d e",sl+1);
    h_hightdce[sl] = new TH1F(hname,htitle,1000,0,4000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoHighTdcP_%d",sl+1);
    htitle.Form("High TDC Tofino slab: %d p",sl+1);
    h_hightdcp[sl] = new TH1F(hname,htitle,1000,0,4000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcEBhaEquAll_%d",sl+1);
    htitle.Form("Tofino Bhabha equalized (all) slab: %d e",sl+1);
    h_adce_bha_equ_all[sl] = new TH1F(hname,htitle,300,-200,1000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcPBhaEquAll_%d",sl+1);
    htitle.Form("Tofino Bhabha equalized (all) slab: %d p",sl+1);
    h_adcp_bha_equ_all[sl] = new TH1F(hname,htitle,300,-200,1000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcBhaSqrAll_%d",sl+1);
    htitle.Form("Tofino Bhabha equalized (all) sqrt(e*p) slab: %d",sl+1);
    h_adc_bha_sqr_all[sl] = new TH1F(hname,htitle,400,0,1600);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcEKaonEquAll_%d",sl+1);
    htitle.Form("Tofino Kaon equalized (all) slab: %d e",sl+1);
    h_adce_kaon_equ_all[sl] = new TH1F(hname,htitle,300,-200,1000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcPKaonEquAll_%d",sl+1);
    htitle.Form("Tofino Kaon equalized (all) slab: %d p",sl+1);
    h_adcp_kaon_equ_all[sl] = new TH1F(hname,htitle,300,-200,1000);
    Printf("Histo name: \"%s\"",hname.Data());

    hname.Form("TofinoAdcKaonSqrAll_%d",sl+1);
    htitle.Form("Tofino Kaon equalized (all) sqrt(e*p) slab: %d",sl+1);
    h_adc_kaon_sqr_all[sl] = new TH1F(hname,htitle,400,0,1600);
    Printf("Histo name: \"%s\"",hname.Data());

  }
}

void DeclareProfiles(Int_t const &run_num_start, Int_t const &run_num_stop){

  TString hname;
  TString htitle;

  Int_t N=run_num_stop-run_num_start+2;

  for (Int_t sl=0; sl<K_N_TOFINO_SLABS; sl++){
    hname.Form("PedE_%d",sl+1);
    htitle.Form("Tofino Pedestal slab: %d e",sl+1);
    p_ped_e[sl] = new TH1F(hname, htitle, N, run_num_start-0.5, run_num_stop+1.5);

    hname.Form("PedP_%d",sl+1);
    htitle.Form("Tofino Pedestal slab: %d p",sl+1);
    p_ped_p[sl] = new TH1F(hname, htitle, N, run_num_start-0.5, run_num_stop+1.5);

    hname.Form("BhaE_%d",sl+1);
    htitle.Form("Tofino Bhabha Peak slab: %d e",sl+1);
    p_bha_e[sl] = new TH1F(hname, htitle, N, run_num_start-0.5, run_num_stop+1.5);

    hname.Form("BhaP_%d",sl+1);
    htitle.Form("Tofino Bhabha Peak slab: %d p",sl+1);
    p_bha_p[sl] = new TH1F(hname, htitle, N, run_num_start-0.5, run_num_stop+1.5);

    hname.Form("EqBhaE_%d",sl+1);
    htitle.Form("Tofino Bhabha Peak (eq.) slab: %d e",sl+1);
    p_EQbha_e[sl] = new TH1F(hname, htitle, N, run_num_start-0.5, run_num_stop+1.5);

    hname.Form("EqBhaP_%d",sl+1);
    htitle.Form("Tofino Bhabha Peak (eq.) slab: %d p",sl+1);
    p_EQbha_p[sl] = new TH1F(hname, htitle, N, run_num_start-0.5, run_num_stop+1.5);

//      hname.Form("EqBhaSqr_%d",sl+1);
//      htitle.Form("Tofino Bhabha Peak slab sqrt(e*p) (eq.): %d",sl+1);
//      p_EQbha_sqr[sl] = new TH1F(hname, htitle, N, run_num_start-0.5, run_num_stop+1.5);

    hname.Form("KaonE_%d",sl+1);
    htitle.Form("Tofino Kaon Peak slab: %d e",sl+1);
    p_kaon_e[sl] = new TH1F(hname, htitle, N, run_num_start-0.5, run_num_stop+1.5);

    hname.Form("KaonP_%d",sl+1);
    htitle.Form("Tofino Kaon Peak slab: %d p",sl+1);
    p_kaon_p[sl] = new TH1F(hname, htitle, N, run_num_start-0.5, run_num_stop+1.5);

    hname.Form("EqKaonE_%d",sl+1);
    htitle.Form("Tofino Kaon Peak (eq.) slab: %d e",sl+1);
    p_EQkaon_e[sl] = new TH1F(hname, htitle, N, run_num_start-0.5, run_num_stop+1.5);

    hname.Form("EqKaonP_%d",sl+1);
    htitle.Form("Tofino Kaon Peak slab (eq.): %d p",sl+1);
    p_EQkaon_p[sl] = new TH1F(hname, htitle, N, run_num_start-0.5, run_num_stop+1.5);

//      hname.Form("EqKaonSqr_%d",sl+1);
//      htitle.Form("Tofino Kaon Peak slab sqrt(e*p) (eq.): %d",sl+1);
//      p_EQkaon_sqr[sl] = new TH1F(hname, htitle, N, run_num_start-0.5, run_num_stop+1.5);

    hname.Form("EqE_%d",sl+1);
    htitle.Form("Tofino Equalization Factor slab: %d e",sl+1);
    p_eq_e[sl] = new TH1F(hname, htitle, N, run_num_start-0.5, run_num_stop+1.5);

    hname.Form("EqP_%d",sl+1);
    htitle.Form("Tofino Equalization Factor slab: %d p",sl+1);
    p_eq_p[sl] = new TH1F(hname, htitle, N, run_num_start-0.5, run_num_stop+1.5);

    hname.Form("EqFactorE_%d",sl+1);
    htitle.Form("Tofino Equalization Factor slab: %d e",sl+1);
    p_EQ_e[sl] = new TH1F(hname, htitle, N, run_num_start-0.5, run_num_stop+1.5);

    hname.Form("EqFactorP_%d",sl+1);
    htitle.Form("Tofino Equalization Factor slab: %d p",sl+1);
    p_EQ_p[sl] = new TH1F(hname, htitle, N, run_num_start-0.5, run_num_stop+1.5);

    hname.Form("OverKaonThreshold_%d",sl+1);
    htitle.Form("Tofino Ratio over kaon threshold slab: %d",sl+1);
    p_OverKaonThreshold[sl] = new TH1F(hname, htitle, N, run_num_start-0.5, run_num_stop+1.5);
  }
}

void ResetHistos(){
  for(Int_t sl=0;sl<K_N_TOFINO_SLABS;sl++ ){
    h_adce[sl]->Reset();
    h_adce_ped[sl]->Reset();
    h_adce_bha[sl]->Reset();
    h_adce_kaon[sl]->Reset();
    h_adce_bha_nosel[sl]->Reset();
    h_adce_kaon_nosel[sl]->Reset();
    h_adce_bhakaon[sl]->Reset();

    h_adce_subtr[sl]->Reset();
    h_adce_ped_subtr[sl]->Reset();
    h_adce_bha_subtr[sl]->Reset();
    h_adce_kaon_subtr[sl]->Reset();
    h_adce_bha_subtr_nosel[sl]->Reset();
    h_adce_kaon_subtr_nosel[sl]->Reset();
    h_adce_bhakaon_subtr[sl]->Reset();

    h_adce_equ[sl]->Reset();
    h_adce_ped_equ[sl]->Reset();
    h_adce_bha_equ[sl]->Reset();
    h_adce_kaon_equ[sl]->Reset();
    h_adce_bha_equ_nosel[sl]->Reset();
    h_adce_kaon_equ_nosel[sl]->Reset();
    h_adce_bhakaon_equ[sl]->Reset();

    h_adcp[sl]->Reset();
    h_adcp_ped[sl]->Reset();
    h_adcp_bha[sl]->Reset();
    h_adcp_kaon[sl]->Reset();
    h_adcp_bha_nosel[sl]->Reset();
    h_adcp_kaon_nosel[sl]->Reset();
    h_adcp_bhakaon[sl]->Reset();

    h_adcp_subtr[sl]->Reset();
    h_adcp_ped_subtr[sl]->Reset();
    h_adcp_bha_subtr[sl]->Reset();
    h_adcp_kaon_subtr[sl]->Reset();
    h_adcp_bha_subtr_nosel[sl]->Reset();
    h_adcp_kaon_subtr_nosel[sl]->Reset();
    h_adcp_bhakaon_subtr[sl]->Reset();

    h_adcp_equ[sl]->Reset();
    h_adcp_ped_equ[sl]->Reset();
    h_adcp_bha_equ[sl]->Reset();
    h_adcp_kaon_equ[sl]->Reset();
    h_adcp_bha_equ_nosel[sl]->Reset();
    h_adcp_kaon_equ_nosel[sl]->Reset();
    h_adcp_bhakaon_equ[sl]->Reset();

    h_adc_bha_sqr[sl]->Reset();
    h_adc_kaon_sqr[sl]->Reset();
    h_adc_bhakaon_sqr[sl]->Reset();
    h_adc_bha_rawsqr[sl]->Reset();
    h_adc_kaon_rawsqr[sl]->Reset();
    h_adc_bhakaon_rawsqr[sl]->Reset();
  }
}


void CloseHdtFile(){
  
  if (fHdtTree){
    delete fHdtTree;
    // fHdtTree->Delete();
    fHdtTree = 0;
  }
  
  if(fHdtFile){
    if(fHdtFile->IsOpen()) fHdtFile->Close();
    delete fHdtFile;
    fHdtFile = 0;
  }
  
}

Int_t OpenHdtFile(const TString &runtype,const Int_t &run_num){
  // return value:
  //              0: ok
  //             -1: file not found

  CloseHdtFile();

  // TString runname = BuildRunName(runtype,run_num, kTRUE);
  TString runname;
  TString fnam;
  if (run_num>=7000){
    runname=BuildRunName(runtype,run_num, kTRUE);
    fnam.Form("%s/%s.hdt.root",ExpandPathName("$FND_HDT").Data(),
              runname.Data());
  } else {
    runname=BuildRunName(runtype,run_num, kFALSE);
    fnam.Form("%s_v9.0/%s.hdt.root",ExpandPathName("$FND_HDT").Data(),
              runname.Data());
  }

  
  if(gSystem->AccessPathName(fnam)){
    gROOT->Warning("OpenHdtFile","file \"%s\" not found",fnam.Data());
    return -1;
  }

  gROOT->Info("OpenHdtFile","opening file \"%s\"",fnam.Data());
  fHdtFile = new TFile(fnam,"OPEN");

  return 0;
}


Int_t LoadHdtTree(){
  // return value:
  //              0: ok
  //             -1: Problems while getting HDT tree from file
  
  if (fHdtTree){
    delete fHdtTree;
    // fHdtTree->Delete();
    fHdtTree = 0;
  }
  
  fHdtTree = (TTree *) ( fHdtFile->Get("FIN_HDT_TREE") );
  if(!fHdtTree){
    gROOT->Warning("LoadHdtTree","Problems while getting HDT tree from file");
    return -1;
  }
  //  fHdtTree->Print();

  fHdtTree->SetBranchAddress("fndhdt",&fCurHdtEv);

  return 0;
}


void FillHistosPed(){
  
  TClonesArray *TofHits = fCurHdtEv->GetTofHits();
  
  Int_t lay, slab;
  TFndHTof *cur_hit = 0;
  for(Int_t ic=0;ic< TofHits->GetEntries(); ic++) {
    cur_hit =  (TFndHTof *) ( TofHits->At(ic) );
    //    cur_hit->PrintHit();
    cur_hit->GetLayer_Slab(lay, slab);
    //     Printf("%d %d %d",lay, cham, wire);
    //     Printf("");
    if(lay==E_FIN_OUTER_LAYER) continue;

    Bool_t both_adc = kFALSE;
    Bool_t both_tdcl = kFALSE;
    Bool_t both_tdch = kFALSE;
    if(cur_hit->GetBothAdc() == TFndHit::E_HIT_BOTH_CH ) both_adc = kTRUE;
    if(cur_hit->GetBothLowTdc() == TFndHit::E_HIT_BOTH_CH ) both_tdcl = kTRUE;
    if(cur_hit->GetBothHighTdc() == TFndHit::E_HIT_BOTH_CH ) both_tdch = kTRUE;

    //     cout << "a: " << cur_hit->GetBothAdc() << " ; l: " << cur_hit->GetBothLowTdc() << " ; h: " << cur_hit->GetBothHighTdc() << endl;
    //     if(both_tdcl) continue; // pedestal
    //     if(!both_tdcl || both_tdch) continue; // bhabha
    //     if(!both_tdcl) continue; // bhabha + kaons
    //    if(!both_tdch) continue; // kaons

    Double_t adce=cur_hit->GetAdcE();
    Double_t adcp=cur_hit->GetAdcP();
    if (adce>1) h_adce[slab-1]->Fill(adce);
    if (adcp>1) h_adcp[slab-1]->Fill(adcp);
    if(!both_tdcl){ // pedestal
      if (adce>1) h_adce_ped[slab-1]->Fill(adce);
      if (adcp>1) h_adcp_ped[slab-1]->Fill(adcp);
    }   
    h_lowtdce[slab-1]->Fill(cur_hit->GetLowTdcE());
    h_lowtdcp[slab-1]->Fill(cur_hit->GetLowTdcP());
    h_hightdce[slab-1]->Fill(cur_hit->GetHighTdcE());
    h_hightdcp[slab-1]->Fill(cur_hit->GetHighTdcP());
  }
}

Int_t GetTrigType(){
  TFndTrig *cur_trig = (TFndTrig *) (fCurHdtEv->GetGts());
  Short_t hyp_bit=cur_trig->GetPU_cont_bit(0, 0);
  Short_t bha_bit=cur_trig->GetPU_cont_bit(0, 1);
  Short_t cosm_bit=cur_trig->GetPU_cont_bit(0, 5);
  if ((CurrentRunNumber>=5872&&CurrentRunNumber<=5880)||
      (CurrentRunNumber>=5887&&CurrentRunNumber<=5905)||
      (CurrentRunNumber>=5912&&CurrentRunNumber<=5935)||
      (CurrentRunNumber>=5941&&CurrentRunNumber<=5957)||
      (CurrentRunNumber>=5998&&CurrentRunNumber<=6014)||
      (CurrentRunNumber>=6020&&CurrentRunNumber<=6014)||
      (CurrentRunNumber>=7290&&CurrentRunNumber<=7316)||
      (CurrentRunNumber>=7319&&CurrentRunNumber<=7329)||
      (CurrentRunNumber>=7332&&CurrentRunNumber<=7346)||
      (CurrentRunNumber>=7349&&CurrentRunNumber<=7364)||
      (CurrentRunNumber>=7367&&CurrentRunNumber<=7382)||
      (CurrentRunNumber>=7385&&CurrentRunNumber<=7400)||
      (CurrentRunNumber>=7403&&CurrentRunNumber<=7405)||
      (CurrentRunNumber>=7403&&CurrentRunNumber<=7405)||
      (CurrentRunNumber>=7407&&CurrentRunNumber<=7421)||
      (CurrentRunNumber>=7423&&CurrentRunNumber<=7433)||
      (CurrentRunNumber>=7436&&CurrentRunNumber<=7445)||
      (CurrentRunNumber>=7447&&CurrentRunNumber<=7453)||
      (CurrentRunNumber>=7458&&CurrentRunNumber<=7463)||
      (CurrentRunNumber>=7465&&CurrentRunNumber<=7475)||
      (CurrentRunNumber>=7477&&CurrentRunNumber<=7486)||
      (CurrentRunNumber>=7488&&CurrentRunNumber<=7497)||
      (CurrentRunNumber>=7499&&CurrentRunNumber<=7517)||
      (CurrentRunNumber>=7519&&CurrentRunNumber<=7528)||
      (CurrentRunNumber>=7530&&CurrentRunNumber<=7541)||
      (CurrentRunNumber>=7543&&CurrentRunNumber<=7552)||
      (CurrentRunNumber>=7554&&CurrentRunNumber<=7563)||
      (CurrentRunNumber>=7565&&CurrentRunNumber<=7574)||
      (CurrentRunNumber>=7576&&CurrentRunNumber<=7585))      
    { hyp_bit=1; bha_bit=0; cosm_bit=0;}
  if ((CurrentRunNumber>=5881&&CurrentRunNumber<=5886)||
      (CurrentRunNumber>=5906&&CurrentRunNumber<=5911)||
      (CurrentRunNumber>=5936&&CurrentRunNumber<=5940)||
      (CurrentRunNumber>=5969&&CurrentRunNumber<=5973)||
      (CurrentRunNumber>=5986&&CurrentRunNumber<=5991)||
      (CurrentRunNumber>=7289&&CurrentRunNumber<=7289)||
      (CurrentRunNumber>=7317&&CurrentRunNumber<=7318)||
      (CurrentRunNumber>=7330&&CurrentRunNumber<=7331)||
      (CurrentRunNumber>=7347&&CurrentRunNumber<=7348)||
      (CurrentRunNumber>=7365&&CurrentRunNumber<=7366)||
      (CurrentRunNumber>=7383&&CurrentRunNumber<=7384)||
      (CurrentRunNumber>=7401&&CurrentRunNumber<=7402)||
      (CurrentRunNumber>=7406&&CurrentRunNumber<=7406)||
      (CurrentRunNumber>=7422&&CurrentRunNumber<=7422)||
      (CurrentRunNumber>=7434&&CurrentRunNumber<=7435)||
      (CurrentRunNumber>=7446&&CurrentRunNumber<=7446)||
      (CurrentRunNumber>=7454&&CurrentRunNumber<=7457)||
      (CurrentRunNumber>=7464&&CurrentRunNumber<=7464)||
      CurrentRunNumber==7476||CurrentRunNumber==7487||
      CurrentRunNumber==7498||CurrentRunNumber==7518||
      CurrentRunNumber==7529||CurrentRunNumber==7542||
      CurrentRunNumber==7553||CurrentRunNumber==7564||
      CurrentRunNumber==7575||CurrentRunNumber==7586)
    { hyp_bit=0; bha_bit=1; cosm_bit=0;}
  if ((CurrentRunNumber>=7587&&CurrentRunNumber<=8335)||
      (CurrentRunNumber>=9027&&CurrentRunNumber<=9084)||
      (CurrentRunNumber>=9394&&CurrentRunNumber<=9585))
  {
    if (CurrentRunNumber==7597||
        CurrentRunNumber==7608||
        CurrentRunNumber==7619||
        CurrentRunNumber==7630||
        CurrentRunNumber==7639||
        CurrentRunNumber==7650||
        CurrentRunNumber==7661||
        CurrentRunNumber==7670||
        CurrentRunNumber==7690||
        CurrentRunNumber==7691||
        CurrentRunNumber==7702||
        CurrentRunNumber==7713||
        CurrentRunNumber==7719||
        CurrentRunNumber==7733||
        CurrentRunNumber==7744||
        CurrentRunNumber==7756||
        CurrentRunNumber==7777||
        CurrentRunNumber==7788||
        CurrentRunNumber==7792||
        CurrentRunNumber==7793||
        (CurrentRunNumber>=7799&&CurrentRunNumber<=7802)||
        CurrentRunNumber==7814||
        CurrentRunNumber==7825||
        CurrentRunNumber==7836||
        CurrentRunNumber==7847||
        CurrentRunNumber==7858||
        CurrentRunNumber==7869||
        CurrentRunNumber==7873||
        CurrentRunNumber==7874||
        CurrentRunNumber==7885||
        CurrentRunNumber==7896||
        CurrentRunNumber==7907||
        CurrentRunNumber==7919||
        CurrentRunNumber==7921||
        CurrentRunNumber==7932||
        CurrentRunNumber==7943||
        CurrentRunNumber==7954||
        CurrentRunNumber==7965||
        CurrentRunNumber==7976||
        CurrentRunNumber==7987||
        CurrentRunNumber==7991||
        CurrentRunNumber==7992||
        CurrentRunNumber==7993||
        CurrentRunNumber==8004||
        CurrentRunNumber==8005||
        CurrentRunNumber==8023||
        CurrentRunNumber==8064||
        CurrentRunNumber==8065||
        CurrentRunNumber==8074||
        CurrentRunNumber==8085||
        CurrentRunNumber==8098||
        CurrentRunNumber==8109||
        CurrentRunNumber==8120||
        CurrentRunNumber==8131||
        CurrentRunNumber==8142||
        CurrentRunNumber==8153||
        CurrentRunNumber==8163||
        CurrentRunNumber==8174||
        CurrentRunNumber==8185||
        CurrentRunNumber==8208||
        CurrentRunNumber==8222||
        CurrentRunNumber==8233||
        CurrentRunNumber==8244||
        CurrentRunNumber==8255||
        CurrentRunNumber==8266||
        CurrentRunNumber==8277||
        CurrentRunNumber==8288||
        CurrentRunNumber==8299||
        CurrentRunNumber==8303||
        CurrentRunNumber==8308||
        CurrentRunNumber==8319||
        CurrentRunNumber==8334||
        CurrentRunNumber==8335||
        CurrentRunNumber==9029||
        CurrentRunNumber==9040||
        CurrentRunNumber==9051||
        CurrentRunNumber==9062||
        CurrentRunNumber==9073||
        CurrentRunNumber==9401||
        CurrentRunNumber==9412||
        CurrentRunNumber==9416||
        CurrentRunNumber==9428||
        CurrentRunNumber==9439||
        CurrentRunNumber==9450||
        CurrentRunNumber==9461||
        CurrentRunNumber==9472||
        CurrentRunNumber==9483||
        CurrentRunNumber==9494||
        CurrentRunNumber==9505||
        CurrentRunNumber==9516||
        CurrentRunNumber==9527||
        CurrentRunNumber==9538||
        CurrentRunNumber==9549||
        CurrentRunNumber==9560||
        CurrentRunNumber==9571||
        CurrentRunNumber==9582)
      { hyp_bit=0; bha_bit=1; cosm_bit=0;}
    else { hyp_bit=1; bha_bit=0; cosm_bit=0;}
  }
  return (Int_t) hyp_bit+2*bha_bit+32*cosm_bit;
}

void GetMultiplicity(Int_t *mult){
  // mult[0] ... low threshold
  // mult[1] ... high threshold
  mult[0]=0; mult[1]=0;

  TClonesArray *TofHits = fCurHdtEv->GetTofHits();
  Int_t lay, slab;
  TFndHTof *cur_hit = 0;

  for(Int_t ic=0;ic< TofHits->GetEntries(); ic++) {
    cur_hit =  (TFndHTof *) ( TofHits->At(ic) );
    cur_hit->GetLayer_Slab(lay, slab);
    if(lay==E_FIN_OUTER_LAYER) continue;

    Bool_t both_adc = kFALSE;
    Bool_t both_tdcl = kFALSE;
    Bool_t both_tdch = kFALSE;
    if(cur_hit->GetBothAdc() == TFndHit::E_HIT_BOTH_CH ) both_adc = kTRUE;
    if(cur_hit->GetBothLowTdc() == TFndHit::E_HIT_BOTH_CH ) both_tdcl = kTRUE;
    if(cur_hit->GetBothHighTdc() == TFndHit::E_HIT_BOTH_CH ) both_tdch = kTRUE;

    Double_t adce=cur_hit->GetAdcE();
    Double_t adcp=cur_hit->GetAdcP();

    if(both_tdch&&adce>0&&adcp>0){ // kaons
      mult[1]++;
    }   
    if(both_tdcl&&adce>0&&adcp>0){ // bhabha+kaons
      mult[0]++;
    }
  }
}
  

void FillHistos(){
  
  TClonesArray *TofHits = fCurHdtEv->GetTofHits();
  
  Int_t lay, slab;
  TFndHTof *cur_hit = 0;
  Int_t mult[2];
  GetMultiplicity(mult);
  for(Int_t ic=0;ic< TofHits->GetEntries(); ic++) {
    cur_hit =  (TFndHTof *) ( TofHits->At(ic) );
    //    cur_hit->PrintHit();
    cur_hit->GetLayer_Slab(lay, slab);
    //     Printf("%d %d %d",lay, cham, wire);
    //     Printf("");
    if(lay==E_FIN_OUTER_LAYER) continue;

    Bool_t both_adc = kFALSE;
    Bool_t both_tdcl = kFALSE;
    Bool_t both_tdch = kFALSE;
    if(cur_hit->GetBothAdc() == TFndHit::E_HIT_BOTH_CH ) both_adc = kTRUE;
    if(cur_hit->GetBothLowTdc() == TFndHit::E_HIT_BOTH_CH ) both_tdcl = kTRUE;
    if(cur_hit->GetBothHighTdc() == TFndHit::E_HIT_BOTH_CH ) both_tdch = kTRUE;

    //     cout << "a: " << cur_hit->GetBothAdc() << " ; l: " << cur_hit->GetBothLowTdc() << " ; h: " << cur_hit->GetBothHighTdc() << endl;
    //     if(both_tdcl) continue; // pedestal
    //     if(!both_tdcl || both_tdch) continue; // bhabha
    //     if(!both_tdcl) continue; // bhabha + kaons
    //    if(!both_tdch) continue; // kaons

    Double_t adce=cur_hit->GetAdcE();
    Double_t adcp=cur_hit->GetAdcP();

    if (adce>1) h_adce_subtr[slab-1]->Fill(adce-ped_e[slab-1]);
    if (adcp>1) h_adcp_subtr[slab-1]->Fill(adcp-ped_p[slab-1]);
    if(!both_tdcl){ // pedestal
      if (adce>1)
        h_adce_ped_subtr[slab-1]->Fill(adce-ped_e[slab-1]);
      if (adcp>1)
        h_adcp_ped_subtr[slab-1]->Fill(adcp-ped_p[slab-1]);
    }   
    if(both_tdcl && !both_tdch){ // bhabha 
      if (adce>1){
        h_adce_bha_nosel[slab-1]->Fill(adce);
        h_adce_bha_subtr_nosel[slab-1]->Fill(adce-ped_e[slab-1]);
        if ((GetTrigType()==2||GetTrigType()==32)&&mult[0]==2){
          h_adce_bha[slab-1]->Fill(adce);
          h_adce_bha_subtr[slab-1]->Fill(adce-ped_e[slab-1]);
        }
      }
      if (adcp>1){
        h_adcp_bha_nosel[slab-1]->Fill(adcp);
        h_adcp_bha_subtr_nosel[slab-1]->Fill(adcp-ped_p[slab-1]);
        if ((GetTrigType()==2||GetTrigType()==32)&&mult[0]==2){
          h_adcp_bha[slab-1]->Fill(adcp);
          h_adcp_bha_subtr[slab-1]->Fill(adcp-ped_p[slab-1]);
        }
      }
    }   
    if(both_tdch){ // kaons
      if (adce>1){
        h_adce_kaon_nosel[slab-1]->Fill(adce);
        h_adce_kaon_subtr_nosel[slab-1]->Fill(adce-ped_e[slab-1]);
        if (GetTrigType()==1&&mult[1]==2){
          h_adce_kaon[slab-1]->Fill(adce);
          h_adce_kaon_subtr[slab-1]->Fill(adce-ped_e[slab-1]);
        }
      }
      if (adcp>1){
        h_adcp_kaon_nosel[slab-1]->Fill(adcp);
        h_adcp_kaon_subtr_nosel[slab-1]->Fill(adcp-ped_p[slab-1]);
        if (GetTrigType()==1&&mult[1]==2){
          h_adcp_kaon[slab-1]->Fill(adcp);
          h_adcp_kaon_subtr[slab-1]->Fill(adcp-ped_p[slab-1]);
        }
      }
    }   
    if(both_tdcl){ // bhabha+kaons
      if (adce>1){
        h_adce_bhakaon[slab-1]->Fill(adce);
        h_adce_bhakaon_subtr[slab-1]->Fill(adce-ped_e[slab-1]);
      }
      if (adcp>1){
        h_adcp_bhakaon[slab-1]->Fill(adcp);
        h_adcp_bhakaon_subtr[slab-1]->Fill(adcp-ped_p[slab-1]);
      }
    }
  }
  if (GetTrigType()==2){
      h_LowBhaMult->Fill(mult[0]);
      h_HighBhaMult->Fill(mult[1]);
  }
  if (GetTrigType()==1){
    h_LowHypMult->Fill(mult[0]);
    h_HighHypMult->Fill(mult[1]);
  }
}

void FillHistosEqu(){
  
  TClonesArray *TofHits = fCurHdtEv->GetTofHits();
  
  Int_t lay, slab;
  TFndHTof *cur_hit = 0;
  Int_t mult[2];
  GetMultiplicity(mult);
  for(Int_t ic=0;ic< TofHits->GetEntries(); ic++) {
    cur_hit =  (TFndHTof *) ( TofHits->At(ic) );
    //    cur_hit->PrintHit();
    cur_hit->GetLayer_Slab(lay, slab);
    //     Printf("%d %d %d",lay, cham, wire);
    //     Printf("");
    if(lay==E_FIN_OUTER_LAYER) continue;

    Bool_t both_adc = kFALSE;
    Bool_t both_tdcl = kFALSE;
    Bool_t both_tdch = kFALSE;
    if(cur_hit->GetBothAdc() == TFndHit::E_HIT_BOTH_CH ) both_adc = kTRUE;
    if(cur_hit->GetBothLowTdc() == TFndHit::E_HIT_BOTH_CH ) both_tdcl = kTRUE;
    if(cur_hit->GetBothHighTdc() == TFndHit::E_HIT_BOTH_CH ) both_tdch = kTRUE;

    //     cout << "a: " << cur_hit->GetBothAdc() << " ; l: " << cur_hit->GetBothLowTdc() << " ; h: " << cur_hit->GetBothHighTdc() << endl;
    //     if(both_tdcl) continue; // pedestal
    //     if(!both_tdcl || both_tdch) continue; // bhabha
    //     if(!both_tdcl) continue; // bhabha + kaons
    //    if(!both_tdch) continue; // kaons

    Double_t adce=cur_hit->GetAdcE();
    Double_t adcp=cur_hit->GetAdcP();
    if (adce>1)
      h_adce_equ[slab-1]->Fill((adce-ped_e[slab-1])*eq_e[slab-1]);
    if (adcp>1)
      h_adcp_equ[slab-1]->Fill((adcp-ped_p[slab-1])*eq_p[slab-1]);
     if(!both_tdcl){ // pedestal
      if (adce>1)
        h_adce_ped_equ[slab-1]->Fill((adce-ped_e[slab-1])*eq_e[slab-1]);
      if (adcp>1)
        h_adcp_ped_equ[slab-1]->Fill((adcp-ped_p[slab-1])*eq_p[slab-1]);
    }   
    if(both_tdcl && !both_tdch){ // bhabha 
      if (adce>1)
        h_adce_bha_equ_nosel[slab-1]->Fill((adce-ped_e[slab-1])*eq_e[slab-1]);
      if (adcp>1)
        h_adcp_bha_equ_nosel[slab-1]->Fill((adcp-ped_p[slab-1])*eq_p[slab-1]);
      if ((GetTrigType()==2||GetTrigType()==32)&&mult[0]==2){
        if (adce>1)
          h_adce_bha_equ[slab-1]->Fill((adce-ped_e[slab-1])*eq_e[slab-1]);
        if (adcp>1)
          h_adcp_bha_equ[slab-1]->Fill((adcp-ped_p[slab-1])*eq_p[slab-1]);
        if (adce>1&&adcp>1){
          h_adc_bha_sqr[slab-1]->Fill(sqrt((adce-ped_e[slab-1])*eq_e[slab-1]
                                      *(adcp-ped_p[slab-1])*eq_p[slab-1]));
          h_adc_bha_rawsqr[slab-1]->Fill(sqrt(adce*adcp));
        }
      }   
    }
    if(both_tdch){ // kaons
      if (adce>1)
        h_adce_kaon_equ_nosel[slab-1]->Fill((adce-ped_e[slab-1])*eq_e[slab-1]);
      if (adcp>1)
        h_adcp_kaon_equ_nosel[slab-1]->Fill((adcp-ped_p[slab-1])*eq_p[slab-1]);
      if (GetTrigType()==1&&mult[1]==2){
        if (adce>1)
          h_adce_kaon_equ[slab-1]->Fill((adce-ped_e[slab-1])*eq_e[slab-1]);
        if (adcp>1)
          h_adcp_kaon_equ[slab-1]->Fill((adcp-ped_p[slab-1])*eq_p[slab-1]);
        if (adce>1&&adcp>1){
          h_adc_kaon_sqr[slab-1]->Fill(sqrt((adce-ped_e[slab-1])*eq_e[slab-1]
                                       *(adcp-ped_p[slab-1])*eq_p[slab-1]));
          h_adc_kaon_rawsqr[slab-1]->Fill(sqrt(adce*adcp));
        }
      }   
    }
    if(both_tdcl){ // bhabha+kaons
      if (adce>1)
        h_adce_bhakaon_equ[slab-1]->Fill((adce-ped_e[slab-1])*eq_e[slab-1]);
      if (adcp>1)
        h_adcp_bhakaon_equ[slab-1]->Fill((adcp-ped_p[slab-1])*eq_p[slab-1]);
      if (adce>1&&adcp>1){
        h_adc_bhakaon_sqr[slab-1]->Fill(sqrt((adce-ped_e[slab-1])*eq_e[slab-1]
                                        *(adcp-ped_p[slab-1])*eq_p[slab-1]));
        h_adc_bhakaon_rawsqr[slab-1]->Fill(sqrt(adce*adcp));
      }
    }   
  }
}

void FillHistosAll(){
  
  TClonesArray *TofHits = fCurHdtEv->GetTofHits();
  
  Int_t lay, slab;
  TFndHTof *cur_hit = 0;
  Int_t mult[2];
  GetMultiplicity(mult);
  for(Int_t ic=0;ic< TofHits->GetEntries(); ic++) {
    cur_hit =  (TFndHTof *) ( TofHits->At(ic) );
    //    cur_hit->PrintHit();
    cur_hit->GetLayer_Slab(lay, slab);
    //     Printf("%d %d %d",lay, cham, wire);
    //     Printf("");
    if(lay==E_FIN_OUTER_LAYER) continue;

    Bool_t both_adc = kFALSE;
    Bool_t both_tdcl = kFALSE;
    Bool_t both_tdch = kFALSE;
    if(cur_hit->GetBothAdc() == TFndHit::E_HIT_BOTH_CH ) both_adc = kTRUE;
    if(cur_hit->GetBothLowTdc() == TFndHit::E_HIT_BOTH_CH ) both_tdcl = kTRUE;
    if(cur_hit->GetBothHighTdc() == TFndHit::E_HIT_BOTH_CH ) both_tdch = kTRUE;

    //     cout << "a: " << cur_hit->GetBothAdc() << " ; l: " << cur_hit->GetBothLowTdc() << " ; h: " << cur_hit->GetBothHighTdc() << endl;
    //     if(both_tdcl) continue; // pedestal
    //     if(!both_tdcl || both_tdch) continue; // bhabha
    //     if(!both_tdcl) continue; // bhabha + kaons
    //    if(!both_tdch) continue; // kaons

    Double_t adce=cur_hit->GetAdcE();
    Double_t adcp=cur_hit->GetAdcP();
    Double_t sqradc;
    sqradc=sqrt((adce-ped_table_e[slab-1])*eq_table_e[slab-1]
                *(adcp-ped_table_p[slab-1])*eq_table_p[slab-1]);

    if(both_tdcl && !both_tdch){ // bhabha 
      if ((GetTrigType()==2||GetTrigType()==32)&&mult[0]==2){
        if (adce>1)
          h_adce_bha_equ_all[slab-1]->Fill((adce-ped_table_e[slab-1])*eq_table_e[slab-1]);
        if (adcp>1)
          h_adcp_bha_equ_all[slab-1]->Fill((adce-ped_table_p[slab-1])*eq_table_p[slab-1]);
        if (adce>1&&adcp>1)
          h_adc_bha_sqr_all[slab-1]->Fill(sqradc);
      }  
    }

    if(both_tdch){ // kaons
      if (GetTrigType()==1&&mult[1]==2){
        if (adce>1)
          h_adce_kaon_equ_all[slab-1]->Fill((adce-ped_table_e[slab-1])*eq_table_e[slab-1]);
        if (adcp>1)
          h_adcp_kaon_equ_all[slab-1]->Fill((adcp-ped_table_p[slab-1])*eq_table_p[slab-1]);
        if (adce>1&&adcp>1){
          h_adc_kaon_sqr_all[slab-1]->Fill(sqradc);
          if (sqradc>KaonThreshold[slab-1]) NevOverThreshold[slab-1]++;
        }
      }   
    }
  }
}

Double_t FindPedestal(TH1F* histo){
  Int_t min_bin = 2;
  Int_t max_bin = 100;
  Double_t max_entry = -1.;
  Int_t max_entry_bin = -1;
  for (Int_t bin=min_bin; bin<max_bin; bin++){
    Double_t entry=histo->GetBinContent(bin);
    if (max_entry<entry){
      max_entry=entry;
      max_entry_bin=bin;
    }
  }
  return histo->GetBinCenter(max_entry_bin);
}

Double_t FindBhaPeak(TH1F* histo, Int_t B){
  // Int_t min_ch = -50;
  // Int_t max_ch = 150;
  //  histo->Fit("gaus","","",min_ch,max_ch);
  // TF1 *gaus=(TF1*)gROOT->FindObject("gaus");
  // return gaus->GetParameter(1);
  Double_t fr[2];
  Double_t sv[4], pllo[4], plhi[4], fp[4], fpe[4];

//    if (B==1){
    fr[0]=-20.;
    fr[1]=400.;
//    } else {
//      fr[0]=0.;
//      fr[1]=1500.;
//    }
    pllo[0]=1.0; pllo[1]=10.; pllo[2]=1.0; pllo[3]=0.1;
//    if (B==1){
    plhi[0]=100.0; plhi[1]=150.; plhi[2]=100000000.0; plhi[3]=30.0;
    sv[0]=1.0; sv[1]=50.0; sv[2]=50000.0; sv[3]=4.0;
//    } else {
//      plhi[0]=1000.0; plhi[1]=1000.; plhi[2]=100000000.0; plhi[3]=100.0;
//      sv[0]=400.0; sv[1]=600.0; sv[2]=50000.0; sv[3]=4.0;
    //  }

   Double_t chisqr;
   Int_t    ndf;
   if (B==1){
     TF1 *fitsnr = langaufit(histo,fr,sv,pllo,plhi,fp,fpe,&chisqr,&ndf);  
     fitsnr->Draw("lsame");
     delete fitsnr;
     return fp[1];
   } else {
     histo->Fit("landau","Q0","",0,3000);
     TF1 *landau=(TF1*)gROOT->FindObject("landau");
     return landau->GetParameter(1);
   }
}

Double_t FindKaonPeak(TH1F* histo){
  cout << "histo:Entry" << histo->GetEntries() << endl;
  if (histo->GetEntries()<100) {return -1;}

  TH1F *histo_clone_old=(TH1F*)gROOT->FindObject("histo_clone");
  if (histo_clone_old) delete histo_clone_old;
  TH1F *histo_clone=(TH1F*)histo->Clone();
  histo_clone->SetName("histo_clone");
  cout << "pointer " << histo << " " << histo_clone << endl;
  histo_clone->Rebin(5);

  histo->Fit("gaus","Q0","",300,1000);
  TF1 *gaus=(TF1*)gROOT->FindObject("gaus");
  Double_t center0=0., center2=0.;
  Double_t center=gaus->GetParameter(1);
  cout << "GAUS" << gaus << endl;

  Int_t count=0;
  TF1* landau=(TF1*)gROOT->FindObject("landau");
  do {
    if (count>0) center2=center0;
    center0=center;
    histo_clone->Fit("landau","Q0","",center*0.75,1000);
    center=landau->GetParameter(1);
    cout << center2 << " --> " << center0 << " -- > " << center << endl;
    if (center2==center){
      cout << "the fit result is oscillating..." << endl;
      center=(center+center0)/2.;
    }
    if (center>1000) center=1000;
    count++;
  } while (fabs(center-center0)>0.1&&count<20);
  return center;
}

void EvaluatePed(){
  for (Int_t sl=0; sl<K_N_TOFINO_SLABS; sl++){
    ped_e[sl]=FindPedestal(h_adce_ped[sl]);
    ped_p[sl]=FindPedestal(h_adcp_ped[sl]);
  }
}

void EvaluateEqu(Int_t B){
  for (Int_t sl=0; sl<K_N_TOFINO_SLABS; sl++){
    if (B==1){
      if (h_adce_bha_subtr[sl]->GetEntries()>500){
        bhapeak_e[sl]=FindBhaPeak(h_adce_bha_subtr[sl],B);
      } else {
        bhapeak_e[sl]=FindBhaPeak(h_adce_bha_subtr_nosel[sl],B);
      }
      if (h_adcp_bha_subtr[sl]->GetEntries()>500){
        bhapeak_p[sl]=FindBhaPeak(h_adcp_bha_subtr[sl],B);
      } else {
        bhapeak_p[sl]=FindBhaPeak(h_adcp_bha_subtr_nosel[sl],B);
      }

      
    } else {
      bhapeak_e[sl]=FindBhaPeak(h_adce_kaon_subtr_nosel[sl],B);
      bhapeak_p[sl]=FindBhaPeak(h_adcp_kaon_subtr_nosel[sl],B);
    }    
  }
  for (Int_t sl=0; sl<K_N_TOFINO_SLABS; sl++){
    eq_e[sl]= 50.0/bhapeak_e[sl];
    eq_p[sl]= 50.0/bhapeak_p[sl];
  }
}

void EvaluateKaonPeak(){
  for (Int_t sl=0; sl<K_N_TOFINO_SLABS; sl++){
    cout << h_adce_kaon_equ[sl]->GetEntries() << endl;
      kaonpeak_e[sl]=FindKaonPeak(h_adce_kaon_equ[sl]);
      kaonpeak_p[sl]=FindKaonPeak(h_adcp_kaon_equ[sl]);
  }
}
void EventLoopPed(const Int_t &nevs){
 
  Int_t entries = fHdtTree->GetEntries();
  
  Int_t cur_ev = 0;
  while(1){
    if(cur_ev >= entries){      
      gROOT->Info("EventLoop","run completed");
      break;
    }
    if(nevs != -1 && cur_ev >= nevs){      
      gROOT->Info("EventLoop","requested number of events reached");
      break;
    }
    //
    fHdtTree->GetEntry(cur_ev);
    //
    FillHistosPed();
    if(cur_ev%3000 == 0 ) cout << "cur_ev: " << cur_ev << endl;
    //     fCurHdtEv->PrintTofHits();
    //
    cur_ev++;
  }
}

void EventLoop(const Int_t &nevs){
  
  Int_t entries = fHdtTree->GetEntries();
  
  Int_t cur_ev = 0;
  while(1){
    if(cur_ev >= entries){      
      gROOT->Info("EventLoop","run completed");
      break;
    }
    if(nevs != -1 && cur_ev >= nevs){      
      gROOT->Info("EventLoop","requested number of events reached");
      break;
    }
    //
    fHdtTree->GetEntry(cur_ev);
    //
    FillHistos();
    if(cur_ev%3000 == 0 ) cout << "cur_ev: " << cur_ev << endl;
    //     fCurHdtEv->PrintTofHits();
    //
    cur_ev++;
  }
}

void EventLoopEqu(const Int_t &nevs){
  
  Int_t entries = fHdtTree->GetEntries();
  
  Int_t cur_ev = 0;
  while(1){
    if(cur_ev >= entries){      
      gROOT->Info("EventLoop","run completed");
      break;
    }
    if(nevs != -1 && cur_ev >= nevs){      
      gROOT->Info("EventLoop","requested number of events reached");
      break;
    }
    //
    fHdtTree->GetEntry(cur_ev);
    //
    FillHistosEqu();
    if(cur_ev%3000 == 0 ) cout << "cur_ev: " << cur_ev << endl;
    //     fCurHdtEv->PrintTofHits();
    //
    cur_ev++;
  }
}

void EventLoopEquAll(const Int_t &nevs){
  
  Int_t entries = fHdtTree->GetEntries();
  
  Int_t cur_ev = 0;
  while(1){
    if(cur_ev >= entries){      
      gROOT->Info("EventLoop","run completed");
      break;
    }
    if(nevs != -1 && cur_ev >= nevs){      
      gROOT->Info("EventLoop","requested number of events reached");
      break;
    }
    //
    fHdtTree->GetEntry(cur_ev);
    //
    FillHistosEqu();
    FillHistosAll();
    if(cur_ev%3000 == 0 ) cout << "cur_ev: " << cur_ev << endl;
    //     fCurHdtEv->PrintTofHits();
    //
    cur_ev++;
  }
}

void SaveHistos(const TString &fname){
  TFile *fOutFile = new TFile(fname,"RECREATE");
  h_LowBhaMult->Write();
  h_LowHypMult->Write();
  h_HighBhaMult->Write();
  h_HighHypMult->Write();

  for(Int_t sl=0;sl<K_N_TOFINO_SLABS;sl++ ){
    h_adce[sl]->Write();
    h_adce_ped[sl]->Write();
    h_adce_bha[sl]->Write();
    h_adce_kaon[sl]->Write();
    h_adce_bha_nosel[sl]->Write();
    h_adce_kaon_nosel[sl]->Write();
    h_adce_bhakaon[sl]->Write();
    h_adce_subtr[sl]->Write();
    h_adce_ped_subtr[sl]->Write();
    h_adce_bha_subtr[sl]->Write();
    h_adce_kaon_subtr[sl]->Write();
    h_adce_bha_subtr_nosel[sl]->Write();
    h_adce_kaon_subtr_nosel[sl]->Write();
    h_adce_bhakaon_subtr[sl]->Write();
    h_adce_equ[sl]->Write();
    h_adce_ped_equ[sl]->Write();
    h_adce_bha_equ[sl]->Write();
    h_adce_kaon_equ[sl]->Write();
    h_adce_bha_equ_nosel[sl]->Write();
    h_adce_kaon_equ_nosel[sl]->Write();
    h_adce_bhakaon_equ[sl]->Write();
    h_adcp[sl]->Write();
    h_adcp_ped[sl]->Write();
    h_adcp_bha[sl]->Write();
    h_adcp_kaon[sl]->Write();
    h_adcp_bha_nosel[sl]->Write();
    h_adcp_kaon_nosel[sl]->Write();
    h_adcp_bhakaon[sl]->Write();
    h_adcp_subtr[sl]->Write();
    h_adcp_ped_subtr[sl]->Write();
    h_adcp_bha_subtr[sl]->Write();
    h_adcp_kaon_subtr[sl]->Write();
    h_adcp_bha_subtr_nosel[sl]->Write();
    h_adcp_kaon_subtr_nosel[sl]->Write();
    h_adcp_bhakaon_subtr[sl]->Write();
    h_adcp_equ[sl]->Write();
    h_adcp_ped_equ[sl]->Write();
    h_adcp_bha_equ[sl]->Write();
    h_adcp_kaon_equ[sl]->Write();
    h_adcp_bha_equ_nosel[sl]->Write();
    h_adcp_kaon_equ_nosel[sl]->Write();
    h_adcp_bhakaon_equ[sl]->Write();
    h_adc_bha_sqr[sl]->Write();
    h_adc_kaon_sqr[sl]->Write();
    h_adc_bhakaon_sqr[sl]->Write();
    h_adc_bha_rawsqr[sl]->Write();
    h_adc_kaon_rawsqr[sl]->Write();
    h_adc_bhakaon_rawsqr[sl]->Write();
    h_lowtdce[sl]->Write();
    h_lowtdcp[sl]->Write();
    h_hightdce[sl]->Write();
    h_hightdcp[sl]->Write();
   }
  
//   TofCalibration->Write();
//   cout << "================> " << TofCalibration->GetPedE(1) << endl;
//   cout << "================> " << TofCalibration->GetPedE(4) << endl;
  fCalibTree->Fill();
  fCalibTree->Write();

  fOutFile->Close();
  delete fOutFile;
  fOutFile = 0;

}

void SaveProfiles(const TString &fname){
  TFile *fOutFile = new TFile(fname,"RECREATE");

  for (Int_t sl=0; sl<K_N_TOFINO_SLABS; sl++){
    p_ped_e[sl]->Write();
    p_ped_p[sl]->Write();
    p_bha_e[sl]->Write();
    p_bha_p[sl]->Write();
    p_kaon_e[sl]->Write();
    p_kaon_p[sl]->Write(); 
    p_eq_e[sl]->Write();
    p_eq_p[sl]->Write();
    p_EQbha_e[sl]->Write();
    p_EQbha_p[sl]->Write();
//      p_EQbha_sqr[sl]->Write();
    p_EQkaon_e[sl]->Write();
    p_EQkaon_p[sl]->Write();
//      p_EQkaon_sqr[sl]->Write();
    p_EQ_e[sl]->Write();
    p_EQ_p[sl]->Write();
    p_OverKaonThreshold[sl]->Write();

    h_adce_bha_equ_all[sl]->Write();
    h_adcp_bha_equ_all[sl]->Write();
    h_adc_bha_sqr_all[sl]->Write();
    h_adce_kaon_equ_all[sl]->Write();
    h_adcp_kaon_equ_all[sl]->Write();
    h_adc_kaon_sqr_all[sl]->Write();
  }

  fOutFile->Close();
  delete fOutFile;
  fOutFile = 0;

}

Int_t GetRunInfo(TString runtype,Int_t runnum){

  TString run_path;
  if (runnum>=7000){
    run_path.Form("%s/%s.hdt.root",ExpandPathName("$FND_HDT").Data(),
                  FROOT::BuildRunName(runtype,runnum,kTRUE).Data());
  } else {
    run_path.Form("%s_v9.0/%s.hdt.root",ExpandPathName("$FND_HDT").Data(),
                  FROOT::BuildRunName(runtype,runnum,kFALSE).Data());
  }

  TFile *hdt_file = new TFile(run_path);
  if (!hdt_file) return -1;

  //  hdt_file->ls();
  TFndGenInfo *geninfo = (TFndGenInfo *)hdt_file->Get("TFndGenInfo");
  
  //  geninfo->Print();
  
  TDatime d = geninfo->GetDateTime();
  
  return (Int_t) d.Convert();

}

void FillTofFromHdt(const TString &runtype,const Int_t &run_num_start,const Int_t &run_num_stop,const TString &root_filename,const Int_t &nevs, Int_t B){
  
  gStyle->SetOptFit(1);
//   TofCalibration = new TofiCalib();
  Int_t start_run=run_num_start, stop_run=run_num_stop;
  fCalibTree = new TTree("CalibTree","Calibreatio Tree");
  fCalibTree->Branch("StartRun",&start_run,"/I");
  fCalibTree->Branch("StartTime",&start_unixtime,"/I");
  fCalibTree->Branch("StopRun",&stop_run,"/I");
  fCalibTree->Branch("StopTime",&stop_unixtime,"/I");
  fCalibTree->Branch("PedE",ped_e,"[12]/F");
  fCalibTree->Branch("PedP",ped_p,"[12]/F");
  fCalibTree->Branch("EqE",eq_e,"[12]/F");
  fCalibTree->Branch("EqP",eq_p,"[12]/F");

  if (root_filename!=""){
    start_unixtime=GetRunInfo(runtype, run_num_start);
    //stop_unixtime=GetRunInfo(runtype, run_num_stop+1);
    stop_unixtime=GetRunInfo(runtype, run_num_stop);
    
    DeclareHistos();
  }

  // --- loop in runs
  for(Int_t cur_run=run_num_start; cur_run <=run_num_stop ; cur_run++){
    cout << "Filling from run \"" << runtype.Data() << "\" " << cur_run << endl;
    CurrentRunNumber = cur_run;
    if(OpenHdtFile(runtype,cur_run) !=0) continue; 
    if(LoadHdtTree() != 0) continue;
    
    EventLoopPed(nevs);
    CloseHdtFile();
  }
  // --- loop in runs
  EvaluatePed();

  // --- loop in runs
  for(Int_t cur_run=run_num_start; cur_run <=run_num_stop ; cur_run++){
    cout << "Filling from run \"" << runtype.Data() << "\" " << cur_run << endl;
    CurrentRunNumber = cur_run;
    if(OpenHdtFile(runtype,cur_run) !=0) continue; 
    if(LoadHdtTree() != 0) continue;
    
    EventLoop(nevs);
    CloseHdtFile();
  }
  // --- loop in runs
  EvaluateEqu(B);

  // --- loop in runs
  for(Int_t cur_run=run_num_start; cur_run <=run_num_stop ; cur_run++){
    cout << "Filling from run \"" << runtype.Data() << "\" " << cur_run << endl;
    CurrentRunNumber = cur_run;
    if(OpenHdtFile(runtype,cur_run) !=0) continue; 
    if(LoadHdtTree() != 0) continue;
    
    if (root_filename!=""){
      EventLoopEqu(nevs);}
    else {
      for (Int_t sl=0; sl<K_N_TOFINO_SLABS; sl++){
        NevOverThreshold[sl]=0;
      }
      EventLoopEquAll(nevs);
      for (Int_t sl=0; sl<K_N_TOFINO_SLABS; sl++){
        cout << NevOverThreshold[sl]  << "***" << h_adc_kaon_sqr[sl]->GetEntries() << endl;
        if (h_adc_kaon_sqr[sl]->GetEntries()>0){
          p_OverKaonThreshold[sl]->Fill(cur_run, NevOverThreshold[sl]/h_adc_kaon_sqr[sl]->GetEntries());
        }
      }
    }
    CloseHdtFile();
  }
  EvaluateKaonPeak();

  // --- loop in runs

  // Store information  
  //   TofCalibration->SetPedE(4,998);
  //   cout << "================> " << TofCalibration->GetPedE(4) << endl;
  // --
  if (root_filename!="") SaveHistos(root_filename);

}

void FillTofFromHdt(const TString &runtype,const Int_t &run_num_start,const Int_t &run_num_stop,const TString &root_filename,const Int_t &nevs){
  Int_t B=1;
  if (runtype=="COSM") {cout << "By default B=0 for COSM run" << endl; B=0;}
  FillTofFromHdt(runtype, run_num_start, run_num_stop, root_filename, nevs, B);
}


void FillTofFromHdt_TS(const TString &runtype,const Int_t &run_num){
  FillTofFromHdt(runtype, run_num, run_num, "", -1, 1);
}

void TofinoStability(const TString &runtype,const Int_t &run_num_start,const Int_t &run_num_stop,const TString &root_filename){

  DeclareHistos();
  DeclareProfiles(run_num_start, run_num_stop);

//    Int_t N=run_num_stop-run_num_start+1;
//    Double_t ped_e_x[12][N], ped_e_y[12][N];
//    Double_t ped_p_x[12][N], ped_p_y[12][N];
//    Double_t bha_e_x[12][N], bha_e_y[12][N];
//    Double_t bha_p_x[12][N], bha_p_y[12][N];
//    Double_t bhabha_e_x[12][N], bhabha_e_y[12][N];
//    Double_t bhabha_p_x[12][N], bhabha_p_y[12][N];
//    Double_t kaon_e_x[12][N], kaon_e_y[12][N];
//    Double_t kaon_p_x[12][N], kaon_p_y[12][N];
//    Double_t eq_e_x[12][N], eq_e_y[12][N];
//    Double_t eq_p_x[12][N], eq_p_y[12][N];

//    Int_t N_bha_e[12], N_bha_p[12];
//    Int_t N_bhabha_e[12], N_bhabha_p[12];
//    Int_t N_kaon_e[12], N_kaon_p[12];

//    for (Int_t sl=0; sl<K_N_TOFINO_SLABS; sl++){
//      N_bha_e[sl]=0;
//      N_bha_p[sl]=0;
//      N_bhabha_e[sl]=0;
//      N_bhabha_p[sl]=0;
//      N_kaon_e[sl]=0;
//      N_kaon_p[sl]=0;
//    }

  for (Int_t run_num=run_num_start; run_num<=run_num_stop; run_num++){
    ResetHistos();
    cout << "****" << h_adce[0]->GetEntries() << endl; 
    FillTofFromHdt_TS(runtype, run_num);
    cout << "****" << h_adce[0]->GetEntries() << endl; 
    for (Int_t sl=0; sl<K_N_TOFINO_SLABS; sl++){
//        ped_e_x[sl][run_num-run_num_start]=run_num;
//        ped_e_y[sl][run_num-run_num_start]=ped_e[sl];
//        ped_p_x[sl][run_num-run_num_start]=run_num;
//        ped_p_y[sl][run_num-run_num_start]=ped_p[sl];
      p_ped_e[sl]->Fill(run_num, ped_e[sl]);
      p_ped_p[sl]->Fill(run_num, ped_p[sl]);

      cout << "Run : " << run_num << " Slab : " << sl << endl;
      cout << "**Number of Bhabha events = " << h_adce_bha[sl]->GetEntries() << endl;
      cout << "**Number of K+K- events = " << h_adce_kaon[sl]->GetEntries() << endl;
      
      if (ped_e[sl]>0&&h_adce_bha[sl]->GetEntries()>500){
//      bha_e_x[sl][N_bha_e[sl]]=run_num;
//      bha_e_y[sl][N_bha_e[sl]]=ped_e[sl]+bhapeak_e[sl];
//      eq_e_x[sl][N_bha_e[sl]]=run_num;
//      eq_e_y[sl][N_bha_e[sl]]=eq_e[sl];
//      N_bha_e[sl]++;
        p_bha_e[sl]->Fill(run_num, ped_e[sl]+bhapeak_e[sl]);
        p_eq_e[sl]->Fill(run_num, ped_e[sl]+bhapeak_e[sl]);
      }

      if (ped_p[sl]>0&&h_adcp_bha[sl]->GetEntries()>500){
//      bha_p_x[sl][N_bha_p[sl]]=run_num;
//      bha_p_y[sl][N_bha_p[sl]]=ped_p[sl]+bhapeak_p[sl];
//      eq_p_x[sl][N_bha_p[sl]]=run_num;
//      eq_p_y[sl][N_bha_p[sl]]=eq_p[sl];
//      N_bha_p[sl]++;
        p_bha_p[sl]->Fill(run_num, ped_p[sl]+bhapeak_p[sl]);
        p_eq_p[sl]->Fill(run_num, ped_p[sl]+bhapeak_p[sl]);
      }

//        if (h_adce_bha[sl]->GetEntries()>2000){
//      bhabha_e_x[sl][N_bhabha_e[sl]]=run_num;
//      bhabha_e_y[sl][N_bhabha_e[sl]]=ped_e[sl]+bhapeak_e[sl];
//      N_bhabha_e[sl]++;

//        }

//        if (h_adcp_bha[sl]->GetEntries()>2000){
//      bhabha_p_x[sl][N_bhabha_p[sl]]=run_num;
//      bhabha_p_y[sl][N_bhabha_p[sl]]=ped_p[sl]+bhapeak_p[sl];
//      N_bhabha_p[sl]++;
//        }

      if (ped_e[sl]>0&&kaonpeak_e[sl]>0){
//      kaon_e_x[sl][N_kaon_e[sl]]=run_num;
//      kaon_e_y[sl][N_kaon_e[sl]]=ped_e[sl]+kaonpeak_e[sl]/eq_e[sl];
        //      cout << "KaonPeak" << sl << " : " << N_kaon_e[sl] << " : " <<
        //        kaon_e_x[sl][N_kaon_e[sl]] << " " << kaon_e_y[sl][N_kaon_e[sl]] << endl;
//      N_kaon_e[sl]++;
        p_kaon_e[sl]->Fill(run_num, ped_e[sl]+kaonpeak_e[sl]/eq_e[sl]);
      }

      if (ped_p[sl]>0&&kaonpeak_p[sl]>0){
//      kaon_p_x[sl][N_kaon_p[sl]]=run_num;
//      kaon_p_y[sl][N_kaon_p[sl]]=ped_p[sl]+kaonpeak_p[sl]/eq_p[sl];
//      N_kaon_p[sl]++;
        p_kaon_p[sl]->Fill(run_num, ped_p[sl]+kaonpeak_p[sl]/eq_p[sl]);
      }
    }
  }

//    TGraph *g_ped_e[12];
//    TGraph *g_ped_p[12];
//    TGraph *g_bha_e[12];
//    TGraph *g_bha_p[12];
//    TGraph *g_bhabha_e[12];
//    TGraph *g_bhabha_p[12];
//    TGraph *g_kaon_e[12];
//    TGraph *g_kaon_p[12];
//    TGraph *g_eq_e[12];
//    TGraph *g_eq_p[12];

//    TString hname;
//    TString htitle;

//    TFile *fOutFile = new TFile(root_filename,"RECREATE");
//    for(Int_t sl=0;sl<K_N_TOFINO_SLABS;sl++ ){
//      g_ped_e[sl] = new TGraph(N,ped_e_x[sl],ped_e_y[sl]);
//      hname.Form("PedE_%d",sl+1);
//      htitle.Form("Tofino Pedestal slab: %d e",sl+1);
//      g_ped_e[sl]->SetName(hname);
//      g_ped_e[sl]->SetTitle(htitle);

//      g_ped_p[sl] = new TGraph(N,ped_p_x[sl],ped_p_y[sl]);
//      hname.Form("PedP_%d",sl+1);
//      htitle.Form("Tofino Pedestal slab: %d p",sl+1);
//      g_ped_p[sl]->SetName(hname);
//      g_ped_p[sl]->SetTitle(htitle);

//      g_bha_e[sl] = new TGraph(N_bha_e[sl],bha_e_x[sl],bha_e_y[sl]);
//      hname.Form("BhaE_%d",sl+1);
//      htitle.Form("Tofino Bhabha Peak slab: %d e",sl+1);
//      g_bha_e[sl]->SetName(hname);
//      g_bha_e[sl]->SetTitle(htitle);

//      g_bha_p[sl] = new TGraph(N_bha_p[sl],bha_p_x[sl],bha_p_y[sl]);
//      hname.Form("BhaP_%d",sl+1);
//      htitle.Form("Tofino Bhabha Peak slab: %d p",sl+1);
//      g_bha_p[sl]->SetName(hname);
//      g_bha_p[sl]->SetTitle(htitle);

//      g_bhabha_e[sl] = new TGraph(N_bhabha_e[sl],bhabha_e_x[sl],bhabha_e_y[sl]);
//      hname.Form("BhabhaE_%d",sl+1);
//      htitle.Form("Tofino Bhabha Peak slab: %d e (for Bhabha run)",sl+1);
//      g_bhabha_e[sl]->SetName(hname);
//      g_bhabha_e[sl]->SetTitle(htitle);

//      g_bhabha_p[sl] = new TGraph(N_bhabha_p[sl],bhabha_p_x[sl],bhabha_p_y[sl]);
//      hname.Form("BhabhaP_%d",sl+1);
//      htitle.Form("Tofino Bhabha Peak slab: %d p (for Bhabha run)",sl+1);
//      g_bhabha_p[sl]->SetName(hname); 
//      g_bhabha_p[sl]->SetTitle(htitle);

//      g_kaon_e[sl] = new TGraph(N_kaon_e[sl],kaon_e_x[sl],kaon_e_y[sl]);
//      hname.Form("KaonE_%d",sl+1);
//      htitle.Form("Tofino Kaon Peak slab: %d e",sl+1);
//      g_kaon_e[sl]->SetName(hname);
//      g_kaon_e[sl]->SetTitle(htitle);

//      g_kaon_p[sl] = new TGraph(N_kaon_p[sl],kaon_p_x[sl],kaon_p_y[sl]);
//      hname.Form("KaonP_%d",sl+1);
//      htitle.Form("Tofino Kaon Peak slab: %d p",sl+1);
//      g_kaon_p[sl]->SetName(hname);
//      g_kaon_p[sl]->SetTitle(htitle);

//      g_eq_e[sl] = new TGraph(N_bha_e[sl],eq_p_x[sl],eq_p_y[sl]);
//      hname.Form("EqE_%d",sl+1);
//      htitle.Form("Tofino Equalization Factor slab: %d e",sl+1);
//      g_eq_e[sl]->SetName(hname);
//      g_eq_e[sl]->SetTitle(htitle);

//      g_eq_p[sl] = new TGraph(N_bha_p[sl],eq_p_x[sl],eq_p_y[sl]);
//      hname.Form("EqP_%d",sl+1);
//      htitle.Form("Tofino Equalization Factor slab: %d p",sl+1);
//      g_eq_p[sl]->SetName(hname);
//      g_eq_p[sl]->SetTitle(htitle);
//      g_ped_e[sl]->Write();
//      g_ped_p[sl]->Write();
//      g_bha_e[sl]->Write();
//      g_bha_p[sl]->Write();
//      g_bhabha_e[sl]->Write();
//      g_bhabha_p[sl]->Write();
//      g_kaon_e[sl]->Write();
//      g_kaon_p[sl]->Write();
//      g_eq_e[sl]->Write();
//      g_eq_p[sl]->Write();
//    }

//    fOutFile->Close();
//    delete fOutFile;
//    fOutFile = 0;

  SaveProfiles(root_filename);
}

void TofinoStability(const TString &runtype,const Int_t &run_num_start,const Int_t &run_num_stop, const TString &calib_filename, const TString &root_filename){
  // Threshold For Kaon Selection (very tentative ver.20070126)
  // Threshold For Kaon Selection (very tentative ver.20070710)
  KaonThreshold[0]=350;
  KaonThreshold[1]=400;
  KaonThreshold[2]=400;
  KaonThreshold[3]=450;
  KaonThreshold[4]=470;
  KaonThreshold[5]=590;
  KaonThreshold[6]=500;
  KaonThreshold[7]=500;
  KaonThreshold[8]=500;
  KaonThreshold[9]=300;
  KaonThreshold[10]=360;
  KaonThreshold[11]=400;
  
  Int_t start_time, stop_time;
  TString temp1, temp2, temp3, temp4, temp5;
  TString temp6, temp7, temp8, temp9, temp10;

  std::fstream cFstrm (calib_filename.Data(), std::ios::in);
  if (cFstrm.fail()) return;

  Int_t line=-1;
  Char_t buff[1000];
  while (!cFstrm.eof()){
    cFstrm.getline(buff, 1000);
    line++;
  }
  cout << "LINE=" << line << endl;
  cFstrm.clear();
  cFstrm.seekg(0, std::ios::beg);
  cout << cFstrm.tellp() << " " << cFstrm.tellg() << endl;

  Int_t start_run[line], stop_run[line];
  Float_t PED_e[line][K_N_TOFINO_SLABS], PED_p[line][K_N_TOFINO_SLABS];
  Float_t EQ_e[line][K_N_TOFINO_SLABS], EQ_p[line][K_N_TOFINO_SLABS];
    
    for (Int_t i=0; i<line; i++){ 
      cFstrm >> start_run[i] >> start_time >> temp1 >> temp2 >> temp3 >> temp4 >> temp5 
             >>  stop_run[i] >> stop_time >> temp6 >> temp7 >> temp8 >> temp9 >> temp10;
      cout << start_run[i] << " " << stop_run[i] << endl;
      for (Int_t sl=0; sl<K_N_TOFINO_SLABS; sl++){
        cFstrm >> PED_e[i][sl];
      }
      for (Int_t sl=0; sl<K_N_TOFINO_SLABS; sl++){
        cFstrm >> PED_p[i][sl];
      }
      for (Int_t sl=0; sl<K_N_TOFINO_SLABS; sl++){
        cFstrm >> EQ_e[i][sl];
      }
      for (Int_t sl=0; sl<K_N_TOFINO_SLABS; sl++){
        cFstrm >> EQ_p[i][sl];
      }
    }

  DeclareHistos();
  DeclareProfiles(run_num_start, run_num_stop);

  for (Int_t run_num=run_num_start; run_num<=run_num_stop; run_num++){
    ResetHistos();
    cout << "****" << h_adce[0]->GetEntries() << endl; 

    // get calibration constants obtained already
    Int_t i;
    for (i=0; i<line; i++){
      cout << i << " " << run_num << " " << start_run[i] << " " << stop_run[i] << endl;
      if (run_num>=start_run[i]&&run_num<=stop_run[i]) break;
    }
    if (i==line) continue;
    cout << i  << endl;
    cout << "Use the calibration table adjusted between run " <<
       start_run[i] << " and " << stop_run[i] << endl;
    for (Int_t sl=0; sl<K_N_TOFINO_SLABS; sl++){
      ped_table_e[sl]=PED_e[i][sl];
      ped_table_p[sl]=PED_p[i][sl];
      eq_table_e[sl]=EQ_e[i][sl];
      eq_table_p[sl]=EQ_p[i][sl];
    }
    FillTofFromHdt_TS(runtype, run_num);
    cout << "****" << h_adce[0]->GetEntries() << endl; 
    for (Int_t sl=0; sl<K_N_TOFINO_SLABS; sl++){
      p_ped_e[sl]->Fill(run_num, ped_e[sl]);
      p_ped_p[sl]->Fill(run_num, ped_p[sl]);

      cout << "Run : " << run_num << " Slab : " << sl << endl;
      cout << "**Number of Bhabha events = " << h_adce_bha[sl]->GetEntries() << endl;
      cout << "**Number of K+K- events = " << h_adce_kaon[sl]->GetEntries() << endl;
      
      if (ped_e[sl]>0&&h_adce_bha[sl]->GetEntries()>500){
        p_bha_e[sl]->Fill(run_num, ped_e[sl]+bhapeak_e[sl]);
        p_EQbha_e[sl]->Fill(run_num, (ped_e[sl]+bhapeak_e[sl]-PED_e[i][sl])*EQ_e[i][sl]);
        p_eq_e[sl]->Fill(run_num, eq_e[sl]);
      }

      if (ped_p[sl]>0&&h_adcp_bha[sl]->GetEntries()>500){
        p_bha_p[sl]->Fill(run_num, ped_p[sl]+bhapeak_p[sl]);
        p_EQbha_p[sl]->Fill(run_num, (ped_p[sl]+bhapeak_p[sl]-PED_p[i][sl])*EQ_p[i][sl]);
        p_eq_p[sl]->Fill(run_num, eq_p[sl]);
      }

      if (ped_e[sl]>0&&kaonpeak_e[sl]>0){
        p_kaon_e[sl]->Fill(run_num, ped_e[sl]+kaonpeak_e[sl]/eq_e[sl]);
        p_EQkaon_e[sl]->Fill(run_num, (ped_e[sl]+kaonpeak_e[sl]/eq_e[sl]-PED_e[i][sl])*EQ_e[i][sl]);
      }

      if (ped_p[sl]>0&&kaonpeak_p[sl]>0){
        p_kaon_p[sl]->Fill(run_num, ped_p[sl]+kaonpeak_p[sl]/eq_p[sl]);
        p_EQkaon_p[sl]->Fill(run_num, (ped_p[sl]+kaonpeak_p[sl]/eq_p[sl]-PED_p[i][sl])*EQ_p[i][sl]);
      }
      p_EQ_e[sl]->Fill(run_num, EQ_e[i][sl]);
      p_EQ_p[sl]->Fill(run_num, EQ_p[i][sl]);

    }
  }
   SaveProfiles(root_filename);
}

void MakePeakTable(const TString &runtype,const Int_t &run_num_start,const Int_t &run_num_stop, const TString &directory){

  DeclareHistos();
  DeclareProfiles(run_num_start, run_num_stop);

  for (Int_t run_num=run_num_start; run_num<=run_num_stop; run_num++){
    ResetHistos();
    cout << "****" << h_adce[0]->GetEntries() << endl; 

    FillTofFromHdt_TS(runtype, run_num);
    cout << "****" << h_adce[0]->GetEntries() << endl; 
    
    TString table_filename="./";
    table_filename += directory;
    table_filename += runtype;
    table_filename += "_";
    table_filename += run_num;
    table_filename += ".txt";

    std::fstream cFstrm (table_filename.Data(), std::ios::out);
    if (cFstrm.fail()) return;

    for (Int_t sl=0; sl<K_N_TOFINO_SLABS; sl++){
      cFstrm << ped_e[sl] << " " << ped_p[sl] << " "
             << h_adce_bha[sl]->GetEntries() << " " << ped_e[sl]+bhapeak_e[sl] << " " 
             << h_adcp_bha[sl]->GetEntries() << " " << ped_p[sl]+bhapeak_p[sl] << " "
             << h_adce_kaon[sl]->GetEntries() << " " << ped_e[sl]+kaonpeak_e[sl]/eq_e[sl] << " "
             << h_adcp_kaon[sl]->GetEntries() << " " << ped_p[sl]+kaonpeak_p[sl]/eq_p[sl] << " "
             << eq_e[sl] << " " << eq_p[sl] << std::endl;
    }  }
}

---