void setGraphOptions(TGraph2D &g)
{
  g.SetTitle("");
  g.SetMarkerColor(1);
  g.SetMarkerStyle(24);
  g.SetMarkerSize(.5);
}

TGraph2D* preparedraw::data(){
  TGraph2D *toreturn = new TGraph2D(Z.size());
  for(int j=0;j<Nt;j++){
    for(int i=0;i<Nl;i++)
      toreturn->SetPoint(i+j*Nl,sstep*sskip*i,times.at(j),abs(Z.at(i+j*Nl)));
  }
  return toreturn;
}

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;
}

  bool PlotBundle::drawHistsAndGraphs(const PlotStyle &plot_style) const {
    if (hasDrawables()) {
      std::vector<DrawableDataObjectDrawOptionPair<TH1*> >::const_iterator hist_it;

      for (hist_it = histograms.begin(); hist_it != histograms.end();
          hist_it++) {
        TH1* hist = hist_it->data_object;
        std::string draw_option(hist_it->draw_option);
        if (hist) {
          if (hist->GetDimension() > 1) {
            hist->GetZaxis()->SetLimits(plot_axis.z_axis_range.low,
                plot_axis.z_axis_range.high);
            hist->GetZaxis()->SetRangeUser(plot_axis.z_axis_range.low,
                plot_axis.z_axis_range.high);
          }
          draw_option.append("SAME");
          hist->Draw(draw_option.c_str());
        }
      }

      std::vector<DrawableDataObjectDrawOptionPair<TGraph*> >::const_iterator graph_it;

      for (graph_it = graphs.begin(); graph_it != graphs.end(); graph_it++) {
        TGraph* graph = graph_it->data_object;
        std::string draw_option(graph_it->draw_option);
        if (graph) {
          draw_option.append("SAME");
          graph->Draw(draw_option.c_str());
        }
      }

      std::vector<DrawableDataObjectDrawOptionPair<TGraph2D*> >::const_iterator graph2d_it;

      for (graph2d_it = graphs2d.begin(); graph2d_it != graphs2d.end();
          graph2d_it++) {
        TGraph2D* graph = graph2d_it->data_object;
        std::string draw_option(graph2d_it->draw_option);
        if (graph) {
          draw_option.append("SAME");
          graph->Draw(draw_option.c_str());
        }
      }
      return true;
    } else {
      std::cout
          << "Dude, you forgot to add a histogram or graph in the plot bundle..."
          << " Please add at least one and make sure it points to an existing object!"
          << std::endl;
      return false;
    }
  }

//g++ allDataPrint.cpp `root-config --cflags --glibs`
int main(int argc, char** argv)
#endif
{
  TCanvas c3("c3","Grafico",640,512);
  TCanvas c1("c1","Confronto",1280,512);
  c1.Divide(2,1);
  
  preparedraw myData (argv[1],
		      // preparedraw::doMax |
		      preparedraw::doFh |
		      preparedraw::doSh |
		      preparedraw::doErr);
  
  TGraph2D *g = myData.data();
  TGraph *gb = myData.firsthalf();//before
  TGraph *ga = myData.secondhalf();//after
  TGraph *gerrs = myData.errs();
  //  TGraph *maxs = myData.maximum();
  
  c3.cd();
  g->GetXaxis()->SetTitle("X");
  g->GetYaxis()->SetTitle("T");
  //g->Draw("cont1");
  g->Draw("pcol");
  //g->Draw();
  //grafo.Draw("surf1");

  cout<<"Disegno i grafici\n";
  TMultiGraph *mg = new TMultiGraph("integrali","Integrali prima e dopo la barriera");
  
  ga->SetLineColor(2);
  mg->Add(gb);
  mg->Add(ga);
  c1.cd(1);
  mg->Draw("apl");
    

  gerrs->SetTitle("Andamento degli errori");
  c1.cd(2);
  gerrs->Draw("apl");
#ifndef __CINT__
  theApp.Run(true);
  return 0;
#endif
}

//g++ allDataPrint.cpp `root-config --cflags --glibs`
int main(int argv, char** argc)
#endif
{
  /*  TCanvas c("c","Real",640,512);
  TGraph2D *grafo = new TGraph2D("out.txt");
  grafo->Draw("pcol");
  TCanvas c2("c2","Imaginary",640,512);
  TGraph2D *Cgrafo = new TGraph2D("iout.txt");
  Cgrafo->Draw("pcol");*/
  TCanvas c3("c3","Norm",640,512);
  TGraph2D *Ngrafo = new TGraph2D(argc[1]);
  Ngrafo->Draw("pcol");
  //grafo.Draw("surf1");
#ifndef __CINT__
  theApp.Run(true);
  return 0;
#endif
}

TGraph2D* GetInterpolatingGraph(TH2F *hold){
  float binsize = hold->GetXaxis()->GetBinWidth(1)/2.;
  TString name = hold->GetName();
  name.ReplaceAll("Org","");
  TGraph2D *g = new TGraph2D(hold);
  g->SetNpx(int(g->GetXmax()-g->GetXmin())/binsize);
  g->SetNpy(int(g->GetYmax()-g->GetYmin())/binsize);
  //cout << "name " << g->GetN() << " " << hold->Integral() << endl;
  return g;
}

TCanvas* graph2dfit_test()
{
	gStyle->SetOptStat(0);
	gStyle->SetOptFit();

	TCanvas *c = new TCanvas("c", "Graph2D example", 0, 0, 600, 800);
	
	Double_t rnd, x, y, z;
	Double_t e = 0.3;
	Int_t nd = 400;
	Int_t np = 10000;

	TRandom r;
	Double_t fl = 6;
	TF2  *f2 = new TF2("f2", "1000*(([0]*sin(x)/x)*([1]*sin(y)/y))+200", -fl, fl, -fl, fl);
	f2->SetParameters(1, 1);
	TGraph2D *dt = new TGraph2D();

	// Fill the 2D graph
	Double_t zmax = 0;
	for (Int_t N = 0; N<nd; N++) {
		f2->GetRandom2(x, y);
		// Generate a random number in [-e,e]
		rnd = 2 * r.Rndm()*e - e;
		z = f2->Eval(x, y)*(1 + rnd);
		if (z>zmax) zmax = z;
		dt->SetPoint(N, x, y, z);
	}


	f2->SetParameters(0.5, 1.5);
	dt->Fit(f2);
	TF2 *fit2 = (TF2*)dt->FindObject("f2");

	f2->SetParameters(1, 1);

	for (Int_t N = 0; N<np; N++) {
		f2->GetRandom2(x, y);
		// Generate a random number in [-e,e]
		rnd = 2 * r.Rndm()*e - e;
		z = f2->Eval(x, y)*(1 + rnd);
		h1->Fill(f2->Eval(x, y) - z);
		z = dt->Interpolate(x, y);
		h2->Fill(f2->Eval(x, y) - z);
		z = fit2->Eval(x, y);
		h3->Fill(f2->Eval(x, y) - z);
	}

	gStyle->SetPalette(1);
	
	f2->SetTitle("Original function with Graph2D points on top");
	f2->SetMaximum(zmax);
	gStyle->SetHistTopMargin(0);
	f2->Draw("surf1");
	dt->Draw("same p0");



	return c;
}

TH2D grMu(TTree *tree, double quantileExpected){

	TGraph2D *gr = new TGraph2D();
	float X, Y, qE;
	double limit;

	std::string x = "trackedParam_proc_scaling_muV";
	std::string y = "trackedParam_proc_scaling_muF";

	tree->SetBranchAddress(Form("%s",x.c_str()),&X);
	tree->SetBranchAddress(Form("%s",y.c_str()),&Y);
	tree->SetBranchAddress("quantileExpected",&qE);
	tree->SetBranchAddress("limit",&limit);

	int pt = 0;
	for (int i=0;i<tree->GetEntries();i++){
	  tree->GetEntry(i);
	  if (TMath::Abs(qE-quantileExpected)>0.001) continue;
	  gr->SetPoint(pt,X,Y,limit); pt++;
	  //std::cout << " Lim " << X << ", " << Y << ", " << limit << std::endl;
	}
	gr->Draw("colz");
	TH2D * expected = (TH2D*) gr->GetHistogram();
	expected->SetTitle("");
	expected->GetZaxis()->SetTitleOffset(1.2);
	expected->GetYaxis()->SetTitle("#it{#mu}_{ggH}");
	expected->GetXaxis()->SetTitle("#it{#mu}_{qqH,VH}");
	expected->GetZaxis()->SetTitle("B(H #rightarrow inv.) - 95% CL upper limit");
	expected->GetXaxis()->SetRangeUser(MUVMIN,MUVMAX);
	expected->GetYaxis()->SetRangeUser(MUFMIN,MUFMAX);
	//expected->GetXaxis()->SetLimits(MUVMIN,MUVMAX);
	//expected->GetYaxis()->SetLimits(MUFMIN,MUFMAX);

	expected->SetMaximum(0.6);
	expected->SetMinimum(0.05);
	return *expected;
}

TH2F* InterpolateThisHistogram(TH2F *hold/*, TH2F* hnew*/){
  float binsize = hold->GetXaxis()->GetBinWidth(1)/2.;
  TString name = hold->GetName();
  name.ReplaceAll("Org","");
  TGraph2D *g = new TGraph2D(hold);
  //cout << g->GetXmax() << " " << g->GetXmin() << " " << g->GetYmax() << " " << g->GetYmin() << " " << binsize << endl;
  g->SetNpx(int(g->GetXmax()-g->GetXmin())/binsize);
  g->SetNpy(int(g->GetYmax()-g->GetYmin())/binsize);
  TH2F *hnew = (TH2F*)g->GetHistogram();
  //TH2F *htemp = (TH2F*)hnew->Clone(name);
  //name.ReplaceAll("YXZ","");
  TH2F *h = new TH2F(name.Data(),hold->GetTitle(),hnew->GetNbinsX(),g->GetXmin()-binsize,g->GetXmax()-binsize,hnew->GetNbinsY(),g->GetYmin()-binsize,g->GetYmax()-binsize);
  for(unsigned int x = 1; x<=hnew->GetNbinsX(); ++x){
    for(unsigned int y = 1; y<=hnew->GetNbinsY(); ++y){
      h->SetBinContent(x,y,hnew->GetBinContent(x,y));
    }
  }
  delete g;
  return h;
}

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::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 w1() {

	TCanvas *c = new TCanvas("c1", "", 0, 0, 700, 600);
	TGraph2D *dt = new TGraph2D();

	ifstream in("/opt/workspace-cpp/simulation/result_rand3.csv");

	string buff;

	getline(in, buff);

	int meanResidenceTime;
	int products;
	float lastProductPrice;
	float shopIncome;

	int n = 0;
	while (getline(in, buff)) {
			for (unsigned int i = 0; i < buff.size(); i++) {
				if (buff[i] == ';') {
					buff[i] = ' ';
				}
			}
			stringstream ss(buff);
			ss >> meanResidenceTime;
			ss >> products;
			ss >> lastProductPrice;
			ss >> shopIncome;
			dt->SetPoint(n, products, meanResidenceTime, shopIncome);
			n++;
	}
	gStyle->SetPalette(1);
	dt->SetTitle("Zaleznosc zysku od liczby produktow oraz czasu przebywania w sklepie");
	dt->GetXaxis()->SetTitle("Liczba produktow");
	dt->GetYaxis()->SetTitle("Sredni czas przebywania w sklepie");
	dt->GetZaxis()->SetTitle("Dochod");
	dt->Draw("surf1");
}

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

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


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


//    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);
   }
   // fit the graph now 
   
   ROOT::Fit::Fitter  fitter;
   
   
   // make the functor objet
   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();
   if (!ok) {
      Error("line3Dfit","Line3D Fit failed");
      return 1;
   }
   
   const ROOT::Fit::FitResult & result = fitter.Result();
   
   std::cout << "Total final distance square " << result.MinFcnValue() << std::endl;
   result.Print(std::cout);
   

   gr->Draw("p0");

   // get fit parameters
   const double * parFit = result.GetParams();

   
   // draw the fitted line
   int n = 1000;
   double t0 = 0;
   double dt = 10;
   TPolyLine3D *l = new TPolyLine3D(n);
   for (int i = 0; i <n;++i) {
      double t = t0+ dt*i/n;
      double x,y,z;
      line(t,parFit,x,y,z);
      l->SetPoint(i,x,y,z);
   }
   l->SetLineColor(kRed);
   l->Draw("same");
   
   // draw original line
   TPolyLine3D *l0 = new TPolyLine3D(n);
   for (int i = 0; i <n;++i) {
      double t = t0+ dt*i/n;
      double x,y,z;
      line(t,p0,x,y,z);
      l0->SetPoint(i,x,y,z);
   }
   l0->SetLineColor(kBlue);
   l0->Draw("same");
   return 0;
}

//.........这里部分代码省略.........
    text->DrawLatex(-3.,0.7,"#kappa_{#lambda} = -2.4, c_{g}  = c_{2} = 0");  
    c1->cd(11);
    //c1_5->SetTicks(0,0);
    //c1_2->SetRightMargin(0.15);
    //c1_2->SetLeftMargin(0.15);
    //c1_2->SetBottomMargin(0.02);
    lm5->Draw("colz1");
    text->DrawLatex(-3.,0.7,"#kappa_{#lambda} = -5, c_{g}  = c_{2} = 0");  
    c1->cd(12);
    //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]);    

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

	po::options_description desc("Allowed Options");
	desc.add_options()
		 ("help,h", "Produce this help message")
		 ("startwl,s",po::value<double>(),"Set the start Wavelength for the Analysis")
		 ("stopwl,p",po::value<double>(),"Set the stop Wavelength for the Analysis")
		 ("non-interactive,n","Runs the program in Noninteractive mode. It quits when it's finished")
		 ("version,v","Prints Version")
	;
		 
	po::variables_map vm;
		 po::store(po::parse_command_line(argc,argv,desc),vm);
		 po::notify(vm);
	if (vm.count("help")) {
		std::cout << desc<< std::endl;
		return 3;
	}
	if (vm.count("version")) {
		std::cout << "VCSEL Laser Analysis Version " << _VERSION << std::endl;
		std::cout << "Using ROOT version " << _ROOT_VERSION << " and Boost version " << _BOOST_VERSION << std::endl;
		return 0;
	}
	
	if (argc < 4) {
		std::cout << desc;
		return 2;
	}	
	double startwl, stopwl;
	startwl = 842.;
	stopwl = 860.;
	bool run = true;
	if (vm.count("startwl")) {
		startwl = vm["startwl"].as<double>();
		NUM_ARGS +=2;
	}
	if (vm.count("stopwl")) {
		double tmp =  vm["stopwl"].as<double>();
		stopwl =tmp;
		NUM_ARGS +=2;
	}
	if (vm.count("non-interactive")) {
		run = false;
		NUM_ARGS++;
	}
	
	
	//checking filetypes must be txt, csv or CSV
	if (!check_extensions(argc, argv)) {
		return 1;
	}
	std::cout <<"startwl: "<< startwl << '\t' << "stopwl: " << stopwl << std::endl;
	
	Double_t max = -210;
	Double_t maxwl = 0;
	int _argc = argc;
	TApplication *t = new TApplication("big",&_argc,argv);
	std::cout << "Running with boost and ROOT" <<std::endl;
	std::vector<double> _x,_y;
	Double_t x[LINES], y[LINES], _inta[LINES], _intb[LINES]; 
	
	Double_t *cmp_int = new Double_t[argc];
	Double_t *argc_ary = new Double_t[argc];
	Double_t *cmp_int_root = new Double_t[argc];
	Double_t *asymmety_ary = new Double_t[argc];
	Double_t *width_ary = new Double_t [argc];
	

	
	TGraph2D *gr = new TGraph2D(LINES*(argc-1));
		//Setting up canvas for plot of all sectrums (is it called spectrums? ;) )
	TCanvas *c1 = new TCanvas("All Plots","All Plots",10,10,3000,1500);
	TH1F *integral_hist = new TH1F("Asymmerty", "Asymmetry", 100,0, 100);
	

	if(!(argc % ROWS)){
		c1->Divide(argc/ROWS,ROWS);
		
	}else{
		c1->Divide(argc/ROWS+(argc %ROWS -1),ROWS);
	}
	
	for (Int_t i = NUM_ARGS +1; i < argc ; i++){
		try{ 
			
			max = -211;
			maxwl = 0;
			argc_ary[i] = i-NUM_ARGS;
			
			std::ifstream in;
			in.seekg(0, std::ios::beg);
				// voodoo keep this;
			char **arg1 = t->Argv() ;
			std::string tmp = arg1[i];
			in.open(tmp.c_str());
			
			std::cout<< "file: " << tmp << std::endl;
			std::string line;
//.........这里部分代码省略.........

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");
}

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 CommandMSUGRA(TString plotName_,Int_t tanBeta_, Bool_t plotLO_){
  gROOT->SetStyle("CMS");//jmt specific
  gROOT->ForceStyle();

  gStyle->SetOptTitle(0);
  gStyle->SetOptStat(0);
  gStyle->SetPalette(1); 
  gStyle->SetTextFont(42);
  gStyle->SetFrameBorderMode(0);

  //convert tanb value to string
  std::stringstream tmp;
  tmp << tanBeta_;
  TString tanb( tmp.str() );
  
  
  // Output file
  std::cout << " create " << plotName_ << std::endl;
  TFile* output = new TFile( plotName_, "RECREATE" );
  if ( !output || output->IsZombie() ) { std::cout << " zombie alarm output is a zombie " << std::endl; }
  

  //set old exclusion Limits
  TGraph* LEP_ch = set_lep_ch(tanBeta_);
  TGraph* LEP_sl = set_lep_sl(tanBeta_);//slepton curve
  TGraph* TEV_sg_cdf = set_tev_sg_cdf(tanBeta_);//squark gluino cdf
  TGraph* TEV_sg_d0 = set_tev_sg_d0(tanBeta_);//squark gluino d0
  //  TGraph* TEV_tlp_cdf = set_tev_tlp_cdf(tanBeta_);//trilepton cdf
  //  TGraph* TEV_tlp_d0 = set_tev_tlp_d0(tanBeta_);//trilepton d0
  TGraph* stau   = set_tev_stau(tanBeta_);//stau 
  TGraph* NoEWSB = set_NoEWSB(tanBeta_); 

  TGraph* TEV_sn_d0_1 = set_sneutrino_d0_1(tanBeta_);
  TGraph* TEV_sn_d0_2 = set_sneutrino_d0_2(tanBeta_);

  int nPoints = nSusyGridPoints();
  double m0[nPoints],m12[nPoints],squarkMass[nPoints],gluinoMass[nPoints];

  susyGrid(m0,m12,squarkMass,gluinoMass);

  TGraph2D* squarkMasses = new TGraph2D("squarkMasses","",nPoints,m0,m12,squarkMass);
  TGraph2D* gluinoMasses = new TGraph2D("gluinoMasses","",nPoints,m0,m12,gluinoMass);

  TH2D* gluinoMassPlot = gluinoMasses->GetHistogram();
  TH2D* squarkMassPlot = squarkMasses->GetHistogram();

  //constant ssqquark and gluino lines
  TF1* lnsq[15];
  TF1* lngl[15];

  TGraph* lnsq_40[15];
  TGraph* lngl_40[15];
  
  TLatex* sq_text[15];
  TLatex* gl_text[15];

  TLatex* sq_40_text[15];
  TLatex* gl_40_text[15];

  for(int i = 1; i < 15; i++){
    //lnsq[i] = constant_squark(tanBeta_,i);
    //sq_text[i] = constant_squark_text(i,*lnsq[i],tanBeta_);
    //lngl[i] = constant_gluino(tanBeta_,i);
    //gl_text[i] = constant_gluino_text(i,*lngl[i]);
    lnsq_40[i] = constant_mass(i*250,squarkMasses);
    lngl_40[i] = constant_mass(i*250,gluinoMasses);
    sq_40_text[i] = constant_squark_text_tanBeta40(i*250,lnsq_40[i]);
    gl_40_text[i] = constant_gluino_text_tanBeta40(i*250,lngl_40[i]);;
  }


  //Legends
  TLegend* legst  = makeStauLegend(0.05,tanBeta_);
  TLegend* legNoEWSB  = makeNoEWSBLegend(0.05,tanBeta_);
  TLegend* legexp = makeExpLegend( *TEV_sg_cdf,*TEV_sg_d0,*LEP_ch,*LEP_sl,*TEV_sn_d0_1,0.035,tanBeta_);
  
 
  //make Canvas
  TCanvas* cvsSys = new TCanvas("cvsnm","cvsnm",0,0,800,600);
  gStyle->SetOptTitle(0);
  cvsSys->SetFillColor(0);
  cvsSys->GetPad(0)->SetRightMargin(0.07);
  cvsSys->Range(-120.5298,26.16437,736.0927,750);
  //  cvsSys->Range(-50.5298,26.16437,736.0927,500);
  cvsSys->SetFillColor(0);
  cvsSys->SetBorderMode(0);
  cvsSys->GetPad(0)->SetBorderMode(0);
  cvsSys->GetPad(0)->SetBorderSize(2);
  cvsSys->GetPad(0)->SetLeftMargin(0.1407035);
  cvsSys->GetPad(0)->SetTopMargin(0.08);
  cvsSys->GetPad(0)->SetBottomMargin(0.13);

  cvsSys->SetTitle("tan#beta="+tanb);
 
  output->cd();
  
//and now  the exclusion limits


  TGra