#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include <TLatex.h>
#include "RooRealVar.h"
#include "RooDataSet.h"
#include <TString.h>
#include "RooFormulaVar.h"
#include "RooAddPdf.h"
#include "RooExponential.h"
#include "RooGaussModel.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "RooPlot.h"
#include "TTree.h"
#include "TCut.h"
#include "TMath.h"
#include "TFile.h"
#include "TCanvas.h"
#include "RooHist.h"
#include "RooGenericPdf.h"
#include "RooTruthModel.h"
#include "RooGaussModel.h"
#include "RooDecay.h"
#include "RooProdPdf.h"
#include "RooEffProd.h"
#include <TRandom3.h>
#include "RooRandom.h"
#include "RooPoisson.h"
#include "RooGamma.h"
#include "RooBernstein.h"
#include "RooGaussian.h"
#include "iostream"
#include "fstream"
#include "RooFitResult.h"
#include "RooCBShape.h"
#include "RooDataSet.h"
#include "TH1D.h"
#include "RooConstVar.h"
#include "RooLognormal.h"
#include "TNtupleD.h"
#include "RooStats/ProfileLikelihoodCalculator.h"

using namespace std;
using namespace RooFit;
using namespace RooStats;

void totalpdftoyprofile(){

  int size=0,seed,unsucess=0;

  TNtupleD *nt = new TNtupleD("nt","","taumean:tauerror:pull");
  cout<<"size for fitting and seed "<<endl;
  cin>>size>>seed;
  for(int fc=0;fc<size;fc++){
    cout<<" \n acting on toy:"<<fc<<"\n\n"<<endl;
    
    RooRealVar* mass=new RooRealVar("mass","M_{#mu^{+}#mu^{-}}[ GeV/c^{2}]",4.9,5.9);
    RooRealVar* treco=new RooRealVar("treco","t_{reco}[ps]",1.2,12.5) ;
    RooRealVar* trecoe=new RooRealVar("trecoe","terr[ps]",0.003,0.26);
    
    //comb
    RooRealVar* B1=new RooRealVar("B1", "B_1comb", 0.516, 0.345, 0.687);
    RooFormulaVar* B2=new RooFormulaVar("B2", "B_2comb", "1.-@0", RooArgList(*B1));
    RooBernstein* mass_comb=new RooBernstein("mass_comb", "mass_comb", *mass, RooArgSet(*B1,*B2));

    //Bsmass
    RooRealVar* meanBs=new RooRealVar("meanBs", "mean",5.3642,5.32, 5.40);
    RooRealVar *sigmacbbs=new RooRealVar("sigmacbbs","m1",0.0489);//,0.004,0.2);                                     
    sigmacbbs->setConstant(true);
    RooRealVar *tailcbbs=new RooRealVar("tailcbbs","",1.7617);//,0.1,10.0);                                            
    tailcbbs->setConstant(true);
    RooRealVar *powcbbs=new RooRealVar("powcbbs","",2.12);//,0,10);                                                 
    powcbbs->setConstant(true);
    RooCBShape *cbsbs=new RooCBShape("cbsbs","Signal Lineshape",*mass,*meanBs,*sigmacbbs,*tailcbbs,*powcbbs);
    //Bdmass
    RooRealVar* meanBd=new RooRealVar("meanBd", "mean",5.2772);//,5.25, 5.30);
    RooRealVar* sigmaBd=new RooRealVar("sigmaBd","sigma1 of Gaussian",0.0379);//,0.0,0.4);
    RooGaussian* gausianBd=new RooGaussian("gausianBd","gaus1", *mass, *meanBd, *sigmaBd);
    RooRealVar *sigmacbbd=new RooRealVar("sigmacbbd","m1",0.0619);//,0.004,0.2);
    RooRealVar *tailcbbd=new RooRealVar("tailcbbd","",1.4646);//,0.1,3.0);
    RooRealVar *powcbbd=new RooRealVar("powcbbd","",2.80);//,0,10);
    RooCBShape *cbsbd=new RooCBShape("cbsbd","Signal Lineshape",*mass,*meanBd,*sigmacbbd,*tailcbbd,*powcbbd);
    RooRealVar* frac_gau_CB_bd=new RooRealVar("frac_gau_CB_bd","fraction", 0.556);//, 0, 1);
    RooAddPdf* mass_pdf_bd=new RooAddPdf("mass_pdf_bd","gausian+CB shape", RooArgList(*gausianBd,*cbsbd), *frac_gau_CB_bd);
    
    //semi
    RooRealVar* meansemi=new RooRealVar("meansemi", "mean",4.868,4.7, 5.289);
    RooRealVar* sigmasemi=new RooRealVar("sigmasemi","sigma1 of Gaussian",0.155);//,0.0,1.0);
    RooGaussian* gausiansemi=new RooGaussian("gausiansemi","gaus1", *mass, *meansemi, *sigmasemi);
    
    RooRealVar* meant=new RooRealVar("meant", "mean", 0);
    RooRealVar* sigmat=new RooRealVar("sigmat","sigma",0.072);
    RooGaussModel* gm=new RooGaussModel("gm","gm", *treco,*meant,*sigmat);//,*trecoe);
  
    RooFormulaVar* effFor=new RooFormulaVar("effFor","-0.1808-0.006736*treco+0.2738/(1+exp(-1.2853*treco))",RooArgSet(*treco));
    
    //    RooTruthModel* CM =new RooTruthModel("CM","",*treco);

    //bd lifetime
    RooRealVar*  TauBd=new RooRealVar("TauBd","",1.45,0,100);
    RooDecay* Bkg_Ctaubd =new RooDecay("Bkg_Ctaubd","",*treco,*TauBd,*gm,RooDecay::SingleSided);
    RooFormulaVar* effForch1 = new RooFormulaVar("effForch1","8.95745e-07-4.95401e-08*treco-2.43181e-06/(1+exp(1.12513 *treco))", RooArgSet(*treco) );
    RooEffProd* CtEffBd =new RooEffProd("CtEffBd","",*Bkg_Ctaubd,*effForch1);

    
    RooRealVar*  Tau=new RooRealVar("Tau","",1.70,0.00,10);//Bs lifetime
    RooDecay* sig_Ctau1 =new RooDecay("sig_Ctau1","",*treco,*Tau,*gm,RooDecay::SingleSided);
    //semi lifetime
    RooRealVar*  Tausemi=new RooRealVar("Tausemi","",1.338,0.0,10.0);
    RooDecay* semi_Ctau =new RooDecay("semi_Ctau","",*treco,*Tausemi,*gm,RooDecay::SingleSided);

    // background lifetime                                                                                                    
    RooRealVar*  TauBkg2=new RooRealVar("TauBkg2","",3.1068,0,10);
    RooRealVar*  TauBkg3=new RooRealVar("TauBkg3","",0.2773,0,5);
    RooDecay* Bkg_Ctau3 =new RooDecay("Bkg_Ctau3","",*treco,*TauBkg3,*gm,RooDecay::SingleSided);
    RooDecay* Bkg_Ctau2 =new RooDecay("Bkg_Ctau2","",*treco,*TauBkg2,*gm,RooDecay::SingleSided);    
    RooRealVar* fg1=new RooRealVar("fg1","",0.588,0.0,1.0);
    RooAddPdf* bkg_Ctau=new RooAddPdf("bkg_Ctauf","",RooArgSet(*Bkg_Ctau2,*Bkg_Ctau3),*fg1);
    RooEffProd* CtEffSig =new RooEffProd("CtEffSig","",*sig_Ctau1,*effFor);

    RooRealVar* nbkg=new RooRealVar("nbkg","Background fraction",120,20,55000);
    RooRealVar* Nsig=new RooRealVar("Nsig","signal fraction",30,0,25000);
    RooRealVar* nbd1=new RooRealVar("nbd1","Bd",5,0,6000);
    RooRealVar* nsemi=new RooRealVar("nsemi","semi",50,0,50000);
    //==============
    RooProdPdf* BsPdf=new RooProdPdf ("BsPdf","",RooArgList(*cbsbs,*CtEffSig));
    RooProdPdf*  BdPdf=new RooProdPdf ("BdPdf","",RooArgList(*mass_pdf_bd,*CtEffBd) );
    RooProdPdf*  BkgPdf=new RooProdPdf ("BkgPdf","",RooArgList(*mass_comb,*bkg_Ctau) );
    RooProdPdf*  semiPdf=new RooProdPdf ("semiPdf","",RooArgList(*gausiansemi,*semi_Ctau) );
    
    cout<<"\n Experiment # "<<fc<<"\n"<<endl;
    int seed1=(fc+1)*seed;
    
    RooAddPdf* TotPdf=new RooAddPdf("TotPdf"," Signal + Bkg Pdf",RooArgList(*BsPdf,*BdPdf,*BkgPdf,*semiPdf),RooArgList(*Nsig,*nbd1,*nbkg,*nsemi) );  
    RooLognormal* bdconstraint=new RooLognormal("bdconstraint","bd constraint",*nbd1,RooConst(5),RooConst(2)) ;
    RooGaussian* bdconstrainttau=new RooGaussian("bdconstrainttau","bd constraint",*TauBd,RooConst(1.45),RooConst(0.9)) ;
    RooGaussian* meanbsconstraint=new RooGaussian("meanbsconstraint","bd constraint",*meanBs,RooConst(5.36),RooConst(0.04)) ;
    RooGaussian* B1constraint=new RooGaussian("B1constraint", "constraint",*B1,RooConst(0.51),RooConst(0.2)) ;

    RooArgList constr;
    constr.add(*bdconstraint);
    constr.add(*meanbsconstraint);                                                                                                   
    constr.add(*B1constraint);                                                                                                                  
    constr.add(*bdconstrainttau);

    RooProdPdf* constrpdf=new RooProdPdf("constrpdf","",constr);
    RooProdPdf* model_con=new RooProdPdf("model_con","model with constraint",RooArgSet(*TotPdf,*constrpdf)) ;                                 
    
    RooDataSet* tmp_evt = bdconstraint->generate(RooArgSet(*nbd1),1);                                                                        
    nbd1->setVal(tmp_evt->get(0)->getRealValue(nbd1->GetName()));                                                                          
    delete tmp_evt;                                                                                                                               

    tmp_evt = bdconstrainttau->generate(RooArgSet(*TauBd),1);                                                                    
    TauBd->setVal(tmp_evt->get(0)->getRealValue(TauBd->GetName()));                                                                            
    delete tmp_evt;                                                                                                                              

    RooDataSet* tmp_evtB1 = B1constraint->generate(RooArgSet(*B1),1);                                                                     
    B1->setVal(tmp_evtB1->get(0)->getRealValue(B1->GetName()));                                                                            
    delete tmp_evtB1;                                                                                                                               
    RooDataSet* tmp_evtbs = meanbsconstraint->generate(RooArgSet(*meanBs),1);                                                              
    meanBs->setVal(tmp_evtbs->get(0)->getRealValue(meanBs->GetName()));                                                                     
    delete tmp_evtbs; 
    /*    
    RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
    RooMsgService::instance().setSilentMode(kTRUE);
    RooMsgService::instance().setStreamStatus(1,false);
    RooMsgService::instance().getStream(1).removeTopic(Integration) ;
    RooMsgService::instance().getStream(1).removeTopic(Minimization) ;
    RooMsgService::instance().getStream(1).removeTopic(Fitting) ;
    RooMsgService::instance().getStream(1).removeTopic(NumIntegration) ;
    RooMsgService::instance().getStream(1).removeTopic(Optimization) ;
    RooMsgService::instance().getStream(1).removeTopic(ObjectHandling) ;
    RooMsgService::instance().getStream(1).removeTopic(Eval) ;
    RooMsgService::instance().Print() ;    
    */

    RooRandom::randomGenerator()->SetSeed(seed1);
    RooDataSet* data=model_con->generate(RooArgSet(*mass,*treco),Extended(true));
    data->Print();
    RooFitResult* fit3d=model_con->fitTo(*data,Extended(kTRUE),Save());//,ConditionalObservables(*trecoe));          
    fit3d->Print("v");  
    int stat=fit3d->status();
    int cov=fit3d->covQual();
    if(stat==0&&cov==3){      
      double result[3];
      result[0]=Tau->getVal();
      result[1]=Tau->getError();
      result[2]=(Tau->getVal()-1.70)/Tau->getError();
      nt->Fill(result);

    }
    TotPdf->Delete();
    model_con->Delete();
    data->Delete();

  }
  TFile *f = new TFile("2Dtoy_consB1bd_30_5.root","recreate");
  nt->Write();

}