///////////////////////////////////////////////////////////////////////////////
//
///////////////////////////////////////////////////////////////////////////////
#include "TH1.h"
#include "TFile.h"
#include "TTree.h"
#include "TBranch.h"
#include "TFolder.h"
#include "TSystem.h"

#include "OfflineMon/TPhysMon.hh"

#include "Stntuple/oracle/TCdfOracle.hh"
#include "Stntuple/obj/TStnRunSummary.hh"

ClassImp(TPhysMon)


TF1*    TPhysMon::fgFuncEp;
//_____________________________________________________________________________
double f_ep(double* X, double* Par) {
  // straw-man parametrization of E/P distribution - function taken from
  // the real histogram

  double x[] = {
    0.025,0.075, 0.125, 0.175,0.225, 0.275, 0.325,  0.375,  0.425,  0.475,
    0.525,0.575, 0.625, 0.675,0.725, 0.775, 0.825,  0.875,  0.925,  0.975,
    1.025,1.075, 1.125, 1.175,1.225, 1.275, 1.325,  1.375,  1.425,  1.475,
    1.525,1.575, 1.625, 1.675,1.725, 1.775, 1.825,  1.875,  1.925,  1.975,
    2.025,2.075, 2.125, 2.175
  };

  double f[] = {  
    0.000,   0.000,    0.000,    0.000,    0.000,
    0.000,   0.000,    0.000,    0.000,    0.000,
    3.000,   4.000,    6.000,    9.000,   15.000,
    40.000, 102.000,  305.000,  862.000, 1391.000,
    1212.000, 805.000,  494.000,  364.000,  283.000,
    211.000, 188.000,  170.000,  136.000,  134.000,
    114.000,  85.000,   82.000,   78.000,   73.000,
    63.000,  62.000,   62.000,   48.000,   44.000,
    10.000,   0.000,    0.000,    0.000
  };

  // 2 parameters: Par[0] : overall scale, Par[1]: offset in X
  // Assumes the same binning!
 
  double ftouse = 0;
  double bin    = 0.05;
  double xx     = X[0]+Par[1];
  
  for (int i=0; i<43; i++) {
    if ( (x[i] <= xx) && (x[i+1] >= xx)) {
      ftouse = f[i]+(xx-x[i])/bin*(f[i+1]-f[i]);
      break;
    }
  }

  double dxx = 0;
  //  if (xx > 1) dxx = 5*(xx-1);
  if (xx > Par[3]) dxx = 7*(xx-Par[3]);

  return Par[0]*(ftouse + Par[2]*(1-TMath::Exp(-dxx)));

}

//_____________________________________________________________________________
TPhysMon::TPhysMon() {
  fgFuncEp = new TF1("func_ep",f_ep,0.85,1.9,4);
  fgFuncEp->SetParNames("scale","offset","constant","thr");
  Clear();
}

//_____________________________________________________________________________
TPhysMon::~TPhysMon() {
  delete fgFuncEp;
}

//_____________________________________________________________________________
void TPhysMon::Clear(const char* Option) {


  fData->fEOverP        = -1.;
  fData->fSigEOverP     = -1;

  fData->fEoverPincMean = -1.;
  fData->fEoverPincN    = -1.;
  fData->fEoverPincErr  = -1.;

  fData->fEoverPNEMean = -1.;
  fData->fEoverPNEN    = -1.;
  fData->fEoverPNEErr  = -1.;

  fData->fEoverPNWMean = -1.;
  fData->fEoverPNWN    = -1.;
  fData->fEoverPNWErr  = -1.;
}

//_____________________________________________________________________________
int truncated_mean(TH1F* Hist, float X0, float Road, int NIter, 
		   float& XMean, float& SigX) {

  double x, nev, sumx, sumx2, sumn, sigx, xmean, xm, road; 

  int    nbins   = Hist->GetNbinsX();

  road  = Road;
  xmean = X0;

  for (int it=0; it<NIter; it++) {

    sumn  = 0;
    sumx  = 0;
    sumx2 = 0;

    for (int i=1; i<=nbins; i++) {
      x   = Hist->GetBinCenter(i);
      if (TMath::Abs(x-xmean) < road) {
	nev = Hist->GetBinContent(i);
	sumn  += nev;
	sumx  += x*nev;
	sumx2 += x*x*nev;
      }
    }

    if (sumn > 0) {
      xm    = sumx/sumn;
      sigx  = TMath::Sqrt((sumx2/sumn-xm*xm)/sumn);
    }
    else {
      xm    = 0;
      sigx  = 0;
    }
//-----------------------------------------------------------------------------
//end of iteration
//-----------------------------------------------------------------------------
//    road = Road+TMath::Abs(xmean-xm);
    xmean = xm;

    printf(" it, xmean, sigx = %3i %10.4f %10.4f\n",it,xmean,sigx);
  }

  XMean = xmean;
  SigX  = sigx;

  return 0;
}

double f_HadE(double* X, double* Par) {
  // straw-man parametrization of CHA Energy 
  // distribution for muons - function taken from the real histogram
  double x[] = {
    0.5,        0.6, 	     0.7, 	  0.8, 	       0.9, 
    1.0,        1.1, 	     1.2, 	  1.3, 	       1.4, 
    1.5,        1.6, 	     1.7, 	  1.8, 	       1.9, 
    2.0,        2.1, 	     2.2, 	  2.3, 	       2.4, 
    2.5,        2.6, 	     2.7, 	  2.8, 	       2.9, 
    3.0,        3.1, 	     3.2, 	  3.3, 	       3.4  
  };

  double f[] = {  
    740.000,    760.000,     800.000,     860.000,     1220.000,
    1600.000,   2100.000,    2500.000,    3200.000,    3600.000,
    3900.000,   4100.000,    4000.000,    3700.000,    3200.000,
    2700.000,   2200.000,    1800.000,    1500.000,    1250.000,
    880.000,    690.000,     600.000,     550.000,     400.000,
    300.000,    280.000,     260.000,     240.000,     200.000
  };

  // 2 parameters: Par[0] : overall scale, Par[1]: offset in X
  // Assumes the same binning!
 
  double ftouse = 0;
  double bin    = 0.1;
  double xx     = X[0]+Par[1];

  double const dWeight = 4.50608509520177904e+03;
  
  // Calculate f( xx) called ftouse
  for (int i=0; i<30; i++) {
    if ( (x[i] <= xx) && (x[i+1] >= xx)) {
      ftouse = f[i]+(xx-x[i])/bin*(f[i+1]-f[i]);
      break;
    }
  }

  return Par[0] * ftouse / dWeight;
}

//_____________________________________________________________________________
double f_EmE(double* X, double* Par) {
  // straw-man parametrization of CEM Energy 
  // distribution for muons - function taken from the real histogram
  double x[] = {
    0.1,        0.12, 	     0.14, 	  0.16,        0.18, 
    0.2,        0.22, 	     0.24, 	  0.26,        0.28, 
    0.3,        0.32, 	     0.34, 	  0.36,        0.38, 
    0.4,        0.42, 	     0.44, 	  0.46,        0.48, 
    0.5,        0.52, 	     0.54, 	  0.56,        0.58,         0.6
  };

  double f[] = {  
    500.000,    680.000,     950.000,     1200.000,    1780.000,
    2730.000,   4000.000,    5140.000,    5960.000,    5880.000,
    5590.000,   4590.000,    4240.000,    3500.000,    2740.000,
    2000.000,   1750.000,    1410.000,    1270.000,    1190.000,
    1000.000,   1000.000,    930.000,     880.000,     770.000,      620.000
  };

  // 2 parameters: Par[0] : overall scale, Par[1]: offset in X
  // Assumes the same binning!
 
  double ftouse = 0;
  double bin    = 0.02;
  double xx     = X[0]+Par[1];

  double const dWeight = 1.22231088516361524e+03;
  
  // Calculate f( xx) called ftouse
  for (int i = 0; i < 26; i++) {
    if ( ( x[i] <= xx) && ( x[i + 1] >= xx)) {
      ftouse = f[i]  +( xx - x[i]) / bin * ( f[i + 1] - f[i]);
      break;
    }
  }

  return Par[0] * ftouse / dWeight;
}

//_____________________________________________________________________________
int fit_njpsi(TH1F* Hist, float& NJPsi, float& M, float& MFitErr, float& SigM){
  
  TF1 *func = new TF1("func","gaus(0)+pol0(3)",2.9,3.3);

  // Set initial values and parameter names based on the histogram contents

  TAxis* ax = Hist->GetXaxis();
  double xmin, xmax, bin, bin1,bin2;

  xmin  = ax->GetXmin();
  xmax  = ax->GetXmax();
  bin   = Hist->GetBinWidth(1);

  bin1  = (3.04-xmin)/bin + 1;
  bin2  = (3.15-xmin)/bin + 1;
//-----------------------------------------------------------------------------
//  assume sigma = 15 MeV ..
//-----------------------------------------------------------------------------
  double m0   = 3.090;
  double sig0 = 0.015;
  double a0   = Hist->Integral(bin1,bin2)*bin/(sqrt(2.*TMath::Pi())*0.015);

  func->SetParameters(a0,m0,sig0,0);
  func->SetParNames("Constant","Mean_value","Sigma","Constant");
  func->SetParLimits(0,0,1000000);
  func->SetParLimits(1,3.05,3.15);
  func->SetParLimits(2,0.011,0.03);
  func->SetParLimits(3,0,1000000);
                                          // Fit histogram
  Hist->Fit("func","r");
    
  double lowlim = func->GetParameter(1)-3*func->GetParameter(2);
  double higlim = func->GetParameter(1)+3*func->GetParameter(2);

                                          // integral of function between fit limits
  double binwid = Hist->GetBinWidth(1);
  double totint = (func->Integral(lowlim,higlim))/binwid;
   
                                          // flat function to integrate background

  TF1 *f_flat =new TF1("f_flat","pol0",2.9,3.3);
  double mypar = func->GetParameter(3);
  f_flat->SetParameters(&mypar);

  double backint = (f_flat->Integral(lowlim,higlim))/binwid;
      
  double sigint = totint - backint;

  printf("Total Integral (%g to %g): %g\n",lowlim,higlim,totint);
  printf("Bkg Integral:                        %g\n",backint);
  printf("Signal Integral (# Jpsis):s          %g\n",sigint);
  
  NJPsi   = sigint;
  M       = func->GetParameter(1);
  MFitErr = func->GetParError (1);
  SigM    = func->GetParameter(2);
    
  func->Delete();
  f_flat->Delete();

  return 0;
}

//_____________________________________________________________________________
int TPhysMon::UpdateNtuple() {

  // for physmon the 1st parameter is the name of the input directory

  TFile* ntuple_file = new TFile(fNtupleFileName.Data(),"recreate");
  TTree* tree        = new TTree("physmon","Physics Monitor");
  TBranch* br        = tree->Branch("PhysMon","TPhysMon", &fData);
//-----------------------------------------------------------------------------
//  read the data files and fill the ntuple
//-----------------------------------------------------------------------------
  char cmd[200];

  TStnRunSummary rs;

  sprintf(cmd,"rm .tmp; ls %s/stra_mon_hist.* > .tmp",fInputDirName.Data());
  gSystem->Exec(cmd);
  FILE* f = fopen(".tmp","r");

  char fn[200];

  TFile*   hfile;
  TH1F*    hist;
  TF1      *ff, *fep;
  TF1*     fHadE = new TF1( "fit1", f_HadE, 0.6,  2.5, 2);;
  TF1*     hEmE  = new TF1( "fit2", f_EmE,  0.12, 0.6, 2);
  TFolder* fol;

  //  double   xmean, sigx;

  fep = TPhysMon::GetFuncEp();
  fep->SetParameters(1,0.0,100.);

  double fint;

  fint = 3583.9380; // fint = fep->Integral(0.85,1.9);

  while (fgets(fn,200,f) > 0) {
    int len = strlen(fn);
    fn[len-1] = 0;
    int run_number = atoi(&fn[len-7]);
    printf(".%s.  %i\n",fn,run_number);
    fData->fRunNumber = run_number;
    Clear();
//-----------------------------------------------------------------------------
//  open ROOT histogram file and get the histograms we need
//-----------------------------------------------------------------------------
    hfile = TFile::Open(fn);
    fol   = (TFolder*) hfile->Get("Ana");
//-----------------------------------------------------------------------------
//  WenuMon histograms
//  1. Missing Et
//-----------------------------------------------------------------------------
    hist  = (TH1F*) fol->FindObject("WenuMon/Hist/met0_0");
    fData->fNWenu = hist->GetEntries();
//-----------------------------------------------------------------------------
//  2. E/P for isolated high-Pt electrons (but not necessarily W-electrons,
//     no MET cut yet), no normalization to chi2 yet ...
//-----------------------------------------------------------------------------
    double n0;
    hist  = (TH1F*) fol->FindObject("WenuMon/Hist/ele_1/ep");
    if (hist) {
      n0 = hist->GetEntries()/fint;
      fep->SetParameters(n0,0.01,n0/10.,0.8);
      hist->Fit(fep->GetName(),"r");

      fData->fEOverP    = 1 + fep->GetParameter(1);
      fData->fSigEOverP = fep->GetParError(1);

      //truncated_mean(hist,1.0,0.2,3,fData->fEOverP,fData->fSigEOverP);
    }

    hist  = (TH1F*) fol->FindObject("WenuMon/Hist/ele_1/ep_pos");
    if (hist) {
      n0 = hist->GetEntries()/fint;
      fep->SetParameters(n0,0.01,n0/10.,0.8);
      hist->Fit(fep->GetName(),"r");

      fData->fEPPos    = 1 + fep->GetParameter(1);
      fData->fSigEPPos = fep->GetParError(1);

      // truncated_mean(hist,1.0,0.2,3,fData->fEPPos,fData->fSigEPPos);
    }

    hist  = (TH1F*) fol->FindObject("WenuMon/Hist/ele_1/ep_neg");
    if (hist) {
      n0 = hist->GetEntries()/fint;
      fep->SetParameters(n0,0.01,n0/10.,0.8);
      hist->Fit(fep->GetName(),"r");

      fData->fEPNeg    = 1 + fep->GetParameter(1);
      fData->fSigEPNeg = fep->GetParError(1);
      // truncated_mean(hist,1.0,0.2,3,fData->fEPNeg,fData->fSigEPNeg);
    }


    hist             = (TH1F*) fol->FindObject("WenuMon/Hist/ee_mass_0");
    fData->fNZee  = hist->GetEntries();
    hist->Fit("gaus");
    ff = hist->GetFunction("gaus");
    fData->fMZee         = ff->GetParameter(1);
    fData->fMZeeFitErr   = ff->GetParError(1);
    fData->fSigMZee      = ff->GetParameter(2);

    hist = (TH1F*) fol->FindObject("WenuMon/Hist/EoverP_inc");
    if (hist) {
      fData->fEoverPincMean = hist->GetMean();
      fData->fEoverPincN = hist->Integral();
      fData->fEoverPincErr = (hist->GetRMS())/sqrt(hist->Integral());
    }

    hist = (TH1F*) fol->FindObject("WenuMon/Hist/EoverP_NE");
    if (hist) {
      fData->fEoverPNEMean = hist->GetMean();
      fData->fEoverPNEN = hist->Integral();
      fData->fEoverPNEErr = (hist->GetRMS())/sqrt(hist->Integral());
    }

    hist = (TH1F*) fol->FindObject("WenuMon/Hist/EoverP_NW");
    if (hist) {
      fData->fEoverPNWMean = hist->GetMean();
      fData->fEoverPNWN = hist->Integral();
      fData->fEoverPNWErr = (hist->GetRMS())/sqrt(hist->Integral());
    }

    hist = (TH1F*) fol->FindObject("WenuMon/Hist/EoverP_SE");
    if (hist) {
      fData->fEoverPSEMean = hist->GetMean();
      fData->fEoverPSEN = hist->Integral();
      fData->fEoverPSEErr = (hist->GetRMS())/sqrt(hist->Integral());
    }
    else {
      fData->fEoverPSEMean = -1.;
      fData->fEoverPSEN    = -1.;
      fData->fEoverPSEErr  = -1.;
    }

    hist = (TH1F*) fol->FindObject("WenuMon/Hist/EoverP_SW");
    if (hist) {
      fData->fEoverPSWMean = hist->GetMean();
      fData->fEoverPSWN = hist->Integral();
      fData->fEoverPSWErr = (hist->GetRMS())/sqrt(hist->Integral());
    }
    else {
      fData->fEoverPSWMean = -1.;
      fData->fEoverPSWN    = -1.;
      fData->fEoverPSWErr  = -1.;
    }

    hist = (TH1F*) fol->FindObject("WenuMon/Hist/ZMass_inc");
    if (hist) {
      fData->fZMassincMean = hist->GetMean();
      fData->fZMassincN = hist->Integral();
      fData->fZMassincErr = (hist->GetRMS())/sqrt(hist->Integral());
    }
    else {
      fData->fZMassincMean = -1.;
      fData->fZMassincN    = -1.;
      fData->fZMassincErr  = -1.;
    }

    hist = (TH1F*) fol->FindObject("WenuMon/Hist/ZMass_CC");
    if (hist) {
      fData->fZMassCCMean = hist->GetMean();
      fData->fZMassCCN = hist->Integral();
      fData->fZMassCCErr = (hist->GetRMS())/sqrt(hist->Integral());
    }
    else {
      fData->fZMassCCMean = -1.;
      fData->fZMassCCN    = -1.;
      fData->fZMassCCErr  = -1.;
    }

    hist = (TH1F*) fol->FindObject("WenuMon/Hist/ZMass_CP");
    if (hist) {
      fData->fZMassCPMean = hist->GetMean();
      fData->fZMassCPN = hist->Integral();
      fData->fZMassCPErr = (hist->GetRMS())/sqrt(hist->Integral());
    }
    else {
      fData->fZMassCPMean = -1.;
      fData->fZMassCPN    = -1.;
      fData->fZMassCPErr  = -1.;
    }

    hist = (TH1F*) fol->FindObject("WenuMon/Hist/ZMass_CPW");
    if (hist) {
      fData->fZMassCPWMean = hist->GetMean();
      fData->fZMassCPWN = hist->Integral();
      fData->fZMassCPWErr = (hist->GetRMS())/sqrt(hist->Integral());
    }
    else {
      fData->fZMassCPWMean = -1.;
      fData->fZMassCPWN    = -1.;
      fData->fZMassCPWErr  = -1.;
    }

    hist = (TH1F*) fol->FindObject("WenuMon/Hist/ZMass_CPE");
    if (hist) {
      fData->fZMassCPEMean = hist->GetMean();
      fData->fZMassCPEN = hist->Integral();
      fData->fZMassCPEErr = (hist->GetRMS())/sqrt(hist->Integral());
    }
    else {
      fData->fZMassCPEMean = -1.;
      fData->fZMassCPEN    = -1.;
      fData->fZMassCPEErr  = -1.;
    }
				// number of COT hits

    hist = (TH1F*) fol->FindObject("WenuMon/Hist/ele_1/trk_ncot_hits_ax");
    fData->fNCotHitsAx    = hist->GetMean();
    fData->fNCotHitsAxErr = hist->GetRMS()/sqrt(hist->GetEntries()+0.0001);

    hist = (TH1F*) fol->FindObject("WenuMon/Hist/ele_1/trk_ncot_hits_st");
    fData->fNCotHitsSt    = hist->GetMean();
    fData->fNCotHitsAxErr = hist->GetRMS()/sqrt(hist->GetEntries()+0.0001);
//-----------------------------------------------------------------------------
//  JpsiMon histograms
//-----------------------------------------------------------------------------
    hist = (TH1F*) fol->FindObject("JpsiMon/Hist/mumu_mass_0");

    fData->fNJPsiEnts = hist->GetEntries();;
    fit_njpsi(hist,
	      fData->fNJPsimumu,
	      fData->fMJPsi,
	      fData->fMJPsiFitErr,
	      fData->fSigMJPsi);

    // J/psi MIP peak
    hist = (TH1F*) fol->FindObject("JpsiMon/Hist/CMUemE");
    if (hist) {
      fData->fJpsiCMUEmMean = hist->GetMean();
      fData->fJpsiCMUEmN = hist->Integral();
      fData->fJpsiCMUEmErr = (hist->GetRMS())/sqrt(hist->Integral());
    }
    else {
      fData->fJpsiCMUEmMean = -1.;
      fData->fJpsiCMUEmN    = -1.;
      fData->fJpsiCMUEmErr  = -1.;
    }
    hist = (TH1F*) fol->FindObject("JpsiMon/Hist/CMUhadE");
    if (hist) {
      fData->fJpsiCMUHadMean = hist->GetMean();
      fData->fJpsiCMUHadN = hist->Integral();
      fData->fJpsiCMUHadErr = (hist->GetRMS())/sqrt(hist->Integral());
    }
    else {
      fData->fJpsiCMUHadMean = -1.;
      fData->fJpsiCMUHadN    = -1.;
      fData->fJpsiCMUHadErr  = -1.;
    }
    hist = (TH1F*) fol->FindObject("JpsiMon/Hist/emE");
    if (hist) {
      n0 = hist->Integral( 6, 31);
      hEmE->SetParameters( n0, 1.e-3);
      hist->Fit( hEmE->GetName(), "r");

      fData->fJpsiEmE        = 0.26 + hEmE->GetParameter( 1);
      fData->fJpsiEmEErr     = hEmE->GetParError( 1);
      printf( "  -+<==>+-  fJpsiEmE: %10.4f + %10.4f\n", fData->fJpsiEmE, fData->fJpsiEmEErr);
    }
    else {
      fData->fJpsiEmE        = -1.;
      fData->fJpsiEmEErr     = -1.;
    }
    hist = (TH1F*) fol->FindObject("JpsiMon/Hist/hadE");
    if (hist) {
      n0 = hist->Integral( 12, 50);
      fHadE->SetParameters( n0, 1.e-3);
      hist->Fit( fHadE->GetName(), "r");

      fData->fJpsiHadE       = 1.6 + fHadE->GetParameter( 1);
      fData->fJpsiHadEErr    = fHadE->GetParError( 1);
      printf( "  -+<==>+- fJpsiHadE: %10.4f + %10.4f\n", fData->fJpsiHadE, fData->fJpsiHadEErr);
    }
    else {
      fData->fJpsiHadE       = -1.;
      fData->fJpsiHadEErr    = -1.;
    }
//-----------------------------------------------------------------------------
//  WmunuMon histograms
//-----------------------------------------------------------------------------
    hist = (TH1F*) fol->FindObject("WmunuMon/Hist/EmEnergy_CMUP");
    if (hist) {
      fData->fEmEnergyCMUPMean = hist->GetMean();
      fData->fEmEnergyCMUPN = hist->Integral();
      fData->fEmEnergyCMUPErr = (hist->GetRMS())/sqrt(hist->Integral());
    }
    else {
      fData->fEmEnergyCMUPMean = -1.;
      fData->fEmEnergyCMUPN    = -1.;
      fData->fEmEnergyCMUPErr  = -1.;
    }

    hist = (TH1F*) fol->FindObject("WmunuMon/Hist/EmEnergy_CMX");
    if (hist) {
      fData->fEmEnergyCMXMean = hist->GetMean();
      fData->fEmEnergyCMXN = hist->Integral();
      fData->fEmEnergyCMXErr = (hist->GetRMS())/sqrt(hist->Integral());
    }
    else {
      fData->fEmEnergyCMXMean = -1.;
      fData->fEmEnergyCMXN    = -1.;
      fData->fEmEnergyCMXErr  = -1.;
    }

    hist = (TH1F*) fol->FindObject("WmunuMon/Hist/HadEnergy_CMUP");
    if (hist) {
      fData->fHadEnergyCMUPMean = hist->GetMean();
      fData->fHadEnergyCMUPN = hist->Integral();
      fData->fHadEnergyCMUPErr = (hist->GetRMS())/sqrt(hist->Integral());
    }
    else {
      fData->fHadEnergyCMUPMean = -1.;
      fData->fHadEnergyCMUPN    = -1.;
      fData->fHadEnergyCMUPErr  = -1.;
    }

    hist = (TH1F*) fol->FindObject("WmunuMon/Hist/HadEnergy_CMX");
    if (hist) {
      fData->fHadEnergyCMXMean = hist->GetMean();
      fData->fHadEnergyCMXN = hist->Integral();
      fData->fHadEnergyCMXErr = (hist->GetRMS())/sqrt(hist->Integral());
    }
    else {
      fData->fHadEnergyCMXMean = -1.;
      fData->fHadEnergyCMXN    = -1.;
      fData->fHadEnergyCMXErr  = -1.;
    }
//-----------------------------------------------------------------------------
//  luminosity histogram: for 1 run/file GetTotal returns the luminosity for 
//  this run
//-----------------------------------------------------------------------------
//     TH1F* lum_hist = (TH1F*) fol->FindObject("int_lumi_tape");
//     fData->fLumiTape = lum_hist->Integral();

    fCdfofprd->GetRunSummary(run_number,&rs);
    fData->fLumiTape = rs.LumiTape();

    fData->fCalendarTime = 0;
    fData->fNZmumu       = 0;
    //    fData->fNJPsimumu    = 0;
    //    fData->fMJPsi        = 0;
    //    fData->fSigMJPsi     = 0;
    tree->Fill();

    hfile->Close();
    delete hfile;
    //    delete fol;
  }
  ntuple_file->Write();
  delete ntuple_file;
  
  return 0;
}