C++ RooDataSet::GetName方法代码示例

void MakePlots(RooWorkspace* ws){

  // Here we make plots of the discriminating variable (invMass) after the fit
  // and of the control variable (isolation) after unfolding with sPlot.
  std::cout << "make plots" << std::endl;

  // make our canvas
  TCanvas* cdata = new TCanvas("sPlot","sPlot demo", 400, 600);
  cdata->Divide(1,3);

  // get what we need out of the workspace
  RooAbsPdf* model = ws->pdf("model");
  RooAbsPdf* zModel = ws->pdf("zModel");
  RooAbsPdf* qcdModel = ws->pdf("qcdModel");

  RooRealVar* isolation = ws->var("isolation");
  RooRealVar* invMass = ws->var("invMass");

  // note, we get the dataset with sWeights
  RooDataSet* data = (RooDataSet*) ws->data("dataWithSWeights");

  // this shouldn't be necessary, need to fix something with workspace
  // do this to set parameters back to their fitted values.
  model->fitTo(*data, Extended() );

  //plot invMass for data with full model and individual componenets overlayed
  //  TCanvas* cdata = new TCanvas();
  cdata->cd(1);
  RooPlot* frame = invMass->frame() ; 
  data->plotOn(frame ) ; 
  model->plotOn(frame) ;   
  model->plotOn(frame,Components(*zModel),LineStyle(kDashed), LineColor(kRed)) ;   
  model->plotOn(frame,Components(*qcdModel),LineStyle(kDashed),LineColor(kGreen)) ;   
    
  frame->SetTitle("Fit of model to discriminating variable");
  frame->Draw() ;
 

  // Now use the sWeights to show isolation distribution for Z and QCD.  
  // The SPlot class can make this easier, but here we demonstrait in more
  // detail how the sWeights are used.  The SPlot class should make this 
  // very easy and needs some more development.

  // Plot isolation for Z component.  
  // Do this by plotting all events weighted by the sWeight for the Z component.
  // The SPlot class adds a new variable that has the name of the corresponding
  // yield + "_sw".
  cdata->cd(2);

  // create weightfed data set 
  RooDataSet * dataw_z = new RooDataSet(data->GetName(),data->GetTitle(),data,*data->get(),0,"zYield_sw") ;

  RooPlot* frame2 = isolation->frame() ; 
  dataw_z->plotOn(frame2, DataError(RooAbsData::SumW2) ) ; 
    
  frame2->SetTitle("isolation distribution for Z");
  frame2->Draw() ;

  // Plot isolation for QCD component.  
  // Eg. plot all events weighted by the sWeight for the QCD component.
  // The SPlot class adds a new variable that has the name of the corresponding
  // yield + "_sw".
  cdata->cd(3);
  RooDataSet * dataw_qcd = new RooDataSet(data->GetName(),data->GetTitle(),data,*data->get(),0,"qcdYield_sw") ;
  RooPlot* frame3 = isolation->frame() ; 
  dataw_qcd->plotOn(frame3,DataError(RooAbsData::SumW2) ) ; 
    
  frame3->SetTitle("isolation distribution for QCD");
  frame3->Draw() ;

  //  cdata->SaveAs("SPlot.gif");

}

void rf403_weightedevts()
{
  // C r e a t e   o b s e r v a b l e   a n d   u n w e i g h t e d   d a t a s e t 
  // -------------------------------------------------------------------------------

  // Declare observable
  RooRealVar x("x","x",-10,10) ;
  x.setBins(40) ;

  // Construction a uniform pdf
  RooPolynomial p0("px","px",x) ;

  // Sample 1000 events from pdf
  RooDataSet* data = p0.generate(x,1000) ;

 

  // C a l c u l a t e   w e i g h t   a n d   m a k e   d a t a s e t   w e i g h t e d 
  // -----------------------------------------------------------------------------------

  // Construct formula to calculate (fake) weight for events
  RooFormulaVar wFunc("w","event weight","(x*x+10)",x) ;

  // Add column with variable w to previously generated dataset
  RooRealVar* w = (RooRealVar*) data->addColumn(wFunc) ;

  // Dataset d is now a dataset with two observable (x,w) with 1000 entries
  data->Print() ;

  // Instruct dataset wdata in interpret w as event weight rather than as observable
  RooDataSet wdata(data->GetName(),data->GetTitle(),data,*data->get(),0,w->GetName()) ;

  // Dataset d is now a dataset with one observable (x) with 1000 entries and a sum of weights of ~430K
  wdata.Print() ;



  // U n b i n n e d   M L   f i t   t o   w e i g h t e d   d a t a 
  // ---------------------------------------------------------------

  // Construction quadratic polynomial pdf for fitting
  RooRealVar a0("a0","a0",1) ;
  RooRealVar a1("a1","a1",0,-1,1) ;
  RooRealVar a2("a2","a2",1,0,10) ;
  RooPolynomial p2("p2","p2",x,RooArgList(a0,a1,a2),0) ;

  // Fit quadratic polynomial to weighted data

  // NOTE: A plain Maximum likelihood fit to weighted data does in general 
  //       NOT result in correct error estimates, unless individual
  //       event weights represent Poisson statistics themselves.
  //       
  // Fit with 'wrong' errors
  RooFitResult* r_ml_wgt = p2.fitTo(wdata,Save()) ;
  
  // A first order correction to estimated parameter errors in an 
  // (unbinned) ML fit can be obtained by calculating the
  // covariance matrix as
  //
  //    V' = V C-1 V
  //
  // where V is the covariance matrix calculated from a fit
  // to -logL = - sum [ w_i log f(x_i) ] and C is the covariance
  // matrix calculated from -logL' = -sum [ w_i^2 log f(x_i) ] 
  // (i.e. the weights are applied squared)
  //
  // A fit in this mode can be performed as follows:

  RooFitResult* r_ml_wgt_corr = p2.fitTo(wdata,Save(),SumW2Error(kTRUE)) ;



  // P l o t   w e i g h e d   d a t a   a n d   f i t   r e s u l t 
  // ---------------------------------------------------------------

  // Construct plot frame
  RooPlot* frame = x.frame(Title("Unbinned ML fit, binned chi^2 fit to weighted data")) ;

  // Plot data using sum-of-weights-squared error rather than Poisson errors
  wdata.plotOn(frame,DataError(RooAbsData::SumW2)) ;

  // Overlay result of 2nd order polynomial fit to weighted data
  p2.plotOn(frame) ;



  // M L  F i t   o f   p d f   t o   e q u i v a l e n t  u n w e i g h t e d   d a t a s e t
  // -----------------------------------------------------------------------------------------
  
  // Construct a pdf with the same shape as p0 after weighting
  RooGenericPdf genPdf("genPdf","x*x+10",x) ;

  // Sample a dataset with the same number of events as data
  RooDataSet* data2 = genPdf.generate(x,1000) ;

  // Sample a dataset with the same number of weights as data
  RooDataSet* data3 = genPdf.generate(x,43000) ;

  // Fit the 2nd order polynomial to both unweighted datasets and save the results for comparison
  RooFitResult* r_ml_unw10 = p2.fitTo(*data2,Save()) ;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  argSet.add(isoFPRW2_EB);


  //binning variables
  //std::cout<<"Binning variables for eventual binned fit"<<std::endl;
  //combinedPfIsoFPR03Phot.setBins(121);
  //etaPhot.setBins(120);
  //mvaIdPhot.setBins(180);
  //isIsolatedGenPhot.setBins(3);
  //ptPhot.setBins(1200);
  //weight_times_isoWeight.setBins(10000000);

  //creating complete dataset
  std::cout<<"Reading trees of MC into a complete general dataset"<<std::endl;
  RooDataSet* allSet;
  if(isFPR){
    if(isoWeight==0)
      allSet = new RooDataSet("allSet", "allSet", argSet, WeightVar("isoFPRW_EB"), RooFit::Import(data));
    if(isoWeight==1)
      allSet = new RooDataSet("allSet", "allSet", argSet, WeightVar("isoFPRW1_EB"), RooFit::Import(data));
    if(isoWeight==2)
      allSet = new RooDataSet("allSet", "allSet", argSet, WeightVar("isoFPRW2_EB"), RooFit::Import(data));
  }
  else{
    if(isoWeight==0)
      allSet = new RooDataSet("allSet", "allSet", argSet, WeightVar("isoW_EB"), RooFit::Import(data));
    if(isoWeight==1)
      allSet = new RooDataSet("allSet", "allSet", argSet, WeightVar("isoW1_EB"), RooFit::Import(data));
    if(isoWeight==2)
      allSet = new RooDataSet("allSet", "allSet", argSet, WeightVar("isoW2_EB"), RooFit::Import(data));
  }
  //allSet.setWeightVar(weight_times_isoWeight);

  std::cout<<"Complete dataset "<<allSet->GetName()<<" created"<<std::endl<<std::endl;
  std::cout<<"******************** "<<std::endl;
  std::cout<<" isWeighted =  " <<  allSet->isWeighted() <<std::endl;
  std::cout<<"******************** "<<std::endl<<std::endl;

  //reducing complete dataset to interesting ones
  
  RooDataSet* d_r = (RooDataSet*)allSet->reduce((cut+" && mvaIdPhot<0.6 && mvaIdPhot>-0.6").c_str());
  d_r->SetName("d_r");
  std::cout<<"Reduced dataset "<<d_r->GetName()<<" created with cut "<<cut<<" && -0.6 < mvaIdPhot < 0.6"<<std::endl;
    std::cout<<"d_r entries: "<<d_r->sumEntries()<<std::endl;
  
  //models for fit
  //gaussian
  RooRealVar gaussmean("gaussmean","gaussmean", -1., -10., 20.);
  RooRealVar gausssigma("gausssigma", "gausssigma", 1., 0., 20.);

  RooGaussian* my_gauss;
  if(isFPR)
    my_gauss = new RooGaussian("my_gauss", "my_gauss", combinedPfIsoFPR03Phot, gaussmean, gausssigma);
  else
    my_gauss = new RooGaussian("my_gauss", "my_gauss", combinedPfIso03Phot, gaussmean, gausssigma);
  

  //crystalBall
  RooRealVar cbmean("cbmean", "cbmean", -1., -5., 15.);
  RooRealVar cbsigma("cbsigma", "cbsigma", 1., 0., 20.);

  RooRealVar cbalpha("cbalpha", "cbaplha", -1., -1000., 0.);
  RooRealVar cbn("cbn","cbn",10., 0., 1000.);

  RooCBShape* my_cb;
  if(isFPR)

//.........这里部分代码省略.........
    wsp->Print();

    RooRealVar* B_postcalib_M = wsp -> var("B_postcalib_M");
    RooRealVar* nsig_sw = wsp-> var("nsig_sw");
    RooRealVar* nbkg_sw = wsp-> var("nbkg_sw");
    RooRealVar* B_M13_Subst3_gamma2pi0 = wsp-> var("B_M13_Subst3_gamma2pi0");
    RooRealVar* B_M023 = wsp -> var("B_M023");
    RooRealVar* K_1_1270_plus_M = wsp -> var("K_1_1270_plus_M");
    RooRealVar* K_1_1270_plus_SMALLESTDELTACHI2 = wsp -> var("K_1_1270_plus_SMALLESTDELTACHI2");
    RooRealVar* gamma_CL = wsp -> var("gamma_CL");
    RooRealVar* piminus_PIDK = wsp -> var("piminus_PIDK");
    RooRealVar* piplus_PIDK = wsp -> var("piplus_PIDK");
    RooRealVar* Kplus_PIDp = wsp -> var("Kplus_PIDp");
    RooRealVar* Kplus_PIDK = wsp -> var("Kplus_PIDK");
    RooRealVar* B_M02 = wsp -> var("B_M02");
    RooRealVar* L_nsig = wsp -> var("L_nsig");
    RooRealVar* L_nbkg = wsp -> var("L_nbkg");

    RooArgSet arg(*B_postcalib_M,*gamma_CL,*B_M13_Subst3_gamma2pi0,*B_M023,*piminus_PIDK,*piplus_PIDK,*Kplus_PIDK,*Kplus_PIDp);
    arg.add(*K_1_1270_plus_M);
    arg.add(*K_1_1270_plus_SMALLESTDELTACHI2);
    arg.add(*B_M02);
    arg.add(*nsig_sw);
    arg.add(*L_nsig);
    arg.add(*nbkg_sw);
    arg.add(*L_nbkg);

    arg.add(*m_Kpipi);


    RooDataSet* DataSWeights = (RooDataSet*) wsp -> data("DataSWeights");
    RooFormulaVar newMass("m_Kpipi", "m_Kpipi", "K_1_1270_plus_M", RooArgList(*(wsp->var("K_1_1270_plus_M"))));
    DataSWeights->addColumn(newMass);
    RooDataSet* splot = new RooDataSet(DataSWeights->GetName(),DataSWeights->GetTitle(),DataSWeights,RooArgSet(arg),"","nsig_sw");


    cout << "\n\n >>>> Defining components and fitting  \n\n" << endl;

    // Defining here pdfs for other resonances to be fitted

    // K1(1270)

    Double_t R = 0.0015; //  was 3.1 GeV-1
    RooRealVar mean_K1_1270("mean_K1_1270","",1272.,1262.,1282.);
    RooRealVar width_K1_1270("width_K1_1270","",90.,70.,110.);

    // K1(1270) -> K rho
    /*RooBreitWigner totalPdf_K1_1270("totalPdf_K1_1270","",*m_Kpipi,mean_K1_1270,width_K1_1270);*/
    /*RooRelBreitWigner totalPdf_K1_1270("totalPdf_K1_1270","totalPdf_K1_1270_rho",*m_Kpipi,mean_K1_1270,width_K1_1270,RooConst(0),RooConst(R),RooConst(770.),RooConst(493.7));*/

    // K1(1270) -> K*0(1430) pi
    /*RooBreitWigner totalPdf_K1_1270_Kst0_1430("totalPdf_K1_1270_Kst0_1430","totalPdf_K1_1270_Kst0_1430",*m_Kpipi,mean_K1_1270,width_K1_1270);*/
    /*RooRelBreitWigner totalPdf_K1_1270_Kst0_1430("totalPdf_K1_1270_Kst0_1430","totalPdf_K1_1270_Kst0_1430",*m_Kpipi,mean_K1_1270,width_K1_1270,RooConst(0),RooConst(R),RooConst(1425.),RooConst(139.6));*/

    // K1(1270) -> K*0(892) pi
    /*RooBreitWigner totalPdf_K1_1270_Kst0_892("totalPdf_K1_1270_Kst0_892","totalPdf_K1_1270_Kst0_892",*m_Kpipi,mean_K1_1270,width_K1_1270);*/
    /*RooRelBreitWigner totalPdf_K1_1270_Kst0_892("totalPdf_K1_1270_Kst0_892","totalPdf_K1_1270_Kst0_892",*m_Kpipi,mean_K1_1270,width_K1_1270,RooConst(0),RooConst(R),RooConst(895.5),RooConst(139.6));*/

    // K1(1400)

    RooRealVar mean_K1_1400("mean_K1_1400","",1403./*,1396.,1410.*/);
    RooRealVar width_K1_1400("width_K1_1400","",174./*,151.,197.*/);
    RooBreitWigner totalPdf_K1_1400("totalPdf_K1_1400","",*m_Kpipi,mean_K1_1400,width_K1_1400);
    /*[>RooRelBreitWigner totalPdf_K1_1400("totalPdf_K1_1400","totalPdf_K1_1400",*m_Kpipi,mean_K1_1400,width_K1_1400,RooConst(0),RooConst(R),RooConst(895.),RooConst(139.6)); //K*(892) pi is 93%<]*/

    // K2*(1430)

void ws_v03()
{
    gROOT->ProcessLine(".x ./mystyle.C");

    TFile *f1 = new TFile("K1_1270/ws_K1_1270.root");
    TFile *f2 = new TFile("K1_1400/ws_K1_1400.root");
    TFile *f3 = new TFile("K2_1430/ws_K2_1430.root");

    RooWorkspace* ws_K1_1270 = (RooWorkspace*) f1->Get("ws_K1_1270");
    RooWorkspace* ws_K1_1400 = (RooWorkspace*) f2->Get("ws_K1_1400");
    RooWorkspace* ws_K2_1430 = (RooWorkspace*) f3->Get("ws_K2_1430");
    ws_K1_1270->Print();
    ws_K1_1400->Print();
    ws_K2_1430->Print();

    // Importing variables from workspaces

    RooRealVar* m_Kpipi = ws_K1_1270 -> var("m_Kpipi");
    RooAbsPdf* totalPdf_K1_1270 = ws_K1_1270 -> pdf("totalPdf_K1_1270");
    RooAbsData* data_K1_1270 = ws_K1_1270 -> data("totalPdf_K1_1270Data");
    RooHistPdf* pdf_K1_1270_to_Krho = ws_K1_1270 -> pdf("histPdf_K1toKrho");

    /*RooRealVar* m_Kpipi = ws_K1_1400 -> var("m_Kpipi");*/
    /*RooAbsPdf* totalPdf_K1_1400 = ws_K1_1400 -> pdf("totalPdf_K1_1400");*/
    /*RooAbsData* data_K1_1400 = ws_K1_1400 -> data("totalPdf_K1_1400Data");*/
    /*RooHistPdf* pdf_K1_1400_to_Krho = ws_K1_1400 -> pdf("histPdf_K1_1400toKrho");*/

    /*RooRealVar* m_Kpipi = ws_K2_1430 -> var("m_Kpipi");*/
    /*RooAbsPdf* totalPdf_K2_1430 = ws_K2_1430 -> pdf("totalPdf_K2_1430");*/
    /*RooAbsData* data_K2_1430 = ws_K2_1430 -> data("totalPdf_K2_1430Data");*/
    /*RooHistPdf* pdf_K2_1430_to_Krho = ws_K2_1430 -> pdf("histPdf_K2_1430toKrho");*/

    // Plotting pdf

    /*TCanvas* c1 = new TCanvas("c1","canvas",20,20,1200,600);*/
    /*c1->Divide(3,1);*/
    /*c1->cd(1);*/
    /*RooPlot* frame_K1_1270 = m_Kpipi -> frame(Bins(200),Title("K1(1270) -> K#pi#pi"));*/
    /*data_K1_1270->plotOn(frame_K1_1270,DrawOption("C"));*/
    /*frame_K1_1270->Draw();*/
    /*c1->cd(2);*/
    /*RooPlot* frame_K1_1400 = m_Kpipi -> frame(Bins(200),Title("K1(1400) -> K#pi#pi"));*/
    /*data_K1_1400->plotOn(frame_K1_1400,DrawOption("C"));*/
    /*frame_K1_1400->Draw();*/
    /*c1->cd(3);*/
    /*RooPlot* frame_K2_1430= m_Kpipi -> frame(Bins(200),Title("K2*(1430) -> K#pi#pi"));*/
    /*data_K2_1430->plotOn(frame_K2_1430,DrawOption("C"));*/
    /*frame_K2_1430->Draw();*/

 ////////////////////////////////////////////////////////////////////////

    TFile *Fworkspace = new TFile("workspace.root");

    RooWorkspace* wsp = (RooWorkspace*) Fworkspace->Get("wsp");
    wsp->Print();

    RooRealVar* B_postcalib_M = wsp -> var("B_postcalib_M");
    RooRealVar* nsig_sw = wsp-> var("nsig_sw");
    RooRealVar* nbkg_sw = wsp-> var("nbkg_sw");
    RooRealVar* B_M13_Subst3_gamma2pi0 = wsp-> var("B_M13_Subst3_gamma2pi0");
    RooRealVar* B_M023 = wsp -> var("B_M023");
    RooRealVar* K_1_1270_plus_M = wsp -> var("K_1_1270_plus_M");
    RooRealVar* K_1_1270_plus_SMALLESTDELTACHI2 = wsp -> var("K_1_1270_plus_SMALLESTDELTACHI2");
    RooRealVar* gamma_CL = wsp -> var("gamma_CL");
    RooRealVar* piminus_PIDK = wsp -> var("piminus_PIDK");
    RooRealVar* piplus_PIDK = wsp -> var("piplus_PIDK");
    RooRealVar* Kplus_PIDp = wsp -> var("Kplus_PIDp");
    RooRealVar* Kplus_PIDK = wsp -> var("Kplus_PIDK");
    RooRealVar* B_M02 = wsp -> var("B_M02");
    RooRealVar* L_nsig = wsp -> var("L_nsig");
    RooRealVar* L_nbkg = wsp -> var("L_nbkg");

    RooArgSet arg(*B_postcalib_M,*gamma_CL,*B_M13_Subst3_gamma2pi0,*B_M023,*piminus_PIDK,*piplus_PIDK,*Kplus_PIDK,*Kplus_PIDp);
    arg.add(*K_1_1270_plus_M);
    arg.add(*K_1_1270_plus_SMALLESTDELTACHI2);
    arg.add(*B_M02);
    arg.add(*nsig_sw);
    arg.add(*L_nsig);
    arg.add(*nbkg_sw);
    arg.add(*L_nbkg);

    arg.add(*m_Kpipi);


    RooDataSet* DataSWeights = (RooDataSet*) wsp -> data("DataSWeights");
    RooFormulaVar newMass("m_Kpipi", "m_Kpipi", "K_1_1270_plus_M", RooArgList(*(wsp->var("K_1_1270_plus_M"))));
    DataSWeights->addColumn(newMass);
    RooDataSet* splot = new RooDataSet(DataSWeights->GetName(),DataSWeights->GetTitle(),DataSWeights,RooArgSet(arg),"","nsig_sw");


    // Defining here pdfs for other resonances to be fitted

// TRY TO ADD 2 SIMPLE BWs FOR K1 1270 AND K1 1400 AS WELL
// also, binned clone and try to fit a toy dataset

        // K1(1270)

    /*Double_t R = 0.0015; //  was 3.1 GeV-1*/
    /*RooRealVar mean_K1_1270("mean_K1_1270","",1272.,1262.,1282.);*/
    /*RooRealVar width_K1_1270("width_K1_1270","",90.,70.,110.);*/
//.........这里部分代码省略.........

void MakePlots(TString filename, TString energy="8TeV", TString lumi=""){
  TString outDir=filename; outDir.ReplaceAll("fitres","img");
  outDir="tmp/";
  //std::map<TString, TH2F *> deltaNLL_map;
  
  /*------------------------------ Plotto */
  TCanvas *c = new TCanvas("c","c");
  
  TFile f_in(filename, "read");
  if(f_in.IsZombie()){
    std::cerr << "File opening error: " << filename << std::endl;
    return;
  }
  
  TList *KeyList = f_in.GetListOfKeys();
  std::cout << KeyList->GetEntries() << std::endl;
  for(int i =0; i <  KeyList->GetEntries(); i++){
    c->Clear();
    TKey *key = (TKey *)KeyList->At(i);
    if(TString(key->GetClassName())!="RooDataSet") continue;
    RooDataSet *dataset = (RooDataSet *) key->ReadObj();
    
    TString constTermName = dataset->GetName();
    TString alphaName=constTermName; alphaName.ReplaceAll("constTerm","alpha");

    TTree *tree = dataset2tree(dataset);
    TGraphErrors bestFit_ = bestFit(tree, alphaName, constTermName);
    TH2F *hist = prof2d(tree, alphaName, constTermName, "nll", "(12,-0.0005,0.0115,29,-0.0025,0.1425)",true);

//     //    deltaNLL_map.insert(std::pair <TString, TH2F *>(keyName,hist));
    hist->SaveAs(outDir+"/deltaNLL-"+constTermName+".root");
    hist->Draw("colz");
    bestFit_.Draw("P same");
    bestFit_.SetMarkerSize(2);

    Int_t iBinX, iBinY;
    Double_t x,y;
    hist->GetBinWithContent2(0,iBinX,iBinY);
    x= hist->GetXaxis()->GetBinCenter(iBinX);
    y= hist->GetYaxis()->GetBinCenter(iBinY);
    TGraph nllBestFit(1,&x,&y);

    nllBestFit.SetMarkerStyle(3);
    nllBestFit.SetMarkerColor(kRed);
    TList* contour68 = contourFromTH2(hist, 0.68);

    hist->Draw("colz");
    hist->GetZaxis()->SetRangeUser(0,50);
    bestFit_.Draw("P same");
    nllBestFit.Draw("P same");
    //contour68->Draw("same");
    c->SaveAs(outDir+"/deltaNLL-"+constTermName+".png");
    hist->SaveAs("tmp/hist-"+constTermName+".root");
    nllBestFit.SaveAs("tmp/nllBestFit.root");
    contour68->SaveAs("tmp/contour68.root");
    delete hist;
    hist = prof2d(tree, alphaName, constTermName, "nll", "(12,-0.0005,0.0115,29,-0.0025,0.1425)");
    RooHistPdf *histPdf = nllToL(hist);
    delete hist;
    RooDataSet *gen_dataset=histPdf->generate(*histPdf->getVariables(),1000000,kTRUE,kFALSE);
    TTree *genTree = dataset2tree(gen_dataset);
    genTree->SaveAs("tmp/genTree-"+constTermName+".root");
    delete gen_dataset;
    delete histPdf;
    
    TGraphErrors toyGraph = g(genTree, constTermName);
    TGraphErrors bestFitGraph = g(tree,alphaName, constTermName);
    TGraphErrors bestFitScanGraph = g(y, x);
    delete genTree;
    delete tree;
    toyGraph.SetFillColor(kGreen);
    toyGraph.SetLineColor(kBlue);
    toyGraph.SetLineStyle(2);
    bestFitGraph.SetLineColor(kBlack);
    bestFitScanGraph.SetLineColor(kRed);
    bestFitScanGraph.SetLineWidth(2);


    
    TMultiGraph g_multi("multigraph","");
    g_multi.Add(&toyGraph,"L3");
    g_multi.Add(&toyGraph,"L");
    g_multi.Add(&bestFitGraph, "L");
    g_multi.Add(&bestFitScanGraph, "L");
   
    g_multi.Draw("A");

    c->Clear();
    g_multi.Draw("A");
    c->SaveAs(outDir+"/smearing_vs_energy-"+constTermName+".png");

    //    TPaveText *pv = new TPaveText(0.7,0.7,1, 0.8);    
//     TLegend *legend = new TLegend(0.7,0.8,0.95,0.92);
//     legend->SetFillStyle(3001);
//     legend->SetFillColor(1);
//     legend->SetTextFont(22); // 132
//     legend->SetTextSize(0.04); // l'ho preso mettendo i punti con l'editor e poi ho ricavato il valore con il metodo GetTextSize()
//   //  legend->SetFillColor(0); // colore di riempimento bianco
//     legend->SetMargin(0.4);  // percentuale della larghezza del simbolo
    //    SetLegendStyle(legend);
//.........这里部分代码省略.........

int main(int argc, char *argv[]){

	OptionParser(argc,argv);
	

	RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR);
	RooMsgService::instance().setSilentMode(true);


  
	system(Form("mkdir -p %s",outdir_.c_str()));

	vector<string> procs;
	split(infilenames_,infilenamesStr_,boost::is_any_of(","));
  
   TPython::Exec("import os,imp,re");
    const char * env = gSystem->Getenv("CMSSW_BASE") ; 
      std::string globeRt = env;
      TPython::Exec(Form("buildSMHiggsSignalXSBR = imp.load_source('*', '%s/src/flashggFinalFit/Signal/python/buildSMHiggsSignalXSBR.py')",globeRt.c_str()));
      TPython::Eval(Form("buildSMHiggsSignalXSBR.Init%dTeV()", 13));
   for (unsigned int i =0 ; i<infilenames_.size() ; i++){
	    int mH  =(int) TPython::Eval(Form("int(re.search('_M(.+?)_','%s').group(1))",infilenames_[i].c_str())); 
	   double WH_XS  =  (double)TPython::Eval(Form("buildSMHiggsSignalXSBR.getXS(%d,'%s')",mH,"WH"));
	   double ZH_XS  =  (double)TPython::Eval(Form("buildSMHiggsSignalXSBR.getXS(%d,'%s')",mH,"ZH"));
     float tot_XS = WH_XS + ZH_XS;
     float wFrac=  WH_XS /tot_XS ;
     float zFrac=  ZH_XS /tot_XS ;
      std::cout << "mass "<< mH << " wh fraction "<< WH_XS /tot_XS << ", zh fraction "<< ZH_XS /tot_XS <<std::endl; 
     TFile *infile =  TFile::Open(infilenames_[i].c_str());
	   string outname  =(string) TPython::Eval(Form("'%s'.split(\"/\")[-1].replace(\"VH\",\"WH_VH\")",infilenames_[i].c_str())); 
     TFile *outfile = TFile::Open(outname.c_str(),"RECREATE") ;
    TDirectory* saveDir = outfile->mkdir("tagsDumper");
    saveDir->cd();

    RooWorkspace *inWS = (RooWorkspace*) infile->Get("tagsDumper/cms_hgg_13TeV");
    RooRealVar *intLumi = (RooRealVar*)inWS->var("IntLumi");
    RooWorkspace *outWS = new RooWorkspace("cms_hgg_13TeV");
    outWS->import(*intLumi);
    std::list<RooAbsData*> data =  (inWS->allData()) ;
    std::cout <<" [INFO] Reading WS dataset contents: "<< std::endl;
        for (std::list<RooAbsData*>::const_iterator iterator = data.begin(), end = data.end(); iterator != end; ++iterator )  {
              RooDataSet *dataset = dynamic_cast<RooDataSet *>( *iterator );
              if (dataset) {
	            string zhname  =(string) TPython::Eval(Form("'%s'.replace(\"wzh\",\"zh\")",dataset->GetName())); 
	            string whname  =(string) TPython::Eval(Form("'%s'.replace(\"wzh\",\"wh\")",dataset->GetName())); 
              RooDataSet *datasetZH = (RooDataSet*) dataset->emptyClone(zhname.c_str(),zhname.c_str());
              RooDataSet *datasetWH = (RooDataSet*) dataset->emptyClone(whname.c_str(),whname.c_str());
                TRandom3 r;
                r.Rndm();
                double x[dataset->numEntries()];
                r.RndmArray(dataset->numEntries(),x);
                for (int j =0; j < dataset->numEntries() ; j++){
                    
                    if( x[j] < wFrac){
                    dataset->get(j);
                    datasetWH->add(*(dataset->get(j)),dataset->weight());
                    } else{
                    dataset->get(j);
                    datasetZH->add(*(dataset->get(j)),dataset->weight());
                    }
                }
              float w =datasetWH->sumEntries();
              float z =datasetZH->sumEntries();
           if(verbose_){
              std::cout << "Original dataset " << *dataset <<std::endl;
              std::cout << "WH       dataset " << *datasetWH <<std::endl;
              std::cout << "ZH       dataset " << *datasetZH <<std::endl;
              std::cout << "********************************************" <<std::endl;
              std::cout << "WH fraction (obs) : WH " << w/(w+z) <<",   ZH "<< z/(w+z) << std::endl;
              std::cout << "WH fraction (exp) : WH " << wFrac <<",   ZH "<< zFrac << std::endl;
              std::cout << "********************************************" <<std::endl;
              std::cout << "" <<std::endl;
              std::cout << "" <<std::endl;
              std::cout << "" <<std::endl;
              std::cout << "********************************************" <<std::endl;
              }
               outWS->import(*datasetWH);
               outWS->import(*datasetZH);
                }
             RooDataHist *datahist = dynamic_cast<RooDataHist *>( *iterator );

              if (datahist) {
	            string zhname  =(string) TPython::Eval(Form("'%s'.replace(\"wzh\",\"zh\")",datahist->GetName())); 
	            string whname  =(string) TPython::Eval(Form("'%s'.replace(\"wzh\",\"wh\")",datahist->GetName())); 
              RooDataHist *datahistZH = (RooDataHist*) datahist->emptyClone(zhname.c_str(),zhname.c_str());
              RooDataHist *datahistWH = (RooDataHist*) datahist->emptyClone(whname.c_str(),whname.c_str());
                TRandom3 r;
                r.Rndm();
                double x[datahist->numEntries()];
                r.RndmArray(datahist->numEntries(),x);
                for (int j =0; j < datahist->numEntries() ; j++){
                    
                    datahistWH->add(*(datahist->get(j)),datahist->weight()*wFrac);
                    datahistZH->add(*(datahist->get(j)),datahist->weight()*zFrac);
                }
              float w =datahistWH->sumEntries();
              float z =datahistZH->sumEntries();
           if(verbose_){
              std::cout << "Original datahist " << *datahist <<std::endl;
              std::cout << "WH       datahist " << *datahistWH <<std::endl;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    /*RooArgSet param(*m_Kpipi);*/
    /*RooDataSet MCdataset("MCdataset","MCdataset",param);*/

    /*Double_t K_1_1270_plus_TRUEM = 0.;*/

    /*t_tree->SetBranchAddress("K_1_1270_plus_TRUEM",&K_1_1270_plus_TRUEM);*/

    /*Int_t n = 0;*/
    /*for (int i=0;i<t_tree->GetEntries();i++)*/
    /*{*/
        /*t_tree->GetEntry(i);*/
        /**m_Kpipi_K1_1270 = K_1_1270_plus_TRUEM;*/

        /*MCdataset.add(param);*/
        /*n++;*/
    /*}*/

    /*pdf_K1_1270_to_Krho.fitTo(MCdataset,Extended(true));*/

    /*TCanvas* canvas = new TCanvas("canvas","MC data for K1(1270) -> #rho K",20,20,800,600);*/
    /*RooPlot* frame_MCdataset = m_Kpipi.frame(Bins(200));*/
    /*pdf_K1_1270_to_Krho.paramOn(frame_MCdataset);*/
    /*MCdataset.plotOn(frame_MCdataset,DrawOption("C"));*/
    /*pdf_K1_1270_to_Krho->plotOn(frame_MCdataset);*/
    /*frame_MCdataset->Draw();*/

 ////////////////////////////////////////////////////////////////////////

    TFile *Fworkspace = new TFile("workspace.root");

    RooWorkspace* wsp = (RooWorkspace*) Fworkspace->Get("wsp");
    wsp->Print();

    RooRealVar* B_postcalib_M = wsp -> var("B_postcalib_M");
    RooRealVar* nsig_sw = wsp-> var("nsig_sw");
    RooRealVar* nbkg_sw = wsp-> var("nbkg_sw");
    RooRealVar* B_M13_Subst3_gamma2pi0 = wsp-> var("B_M13_Subst3_gamma2pi0");
    RooRealVar* B_M023 = wsp -> var("B_M023");
    RooRealVar* K_1_1270_plus_M = wsp -> var("K_1_1270_plus_M");
    RooRealVar* K_1_1270_plus_SMALLESTDELTACHI2 = wsp -> var("K_1_1270_plus_SMALLESTDELTACHI2");
    RooRealVar* gamma_CL = wsp -> var("gamma_CL");
    RooRealVar* piminus_PIDK = wsp -> var("piminus_PIDK");
    RooRealVar* piplus_PIDK = wsp -> var("piplus_PIDK");
    RooRealVar* Kplus_PIDp = wsp -> var("Kplus_PIDp");
    RooRealVar* Kplus_PIDK = wsp -> var("Kplus_PIDK");
    RooRealVar* B_M02 = wsp -> var("B_M02");
    RooRealVar* L_nsig = wsp -> var("L_nsig");
    RooRealVar* L_nbkg = wsp -> var("L_nbkg");

    RooArgSet arg(*B_postcalib_M,*gamma_CL,*B_M13_Subst3_gamma2pi0,*B_M023,*piminus_PIDK,*piplus_PIDK,*Kplus_PIDK,*Kplus_PIDp);
    arg.add(*K_1_1270_plus_M);
    arg.add(*K_1_1270_plus_SMALLESTDELTACHI2);
    arg.add(*B_M02);
    arg.add(*nsig_sw);
    arg.add(*L_nsig);
    arg.add(*nbkg_sw);
    arg.add(*L_nbkg);

    arg.add(*m_Kpipi);


    RooDataSet* DataSWeights = (RooDataSet*) wsp -> data("DataSWeights");
    RooFormulaVar newMass("m_Kpipi", "m_Kpipi", "K_1_1270_plus_M", RooArgList(*(wsp->var("K_1_1270_plus_M"))));
    DataSWeights->addColumn(newMass);
    RooDataSet* splot = new RooDataSet(DataSWeights->GetName(),DataSWeights->GetTitle(),DataSWeights,RooArgSet(arg),"","nsig_sw");

    // here add pdfs and FIT sum of pdf to splot

    RooRealVar K1_1270_y("K1_1270_y","K1_1270_y",1000.,0.,10000.);
    RooRealVar K1_1400_y("K1_1400_y","K1_1400_y",100.,0.,10000.);
    /*RooRealVar K2_1430_y("K2_1430_y","K2_1430_y",300.,0.,10000.);*/
    RooFormulaVar K2_1430_y("K2_1430_y","K2_1430_y","K1_1270_y/3.",K1_1270_y); // K2_1430 yield is fixed to 1/3 of the K1_1270 yield (as found from Belle ... TO CHECK)

    RooArgList shapes;
    RooArgList yields;
    shapes.add(*totalPdf_K1_1270);
    shapes.add(*totalPdf_K1_1400);
    shapes.add(*totalPdf_K2_1430);
    yields.add(K1_1270_y);
    yields.add(K1_1400_y);
    yields.add(K2_1430_y);

    RooAddPdf PDF("PDF","total Pdf for the resonances considered", shapes,yields);

    PDF->fitTo(*splot,Extended(),SumW2Error(kTRUE),Range(1000,2000));

    // Plotting

    TCanvas* canvas_sPlot = new TCanvas("canvas_sPlot","sPlot with weights",40,20,800,600);
    RooPlot* frame_splot = m_Kpipi->frame(1000,2000,80);
    splot->plotOn(frame_splot);

    PDF->paramOn(frame_splot);
    PDF->plotOn(frame_splot,Components(*totalPdf_K1_1270),LineColor(kRed),LineStyle(kDashed));
    PDF->plotOn(frame_splot,Components(*totalPdf_K1_1400),LineColor(1),LineStyle(kDashed));
    PDF->plotOn(frame_splot,Components(*totalPdf_K2_1430),LineColor(51),LineStyle(kDashed));
    PDF->plotOn(frame_splot);
    frame_splot->Draw();

}

TTree *ToyTree(TString dirname="test/dato/fitres/Hgg_Et-toys/0.01-0.00", TString fname="outProfile-scaleStep2smearing_7-Et_25-trigger-noPF-EB.root", TString opt="", int nSmooth=10){
  TString outDir=dirname; outDir.ReplaceAll("fitres","img");
  outDir="tmp/";
  //std::map<TString, TH2F *> deltaNLL_map;

  //bool smooth=false;
  //if(opt.Contains("smooth")) smooth=true;
  

  /*------------------------------ Plotto */
  TCanvas c("ctoy","c");
  
  
  TTree *toys = new TTree("toys","");
  toys->SetDirectory(0);
  Double_t constTerm_tree, constTermTrue_tree;
  Double_t alpha_tree, alphaTrue_tree;
  char catName[100]; 
  Int_t catIndex;
  toys->Branch("constTerm", &constTerm_tree, "constTerm/D");
  toys->Branch("alpha", &alpha_tree, "alpha/D");
  toys->Branch("constTermTrue", &constTermTrue_tree, "constTermTrue/D");
  toys->Branch("alphaTrue", &alphaTrue_tree, "alphaTrue/D");

  toys->Branch("catName", catName, "catName/C");
  toys->Branch("catIndex", &catIndex, "catIndex/I");
  std::map<TString, Int_t> catIndexMap;

  ///1/
  for(int itoy =2; itoy <= 50; itoy++){
    TString filename=dirname+"/"; filename+=itoy; filename+="/"+fname;
    TString fout=dirname+"/"; fout+=itoy; fout+="/";    
    TFile f_in(filename, "read");
    if(f_in.IsZombie()){
      std::cerr << "File opening error: " << filename << std::endl;
      continue; //return NULL;
    }
    //std::cout << filename << std::endl;
    TList *KeyList = f_in.GetListOfKeys();
    //std::cout << KeyList->GetEntries() << std::endl;
    for(int i =0; i <  KeyList->GetEntries(); i++){
      c.Clear();
      TKey *key = (TKey *)KeyList->At(i);
      if(TString(key->GetClassName())!="RooDataSet") continue;
      RooDataSet *dataset = (RooDataSet *) key->ReadObj();
    
      TString constTermName = dataset->GetName();
      TString alphaName=constTermName; alphaName.ReplaceAll("constTerm","alpha");
      if(constTermName.Contains("scale")) continue;
      if(constTermName.Contains("alpha")) continue;
      if(constTermName.Contains("1.4442-gold")) continue;
      TTree *tree = dataset2tree(dataset);

      TGraph *rhoGraph = GetRho(tree, alphaName, constTermName);

      rhoGraph->SaveAs(fout+"rhoGraph-"+constTermName+".root");


      TGraphErrors bestFit_ = bestFit(tree, alphaName, constTermName);
      //TString binning="(241,-0.0005,0.2405,61,-0.00025,0.03025)"; //"(40,0.00025,0.02025,61,-0.0022975,0.1401475)";
      TString binning="(241,-0.0005,0.2405,301,-0.00005,0.03005)"; 

      TH2F *hist = prof2d(tree, constTermName, alphaName, "nll", binning, true, nSmooth, opt);
      //hist->SaveAs("myhist.root");
      
      Int_t iBinX, iBinY;
      hist->GetBinWithContent2(0.0002,iBinX,iBinY,1,-1,1,-1,0.0000001);
	
      //      if(iBinX!=0 && iBinY!=0 && iBinX < 41 && iBinY < 62){
      {
	TString catName_=constTermName; catName_.ReplaceAll("constTerm_",""); catName_.ReplaceAll("-","_");
	if(catIndexMap.count(catName_)==0) catIndexMap.insert(std::pair<TString,Int_t>(catName_,catIndexMap.size()));
	catIndex=catIndexMap[catName_];	
	constTerm_tree =  hist->GetYaxis()->GetBinCenter(iBinY);
	alpha_tree = hist->GetXaxis()->GetBinCenter(iBinX);
	sprintf(catName,"%s", catName_.Data());
	bestFit_.GetPoint(0, constTermTrue_tree,alphaTrue_tree);
// 	std::cout << constTerm_tree << " " << constTermTrue_tree 
// 		  << "\t" << alpha_tree << " " << alphaTrue_tree 
// 		  << std::endl;

	if(opt.Contains("scandiff")){
	  constTermTrue_tree = getMinimumFromTree(tree, "nll",TString(constTermName).ReplaceAll("-","_"));
	} else       if(opt.Contains("scan")){
	  constTerm_tree = getMinimumFromTree(tree, "nll",TString(constTermName).ReplaceAll("-","_"));
	}
	//std::cout << iBinX << "\t" << iBinY << "\t" << constTerm_tree - getMinimumFromTree(tree, "nll",TString(constTermName).ReplaceAll("-","_")) << std::endl;
	
	toys->Fill();
//       }else{
// 	hist->SaveAs("myhist.root");
// 	exit(0);
      }

      
      delete tree;
      delete hist;
    
    }
    f_in.Close();
//.........这里部分代码省略.........

void MakePlots(TString filename, float zmax=30, int nSmooth=10, TString opt="", TString energy="8TeV", TString lumi=""){
  TString outDir=filename; outDir.Remove(outDir.Last('/')); outDir+="/img/"+opt;
  //outDir="tmp/k5b/";
  //std::map<TString, TH2F *> deltaNLL_map;
  
  /*------------------------------ Plotto */
  TCanvas *c = new TCanvas("c","c");
  
  TFile f_in(filename, "read");
  if(f_in.IsZombie()){
    std::cerr << "File opening error: " << filename << std::endl;
    return;
  }
  
  TList *KeyList = f_in.GetListOfKeys();
  std::cout << KeyList->GetEntries() << std::endl;
  for(int i =0; i <  KeyList->GetEntries(); i++){
    c->Clear();
    TKey *key = (TKey *)KeyList->At(i);
    if(TString(key->GetClassName())!="RooDataSet") continue;
    RooDataSet *dataset = (RooDataSet *) key->ReadObj();
    if(dataset==NULL){
      std::cerr << "[WARNING] No dataset for " << key->GetName() << "\t" << key->GetTitle() << std::endl;
      continue;
    }
    TString constTermName = dataset->GetName();
    TString alphaName=constTermName; alphaName.ReplaceAll("constTerm","alpha");

    if(constTermName.Contains("absEta_1_1.4442-gold")) continue;
    if(constTermName.Contains("rho") || constTermName.Contains("phi")) continue;
        if(constTermName.Contains("scale")) continue;
    TTree *tree = dataset2tree(dataset);
    TGraphErrors bestFit_ = bestFit(tree, alphaName, constTermName);
    //    TString binning="(241,-0.0005,0.2405,60,0.00025,0.03025)"; 
    TString binning="(241,-0.0005,0.2405,301,-0.00005,0.03005)"; 
    
    TH2F *hist = prof2d(tree, constTermName, alphaName, "nll", binning, true,nSmooth, opt);
//     std::cout << "Bin width = " << hist->GetXaxis()->GetBinWidth(10) << "\t" << hist->GetYaxis()->GetBinWidth(10) << std::endl; 
//     std::cout << "Bin 1 center = " << hist->GetXaxis()->GetBinCenter(1) << "\t" << hist->GetYaxis()->GetBinCenter(1) << std::endl; 
//     std::cout << "Bin 10 center = " << hist->GetXaxis()->GetBinCenter(10) << "\t" << hist->GetYaxis()->GetBinCenter(10) << std::endl; 
//     return;
    hist->Draw("colz");
    hist->GetZaxis()->SetRangeUser(0,zmax);
    hist->GetXaxis()->SetRangeUser(0,0.15);
    hist->GetYaxis()->SetRangeUser(0,0.018);

    hist->GetXaxis()->SetTitle("#Delta S");
    hist->GetYaxis()->SetTitle("#Delta C");

    Int_t iBinX, iBinY;
    Double_t x,y;
    hist->GetBinWithContent2(0.0002,iBinX,iBinY,1,-1,1,-1,0.0000001);
    x= hist->GetXaxis()->GetBinCenter(iBinX);
    y= hist->GetYaxis()->GetBinCenter(iBinY);
    std::cout << "Best Fit: " << x << "\t" << y << std::endl;
    TGraph nllBestFit(1,&x,&y);

    TString fileName=outDir+"/"+constTermName;
    fileName+="-"; fileName+=nSmooth;
    
    nllBestFit.SetMarkerStyle(3);
    nllBestFit.SetMarkerColor(kRed);
    nllBestFit.Draw("P same");

    std::cout << fileName << std::endl;
    ofstream fout(fileName+".dat", ios_base::app);
    fout << constTermName << "\t" << x << "\t" << y << std::endl;


    c->SaveAs(fileName+".png");
    c->SaveAs(fileName+".eps");
    if(fileName.Contains("constTerm")) c->SaveAs(fileName+".C");

    fileName+="-zoom";
    hist->GetZaxis()->SetRangeUser(0,1);
    //hist->GetXaxis()->SetRangeUser(0.00,0.12);
    //hist->GetYaxis()->SetRangeUser(0,0.005);
    c->SaveAs(fileName+".png");
    c->SaveAs(fileName+".eps");


//     hist->SaveAs(outDir+"/deltaNLL-"+constTermName+".root");
//     hist->Draw("colz");
//     bestFit_.Draw("P same");
//     bestFit_.SetMarkerSize(2);


    
//     nllBestFit.SetMarkerStyle(3);
//     nllBestFit.SetMarkerColor(kRed);
//     TList* contour68 = contourFromTH2(hist, 0.68);

//     hist->Draw("colz");
//     hist->GetZaxis()->SetRangeUser(0,zmax);
//     //bestFit_.Draw("P same");
//     nllBestFit.Draw("P same");
//     //contour68->Draw("same");
    delete hist;
    RooAbsPdf *histPdf = NULL;
    if(!opt.Contains("keys")){
//.........这里部分代码省略.........

int main(int argc, char *argv[]){

	OptionParser(argc,argv);
	

	RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR);
	RooMsgService::instance().setSilentMode(true);


  
	system(Form("mkdir -p %s",outdir_.c_str()));

	vector<string> procs;
	split(infilenames_,infilenamesStr_,boost::is_any_of(","));
  
   TPython::Exec("import os,imp,re");
   for (unsigned int i =0 ; i<infilenames_.size() ; i++){
     TFile *infile =  TFile::Open(infilenames_[i].c_str());
	   string outname  =(string) TPython::Eval(Form("'%s'.split(\"/\")[-1].replace('root','_reduced.root')",infilenames_[i].c_str())); 
     TFile *outfile = TFile::Open(outname.c_str(),"RECREATE") ;
    TDirectory* saveDir = outfile->mkdir("tagsDumper");
    saveDir->cd();

    RooWorkspace *inWS = (RooWorkspace*) infile->Get("tagsDumper/cms_hgg_13TeV");
    RooRealVar *intLumi = (RooRealVar*)inWS->var("IntLumi");
    RooWorkspace *outWS = new RooWorkspace("cms_hgg_13TeV");
    outWS->import(*intLumi);
    std::list<RooAbsData*> data =  (inWS->allData()) ;
    std::cout <<" [INFO] Reading WS dataset contents: "<< std::endl;
        for (std::list<RooAbsData*>::const_iterator iterator = data.begin(), end = data.end(); iterator != end; ++iterator )  {
              RooDataSet *dataset = dynamic_cast<RooDataSet *>( *iterator );
              if (dataset) {
              RooDataSet *datasetReduced = (RooDataSet*) dataset->emptyClone(dataset->GetName(),dataset->GetName());
                TRandom3 r;
                r.Rndm();
                double x[dataset->numEntries()];
                r.RndmArray(dataset->numEntries(),x);
                int desiredEntries = floor(0.5+ dataset->numEntries()*fraction_);
                int modFraction = floor(0.5+ 1/fraction_);
                int finalEventCount=0;
                for (int j =0; j < dataset->numEntries() ; j++){
                    if( j%modFraction==0){
                      finalEventCount++;
                    }
                 }
                float average_weight= dataset->sumEntries()/finalEventCount;
                for (int j =0; j < dataset->numEntries() ; j++){
                    if( j%modFraction==0){
                    dataset->get(j);
                    datasetReduced->add(*(dataset->get(j)),average_weight);
                    }
                }
              float entriesIN =dataset->sumEntries();
              float entriesOUT =datasetReduced->sumEntries();
           if(verbose_){
              std::cout << "Original dataset " << *dataset <<std::endl;
              std::cout << "Reduced       dataset " << *datasetReduced <<std::endl;
              std::cout << "********************************************" <<std::endl;
              std::cout << "fraction (obs) : " << entriesOUT/entriesIN << std::endl;
              std::cout << "fraction (exp) : " << fraction_ << std::endl;
              std::cout << "********************************************" <<std::endl;
              std::cout << "" <<std::endl;
              std::cout << "" <<std::endl;
              std::cout << "" <<std::endl;
              std::cout << "********************************************" <<std::endl;
              }
               outWS->import(*datasetReduced);
                }
                
             RooDataHist *datahist = dynamic_cast<RooDataHist *>( *iterator );

              if (datahist) {
              RooDataHist *datahistOUT = (RooDataHist*) datahist->emptyClone(datahist->GetName(),datahist->GetName());
                TRandom3 r;
                r.Rndm();
                for (int j =0; j < datahist->numEntries() ; j++){
                    
                    datahistOUT->add(*(datahist->get(j)),datahist->weight());
                }
              float w =datahistOUT->sumEntries();
              float z =datahist->sumEntries();
           if(verbose_){
              std::cout << "Original datahist " << *datahist <<std::endl;
              std::cout << "Reduced  datahist " << *datahistOUT<<std::endl;
              std::cout << "********************************************" <<std::endl;
              std::cout << "WH fraction (obs) : " << w/(z) <<std::endl;
              std::cout << "WH fraction (exp) : " << fraction_ << std::endl;
              std::cout << "********************************************" <<std::endl;
              std::cout << "" <<std::endl;
              std::cout << "" <<std::endl;
              std::cout << "" <<std::endl;
              std::cout << "********************************************" <<std::endl;
              }
               outWS->import(*datahistOUT);
                }
                  }
   saveDir->cd();
   outWS->Write();
   outfile->Close();
   infile->Close();
//.........这里部分代码省略.........