//////////////////////////////////////////////////////////////////////////
//
// PEN Maximum Likelihood RooFit Code
// 
// https://twiki.cern.ch/twiki/bin/view/RooStats/WebHome
//
// Define basic PEN p.d.f's
//
// pdf list:
//
// 04/2009 - Emil Frlez 
// 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "TFile.h"
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooNovosibirsk.h"
#include "RooExponential.h" 
#include "TCanvas.h"
#include "RooPlot.h"
#include <iostream.h>
#include "RooClassFactory.h"
#include "TROOT.h"

#ifndef __CINT__
#include "MichelDecayPdf.h"
#include "MichelDecayPdf.cxx"
#endif

using namespace RooFit ;

// fit model, binning, number of samples, number of events, iteration number (0 for gen_writeout, pdf generation, -1 for histo fitting with different parms)

int simplePEN(int n=4, int tbin=280, int ebin=200, int zbin=150, int nsample=10, long int nev=50000, int iseed=6159, int niter=1)
//int simplePEN(int n=7, int tbin=0, int ebin=0, , int zbin=150,int nsample=10, int nev=1E6, int iseed=223, int niter=1)
//int simplePEN(int n=100, int tbin=2800, int ebin=2000, int zbin=150, int nsample=10, int nev=1E4, int iseed=223, int niter=1)
//int simplePEN(int n=101, int tbin=280, int ebin=200, int zbin=150, int nsample=10, int nev=1E5, int iseed=223, int niter=1)
{

  gROOT->Reset();
  gROOT->SetStyle("Plain");

  if (n < 100)TCanvas* c = new TCanvas("MCstudy","MCstudy",1100,900) ;
  if (n == 100)TCanvas* c = new TCanvas("MCstudy","MCstudy",900,1100) ;
  if (n == 101)TCanvas* c = new TCanvas("MCstudy","MCstudy",1100,900) ;

  TLatex lat ;
  lat.SetTextFont(42); // helvetica-medium-r-normal
  lat.SetTextSize(0.025);
  lat.SetTextAlign(12);

  // TFile f("run84101.root") ;
  // t1=(TTree *) gROOT->FindObject("usertree") ;

  // output text file: mean fitted events +- rms
  if (tbin > 0) ofstream out(Form("datasimple/nsigb%d.dat",niter)) ;
  if (tbin == 0) ofstream out(Form("datasimple/nsigu%d.dat",niter)) ;
  //ofstream out(Form("dataout/nsig.dat")) ;
  ofstream out280("ht280.dat") ;
  ofstream out2800("ht2800.dat") ;

  //  Output root file
  if (tbin > 0) TFile *file  = new TFile(Form("rootsimple/t%db_%.0eev.root",tbin,(float) nev),"RECREATE");
  if (tbin == 0) TFile *file  = new TFile(Form("rootsimple/tunb_%.0eev.root",(float) nev),"RECREATE");

   //t1->Draw("clmp.n_clump") ;
   //c->Update();

  // S e t u p   t i m e   m o d e l s
  // ---------------------------------

  // 
  Int_t binned=0 ;

  RooRealVar t("t","t",-80,200,"ns") ; // time scale: -50+200 ns
  RooRealVar dt("dt","per-event-error on t",1,0,5) ; // this has to equal decay time rms
  dt.setConstant(kTRUE) ;

  RooRealVar toff_mich("toff_mich","toff_mich",0.1,-10,10) ; // RooFactory t_mich time offset: -3+3 ns
  toff_mich.setConstant(kTRUE) ;                            

  RooRealVar bias_pt("bias_pt","bias",0.1,-10,10) ; // decay time offset for p2e
  bias_pt.setConstant(kTRUE) ; 
  RooRealVar bias_mt("bias_mt","bias",0.1,-10,10) ; // decay time offset for Michel
  bias_mt.setConstant(kTRUE) ; 

  RooRealVar s_pt("s_pt","per-event error scale factor",1,0.1,10) ; // TDC rms for p2e
  s_pt.setConstant(kTRUE) ;
  RooGaussModel g_pt("g_pt","gauss model scaled by per-event error",t,bias_pt,s_pt,dt) ;
  RooRealVar s_mt("s_mt","per-event error scale factor",1,0.1,10) ;  // TDC rms for Michel
  s_mt.setConstant(kTRUE) ;
  RooGaussModel g_mt("g_mt","gauss model scaled by per-event error",t,bias_mt,s_mt,dt) ;

  RooRealVar tau_pi("tau_pi","pion lifetime",26.033,24.,28.) ; // parameter variable
  tau_pi.setConstant(kTRUE) ; 
  RooRealVar tau_mu("tau_mu","muon lifetime",2197,2190,2300) ; // parameter variable
  tau_mu.setConstant(kTRUE) ; 

  // Build p2e smeared timing p.d.f
  RooDecay t_p2e("t_p2e","pi+ decay",t,tau_pi,g_pt,RooDecay::SingleSided) ; // use RooDecay or 
  //#ifdef __CINT__                                                         // custom alternative below! 
  //   gROOT->ProcessLineSync(".x ts_p2e.cxx+") ;
  //#endif
  // load the t_p2e Pdf class
  //gROOT->ProcessLineSync(".L ts_p2e.cxx+") ; // compile and load
  //ts_p2e t_p2e("t_p2e","pi+ decay",t,tau_pi,toff_mich,g_pt,ts_p2e::SingleSided) ;

  // Build Michel smeared timing p.d.f
  RooDecay t_mu("t_mu","mu+ decay",t,tau_mu,g_mt,RooDecay::SingleSided) ;

  // Detector time resolution (response) function for Michel decay
  RooRealVar tmg("tmg","tmg",0) ; // mean
  tmg.setConstant(kTRUE) ; 
  RooRealVar tsg("tsg","tsg",1,0.1,3) ; // rms in ns
  tsg.setConstant(kTRUE) ; 
  RooGaussian tgauss("tgauss","tgauss",t,tmg,tsg) ;
  // Define sampling frequency
  t.setBins(10000,"fft") ; // use binning only with FFT

  // Build pi-mu-e smeared timing p.d.f using RooClassFactory (make template and edit)
  //RooClassFactory::makePdf("MichelDecayPdf","t,tau_mu,tau_pi,toff_mich") ;
  //RooClassFactory::makePdf("PrompttPdf","t,trms_pr,toff_mich") ;
  //#ifdef __CINT__
  //   gROOT->ProcessLineSync(".x MichelDecayPdf.cxx+") ;
  //#endif

  // load the MichelDecayPdf class
  gROOT->ProcessLineSync(".L MichelDecayPdf.cxx+") ; // compile and load
  // Instance of the Michel Decay Class
  MichelDecayPdf t_mic("t_mic","pi+->mu+->e+ chain decay",t,tau_mu,tau_pi,toff_mich) ;
  // smeared Michel Decay Class - FFT
  RooFFTConvPdf t_mich("t_mich","Michel decay time t_mich (x) gauss",t,t_mic,tgauss) ;
  // smeared Michel Decay Class - numerical convolution
  //RooNumConvPdf t_mich("t_mich","Michel decay time t_mich (x) gauss",t,t_mic,tgauss) ;
  //t_mich.setConvolutionWindow(0,5,1) ;

  // accidental events, flat time distribution
  RooPolynomial t_acc("t_acc","accidental time PDF",t,RooArgList()) ; // flat background

  // prompt events
  RooRealVar trms_pr("trms_pr","trms_pr",1,0.5,2) ; // mean
  trms_pr.setConstant(kTRUE) ; 
  gROOT->ProcessLineSync(".L PrompttPdf.cxx+") ; // compile and load
  PrompttPdf t_prom("t_prom","Prompt time",t,trms_pr,toff_mich) ;

  // decay time model: p2e+Michel+accidental+DIF
  RooRealVar nsig("nsig","signal events",nev,1,1E9) ;             // p2e signal fraction
  RooRealVar nbkg1("nbkg1","Michel background events",nev,1,1E9) ; // Michel fraction
  RooRealVar nbkg2("nbkg2","acc background events",nev,1,1E9) ;     // accidental fraction
  RooRealVar nbkg3("nbkg3","DIF background events",nev/100,1,1E9) ;     // DIF fraction
  RooRealVar nbkg4("nbkg4","Prompt background events",nev/50,1,1E9) ;  // Prompt fraction
  //nbkg2.setConstant(kTRUE) ;  

  Double_t true_nsig=(double) nev ;              // p2e signal
  Double_t true_nbkg1=1.50*true_nsig ;  // Michel
  Double_t true_nbkg2=1.00*true_nsig ;  // Accidentals
  Double_t true_nbkg3=0.20*true_nsig ;   // DIF
  Double_t true_nbkg4=0.05*true_nsig ;  // prompt

  long int nevtot ;

  nsig.setVal(true_nsig) ;
  nbkg1.setVal(true_nbkg1) ;
  nbkg2.setVal(true_nbkg2) ;
  nbkg3.setVal(true_nbkg3) ;
  nbkg4.setVal(true_nbkg4) ;
  //nev=(int) (true_nsig+true_nbkg1) ;
  //nev=(int) (true_nsig+true_nbkg1+true_nbkg2) ;
  //nev=(int) (true_nsig+true_nbkg1+true_nbkg2+true_nbkg3) ;
  nev=(int) (true_nsig+true_nbkg1+true_nbkg2+true_nbkg3+true_nbkg4) ; // p2e+Michel+acc+dif+prompt

  nevtot=nev ;

  //RooAddPdf tmodel("tmodel","tmodel",RooArgList(t_p2e,t_mich),RooArgList(nsig,nbkg1)) ; 
  //RooAddPdf tmodel("tmodel","tmodel",RooArgList(t_p2e,t_mich,t_acc),RooArgList(nsig,nbkg1,nbkg2)) ;
  //RooAddPdf tmodel("tmodel","tmodel",RooArgList(t_p2e,t_mich,t_acc,t_mu),RooArgList(nsig,nbkg1,nbkg2,nbkg3)) ;
  RooAddPdf tmodel("tmodel","tmodel",RooArgList(t_p2e,t_mich,t_acc,t_mu,t_prom),RooArgList(nsig,nbkg1,nbkg2,nbkg3,nbkg4)) ;

  // calculate the execution time
   Double_t time_exec ;
  // if (n == 1) time_exec = 
  

  // S e t u p   e n e r g y   m o d e l s
  // --------------------------------------

   RooRealVar e("e","e",0,100,"MeV") ;                                     // CsI energy scale: 0-100 MeV

  // p2e CsI line shape
  RooRealVar peak_p2e("peak_p2e","peak of p2e",66.23,1,100) ;          // p2e peak position
  RooRealVar width_p2e("width_p2e","width of p2e",3.49,0.5,10) ;       // p2e width, SD
  RooRealVar tail_p2e("tail_p2e","tail of p2e",-2.18,-5,5) ;           // p2e low-energy (negative) tail
  RooRealVar tailcb_p2e("tailcb_p2e","tailcb of p2e",1.085,0.2,5) ;    // p2e low-energy (negative) tail
  RooRealVar n_p2e("n_p2e","exponent of p2e",7.34,0,10) ;              // p2e exponent in CB function 
  peak_p2e.setConstant(kTRUE) ;
  width_p2e.setConstant(kTRUE) ;
  tail_p2e.setConstant(kTRUE) ;
  tailcb_p2e.setConstant(kTRUE) ;
  n_p2e.setConstant(kTRUE) ;

  // 1st option: modified Gaussian p.d.f in terms of e, peak_p2e, width_p2e, tail_p2e
  //RooNovosibirsk e_p2e("e_p2e","modified gaussian PDF",e,peak_p2e,width_p2e,tail_p2e) ;  
  // 2nd option: Crystal Ball function, already smeared
  RooCBShape e_p2es("e_p2es","crystal ball PDF",e,peak_p2e,width_p2e,tailcb_p2e,n_p2e) ;

  // CsI Response Function, Gaussian
  RooRealVar cmg("cmg","cmg",0,-0.01,0.01) ;   // mean
  RooRealVar csg("csg","csg",2.0,0.5,5) ;  // width, SD
  RooGaussian cgauss("cgauss","cgauss",e,cmg,csg) ;  
  // Define sampling frequency
  //e.setBins(10000,"fft") ;
  // Construct convolution, smeared modified Gaussian p.d.f for p2e line shape
  //RooFFTConvPdf e_p2es("e_p2es","novosibirsk*gauss",e,e_p2e,cgauss) ;

  // Michel energy CsI line shape
  // empirical mp1, mp2, mp3, mp4 (1, 39, 48 and 3.5 MeV) 
  RooRealVar mp1("mp1","mp1",1.48,0.1,2) ;
  mp1.setConstant(kTRUE) ;
  RooRealVar mp2("mp2","mp2",46.08,38,50) ;
  mp2.setConstant(kTRUE) ; 
  RooRealVar mp3("mp3","mp3",50.22,47,55) ;
  mp3.setConstant(kTRUE) ; 
  RooRealVar mp4("mp4","mp4",4.07,3,5) ;
  mp4.setConstant(kTRUE) ; 
  RooGenericPdf e_mich("e_mich","abs((e-mp1)**2*(3-2*e/mp2)*(1-TMath::Erf((e-mp3)/mp4)))",RooArgSet(e,mp1,mp2,mp3,mp4)) ;

  // Accidental CsI energy spectrum
  RooGenericPdf e_acc("e_acc","abs((e-mp1)**2*(3-2*e/mp2)*(1-TMath::Erf((e-mp3)/mp4)))",RooArgSet(e,mp1,mp2,mp3,mp4)) ;

  // DIF CsI energy spectrum
  RooGenericPdf e_dif("e_dif","abs((e-mp1)**2*(3-2*e/mp2)*(1-TMath::Erf((e-mp3)/mp4)))",RooArgSet(e,mp1,mp2,mp3,mp4)) ;


  RooRealVar pp1("pp1","pp1",100,1,500) ;
  mp4.setConstant(kTRUE) ;
  // Prompt (proton) CsI energy spectrum
  RooGenericPdf e_prom("e_prom","pp1",RooArgSet(e,pp1)) ;

  // Target energy lineshape
  RooRealVar et("et","et",0,15,"MeV") ;  // target energy scale: 0-15 MeV

  RooRealVar tgt_mean("tgt_mean","target mean energy",3,0,10) ;
  RooRealVar tgt_sigma("tgt_sigma","target energy sigma",0.5,0.1,10) ;
  // Build target energy p.d.f in terms of tgt_mean, tgt_sigma
  RooLandau e_tgt("e_tgt","target energy PDF",et,tgt_mean,tgt_sigma) ;  
  // Target Response Function
  RooRealVar emg("emg","emg",0) ;
  RooRealVar esg("esg","esg",0.5,0.1,3) ;
  RooGaussian egauss("egauss","egauss",et,emg,esg) ;
  // Define sampling frequency
  et.setBins(10000,"fft") ;
  // Construct convolution: option 1: FFT, option 2: numerical convolution
  RooFFTConvPdf e_tgts("e_tgts","landau*gauss",et,e_tgt,egauss) ;
  //RooNumConvPdf e_tgts("e_tgts","landau*gauss",et,e_tgt,tgauss) ;

  // composite CsI energy line shape model: p2e+Michel+accidental+DIF+prompt+scattered e+

  //RooAddPdf emodel("emodel","emodel",RooArgList(e_p2es,e_mich),RooArgList(nsig,nbkg1)) ;
  //RooAddPdf emodel("emodel","emodel",RooArgList(e_p2es,e_mich,e_acc),RooArgList(nsig,nbkg1,nbkg2)) ;
  //RooAddPdf emodel("emodel","emodel",RooArgList(e_p2es,e_mich,e_acc,e_dif),RooArgList(nsig,nbkg1,nbkg2,nbkg3)) ;
  RooAddPdf emodel("emodel","emodel",RooArgList(e_p2es,e_mich,e_acc,e_dif,e_prom),RooArgList(nsig,nbkg1,nbkg2,nbkg3,nbkg4)) ;

  // S e t u p   v e r t e x   m o d e l s
  // --------------------------------------

  // axial (z) stop distribution, p2e z decay vertex, x,y: +-14 mm, z:-30+15 mm
  RooRealVar z("z","z",-100,25,"mm") ;
  RooRealVar x("x","x",-15,15,"mm") ;
  RooRealVar y("y","y",-15,15,"mm") ;

  // p2e
  RooRealVar zm_p2e("zm_p2e","zm_p2e",1,-1,14) ;
  RooRealVar zs1_p2e("zs1_p2e","zs1_p2e",6.5,0.1,10) ;
  RooRealVar zs2_p2e("zs2_p2e","zs2_p2e",2.5,0.1,10) ;
  zm_p2e.setConstant(kTRUE) ; 
  zs1_p2e.setConstant(kTRUE) ; 
  zs2_p2e.setConstant(kTRUE) ; 

  // p2e z vertices
  //RooAbsPdf* z_p2=RooClassFactory::makePdfInstance("z_p2e",
  //"0.826*exp(-0.5*(z-zm_p2e)*(z-zm_p2e)/(zs1_p2e*zs1_p2e))+exp(-0.5*(z-zm_p2e)*(z-zm_p2e)/(zs2_p2e*zs2_p2e))",RooArgSet(z,zm_p2e,zs1_p2e,zs2_p2e)) ;
  //RooClassFactory::makePdf("z_p2ePdf","z,zm_p2e,zs1_p2e,zs2_p2e") ;
  //RooClassFactory::makePdf("z_michPdf","z,zm_mich,zs1_mich,zs2_mich") ;

  // mich
  RooRealVar zm_mich("zm_mich","zm_mich",1,-1,14) ;
  RooRealVar zs1_mich("zs1_mich","zs1_mich",7.0,0.1,10) ;
  RooRealVar zs2_mich("zs2_mich","zs2_mich",3.0,0.1,10) ;
  zm_mich.setConstant(kTRUE) ; 
  zs1_mich.setConstant(kTRUE) ; 
  zs2_mich.setConstant(kTRUE) ; 

  // Michel z vertices
  //RooAbsPdf* z_mic=RooClassFactory::makePdfInstance("z_mich",
  //"0.826*exp(-0.5*(z-zm_mich)*(z-zm_mich)/(zs1_mich*zs1_mich))+exp(-0.5*(z-zm_mich)*(z-zm_mich)/(zs2_mich*zs2_mich))",RooArgSet(z,zm_mich,zs1_mich,zs2_mich)) ;

  // acc
  RooRealVar zm_acc("zm_acc","zm_acc",1,-1,14) ;
  RooRealVar zs1_acc("zs1_acc","zs1_acc",7.0,0.1,10) ;
  RooRealVar zs2_acc("zs2_acc","zs2_acc",3.0,0.1,10) ;
  zm_acc.setConstant(kTRUE) ; 
  zs1_acc.setConstant(kTRUE) ; 
  zs2_acc.setConstant(kTRUE) ; 

  // accidental z vertices
  //RooAbsPdf* z_ac=RooClassFactory::makePdfInstance("z_acc",
  //"0.826*exp(-0.5*(z-zm_acc)*(z-zm_acc)/(zs1_acc*zs1_acc))+exp(-0.5*(z-zm_acc)*(z-zm_acc)/(zs2_acc*zs2_acc))",RooArgSet(z,zm_acc,zs1_acc,zs2_acc)) ;

  // Input histograms for z_dif pdf
  TFile input("/net/group1/data2/ef2p/pen08/maxlike/histos/z_dif.root"); 
  TH1F *hz_dif = (TH1F *) input.Get("z_dif");
  RooDataHist z_dif_data("z_dif_data","histogram of z dif",z,hz_dif) ;
  RooHistPdf z_dif("z_dif","z_dif G4 data",z,z_dif_data) ;

  // Input histograms for z_prom pdf
  TH1F *hz_prom = (TH1F *) input.Get("z_dif");
  RooDataHist z_prom_data("z_prom_data","histogram of z prom",z,hz_prom) ;
  RooHistPdf z_prom("z_prom","z_prom G4 data",z,z_prom_data) ;

  gROOT->ProcessLineSync(".L RooZ_p2ePdf.cxx+") ; // compile and load
  gROOT->ProcessLineSync(".L RooZ_michPdf.cxx+") ; // compile and load
  gROOT->ProcessLineSync(".L RooZ_accPdf.cxx+") ; // compile and load
  RooZ_p2ePdf z_p2e("z_p2e","z_p2e",z,zm_p2e,zs1_p2e,zs2_p2e) ;
  RooZ_michPdf z_mich("z_mich","z_mich",z,zm_mich,zs1_mich,zs2_mich) ;
  RooZ_accPdf z_acc("z_acc","z_acc",z,zm_acc,zs1_acc,zs2_acc) ;

  //RooAddPdf zmodel("zmodel","zmodel",RooArgList(z_p2e,z_mich),RooArgList(nsig,nbkg1)) ;
  //RooAddPdf zmodel("zmodel","zmodel",RooArgList(z_p2e,z_mich,z_acc),RooArgList(nsig,nbkg1,nbkg2)) ;
  //RooAddPdf zmodel("zmodel","zmodel",RooArgList(z_p2e,z_mich,z_acc,z_dif),RooArgList(nsig,nbkg1,nbkg2,nbkg3)) ;
  RooAddPdf zmodel("zmodel","zmodel",RooArgList(z_p2e,z_mich,z_acc,z_dif,z_prom),RooArgList(nsig,nbkg1,nbkg2,nbkg3,nbkg4)) ;
  
  // lateral (x,y) stop distribution
  RooRealVar xm_p2e("xm_p2e","xm_p2e",0,-15,15) ;
  RooRealVar xs_p2e("xs_p2e","xs_p2e",2,0.1,5) ;
  RooGaussian x_p2e("x_p2e","x_p2e",x,xm_p2e,xs_p2e) ;     // p2e and Michel x vertices
  RooGaussian x_mich("x_mich","x_mich",x,xm_p2e,xs_p2e) ;

  RooRealVar ym_p2e("ym_p2e","ym_p2e",0,-15,15) ;
  RooRealVar ys_p2e("ys_p2e","ys_p2e",2,0.1,5) ;
  RooGaussian y_p2e("y_p2e","y_p2e",y,ym_p2e,ys_p2e) ;     // p2e and Michel y vertices
  RooGaussian y_mich("y_mich","y_mich",y,ym_p2e,ys_p2e) ;

  //RooProdPdf xy_p2e("xy_p2e","xy_p2e",RooArgSet(x_p2e,y_p2e)) ;       // p2e and Michel 2D xy vertices 
  //RooProdPdf xy_mich("xy_mich","xy_mich",RooArgSet(x_mich,y_mich)) ; 

  //RooProdPdf xyz_p2e("xyz_p2e","xyz_p2e",RooArgSet(x_p2e,y_p2e,z_p2e)) ;      // p2e and Michel 3D xy vertices 
  //RooProdPdf xyz_mich("xyz_mich","xyz_mich",RooArgSet(x_mich,y_mich,z_mich)) ;

  //RooAddPdf zmodel("zmodel","zmodel",RooArgList(z_p2e,z_mich),RooArgList(nsig,nbkg1)) ;  

  // vertex distribution model: p2e+Michel+accidental+DIF

  // composite multi-dimensional model

  RooProdPdf t_e_model("t_e_model","t_e_model",RooArgSet(tmodel,emodel)) ;
  //RooProdPdf t_e_v_model("t_e_v_model","t_e_v_model",RooArgSet(tmodel,emodel,zmodel)) ;

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

  switch (n) {

  case 1:

  gROOT->ProcessLine(".! date") ;

  // Sample nev events in t from t_p2e, including time offset and resolution
  RooDataSet* t_p2e_data=t_p2e.generate(t,nevtot) ;

  // Fit t_p2e to t_p2e_data
  t_p2e.fitTo(*t_p2e_data) ;

  // Construct and plot frame in t
  RooPlot* t_p2e_frame = t.frame(Title("p2e decay time p.d.f.")) ;

  // Plot and draw p2e decay time t_p2e_data and fit
  t_p2e_data->plotOn(t_p2e_frame) ;
  t_p2e.plotOn(t_p2e_frame) ;
  t_p2e_frame->Draw() ; 

  Double_t chi2_fit=t_p2e_frame->chiSquare() ;
  printf("\n\n CHI2 VALUE: %.2f\n\n",chi2_fit) ;

  gROOT->ProcessLine(".! date") ;

  break ;
  // -----------------------

  case 2:

  gROOT->ProcessLine(".! date") ;

  // Sample nev events in t from t_mu, including time offset and resolution
  RooDataSet* t_mu_data=t_mu.generate(t,nevtot) ;

  // Fit t_p2e to t_p2e_data
  t_mu.fitTo(*t_mu_data) ;

  // Construct and plot frame in t
  RooPlot* t_mu_frame = t.frame(Title("mu+ decay time p.d.f.")) ;

  // Plot and draw mu+ decay time t_mu_data and fit
  t_mu_data->plotOn(t_mu_frame) ;
  t_mu.plotOn(t_mu_frame) ;
  t_mu_frame->Draw() ; 

  Double_t chi2_fit=t_mu_frame->chiSquare() ;
  printf("\n\n CHI2 VALUE: %.2f\n\n",chi2_fit) ;

  gROOT->ProcessLine(".! date") ;

  break ;
  // -----------------------

  case 3:

  gROOT->ProcessLine(".! date") ;

  // Sample nev events in t from t_mich
  RooDataSet* t_mich_data=t_mich.generate(t,nevtot) ;

  // Fit t_p2e to t_p2e_data
  t_mich.fitTo(*t_mich_data) ;

  // Construct and plot frame in t
  RooPlot* t_mich_frame = t.frame(Title("pi+->mu+->e+ decay time p.d.f.")) ;

  // Plot and draw Michel decay time t_mich_data and fit
  t_mich_data->plotOn(t_mich_frame) ;
  t_mich.plotOn(t_mich_frame) ;
  t_mich_frame->Draw() ; 

  Double_t chi2_fit=t_mich_frame->chiSquare() ;
  printf("\n\n CHI2 VALUE: %.2f\n\n",chi2_fit) ;

  gROOT->ProcessLine(".! date") ;

  break ;
  // -----------------------

  case 4:

   gStyle->SetPadLeftMargin(0.28);
   gStyle->SetPadRightMargin(0.08);
   gStyle->SetTitleXOffset(1.1);
   gStyle->SetTitleYOffset(1.3);

    // timing the generation-and-fit
   gROOT->ProcessLine(".! date") ;

  // Construct and plot frame in t
  //RooPlot* tmodel_frame = t.frame(Title("p2e+Michel+accidentals+DIF+Prompt decay time p.d.f.")) ;

  RooRandom::randomGenerator()->SetSeed(iseed) ; 

  if (niter == 0) { RooDataSet* tmodel_data=tmodel.generate(t,nevtot) ; TFile *ft  = new TFile("tmodel_hist.root","RECREATE"); 
    TH1* ht280=tmodel_data->createHistogram("ht280",t,Binning(tbin,-80,200)) ;
    TH1* ht2800=tmodel_data->createHistogram("ht2800",t,Binning(10*tbin,-80,200)) ;
    TH1* ht28000=tmodel_data->createHistogram("ht28000",t,Binning(100*tbin,-80,200)) ;
    TH1* ht280000=tmodel_data->createHistogram("ht280000",t,Binning(1000*tbin,-80,200)) ;
    //for (int k=1; k<tbin+1; k++) { out280 << ht280->GetBinContent(k) << " " << ht280->GetBinError(k) << endl; }
    //for (int k=1; k<10*tbin+1; k++) { out2800 << ht2800->GetBinContent(k) << " " << ht2800->GetBinError(k) << endl; }
    ft->Write() ; ft->Close() ; return ; }

  // Sample nev events in t from tmodel
    if ( tbin == 0 ) {
      RooDataSet* tmodel_data=tmodel.generate(t,nevtot) ;
    // Fit t_p2e to t_p2e_data
    //tmodel.fitTo(*tmodel_data,Extended(),Save()) ;
    tmodel.fitTo(*tmodel_data,Save()) ;
    // Plot and draw Michel decay time t_mich_data and fit
    tmodel_data->plotOn(tmodel_frame) ; }
    // Alternative: binned event sample in t from tmodel, 1st set numkber of bins
    if ( tbin > 0 ) {
    t.setBins(tbin) ; 
    if (niter == 1 ) RooDataHist* tmodel_hist = tmodel.generateBinned(t,nevtot) ;                                // tmodel == tmodel_hist
    if (niter == 2 ) { TFile *ft  = new TFile("tmodel_hist_1.root"); TH1* ht = (TH1*) gDirectory->Get("ht2800__t"); // tmodel != tmodel_hist
      //for (int k=1; k<tbin+1; k++) out1 << ht->GetBinContent(k) << " " << ht->GetBinError(k) << endl;
      RooDataHist* tmodel_hist= new RooDataHist("tmodel_hist","tmodel_hist",t,Import(*ht)) ; }   // _0 3E8 evts, _1 1.875E8 evts
    tmodel.fitTo(*tmodel_hist,Extended(kTRUE)) ;                                              //    _0 8e7 p2e,  _1 5E7 p2e

    // Construct and plot frame in t
    RooPlot* tmodel_frame = t.frame(Title("p2e+Michel+accidentals+DIF+Prompt decay time p.d.f.")) ;

    tmodel_hist->plotOn(tmodel_frame) ; }
    tmodel.plotOn(tmodel_frame,LineColor(kBlack)) ;
    tmodel.plotOn(tmodel_frame,Components(t_acc),LineStyle(kDashed),LineColor(kRed)) ;
    tmodel.plotOn(tmodel_frame,Components(t_mich),LineStyle(kDashed),LineColor(kGreen)) ;
    tmodel.plotOn(tmodel_frame,Components(t_mu),LineStyle(kDashed),LineColor(kBlue)) ;
    tmodel.plotOn(tmodel_frame,Components(t_prom),LineStyle(kDashed),LineColor(6)) ;
    tmodel.plotOn(tmodel_frame,Components(t_p2e),LineStyle(kDashed),LineColor(28)) ;

    lat.DrawLatex(35,0.034*true_nsig,Form("tmodel: %.3f ns/bin, %.0e p2e events",(float) 280/tbin,true_nsig));
    lat.DrawLatex(35,0.032*true_nsig,Form("Fit Error: (N_{fit}-N_{true})/N_{true} (p2e)= %.3f perc.",fabs(100*(nsig.getVal()-true_nsig)/true_nsig)));
    lat.DrawLatex(35,0.030*true_nsig,Form("#sigma_{t}=1 ns, t_{off}=0.0 ns"));

    tmodel_frame->Draw() ; 

    Double_t BR=nsig.getVal()/(nbkg1.getVal()+nbkg2.getVal()) ;
    printf (" \n\n True p2e: %f, Fitted p2e: %f, p2e Event# Uncertainty (perc.): %f \n\n",true_nsig,nsig.getVal(),
            fabs(100*(nsig.getVal()-true_nsig)/true_nsig)) ;

    Double_t chi2_fit=0 ; //=tmodel_frame->chiSquare() ;
    //printf("\n\n CHI2 VALUE: %.2f\n\n",chi2_fit) ;

    RooHist* histo = (RooHist*) tmodel_frame->findObject("tmodel_hist") ;
    RooCurve* func = (RooCurve*) tmodel_frame->findObject("tmodel") ;
    for (Int_t i=0 ; i<histo->GetN() ; i++) {
    Double_t xdata,ydata ;
    histo->GetPoint(i,xdata,ydata) ;
    Double_t yfunc = curve->interpolate(xdata) ;
    if (ydata>0) chi2_fit=chi2_fit + (ydata-yfunc)*(ydata-yfunc)/ydata ; }
    printf("\n\n CHI2 VALUE: %.2f\n\n",chi2_fit) ;

    gROOT->ProcessLine(".! date") ;

    // comment out for condor_q
    if (tbin > 0)   c->Print(Form("figsimple/tmodel%db_%.0eev.pdf",tbin,true_nsig));
    if (tbin == 0)  c->Print(Form("figsimple/tmodelunb_%.0eev.pdf",true_nsig));

    //file->Write();
  
  out << nsig.getVal() << " " << nsig.getError() << endl; 

  break ;
  // -----------------------

  case 5:

   gStyle->SetPadLeftMargin(0.28);
   gStyle->SetPadRightMargin(0.08);
   gStyle->SetTitleXOffset(1.1);
   gStyle->SetTitleYOffset(1.3);

  // Sample nev events in e from e_p2es
  RooDataSet* e_p2es_data=e_p2es.generate(e,nevtot) ;

  // Fit e_p2e to e_p2es_data
  e_p2es.fitTo(*e_p2es_data) ;

  // Construct and plot frame in e
  RooPlot* e_p2es_frame = e.frame(Title("p2e energy line shape p.d.f.")) ;

  // Plot and draw p2e line shape e_p2es_data and fit
  e_p2es_data->plotOn(e_p2es_frame) ;
  e_p2es.plotOn(e_p2es_frame) ;
  e_p2es_frame->Draw() ; 

  break ;
  // -----------------------

  case 6:

   gStyle->SetPadLeftMargin(0.28);
   gStyle->SetPadRightMargin(0.08);
   gStyle->SetTitleXOffset(1.1);
   gStyle->SetTitleYOffset(1.3);

  // Sample 10000 events in e from e_mich
  RooDataSet* e_mich_data=e_mich.generate(e,nevtot) ;

  // Fit e_mich to data histogram
  e_mich.fitTo(*e_mich_data) ;

  // Construct and plot frame in e
  RooPlot* e_mich_frame = e.frame(Title("Michel energy line shape p.d.f.")) ;

  // Plot and draw Michel line shape e_mich_data and fit
  e_mich_data->plotOn(e_mich_frame) ;
  e_mich.plotOn(e_mich_frame) ;
  e_mich_frame->Draw() ; 

  break ;
  // -----------------------

  case 7:

   gStyle->SetPadLeftMargin(0.28);
   gStyle->SetPadRightMargin(0.08);
   gStyle->SetTitleXOffset(1.1);
   gStyle->SetTitleYOffset(1.3);

    // timing the generation-and-fit
  gROOT->ProcessLine(".! date") ;

  // Construct and plot frame in t
  RooPlot* emodel_frame = e.frame(Title("p2e+Michel+accidentals+DIF+Prompt CsI energy p.d.f.")) ;

  RooRandom::randomGenerator()->SetSeed(iseed) ; 

  // Sample nev events in e from emodel
    if ( ebin == 0 ) {
    RooDataSet* emodel_data=emodel.generate(e,nevtot) ;
    // Fit e_p2e to e_p2e_data
    emodel.fitTo(*emodel_data,Save()) ; 
    // Plot and draw Michel decay time t_mich_data and fit
    emodel_data->plotOn(emodel_frame) ; }
    // Alternative: binned event sample in t from tmodel, 1st set numkber of bins
    if ( ebin > 0 ) {
    e.setBins(ebin) ; 
    RooDataHist* emodel_hist = emodel.generateBinned(e,nevtot) ; 
    emodel.fitTo(*emodel_hist,Save()) ; 
    emodel_hist->plotOn(emodel_frame) ; }

    emodel.plotOn(emodel_frame,LineColor(kBlack)) ;
    emodel.plotOn(emodel_frame,Components(e_acc),LineStyle(kDashed),LineColor(kRed)) ;
    emodel.plotOn(emodel_frame,Components(e_mich),LineStyle(kDashed),LineColor(kGreen)) ;
    emodel.plotOn(emodel_frame,Components(e_dif),LineStyle(kDashed),LineColor(kBlue)) ;
    emodel.plotOn(emodel_frame,Components(e_prom),LineStyle(kDashed),LineColor(6)) ;
    emodel.plotOn(emodel_frame,Components(e_p2es),LineStyle(kDashed),LineColor(28)) ;

    emodel_frame->Draw() ; 

    lat.DrawLatex(50,0.060*true_nsig,Form("emodel: %.2f MeV/bin, %.0e p2e events",(float) 200/ebin,true_nsig));
    lat.DrawLatex(50,0.057*true_nsig,Form("Fit Error: (N_{fit}-N_{true})/N_{true} (p2e)= %.3f perc.",fabs(100*(nsig.getVal()-true_nsig)/true_nsig)));
    lat.DrawLatex(50,0.054*true_nsig,Form("#sigma_{E}=3 MeV, E_{off}=0.0 MeV"));

    Double_t BR=nsig.getVal()/(nbkg1.getVal()+nbkg2.getVal()) ;
    //printf (" \n\n p2e Event# Uncertainty (perc.): %f \n\n",fabs(100*(nsig.getVal()-true_nsig)/true_nsig)) ;

    gROOT->ProcessLine(".! date") ;

    // comment out for condor_q
    if (ebin > 0)   c->Print(Form("figsimple/emodel%db_%.0eev.pdf",ebin,true_nsig));
    if (ebin == 0)  c->Print(Form("figsimple/emodelunb_%.0eev.pdf",true_nsig));

    //file->Write();
  
  out << nsig.getVal() << " " << nsig.getError() << endl;

  break ;
  // -----------------------

  case 8:
  // Sample 10000 events in et from e_tgts
  RooDataSet* e_tgts_data=e_tgts.generate(et,10000) ;

  // Fit me_tgts to data
  e_tgts.fitTo(*e_tgts_data) ;

  // Construct and plot frame in et
  RooPlot* e_tgts_frame = et.frame(Title("target energy p.d.f.")) ;

  // Plot and draw target energy in frame 
  e_tgts_data->plotOn(e_tgts_frame) ;
  e_tgt.plotOn(e_tgts_frame) ;
  e_tgts.plotOn(e_tgts_frame,LineStyle(kDashed),LineColor(kRed)) ;
  e_tgts_frame->Draw() ; 

  break ;
  // -----------------------

  case 9:
  // Sample 10000 events in z from z_p2e
  RooDataSet* z_p2e_data=z_p2e.generate(z,10000) ;

  // Fit z_p2e to data
  z_p2e.fitTo(*z_p2e_data) ;
 
  // Construct and plot frame in z
  RooPlot* z_p2e_frame = z.frame(Title("z p2e vertex p.d.f.")) ;

  // Plot and draw z vertex in frame 
  z_p2e_data->plotOn(z_p2e_frame) ;
  z_p2e.plotOn(z_p2e_frame) ;
  z_p2e_frame->Draw() ; 

  break ;
  // -----------------------

  case 10:
  // Sample 10000 events in z from z_mich
  RooDataSet* z_mich_data=z_mich.generate(z,10000) ;

  // Fit z_mich to data
  z_mich.fitTo(*z_mich_data) ;
 
  // Construct and plot frame in z
  RooPlot* z_mich_frame = z.frame(Title("z Michel vertex p.d.f.")) ;

  // Plot and draw z vertex in frame 
  z_mich_data->plotOn(z_mich_frame) ;
  z_mich.plotOn(z_mich_frame) ;
  z_mich_frame->Draw() ; 

  break ;
  // -----------------------

  case 12:
  // Sample 10000 events in x,y from xy_p2e
  RooDataSet* xy_p2e_data=xy_p2e.generate(RooArgSet(x,y),10000) ;

  // Fit xy_p2e to data
  xy_p2e.fitTo(*xy_p2e_data) ;
 
  // create 2D histogram in x,y
  TH2* xy_p2e_pdf=xy_p2e.createHistogram("x,y",20,20) ;

  // Construct and plot frame in x,y
  RooPlot* xy_p2e_frame = x.frame(Title("(x,y) p2e vertex p.d.f.")) ;

  // Plot and draw (x,y) pdf in frame 
  //xy_p2e_data->plotOn(xy_p2e_frame) ;
  //xy_p2e.plotOn(xy_p2e_frame) ;
  //xy_p2e_frame->Draw() ; 
  xy_p2e_pdf->Draw("box") ;

  break ;
  // -----------------------

  case 13:
  // Sample 10000 events in x,y from xy_mich
  RooDataSet* xy_mich_data=xy_mich.generate(RooArgSet(x,y),10000) ;

  // Fit xy_mich to data
  xy_mich.fitTo(*xy_mich_data) ;
 
  // create 2D histogram in x,y
  TH2* xy_mich_pdf=xy_mich.createHistogram("x,y",20,20) ;

  // Construct and plot frame in x,y
  RooPlot* xy_mich_frame = x.frame(Title("(x,y) Michel vertex p.d.f.")) ;

  // Plot and draw (x,y) pdf in frame 
  //xy_mich_data->plotOn(xy_mich_frame) ;
  //xy_mich.plotOn(xy_mich_frame) ;
  //xy_mich_frame->Draw() ; 
  xy_mich_pdf->Draw("lego") ;

  break ;
  // -----------------------

  case 12:
  // Sample 10000 events in x,y from xy_p2e
  RooDataSet* xy_p2e_data=xy_p2e.generate(RooArgSet(x,y),10000) ;

  // Fit xy_p2e to data
  xy_p2e.fitTo(*xy_p2e_data) ;
 
  // create 2D histogram in x,y
  TH2* xy_p2e_pdf=xy_p2e.createHistogram("x,y",20,20) ;

  // Construct and plot frame in x,y
  RooPlot* xy_p2e_frame = x.frame(Title("(x,y) p2e vertex p.d.f.")) ;

  // Plot and draw (x,y) pdf in frame 
  //xy_p2e_data->plotOn(xy_p2e_frame) ;
  //xy_p2e.plotOn(xy_p2e_frame) ;
  //xy_p2e_frame->Draw() ; 
  xy_p2e_pdf->Draw("box") ;

  break ;
  // -----------------------

  case 101:

   gStyle->SetPadLeftMargin(0.50);
   gStyle->SetPadRightMargin(0.02);
   gStyle->SetTitleXOffset(1.1);
   gStyle->SetTitleYOffset(1.5);

    // timing the generation-and-fit
  gROOT->ProcessLine(".! date") ;

  // Construct and plot frame in t
  RooPlot* zmodel_frame = z.frame(Title("p2e+Michel+accidentals+DIF+Prompt z vertex p.d.f.")) ;

  RooRandom::randomGenerator()->SetSeed(iseed) ; 

  // Sample nev events in e from emodel
    if ( zbin == 0 ) {
    RooDataSet* zmodel_data=zmodel.generate(z,nevtot) ;
    // Fit z_p2e to z_p2e_data
    zmodel.fitTo(*zmodel_data,Extended(),Save()) ; 
    // Plot and draw Michel decay time t_mich_data and fit
    zmodel_data->plotOn(zmodel_frame) ; }
    // Alternative: binned event sample in t from tmodel, 1st set numkber of bins
    if ( zbin > 0 ) {
    z.setBins(tbin) ; 
    RooDataHist* zmodel_hist = zmodel.generateBinned(z,nevtot) ; 
    zmodel.fitTo(*zmodel_hist,Save()) ; 
    zmodel_hist->plotOn(zmodel_frame) ; }

    zmodel.plotOn(zmodel_frame,LineColor(kBlack)) ;
    zmodel.plotOn(zmodel_frame,Components(z_p2e),LineStyle(kDashed),LineColor(28)) ;
    zmodel.plotOn(zmodel_frame,Components(z_mich),LineStyle(kDashed),LineColor(kGreen)) ;
    zmodel.plotOn(zmodel_frame,Components(z_acc),LineStyle(kDashed),LineColor(kRed)) ;
    zmodel.plotOn(zmodel_frame,Components(z_dif),LineStyle(kDashed),LineColor(kBlue)) ;
    zmodel.plotOn(zmodel_frame,Components(z_prom),LineStyle(kDashed),LineColor(6)) ;

    zmodel_frame->Draw() ; 

    lat.DrawLatex(50,0.060*true_nsig,Form("zmodel: %.2f mm/bin, %.0e p2e events",(float) 200/zbin,true_nsig));
    lat.DrawLatex(50,0.057*true_nsig,Form("Fit Error: (N_{fit}-N_{true})/N_{true} (p2e)= %.3f perc.",fabs(100*(nsig.getVal()-true_nsig)/true_nsig)));
    lat.DrawLatex(50,0.054*true_nsig,Form("#sigma_{E}=3 mm, z_{off}=0.0 mm"));

    Double_t BR=nsig.getVal()/(nbkg1.getVal()+nbkg2.getVal()) ;
    //printf (" \n\n p2e Event# Uncertainty (perc.): %f \n\n",fabs(100*(nsig.getVal()-true_nsig)/true_nsig)) ;

    gROOT->ProcessLine(".! date") ;

    // comment out for condor_q
    if (zbin > 0)   c->Print(Form("figsimple/zmodel%db_%.0eev.pdf",zbin,true_nsig));
    if (zbin == 0)  c->Print(Form("figsimple/zmodelunb_%.0eev.pdf",true_nsig));

    //file->Write();
  
  out << nsig.getVal() << " " << nsig.getError() << endl;

  break ;
  // -----------------------

  case 100:

   gStyle->SetPadLeftMargin(0.28);
   gStyle->SetPadRightMargin(0.08);
   gStyle->SetTitleXOffset(1.1);
   gStyle->SetTitleYOffset(1.3);

   RooRandom::randomGenerator()->SetSeed(iseed) ; 

  // timing the generation-and-fit
  gROOT->ProcessLine(".! date") ;

  // Construct and plot frame in t
  RooPlot* tmodel_frame = t.frame(Title("p2e+Michel+accidentals+DIF+prompt e+ Decay Time p.d.f.")) ;
 // Construct and plot frame in e
  RooPlot* emodel_frame = e.frame(Title("p2e+Michel+accidentals+DIF+prompt e+ CsI Energy p.d.f.")) ;

  if ( tbin*ebin == 0 ) {
  // Sample nevtot events in (t,e) from t_e_model
  RooDataSet* t_e_model_data=t_e_model.generate(RooArgSet(t,e),nevtot) ;

  // Fit emodel to emodel_data
  t_e_model.fitTo(*t_e_model_data,Save()) ; }

  if ( tbin*ebin > 0 ) {
    t.setBins(tbin) ; 
    e.setBins(ebin) ; 
    RooDataHist* t_e_model_hist=t_e_model.generateBinned(RooArgSet(t,e),nevtot) ;
    t_e_model.fitTo(*t_e_model_hist,Save()) ; }

  c->Divide(1,2) ;
  c->cd(1) ;

  // Plot and draw Michel decay time t_mich_data and fit
  if ( tbin*ebin == 0 ) t_e_model_data->plotOn(tmodel_frame) ;
  if ( tbin*ebin > 0 )  t_e_model_hist->plotOn(tmodel_frame) ;
  t_e_model.plotOn(tmodel_frame,LineColor(kBlack)) ;
  t_e_model.plotOn(tmodel_frame,Components(t_acc),LineStyle(kDashed),LineColor(kRed)) ;
  t_e_model.plotOn(tmodel_frame,Components(t_mich),LineStyle(kDashed),LineColor(kGreen)) ;
  t_e_model.plotOn(tmodel_frame,Components(t_mu),LineStyle(kDashed),LineColor(kBlue)) ;
  t_e_model.plotOn(tmodel_frame,Components(t_prom),LineStyle(kDashed),LineColor(6)) ;
  t_e_model.plotOn(tmodel_frame,Components(t_p2e),LineStyle(kDashed),LineColor(28)) ;

  tmodel_frame->Draw() ; 

  c->cd(2) ;

  // Plot and draw e+ CsI energy lineshape and fit
  if ( tbin*ebin == 0 ) t_e_model_data->plotOn(emodel_frame) ;
  if ( tbin*ebin > 0 )  t_e_model_hist->plotOn(emodel_frame) ;
  t_e_model.plotOn(emodel_frame) ;
  t_e_model.plotOn(emodel_frame,Components(e_acc),LineStyle(kDashed),LineColor(kRed)) ;
  t_e_model.plotOn(emodel_frame,Components(e_mich),LineStyle(kDashed),LineColor(kGreen)) ;
  t_e_model.plotOn(emodel_frame,Components(e_dif),LineStyle(kDashed),LineColor(kBlue)) ;
  t_e_model.plotOn(emodel_frame,Components(e_prom),LineStyle(kDashed),LineColor(6)) ;
  t_e_model.plotOn(emodel_frame,Components(e_p2es),LineStyle(kDashed),LineColor(28)) ;

  emodel_frame->Draw() ; 

  gROOT->ProcessLine(".! date") ;

    // comment out for condor_q
  if (ebin > 0 && tbin >0)   c->Print(Form("figsimple/t_e_model%db_%db_%.0eev.pdf",tbin,ebin,true_nsig));
  if (ebin == 0 && tbin ==0)  c->Print(Form("figsimple/t_e_modelunb_%.0eev.pdf",true_nsig));


  break ;
  // -----------------------

  }
}