void Test(){
////////////////////////////////////////////////////////////
// ********** Defining Constants and Variables ********** //
////////////////////////////////////////////////////////////


  std::cout << "*** Defining Constants and Variables ***" << std::endl;

  // ********** Detector Observables ********** //
  ////////////////////////////////////////////////

  // *** B0 Mass ***
  B0m = new RooRealVar("B0m","Reconstructed B0 mass",5346,5529,"MeV"); // Bs toy analysis

  // *** Bs Lifetime ***
  Bst = new RooRealVar("Bst","Reconstructed Bs lifetime",-0.5,10,"ps");

 
  // ********** Used Constants: B0 Mass ********** //
  ///////////////////////////////////////////////////
 
  // *** signal sample ***
  B0mass_Bs = new RooRealVar("B0mass_Bs","B0 Mass of the Bs meson",5366.3,"MeV"); // PDG 06-01-2009
  B0mass_sigma_G = new RooRealVar("B0mass_sigma_G","Width of the Gaussians describing the reconstructed B mass",10,"MeV"); // From Ntuple fit

  // *** Incl background samples ***
  B0mass_lambda = new RooRealVar("B0mass_lambda","Decay constant of the exponential describing the reconstructed B mass",-0.00033,"MeV"); // From Real Data fit


  // ********** Used Constants: Time ********** //
  ////////////////////////////////////////////////

  // *** Lifetime ***
  tau_Bs = new RooRealVar("tau_Bs","Average lifetime for Bs",1.472,"ps"); // PDG 06-01-2009
  dG_Bs  = new RooRealVar("dG_Bs","Lifetime difference for Bs",0.062,"ps^-1"); // HFAG 2009
  dM_Bs  = new RooRealVar("dM_Bs","Mass difference between light and heavy weak-eigenstate for Bs",17.77,"ps^-1");// PDG 06-01-2009

  // *** CP Observables *** //
  ////////////////////////////

  Adir_Bs = new RooRealVar("Adir_Bs","Direct CP asymmetry for Bs", -0.12,-1.5,1.5);
  Amix_Bs = new RooRealVar("Amix_Bs","Mixing induced CP asymmetry for Bs",-0.53,-1.5,1.5);
  ADG_Bs  = new RooRealVar("ADG_Bs" ,"Delta Gamma CP asymmetry for Bs", 0.84,-1.5,1.5);


  // *** Decay equations *** //
  /////////////////////////////
  // *** Mistagging ***
  // |A(t)|_obs^2 = w|Abar(t)|_phys^2 + (1-w)|A(t)|_phys^2
  mistag = new RooRealVar("mistag","Mistag fraction of original B and Bbar",0.3651);//Roadmap for Bd->J/psi phi Table 15
  Int_t tageff = 46; //Roadmap for Bd->J/psi phi Table 15 (45.61%)
  tag = new RooCategory("tag","Flavour tag of the B meson");
  tag->defineType("Bs",1);
  tag->defineType("Bsbar",-1);
  tag->defineType("untag",0);

  // *** Coefficients for the decay equations ***
  // Overall normalization is done automatically
  ACosh_Bs = new RooRealVar("ACosh_Bs","Parameter of the Cosh in decay equation for Bs",1);
  DSinh_Bs = new RooFormulaVar("DSinh_Bs","Parameter of the Sinh in decay equation for Bs","@0",RooArgList(*ADG_Bs));
  CCos_Bs = new RooFormulaVar("CCos_Bs","Parameter of the Cos in decay equation for Bs","@2*@0*(1-2*@1)",RooArgList(*Adir_Bs,*mistag,*tag));
  SSin_Bs = new RooFormulaVar("SSin_Bs","Parameter of the Sin in decay equation for Bs","@2*@0*(1-2*@1)",RooArgList(*Amix_Bs,*mistag,*tag)); 

  // *** Constants for Proper Time Error smearing ***
  resMean = new RooRealVar("resMean","resolution mean",0.,"ps");
  resDecay = new RooRealVar("resDecay","resolution width for signal",0.054,"ps"); // From NTuple fit
  resPrompt_III = new RooRealVar("resPrompt_III","resolution width for prompt Combi I",0.0661,"ps"); // From Real Data fit

  // ********** Number of Events ********** //
  ////////////////////////////////////////////

  // *** number of events for the fitter ***
  nevents_fit_Bs = new RooRealVar("nevents_fit_Bs","Number of Bs events - for fitting",2494);
  nevents_fit_Jpsi = new RooRealVar("nevents_fit_Jpsi","Number of Jpsi events - for fitting",65266);
  Int_t nevents_toy_input = nevents_fit_Bs->getVal() + nevents_fit_Jpsi->getVal();


////////////////////////////////////////////////
// ********** Building B0 Mass PDF ********** //
////////////////////////////////////////////////

  std::cout << "*** Building B0 Mass PDF ***" << std::endl;

  // *** Gauss ***
  B0mass_G_Bs = new RooGaussian("B0mass_G_Bs","Reconstructed Bs mass",*B0m,*B0mass_Bs,*B0mass_sigma_G);

  // *** Exponential Background ***
  B0mass_Exp_bkg = new RooExponential("B0mass_Exp_bkg","Reconstructed B mass for background",*B0m,*B0mass_lambda); // Inclusive background samples


////////////////////////////////////////////////////
// ********** Building Bs Lifetime PDF ********** //
////////////////////////////////////////////////////

  BsLifetime_Res_Decay = new RooGaussModel("BsLifetime_Res_Decay","Reconstructed Bs lifetime for Decay",*Bst,*resMean,*resDecay);
  BsLifetime_Gauss_Prompt_III = new RooGaussian("BsLifetime_Gauss_Prompt_I","Reconstructed Bs lifetime for Combinatorics",*Bst,*resMean,*resPrompt_III);

  decay_Bs = new RooBDecay("decay_Bs","Decay Bs->J/#psi Ks",*Bst, *tau_Bs,*dG_Bs,*ACosh_Bs,*DSinh_Bs,*CCos_Bs,*SSin_Bs,*dM_Bs,*BsLifetime_Res_Decay,RooBDecay::SingleSided);// B signal



//////////////////////////////////////////////
// ********** Building Total PDF ********** //
//////////////////////////////////////////////

  Sig_Bs = new RooProdPdf("Sig_Bs","Total PDF - Signal Bs - LL",RooArgList(*B0mass_G_Bs,*decay_Bs));

  Jpsi = new RooProdPdf("Jpsi","Total PDF - Jpsi - LL",RooArgList(*B0mass_Exp_bkg,*BsLifetime_Gauss_Prompt_III));

  total_pdf = new RooAddPdf("total_pdf","Total PDF - Long Long",RooArgList(*Sig_Bs,*Jpsi),RooArgList(*nevents_fit_Bs,*nevents_fit_Jpsi),kFALSE);



///////////////////////////////////////////
// ********** Generating data ********** //
///////////////////////////////////////////

  std::cout << "*** Starting to generate data from PDF ***" << std::endl;

  // *** Identifying and starting Random Generator ***
  RooRandom::randomGenerator()->GetName();
  RooRandom::randomGenerator()->GetSeed();
  RooRandom::randomGenerator()->SetSeed(123);
  RooRandom::randomGenerator()->GetSeed();

  std::cout << "*** Using Random Generator:  " << RooRandom::randomGenerator()->GetName() << " ***" << std::endl;
  std::cout << "*** Starting with given seed = " << RooRandom::randomGenerator()->GetSeed() << " ***" << std::endl;

  // ********** Generating ********** //
  //////////////////////////////////////

  std::cout << "*** Generating Proto Data for RooCategory ***" << std::endl;
  dataset_proto = new RooDataSet("dataset_proto","Proto Data for RooCatergory",*tag);

  for(Int_t i=0; i<tageff; i++){
    tag->setLabel("Bs");
    dataset_proto->add(*tag);
  }
  for(Int_t i=0; i<tageff; i++){
    tag->setLabel("Bsbar");
    dataset_proto->add(*tag);
  }
  for(Int_t i=0; i<2*(100-tageff); i++){
    tag->setLabel("untag");
    dataset_proto->add(*tag);
  }
  std::cout << "Relative distribution of tagging categories" << std::endl;
  Roo1DTable* protoTable = dataset_proto->table(*tag) ;
  protoTable->Print("v") ;

  std::cout << "*** Generating events ***" << std::endl;
  //dataset = total_pdf->generate(RooArgSet(*B0m,*Bst,*tag), nevents_toy_input, ProtoData(*dataset_proto,kTRUE));// Original Code - Works with Root 5.24 on Stoomboot
  dataset = total_pdf->generate(RooArgSet(*B0m,*Bst), nevents_toy_input, ProtoData(*dataset_proto,kTRUE));

  //dataset = total_pdf->generate(RooArgSet(*B0m,*Bst,*tag), nevents_toy_input); // No ProtoData -> each tagging category gets 33.3%

// Note:
// Something seems wrong with the use of ProtoData in Root 5.28
// The data that is associated with the RooBDecay, which is the only one that explicitly depends on the RooCategory gets messed up


////////////////////////////////////
// ********** Plotting ********** //
////////////////////////////////////

  // ********** B0 Mass frame ********** //
  /////////////////////////////////////////

  std::cout << "*** Plotting PDF components: B0 mass ***" << std::endl;

  // *** Datapoints and distribution ***
  RooPlot *B0mframe = B0m->frame(Range(5349,5529),Bins(36)); // Full
  B0mframe->SetYTitle("# Events/(5 MeV/c^{2})");

  B0mframe->SetXTitle("B Mass (MeV/c^{2})");
  B0mframe->SetTitle("Decay of B_{d/s}#rightarrow J/#psi K_{S}");
  B0mframe->SetTitleOffset(1.3,"y");
  B0mframe->SetTitleSize(0.07,"x");
  B0mframe->SetTitleSize(0.07,"y");
  B0mframe->SetTitleFont(132,"y");
  B0mframe->SetTitleFont(132,"x");
  B0mframe->SetLabelSize(0.05,"x");
  B0mframe->SetLabelSize(0.05,"y");
  B0mframe->SetLabelFont(132,"x");
  B0mframe->SetLabelFont(132,"y");
  B0mframe->GetXaxis()->SetNdivisions(505,kTRUE);

  // *** All tracks ***
  dataset->plotOn(B0mframe,LineColor(1),MarkerStyle(8),MarkerSize(0.8));
  total_pdf->plotOn(B0mframe,LineColor(2),LineWidth(2),LineStyle(11));
  total_pdf->plotOn(B0mframe,Components(RooArgSet(*Jpsi)),LineColor(3),LineWidth(2));

  // ********** Bs Lifetime frame ********** //
  /////////////////////////////////////////////

  std::cout << "*** Plotting PDF components: Bs Lifetime ***" << std::endl;

  // *** Datapoints and distribution ***
  RooPlot *Bstframe = Bst->frame(Range(-0.5,8),Bins(85)); 
  Bstframe->SetYTitle("# Events/(0.1 ps)");
  Bstframe->SetXTitle("B_{s}^{0} Lifetime (ps)");
  Bstframe->SetTitle("Decay of B_{s}#rightarrow J/#psi K_{S}");

  Bstframe->SetTitleOffset(1.3,"y");
  Bstframe->SetTitleSize(0.07,"x");
  Bstframe->SetTitleSize(0.07,"y");
  Bstframe->SetTitleFont(132,"y");
  Bstframe->SetTitleFont(132,"x");
  Bstframe->SetLabelSize(0.05,"x");
  Bstframe->SetLabelSize(0.05,"y");
  Bstframe->SetLabelFont(132,"x");
  Bstframe->SetLabelFont(132,"y");

  // *** Bs only ***
  dataset->plotOn(Bstframe,Cut("tag==1"),MarkerStyle(8),MarkerSize(0.8));
  total_pdf->plotOn(Bstframe,Slice(*tag,"Bs"),Components(RooArgSet(*Sig_Bs)),LineColor(4),LineWidth(2));
  total_pdf->plotOn(Bstframe,Slice(*tag,"Bs"),Components(RooArgSet(*Jpsi)),LineColor(2),LineWidth(2),LineStyle(11));

  Bstframe->SetMinimum(1);
  Bstframe->SetMaximum(100000);


  // ********** Draw ********** //
  ////////////////////////////////

  std::cout << "*** Plotting PDF components on Canvas ***" << std::endl;

  TCanvas *B0doek = new TCanvas("B0doek","B0doek",900,800);
  B0doek->SetLeftMargin(0.2);
  B0doek->SetBottomMargin(0.15);
  B0doek->cd(1);
  B0mframe->Draw();

  TCanvas *Bst_doek = new TCanvas("Bst_doek","Bst_doek",900,800);
  Bst_doek->SetLeftMargin(0.2);
  Bst_doek->SetBottomMargin(0.15);
  Bst_doek->cd(1);
  gPad->SetLogy();
  Bstframe->Draw();

}