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

   // implementation of the function to be minimized
   double operator() (const double * par) {

      assert(fGraph    != 0);
      double * x = fGraph->GetX();
      double * y = fGraph->GetY();
      double * z = fGraph->GetZ();
      int npoints = fGraph->GetN();
      double sum = 0;
      for (int i  = 0; i < npoints; ++i) {
         double d = distance2(x[i],y[i],z[i],par);
         sum += d;
#ifdef DEBUG
         if (first) std::cout << "point " << i << "\t"
            << x[i] << "\t"
            << y[i] << "\t"
            << z[i] << "\t"
            << std::sqrt(d) << std::endl;
#endif
      }
      
      // if (first)
      // 	std::cout << "Total distance square w/ initial params: " << sum << std::endl;
      // 	first = false;
      return sum;
   }

void Analyzer::ComputeMin(){
	g2=new TGraph2D(); //TODO
	
	alpha=1.0;beta=0;
	Loop(t_data,1);
	if(varName=="QGLMLP")
		Loop(t_mc,4);
	//scan
	alpha=1.0;beta=0;
	for(float ai=0.7; ai<=1.1; ai+=0.02)
		{
		Reset(h_mc);	
		alpha=ai;
		Loop(t_mc,2);
		h_mc->Scale(h_data->Integral()/h_mc->Integral());
		g2->SetPoint(g2->GetN(),alpha,beta, h_data->Chi2Test(h_mc,opt.c_str())  );	
		}
	alpha=1.0;beta=0;
	for(float bi=-0.5; bi<=0.5; bi+=0.01)
		{
		Reset(h_mc);	
		beta=bi;
		Loop(t_mc,2);
		h_mc->Scale(h_data->Integral()/h_mc->Integral());
		g2->SetPoint(g2->GetN(),alpha,beta, h_data->Chi2Test(h_mc,opt.c_str())  );	
		}
	
	//Find min0;min1
	float min0=1,min1=0;
	pair<float,float> R=MinG(g2);
	min0=R.first;min1=R.second;

	for(int i=-nstep;i<=nstep;i++)
	for(int j=-nstep;j<=nstep;j++)
        	{
        	alpha=min0+i*stp0;
        	beta=min1+j*stp1;
		Loop(t_mc,2);
		h_mc->Scale(h_data->Integral()/h_mc->Integral());
		g2->SetPoint(g2->GetN(),alpha,beta, h_data->Chi2Test(h_mc,opt.c_str())  );
		}
	//double min0,min1;	
	R=MinG(g2);
	printf("a=%.3f;b=%.3f;lmin=%.3f;lmax=%.3f;break;\n",R.first,R.second,lmin,lmax);
	return;
	
}

int
plotting_output()
{

	vector< int > x, y;
	vector< double > laplace;
	string mapping_tower_file = "./numerical_solution.txt";


	/* Stream to read table from file */
	ifstream istream_mapping;

	/* Open the datafile, if it won't open return an error */
	if (!istream_mapping.is_open())
	{
		istream_mapping.open( mapping_tower_file.c_str() );
		if(!istream_mapping)
		{
			cerr << "CaloTowerGeomManager::ReadGeometryFromTable - ERROR Failed to open mapping file " << mapping_tower_file << endl;
			exit(1);
		}
	}

	string line_mapping;
	int x_i, y_i;
	double laplace_i;

	while ( getline( istream_mapping, line_mapping ) )
	{
		istringstream iss(line_mapping);
		iss >> x_i >> y_i >> laplace_i;
		x.push_back( x_i );
		y.push_back( y_i );
		laplace.push_back( laplace_i );
	}

	TCanvas *c1 = new TCanvas();
	TGraph2D *gr = new TGraph2D();
	for(int i = 0; i < x.size(); i++)
	{
		gr->SetPoint( gr->GetN(), x.at(i), y.at(i), laplace.at(i) );
	}
	gr->Draw("surf1");
	cout << "graph drawn " << endl;

	TCanvas *c1 = new TCanvas();
	gr->Draw("P");





	return 0;
}

pair<float, float> Analyzer::MinG(TGraph2D *g,double *min0,double*min1){
	pair<float,float> R(-99,-99);
	if(g==NULL) return R;
	if(g->GetN()==0)return R;
	double *x1,*y1,*z1;
	x1=g->GetX();
	y1=g->GetY();
	z1=g->GetZ();
	
	float a=z1[0];R=pair<float,float>(x1[0],y1[0]);
	for(int i=0;i<g->GetN();i++){if((z1[i]<a)||(a<0)){a=z1[i];    R=pair<float,float>(x1[i],y1[i]); }}
	
	if(min0!=NULL) *min0=a;
	if(min1!=NULL) *min1=z1[g2->GetN()-1]	;
	return R;
}

void Analyzer::ComputeMinFast(){
	g2=new TGraph2D(); 
	
	alpha=1.0;beta=0;
	Loop(t_data,1);
	for(int j=0;j<=h_data->GetNbinsX()+1;j++)if(h_data->GetBinError(j)==0)h_data->SetBinError(j,1);
	if(varName=="QGLMLP")
		Loop(t_mc,4);
	//scan
	//reset Fast
	ResetFast();
	
	alpha=1.0;beta=0;
	for(float ai=aMin; ai<=aMax; ai+=0.02)
		{
		Reset(h_mc);	
		alphaFast.push_back(ai);	
		betaFast.push_back(0);	
		h_mcFast.push_back(new TH1F( Form("hmc_%d",int(h_mcFast.size())),"hmc",nBins,xMin,xMax) );
		}
	alpha=1.0;beta=0;
	for(float bi=bMin; bi<=bMax; bi+=0.01)
		{
		Reset(h_mc);	
		alphaFast.push_back(1.0);	
		betaFast.push_back(bi);	
		h_mcFast.push_back(new TH1F( Form("hmc_%d",int(h_mcFast.size())),"hmc",nBins,xMin,xMax) );
		}

	for(int j=0;j< int(h_mcFast.size());j++) h_mcFast[j]->Sumw2();	
	Loop(t_mc,16);
		for(int i=0 ;i<int(alphaFast.size());i++)
			{
			for(int j=0;j<=h_mcFast[i]->GetNbinsX()+1;j++)if(h_mcFast[i]->GetBinError(j)==0)h_mcFast[i]->SetBinError(j,1);
			h_mcFast[i]->Scale(h_data->Integral()/h_mcFast[i]->Integral());
			g2->SetPoint(g2->GetN(),alphaFast[i],betaFast[i], h_data->Chi2Test(h_mcFast[i],opt.c_str())  );
			}
	//Find min0;min1
	float min0=1,min1=0;
	pair<float,float> R=MinG(g2);
	min0=R.first;min1=R.second;
		
	ResetFast();

	delete g2;
	g2=new TGraph2D();

	for(int i=-nstep;i<=nstep;i++)
	for(int j=-nstep;j<=nstep;j++)
        	{
        	alphaFast.push_back(min0+i*stp0 );
        	betaFast.push_back(min1+j*stp1 );
		h_mcFast.push_back(new TH1F( Form("hmc_%d",int(h_mcFast.size())),"hmc",nBins,xMin,xMax) );
		//g2->SetPoint(g2->GetN(),alpha,beta, h_data->Chi2Test(h_mc,opt.c_str())  );
		}
	alphaFast.push_back(min0);
	betaFast.push_back(0);
	h_mcFast.push_back(new TH1F( Form("hmc_%d",int(h_mcFast.size())),"hmc",nBins,xMin,xMax) );

	alphaFast.push_back(min1);
	betaFast.push_back(bMax);
	h_mcFast.push_back(new TH1F( Form("hmc_%d",int(h_mcFast.size())),"hmc",nBins,xMin,xMax) );
	
	alphaFast.push_back(1);
	betaFast.push_back(0);
	h_mcFast.push_back(new TH1F( Form("hmc_%d",int(h_mcFast.size())),"hmc",nBins,xMin,xMax) );
	
	for(int j=0;j<int(h_mcFast.size());j++) h_mcFast[j]->Sumw2();	
	Loop(t_mc,16);
		for(int i=0 ;i<int(alphaFast.size());i++)
			{
			for(int j=0;j<=h_mcFast[i]->GetNbinsX()+1;j++)if(h_mcFast[i]->GetBinError(j)==0)h_mcFast[i]->SetBinError(j,1);
			h_mcFast[i]->Scale(h_data->Integral()/h_mcFast[i]->Integral());
			g2->SetPoint(g2->GetN(),alphaFast[i],betaFast[i], h_data->Chi2Test(h_mcFast[i],opt.c_str())  );
			}
	double m0,m1;
	R=MinG(g2,&m0,&m1);
	printf("a=%.3f;b=%.3f;lmin=%.3f;lmax=%.3f;break;//chi2=%.3lf; chi2_0=%.3lf\n",R.first,R.second,lmin,lmax,m0,m1);
	{
	TFile *out=TFile::Open("output.root","UPDATE");out->cd();
	for(int i=0;i<int(h_mcFast.size());i++)
		{
		h_mcFast[i]->SetName(Form("%s_alpha%.2f_beta%.2f_lmin%.3f_lmax%.3f_pt%0f_%.0f_rho%.0f_%.0f_eta%.0f_%.0f",varName.c_str(),alphaFast[i],betaFast[i],lmin,lmax,PtMin,PtMax,RhoMin,RhoMax,EtaMin,EtaMax));
		h_mcFast[i]->Write();
		}
	h_data->SetName(Form("%s_data_pt%0f_%.0f_rho%.0f_%.0f_eta%.0f_%.0f",varName.c_str(),PtMin,PtMax,RhoMin,RhoMax,EtaMin,EtaMax));
	h_data->Write();
	}
	return;
}

//.........这里部分代码省略.........
        laplace[1][0] = v_left / 2;
        laplace[2][0] = v_left / 2;
        laplace[3][0] = v_left / 2;
        laplace[1][40] = v_left;
        laplace[2][40] = v_left / 2;
        laplace[3][40] = v_left / 2;
        laplace[4][40] = v_left / 4;
        cout << "    (additional Constraints applied)" << endl;
    }

    //Write out the initial matrix
    if(write_matrix_init == true)
    {
        for(int i = 0; i < x_max+2; i++)
        {
            cout << " | ";
            for(int j = 0; j < y_max+2; j++)
            {
                cout << laplace[i][j] << " | ";
            }
            cout << endl;
        }
        cout << endl << endl;
    }


    //Iterate over the inner rows/columns of the matrix, averaging at each point.
    for ( int k = 0; k < iterations; k++)
    {
        //Loop over all iterations, averaging each point with those around them... skipping over boundary conditions
        for(int i = 1; i < x_max+1; i++)
        {
            for(int j = 1; j < y_max+1; j++)
            {
                laplace[i][j] = 0.25*( laplace[i-1][j] + laplace[i+1][j] + laplace[i][j-1] + laplace[i][j+1] );
//				laplace[i][j] = 0.125*( laplace[i-1][j] + laplace[i+1][j] + laplace[i][j-1] + laplace[i][j+1] + laplace[i-1][j-1] + laplace[i-1][j+1] + laplace[i+1][j-1] + laplace [i+1][j+1] );
            }
        }
    }
    cout << "Computations Completed. Plotting..." << endl;

    //Write out the final matrix
    if(write_matrix_final == true)
    {
        for(int i = 0; i < x_max+2; i++)
        {
            cout << " | ";
            for(int j = 0; j < y_max+2; j++)
            {
                cout << laplace[i][j] << " | ";
            }
            cout << endl;
        }
        cout << endl << endl;
    }

    //Write the output to a file
    if(output)
    {
        ofstream myfile;
        myfile.open("numerical_solution.txt");

        for(int i = 1; i < x_max+1; i++)
        {
            for(int j = 1; j < y_max+1; j++)
            {
                myfile << i << " " << j << " " << laplace[i][j] << endl;
//				myfile << j << " " << i << " " << laplace[i][j] << endl;
            }
        }
        myfile.close();
        cout << "Output file created." << endl;
    }


    //Plot the final output
    if(plot)
    {
        //Create the output graph
        TCanvas *c1 = new TCanvas();
        TGraph2D *gr = new TGraph2D();
        for(int i = 1; i < x_max+1; i++)
        {
            for(int j = 1; j < y_max+1; j++)
            {
                gr->SetPoint( gr->GetN(), i, j, laplace[i][j] );
            }
        }
        //gr->Draw("surf1");
        gr->Draw("COLZ");

        TCanvas *c2 = new TCanvas();
        gr->Draw("surf1");

        //	TH2 *h2 = gr->Project("xy");
        //	h2->Draw();
    }

    return 0;
}

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 fillGraphsFromFiles( const TString& par1,
			  const TString& par2,
			  const vector<TString>& fnames,
			  vector<string>&  keys,
			  map<string,TGraph2D *>& m_graphs)
{
  keys.push_back("-2s");
  keys.push_back("-1s");
  keys.push_back("median");
  keys.push_back("+1s");
  keys.push_back("+2s");
  keys.push_back("obs");

  for (int i=0; i<6; i++) {
    m_graphs[keys[i]] = new TGraph2D();
    m_graphs[keys[i]]->SetName(Form("graph2D%s",keys[i].c_str()));
  }

  int nobs=0, nexp=0;
  for( size_t i=0; i<fnames.size(); i++) {
    
    double par1val = extractParValue(par1,fnames[i]);
    double par2val = extractParValue(par2,fnames[i]);

    //cout << par1val << "\t" << par2val << endl;

    if (par1val == -9e99 ||
	par2val == -9e99)
      continue;
    
    TFile *f = new TFile(fnames[i]);

    double obs,median,s68hi,s68lo,s95hi,s95lo;
    int num = getBands(f,1,0,obs,median,s68hi,s68lo,s95hi,s95lo);

    switch (num) {
    case 0: break;
    case 1:
      //cout << "SetPoint(i="<<nobs<<",par1="<<par1val<<",par2="<<par2val<<",obs="<<obs<<");"<<endl;
      m_graphs["obs"]->SetPoint(nobs,par1val,par2val,obs);
      nobs++;
      break;
    default:
      m_graphs["+2s"]->SetPoint(nexp,par1val,par2val,s95hi);
      m_graphs["-2s"]->SetPoint(nexp,par1val,par2val,s95lo);
      m_graphs["+1s"]->SetPoint(nexp,par1val,par2val,s68hi);
      m_graphs["-1s"]->SetPoint(nexp,par1val,par2val,s68lo);
      m_graphs["median"]->SetPoint(nexp,par1val,par2val,median);
      nexp++;
      break;
    }
        
    f->Close();

    delete f;

    //if (!(i%10)) cout << i << " " << std::flush;


  } // file loop

#if 0
  TGraph2D *gobs = m_graphs["obs"];
  cout << "obs has " << gobs->GetN() << " points" << endl;
  cout << "npx = " << gobs->GetNpx() << endl;
  cout << "npy = " << gobs->GetNpy() << endl;
  cout << "xmin = " << gobs->GetXmin() << endl;
  cout << "xmax = " << gobs->GetXmax() << endl;
  cout << "ymin = " << gobs->GetYmin() << endl;
  cout << "ymax = " << gobs->GetYmax() << endl;
  cout << "zmin = " << gobs->GetZmin() << endl;
  cout << "zmax = " << gobs->GetZmax() << endl;

  double *xvec = gobs->GetX();
  double *yvec = gobs->GetY();
  double *zvec = gobs->GetZ();
  for (int i=0; i<gobs->GetN(); i++)
    printf("%lf\t%lf\t%lf\n",xvec[i],yvec[i],zvec[i]);
#endif
}                                                 // fillGraphsFromFiles

void atgcplotLimit()
{
  atgcstyle();

  gStyle->SetPalette(1);

  vector<TString> fnames;
  //getFileNames(fileglob, fnames);

  fnames.push_back("higgsCombineTest.MultiDimFit.mH120_expected.root");
  fnames.push_back("higgsCombineTest.MultiDimFit.mH120_measured.root");

  assert(fnames.size());

  TString   par1;
  TString   par2;

  // get names of coupling parameters from root filename
  //
  //if (fnames[0].Contains("lz",TString::kIgnoreCase)) {       par1 = TString("lZ");
  //  if (fnames[0].Contains("dkg",TString::kIgnoreCase))      par2 = TString("dkg");
  //  else if (fnames[0].Contains("dg1",TString::kIgnoreCase)) par2 = TString("dg1");
  //} else if  (fnames[0].Contains("dkg",TString::kIgnoreCase) &&
  //	      fnames[0].Contains("dg1",TString::kIgnoreCase)) {
  //  par1 = TString("dkg");
  //  par2 = TString("dg1");
  //} 

  par1 = TString("lZ");
  par2 = TString("dkg");

  if (!par1.Length() || !par2.Length() ) {
    cerr << "Unknown coupling parameters in name " << fnames[0] << endl;
    exit(-1);
  }

  TString method("Other");
  if (fnames[0].Contains("Asymptotic"))
    method = TString("asympCLs");
  else if (fnames[0].Contains("MultiDimFit"))
    method = TString("deltaNLL");

  // try this:
  //  method = TString("asympCLs");
  

  cout << "Plotting " << par2 << " versus " << par1 << ", method = " << method << endl;

  vector<string> keys;
  map<string,double> m_contourlevels;
  map<string,TGraph2D *> m_graphs;

  if (method.EqualTo("asympCLs")) {
    fillGraphsFromFilesAsymp(par1,par2,fnames,keys,m_graphs);
    for (size_t i=0; i<keys.size(); i++)
      m_contourlevels[keys[i]] = 1;
  }
  if (method.EqualTo("deltaNLL")) {
    fillGraphsFromFilesDeltaNLL(par1,par2,fnames,keys,m_graphs);
    m_contourlevels["exp68"] = 2.3;
    m_contourlevels["exp95"] = 5.99;
    m_contourlevels["exp99"] = 9.21;
    m_contourlevels["obs95"] = 5.99;
  }
  else if (fnames.size() == 1) {  
    fillGraphsFromTextTables          (fnames[0],keys,m_graphs);
    for (size_t i=0; i<keys.size(); i++)
      m_contourlevels[keys[i]] = 0.05;
#if 0
    m_graphs["-2s"]->Draw("COLZ TEXT");
#else
    TGraph2D *gr = m_graphs["-2s"];
    TH2D *h2 = new TH2D("h2d","h2d",
			parbins(par1),parmin(par1)-parinc(par1)/2,parmax(par1)+parinc(par1)/2,
			parbins(par2),parmin(par2)-parinc(par2)/2,parmax(par2)+parinc(par2)/2);
    h2->FillN(gr->GetN(),gr->GetX(),gr->GetY(),gr->GetZ());
    h2->GetXaxis()->SetTitle(par2latex(par1));
    h2->GetYaxis()->SetTitle(par2latex(par2));
    h2->Draw("COLZ TEXT");
#endif
  }  else {
    fillGraphsFromFiles        (par1,par2,fnames,keys,m_graphs);
    for (size_t i=0; i<keys.size(); i++)
      m_contourlevels[keys[i]] = 1;
  }

  // try this
  /*
    fillGraphsFromFiles        (par1,par2,fnames,keys,m_graphs);
    for (size_t i=0; i<keys.size(); i++)
      m_contourlevels[keys[i]] = 1;
  */
  // return;

  // for limit in limits:
  //   limits[limit]->Print()


#if 0
  // for a first look
//.........这里部分代码省略.........

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);
//.........这里部分代码省略.........