C++ TGraph2D::SetPoint方法代码示例

void fillGraphsFromFilesDeltaNLL( const TString& par1name,
				  const TString& par2name,
				  const vector<TString>& fnames,
				  vector<string>&  keys,
				  map<string,TGraph2D *>& m_graphs)
{

  std::cout << "fillGraphsFromFilesDeltaNLL 1" << std::endl;

  keys.push_back("exp68");
  keys.push_back("exp95");
  keys.push_back("exp99");

  // uncommented below to plot observed!
  keys.push_back("obs95");

  TGraph2D *grobs = new TGraph2D();
  TGraph2D *grexp = new TGraph2D();

  m_graphs["obs95"] = grobs;
  m_graphs["exp95"] = grexp;

  grobs->SetName("graph2Dobs95");
  grexp->SetName("graph2Dexp95");

  Int_t nobs=0, nexp=0;

  for( size_t i=0; i<fnames.size(); i++) {
    
    TFile *f = new TFile(fnames[i]);
    TTree *t = (TTree *) f->Get("limit");

    if (!t) { 
      std::cerr<<"TFile "<<f->GetName()<<" does not contain the tree"<<std::endl;
      return;
    }
    cout << fnames[i] << " has limit tree with " << t->GetEntries() << " entries." << endl;

    Float_t deltaNLL, par1, par2;
    Int_t iToy;

    t->SetBranchAddress("iToy", &iToy);
    t->SetBranchAddress("deltaNLL", &deltaNLL);
    t->SetBranchAddress(par1name, &par1);
    t->SetBranchAddress(par2name, &par2);

    for (size_t j = 0, n = t->GetEntries(); j < n; ++j) {
      t->GetEntry(j);
      printf ("%d\r",j);
      if( !iToy){
	//	cout <<"!iToy" << endl;
	grobs->SetPoint(nobs++,par1,par2,2*deltaNLL);
	//	cout <<"grobs->SetPoint("<<nobs++<<","<<par1<<","<< par2<< ","<< 2*deltaNLL << endl;
      }
      else if (iToy == -1) {
	//	cout <<"iToy == -1" << endl;
	grexp->SetPoint(nexp++,par1,par2,2*deltaNLL);
      }
      else {
	cerr << "Unexpected value for iToy, = " << iToy << endl;
	exit(-1);
      }
    } // tree entry loop

    f->Close();
    delete f;

  } // file loop
  cout << endl;

  m_graphs["exp68"] = (TGraph2D*)grexp->Clone("graph2Dexp68");
  m_graphs["exp99"] = (TGraph2D*)grexp->Clone("graph2Dexp99");

#if 0
  TCanvas *canv = new TCanvas("tester","tester",500,500);
  cout << grexp->GetN()<<" points. " <<endl;
  grexp->Draw("TRI"); // cont 5z list");
#endif
}                                         // fillGraphsFromFilesDeltaNLL

void view_SMEvents_3D_from_Hits() {
	/*** Displays an 3D occupancy plot for each SM Event. (stop mode event)

	Can choose which SM event to start at. (find "CHOOSE THIS" in this script)
	Input file must be a Hits file (_interpreted_Hits.root file).
	***/
	gROOT->Reset();

	// Setting up file, treereader, histogram
	TFile *f = new TFile("/home/pixel/pybar/tags/2.0.2_new/pyBAR-master/pybar/module_202_new/101_module_202_new_stop_mode_ext_trigger_scan_interpreted_Hits.root");


	if (!f) { // if we cannot open the file, print an error message and return immediately
		cout << "Error: cannot open the root file!\n";
		//return;
	}

	TTreeReader *reader = new TTreeReader("Table", f);

	TTreeReaderValue<UInt_t> h5_file_num(*reader, "h5_file_num");
	TTreeReaderValue<Long64_t> event_number(*reader, "event_number");
	TTreeReaderValue<UChar_t> tot(*reader, "tot");
	TTreeReaderValue<UChar_t> relative_BCID(*reader, "relative_BCID");
	TTreeReaderValue<Long64_t> SM_event_num(*reader, "SM_event_num");
	TTreeReaderValue<Double_t> x(*reader, "x");
	TTreeReaderValue<Double_t> y(*reader, "y");
	TTreeReaderValue<Double_t> z(*reader, "z");

	// Initialize the canvas and graph
	TCanvas *c1 = new TCanvas("c1","3D Occupancy for Specified SM Event", 1000, 10, 900, 550);
	c1->SetRightMargin(0.25);
	TGraph2D *graph = new TGraph2D();

	// Variables used to loop the main loop
	bool endOfReader = false; // if reached end of the reader
	bool quit = false; // if pressed q
	int smEventNum = 1; // the current SM-event CHOOSE THIS to start at desired SM event number
	
	// Main Loop (loops for every smEventNum)
	while (!endOfReader && !quit) {
		// Variables used in this main loop
		int startEntryNum = 0;
		int endEntryNum = 0;
		string histTitle = "3D Occupancy for SM Event ";
		string inString = "";
		bool fitFailed = false; // true if the 3D fit failed
		bool lastEvent = false;

		// Declaring some important output values for the current graph and/or line fit
		int numEntries = 0;
		double sumSquares = 0;

		// Get startEntryNum and endEntryNum
		startEntryNum = getEntryNumWithSMEventNum(reader, smEventNum);
		endEntryNum = getEntryNumWithSMEventNum(reader, smEventNum + 1);

		if (startEntryNum == -2) { // can't find the smEventNum
			cout << "Error: There should not be any SM event numbers that are missing." << "\n";
		} else if (startEntryNum == -3) { 
			endOfReader = true;
			break;
		} else if (endEntryNum == -3) { // assuming no SM event nums are skipped
			endEntryNum = reader->GetEntries(false);
			lastEvent = true;
		}

		// Fill TGraph with points and set title and axes
		graph = new TGraph2D(); // create a new TGraph to refresh

		reader->SetEntry(startEntryNum);
		for (int i = 0; i < endEntryNum - startEntryNum; i++) {
			graph->SetPoint(i, (*x - 0.001), (*y + 0.001), (*z - 0.001));
			endOfReader = !(reader->Next());
		}

		histTitle.append(to_string(smEventNum));
		graph->SetTitle(histTitle.c_str());
		graph->GetXaxis()->SetTitle("x (mm)");
		graph->GetYaxis()->SetTitle("y (mm)");
		graph->GetZaxis()->SetTitle("z (mm)");

		graph->GetXaxis()->SetLimits(0, 20); // ROOT is buggy, x and y use setlimits()
		graph->GetYaxis()->SetLimits(-16.8, 0); // but z uses setrangeuser()
		graph->GetZaxis()->SetRangeUser(0, 40.96);
		c1->SetTitle(histTitle.c_str());

		// 3D Fit, display results, draw graph and line fit, only accept "good" events, get input
		if (!endOfReader || lastEvent) {
			// Display some results
			numEntries = graph->GetN();
			cout << "Current SM Event Number: " << smEventNum << "\n";
			cout << "Number of entries:       " << numEntries << "\n";

			// Starting the fit. First, get decent starting parameters for the fit - do two 2D fits (one for x vs z, one for y vs z)
			TGraph *graphZX = new TGraph();
			TGraph *graphZY = new TGraph();
			reader->SetEntry(startEntryNum);
			for (int i = 0; i < endEntryNum - startEntryNum; i++) {
				graphZX->SetPoint(i, (*z - 0.001), (*x + 0.001));
				graphZY->SetPoint(i, (*z - 0.001), (*y + 0.001));
//.........这里部分代码省略.........

void makeLikelihood(){

   gSystem->Load("libHiggsAnalysisCombinedLimit.so");
   //TFile *fi = TFile::Open("lduscan_neg_ext/3D/lduscan_neg_ext_3D.root");
   TFile *fi = TFile::Open("allthepoints.root");
   TTree *tree = (TTree*)fi->Get("limit");
   //TTree *tree = new TTree("tree_vals","tree_vals");  

   // ------------------------------ THIS IS WHERE WE BUILD THE SPLINE ------------------------ //
   // Create 2 Real-vars, one for each of the parameters of the spline 
   // The variables MUST be named the same as the corresponding branches in the tree
   RooRealVar ldu("lambda_du","lambda_du",0.1,-2.5,2.5); 
   RooRealVar lVu("lambda_Vu","lambda_Vu",0.1,0,3);
   RooRealVar kuu("kappa_uu","kappa_uu",0.1,0,3);
   
   RooSplineND *spline = new RooSplineND("spline","spline",RooArgList(ldu,lVu,kuu),tree,"deltaNLL",1.,true,"deltaNLL<1000 && TMath::Abs(quantileExpected)!=1 && TMath::Abs(quantileExpected)!=0");
   // ----------------------------------------------------------------------------------------- //
   
   //TGraph *gr = spline->getGraph("x",0.1); // Return 1D graph. Will be a slice of the spline for fixed y generated at steps of 0.1
   fOut = new TFile("outplots-2hdm-neg-fine-mssm-final-try2.root","RECREATE");
   // Plot the 2D spline 
   TGraph2D *gr  		= new TGraph2D();  gr->SetName("type1");
   TGraph2D *gr2 		= new TGraph2D();  gr2->SetName("type2");
   TGraph2D *gr_ldu 		= new TGraph2D(); gr_ldu->SetName("ldu");
   TGraph2D *gr_lVu 		= new TGraph2D(); gr_lVu->SetName("lVu");
   TGraph2D *gr_kuu 		= new TGraph2D(); gr_kuu->SetName("kuu");
   TGraph2D *gr2_ldu		= new TGraph2D(); gr2_ldu->SetName("ldu_2");
   TGraph2D *gr2_lVu 		= new TGraph2D(); gr2_lVu->SetName("lVu_2");
   TGraph2D *gr2_kuu 		= new TGraph2D(); gr2_kuu->SetName("kuu_2");
   TGraph2D *gr_t1_lVu_V_kuu    = new TGraph2D(); gr_t1_lVu_V_kuu->SetName("t1_lVu_V_kuu");

   TGraph2D *gr_mssm_ldu 	= new TGraph2D(); gr_mssm_ldu->SetName("mssm_ldu");
   TGraph2D *gr_mssm_lVu 	= new TGraph2D(); gr_mssm_lVu->SetName("mssm_lVu");
   TGraph2D *gr_mssm_kuu 	= new TGraph2D(); gr_mssm_kuu->SetName("mssm_kuu");

   TGraph2D *g_FFS = new TGraph2D(); g_FFS->SetName("ffs_kuu1");
   double x,y,z;
   double mintF = 10000;
   int pt=0 ;
   /*
   for (double x=-1.6;x<=1.6;x+=0.05){
     for (double y=0.5;y<=1.5;y+=0.05){
	ldu.setVal(x);
	lVu.setVal(y);
	kuu.setVal(1);
	double dnll2 = 2*spline->getVal();
	if (dnll2 < mintF) mintF = dnll2;
	g_FFS->SetPoint(pt,x,y,dnll2);
	pt++;
     }
   }
   */

   TGraph2D *gcvcf = new TGraph2D(); gcvcf->SetName("cvcf");
   TGraph2D *gcvcf_kuu = new TGraph2D(); gcvcf_kuu->SetName("cvcf_kuu");
   TGraph2D *gcvcf_lVu = new TGraph2D(); gcvcf_lVu->SetName("cvcf_lVu");
   double mintkvkf = 10000;
   int pt=0 ;
   if (doXCHECK){
   // Sanity check, for ldu = 1, we should resolve kv kf ?
   //
   for (double cv=0.5;cv<=1.4;cv+=0.05){
     for (double cf=0.3;cf<=1.7;cf+=0.05){
	ldu.setVal(1.);
	lVu.setVal(cv/cf);
	kuu.setVal(cf*cf/gamma(cv,cf,cf));
	double dnll2 = 2*spline->getVal();
	if (dnll2 < mintkvkf) mintkvkf = dnll2;
	gcvcf->SetPoint(pt,cv,cf,dnll2);
	gcvcf_lVu->SetPoint(pt,cv,cf,lVu.getVal());
	gcvcf_kuu->SetPoint(pt,cv,cf,kuu.getVal());
	pt++;
     }
   }
   std::cout << " Min cV-cF = " << mintkvkf << std::endl;
   for (int p=0;p<gcvcf->GetN();p++){
        double z = (gcvcf->GetZ())[p] - mintkvkf;
        double x = (gcvcf->GetX())[p];
        double y = (gcvcf->GetY())[p];
	gcvcf->SetPoint(p,x,y,z);
   }
   }


   double Vldu, VlVu, Vkuu;

   int pt = 0;
   double mint2 = 10000;
   double mint1 = 10000;

   if (doTHDM){
   for (double scbma=-1;scbma<1;scbma+=0.01){
     for (double b=0.01;b<1.45;b+=0.01){
        double tanb = TMath::Tan(b);
     	if (tanb>1. ) b+=0.05;
	double cbma;
	if (scbma < 0) cbma = -1*scbma*scbma;
	else cbma = scbma*scbma;
	// Type 1 
	type1(cbma, tanb, &Vldu, &VlVu, &Vkuu);
//.........这里部分代码省略.........

void plotContsSingle(TFile *fOUT, std::string dirname, std::string fin, float X, int keepMDM=10){

   //gSystem->Load("libHiggsAnalysisCombinedLimit.so");
   gROOT->SetBatch(1);

   TFile *fiSignals = TFile::Open("signalsVA.root");
   RooWorkspace *workspace = (RooWorkspace*)fiSignals->Get("combinedws");
   TFile *fiSignalsPS = TFile::Open("signalsPS.root");
   RooWorkspace *workspacePS = (RooWorkspace*)fiSignalsPS->Get("combinedws");

   //TFile *fi = TFile::Open("limits-output.root");
   //TFile *fi = TFile::Open("signal-scans.root");
   TFile *fi = TFile::Open(fin.c_str());
   TTree *tree = (TTree*)fi->Get("limit");

   double mh;
   double limit;
   float quantile;
   tree->SetBranchAddress("mh",&mh);
   tree->SetBranchAddress("limit",&limit);
   tree->SetBranchAddress("quantileExpected",&quantile);

   int nvt = tree->GetEntries();
   
   TGraph2D *grV = new TGraph2D(); grV->SetName("vector");
   TGraph2D *grA = new TGraph2D(); grA->SetName("axial");
   TGraph2D *grS = new TGraph2D(); grS->SetName("scalar");
   TGraph2D *grP = new TGraph2D(); grP->SetName("pseudoscalar");


   TGraph2D *grVs = new TGraph2D(); grVs->SetName("vector_signal");
   TGraph2D *grAs = new TGraph2D(); grAs->SetName("axial_signal");
   TGraph2D *grSs = new TGraph2D(); grSs->SetName("scalar_signal");
   TGraph2D *grPs = new TGraph2D(); grPs->SetName("pseudoscalar_signal");

   int ptV=0;
   int ptA=0;
   int ptS=0;
   int ptP=0;

   TGraph *grV_mMED = new TGraph(); grV_mMED->SetName(Form("vector_mMED_mDM%d",keepMDM));
   TGraph *grA_mMED = new TGraph(); grA_mMED->SetName(Form("axial_mMED_mDM%d",keepMDM));
   TGraph *grS_mMED = new TGraph(); grS_mMED->SetName(Form("scalar_mMED_mDM%d",keepMDM));
   TGraph *grP_mMED = new TGraph(); grP_mMED->SetName(Form("pseudoscalar_mMED_mDM%d",keepMDM));
   int ptVm=0;
   int ptAm=0;
   int ptSm=0;
   int ptPm=0;

   grV_mMED->SetLineWidth(2); grV_mMED->SetMarkerSize(1.0);
   grA_mMED->SetLineWidth(2); grA_mMED->SetMarkerSize(1.0);
   grS_mMED->SetLineWidth(2); grS_mMED->SetMarkerSize(1.0);
   grP_mMED->SetLineWidth(2); grP_mMED->SetMarkerSize(1.0);

   for (int i=0; i<nvt;i++){
     //if ( ! ( i%6==X ) ) continue; // 2 or 5
     tree->GetEntry(i);
     if (quantile!=X) continue;
      
     int cd = code(mh);

     float mmed = MMED(mh,cd);
     float mdm  = MDM(mh,cd);

     // onshell crazyness?
     if ((int)mmed==2*( (int)mdm) ) continue;
     //
     //std::cout << " int X = " << i << std::endl;
     //std::cout << mh << ", " << cd << ", " << mmed << ", " << mdm <<  ", " << limit << std::endl;  
     //std::cout << mh << ", " << Form("monojet_signal_signal_%3d%04d%04d",cd,(int)mmed,(int)mdm) << std::endl; 
     //exit();
     
     //std::vector<std::pair<double,double>> pointsV_mMED;
     //std::vector<std::pair<double,double>> pointsA_mMED;
     //std::vector<std::pair<double,double>> pointsS_mMED;
     //std::vector<std::pair<double,double>> pointsP_mMED;

     RooDataHist *dh; 
     if      (cd==801 || cd==800) dh = (RooDataHist*) workspace->data(Form("monojet_signal_signal_%3d%04d%04d",cd,(int)mmed,(int)mdm));
     else if (cd==805 || cd==806) dh = (RooDataHist*) workspacePS->data(Form("monojet_signal_signal_%3d%04d%04d",cd,(int)mmed,(int)mdm));
     double nsignal = 0;
     if (dh) {
     	nsignal = dh->sumEntries();
     }
     if (cd==800) {
	grVs->SetPoint(ptV,mmed,mdm,nsignal);
     	grV->SetPoint(ptV,mmed,mdm,limit);
	ptV++;
	if ( (int)mdm == keepMDM ) { 
		grV_mMED->SetPoint(ptVm,mmed,limit);
		//pointsV_mMED.push_back(std::mk_pair<double,double>(mmed,limit));
		ptVm++;
	}
     } else if (cd==801){
	grAs->SetPoint(ptA,mmed,mdm,nsignal);
     	grA->SetPoint(ptA,mmed,mdm,limit);
	ptA++;
	if ( (int)mdm == keepMDM ) { 
		grA_mMED->SetPoint(ptAm,mmed,limit);
		ptAm++;
//.........这里部分代码省略.........

void PlotSVMopt(){

 //  // TString FileData = "SVMoptData.dat";
//   TFile *f = new TFile("SVMoptTree.root","RECREATE");
//   TTree *T = new TTree("SVMoptTree","data from ascii file");
//   Long64_t nlines = T->ReadFile("SVMoptData.dat","Gamma:C:ROCint:SigAt1Bkg_test:SigAt1Bkg_train");
//   printf(" found %lld pointsn",nlines);
//   T->Write();

// TGraph2D *grROCint = new TGraph2D();
//  for(int i=0;i<nlines;i++){

//    grROCint->SetPoint(i,);
//  }

  TGraph2D *grROCint = new TGraph2D();
  

TGraph2D *grSigAt1Bkg_test = new TGraph2D();
TGraph2D *grOverTrain = new TGraph2D();

 TString dir = gSystem->UnixPathName(__FILE__);
 dir.ReplaceAll("PlotSVMopt.C","");
 dir.ReplaceAll("/./","/");
 ifstream in;
 // in.open(Form("%sSVMoptData.dat",dir.Data()));
 in.open(Form("%sSVMoptData_NEW2.dat",dir.Data()));

 Float_t Gamma,C,ROCint, SigAt1Bkg_test, SigAt1Bkg_train;
 Int_t nlines = 0;
 while (1) {
   in >> Gamma >> C >> ROCint >> SigAt1Bkg_test >> SigAt1Bkg_train;
   if (!in.good()) break;
   // if (in.good()){
      if (nlines < 5) printf("Gamma=%8f, C=%8f, ROCint=%8f\n",Gamma, C, ROCint);
      grROCint->SetPoint(nlines,Gamma,C,ROCint);
      grSigAt1Bkg_test->SetPoint(nlines,Gamma,C,SigAt1Bkg_test);
      grOverTrain->SetPoint(nlines,Gamma,C,fabs(SigAt1Bkg_train-SigAt1Bkg_test));
      nlines++;
      //   }
   }
   printf(" found %d points\n",nlines);
   in.close();
   TCanvas *cSVMopt = new TCanvas("cSVMopt","SVM model choice",600,600);
   cSVMopt->Divide(2,2);
   cSVMopt->cd(1);
   cSVMopt_1->SetLogx();
   cSVMopt_1->SetLogy();
   grROCint->GetXaxis()->SetLabelSize(0.04);
   grROCint->GetYaxis()->SetLabelSize(0.04);
   grROCint->GetZaxis()->SetLabelSize(0.04);
   grROCint->GetXaxis()->SetTitle("#gamma");
   grROCint->GetYaxis()->SetTitle("C");
   grROCint->GetZaxis()->SetTitle("ROC integral");
   grROCint->SetTitle("ROC integral");
   grROCint->SetMaximum(1.0);
   grROCint->SetMinimum(0.9);
   grROCint->Draw("COLZ");
   cSVMopt->cd(2);
   cSVMopt_2->SetLogx();
   cSVMopt_2->SetLogy();
   grSigAt1Bkg_test->SetTitle("Signal at 1% Bkg level (test)");
   //   grSigAt1Bkg_test->Draw("surf1");
   grSigAt1Bkg_test->Draw("COLZ");
   cSVMopt->cd(3);
   cSVMopt_3->SetLogx();
   cSVMopt_3->SetLogy();
   grOverTrain->SetTitle("Overtraining");
   grOverTrain->Draw("COLZ");
   //cSVMopt->SaveAs("SVMoptC1.root");
}

//.........这里部分代码省略.........
			
			std::cout <<"size of int " <<  intb.size()<< std::endl;
			for (size_t it = 0; it < intb.size() - 1; it++){
				double y_val = (intb[it]+intb[it+1])/2;
				assert (y_val >= 0);
				double area = 0.002*y_val;
				if(area > 0 )
					s_integral += area;
			}
			
			
			
			
			std::cout << "Simpson integral: " <<s_integral <<std::endl;
			integral_hist->Fill(s_integral);
			cmp_int[i] = s_integral;
			Int_t lines = (Int_t)intb.size();
			TGraph *r_integral = new TGraph(lines, _inta, _intb);
			
			std::cout << "ROOT integral: " << r_integral->Integral() << std::endl;
			cmp_int_root[i] = r_integral->Integral();
			
				//expanding
				//expand(y, THRS_EXPAND, RATIO_EXPAND, LINES);
			
			
				//Filling TGraph2D
			
			for(Int_t j = 0; j <LINES ; j++){
				if (y[j] > max){
					max = y[j];
					maxwl = x[j];
				}
				gr->SetPoint(j+i*LINES, x[j],i,y[j]);
			}
			
			
			in.seekg(0, std::ios::beg);
			in.close();
			
				//Plotting each spectrum
			
			TGraph *_gr = new TGraph(LINES,x,y);
			_gr->GetHistogram()->GetXaxis()->SetTitle("#lambda in nm");
			_gr->GetHistogram()->GetYaxis()->SetTitle("Intensity in dB");
			c1->cd(i-NUM_ARGS);
			_gr->Draw("AP");
			_gr->GetYaxis()->SetRangeUser(-80.,-10.);
			_gr->GetXaxis()->SetRangeUser(startwl,stopwl);
			_gr->SetTitle(tmp.c_str());
			c1->Update();
			
			
				//Calculating asymmetry
			std::cout << "maximum: " << max << std::endl;
			double leftlimit, rightlimit = 1;
			leftlimit = findlower(x,y, max);
			rightlimit = findupper(x,y, max);
			if (leftlimit != 1 && rightlimit != 1){
				width_ary[i] = (leftlimit +rightlimit)/2;
			}else{
				width_ary[i] = maxwl;
			}
			double calced_asy = (maxwl-leftlimit)/(rightlimit-maxwl);
			asymmety_ary[i-NUM_ARGS] = calced_asy;
			

//.........这里部分代码省略.........
    //c1_5->SetTicks(0,0);
    //c1_2->SetRightMargin(0.15);
    //c1_2->SetLeftMargin(0.15);
    //c1_2->SetBottomMargin(0.02);
    lm10->Draw("colz1");
    //text->DrawLatex(-3,1,"SM plane in log scale");
    //text->SetTextSize(0.08);
    //text->SetTextColor(0);
    text->DrawLatex(-3.,0.7,"#kappa_{#lambda} = -10, c_{g}  = c_{2} = 0");  
    c1->SaveAs("C2Fit.pdf");
    c1->Close();
     */
    //////////////////////////////////////////////////
    //
    // do histrograms with errors
    //
    // plot (point - fit)/fit between int nmintest, int nmaxtest
    // do by the planes
    //////////////////////////////////////////////////
    // take the fit
    // need to be done by planes
    //c1->Clear();
    // a
    double SMxs =  0.013531; // 1 0.017278;// 14 0.0041758;// 8tev 0.013531; // 13 tev 0.017278;// 0.0041758;
    TGraph2D *g2 = new TGraph2D(117);//(118);
    g2->SetMarkerStyle(20);
    g2->SetMarkerSize(2);
    g2->SetTitle("0");
    g2->SetTitle("#kappa_{t} = #kappa_{#lambda} = 1 , c_{2} = 0 ; c_{g} ; c_{2g}");
    int j=0;
    for (unsigned int ij = 0; ij < nmaxx ; ij++) if( par1[ij] ==1 && par0[ij] ==1 && par2[ij]==0 && cross_section[ij] >0.0001) if(ij!=301) {
        double fit = SMxs*(fg2->Eval(par3[ij], par4[ij]));
        cout<<j<<" "<< par3[ij]<<" "<< par4[ij]<<" "<<fit <<" "<< cross_section[ij]<<" diff: " <<(fit - cross_section[ij])/fit<< endl;
        g2->SetPoint(j, par3[ij], par4[ij], 100*(fit - cross_section[ij])/fit); 
        j++;
        //Differences2->Fill(par3[i], par4[i], (fit - cross_section[i])/fit); 
    }
    // b 
    ////////////////////////////////
    int ktb=1.0;
    int klb=1.0;
    // cg ===> x ==> c2
    // c2g ===> y ==> kt ==> cg = c2g
    TF2 *pb = new TF2("pb","([0]*[15]**4 + [1]*x**2 + [2]*[15]**2*[16]**2 + [3]*y**2*[16]**2 +  [4]*y**2 + [5]*x*[15]**2 + [6]*[15]*[16]*[15]**2 + [7]*[15]*[16]*x + [8]*y*[16]*x - [9]*x*y + [10]*y*[16]*[15]**2 - [11]*y*[15]**2 + [12]*[16]*y*[15]*[16] - [13]*y*[15]*[16] - [14]*y*y*[16])/[17]",-3,3,-1,1);
    pb->SetParameter(0,a[0]);
    pb->SetParameter(1,a[1]);
    pb->SetParameter(2,a[2]);    
    pb->SetParameter(3,a[3]);
    pb->SetParameter(4,a[4]);
    pb->SetParameter(5,a[5]); 
    pb->SetParameter(6,a[6]);
    pb->SetParameter(7,a[7]);
    pb->SetParameter(8,a[8]); 
    pb->SetParameter(9,a[9]);
    pb->SetParameter(10,a[10]);
    pb->SetParameter(11,a[11]); 
    pb->SetParameter(12,a[12]);    
    pb->SetParameter(13,a[13]);
    pb->SetParameter(14,a[14]);
    pb->SetTitle("#kappa_{t} = #kappa_{#lambda} = 1 , c_{2g} = - c_{g} ; c_{2} ; c_{g}");
    pb->SetParameter(15,ktb); //==> c2g ==>kt
    pb->SetParameter(16,klb);
    pb->SetContour(200);
    //l5->SetParameter(17,cg); //==> cg
    pb->SetParameter(17,norm);//0.013531
    //pb->SetMinimum(0);

Int_t line3dfit_copy()
{
   gStyle->SetOptStat(0);
   gStyle->SetOptFit();


   //double e = 0.1;
   // Int_t nd = 5;


//    double xmin = 0; double ymin = 0;
//    double xmax = 10; double ymax = 10;

   TGraph2D * gr = new TGraph2D();

   // Fill the 2D graph
   // double p0[4] = {10,20,1,2};

   // generate graph with the 3d points
   // for (Int_t N=0; N<nd; N++) {
   //    double x,y,z = 0;
   //    // Generate a random number
   //    double t = gRandom->Uniform(0,10);
   //    line(t,p0,x,y,z);
   //    double err = 1;
   //  // do a gaussian smearing around the points in all coordinates
   //    x += gRandom->Gaus(0,err);
   //    y += gRandom->Gaus(0,err);
   //    z += gRandom->Gaus(0,err);
   //    gr->SetPoint(N,x,y,z);
   //    //dt->SetPointError(N,0,0,err);
   // }


   // gr->SetPoint(0, 3, 3, 5);
   // gr->SetPoint(1, 3, 3, 3);
   // gr->SetPoint(2, 4, 3.5, 5);
   // gr->SetPoint(3, 6.2, 7.3, 5.8);
   // gr->SetPoint(4, 9, 8, 7);
   // gr->SetPoint(5, 5, 5, 5);

   gr->SetPoint(0, 19, -4.25, 1.92);
   gr->SetPoint(1, 19, -4.30, 1.92);
   gr->SetPoint(2, 19, -4.35, 1.92);
   gr->SetPoint(3, 19, -4.40, 1.92);
   gr->SetPoint(4, 19, -4.45, 1.92);
   gr->SetPoint(5, 19, -4.50, 1.92);
   gr->SetPoint(6, 19, -4.55, 1.92);
   gr->SetPoint(7, 19, -4.60, 1.92);
   gr->SetPoint(8, 18.75, -4.30, 1.92);
   gr->SetPoint(9, 18.75, -4.35, 1.92);
   gr->SetPoint(10, 18.75, -4.40, 1.92);
   gr->SetPoint(11, 18.75, -4.45, 1.92);
   gr->SetPoint(12, 18.75, -4.50, 1.92);
   gr->SetPoint(13, 18.75, -4.55, 1.92);
   gr->SetPoint(14, 19, -4.20, 2.08);
   gr->SetPoint(15, 19, -4.65, 2.08);
   gr->SetPoint(16, 19, -4.70, 2.08);
   gr->SetPoint(17, 18.75, -4.25, 2.08);
   gr->SetPoint(18, 18.75, -4.60, 2.08);
   // gr->SetPoint(19, 19.001, -4.151, 2.241);
   // gr->SetPoint(20, 18.751, -4.201, 2.241);
   gr->SetPoint(19, 19., -4.15, 2.24);
   gr->SetPoint(20, 18.5723, -4.20, 2.24);

   // gr->SetMarkerStyle(8);
   // gr->Draw("p");
   
   // gr->Draw("p0");
   // TFitResultPtr fit = gr->Fit("pol1", "WS");
   // // fit->Print("V");
   // Double_t p0 = fit->Value(0);
   // Double_t p1 = fit->Value(1);

   // // draw the line
   // TPolyLine3D *l = new TPolyLine3D(2);
   // double dz = 8;
   // l->SetPoint(0,0,0,p0);
   // l->SetPoint(1,dz,0,dz * p1);

   // l->SetLineColor(kRed);
   // l->Draw("same");

   // fit the graph now, and make the functor objet

   ROOT::Fit::Fitter  fitter;

   SumDistance2 sdist(gr);
#ifdef __CINT__
   ROOT::Math::Functor fcn(&sdist,4,"SumDistance2");
#else
   ROOT::Math::Functor fcn(sdist,4);
#endif
   // set the function and the initial parameter values
   double pStart[4] = {1,1,1,1};
   fitter.SetFCN(fcn,pStart);
   // set step sizes different than default ones (0.3 times parameter values)
   for (int i = 0; i < 4; ++i) fitter.Config().ParSettings(i).SetStepSize(0.01);

   bool ok = fitter.FitFCN();
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    fcint = fc.GetInterval();  // that was easy.

    RooFitResult* fit = modelWithConstraints->fitTo(*data, Save(true));

    // Third, use a Calculator based on Markov Chain monte carlo
    // Before configuring the calculator, let's make a ProposalFunction
    // that will achieve a high acceptance rate
    ProposalHelper ph;
    ph.SetVariables((RooArgSet&)fit->floatParsFinal());
    ph.SetCovMatrix(fit->covarianceMatrix());
    ph.SetUpdateProposalParameters(true);
    ph.SetCacheSize(100);
    ProposalFunction* pdfProp = ph.GetProposalFunction();  // that was easy

    MCMCCalculator mc(*data, modelConfig);
    mc.SetNumIters(20000); // steps to propose in the chain
    mc.SetTestSize(.05); // 95% CL
    mc.SetNumBurnInSteps(40); // ignore first N steps in chain as "burn in"
    mc.SetProposalFunction(*pdfProp);
    mc.SetLeftSideTailFraction(0.5);  // find a "central" interval
    MCMCInterval* mcInt = (MCMCInterval*)mc.GetInterval();  // that was easy


    // Get Lower and Upper limits from Profile Calculator
    cout << "Profile lower limit on s = " << ((LikelihoodInterval*) lrinterval)->LowerLimit(*s) << endl;
    cout << "Profile upper limit on s = " << ((LikelihoodInterval*) lrinterval)->UpperLimit(*s) << endl;

    // Get Lower and Upper limits from FeldmanCousins with profile construction
    if (fcint != NULL) {
        double fcul = ((PointSetInterval*) fcint)->UpperLimit(*s);
        double fcll = ((PointSetInterval*) fcint)->LowerLimit(*s);
        cout << "FC lower limit on s = " << fcll << endl;
        cout << "FC upper limit on s = " << fcul << endl;
        TLine* fcllLine = new TLine(fcll, 0, fcll, 1);
        TLine* fculLine = new TLine(fcul, 0, fcul, 1);
        fcllLine->SetLineColor(kRed);
        fculLine->SetLineColor(kRed);
        fcllLine->Draw("same");
        fculLine->Draw("same");
        dataCanvas->Update();
    }

    // Plot MCMC interval and print some statistics
    MCMCIntervalPlot mcPlot(*mcInt);
    mcPlot.SetLineColor(kMagenta);
    mcPlot.SetLineWidth(2);
    mcPlot.Draw("same");

    double mcul = mcInt->UpperLimit(*s);
    double mcll = mcInt->LowerLimit(*s);
    cout << "MCMC lower limit on s = " << mcll << endl;
    cout << "MCMC upper limit on s = " << mcul << endl;
    cout << "MCMC Actual confidence level: "
         << mcInt->GetActualConfidenceLevel() << endl;

    // 3-d plot of the parameter points
    dataCanvas->cd(2);
    // also plot the points in the markov chain
    RooDataSet * chainData = mcInt->GetChainAsDataSet();

    assert(chainData);
    std::cout << "plotting the chain data - nentries = " << chainData->numEntries() << std::endl;
    TTree* chain =  RooStats::GetAsTTree("chainTreeData","chainTreeData",*chainData);
    assert(chain);
    chain->SetMarkerStyle(6);
    chain->SetMarkerColor(kRed);

    chain->Draw("s:ratioSigEff:ratioBkgEff","nll_MarkovChain_local_","box"); // 3-d box proportional to posterior

    // the points used in the profile construction
    RooDataSet * parScanData = (RooDataSet*) fc.GetPointsToScan();
    assert(parScanData);
    std::cout << "plotting the scanned points used in the frequentist construction - npoints = " << parScanData->numEntries() << std::endl;
    // getting the tree and drawing it -crashes (very strange....);
    // TTree* parameterScan =  RooStats::GetAsTTree("parScanTreeData","parScanTreeData",*parScanData);
    // assert(parameterScan);
    // parameterScan->Draw("s:ratioSigEff:ratioBkgEff","","goff");
    TGraph2D *gr = new TGraph2D(parScanData->numEntries());
    for (int ievt = 0; ievt < parScanData->numEntries(); ++ievt) {
        const RooArgSet * evt = parScanData->get(ievt);
        double x = evt->getRealValue("ratioBkgEff");
        double y = evt->getRealValue("ratioSigEff");
        double z = evt->getRealValue("s");
        gr->SetPoint(ievt, x,y,z);
        // std::cout << ievt << "  " << x << "  " << y << "  " << z << std::endl;
    }
    gr->SetMarkerStyle(24);
    gr->Draw("P SAME");


    delete wspace;
    delete lrinterval;
    delete mcInt;
    delete fcint;
    delete data;

    // print timing info
    t.Stop();
    t.Print();
}