C++ RooRealVar::setRange方法代码示例

pair<double,double> bkgEvPerGeV(RooWorkspace *work, int m_hyp, int cat, int spin=false){
  
  RooRealVar *mass = (RooRealVar*)work->var("CMS_hgg_mass");
  if (spin) mass = (RooRealVar*)work->var("mass");
  mass->setRange(100,180);
  RooAbsPdf *pdf = (RooAbsPdf*)work->pdf(Form("pdf_data_pol_model_8TeV_cat%d",cat));
  RooAbsData *data = (RooDataSet*)work->data(Form("data_mass_cat%d",cat));
  RooPlot *tempFrame = mass->frame();
  data->plotOn(tempFrame,Binning(80));
  pdf->plotOn(tempFrame);
  RooCurve *curve = (RooCurve*)tempFrame->getObject(tempFrame->numItems()-1);
  double nombkg = curve->Eval(double(m_hyp));
 
  RooRealVar *nlim = new RooRealVar(Form("nlim%d",cat),"",0.,0.,1.e5);
  //double lowedge = tempFrame->GetXaxis()->GetBinLowEdge(FindBin(double(m_hyp)));
  //double upedge  = tempFrame->GetXaxis()->GetBinUpEdge(FindBin(double(m_hyp)));
  //double center  = tempFrame->GetXaxis()->GetBinUpCenter(FindBin(double(m_hyp)));

  nlim->setVal(nombkg);
  mass->setRange("errRange",m_hyp-0.5,m_hyp+0.5);
  RooAbsPdf *epdf = 0;
  epdf = new RooExtendPdf("epdf","",*pdf,*nlim,"errRange");
		
  RooAbsReal *nll = epdf->createNLL(*data,Extended(),NumCPU(4));
  RooMinimizer minim(*nll);
  minim.setStrategy(0);
  minim.setPrintLevel(-1);
  minim.migrad();
  minim.minos(*nlim);
  
  double error = (nlim->getErrorLo(),nlim->getErrorHi())/2.;
  data->Print(); 
  return pair<double,double>(nombkg,error); 
}

void THSEventsPDF::adjustBinning(Int_t* offset1) const
{
   RooRealVar* xvar = fx_off ;
  if (!dynamic_cast<RooRealVar*>(xvar)) {
    coutE(InputArguments) << "RooDataHist::adjustBinning(" << GetName() << ") ERROR: dimension " << xvar->GetName() << " must be real" << endl ;
    assert(0) ;
  }
  Double_t xlo = xvar->getMin() ;
  Double_t xhi = xvar->getMax() ;
  //adjust bin range limits with new scale parameter
  //cout<<scale<<" "<<fMean<<" "<<xlo<<" "<<xhi<<endl;
  xlo=(xlo-fMean)/scale+fMean;
  xhi=(xhi-fMean)/scale+fMean;
  if(xvar->getBinning().lowBound()==xlo&&xvar->getBinning().highBound()==xhi) return;
  xvar->setRange(xlo,xhi) ;
  // Int_t xmin(0) ;
  // cout<<"THSEventsPDF::adjustBinning( "<<xlo <<" "<<xhi<<endl;
  //now adjust fitting range to bin limits??Possibly not
  if (fRHist->GetXaxis()->GetXbins()->GetArray()) {

    RooBinning xbins(fRHist->GetNbinsX(),fRHist->GetXaxis()->GetXbins()->GetArray()) ;

    Double_t tolerance = 1e-6*xbins.averageBinWidth() ;
    
    // Adjust xlo/xhi to nearest boundary
    Double_t xloAdj = xbins.binLow(xbins.binNumber(xlo+tolerance)) ;
    Double_t xhiAdj = xbins.binHigh(xbins.binNumber(xhi-tolerance)) ;
    xbins.setRange(xloAdj,xhiAdj) ;

    xvar->setBinning(xbins) ;
    if (fabs(xloAdj-xlo)>tolerance||fabs(xhiAdj-xhi)<tolerance) {
      coutI(DataHandling) << "RooDataHist::adjustBinning(" << GetName() << "): fit range of variable " << xvar->GetName() << " expanded to nearest bin boundaries: [" 
			  << xlo << "," << xhi << "] --> [" << xloAdj << "," << xhiAdj << "]" << endl ;
    }


  } else {

    RooBinning xbins(fRHist->GetXaxis()->GetXmin(),fRHist->GetXaxis()->GetXmax()) ;
    xbins.addUniform(fRHist->GetNbinsX(),fRHist->GetXaxis()->GetXmin(),fRHist->GetXaxis()->GetXmax()) ;

    Double_t tolerance = 1e-6*xbins.averageBinWidth() ;

    // Adjust xlo/xhi to nearest boundary
    Double_t xloAdj = xbins.binLow(xbins.binNumber(xlo+tolerance)) ;
    Double_t xhiAdj = xbins.binHigh(xbins.binNumber(xhi-tolerance)) ;
    xbins.setRange(xloAdj,xhiAdj) ;
    xvar->setRange(xloAdj,xhiAdj) ;
    //xvar->setRange(xlo,xhi) ;
 
  }
  return;
}

//
// set value and range for a variable in the workspace
//
void setValRange (RooWorkspace* workspace, const char* name, double val, double vmin, double vmax)
{
  RooRealVar* var = workspace->var(name);
  if ( var ) {
    if ( vmax>vmin )  var->setRange(vmin,vmax);
    var->setVal(val);
  }
}

MakeBiasStudy::MakeBiasStudy() {
  int Nmodels = 8;
  //RooRealVar* mass = ws->var("mass");
  RooRealVar *mass = new RooRealVar("mass","mass", 100,180);
  RooRealVar *nBkgTruth = new RooRealVar("TruthNBkg","", 0,1e9);
  //  RooAbsData* realData = ws->data("Data_Combined")->reduce( Form("evtcat==evtcat::%s",cat.Data()) );

  double Bias[Nmodels][Nmodels];
  double BiasE[Nmodels][Nmodels];
  MakeAICFits MakeAIC_Fits;
  for(int truthType = 0; truthType < Nmodels; truthType++){

    RooAbsPdf *truthPdf = MakeAIC_Fits.getBackgroundPdf(truthType,mass);
    RooExtendPdf *truthExtendedPdf = new RooExtendPdf("truthExtendedPdf","",*truthPdf,*nBkgTruth);
    //truthExtendedPdf.fitTo(*realData,RooFit::Strategy(0),RooFit::NumCPU(NUM_CPU),RooFit::Minos(kFALSE),RooFit::Extended(kTRUE));
    //truthExtendedPdf.fitTo(*realData,RooFit::Strategy(2),RooFit::NumCPU(NUM_CPU),RooFit::Minos(kFALSE),RooFit::Extended(kTRUE));
    
    double BiasWindow = 2.00;
    mass->setRange("biasRegion", mh-BiasWindow, mh+BiasWindow);
    double TruthFrac = truthExtendedPdf->createIntegral(mass,RooFit::Range("biasRegion"),RooFit::NormSet(*mass))->getVal();
    double NTruth = TruthFrac * nBkgTruth->getVal();
    double NTruthE = TruthFrac * nBkgTruth->getError();
    
    RooDataSet* truthbkg = truthPdf->generate(RooArgSet(*mass),nBkgTruth);

    for(int modelType = 0; modelType < Nmodels; modelType++){
      RooAbsPdf* ModelShape = MakeAIC_Fits.getBackgroundPdf(modelType,mass);
      RooRealVar *nBkgFit = new RooRealVar("FitNBkg", "", 0, 1e9);
      RooExtendPdf ModelExtendedPdf = new RooExtendPdf("ModelExtendedPdf", "",*ModelShape, *nBkgFit);
      ModelExtendedPdf.fitTo(truthbkg, RooFit::Strategy(0),RooFit::NumCPU(NUM_CPU),RooFit::Minos(kFALSE),RooFit::Extended(kTRUE));
      ModelExtendedPdf.fitTo(truthbkg, RooFit::Strategy(2),RooFit::NumCPU(NUM_CPU),RooFit::Minos(kFALSE),RooFit::Extended(kTRUE));

      double FitFrac = ModelExtendedPdf.createIntegral(mass,RooFit::Range("biasRegion"),RooFit::NormSet(mass))->getVal();
      double NFit = FitFrac * nBkgFit->getVal();
      double NFitE = FitFrac * nBkgFit->getError();

      Bias[truthType][modelType] = fabs(NFit - NTruth);
      BiasE[truthType][modelType] = fabs(NFitE - NTruthE);
    }

    
  }
  
  for(int i = 0; i < Nmodels; i++) { 
    std::cout <<  "===== Truth Model : " << MakeBiasStudy::Category(i) << " ===== " << std::endl; 
    for (int j = 0; j < Nmodels; j++) { 
      std::cout << "Fit Model: " << MakeBiasStudy::Category(j) << "  , Bias = " << Bias[i][j] << " +/- " << BiasE[i][j] << std::endl; 
    }   
  }
  

}

void FitterUtils::initiateParams(int nGenSignalZeroGamma, int nGenSignalOneGamma, int nGenSignalTwoGamma, RooRealVar const& expoConstGen, RooRealVar& nSignal, RooRealVar& nPartReco, 
      RooRealVar& nComb, RooRealVar& fracZero, RooRealVar& fracOne, RooRealVar& expoConst, RooRealVar&  nJpsiLeak, bool constPartReco, RooRealVar const& fracPartRecoSigma)
{
   TRandom rand;
   rand.SetSeed();

   int nGenSignal = nGenSignalZeroGamma + nGenSignalOneGamma + nGenSignalTwoGamma;

   double nGenSignal2;
   double nGenPartReco2;
   if(!constPartReco)
   {
      nGenSignal2 = rand.Uniform(nGenSignal-5*sqrt(nGenSignal), nGenSignal+5*sqrt(nGenSignal));
      nGenPartReco2 = rand.Uniform(nGenPartReco-5*sqrt(nGenPartReco), nGenPartReco+5*sqrt(nGenPartReco));
   }
   if(constPartReco)
   { 
      double nGenSigPartReco( nGenSignal+nGenPartReco );
      double nGenSigPartReco2( rand.Uniform( nGenSigPartReco-5*sqrt(nGenSigPartReco), nGenSigPartReco+5*sqrt(nGenSigPartReco) ) );
      double fracPartReco1( nGenPartReco/(1.*nGenSignal));
      double fracPartReco2( rand.Uniform(fracPartReco1-5*fracPartRecoSigma.getVal(), fracPartReco1+5*fracPartRecoSigma.getVal()) ); 

      nGenPartReco2 = fracPartReco2*nGenSigPartReco2 / (1+fracPartReco2); 
      nGenSignal2 = nGenSigPartReco2 / (1+fracPartReco2); 
   }
   double nGenComb2 = rand.Uniform(nGenComb-5*sqrt(nGenComb), nGenComb+5*sqrt(nGenComb));
   double nGenJpsiLeak2 = rand.Uniform(nGenJpsiLeak-5*sqrt(nGenJpsiLeak), nGenJpsiLeak+5*sqrt(nGenJpsiLeak));


   nSignal.setVal(nGenSignal2);
   nSignal.setRange(TMath::Max(0.,nGenSignal2-10.*sqrt(nGenSignal)) , nGenSignal2+10*sqrt(nGenSignal));

   nPartReco.setVal(nGenPartReco2);
   nPartReco.setRange(TMath::Max(0.,nGenPartReco2-10.*sqrt(nGenPartReco)), nGenPartReco2+10*sqrt(nGenPartReco));


   nComb.setVal(nGenComb2);
   nComb.setRange(TMath::Max(0.,nGenComb2-10.*sqrt(nGenComb)), nGenComb2+10*sqrt(nGenComb));

   nJpsiLeak.setVal(nGenJpsiLeak2);
   nJpsiLeak.setRange(TMath::Max(0., nGenJpsiLeak2-10*sqrt(nGenJpsiLeak)), nGenJpsiLeak2+10*sqrt(nGenJpsiLeak));

   double fracGenZero(nGenSignalZeroGamma/(1.*nGenSignal));
   double fracGenOne(nGenSignalOneGamma/(1.*nGenSignal));

   fracZero.setVal(rand.Gaus(fracGenZero, sqrt(nGenSignalZeroGamma)/(1.*nGenSignal))) ;
   fracZero.setRange(0., 1.);
   fracOne.setVal(rand.Gaus(fracGenOne, sqrt(nGenSignalOneGamma)/(1.*nGenSignal))) ;
   fracOne.setRange(0., 1.);

   expoConst.setVal(rand.Uniform( expoConstGen.getVal() - 5*expoConstGen.getError(), expoConstGen.getVal() + 5*expoConstGen.getError() ) );
   expoConst.setRange( expoConstGen.getVal() - 10*expoConstGen.getError(), expoConstGen.getVal() + 10*expoConstGen.getError() );
}

double NormalizedIntegral(RooAbsPdf & function, RooRealVar & integrationVar, double lowerLimit, double upperLimit){

  integrationVar.setRange("integralRange",lowerLimit,upperLimit);
  RooAbsReal* integral = function.createIntegral(integrationVar,NormSet(integrationVar),Range("integralRange"));


  double normlizedIntegralValue = integral->getVal();

  //  cout<<normlizedIntegralValue<<endl;


  return normlizedIntegralValue;


}

void FitterUtilsSimultaneousExpOfPolyTimesX::initiateParams(int nGenSignalZeroGamma, int nGenSignalOneGamma, int nGenSignalTwoGamma,
      RooRealVar& nKemu, RooRealVar& nSignal, RooRealVar& nPartReco,
      RooRealVar& nComb, RooRealVar& fracZero, RooRealVar& fracOne,
      RooRealVar&  nJpsiLeak, bool constPartReco, RooRealVar const& fracPartRecoSigma,
      RooRealVar& l1Kee, RooRealVar& l2Kee, RooRealVar& l3Kee, RooRealVar& l4Kee, RooRealVar& l5Kee,
      RooRealVar& l1Kemu, RooRealVar& l2Kemu, RooRealVar& l3Kemu, RooRealVar& l4Kemu, RooRealVar& l5Kemu,
      RooRealVar const& l1KeeGen, RooRealVar const& l2KeeGen, RooRealVar const& l3KeeGen, RooRealVar const& l4KeeGen, RooRealVar const& l5KeeGen 

      )
{
  FitterUtilsExpOfPolyTimesX::initiateParams(nGenSignalZeroGamma, nGenSignalOneGamma, nGenSignalTwoGamma,
           nSignal, nPartReco, nComb, fracZero, fracOne, nJpsiLeak, constPartReco, fracPartRecoSigma,
           l1Kee, l2Kee, l3Kee, l4Kee, l5Kee,
           l1KeeGen, l2KeeGen, l3KeeGen, l4KeeGen, l5KeeGen ); 



  TRandom rand;
  rand.SetSeed();

  nKemu.setVal(rand.Uniform(nGenKemu-5*sqrt(nGenKemu), nGenKemu+5*sqrt(nGenKemu)));
  nKemu.setRange(nGenKemu-10*sqrt(nGenKemu), nGenKemu+10*sqrt(nGenKemu));

  l1Kemu.setVal(rand.Uniform( l1KeeGen.getVal() - 5*l1KeeGen.getError(), l1KeeGen.getVal() + 5*l1KeeGen.getError() ) );
  l1Kemu.setRange( l1KeeGen.getVal() - 10*l1KeeGen.getError(), l1KeeGen.getVal() + 10*l1KeeGen.getError() );

  l2Kemu.setVal(rand.Uniform( l2KeeGen.getVal() - 5*l2KeeGen.getError(), l2KeeGen.getVal() + 5*l2KeeGen.getError() ) );
  l2Kemu.setRange( l2KeeGen.getVal() - 10*l2KeeGen.getError(), l2KeeGen.getVal() + 10*l2KeeGen.getError() );

  l3Kemu.setVal(rand.Uniform( l3KeeGen.getVal() - 5*l3KeeGen.getError(), l3KeeGen.getVal() + 5*l3KeeGen.getError() ) );
  l3Kemu.setRange( l3KeeGen.getVal() - 10*l3KeeGen.getError(), l3KeeGen.getVal() + 10*l3KeeGen.getError() );

  l4Kemu.setVal(rand.Uniform( l4KeeGen.getVal() - 5*l4KeeGen.getError(), l4KeeGen.getVal() + 5*l4KeeGen.getError() ) );
  l4Kemu.setRange( l4KeeGen.getVal() - 10*l4KeeGen.getError(), l4KeeGen.getVal() + 10*l4KeeGen.getError() );

  l5Kemu.setVal(rand.Uniform( l5KeeGen.getVal() - 5*l5KeeGen.getError(), l5KeeGen.getVal() + 5*l5KeeGen.getError() ) );
  l5Kemu.setRange( l5KeeGen.getVal() - 10*l5KeeGen.getError(), l5KeeGen.getVal() + 10*l5KeeGen.getError() );

}

int main(int argc, char **argv)
{
	bool printeff = true;
	string fc = "none";
	
	gROOT->ProcessLine(".x lhcbStyle.C");

	if(argc > 1)
	{
		for(int a = 1; a < argc; a++)
		{
			string arg = argv[a];
			string str = arg.substr(2,arg.length()-2);

			if(arg.find("-E")!=string::npos) fc = str;
			if(arg=="-peff") printeff = true;
		}
	}
	
	int nexp = 100;
	int nbins = 6;
	double q2min[] = {8.,15.,11.0,15,16,18};
	double q2max[] = {11.,20.,12.5,16,18,20};

	TString datafilename = "/afs/cern.ch/work/p/pluca/weighted/Lmumu/candLb.root";
	TreeReader * data = new TreeReader("candLb2Lmumu");
	data->AddFile(datafilename);
	TreeReader * datajpsi = new TreeReader("candLb2JpsiL");
	datajpsi->AddFile(datafilename);

	TFile * histFile = new TFile("Afb_bkgSys.root","recreate");

	string options = "-quiet-noPlot-lin-stdAxis-XM(#Lambda#mu#mu) (MeV/c^{2})-noCost-noParams";
	Analysis::SetPrintLevel("s");

	RooRealVar * cosThetaL = new RooRealVar("cosThetaL","cosThetaL",0.,-1.,1.);
	RooRealVar * cosThetaB = new RooRealVar("cosThetaB","cosThetaB",0.,-1.,1.);
	RooRealVar * MM = new RooRealVar("Lb_MassConsLambda","Lb_MassConsLambda",5621.,5400.,6000.);
	MM->setRange("Signal",5600,5640);
	RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR);

	//TGraphAsymmErrors * fL_vs_q2 = new TGraphAsymmErrors();
	//TCanvas * ceff = new TCanvas();

	RooCategory * samples = new RooCategory("samples","samples");
	samples->defineType("DD");
	samples->defineType("LL");

	RooRealVar * afb = new RooRealVar("afb","afb",0.,-0.75,0.75);
	RooRealVar * fL = new RooRealVar("fL","fL",0.6,0.,1.);
	TString afbLpdf = "((3./8.)*(1.-fL)*(1 + TMath::Power(cosThetaL,2)) + afb*cosThetaL + (3./4.)*fL*(1 - TMath::Power(cosThetaL,2)))";
	RooRealVar * afbB = new RooRealVar("afbB","afbB",0.,-0.5,0.5);
	TString afbBpdf = "(1 + 2*afbB*cosThetaB)";
	RooAbsPdf * teoPdf = new RooGenericPdf("teoPdf",afbLpdf,RooArgSet(*cosThetaL,*afb,*fL));
	RooAbsPdf * teoPdfB = new RooGenericPdf("teoPdfB",afbBpdf,RooArgSet(*cosThetaB,*afbB));

	TreeReader * mydata = datajpsi;
	Str2VarMap jpsiParsLL = getJpsiPars("LL", CutsDef::LLcut, histFile);
	Str2VarMap jpsiParsDD = getJpsiPars("DD", CutsDef::DDcut, histFile);

	vector<TH1 *> fLsysh, afbsysh, afbBsysh, fLsysh_frac, afbsysh_frac, afbBsysh_frac;

	for(int i = 0; i < nbins; i++)
	{
		TString q2name = ((TString)Form("q2_%4.2f_%4.2f",q2min[i],q2max[i])).ReplaceAll(".","");
		if(i>0) { mydata = data; MM->setRange(5400,6000); }
		else { q2name = "jpsi"; MM->setRange(5500,5850); }
		TString curq2cut = Form("TMath::Power(J_psi_1S_MM/1000,2) >= %e && TMath::Power(J_psi_1S_MM/1000,2) < %e",q2min[i],q2max[i]);	
		
		cout << "------------------- q2 bin: " << q2min[i] << " - " << q2max[i] << " -----------------------" << endl;

		/**               GET AND FIT EFFICIENCIES                  **/

		RooAbsPdf * effDDpdf = NULL, * effLLpdf = NULL, * effLLBpdf = NULL, * effDDBpdf = NULL;	
		getEfficiencies(q2min[i],q2max[i],&effLLpdf,&effDDpdf,&effLLBpdf,&effDDBpdf,printeff);
		cout << "Efficiencies extracted" << endl;
		histFile->cd();


		/**                    FIT AFB                  **/


		afb->setVal(0);
		afbB->setVal(-0.37);
		fL->setVal(0.6);

		RooAbsPdf * corrPdfLL = new RooProdPdf("sigPdfLL"+q2name,"corrPdfLL",*teoPdf,*effLLpdf);
		RooAbsPdf * corrPdfDD = new RooProdPdf("sigPdfDD"+q2name,"corrPdfDD",*teoPdf,*effDDpdf);
		RooAbsPdf * corrPdfLLB = new RooProdPdf("sigPdfLLB"+q2name,"corrPdfLLB",*teoPdfB,*effLLBpdf);
		RooAbsPdf * corrPdfDDB = new RooProdPdf("sigPdfDDB"+q2name,"corrPdfDDB",*teoPdfB,*effDDBpdf);

		TCut baseCut = "";
		TCut cutLL = CutsDef::LLcut + (TCut)curq2cut + baseCut;
		TCut cutDD = CutsDef::DDcut + (TCut)curq2cut + baseCut;

		histFile->cd();
		double fracDDv[2], fracLLv[2];
		double nsigDD, nsigLL;
		RooDataSet * dataLL = getDataAndFrac("LL",q2name,mydata,cutLL,MM,&fracLLv[0],jpsiParsLL,&nsigLL);
		RooDataSet * dataDD = getDataAndFrac("DD",q2name,mydata,cutDD,MM,&fracDDv[0],jpsiParsDD,&nsigDD);
//.........这里部分代码省略.........

void FitterUtilsSimultaneousExpOfPolyTimesX::generate(bool wantPlots, string plotsfile)
{
   FitterUtilsExpOfPolyTimesX::generate(wantPlots, plotsfile);
   TFile fw(workspacename.c_str(), "UPDATE");
   RooWorkspace* workspace = (RooWorkspace*)fw.Get("workspace");

   RooRealVar *B_plus_M = workspace->var("B_plus_M");
   RooRealVar *misPT = workspace->var("misPT");
   RooDataSet* dataSetCombExt = (RooDataSet*)workspace->data("dataSetCombExt");
   RooDataSet* dataSetComb = (RooDataSet*)workspace->data("dataSetComb");
//   RooRealVar *l1KeeGen = workspace->var("l1KeeGen");  
//   RooRealVar *l2KeeGen = workspace->var("l2KeeGen");  
//   RooRealVar *l3KeeGen = workspace->var("l3KeeGen");  
//   RooRealVar *l4KeeGen = workspace->var("l4KeeGen");  
//   RooRealVar *l5KeeGen = workspace->var("l5KeeGen");  
//
//
//   RooExpOfPolyTimesX kemuPDF("kemuPDF", "kemuPDF",  *B_plus_M, *misPT,  *l1KeeGen, *l2KeeGen, *l3KeeGen, *l4KeeGen, *l5KeeGen);
//
//   RooAbsPdf::GenSpec* GenSpecKemu = kemuPDF.prepareMultiGen(RooArgSet(*B_plus_M, *misPT), RooFit::Extended(1), NumEvents(nGenKemu));
//
//   cout<<"Generating Kemu"<<endl;
//   RooDataSet* dataGenKemu = kemuPDF.generate(*GenSpecKemu);//(argset, 100, false, true, "", false, true);
//   dataGenKemu->SetName("dataGenKemu"); dataGenKemu->SetTitle("dataGenKemu");
//
//
//   RooWorkspace* workspaceGen = (RooWorkspace*)fw.Get("workspaceGen");
//   workspaceGen->import(*dataGenKemu);
//
//   workspaceGen->Write("", TObject::kOverwrite);
//   fw.Close();
//   delete dataGenKemu;
//   delete GenSpecKemu;


   TVectorD rho(2);
   rho[0] = 2.5;
   rho[1] = 1.5;
   misPT->setRange(-2000, 5000);
   RooNDKeysPdf kemuPDF("kemuPDF", "kemuPDF", RooArgList(*B_plus_M, *misPT), *dataSetCombExt, rho, "ma",3, true);
   misPT->setRange(0, 5000);


   RooAbsPdf::GenSpec* GenSpecKemu = kemuPDF.prepareMultiGen(RooArgSet(*B_plus_M, *misPT), RooFit::Extended(1), NumEvents(nGenKemu));

   cout<<"Generating Kemu"<<endl;
   RooDataSet* dataGenKemu = kemuPDF.generate(*GenSpecKemu);//(argset, 100, false, true, "", false, true);
   dataGenKemu->SetName("dataGenKemu"); dataGenKemu->SetTitle("dataGenKemu");


   RooWorkspace* workspaceGen = (RooWorkspace*)fw.Get("workspaceGen");
   workspaceGen->import(*dataGenKemu);

   if(wantPlots) PlotShape(*dataSetComb, *dataGenKemu, kemuPDF, plotsfile, "cKemuKeys", *B_plus_M, *misPT);

   fw.cd();
   workspaceGen->Write("", TObject::kOverwrite);
   fw.Close();
   delete dataGenKemu;
   delete GenSpecKemu;
}

void setup(ModelConfig* mcInWs) {
  RooAbsPdf* combPdf = mcInWs->GetPdf();

  RooArgSet mc_obs = *mcInWs->GetObservables();
  RooArgSet mc_globs = *mcInWs->GetGlobalObservables();
  RooArgSet mc_nuis = *mcInWs->GetNuisanceParameters();

  // pair the nuisance parameter to the global observable
  RooArgSet mc_nuis_tmp = mc_nuis;
  RooArgList nui_list;
  RooArgList glob_list;
  RooArgSet constraint_set_tmp(*combPdf->getAllConstraints(mc_obs, mc_nuis_tmp, false));
  RooArgSet constraint_set;
  int counter_tmp = 0;
  unfoldConstraints(constraint_set_tmp, constraint_set, mc_obs, mc_nuis_tmp, counter_tmp);

  TIterator* cIter = constraint_set.createIterator();
  RooAbsArg* arg;
  while ((arg = (RooAbsArg*)cIter->Next())) {
    RooAbsPdf* pdf = (RooAbsPdf*)arg;
    if (!pdf) continue;

    // pdf->Print();

    TIterator* nIter = mc_nuis.createIterator();
    RooRealVar* thisNui = NULL;
    RooAbsArg* nui_arg;
    while ((nui_arg = (RooAbsArg*)nIter->Next())) {
      if (pdf->dependsOn(*nui_arg)) {
        thisNui = (RooRealVar*)nui_arg;
        break;
      }
    }
    delete nIter;

    // need this incase the observable isn't fundamental. 
    // in this case, see which variable is dependent on the nuisance parameter and use that.
    RooArgSet* components = pdf->getComponents();
    // components->Print();
    components->remove(*pdf);
    if (components->getSize()) {
      TIterator* itr1 = components->createIterator();
      RooAbsArg* arg1;
      while ((arg1 = (RooAbsArg*)itr1->Next())) {
        TIterator* itr2 = components->createIterator();
        RooAbsArg* arg2;
        while ((arg2 = (RooAbsArg*)itr2->Next())) {
          if (arg1 == arg2) continue;
          if (arg2->dependsOn(*arg1)) {
            components->remove(*arg1);
          }
        }
        delete itr2;
      }
      delete itr1;
    }

    if (components->getSize() > 1) {
      cout << "ERROR::Couldn't isolate proper nuisance parameter" << endl;
      return;
    }
    else if (components->getSize() == 1) {
      thisNui = (RooRealVar*)components->first();
    }

    TIterator* gIter = mc_globs.createIterator();
    RooRealVar* thisGlob = NULL;
    RooAbsArg* glob_arg;
    while ((glob_arg = (RooAbsArg*)gIter->Next())) {
      if (pdf->dependsOn(*glob_arg)) {
        thisGlob = (RooRealVar*)glob_arg;
        break;
      }
    }
    delete gIter;

    if (!thisNui || !thisGlob) {
      cout << "WARNING::Couldn't find nui or glob for constraint: " << pdf->GetName() << endl;
      //return;
      continue;
    }

    // cout << "Pairing nui: " << thisNui->GetName() << ", with glob: " << thisGlob->GetName() << ", from constraint: " << pdf->GetName() << endl;

    nui_list.add(*thisNui);
    glob_list.add(*thisGlob);

    if (string(pdf->ClassName()) == "RooPoisson")  {
      double minVal = max(0.0, thisGlob->getVal() - 8*sqrt(thisGlob->getVal()));
      double maxVal = max(10.0, thisGlob->getVal() + 8*sqrt(thisGlob->getVal()));
      thisNui->setRange(minVal, maxVal);
      thisGlob->setRange(minVal, maxVal);
    }
    else if (string(pdf->ClassName()) == "RooGaussian") {
      thisNui->setRange(-7, 7);
      thisGlob->setRange(-10, 10);
    }

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

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

 doofit::builder::EasyPdf *epdf = new doofit::builder::EasyPdf();

    

 epdf->Var("sig_yield");
 epdf->Var("sig_yield").setVal(153000);
 epdf->Var("sig_yield").setConstant(false);
 //decay time
 epdf->Var("obsTime");
 epdf->Var("obsTime").SetTitle("t_{#kern[-0.2]{B}_{#kern[-0.1]{ d}}^{#kern[-0.1]{ 0}}}");
 epdf->Var("obsTime").setUnit("ps");
 epdf->Var("obsTime").setRange(0.,16.);

 // tag, respectively the initial state of the produced B meson
 epdf->Cat("obsTag");
 epdf->Cat("obsTag").defineType("B_S",1);
 epdf->Cat("obsTag").defineType("Bbar_S",-1);

  //finalstate
  epdf->Cat("catFinalState");
  epdf->Cat("catFinalState").defineType("f",1);
  epdf->Cat("catFinalState").defineType("fbar",-1);

  epdf->Var("obsEtaOS");
  epdf->Var("obsEtaOS").setRange(0.0,0.5);


  std::vector<double> knots;
    knots.push_back(0.07);
    knots.push_back(0.10);
    knots.push_back(0.138);
    knots.push_back(0.16);
    knots.push_back(0.23);
    knots.push_back(0.28);
    knots.push_back(0.35);
    knots.push_back(0.42);
    knots.push_back(0.44);
    knots.push_back(0.48);
    knots.push_back(0.5);

  // empty arg list for coefficients
  RooArgList* list = new RooArgList();

  // create first coefficient
  RooRealVar* coeff_first = &(epdf->Var("parCSpline1"));
  coeff_first->setRange(0,10000);
  coeff_first->setVal(1);
  coeff_first->setConstant(false);
  list->add( *coeff_first );

  for (unsigned int i=1; i <= knots.size(); ++i){
    std::string number = boost::lexical_cast<std::string>(i);
    RooRealVar* coeff = &(epdf->Var("parCSpline"+number));
    coeff->setRange(0,10000);
    coeff->setVal(1);
    coeff->setConstant(false);
    list->add( *coeff );
    }
  
  // create last coefficient
  RooRealVar* coeff_last = &(epdf->Var("parCSpline"+boost::lexical_cast<std::string>(knots.size())));
  coeff_last->setRange(0,10000);
  coeff_last->setVal(1);
  coeff_last->setConstant(false);
  list->add( *coeff_last );
  list->Print();
  // define Eta PDF
  doofit::roofit::pdfs::DooCubicSplinePdf splinePdf("splinePdf",epdf->Var("obsEtaOS"),knots,*list,0,0.5);
  
  //Berechne die Tagging Assymetrie
  epdf->Var("p0");
  epdf->Var("p0").setVal(0.369);
  epdf->Var("p0").setConstant(true);

  epdf->Var("p1");
  epdf->Var("p1").setVal(0.952);
  epdf->Var("p1").setConstant(true);

  epdf->Var("delta_p0");
  epdf->Var("delta_p0").setVal(0.019);
  epdf->Var("delta_p0").setConstant(true);

  epdf->Var("delta_p1");
  epdf->Var("delta_p1").setVal(-0.012);
  epdf->Var("delta_p1").setConstant(true);

  epdf->Var("etamean");
  epdf->Var("etamean").setVal(0.365);
  epdf->Var("etamean").setConstant(true);

  epdf->Formula("omega","@0 [email protected]/2 +(@[email protected]/2)*(@[email protected])", RooArgList(epdf->Var("p0"),epdf->Var("delta_p0"),epdf->Var("p1"),epdf->Var("delta_p1"),epdf->Var("obsEtaOS"),epdf->Var("etamean")));
  epdf->Formula("omegabar","@0 [email protected]/2 +(@[email protected]/2)*(@[email protected])", RooArgList(epdf->Var("p0"),epdf->Var("delta_p0"),epdf->Var("p1"),epdf->Var("delta_p1"),epdf->Var("obsEtaOS"),epdf->Var("etamean")));
      


  //Koeffizienten
  DecRateCoeff *coeff_c = new DecRateCoeff("coef_cos","coef_cos",DecRateCoeff::CPOdd,epdf->Cat("catFinalState"),epdf->Cat("obsTag"),epdf->Var("C_f"),epdf->Var("C_fbar"),epdf->Var("obsEtaOS"),splinePdf,epdf->Var("tageff"),epdf->Real("omega"),epdf->Real("omegabar"),epdf->Var("asym_prod"),epdf->Var("asym_det"),epdf->Var("asym_tageff"));
  DecRateCoeff *coeff_s = new DecRateCoeff("coef_sin","coef_sin",DecRateCoeff::CPOdd,epdf->Cat("catFinalState"),epdf->Cat("obsTag"),epdf->Var("S_f"),epdf->Var("S_fbar"),epdf->Var("obsEtaOS"),splinePdf,epdf->Var("tageff"),epdf->Real("omega"),epdf->Real("omegabar"),epdf->Var("asym_prod"),epdf->Var("asym_det"),epdf->Var("asym_tageff"));
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  cout << "  a11_bkgd = " << CMS_qqzzbkg_a11.getVal() << endl;
  cout << "  a12_bkgd = " << CMS_qqzzbkg_a12.getVal() << endl;
  cout << "  a13_bkgd = " << CMS_qqzzbkg_a13.getVal() << endl;
  cout << "}" << endl;
  cout << "---------------------------" << endl;


  of << "qqZZshape a0_bkgd   " << CMS_qqzzbkg_a0.getVal() << endl;
  of << "qqZZshape a1_bkgd   " << CMS_qqzzbkg_a1.getVal() << endl;
  of << "qqZZshape a2_bkgd   " << CMS_qqzzbkg_a2.getVal() << endl;
  of << "qqZZshape a3_bkgd   " << CMS_qqzzbkg_a3.getVal() << endl;
  of << "qqZZshape a4_bkgd   " << CMS_qqzzbkg_a4.getVal() << endl;
  of << "qqZZshape a5_bkgd   " << CMS_qqzzbkg_a5.getVal() << endl;
  of << "qqZZshape a6_bkgd   " << CMS_qqzzbkg_a6.getVal() << endl;
  of << "qqZZshape a7_bkgd   " << CMS_qqzzbkg_a7.getVal() << endl;
  of << "qqZZshape a8_bkgd   " << CMS_qqzzbkg_a8.getVal() << endl;
  of << "qqZZshape a9_bkgd   " << CMS_qqzzbkg_a9.getVal() << endl;
  of << "qqZZshape a10_bkgd  " << CMS_qqzzbkg_a10.getVal() << endl;
  of << "qqZZshape a11_bkgd  " << CMS_qqzzbkg_a11.getVal() << endl;
  of << "qqZZshape a12_bkgd  " << CMS_qqzzbkg_a12.getVal() << endl;
  of << "qqZZshape a13_bkgd  " << CMS_qqzzbkg_a13.getVal() << endl;
  of << endl << endl;
  of.close();

  cout << endl << "Output written to: " << outfile << endl;
  
    
  double qqzznorm;
  if (channel == 1) qqzznorm = 20.5836;
  else if (channel == 2) qqzznorm = 13.8871;
  else if (channel == 3) qqzznorm = 32.9883;
  else { cout << "disaster!" << endl; }

  ZZMass->setRange("fullrange",100.,1000.);
  ZZMass->setRange("largerange",100.,600.);
  ZZMass->setRange("zoomrange",100.,200.);
    
  double rescale = qqzznorm/totalweight;
  double rescale_z = qqzznorm/totalweight_z;
  cout << "rescale: " << rescale << ", rescale_z: " << rescale_z << endl;


  // Plot m4l and
  RooPlot* frameM4l = ZZMass->frame(Title("M4L"),Range(100,600),Bins(250)) ;
  set->plotOn(frameM4l, MarkerStyle(20), Rescale(rescale)) ;
  
  //set->plotOn(frameM4l) ;
  RooPlot* frameM4lz = ZZMass->frame(Title("M4L"),Range(100,200),Bins(100)) ;
  set->plotOn(frameM4lz, MarkerStyle(20), Rescale(rescale)) ;


  int iLineColor = 1;
  string lab = "blah";
  if (channel == 1) { iLineColor = 2; lab = "4#mu"; }
  if (channel == 3) { iLineColor = 4; lab = "2e2#mu"; }
  if (channel == 2) { iLineColor = 6; lab = "4e"; }

  bkg_qqzz->plotOn(frameM4l,LineColor(iLineColor),NormRange("largerange")) ;
  bkg_qqzz->plotOn(frameM4lz,LineColor(iLineColor),NormRange("zoomrange")) ;
    
//second shape to compare with (if previous comparison code unceommented)
  //bkg_qqzz_bkgd->plotOn(frameM4l,LineColor(1),NormRange("largerange")) ;
  //bkg_qqzz_bkgd->plotOn(frameM4lz,LineColor(1),NormRange("zoomrange")) ;
    
  
  double normalizationBackground_qqzz = bkg_qqzz->createIntegral( RooArgSet(*ZZMass), Range("fullrange") )->getVal();

//void RunToyScan5(TString fileName, double startVal, double stopVal, TString outFile) {
void frequentist(TString fileName) {
  cout << "Starting frequentist " << time(NULL) << endl;
  double startVal = 0;
  double stopVal = 200;
  TString outFile = "";

  int nToys = 1 ;
  int nscanpoints = 2 ;

  /*
  gROOT->LoadMacro("RooBetaPdf.cxx+") ;
  gROOT->LoadMacro("RooRatio.cxx+") ;
  gROOT->LoadMacro("RooPosDefCorrGauss.cxx+") ;
  */

  // get relevant objects out of the "ws" file

  TFile *file = TFile::Open(fileName);
  if(!file){
    cout <<"file not found" << endl;
    return;
  } 

  RooWorkspace* w = (RooWorkspace*) file->Get("workspace");
  if(!w){
    cout <<"workspace not found" << endl;
    return;
  }

  ModelConfig* mc = (ModelConfig*) w->obj("S+B_model");
  RooAbsData* data = w->data("data");

  if( !data || !mc ){
    w->Print();
    cout << "data or ModelConfig was not found" <<endl;
    return;
  }

  RooRealVar* myPOI = (RooRealVar*) mc->GetParametersOfInterest()->first();
  myPOI->setRange(0, 1000.);

  ModelConfig* bModel = (ModelConfig*) w->obj("B_model");
  ModelConfig* sbModel = (ModelConfig*) w->obj("S+B_model");

  ProfileLikelihoodTestStat profll(*sbModel->GetPdf());
  profll.SetPrintLevel(2);
  profll.SetOneSided(1);
  TestStatistic * testStat = &profll;

  HypoTestCalculatorGeneric *  hc = 0;
  hc = new FrequentistCalculator(*data, *bModel, *sbModel);
  
  ToyMCSampler *toymcs = (ToyMCSampler*)hc->GetTestStatSampler();
  toymcs->SetMaxToys(10000);
  toymcs->SetNEventsPerToy(1);
  toymcs->SetTestStatistic(testStat);


  ((FrequentistCalculator *)hc)->SetToys(nToys,nToys);
  
  HypoTestInverter calc(*hc);
  calc.SetConfidenceLevel(0.95);
  calc.UseCLs(true);
  //calc.SetVerbose(true);
  calc.SetVerbose(2);

  cout << "About to set fixed scan " << time(NULL) << endl;
  calc.SetFixedScan(nscanpoints,startVal,stopVal);
  cout << "About to do inverter " << time(NULL) << endl;
  HypoTestInverterResult * res_toysCLs_calculator = calc.GetInterval();

  cout << "CLs = " << res_toysCLs_calculator->UpperLimit() 
	    << "   CLs_exp = " << res_toysCLs_calculator->GetExpectedUpperLimit(0) 
	    << "   CLs_exp(-1s) = " << res_toysCLs_calculator->GetExpectedUpperLimit(-1) 
	    << "   CLs_exp(+1s) = " << res_toysCLs_calculator->GetExpectedUpperLimit(1) << endl ;

  /*
  // dump results string to output file
  ofstream outStream ;
  outStream.open(outFile,ios::app) ;
  
  outStream << "CLs = " << res_toysCLs_calculator->UpperLimit() 
	    << "   CLs_exp = " << res_toysCLs_calculator->GetExpectedUpperLimit(0) 
	    << "   CLs_exp(-1s) = " << res_toysCLs_calculator->GetExpectedUpperLimit(-1) 
	    << "   CLs_exp(+1s) = " << res_toysCLs_calculator->GetExpectedUpperLimit(1) << endl ;
  
  outStream.close() ;
  */


  cout << "End of frequentist " << time(NULL) << endl;
  return ;

}

//.........这里部分代码省略.........
	smallTree->Branch( "L2PT", &l2pt, "L2PT/D" );
	smallTree->Branch( "L2ETA", &l2eta, "L2ETA/D" );
	smallTree->Branch( "L2PHI", &l2phi, "L2PHI/D" );
	smallTree->Branch( "ZMASS", &zmass, "ZMASS/D" );
	smallTree->Branch( "ZPT", &zpt, "ZPT/D" );
	smallTree->Branch( "ZETA", &zeta, "ZETA/D" );
	smallTree->Branch( "MT", &mt, "MT/D" );
	smallTree->Branch( "NSOFTJET", &nsoftjet, "NSOFTJET/I" );
	smallTree->Branch( "NHARDJET", &nhardjet, "NHARDJET/I" );
	smallTree->Branch( "MAXJETBTAG", &maxJetBTag, "MAXJETBTAG/D" );
	smallTree->Branch( "MINDPJETMET", &minDeltaPhiJetMet, "MINDPJETMET/D" );
	smallTree->Branch( "DETAJJ", &detajj, "DETAJJ/D" );
	smallTree->Branch( "MJJ", &mjj, "MJJ/D" );
	smallTree->Branch( "NVTX", &nvtx, "NVTX/I" );
	smallTree->Branch( "nInter" , &ni, "nInter/I" );
	smallTree->Branch( "CATEGORY", &category, "CATEGORY/I" );
	smallTree->Branch( "Weight" , &weight, "Weight/D" );
	smallTree->Branch( "HMASS", &hmass, "HMASS/D" );
	smallTree->Branch( "HWEIGHT", &hweight, "HWEIGHT/D" );

	bool isData = opt.checkBoolOption("DATA");

	unsigned long nentries = tree->GetEntries();

	RooDataSet * events = nullptr;

	PhotonPrescale photonPrescales;

	vector<int> thresholds;
	if (type == PHOT) {
		if (w == nullptr)
			throw string("ERROR: No mass peak pdf!");
		RooRealVar * zmass = w->var("mass");
		zmass->setRange(76.0, 106.0);
		RooAbsPdf * pdf = w->pdf("massPDF");
		events = pdf->generate(*zmass, nentries);

		photonPrescales.addTrigger("HLT_Photon36_R9Id90_HE10_Iso40_EBOnly", 36, 3, 7);
		photonPrescales.addTrigger("HLT_Photon50_R9Id90_HE10_Iso40_EBOnly", 50, 5, 8);
		photonPrescales.addTrigger("HLT_Photon75_R9Id90_HE10_Iso40_EBOnly", 75, 7, 9);
		photonPrescales.addTrigger("HLT_Photon90_R9Id90_HE10_Iso40_EBOnly", 90, 10, 10);
	}

	TH1D ptSpectrum("ptSpectrum", "ptSpectrum", 200, 55, 755);
	ptSpectrum.Sumw2();

	unordered_set<EventAdr> eventsSet;
	for ( unsigned long iEvent = 0; iEvent < nentries; iEvent++ ) {
//		if (iEvent < 6060000)
//			continue;

		if ( iEvent % 10000 == 0) {
			cout << string(40, '\b');
			cout << setw(10) << iEvent << " / " << setw(10) << nentries << " done ..." << std::flush;
		}

		tree->GetEntry( iEvent );

		run = -999;
		lumi = -999;
		event = -999;
		pfmet = -999;
		nele = -999;
		nmu = -999;
		nsoftmu = -999;
		l1pt = -999;

TH1D* runSig(RooWorkspace* ws,
             const char* modelConfigName = "ModelConfig",
             const char* dataName = "obsData",
             const char* asimov1DataName = "asimovData_1",
             const char* conditional1Snapshot = "conditionalGlobs_1",
             const char* nominalSnapshot = "nominalGlobs")
{
    string defaultMinimizer    = "Minuit";     // or "Minuit"
    int defaultStrategy        = 2;             // Minimization strategy. 0-2. 0 = fastest, least robust. 2 = slowest, most robust

    double mu_profile_value = 1; // mu value to profile the obs data at wbefore generating the expected
    bool doUncap            = 1; // uncap p0
    bool doInj              = 0; // setup the poi for injection study (zero is faster if you're not)
    bool doMedian           = 1; // compute median significance
    bool isBlind            = 0; // Dont look at observed data
    bool doConditional      = !isBlind; // do conditional expected data
    bool doObs              = !isBlind; // compute observed significance

    TStopwatch timer;
    timer.Start();

    if (!ws)
    {
        cout << "ERROR::Workspace is NULL!" << endl;
        return NULL;
    }
    ModelConfig* mc = (ModelConfig*)ws->obj(modelConfigName);
    if (!mc)
    {
        cout << "ERROR::ModelConfig: " << modelConfigName << " doesn't exist!" << endl;
        return NULL;
    }
    RooDataSet* data = (RooDataSet*)ws->data(dataName);
    if (!data)
    {
        cout << "ERROR::Dataset: " << dataName << " doesn't exist!" << endl;
        return NULL;
    }

    mc->GetNuisanceParameters()->Print("v");

    //RooNLLVar::SetIgnoreZeroEntries(1);
    ROOT::Math::MinimizerOptions::SetDefaultMinimizer(defaultMinimizer.c_str());
    ROOT::Math::MinimizerOptions::SetDefaultStrategy(defaultStrategy);
    ROOT::Math::MinimizerOptions::SetDefaultPrintLevel(-1);
    //  cout << "Setting max function calls" << endl;
    //ROOT::Math::MinimizerOptions::SetDefaultMaxFunctionCalls(20000);
    //RooMinimizer::SetMaxFunctionCalls(10000);

    ws->loadSnapshot("conditionalNuis_0");
    RooArgSet nuis(*mc->GetNuisanceParameters());

    RooRealVar* mu = (RooRealVar*)mc->GetParametersOfInterest()->first();

    RooAbsPdf* pdf_temp = mc->GetPdf();

    string condSnapshot(conditional1Snapshot);
    RooArgSet nuis_tmp2 = *mc->GetNuisanceParameters();
    RooNLLVar* obs_nll = doObs ? (RooNLLVar*)pdf_temp->createNLL(*data, Constrain(nuis_tmp2)) : NULL;

    RooDataSet* asimovData1 = (RooDataSet*)ws->data(asimov1DataName);
    if (!asimovData1)
    {
        cout << "Asimov data doesn't exist! Please, allow me to build one for you..." << endl;
        string mu_str, mu_prof_str;

        asimovData1 = makeAsimovData(mc, doConditional, ws, obs_nll, 1, &mu_str, &mu_prof_str, mu_profile_value, true);
        condSnapshot="conditionalGlobs"+mu_prof_str;

        //makeAsimovData(mc, true, ws, mc->GetPdf(), data, 0);
        //ws->Print();
        //asimovData1 = (RooDataSet*)ws->data("asimovData_1");
    }

    if (!doUncap) mu->setRange(0, 40);
    else mu->setRange(-40, 40);

    RooAbsPdf* pdf = mc->GetPdf();
    RooArgSet nuis_tmp1 = *mc->GetNuisanceParameters();
    RooNLLVar* asimov_nll = (RooNLLVar*)pdf->createNLL(*asimovData1, Constrain(nuis_tmp1));

    //do asimov
    mu->setVal(1);
    mu->setConstant(0);
    if (!doInj) mu->setConstant(1);

    int status,sign;
    double med_sig=0,obs_sig=0,asimov_q0=0,obs_q0=0;

    if (doMedian)
    {
        ws->loadSnapshot(condSnapshot.c_str());
        if (doInj) ws->loadSnapshot("conditionalNuis_inj");
        else ws->loadSnapshot("conditionalNuis_1");
        mc->GetGlobalObservables()->Print("v");
        mu->setVal(0);
        mu->setConstant(1);
        status = minimize(asimov_nll, ws);
        if (status >= 0) cout << "Success!" << endl;

//.........这里部分代码省略.........