C++ TGraphErrors::SetLineWidth方法代码示例

void ExpManager::GetSim1DGraphZ(TString NameTitle, double xmin, double xmax, double ymin, double ymax, TString grid ){   // do the same thing as for the polar interpolation

    TGraphErrors *fGraph = new TGraphErrors(); //= new TGraph2DErrors(np, x_array, y_array, bz_array, ex, ey, ez);
    fGraph->SetTitle("Simulated Data;Z (mm);Magnetic Field (mT)");
    fGraph->SetMarkerSize(1.2);
    fGraph->SetMarkerStyle(20);
    fGraph->SetMarkerColor(kBlue);
    fGraph->SetLineColor(kBlue);
    fGraph->SetLineWidth(2);
       
    int graph_counter = 0 ; 
   for (unsigned i=0; i< fExpY.size(); i++)   {
    //cout <<  " X "  << fExpX.at(i) ;  
        if( (fExpX.at(i) >= xmin && fExpX.at(i) <= xmax) && (fExpY.at(i) >= ymin && fExpY.at(i) <= ymax) &&  fGrid.at(i)==grid ){
            fGraph->SetPoint(graph_counter,fExpZ.at(i),fSimB.at(i));    
            fGraph->SetPointError(graph_counter,fExpZErr.at(i),0);   
            graph_counter++;
        } 
   }

    fGraph->SetTitle(NameTitle+Form(" Simulated Data : %.2f < X < %.2f mm  __  %.2f < Y < %.2f mm;Z (mm);Magnetic Field (mT)",xmin,xmax,ymin,ymax));
    fGraph->SetName(NameTitle+Form("_Sim_X_%.2f_%.2fmm_Y_%.2f_%.2fmm",xmin,xmax,ymin,ymax));
    fGraph->Write();

}

//______________________________________________________________________________
// Construct a graph from vectors and y error vector
TGraphErrors *LoadGraphFromVectorsWithError(vector<double> xVector, vector<double> yVector, vector<double> yErrorVector, string xTitle, string yTitle)
{
	int n = xVector.size();

	if ((xVector.size() == yVector.size()) &&
		(yVector.size() == yErrorVector.size()))
	{
		//Create a graph
		TGraphErrors *gr = new TGraphErrors(n, &xVector[0], &yVector[0], 0, &yErrorVector[0]);
		gr->SetTitle("");
		gr->SetMarkerStyle(20);
		gr->SetMarkerSize(1.2);
		gr->SetLineWidth(2);
		gr->GetXaxis()->SetTitle(xTitle.c_str());
		gr->GetXaxis()->CenterTitle();
		gr->GetYaxis()->SetTitle(yTitle.c_str());
		gr->GetYaxis()->CenterTitle();
		return gr;
		delete gr;
	}
	else
	{
		TGraphErrors *gr0 = new TGraphErrors();
		return gr0;
		delete gr0;
	}
}

void plot_mc(string name="g", double q2=1.9, Int_t n=0)
{
//    SetAtlasStyle();
    TGraphErrors *p;

    TString nn;
    nn.Form("ave_%s_vs_x_for_Q2=%g",name.c_str(),q2);
    gDirectory->GetObject(nn,p);

//    p->Print();

    Double_t ratsize = 0.0;

    TCanvas *c = new TCanvas("PDF","pdf",600,600);
    TPad *pad1 = new TPad("pad1","pad1",0.,ratsize,1.,1.);
    pad1->SetLogx();
    pad1->Draw();
    pad1->cd();
    p->GetXaxis()->Set(101,0.0001,1.);
    p->SetFillColor(kRed-2);
    p->SetFillStyle(3001);
    p->SetLineWidth(1);
    p->SetLineColor(kRed);
    p->GetYaxis()->SetTitle(name.c_str());

    p->GetYaxis()->SetTitleSize(0.06);
    p->GetYaxis()->SetTitleOffset(1.);
    p->GetXaxis()->Set(101,0.0001,1.);
    p->Draw("ALE3");

    Double_t av = 0;
    Double_t av2 = 0;

    for (Int_t i = 1; i<n+1 ; i++) {
        nn.Form("%s_vs_x_for_Q^{2}=%g_%i",name.c_str(),q2,i);
        TGraph *pp = NULL;
        gDirectory->GetObject(nn,pp);
        if (pp != NULL) {
            pp->SetLineColor(kGreen);
            pp->Draw("L");

            av  += pp->GetY()[0];
            av2 += pp->GetY()[0]*pp->GetY()[0];
            cout << i << " "<<pp->GetY()[0] << endl;


        }
    }

    av  /= n;
    av2 = sqrt(av2/n - av*av);
//    cout << n << " "<<av << " "<< av2<<endl;


    p->Draw("E3C");


}

void ExpManager::DrawMap(TString NameTitle, double xmin, double xmax, double ymin, double ymax, double zmin, double zmax) {

    TMultiGraph *mg = new TMultiGraph();
    
    //Draw a cross
    TGraphErrors *frame = new TGraphErrors(); //= new TGraph2DErrors(np, x_array, y_array, bz_array, ex, ey, ez);
    frame->SetPoint(0,+0,+100);
    frame->SetPoint(1,+0,-100);
    frame->SetPoint(2,0,0); 
    frame->SetPoint(3,-100,0); 
    frame->SetPoint(4,+100,0);
    frame->SetPoint(5,+0,+0); 
    frame->SetPoint(6,+0,-100);
    frame->SetMarkerColor(kWhite);
    frame->SetDrawOption("ap");
    
    //Draw the map 
    TGraphErrors *fGraph = new TGraphErrors(); //= new TGraph2DErrors(np, x_array, y_array, bz_array, ex, ey, ez);
    fGraph->SetMarkerSize(1.2);
    fGraph->SetMarkerStyle(20);
    fGraph->SetMarkerColor(kBlue);
    fGraph->SetLineColor(kBlue);
    fGraph->SetLineWidth(1);
    //fGraph->SetDrawOption("ap");
    
    int graph_counter = 0 ; 
   for (unsigned i=0; i< fExpY.size(); i++)   {
        if( (fExpX.at(i) >= xmin && fExpX.at(i) <= xmax) && (fExpY.at(i) >= ymin && fExpY.at(i) <= ymax) && (fExpZ.at(i) >= zmin && fExpZ.at(i) <= zmax) ) {
            fGraph->SetPoint(graph_counter,fExpX.at(i),fExpY.at(i));    
            fGraph->SetPointError(graph_counter,fExpXErr.at(i),fExpYErr.at(i)); 
            graph_counter++;
        } 
   }

    //fGraph->Write();
    mg->SetTitle(NameTitle+Form(" Map : %.2f < X < %.2f mm    %.2f < Y < %.2f mm    %.2f < Z < %.2f mm;X (mm);Y (mm)",xmin,xmax,ymin,ymax,zmin,zmax));
    mg->SetName(NameTitle+Form("_Map__X_%.2f_%.2fmm__Y_%.2f_%.2fmm__Z_%.2f_%.2fmm",xmin,xmax,ymin,ymax,zmin,zmax));
    mg->Add(frame); 
    mg->Add(fGraph);
    mg->Write();
}

void PlotParValVsLabelWithErr(TFile* f, TTree* tr, TString strMillepedeRes, TString strOutdir)
{
 
  f->cd();
  TString canvName="c_";
  canvName+=strMillepedeRes;
  canvName+="_";
  canvName+=StrPlotType(PARSwithERRvsLABEL);
  canvName.ReplaceAll(".res","");

  TCanvas* canv = new TCanvas(canvName,canvName,900,600);
  canv->Divide(3,2);

  for (int ind=1; ind<=6; ind++){
    canv->cd(ind);
    TPad* pad = (TPad*)canv->GetPad(ind);
    TString strCut="((label%20-1)%9+1)==";
    strCut+=ind;
    int n = tr->Draw("label%700000:10000*parVal:10000*parErr:0.01*(label%700000)",strCut,"goff");
    TGraphErrors *gr = new TGraphErrors(n,tr->GetV1(),tr->GetV2(),tr->GetV4(),tr->GetV3());
    gr->SetMarkerStyle(20);
    gr->SetLineWidth(2);
    for (int i=0; i<n; i++){
      std::cout<<tr->GetV1()[i]<<" "<<tr->GetV2()[i]<<"+-"<<tr->GetV3()[i]<<std::endl;
    }
    gr->SetTitle(StrPar(ind)+TString(", 10000*(par+-err)"));
    gr->GetXaxis()->SetTitle("label%700000");
    gr->Draw("AP");
  }// end of loop over ind
  canvName+=".png";
  TString saveName=strOutdir+canvName;
  canv->SaveAs(saveName);
  saveName.ReplaceAll(".png",".pdf");
  canv->SaveAs(saveName);

}// end of PlotParValVsLabelWithErr

void getY0()
{
  gROOT->Reset();
  gROOT->LoadMacro("fitfun.C");   // fit function for localY, local Z

  const Int_t N = 120;

  //histograms to check the par2 and par3 distributions
  TH1D *hPar2 = new TH1D("Par2","",500,0., 5.);
  TH1D *hPar3 = new TH1D("Par3","",500,0., 1.);

  TFile *fin = new TFile("align.root");
  TH2D *yLocal2D = (TH2D *)fin->Get("yLocal_all");
  yLocal2D->RebinY(5);
  TH1D *yLocal[N];
  for(int i=0;i<N;i++) {
    char hisname[100];
    sprintf(hisname,"Tray_%d",i+1);
    yLocal[i] = yLocal2D->ProjectionY(hisname,i+1,i+1);
    yLocal[i]->Sumw2();
  }

  TF1 *fitfun = new TF1("fitfun",fitfun,-5.,5.,5);

  double t[N];
  double y0[N], ye[N];

  gStyle->SetOptStat(101);
  gStyle->SetOptFit(100);
  gStyle->SetTextSize(0.06);
  gStyle->SetTextFont(42);
  TCanvas *c1 = new TCanvas("c1","c1",0,0,800,600);
  c1->SetFillColor(10);
  c1->SetBorderMode(0);
  c1->SetBorderSize(2);
  c1->Divide(4,3);

  TPostScript *ps = new TPostScript("fig/yLocalFit.ps",112);

  for(int i=0;i<N;i++) {
    t[i] = i+1;  // tray number

    if(i%12==0) ps->NewPage();   // 12 trays per page
    c1->cd(i%12+1);

    double par[5], err[5];
    fitfun->SetParameters(0., yLocal[i]->GetBinContent(50), 0.1, 1., 0.);
    fitfun->SetLineWidth(2);
    yLocal[i]->Draw("e");
    double entries = yLocal[i]->Integral(1,100);
    if(entries>10) {
      fitfun->SetParLimits(2, 1.55, 1.7);
      fitfun->SetParLimits(3, 0.1, 0.3);
      yLocal[i]->Fit("fitfun","R");
      fitfun->GetParameters(&par[0]);
      err[0] = fitfun->GetParError(0);
    } else {
      par[0] = 0.;
      err[0] = 0.;
    }

    hPar2->Fill(par[2]);
    hPar3->Fill(par[3]);

    char text[100];
    sprintf(text, "Y0 = %5.3f #pm %5.3f", par[0], err[0]);

    TLatex *tex = new TLatex(-3, par[1]/3., text);
    tex->SetTextSize(0.07);
    tex->SetTextFont(12);
    tex->Draw("same");

    y0[i] = par[0];
    ye[i] = err[0];
    
    c1->Update();
  }

  ps->NewPage();

  gStyle->SetOptStat(0);
  gStyle->SetOptFit(0);
  gStyle->SetEndErrorSize(0.01);
   
  TCanvas *c1 = new TCanvas("c1","c1",0,0,800,600);
  c1->SetFillColor(10);
  c1->SetBorderMode(0);
  c1->SetBorderSize(2);
  c1->SetLeftMargin(0.14);
  c1->SetBottomMargin(0.15);
  c1->SetTopMargin(0.02);
  c1->SetRightMargin(0.04);
  c1->Draw();

  double x1 = 0.5;
  double x2 = 120.5;
  double y1 = -2.;
  double y2 = 2.;
  TH1D *h0 = new TH1D("h0","",1, x1, x2);
  h0->SetMinimum(y1);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    can->cd(iw+1);
    if(iw==3||iw==7||iw==11) gPad->SetRightMargin(0.02);
    char ch1[8];
    sprintf(ch1,"h1_%d",iw);
    h1[iw] = new TH1D(ch1,"",500,hxmin[iw],hxmax[iw]);
    h1[iw]->SetMinimum(hymin[iw]); h1[iw]->SetMaximum(hymax[iw]);
    h1[iw]->SetXTitle(atit[iw][0]); h1[iw]->SetYTitle(atit[iw][1]);
    h1[iw]->GetXaxis()->CenterTitle(1);
    h1[iw]->GetYaxis()->CenterTitle(1);

    // futz with the axes
    h1[iw]->GetYaxis()->SetNdivisions(606);
    h1[iw]->GetXaxis()->SetNdivisions(606);

      h1[iw]->GetYaxis()->SetTitleSize(textsize);
      h1[iw]->GetYaxis()->SetTitleOffset(3.5);
      h1[iw]->GetYaxis()->SetLabelSize(textsize);
      h1[iw]->GetXaxis()->SetTitleSize(textsize);
      h1[iw]->GetXaxis()->SetTitleOffset(3.0);
      h1[iw]->GetXaxis()->SetLabelSize(textsize);
     // h1[iw]->GetXaxis()->SetLabelOffset(1.2);

    h1[iw]->Draw();
  }
  //
  // Draw! 
  //
  cout << "Now Draw!" << endl;
  TGraphErrors *ge;
  for (int iw=0;iw<nw;iw++)
  {
    can->cd(iw+1);

    //
    if (iw==0)
    {
      tex=new TLatex(0.7,0.05,"HIN11005 is 2.76TeV, EP");
      tex->SetTextSize(textsize*0.9);
      //tex->Draw();
      tex=new TLatex(1.5,0.07,"|#eta|<0.8");
      tex->SetTextSize(textsize);
      //tex->Draw();
    }
    if (iw < 4) tex=new TLatex(ptit[iw][0]*1.2,ptit[iw][1]*1.27,htit[iw]);
    else tex=new TLatex(ptit[iw][0]*1.2,ptit[iw][1]*1.27,htit[iw]);
    if ( iw == 0)   tex->SetTextSize(textsize);
    else if (iw == 8) tex->SetTextSize(textsize);
    else  tex->SetTextSize(textsize);
    tex->Draw();
    //
    // Legend
    //
    TLegend *leg = new TLegend(lxmin[iw]*3.99,lymin[iw]*1.1,lxmax[iw]*1.99,lymax[iw]*1.2);
    leg->SetFillColor(0);
    leg->SetTextSize(textsize);
    //
    // Read data
    //
    for (int im=0;im<1;im++)
    {
      int j=im*12+iw;
      int col=2;
      //if (ndp[j]==0) continue;
      ge=new TGraphErrors(16,&mxa[j][0],&mya[j][0],&mxe[j][0],&mye[j][0]);
      drawSysBoxValue(ge,16,0.2, myesys[j]);
      ge->SetTitle("");
      ge->SetMarkerStyle(25);
      //      ge->SetMarkerStyle(msty[im]);
      ge->SetMarkerSize(msiz[im]);
      ge->SetMarkerColor(col);
      ge->SetLineWidth(1.2);
      ge->SetLineColor(col);
      ge->Draw("pe");
      leg->AddEntry(ge,"v_{5}{#Psi_{23}} Steve","p");
    }
    for (int im=0;im<1;im++)
    {
      int j=im*12+iw;
      int col=1;
      //if (ndp[j]==0) continue;
      ge=new TGraphErrors(19,&xa[j][0],&ya[j][0],&xe[j][0],&ye[j][0]);
      //drawSysBoxValue(ge,16,0.2, yesys[j]);
      ge->SetTitle("");
      ge->SetMarkerStyle(msty[im][1]);
      //      ge->SetMarkerStyle(msty[im]);
      ge->SetMarkerSize(msiz[im]);
      ge->SetMarkerColor(col);
      ge->SetLineWidth(1.2);
      ge->SetLineColor(col);
      ge->Draw("pe");
      leg->AddEntry(ge,"v_{5}{#Psi_{23}} Shengquan","p");
    }
    if (iw==0) leg->Draw();

  }
  cout << "end of process" << endl;
  can->cd();
  can->Print("./figures/v523_pt_ShengquanSteveb.png");
  can->Print("./figures/v523_pt_ShengquanSteveb.pdf");
}

//.........这里部分代码省略.........
	mgClipOutputFWHM->GetXaxis()->SetTitle("V clipping (mV)");
	mgClipOutputFWHM->GetYaxis()->SetTitle("Output FWHM (ns)");
	cmgClipOutputFWHM->Update();
		
	legend = new TLegend(0.6,0.67,0.89,0.86,"V input");
	legend->AddEntry(gClip1Point9V, "1.9 V", "lp");
	legend->AddEntry(gClip960mV, "960 mV", "lp");
	legend->AddEntry(gExpClip960mV, "Exp 960 mV", "lp");
	legend->AddEntry(gExpClip1Point9V, "Exp 1.9 V", "lp");
	legend->SetTextSize(0.04);
	legend->SetMargin(0.5);
	legend->Draw();
	
	// Hysteresis plot: V output (t) in function of V input (t) for several clipping values
	
	// variables used in the analysis
	
	Long64_t iInputMax, iOutputMax;
	Float_t xInputMax, xOutputMax, xInputHalf, xOutputHalf;
	Double_t InputMax, OutputMax, InputHalf, OutputHalf;
	
	Long64_t firstIndex = 0;
	Long64_t lastIndex = 0;
	Long64_t inputGraphPoints = 0;
	Long64_t outputGraphPoints = 0;
	
	// hard coded values to cut the x axis of both waves
	// Input wave
	
	inputWaveName += " graph";
	
	TGraphErrors *gInput = listOfGraphs->FindObject(inputWaveName);
	gInput->SetLineColor(kRed);
	gInput->SetLineWidth(2);
	
	xInput = gInput->GetX();
	yInput = gInput->GetY();
	inputGraphPoints = gInput->GetN();
	
	cout << inputGraphPoints << endl;
	
	// Invert the input wave
	
	for(Int_t i = 0; i < inputGraphPoints; i++) {
		yInput[i] = -(yInput[i]);
	}
	
	// find the x at which the graph reaches the max value
	
	iInputMax = TMath::LocMax(inputGraphPoints, yInput);
	xInputMax = xInput[iInputMax];
	cout << "iInputMax = " << iInputMax << endl;
	cout << "xInputMax = " << xInputMax << endl;
	InputMax = gInput->Eval(xInput[iInputMax]);
	cout << "InputMax = " << InputMax << endl;
	
	// Output wave
	
	outputWaveName += " graph";
	
	TGraphErrors *gOutput = listOfGraphs->FindObject(outputWaveName);
	gOutput->SetLineWidth(2);
	
	xOutput = gOutput->GetX();
	yOutput = gOutput->GetY();
	outputGraphPoints = gOutput->GetN();

//.........这里部分代码省略.........
  h0->GetXaxis()->SetLabelSize(0.04);
  h0->GetXaxis()->SetLabelFont(42);
  h0->GetXaxis()->SetTitleFont(42);
  h0->GetYaxis()->SetNdivisions(505);
  h0->GetYaxis()->CenterTitle();
  h0->GetYaxis()->SetTitle("Anisotropy Parameter, v_{2}");
  h0->GetYaxis()->SetTitleOffset(.9);
  h0->GetYaxis()->SetTitleSize(0.05);
  h0->GetYaxis()->SetLabelOffset(0.005);
  h0->GetYaxis()->SetLabelSize(0.04);
  h0->GetYaxis()->SetLabelFont(42);
  h0->GetYaxis()->SetTitleFont(42);
  h0->Draw("c");

  
  double pt[2] = {1.766, 4.459};
  double pt1Bin[1] = {1.866};
  double v2[2] = {0.08329931, 0.126231279};
  double eptLow[2] = {.766, .459};
  double eptHigh[2] = {2.234, 2.341};
  double eptSyst[2] = { 0.0, 0.0};
  double eptSyst1Bin[1] = { 0.0};
  double ev2Low[2] = { 0.00736501, 0.012241649};
  double ev2High[2] = { 0.00736501, 0.012241649};
  double ev2Syst[2] = {.0124948965, .02145931743 };
  
  
  TGraphAsymmErrors *AlexData = new TGraphAsymmErrors(2, pt, v2, eptLow, eptHigh, ev2Low, ev2High);
  
  AlexData->SetMarkerColor(2);
  AlexData->SetMarkerStyle(21);
  AlexData->SetMarkerSize(2);
  AlexData->SetLineColor(2);
  AlexData->SetLineWidth(3);
  
  TGraphErrors *AlexDataSyst = new TGraphErrors(2, pt, v2, eptSyst, ev2Syst);
  
  AlexDataSyst->SetMarkerColor(2);
  AlexDataSyst->SetMarkerStyle(21);
  AlexDataSyst->SetMarkerSize(1.1);
  AlexDataSyst->SetLineColor(2);
  AlexDataSyst->SetLineWidth(33);
  
  
  
  // 8 data points
   
 
  
  double weights[6] = {0.4034, .3058, .1696, .0778, .0316, .0119};

  double avgV2[1] = {0};
  double avgV2StatErrorSquared[1] = {0};
  double avgV2SystErrorSquared[1] = {0};
  double avgV2StatError[1] = {0};
  double avgV2SystError[1] = {0};
  
  for(int i = 0; i < 6; i++){

    cout << "V2: " << y_data_cen[1][i] << endl;
    
    avgV2[0] = avgV2[0] + (weights[i]*y_data_cen[1][i]);
    
    avgV2StatErrorSquared[0] = avgV2StatErrorSquared[0] + (weights[i]*weights[i])*(ye_data_cen[1][i]*ye_data_cen[1][i]);
    avgV2SystErrorSquared[0] = avgV2SystErrorSquared[0] + (weights[i]*weights[i])*(yes_data_cen[1][i]*yes_data_cen[1][i]);
    

//.........这里部分代码省略.........
      stringstream filePath;
      filePath << "pmtratestudy/run" << data[0] << "*.root";
      cout << "opening file at " << filePath.str() << endl;
      cout << "file counter: " << fileCounter << " channel=" << ch << endl;
      //TFile* f = new TFile(filePath.str().c_str());
      //TTree* t = (TTree *)f->Get("eventtree");
      TChain* t = new TChain("eventtree");
      t->Add( filePath.str().c_str() );
      out->cd();

      x[fileCounter] = data[1];

      int nfires[NCH] = {0};
      int samples = 0;
      float chmax = 0.0;
      t->SetBranchAddress("nfires", &nfires);
      t->SetBranchAddress("samples", &samples);
      t->SetBranchAddress("chmax", &chmax);
      
      h->Reset();
      
      int nentries = t->GetEntries();
      for (int entry = 0; entry < nentries; ++entry) {
        t->GetEntry(entry);
        if (chmax < 100.0) {
          h->Fill(nfires[ch]);
        }
      }


      
      pois->SetParameter(0,1);
      pois->SetParameter(1, h->GetMean());

      h->Fit(pois,"RQ","",0,50);
      //TF1 *myfit = (TF1 *)h->GetFunction("pois");
      TF1 *myfit = (TF1 *)pois;
      Double_t lambda = myfit->GetParameter(1);  
      h->Draw();
      stringstream histFileName;
      histFileName << "hist/h" << data[0] << "_ch" << ch << ".png";
      c1->SaveAs(histFileName.str().c_str());
      //Graph with poisson method
#if 1
      y[fileCounter] = lambda / ((samples - 1) * 15.625E-6);
      errY[fileCounter] = myfit->GetParError(1) / ((samples - 1) * 15.625E-6);
#endif
      //Graph with mean method
#if 0
      y[fileCounter] = h->GetMean() / ((samples - 1) * 15.625E-6);
      errY[fileCounter] = h->GetMeanError() / ((samples - 1) * 15.625E-6);
#endif
      cout << x[fileCounter] << ", " << y[fileCounter] 
           << " | " << (samples - 1) << endl;
      delete t;
      //f->Close();
      fileCounter++;
    } 

    ppp->SetParameter(0,1);
    ppp->SetParameter(1,0.4);
     
    TGraphErrors* gr = new TGraphErrors(NRUNS, x, y, errX, errY);
    gr->SetLineColor(color[ch % NCOLORS]);
    cout << "color: " << color[ch % NCOLORS] << endl;
    gr->SetLineWidth(2);
    gr->SetMarkerStyle(7);
    gr->Fit("ppp","R0","Q0",0.045,2.0);
    TF1 *afit = (TF1 *)gr->GetFunction("ppp");
    Double_t aRate = 1/afit->GetParameter(1);  
    if (aRate > 0) {
      hRate->Fill(aRate);
    }
    gr->GetXaxis()->SetTitle("Run Date");
    gr->GetYaxis()->SetTitle("Rate [kHz]");

    stringstream entryName, fileName;
    entryName << "Channel" << ch;
    gr->SetTitle(entryName.str().c_str());
    fileName << "plots/" << ch << ".png";
    legend->AddEntry(gr, entryName.str().c_str());
    gr->Draw("alp");
    c1->SaveAs(fileName.str().c_str());
    mg->Add(gr);
    cout << "added plot to mg\n";
    fin.clear();
    fin.seekg(0, ios::beg);
  } // loop over channel
  hRate->Draw();
  hRate->Fit("gaus");
  c1->SaveAs("hrate.pdf");
  mg->Draw("alp");
  mg->GetXaxis()->SetTitle("Days since first run");
  mg->GetYaxis()->SetTitle("Rate [kHz]");
  mg->SetTitle("All channels: Rate vs. Days since first Run");

  legend->Draw();
  c1->SaveAs("mg.pdf");
  
}

void make1DLimitHH(TString combine_dir,bool blind=true){
  //TString combine_dir = "test_runSusyHgg/signalInj_sms_ChiHH_0_175/";

  TGraph obser( (275-125)/25 );
  TGraph graph( (275-125)/25 );
  TGraphAsymmErrors error( (275-125)/25 );

    for(int m=125;m<501;m+=25) {
      int i=(m-125)/25;
      TFile limit_file(Form("%s/higgsCombineChiHH_0_%d.Asymptotic.mH120.root",combine_dir.Data(),m) );
      TTree *limit_tree = (TTree*)limit_file.Get("limit");
      TTreeFormula limit_form("get_limit","limit",limit_tree);

      limit_tree->GetEntry(1);
      float down = limit_form.EvalInstance();
      limit_tree->GetEntry(2);
      float exp = limit_form.EvalInstance();
      limit_tree->GetEntry(3);
      float up = limit_form.EvalInstance();
      limit_tree->GetEntry(5);
      float obs = limit_form.EvalInstance();

      graph.SetPoint(i,float(m),exp);
      error.SetPoint(i,float(m),exp);
      error.SetPointError(i,0,0,exp-down,up-exp);

      obser.SetPoint(i,float(m),obs);
    }

    TGraphErrors* theo = getTheoXSec("/home/amott/HggApp/SusyHgg/xsecs/Higgsino.txt");

    TCanvas cv;
    cv.SetLogy();

    theo->SetMaximum(1e2);
    theo->SetMinimum(1e-2);
    theo->GetYaxis()->SetTitle("95% CL #sigma upper limit (pb)");
    theo->GetXaxis()->SetTitle("m_{chargino}");
    theo->SetFillColor(kBlue);
    theo->SetLineStyle(kDotted);
    theo->SetLineWidth(2.0);
    
    error.SetMaximum(1e2);
    error.SetMinimum(1e-2);
    error.GetYaxis()->SetTitle("95% CL #sigma upper limit (pb)");
    error.GetXaxis()->SetTitle("m_{chargino}");
    error.SetFillColor(kGreen);
    error.Draw("A3");

    theo->Draw("3C");

    graph.SetLineStyle(kDashed);
    graph.SetLineWidth(2);
    graph.Draw("C");

    obser.SetLineStyle(1);
    obser.SetLineWidth(2);
    if(!blind) obser.Draw("C");

    TLegend leg(0.7,0.7,0.85,0.85);
    leg.SetFillColor(0);
    leg.SetBorderSize(0);
    leg.AddEntry(&graph,"expected","l");
    leg.AddEntry(&error,"expected #pm1#sigma","F");
    leg.AddEntry(theo,"theoretical","f");
    if(!blind)     leg.AddEntry(&obser,"observed","l");

    leg.Draw("SAME");

    TLatex prelim(0.65,0.96,"CMS Preliminary");
    prelim.SetNDC();
    prelim.SetTextSize(0.045);
    prelim.Draw();

    TLatex lbl(0.5,0.86,"#sqrt{s} = 8 TeV  #int L dt = 19.78 fb^{-1}");
    lbl.SetNDC();
    lbl.SetTextSize(0.045);
    lbl.Draw();

    cv.SaveAs(combine_dir+"expected_exclusion_HH_1D.png");


}

//.........这里部分代码省略.........
		
    fitFunc->SetParameter(6, 0.05*hNsigma->GetMaximum());
    fitFunc->SetParameter(7, 4.); //was at 5. but it seems like 4. works better
    fitFunc->SetParameter(8, 3.);
		
		hNsigma->Fit("fitFunc");
    
    float  ampPi = fitFunc->GetParameter(0);
    float centPi = fitFunc->GetParameter(1);
    float  sigPi = fitFunc->GetParameter(2);
    piFunc->SetParameters(ampPi,centPi,sigPi);
    piFunc->SetLineColor(kViolet+1);
    
    float  ampKP = fitFunc->GetParameter(3);
    float centKP = fitFunc->GetParameter(4);
    float  sigKP = fitFunc->GetParameter(5);
    kpFunc->SetParameters(ampKP,centKP,sigKP);
    kpFunc->SetLineColor(kGreen+2);
    
    float  ampE  = fitFunc->GetParameter(6);
    float centE  = fitFunc->GetParameter(7);
    float  sigE  = fitFunc->GetParameter(8);
    eFunc->SetParameters(ampE,centE,sigE);
    eFunc->SetLineColor(kBlue);
    
    piFunc->Draw("same");
    kpFunc->Draw("same");
    eFunc->Draw("same");
    
    TLine *l1 = new TLine( 2.5,0, 2.5,0.9*hNsigma->GetMaximum());
    TLine *l2 = new TLine(-1.0,0,-1.0,0.9*hNsigma->GetMaximum());
    l1->SetLineColor(kOrange+2);
    l2->SetLineColor(kOrange+2);
    l1->SetLineWidth(2);
    l2->SetLineWidth(2);
    
    l1->Draw();
    l2->Draw();
		
		//*
    float  totYield = fitFunc->Integral(-1.,2.5);
    float pionYield =  piFunc->Integral(-1.,2.5);
    float   kpYield =  kpFunc->Integral(-1.,2.5);
    float    eYield =   eFunc->Integral(-1.,2.5);
		//*/
		/*
     float  totYield = fitFunc->Integral(-2.,2.);
     float pionYield =  piFunc->Integral(-2.,2.);
     float   kpYield =  kpFunc->Integral(-2.,2.);
     float    eYield =   eFunc->Integral(-2.,2.);
     //*/
    
    
    
    totalYieldWhole += totYield;
    pionYieldWhole += pionYield;
    kpYieldWhole += kpYield;
    electYieldWhile += eYield;
    
    
    printf("\nTotal Yield\tPion Yield\tK/p Yield\tElect Yield\n");
    printf("%e\t%e\t%e\t%e\n",totYield,pionYield,kpYield,eYield);
    printf("Pion Fraction = %e\n",(pionYield/totYield));
    
		
		totalYeild[ihist] = totYield;

void genieta16calib()
{
  glob_t globbuf;

  int stat = glob (fileglob.c_str(), GLOB_MARK, NULL, &globbuf);
  if (stat) {
    switch (stat) {
    case GLOB_NOMATCH: cerr << "No file matching glob pattern "; break;
    case GLOB_NOSPACE: cerr << "glob ran out of memory "; break;
    case GLOB_ABORTED: cerr << "glob read error "; break;
    default: cerr << "unknown glob error stat=" << stat << " "; break;
    }
    cerr << fileglob << endl;
    exit(-1);
  }
  if (gl_verbose)
    cout<<globbuf.gl_pathc<<" files match the glob pattern"<<endl;
      
  for (size_t i=0; i<std::max((size_t)1,globbuf.gl_pathc); i++) {
    TVectorD vxday,vyfcamp,vzlumi;

    string path = globbuf.gl_pathv[i];

    // pick out day of year and mean(sum(fC))
    loadVectorsFromFile(path.c_str(),"%*lf %lf %lf %*lf %*s %*s %lf",vxday,vzlumi,vyfcamp);
    TGraph *pwg = new TGraph(vxday,vyfcamp);
    pwg->SetNameTitle(path.c_str(),path.c_str());

    //pwg->Print();

    v_graphs.push_back(pwg);

    optimizenorm(i,path);

    //pwg->Print();
  }

  // calc average and rms
  int ngraphs = (int)v_graphs.size();
  int npts    = v_graphs[0]->GetN();
  cout << ngraphs << " graphs, " << npts << " points per graph" << endl;
  TVectorD vx(npts),vyavg(npts),vyrms(npts), vxerr(npts);
  for (int i=0; i<npts; i++) {
    double x,y;
    v_graphs[0]->GetPoint(i,x,y);
    vx[i] = x;
    vxerr[i]=0;
    double yavg = y;
    for (int j=1; j<ngraphs; j++) {
      v_graphs[j]->GetPoint(i,x,y);
      if (x!=vx[i]) {
	cerr << x << " != " << vx[i] << endl;
	exit(-1);
      }
      yavg += y;
    }
    yavg /= (double)ngraphs;
    vyavg[i] = yavg;
    double var = 0.0;
    for (int j=0; j<ngraphs; j++) {
      v_graphs[j]->GetPoint(i,x,y);
      var += (y-yavg)*(y-yavg);
    }
    vyrms[i] = sqrt(var/(double)ngraphs);
  }

  TCanvas *c1 = new TCanvas("c1","c1",1500,600);
  v_graphs[0]->Draw("ALP");
  v_graphs[0]->GetHistogram()->GetYaxis()->SetRangeUser(0.0,2.0);
  v_graphs[0]->SetTitle("Average laser response from i#eta=-16 depth 3 over 2016, selected channels; Day # ; Arbitrary norm");
  for (int i=1; i<ngraphs; i++) {
    v_graphs[i]->Draw("LP same");
    v_graphs[i]->SetLineColor(15);
  }

  TGraphErrors *gravg = new TGraphErrors(vx,vyavg,vxerr,vyrms);

  gravg->Draw("LPE same");
  gravg->SetLineWidth(2);
  gPad->SetRightMargin(0.03);
  gPad->SetLeftMargin(0.08);
  gPad->Update();
  gPad->SaveAs("ieta-16calib.png");

  gravg->Print();
}

void MakePlots(TString filename, TString energy="8TeV", TString lumi=""){
  TString outDir=filename; outDir.ReplaceAll("fitres","img");
  outDir="tmp/";
  //std::map<TString, TH2F *> deltaNLL_map;
  
  /*------------------------------ Plotto */
  TCanvas *c = new TCanvas("c","c");
  
  TFile f_in(filename, "read");
  if(f_in.IsZombie()){
    std::cerr << "File opening error: " << filename << std::endl;
    return;
  }
  
  TList *KeyList = f_in.GetListOfKeys();
  std::cout << KeyList->GetEntries() << std::endl;
  for(int i =0; i <  KeyList->GetEntries(); i++){
    c->Clear();
    TKey *key = (TKey *)KeyList->At(i);
    if(TString(key->GetClassName())!="RooDataSet") continue;
    RooDataSet *dataset = (RooDataSet *) key->ReadObj();
    
    TString constTermName = dataset->GetName();
    TString alphaName=constTermName; alphaName.ReplaceAll("constTerm","alpha");

    TTree *tree = dataset2tree(dataset);
    TGraphErrors bestFit_ = bestFit(tree, alphaName, constTermName);
    TH2F *hist = prof2d(tree, alphaName, constTermName, "nll", "(12,-0.0005,0.0115,29,-0.0025,0.1425)",true);

//     //    deltaNLL_map.insert(std::pair <TString, TH2F *>(keyName,hist));
    hist->SaveAs(outDir+"/deltaNLL-"+constTermName+".root");
    hist->Draw("colz");
    bestFit_.Draw("P same");
    bestFit_.SetMarkerSize(2);

    Int_t iBinX, iBinY;
    Double_t x,y;
    hist->GetBinWithContent2(0,iBinX,iBinY);
    x= hist->GetXaxis()->GetBinCenter(iBinX);
    y= hist->GetYaxis()->GetBinCenter(iBinY);
    TGraph nllBestFit(1,&x,&y);

    nllBestFit.SetMarkerStyle(3);
    nllBestFit.SetMarkerColor(kRed);
    TList* contour68 = contourFromTH2(hist, 0.68);

    hist->Draw("colz");
    hist->GetZaxis()->SetRangeUser(0,50);
    bestFit_.Draw("P same");
    nllBestFit.Draw("P same");
    //contour68->Draw("same");
    c->SaveAs(outDir+"/deltaNLL-"+constTermName+".png");
    hist->SaveAs("tmp/hist-"+constTermName+".root");
    nllBestFit.SaveAs("tmp/nllBestFit.root");
    contour68->SaveAs("tmp/contour68.root");
    delete hist;
    hist = prof2d(tree, alphaName, constTermName, "nll", "(12,-0.0005,0.0115,29,-0.0025,0.1425)");
    RooHistPdf *histPdf = nllToL(hist);
    delete hist;
    RooDataSet *gen_dataset=histPdf->generate(*histPdf->getVariables(),1000000,kTRUE,kFALSE);
    TTree *genTree = dataset2tree(gen_dataset);
    genTree->SaveAs("tmp/genTree-"+constTermName+".root");
    delete gen_dataset;
    delete histPdf;
    
    TGraphErrors toyGraph = g(genTree, constTermName);
    TGraphErrors bestFitGraph = g(tree,alphaName, constTermName);
    TGraphErrors bestFitScanGraph = g(y, x);
    delete genTree;
    delete tree;
    toyGraph.SetFillColor(kGreen);
    toyGraph.SetLineColor(kBlue);
    toyGraph.SetLineStyle(2);
    bestFitGraph.SetLineColor(kBlack);
    bestFitScanGraph.SetLineColor(kRed);
    bestFitScanGraph.SetLineWidth(2);


    
    TMultiGraph g_multi("multigraph","");
    g_multi.Add(&toyGraph,"L3");
    g_multi.Add(&toyGraph,"L");
    g_multi.Add(&bestFitGraph, "L");
    g_multi.Add(&bestFitScanGraph, "L");
   
    g_multi.Draw("A");

    c->Clear();
    g_multi.Draw("A");
    c->SaveAs(outDir+"/smearing_vs_energy-"+constTermName+".png");

    //    TPaveText *pv = new TPaveText(0.7,0.7,1, 0.8);    
//     TLegend *legend = new TLegend(0.7,0.8,0.95,0.92);
//     legend->SetFillStyle(3001);
//     legend->SetFillColor(1);
//     legend->SetTextFont(22); // 132
//     legend->SetTextSize(0.04); // l'ho preso mettendo i punti con l'editor e poi ho ricavato il valore con il metodo GetTextSize()
//   //  legend->SetFillColor(0); // colore di riempimento bianco
//     legend->SetMargin(0.4);  // percentuale della larghezza del simbolo
    //    SetLegendStyle(legend);
//.........这里部分代码省略.........

void make1DLimit(TString combine_dir,TString type= "WH",bool blind=true){
  //TString combine_dir = "test_runSusyHgg/signalInj_sms_ChiWH_0_175/";
  //WH
  wh_limits.push_back(wh_125);
  wh_limits.push_back(wh_150);
  wh_limits.push_back(wh_175);
  wh_limits.push_back(wh_200);
  //HH
  hh_limits.push_back(hh_125);
  hh_limits.push_back(hh_150);
  hh_limits.push_back(hh_175);
  hh_limits.push_back(hh_200);

  TGraph obser( (200-125)/25 );
  TGraph graph( (200-125)/25 );
  TGraphAsymmErrors error( (200-125)/25 );
  TGraphAsymmErrors error2S( (200-125)/25 );

  TGraph obser_r( (200-125)/25 );
  TGraph graph_r( (200-125)/25 );
  TGraphAsymmErrors error_r( (200-125)/25 );
  TGraphAsymmErrors error_r2S( (200-125)/25 );
  
  TGraphErrors* theo = 0;
  if(type=="WH") theo = getTheoXSec("xsecs/CharginoNeutralino.txt");
  else theo = getTheoXSec("xsecs/Higgsino_ElectroHiggs.txt");
  //else theo = getTheoXSec("/home/amott/HggApp/SusyHgg/xsecs/Higgsino.txt");

  for(int m=125;m<=200;m+=25) {
    int i=(m-125)/25;
    TFile limit_file(Form("%s/higgsCombineChi%s_0_%d.Asymptotic.mH120.root",combine_dir.Data(),type.Data(),m) );
    TTree *limit_tree = (TTree*)limit_file.Get("limit");
    TTreeFormula limit_form("get_limit","limit",limit_tree);
    
    float down_2s = -1;
    float down = -1;
    float exp = -1;
    float up = -1;
    float up_2s = -1;
    float obs = -1;

    if( type == "WH" )
      {
	down_2s = wh_limits.at(i)[0];
	down = wh_limits.at(i)[1];
	exp = wh_limits.at(i)[2];
	up = wh_limits.at(i)[3];
	up_2s = wh_limits.at(i)[4];
	obs = wh_limits.at(i)[5];
      }
    else if ( type == "HH")
      {
	down_2s = hh_limits.at(i)[0];
	down = hh_limits.at(i)[1];
	exp = hh_limits.at(i)[2];
	up = hh_limits.at(i)[3];
	up_2s = hh_limits.at(i)[4];
	obs = hh_limits.at(i)[5];
      }
    else
      {
	std::cerr << "UNRECOGNIZED OPTION!!! QUITTING" << std::endl;
      }
    
    if(i==0) m+=5; //first point is actually at m=130
    graph.SetPoint(i,float(m), exp);
    error.SetPoint(i,float(m), exp);
    error2S.SetPoint(i, float(m), exp);
    error.SetPointError(i, 0, 0, exp-down, up-exp);
    error2S.SetPointError(i, 0 , 0 , exp-down_2s, up_2s-exp);
  
    graph_r.SetPoint(i,float(m),exp/theo->Eval(m));
    error_r.SetPoint(i,float(m),exp/theo->Eval(m));
    error_r2S.SetPoint(i,float(m),exp/theo->Eval(m));
    error_r.SetPointError(i,0,0,(exp-down)/theo->Eval(m),(up-exp)/theo->Eval(m));
    error_r2S.SetPointError(i, 0, 0, (exp-down_2s)/theo->Eval(m), (up_2s-exp)/theo->Eval(m) );
    
  
    obser.SetPoint(i,float(m),obs);
    obser_r.SetPoint(i,float(m),obs/theo->Eval(m));
    if(i==0) m-=5;
  }


    TCanvas cv;
    cv.SetLogy();
    cv.SetGrid(1,1);
    theo->SetMaximum(1e2);
    theo->SetMinimum(1e-2);
    theo->GetYaxis()->SetLabelSize(0.05);
    theo->GetYaxis()->SetTitleSize(0.06);
    theo->GetYaxis()->SetTitleOffset(0.8);
    theo->GetYaxis()->SetTitle("95% CL #sigma upper limit (pb)");
    theo->GetXaxis()->SetTitle("m_{chargino} (GeV)");
    if(type=="HH") theo->GetXaxis()->SetTitle("m_{neutralino} (GeV)");
    theo->SetFillColor(kBlue);
    theo->SetLineStyle(kDotted);
    theo->SetLineWidth(2);
    
    error.SetMaximum(1e2);
//.........这里部分代码省略.........