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

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