C++ TRandom3::Poisson方法代码示例

// ------------------------------------------------------------------------------------------------------//
// ------------------------------------------------------------------------------------------------------//
void bias_corr_stat_error(double Sa, double Sb, double Sc, double Sd,
		          double Qa, double Qb, double Qc, double Qd, 
		          double Scale, double *ret){
 
      	double ea = (Sa)/(Sa+Sb);
      	double ed = (Sd)/(Sd+Sc);
      	double ew = (Sa+Sb)/(Sa+Sb+Sd+Sc);
 	double st = Sa+Sb+Sc+Sd;
 
 	double na = (Qa + Sa);
      	double nb = (Qb + Sb);
      	double nc = (Qc + Sc);     		
      	double nd = (Qd + Sd);
 	
 	double S[4] = {0,0,0,0};
 	signal(na,nb,nc,nd,ea,ed,ew,S);
 	
 	double s1 = S[0];
 	double s2 = S[2];
 	
 	double cen1 = st/s1;
 	double cen2 = st/s2;
 	
      	TH1F *tmp_u = new TH1F("tmp_u","tmp_u",100,cen1*(1-0.01),cen1*(1+0.01)); 
      	TH1F *tmp_d = new TH1F("tmp_d","tmp_d",100,cen2*(1-0.01),cen2*(1+0.01));
      	
      	TRandom3 *rnd = new TRandom3();
      	      	
      	for(int i=0;i<50000; i++){
      	
      		double qa = Scale*rnd->Poisson((1./Scale)*Qa);
      		double qb = Scale*rnd->Poisson((1./Scale)*Qb);
      		double qc = Scale*rnd->Poisson((1./Scale)*Qc);
      		double qd = Scale*rnd->Poisson((1./Scale)*Qd);
      		
      		na = (qa + Sa);
      		nb = (qb + Sb);
      		nc = (qc + Sc);     		
      		nd = (qd + Sd);
      		
      		signal(na,nb,nc,nd,ea,ed,ew,S);

      		tmp_u->Fill(st/S[0]);      		     
      		tmp_d->Fill(st/S[2]);      		               	
          	
        }
        
	tmp_u->Fit("gaus","Q","Q");
	tmp_d->Fit("gaus","Q","Q");

	double val1 = tmp_u->GetFunction("gaus")->GetParameter(2);
	double val2 = tmp_d->GetFunction("gaus")->GetParameter(2);
	delete tmp_u;
	delete tmp_d;
	
	ret[0] = val1;
	ret[1] = val2;
}

int main(int argc, char *argv[])
{
    std::auto_ptr<TRint> application(new TRint("histInMemory", 0, 0));

    TH1F *hist1 = new TH1F("hist1", "hist1", 50, 0, 10);
    TH1F *hist2 = new TH1F("hist2", "hist2", 50, 0, 10);
    TH1F *hist3 = new TH1F("hist3", "hist3", 50, 0, 10);
    TH1F *hist4 = new TH1F("hist4", "hist4", 50, 0, 10);
    TH1F *hist5 = new TH1F("hist5", "hist5", 50, 0, 10);

    {
        TRandom3 *random = new TRandom3(1);
        for(int i = 0; 10000 > i; ++i)
        {
            hist1->Fill(random->Gaus(5, 1));
            hist2->Fill(random->Poisson(5));
            hist3->Fill(random->Uniform(0, 10));
            hist4->Fill(random->Gaus(4, 1));
            hist5->Fill(random->Poisson(3));
        }
        delete random;
    }

    TFile *output = new TFile("output.root", "RECREATE");
    TDirectory *group1 = output->mkdir("group1");
    group1->cd();

    hist1->Write();
    hist2->Write();

    TDirectory *group2 = output->mkdir("group2");
    group2->cd();

    hist3->Write();
    hist4->Write();
    hist5->Write();

    delete output;

    application->Run(true);

    delete hist1;
    delete hist2;
    delete hist3;
    delete hist4;
    delete hist5;

    return 0;
}

// ------------------------------------------------------------------------------------------------------//
// ------------------------------------------------------------------------------------------------------//
void frac_error(double a, double b, double *ret){

	double c = (b-a);
	gStyle->SetOptFit(0011);
	gStyle->SetOptStat(0); 
	TCanvas *can = new TCanvas("can","can",900,600);

	double frac = a/b;
	TH1F *h_    = new TH1F("h_","h_",100,frac*0.98,frac*1.02);
	TRandom3 *r  = new TRandom3();
	
	for(int i=0;i<50000;i++){
	 
	  double ap = r->Poisson(a);
	  double cp = r->Poisson(c);
	  h_->Fill((double)ap/(ap+cp));
	  
	}
	
	h_->Fit("gaus","Q","Q");
		
	h_->Draw();
	can->SaveAs("statEP.png");
	double val = h_->GetFunction("gaus")->GetParameter(2);
	delete h_;
	
	ret[0] = val;
}

int main()
{
	// Manual input:
	if (true)
	{
		vector<KLV> gen(4); // pt eta phi m
		gen[1].p4.SetCoordinates(30,  1.0, 0.5, 0);
		gen[3].p4.SetCoordinates(25, -1.0, 2.5, 0);
		gen[2].p4.SetCoordinates(15,  2.1, 1.1, 0);
		gen[0].p4.SetCoordinates(10,  2.3, 1.0, 0);
		vector<KLV> rec(3);
		rec[2].p4.SetCoordinates(20,  1.3, 0.4, 0);
		rec[0].p4.SetCoordinates(18, -0.9, 2.2, 0);
		rec[1].p4.SetCoordinates(13,  2.2, 1.1, 0);
		test(gen, rec);
	}

	// Automatic input
	TRandom3 rnd;
	for (size_t i = 0; i < 4; ++i)
	{
		vector<KLV> gen, rec;
		// GEN stuff
		for (size_t j = 0; j < (size_t)rnd.Poisson(10); ++j)
		{
			gen.push_back(KLV());
			gen.back().p4.SetCoordinates(rnd.Uniform(5, 100), rnd.Gaus(0, 5), rnd.Uniform(-M_PI, +M_PI), rnd.Gaus(0, 10));
			// RECO efficiency
			if (rnd.Uniform(0, 1) > 0.05)
			{
				rec.push_back(KLV());
				rec.back().p4.SetCoordinates(
					gen.back().p4.pt() * fabs(rnd.Gaus(0.8, 0.1)),
					gen.back().p4.eta() * fabs(rnd.Gaus(1, 0.1)),
					gen.back().p4.phi() + rnd.Uniform(0, 0.1),
					gen.back().p4.mass() * rnd.Gaus(1, 0.1));
			}
		}
		// RECO noise
		for (size_t j = 0; j < (size_t)rnd.Poisson(5); ++j)
		{
			rec.push_back(KLV());
			rec.back().p4.SetCoordinates(rnd.Uniform(1, 10), rnd.Gaus(0, 5), rnd.Uniform(-M_PI, +M_PI), rnd.Gaus(0, 10));
		}
		test(gen, rec);
	}
}

Int_t TwoBody::CreateDimuonToyMc( void ){
  //
  // generate a toy di-muon dataset with systematics
  // set mData accordingly
  //

  // generate expected number of events from its uncertainty
  //RooDataSet * _ds = ws->pdf("syst_nbkg_dimuon")->generate(*ws->var("beta_nbkg_dimuon"), 1);
  //Double_t _ntoy = ((RooRealVar *)(_ds->get(0)->first()))->getVal() * (ws->var("nbkg_est_dimuon")->getVal());
  //delete _ds;

  Double_t _beta = GetRandom("syst_nbkg_dimuon", "beta_nbkg_dimuon");
  //  Double_t _kappa = ws->var("nbkg_kappa_dimuon")->getVal();
  Double_t _nbkg_est = ws->var("nbkg_est_dimuon")->getVal();
  //Double_t _ntoy = pow(_kappa,_beta) * _nbkg_est;
  Double_t _ntoy = _beta * _nbkg_est;
 
  Int_t _n = r.Poisson(_ntoy);
  // all nuisance parameters:
  //   beta_nsig_dimuon, 
  //   beta_nbkg_dimuon,
  //   lumi_nuis

  // create dataset
  RooRealVar * _mass = ws->var("mass");
  RooArgSet _vars(*_mass);

  RooAbsPdf * _pdf = ws->pdf("bkgpdf_dimuon");

  RooAbsPdf::GenSpec * _spec = _pdf->prepareMultiGen(_vars,
						     Name("toys"),
						     NumEvents(_n),
						     Extended(kFALSE),
						     Verbose(kTRUE));

  //RooPlot* xframe = _mass->frame(Title("Gaussian p.d.f.")) ;
  //realdata->plotOn(xframe,LineColor(kRed),MarkerColor(kRed));

  delete data;
  data = _pdf->generate(*_spec); // class member
  delete _spec;

  //data->plotOn(xframe);
  //TCanvas* c = new TCanvas("test","rf101_basics",800,400) ;
  //gPad->SetLeftMargin(0.15) ; xframe->GetYaxis()->SetTitleOffset(1.6) ; xframe->Draw() ;
  //c->SaveAs("test.pdf");

  Int_t n_generated_entries = (Int_t)(data->sumEntries());

  // debug
  std::cout << "!!!!!!!!!!!!!! _beta = " << _beta << std::endl;
  //std::cout << "!!!!!!!!!!!!!! _kappa = " << _kappa << std::endl;
  std::cout << "!!!!!!!!!!!!!! _nbkg_est = " << _nbkg_est << std::endl;
  std::cout << "!!!!!!!!!!!!!! _ntoy     = " << _ntoy << std::endl;
  std::cout << "!!!!!!!!!!!!!! _n        = " << _n    << std::endl;
  std::cout << "!!!!!!!!!!!!!! n_generated_entries = " << n_generated_entries    << std::endl;
  return n_generated_entries;
}

// Read in an event spectrum from a text file containing
// spectral data, and then generate some random events.
// This assumes that the input data file contains two columns,
// the first with (evenly spaced!) energy bins in units of eV, and the second
// with a number of events per kg per day per keV
void FillRecoilSpectrumFromFile(TH1F* recoilHisto, Double_t time, Double_t mass, const char* filename, Int_t seed)
{
	FILE* dataFile;
	dataFile = fopen(filename,"r");
	if (dataFile == NULL) {
		cout << "The file " << filename << " is not found!!! Quitting now." << endl;
		exit(1);
	}
	
	// Count the number of lines to set the array sizes
	fstream lineCountStream;
	lineCountStream.open(filename, fstream::in);
	Int_t lineCount = 0;
	while (lineCountStream.peek() != EOF) {
		lineCountStream.ignore(128, '\n');
		lineCount++;
	}
	lineCountStream.close();
	cout << "The file has " << lineCount << " lines." << endl;
	
	// Should add a section here to double check the formatting of the file
	
	// Warning message about units
	cout << "Reading data and filling histogram from " << filename << "..." << endl;
	cout << "I am assuming that the first column in the file is in units of MeV "
			 << "and that the second column is in units of 1/(kg day keV)!" << endl << endl;
	
	// Read the data
	const Int_t nLines = lineCount;
	Float_t energies[nLines];
	Float_t rates[nLines];
	for (Int_t thisLine = 0; thisLine < nLines; thisLine++)
		fscanf(dataFile, "%E %E", &energies[thisLine], &rates[thisLine]);
	fclose(dataFile);
	
	// Set up the source histogram
	TH1F SpectralDataHistogram("SpectralDataHistogram", "Spectrum of Events from File", nLines-1, energies);
	for (Int_t thisBin = 0; thisBin < (nLines - 1); thisBin++) {
		SpectralDataHistogram.SetBinContent(thisBin,rates[thisBin]);
	}
	
	// Fill the recoil histogram with the correct number
	// of events, randomly chosen from the source histogram.
	TRandom3* randGen = new TRandom3(seed);
	Float_t meanNEvents = time * mass * SpectralDataHistogram.ComputeIntegral();
	Int_t nEvents = randGen->Poisson(meanNEvents);
	cout << "nEvents = " << nEvents << endl;
	for (Int_t thisEvent; thisEvent<nEvents; thisEvent++)
		recoilHisto->Fill(SpectralDataHistogram.GetRandom());
}

// Generate the number of events for this run.
// This is just Poisson statistics and the cross section.
Int_t GenerateNumOfNuRecoils(Double_t time, Double_t detMass, Double_t distance, Double_t activity, Int_t nNeutrons, Int_t nProtons, TF1* RecoilSpectrum, Double_t SpectrumMin, Double_t SpectrumMax, UInt_t seed)
{
	// Set up the constants
	Double_t AvogadroConst = 6.022e23;
	Double_t NucleonMass = (nNeutrons * 939.565) + (nProtons * 938.272); // [MeV]
	Double_t molarMass = nNeutrons + nProtons;
	
	// Calculate mean number of events
	Double_t SpectrumWeightedCrossSection = RecoilSpectrum->Integral(SpectrumMin, SpectrumMax);
	Double_t MeanEvents = time * AvogadroConst * activity * detMass * SpectrumWeightedCrossSection / (4 * TMath::Pi() * molarMass * pow(distance,2));
	
	// Generate a random number of events from Poisson distribution
	TRandom3* randGen = new TRandom3(seed);
	Int_t nEvt = randGen->Poisson(MeanEvents);
	return nEvt;
}

void bins_chi2_n() {

//-------------------------------------------------
//Define here the luminosities for period 1 and 2
//-------------------------------------------------
    float l1 = 3.99;
    float l2 = 3.66;
    float ltot = l1+l2;


//-------------------------------------------------
//Define here the number of boxes (indepednant) and the yields
//-------------------------------------------------
//observed 27 - 36 | CR 88 - 54
///*
    int nbins = 2;
    float n1[] = {27,88};
    float n2[] = {36,54};
//*/

    //observed 31 - 13 | CR 59 - 43
    /*
    int nbins = 2;
    float n1[] = {31,59};
    float n2[] = {13,43};
    */
    float *ntot = new float[nbins];
    float *pred1 = new float[nbins];
    float *pred2 = new float[nbins];
    for(int i=0; i<nbins; i++) {
        ntot[i] = n1[i]+n2[i];
        pred1[i] = ntot[i]/ltot*l1;
        pred2[i] = ntot[i]/ltot*l2;
    }


    double prob = 1;
//max should be larger that the maximum yield
    int max = 200;
    for(int i=0; i<nbins; i++) {
        if(max<n1[i]) max = n1[i]*2;
        if(max<n2[i]) max = n2[i]*2;

    }
    double chi2=0;
    double like0 = 1;
    for(int i=0; i<nbins; i++) {
        TF1 f1("poiss1","TMath::PoissonI(x,[1])",0,max);
        TF1 f2("poiss1","TMath::PoissonI(x,[1])",0,max);
        f1.SetParameter(1,pred1[i]);
        f2.SetParameter(1,pred2[i]);
        if(f1.Integral(n1[i],max)<0.5)
            prob*=f1.Integral(n1[i],max);
        else
            prob*=f1.Integral(0,n1[i]);
        if(f2.Integral(n2[i],max)<0.5)
            prob*=f2.Integral(n2[i],max);
        else
            prob*=f2.Integral(0,n2[i]);

        //likelihood

        like0*=TMath::PoissonI(n1[i],pred1[i]);
        like0*=TMath::PoissonI(n2[i],pred2[i]);
        //chi2-test

        chi2+=pow((n1[i]-pred1[i]),2)/pred1[i];
        chi2+=pow((n2[i]-pred2[i]),2)/pred2[i];
        //prob *= f1.Integral(n1[i],max)*f2.Integral(n2[i],max);
    }
    cout<<"prob = "<<prob<<endl;
    cout<<"likelihood = "<<like0<<endl;
    cout<<"chi2 = "<<chi2<<" "<<TMath::Prob(chi2,1)<<endl;

//MC toys
    TRandom3 rand;
    int ntoys = 1000000;
    vector<double> vprob;
    for(int i=0; i<ntoys; i++) {
        double p = 1;
        double like = 1;
        for(int b=0; b<nbins; b++) {
            float v1 = rand.Poisson(pred1[b]);
            float v2 = rand.Poisson(pred2[b]);
            ///Compute prob
            like*=TMath::PoissonI(v1,pred1[b]);
            like*=TMath::PoissonI(v2,pred2[b]);

        }
        //vprob.push_back(prob);
        //cout<<like<<endl;
        vprob.push_back(like);
    }
//sort the vector
//std::sort(vprob.begin(),vprob.end(),myfunction);
//compute prob
    int count = 0;
    for(unsigned int x=0; x<vprob.size(); x++) {
        if(like0<vprob[x]) count++;
        // break;
//.........这里部分代码省略.........

void createTBevents(int input){

  printf("Starting Simulation of data\n");

  //creating the output file
  char outputFileName[100] = {"OutputFile.root"};

  printf("Creating output file: %s \n",outputFileName);
  TFile * outputFile = new TFile(outputFileName,"RECREATE");


  //Counter for event number
  unsigned int eventNr;

  //Counter for total number of hits
  unsigned int hitsTotal = 0;

  
  short int col, row, adc;
  short int ladder = 2;
  short int mod = 3;
  short int disk = 2;
  short int blade = 2;
  short int panel = 2;
  
  //create the tree to store the data
  TTree *bpixTree[3];

  char title[30];
  for (int i=1; i<4; i++){
    sprintf(title,"BPIX_Digis_Layer%1d",i);
    bpixTree[i-1]= new TTree(title,title);
    bpixTree[i-1]->Branch("Event", &eventNr, "Event/i");           
    bpixTree[i-1]->Branch("Ladder", &ladder, "Ladder/S");      
    bpixTree[i-1]->Branch("Module", &mod, "Module/S");      
    bpixTree[i-1]->Branch("adc", &adc, "adc/S");       
    bpixTree[i-1]->Branch("col", &col, "col/S");       
    bpixTree[i-1]->Branch("row", &row, "row/S");
  }

  

  //Maximum number of events. Events does not correspond with Hits
  unsigned int maxEventNr = input;

  //Number of Hits per Event
  //This should be randomized later and be dependant on the rate
  double meanHitsPerEvent = 2;

  //Maximum particle flux [MHz cm^-2]
  int maxParticleFlux = 500;

  //number of hits in current event
  int hitsInEvent = -1; 

  //create a random number generator
  TRandom3 * random = new TRandom3();
  TRandom3 * randomrow = new TRandom3();
  TRandom3 * randomcol = new TRandom3();
  TRandom3 * randomadc = new TRandom3();

  //using custom function to distribute values
  //values used from http://ntucms1.cern.ch/pixel_dev/flux/v8/016031/fitspot_bin_11.pdf
  
  TF1 * fx = new TF1("xfunc", "[0]*exp(2.59349*exp(-0.5*((-x+3.24273/.15+[1]/.15)/7.07486*.15)**2)+2.07765*exp(-0.5*((-x+9.33060e-01/.15+[1]/.15)/2.24067*.15)**2)-4.21847)",0, 52);
  
  fx->SetParameters(1,337.0);
  fx->SetParameters(2,1.74);
  
  TF1 * fy = new TF1("yfunc", AsyGaus,0,80,4);
  
  fy->SetParNames("mean","sigma1","sigma2","amplitude");
  fy->SetParameter("mean",43);
  fy->SetParameter("sigma1",11.4);
  fy->SetParameter("sigma2",15.0);
  fy->SetParameter("amplitude",347.0);
  
  TF1 * fadc = new TF1("adcfunc", langaufun,0,400,4);
  fadc->SetParNames("scale","mpv","area","sigma");
  fadc->SetParameter("scale",19);
  fadc->SetParameter("mvp",220);
  fadc->SetParameter("area",10000);
  fadc->SetParameter("sigma",30);




  while (eventNr < maxEventNr){
    //printf("eventNr: %d \n",eventNr);
    //Function used for fitting according to Xin an Stefano
    
    //Start by generating the number of hits per event
    //following a poisson distribution
    random->SetSeed(0);
    hitsInEvent = random->Poisson(meanHitsPerEvent);
    //printf("hitsInEvent %d \n", hitsInEvent);

    hitsTotal += hitsInEvent;

    if(hitsInEvent < 0){
//.........这里部分代码省略.........

//-------------------------------------------------------------
//Main macro for generating data and fitting
//=============================================================  
void FitMassPhotonResolutionSystematics(const string workspaceFile = "/afs/cern.ch/work/d/daan/public/releases/CMSSW_5_3_9_patch3/src/CMSAna/HHToBBGG/data/FitWorkspace_asdf.root", const string outputTree = "/afs/cern.ch/work/d/daan/public/releases/CMSSW_5_3_9_patch3/src/CMSAna/HHToBBGG/data/MassFitResults/MassFitTwoD_asdf.root", Int_t plotOption = 0, Int_t storeOption = 1, Int_t SystematicsUpOrDown = 1, Int_t NToys = 5000) {

  TRandom3 *randomnumber = new TRandom3(1200);
  TFile *wsFile = new TFile (workspaceFile.c_str(), "READ");
  RooWorkspace *ws = (RooWorkspace*)wsFile->Get("MassFitWorkspace");
  
  //Import variables from workspace
  RooAbsPdf *model2Dpdf = ws->pdf("model2Dpdf");
  RooRealVar *massBjet = ws->var("massBjet");
  RooRealVar *massPho = ws->var("massPho");
  RooRealVar *nsig = ws->var("N (Sig)"); RooRealVar constNsig(*nsig);
  RooRealVar *nres = ws->var("N (ResBkg)"); RooRealVar constNres(*nres);
  RooRealVar *nnonres = ws->var("N (NonResBkg)"); RooRealVar constNnonres(*nnonres);
  RooRealVar *expRateBjet = ws->var("expRateBjet"); RooRealVar constexpBjet(*expRateBjet);
  RooRealVar *expRatePho = ws->var("expRatePho"); RooRealVar constexpPho(*expRatePho);

  //Variables to set constant
  RooRealVar *sigMeanBjet = ws->var("sigMeanBjet"); sigMeanBjet->setConstant();
  RooRealVar *sigSigmaBjet = ws->var("sigSigmaBjet"); sigSigmaBjet->setConstant();
  RooRealVar *sigAlpha = ws->var("sigAlpha"); sigAlpha->setConstant();
  RooRealVar *sigPower = ws->var("sigPower"); sigPower->setConstant();
  RooRealVar *sigmaPho = ws->var("sigmaPho"); sigmaPho->setConstant();

  RooRealVar *resMeanBjet = ws->var("resMeanBjet"); resMeanBjet->setConstant();
  RooRealVar *resSigmaBjet = ws->var("resSigmaBjet"); resSigmaBjet->setConstant();
  RooRealVar *resAlpha = ws->var("resAlpha"); resAlpha->setConstant();
  RooRealVar *resPower = ws->var("resPower"); resPower->setConstant();
  RooRealVar *resExpo = ws->var("resExpo"); resExpo->setConstant();
  RooRealVar *nbbH = ws->var("nbbH"); nbbH->setConstant();
  RooRealVar *nOthers = ws->var("nOthers"); nOthers->setConstant();
  
  double inputSigmaPho = sigmaPho->getVal();

  //Create TTree to store the resulting yield data
  TFile *f = new TFile(outputTree.c_str(), "RECREATE");
  TTree *resultTree = new TTree("resultTree", "Parameter results from fitting");
  Float_t nsigOut, nresOut, nnonresOut;
  Float_t nsigStd, nresStd, nnonresStd;
  
  resultTree->Branch("nsigOut",&nsigOut, "nsigOut/F");
  resultTree->Branch("nresOut",&nresOut, "nresOut/F");
  resultTree->Branch("nnonresOut",&nnonresOut, "nnonresOut/F");
  resultTree->Branch("nsigStd",&nsigStd, "nsigStd/F");
  resultTree->Branch("nresStd",&nresStd, "nresStd/F");
  resultTree->Branch("nnonresStd",&nnonresStd, "nnonresStd/F");
  
  //Generate Toy MC experiment data and fits
  for(UInt_t t=0; t < NToys; ++t) {
    cout << "Toy #" << t << endl;
    nsig->setVal(constNsig.getVal()); nres->setVal(constNres.getVal()); nnonres->setVal(constNnonres.getVal());
    expRateBjet->setVal(constexpBjet.getVal()); expRatePho->setVal(constexpPho.getVal());
   
    //set jet energy resolutions to nominal
    sigmaPho->setVal(inputSigmaPho);

    cout << "Before: " << sigmaPho->getVal() << " | ";

    Float_t ran = randomnumber->Poisson(325.);
    RooDataSet *pseudoData2D = model2Dpdf->generate(RooArgList(*massBjet,*massPho), ran);

    //move jet energy resolution up/down
    if (SystematicsUpOrDown == 1) {
      sigmaPho->setVal(inputSigmaPho*1.15);
    } else if (SystematicsUpOrDown == -1) {
      sigmaPho->setVal(inputSigmaPho/1.15);
    }

    cout << "After: " << sigmaPho->getVal() << " \n";

    RooFitResult *fitResult2D = model2Dpdf->fitTo(*pseudoData2D, RooFit::Save(), RooFit::Extended(kTRUE), RooFit::Strategy(2));
//     if (t == 1763) {
//     	ws->import(*pseudoData2D, kTRUE);
//     	ws->import(*pseudoData2D, Rename("pseudoData2D"));
//     }
//     if (plotOption == 1) MakePlots(ws, fitResult2D);
    

    cout << "DEBUG: " << constexpBjet.getVal() << " , " << constexpPho.getVal() << " | " << expRateBjet->getVal() << " " << expRatePho->getVal() << "\n";


    //Store fit parameters into ROOT file
    if (storeOption == 1) {
      //Save variables into separate branches of root tree
  		nsigOut = nsig->getVal();
  		nresOut = nres->getVal();
  		nnonresOut = nnonres->getVal();
  		nsigStd = nsig->getPropagatedError(*fitResult2D);
  		nresStd = nres->getPropagatedError(*fitResult2D);
  		nnonresStd = nnonres->getPropagatedError(*fitResult2D);
  		//cout << "\n\n\n\n\n\n" << nsigOut << " | " << nresOut << " | " << nnonresOut << " | " << ran  << "\n\n\n\n\n" << endl;
  		resultTree->Fill();
    }
  }
  
  //Write to the TTree and close it
  resultTree->Write();
  f->Close();
//.........这里部分代码省略.........

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

  // get the same-sign distributions from bkg
  plotContainer SS_tightIso_bkg_plots ("SS_tightIso_bkg", variablesList, variables2DList, selections_SS_tightIso, bkgSamplesList, 0) ;
  counters SS_tightIso_bkgCount = fillHistos (bkgSamples, SS_tightIso_bkg_plots, 
              variablesList, variables2DList,
              selections_SS_tightIso,
              lumi,
              vector<float> (0),
              false, false, maxEvtsMC) ;
  SS_tightIso_bkg_plots.AddOverAndUnderFlow () ;


  plotContainer SS_QCD_tightIso ("SS_tightIso_QCD", variablesList, variables2DList, selections_SS_tightIso, QCDsample, 0) ;
//  vector <TH1F *> SS_QCD ;
  vector<vector<float>> QCDyieldSSregionTightIso (variablesList.size(), vector<float>(selections_SS_tightIso.size()) ); // var, cut


  cout << "--- MAIN reading DATA and filling SS histos with non-relaxed iso" << endl ;

  // get the same-sign distributions from data
  //plotContainer SS_tightIso_DATA_array[nruns];
  TRandom3 *g = new TRandom3();

    for (unsigned int icut = 0 ; icut < selections_SS_tightIso.size () ; ++icut)
    {
      TH1F *h =new TH1F(selections_SS_tightIso.at (icut).first.Data (),selections_SS_tightIso.at (icut).first.Data (),10000,1,10001);
      QCDYields.push_back(h);   
    }

    plotContainer SS_tightIso_DATA_plots ("SS_tightIso_DATA", variablesList, variables2DList, selections_SS_tightIso, DATASamplesList, 2) ;
    counters SS_tightIso_DATACount = fillHistos (DATASamples, SS_tightIso_DATA_plots, 
      variablesList, variables2DList,
      selections_SS_tightIso,
      lumi,
      vector<float> (0),
      true, false) ;
    SS_tightIso_DATA_plots.AddOverAndUnderFlow () ;  
    //SS_tightIso_DATA_array[irun] = SS_tightIso_DATA_plots;
  //}

  //cout << "--- MAIN preparing to loop on variables and selections to calc SS QCD yield with non-relaxed iso" << endl ;

  // the index in the stack is based on variable ID (iv) and selection ID (isel):
  // iHisto = iv + nVars * isel
  
  //vector<string> QCDsample ;
  //QCDsample.push_back ("QCD") ;
 
    for (unsigned int icut = 0 ; icut < selections_SS_tightIso.size () ; ++icut)
    {
      for (unsigned int ivar = 0 ; ivar < variablesList.size () ; ++ivar)
      {
        //plotContainer SS_tightIso_DATA_plots = SS_tightIso_DATA_array[irun]
        THStack * D_stack = SS_tightIso_DATA_plots.makeStack (variablesList.at (ivar),
          selections_SS_tightIso.at (icut).first.Data ()) ;
        TH1F * tempo = (TH1F *) D_stack->GetStack ()->Last () ;
        TH1F * orig = (TH1F*)tempo->Clone("original");
        THStack * b_stack = SS_tightIso_bkg_plots.makeStack (variablesList.at (ivar),
          selections_SS_tightIso.at (icut).first.Data ()) ;
        TH1F * h_bkg = (TH1F *) b_stack->GetStack ()->Last () ;
        int nbins = tempo->GetNbinsX();

        for(int irun =0;irun<nruns;irun++){
          //TString name; name.Form("%s%d",tempo->GetName (),irun) ;
          //tempo->SetName(name);
          //tempo->SetBinErrorOption(TH1::kPoisson);
          for(int ibin=1; ibin<=nbins;ibin++){
            tempo->SetBinContent(ibin,g->Poisson(orig->GetBinContent(ibin)));
            //cout<<tempo->GetBinContent(ibin)<<" "<<orig->GetBinContent(ibin)<<endl;
          }
        
        //name = TString ("DDQCD_tightIso_") + name ;
        //TH1F * dummy = (TH1F *) tempo->Clone (name) ;
        tempo->Add (h_bkg, -1) ;
        //cout<<"integral "<<tempo->Integral()<<endl;
        QCDYields.at(icut)->Fill(tempo->Integral());
      }
    }
  }

  //}
  // get the SS-to-OS scale factor and scale the QCD distributions
  // ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----


  // Save the histograms
  // ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----

  cout << "--- MAIN before saving" << endl ;
  TString outString ;
  outString.Form ("outPlotter.root") ;
  TFile * fOut = new TFile (outString.Data (), "RECREATE") ;

  fOut->cd();
  for(uint i=0; i<QCDYields.size();i++)QCDYields.at(i)->Write();
  fOut->Close () ;

  
  return 0 ;
}  

//.........这里部分代码省略.........
  TCanvas *cv = new TCanvas("cv","cv",800,600);

  RooPlot *mframeFail_default = m_default->frame(Bins(Int_t(130-60)/2));
  modelFail_default->plotOn(mframeFail_default);
  modelFail_default->plotOn(mframeFail_default,Components("ebackgroundFail"),LineStyle(kDashed),LineColor(kRed));
  modelFail_default->plotOn(mframeFail_default,Components("bkgexpFail"),LineStyle(kDashed),LineColor(kGreen+2));
  mframeFail_default->GetYaxis()->SetTitle("");
  mframeFail_default->GetYaxis()->SetTitleOffset(1.2);
  mframeFail_default->GetXaxis()->SetTitle("m_{ee#gamma} [GeV/c^{2}]");
  mframeFail_default->GetXaxis()->SetTitleOffset(1.05);
  mframeFail_default->SetTitle("");
  mframeFail_default->Draw();

  cv->SaveAs("DefaultModel.gif");

  RooPlot *mframeFail = m_default->frame(Bins(Int_t(130-60)/2));
  modelFail->plotOn(mframeFail);
  modelFail->plotOn(mframeFail,Components("ebackgroundFail_shifted"),LineStyle(kDashed),LineColor(kRed));
  modelFail->plotOn(mframeFail,Components("bkgexpFail_shifted"),LineStyle(kDashed),LineColor(kGreen+2));
  mframeFail->GetYaxis()->SetTitle("");
  mframeFail->GetYaxis()->SetTitleOffset(1.2);
  mframeFail->GetXaxis()->SetTitle("m_{ee#gamma} [GeV/c^{2}]");
  mframeFail->GetXaxis()->SetTitleOffset(1.05);
  mframeFail->SetTitle("");
  mframeFail->Draw();
  cv->SaveAs(Form("ShiftedModel_%d.gif",Option));


  //*************************************************************************************
  //Do Toys
  //*************************************************************************************
  for(uint t=0; t < NToys; ++t) {

    RooDataSet *pseudoData_pass    = modelPass_default->generate(*m_default, randomnumber->Poisson(npass));
    RooDataSet *pseudoData_fail  = 0;
    pseudoData_fail    = modelFail->generate(*m_default, randomnumber->Poisson(nfail));
    RooDataSet *pseudoDataCombined = new RooDataSet("pseudoDataCombined","pseudoDataCombined",RooArgList(*m_default),
                                                    RooFit::Index(sample),
                                                    RooFit::Import("Pass",*pseudoData_pass),
                                                    RooFit::Import("Fail",*pseudoData_fail));

    pseudoDataCombined->write(Form("toy%d.txt",t));

    RooFitResult *fitResult=0;
    fitResult = totalPdf->fitTo(*pseudoDataCombined,
                                RooFit::Extended(),
                                RooFit::Strategy(2),
                                //RooFit::Minos(RooArgSet(eff)),
                                RooFit::Save());

    cout << "\n\n";
    cout << "Eff Fit: " << eff->getVal() << " -" << fabs(eff->getErrorLo()) << " +" << eff->getErrorHi() << endl;

    //Fill Tree
    varEff = eff->getVal();
    varEffErrL = fabs(eff->getErrorLo());
    varEffErrH = eff->getErrorHi();
    outTree->Fill();


//   //*************************************************************************************
//   //Plot Toys
//   //*************************************************************************************
//   TCanvas *cv = new TCanvas("cv","cv",800,600);
//   char pname[50];
//   char binlabelx[100];

TestProblem
AtlasDiJetMass(const int Nt,
               const int Nr,
               double tbins[],
               double rbins[],
               const double apar,
               const double bpar,
               const int nEvents,
               const double evtWeight)
{
  // From "Fully Bayesian Unfolding" by G. Choudalakis.
  // see arXiv:1201.4612v4
  // Recommended binning:
  // double bins[Nt+1] = {0};
  // for (int j=0; j<=Nt; j++)
  //   bins[j] = 500*TMath::Exp(0.15*j);

  static int id = 0; id++;
  TestProblem t;
  TRandom3 ran;

  // Mass distribution dN/dM
  TF1 *mass = new TF1("mass_dist",
                      "TMath::Power(1.-x/7000,6.0)/TMath::Power(x/7000,4.8)",
                      tbins[0], tbins[Nt]);

  // Generate test problem by MC convolution
  TH1D *hT   = new TH1D("hT",   "Truth mass dist. #hat{T}", Nt, tbins);
  TH1D *hTmc = new TH1D("hTmc", "MC Truth mass dist. #tilde{T}", Nt, tbins);
  TH1D *hD   = new TH1D("hD", "Measured mass dist.", Nr, rbins);
  TH2D *hM   = new TH2D("hM", "Migration matrix", Nr, rbins, Nt, tbins);
  hM->SetTitle(Form("%d x %d migration matrix M_{tr};"
                    "mass (observed);mass (true)", Nr, Nt));
  hT->SetLineWidth(2);
  hD->SetLineWidth(2);

  std::cout << Form("Generating test problem...") << std::flush;
  // Fill histos with MC events.
  // The response matrix gets more statistics than the data.
  double xt, xm, sigma;
  for (int i=0; i<nEvents; i++)
  {
    xt = mass->GetRandom();
    sigma = apar*TMath::Sqrt(xt) + bpar*xt;
    xm = xt + ran.Gaus(0, sigma);
    hM->Fill(xm, xt);
    hTmc->Fill(xt, evtWeight);
    hD->Fill(xm, evtWeight);
  }

  // Simulate Poisson fluctuations in real data (integer content, empty bins)
  for (int r=1; r<=Nr; r++)
    hD->SetBinContent(r, ran.Poisson(hD->GetBinContent(r)));

  // The true truth \hat{T}
  double totmass = mass->Integral(tbins[0],tbins[Nt]);
  for (int j=1; j<=Nt; j++)
  {
    hT->SetBinContent(j, evtWeight*nEvents*mass->Integral(tbins[j-1],
                      tbins[j])/totmass);
    hT->SetBinError(j, 0);
  }
  cout << "Done." << endl;

  // Projection of migration matrix to truth axis. hMt is the
  // numerator for efficiency. Bin contents should be counts
  // here. Elements are normalized to contain probabilities only after
  // projection & division by T-tilde.
  TH1D *hMt  = hM->ProjectionY("hMt",1,Nr);
  TH1D *heff = (TH1D *)hMt->Clone("heff");
  heff->Divide(hTmc);
  heff->Scale(evtWeight);
  hM->Scale(1./nEvents);
  hMt->Scale(1./nEvents);

  t.Response  = hM;
  t.xTruth    = hT;
  t.xTruthEst = hTmc;
  t.xIni      = hMt;
  t.bNoisy    = hD;
  t.eff       = heff;

  return t;
}

//.........这里部分代码省略.........
  
  // count events
  if( mode == "exclusion" )
  {
    RooAbsReal* rooTotIntegral = rooTotPdf -> createIntegral(x,NormSet(x),Range("signal"));
    NBKGToy_signal_fit = rooTotIntegral->getVal() * rooNBKGTot->getVal();
    
    std::cout << ">>>>>> BKG toy events (true) in signal region in " << lumi << "/pb: " << NBKGToy_signal << std::endl;  
    std::cout << ">>>>>> BKG toy events (fit)  in signal region in " << lumi << "/pb: " << NBKGToy_signal_fit << std::endl;      
  }
  
  if( mode == "discovery" )
  {
    std::cout << ">>>>>> BKG toy events (true) in " << lumi << "/pb: " << NBKGToy << std::endl;  
    std::cout << ">>>>>> BKG toy events (fit)  in " << lumi << "/pb: " << rooNBKGTot->getVal() << std::endl;  
    std::cout << ">>>>>> SIG toy events (true) in " << lumi << "/pb: " << NSIGToy << std::endl;  
    std::cout << ">>>>>> SIG toy events (fit)  in " << lumi << "/pb: " << rooNSIG->getVal() << std::endl;  
  }
  
  if( drawPlots )
  {
    TCanvas* c3 = new TCanvas("TOY","TOY");
    c3 -> SetGridx();
    c3 -> SetGridy();
    
    RooPlot* rooTOYPlot = x.frame();
    rooBKGToyDataSet -> plotOn(rooTOYPlot,MarkerSize(0.7));
    rooTotPdf -> plotOn(rooTOYPlot, LineColor(kRed));
    rooTotPdf -> plotOn(rooTOYPlot, Components("rooSIGPdf"), LineColor(kRed));
    rooTOYPlot->Draw();
    
    char pngFileName[50];
    sprintf(pngFileName,"BKGToyFit_H%d_%s.png",mH,mode.c_str());
    c3 -> Print(pngFileName,"png");
  }
  
  
  
  
  
  
  //-------------------------
  // generate toy experiments
  
  TH1F* h_BKGRes = new TH1F("h_BKGRes","",200,-400,400);
  TH1F* h_SIGRes = new TH1F("h_SIGRes","",200,-400,400);
  
  TRandom3 B;
  TRandom3 S;
  for(int j = 0; j < nToys; ++j)
  {
    if( j%100 == 0 )
      std::cout << ">>>>>> generating toy experiment " << j << std::endl;
    
    NBKGToy = B.Poisson(BKGTotShapeHisto->Integral());
    NSIGToy = 0;
    if( mode == "discovery" ) NSIGToy = S.Poisson(SIGShapeHisto->Integral());
    
    rooBKGToyDataSet = rooBKGTotPdf->generate(RooArgSet(x),NBKGToy);
    rooSIGToyDataSet = rooSIGPdf->generate(RooArgSet(x),NSIGToy);
    rooBKGToyDataSet -> append(*rooSIGToyDataSet);
    
    NBKGToy_signal = rooBKGToyDataSet->sumEntries(signalCut);
    NBKGToy_signal_fit = 0.;
    
    
    
    // fit
    if( mode == "exclusion" ) rooTotPdf -> fitTo(*rooBKGToyDataSet,Extended(kTRUE),PrintLevel(-1),Range("low,high"));
    if( mode == "discovery" ) rooTotPdf -> fitTo(*rooBKGToyDataSet,Extended(kTRUE),PrintLevel(-1));    
    
    
    
    // count events
    if( mode == "exclusion" )
    {
      RooAbsReal* rooTotIntegral = rooTotPdf -> createIntegral(x,NormSet(x),Range("signal"));
      NBKGToy_signal_fit = rooTotIntegral->getVal() * rooNBKGTot->getVal();
      
      h_BKGRes -> Fill(NBKGToy_signal_fit - NBKGToy_signal);
      h_SIGRes -> Fill(0.);
    }
    
    if( mode == "discovery" )
    {
      h_BKGRes -> Fill(rooNBKGTot->getVal() - NBKGToy);
      h_SIGRes -> Fill(rooNSIG->getVal() - NSIGToy);
    }
    
  }
  
  
  
  outFile -> cd();
  
  h_BKGRes -> Write();
  h_SIGRes -> Write();
  
  outFile -> Close();
}