void bug4(){
  using namespace RooFit;

  RooRealVar *m  = new RooRealVar( "m","invariant mass" , 6100., 6700., "MeV/c^{2}");
  RooRealVar *sm = new RooRealVar("sm","mass resolution",    0.,  100., "MeV/c^{2}");
  RooRealVar *t  = new RooRealVar( "t","proper time",        0.,   5.0, "ps"       );
  RooRealVar *st = new RooRealVar("st","time resolution",    0.001,  0.08, "ps"       );
  st->setBins(100) ;

  RooDataSet* data = RooDataSet::read("asciiData.file", RooArgList(*m, *sm, *t, *st) );

  RooPlot* tframe = t->frame();
  data->plotOn(tframe);
  tframe->Draw();

  // the set of observables...
  RooArgSet* observables = new RooArgSet(*m,*sm,*t,*st);

  //Define the PDF  
  //--------------------------------------------------------------------------------
  //Mass
  RooRealVar *SigMass = new RooRealVar("SigMass","Signal mass", 6400, 6100, 6700, "MeV/c^{2}" );
  RooRealVar *MassResScale = new RooRealVar("MassResScale", "mass resolution scale",1.0, 0.1, 3.0 );
  RooFormulaVar *SigMassRes = new RooFormulaVar("SigMassRes", "Sig mass resolution",
						   "(@0*@1)", RooArgList(*MassResScale,*sm) );

  RooAbsPdf  *SigMassPdf = new RooGaussian("SigMassPdf","MassPDF of signal",*m,*SigMass,*SigMassRes);

  //
  //...propertime distributions
  //=================================================================================
  //Scale factor of Resolutions
  RooRealVar *TauResScale = new RooRealVar("TauResScale","time resolution scale",0.9936, 0.1, 3);  

  //
  //Singal
  //------------------------------------------------
  RooRealVar *SigTau = new RooRealVar("SigTau","Signal lifetime",0.10,0.001,10.0,"ps");
  RooRealVar *SigTauResBias = new RooRealVar("SigTauResBias","bias scale",0,-5,5);
  RooResolutionModel *SigTauRes = new RooGaussModel("SigTauRes","gauss resolution model",*t,
                                                   *SigTauResBias,*TauResScale,*st,*st);
  RooAbsPdf  *SigTauPdf = new RooDecay("SigTauPdf","unbiased proper time PDF",*t,
                                      *SigTau,*SigTauRes,RooDecay::SingleSided);

  
  // Take into account acceptance
  // signal acceptance model
  RooRealVar *Acc_a  = new RooRealVar("Acc_a","acceptance a", 3.189  );
  RooRealVar *Acc_n  = new RooRealVar("Acc_n","acceptance n", 3.14   );
  RooRealVar *Acc_c  = new RooRealVar("Acc_c","acceptance c", 0.1515 );
 
  RooAbsReal *Acc = new RooFormulaVar("Acc","(@0>0)*@3*(@1*@0)**@2/( 1+(@1*@0)**@2 )", RooArgList(*t,*Acc_a,*Acc_n,*Acc_c) );

  RooAbsPdf *SigTauAccPdf = new RooEffProd("SigTauAccPdf","biased proper time PDF", *SigTauPdf, *Acc);

  // product: unconditional
  RooAbsPdf *SigPdf = new RooProdPdf("SigPdf" ,"signal PDF", RooArgList(*SigMassPdf, *SigTauAccPdf) );

  SigTauResBias->setConstant(kTRUE);
  TauResScale->setConstant(kTRUE);

  // fitting


  RooArgSet* param = SigPdf->getParameters(*observables) ;
  if (0) {
    RooFitResult *s = SigPdf->fitTo(*data,Minos(false),Strategy(1),Save(true),Verbose(true),ConditionalObservables(RooArgSet(*sm,*st)),NumCPU(8));
    s->Print("v");
    param->writeToFile("bug4_params.txt") ;
  } else {
    param->readFromFile("bug4_params.txt") ;   
  }


  // Generate new (m,t) data from p.d.f taking (sm,st) from original dataset
  RooDataSet* protodata = data->reduce(RooArgSet(*sm,*st)) ;
  data2 = SigPdf->generate(RooArgSet(*m,*t),data->numEntries(),ProtoData(*protodata)) ;  

  SigPdf->fitTo(*data,Minos(false),Strategy(1),Save(true),Verbose(true),ConditionalObservables(RooArgSet(*sm,*st)),NumCPU(8));
  
  // make plots
  TCanvas *c = new TCanvas();
  c->Divide(1,3);
 
  RooPlot *p1 = m->frame(Bins(100),Title("invariant mass"));
  data->plotOn(p1);
  data2->plotOn(p1,MarkerColor(kRed),LineColor(kRed));

  SigPdf->plotOn(p1,ProjWData(*sm,*data,kTRUE));
  c->cd(1); 
  p1->Draw();
  
  RooDataHist bdata_st("bdata_st","bdata_st",*st,*data) ;
  RooPlot *p2 = t->frame(Bins(100),Title("proper time"));
  data->plotOn( p2 );
  data2->plotOn(p2,MarkerColor(kRed),LineColor(kRed));
  TStopwatch tw ;
  tw.Start() ;
  SigPdf->plotOn( p2, ProjWData(*st,*data,kTRUE),NumCPU(2,kTRUE)); 
  tw.Stop() ;
  tw.Print() ;
 
  c->cd(2); 
  p2->Draw();


  RooPlot *p3 = t->frame(Bins(100),Title("proper time"));
  Acc->plotOn(p3) ;
  c->cd(3) ;
  p3->Draw() ;

}