//--------------------------------------------------------------------------------------------------
void overlayFrame(TString text)
{
  // Overlay a linear frame from user coordinates (0 - 1, 0 - 1) and put the frame text

  // create new transparent pad for the text
  TPad *transPad = new TPad("transPad","Transparent Pad",0,0,1,1);
  transPad->SetFillStyle(4000);
  transPad->Draw();
  transPad->cd();

  // overlay the text in a well defined frame
  TText *plotText = new TText(0.01,0.01,text.Data());
  plotText->SetTextSize(0.04);
  plotText->SetTextColor(kBlue);
  plotText->Draw();

  return;
}

void first() {

   TCanvas *nut = new TCanvas("nut", "FirstSession",100,10,700,900);
   nut->Range(0,0,20,24);
   nut->SetFillColor(10);
   nut->SetBorderSize(2);

   TPaveLabel *pl = new TPaveLabel(3,22,17,23.7,
      "My first ROOT interactive session","br");
   pl->SetFillColor(18);
   pl->Draw();

   TText t(0,0,"a");
   t.SetTextFont(62);
   t.SetTextSize(0.025);
   t.SetTextAlign(12);
   t.DrawText(2,20.3,"ROOT is based on CINT, a powerful C/C++ interpreter.");
   t.DrawText(2,19.3,"Blocks of lines can be entered within {...}.");
   t.DrawText(2,18.3,"Previous typed lines can be recalled.");

   t.SetTextFont(72);
   t.SetTextSize(0.026);
   t.DrawText(3,17,"Root >  float x=5; float y=7;");
   t.DrawText(3,16,"Root >  x*sqrt(y)");
   t.DrawText(3,14,
      "Root >  for (int i=2;i<7;i++) printf(\"sqrt(%d) = %f\\n\",i,sqrt(i));");
   t.DrawText(3,10,"Root >  TF1 f1(\"f1\",\"sin(x)/x\",0,10)");
   t.DrawText(3, 9,"Root >  f1.Draw()");
   t.SetTextFont(81);
   t.SetTextSize(0.018);
   t.DrawText(4,15,"(float) 13.2288f");
   t.DrawText(4,13.3,"sqrt(2) = 1.414214");
   t.DrawText(4,12.7,"sqrt(3) = 1.732051");
   t.DrawText(4,12.1,"sqrt(4) = 2.000000");
   t.DrawText(4,11.5,"sqrt(5) = 2.236068");
   t.DrawText(4,10.9,"sqrt(6) = 2.449490");

   TPad *pad = new TPad("pad","pad",.2,.05,.8,.35);
   pad->Draw();
   pad->cd();
   pad->SetGrid();
   TF1 *f1 = new TF1("f1","sin(x)/x",0,10);
   f1->Draw();
}

TH1F* GetHisto(TFile* fin, string region, string process, string varname, float& norm, bool do_norm, float input_norm)
{
    string cname = CHANNEL_NAME+string("/")+region+"/"+varname;
    TCanvas* c = (TCanvas*) fin->Get(cname.c_str());
    string hname = "v:"+varname+"|p:"+process+"|r:"+region+string("|c:")+CHANNEL_NAME+string("|t:1DEntries");
    TH1F* h = 0;
    if(VERBOSE>0){
 	 cerr<<"cname :"<<cname<<endl;
   	 cerr<<"histo name: "<<hname<<endl;
   	 cerr<<"pointer: "<<c<<endl;
    } 
    TList* l = c->GetListOfPrimitives();
    TPad* pad = (TPad*) l->At(0);
    THStack* stack = (THStack*) pad->GetPrimitive("");
    h = (TH1F*) stack->GetHists()->FindObject(hname.c_str());
    if(do_norm) h->Scale(input_norm/h->Integral());
    norm = h->Integral();
    return (TH1F*) h->Clone();
}

// Render a TF2 looking like a rose.
// Author: Timur Pocheptsov
void glrose()
{
  //Define and set user's palette,
   //use polar system.
   const Int_t paletteSize = 10;
   Float_t rgb[paletteSize * 3] =
      {0.80f, 0.55f, 0.40f,
       0.85f, 0.60f, 0.45f,
       0.90f, 0.65f, 0.50f,
       0.95f, 0.70f, 0.55f,
       1.f,   0.75f, 0.60f,
       1.f,   0.80f, 0.65f,
       1.f,   0.85f, 0.70f,
       1.f,   0.90f, 0.75f,
       1.f,   0.95f, 0.80f,
       1.f,   1.f,   0.85f};

   Int_t palette[paletteSize] = {0};

   for (Int_t i = 0; i < paletteSize; ++i)
      palette[i] = TColor::GetColor(rgb[i * 3], rgb[i * 3 + 1], rgb[i * 3 + 2]);

   gStyle->SetPalette(paletteSize, palette);

   gStyle->SetCanvasPreferGL(1);
   TCanvas *cnv = new TCanvas("glc", "Surface sample", 200, 10, 600, 550);

   TPaveLabel *title = new TPaveLabel(0.04, 0.86, 0.96, 0.98,
                           "\"glsurf2pol\" option + user defined palette.");
   title->SetFillColor(32);
   title->Draw();

   TPad *rosePad  = new TPad("box", "box", 0.04, 0.04, 0.96, 0.8);
   rosePad->Draw();


   TF2 *fun = new TF2("a", "cos(y)*sin(x)+cos(x)*sin(y)", -6, 6, -6, 6);
   fun->SetContour(paletteSize);
   fun->SetNpx(30);
   fun->SetNpy(30);
   rosePad->cd();
   fun->Draw("glsurf2pol");
}

void draw_activation(TCanvas* c, Double_t cx, Double_t cy, 
                     Double_t radx, Double_t rady, Int_t whichActivation)
{
   TImage *activation = NULL;

   switch (whichActivation) {
   case 0:
      activation = TImage::Open("sigmoid-small.png");
      break;
   case 1:
      activation = TImage::Open("line-small.png");
      break;
   default:
      cout << "Activation index " << whichActivation << " is not known." << endl;
      cout << "You messed up or you need to modify network.C to introduce a new "
           << "activation function (and image) corresponding to this index" << endl;
   }

   if (activation == NULL) {
      cout << "Could not create an image... exit" << endl;
      return;
   }
  
   activation->SetConstRatio(kFALSE);

   radx *= 0.7;
   rady *= 0.7;
   TString name = Form("activation%f%f", cx, cy);
   TPad* p = new TPad(name+"", name+"", cx-radx, cy-rady, cx+radx, cy+rady);

   p->Draw();
   p->cd();

   activation->Draw();
   c->cd();
}

void root_util::set_box(Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax, 
                              Int_t lstyle, Double_t lwidth, Double_t lcolor,  Double_t fcolor, 
                              Int_t fptn, Char_t *loctype, Int_t gnum, Int_t xonly, Int_t yonly){

    //cout << pad_size <<endl;
    Int_t clen = clst.size();
    //cout << clen<<endl;
    Int_t cnum = gnum / pad_size;
    cout << "fcolor"<< fcolor<<endl;
    Int_t pad_num = gnum % pad_size + 1;
    //cout << pad_num<<endl;
    //cout << xmin <<endl;
    //cout << xmax <<endl;
    //cout << ymin <<endl;
    //cout << ymax <<endl;
    TPad *b;

    //b = new TEllipse(xmin, ymin, xmin, ymax);
    b = new TPad("box","box", xmin, ymin, xmax, ymax);
    b->SetFillColor(fcolor);
    clst[cnum]->cd(pad_num);
    b->Draw();
//    clst[cnum]->UseCurrentStyle();
}

TPad* plot_ratio(TCanvas *canv, bool up){
  canv->SetFillColor      (0);
  canv->SetBorderMode     (0);
  canv->SetBorderSize     (10);
  // Set margins to reasonable defaults
  canv->SetLeftMargin     (0.18);
  canv->SetLeftMargin     (0.17);
  canv->SetRightMargin    (0.05);
  canv->SetTopMargin      (0.12);
  canv->SetBottomMargin   (0.18);
  // Setup a frame which makes sense
  canv->SetFrameFillStyle (0);
  canv->SetFrameLineStyle (0);
  canv->SetFrameBorderMode(0);
  canv->SetFrameBorderSize(10);
  canv->SetFrameFillStyle (0);
  canv->SetFrameLineStyle (0);
  canv->SetFrameBorderMode(0);
  canv->SetFrameBorderSize(10);      

  canv->cd();
  TPad *pad = 0;
  if (up){
    pad = new TPad("upper","pad",0, 0.26 ,1 ,1);
    pad->SetBottomMargin(0.05);
    pad->SetTopMargin(0.09);
    pad->Draw();
    pad->cd();
    return pad;
  }
  else {
    pad = new TPad("lower","pad",0, 0,1 ,0.26);  
    pad->SetTopMargin(0.05);
    pad->SetBottomMargin(0.24);
    pad->Draw();
    pad->cd();
    return pad;
  }

};

void DrawALICELogo(Bool_t prel, Float_t x1, Float_t y1, Float_t x2, Float_t y2)
{
    // correct for aspect ratio of figure plus aspect ratio of pad (coordinates are NDC!)
    x2 = x1 + (y2 - y1) * (466. / 523) * gPad->GetWh() * gPad->GetHNDC() / (gPad->GetWNDC() * gPad->GetWw());
    // Printf("%f %f %f %f", x1, x2, y1, y2);

    TPad *myPadLogo = new TPad("myPadLogo", "Pad for ALICE Logo", x1, y1, x2, y2);
    myPadLogo->SetLeftMargin(0);
    myPadLogo->SetTopMargin(0);
    myPadLogo->SetRightMargin(0);
    myPadLogo->SetBottomMargin(0);
    myPadLogo->Draw();
    myPadLogo->cd();
    TASImage *myAliceLogo = new TASImage((prel) ? "~/alice_logo_preliminary.eps" : "~/alice_logo_performance.eps");
    myAliceLogo->Draw();
}

void DrawALICELogo(Float_t x1, Float_t y1, Float_t x2, Float_t y2)
{
// Correct for aspect ratio of figure plus aspect ratio of pad.
// Coordinates are NDC!

  x2 = x1 + (y2 - y1)*0.891*gPad->GetCanvas()->GetWindowHeight()*gPad->GetHNDC() / (gPad->GetWNDC() * gPad->GetCanvas()->GetWindowWidth());
  
  TPad *myPadLogo = new TPad("myPadLogo","Pad for ALICE Logo", x1, y1, x2, y2);
  myPadLogo->SetLeftMargin(0);
  myPadLogo->SetTopMargin(0);
  myPadLogo->SetRightMargin(0);
  myPadLogo->SetBottomMargin(0);
  myPadLogo->Draw();
  myPadLogo->cd();
  TASImage *myAliceLogo = 
  new TASImage("alice_logo_preliminary.eps");
  myAliceLogo->Draw("same");
}

void DrawGraph(string name, TGraphErrors* graph)
{
    TCanvas* canvas = new TCanvas(name.c_str(),"",1200,900);
    canvas->SetFillColor(kWhite);
    canvas->cd();
    TPad* pad = new TPad();
    // left right bottom top
    pad->SetMargin(0.15,0.05,0.15,0.05);
    pad->SetTicks(1,1);
    pad->SetFillColor(kWhite);
    pad->Draw();
    pad->cd();
    TGraphErrors* oldgraph = graph;
    graph = new TGraphErrors(*graph);
    delete oldgraph;
    float textsize = 0.055;
    graph->SetMinimum(0);
    graph->SetMaximum(100);
    graph->SetTitle("");
    graph->GetXaxis()->SetTitle("#it{m}(#it{K^{#plus}K^{#minus}}) [MeV/#it{c}^{2}]");
    graph->GetYaxis()->SetTitle("Efficiency [%]");
    graph->SetLineColor(kBlack);
    graph->SetLineWidth(1);
    graph->SetMarkerSize(1);
    graph->SetMarkerStyle(8);
    graph->SetMarkerColor(kBlack);
    graph->Draw("AP");
    graph->GetXaxis()->SetTitleSize(textsize*1.2);
    graph->GetYaxis()->SetTitleSize(textsize*1.2);
    graph->GetXaxis()->SetLabelSize(textsize);
    graph->GetYaxis()->SetLabelSize(textsize);
    graph->GetXaxis()->SetTitleFont(132);
    graph->GetYaxis()->SetTitleFont(132);
    graph->GetXaxis()->SetLabelFont(132);
    graph->GetYaxis()->SetLabelFont(132);
    graph->GetYaxis()->CenterTitle();
    double _blurbx        = 0.75;
    double _blurby        = 0.80;
    string _blurbtext     = "#splitline{LHCb}{#scale[0.75]{Preliminary}}";
    TLatex* _blurb = new TLatex(_blurbx,_blurby,_blurbtext.c_str());
    _blurb->SetTextFont(132);
    _blurb->SetTextColor(1);
    _blurb->SetNDC(kTRUE);
    _blurb->SetTextAlign(11);
    _blurb->SetTextSize(0.07);
    _blurb->Draw();
    canvas->SaveAs((name+".pdf").c_str());
    canvas->Write();
}

//___________________________ 
void RangeClicked() {

	TPad *pad = (TPad *) gCanvas->GetSelectedPad();
	if (!pad) return;

	int eventtype = pad->GetEvent(); 
	
	if ( pad->GetEvent() == 12 && !gMinimumIsSet && !gMaximumIsSet) {   //   kButton1Down   //kButton1Double=61    //middle button 1 click = 12

		//Get the abscissa  
		gMin = pad->GetEventX();
		gMin = pad->AbsPixeltoX(gMin);
		printf(" [ %.3f , \n",gMin);
		gMinimumIsSet = true ;
		gMaximumIsSet = false ; 
		return ; // after this return the condition of the block below will be satisfied  
		} 
		
	if ( pad->GetEvent() == 12 && gMinimumIsSet && !gMaximumIsSet) {

		//Get the abscissa  
		gMax = pad->GetEventX();
		gMax = pad->AbsPixeltoX(gMax);
		
		//print the values
		if(gMin>=gMax) {	cout << " WARNING : Min >= Max , "<< " choose a value > "<< gMin << " \n" << endl ; return ; }
		else printf(" [ %.3f , %.3f ]\n",gMin,gMax);
		
		//reset 
		gMinimumIsSet=false ;
		gMaximumIsSet=false ;
		return ; 
		}
	
	return ;

}

void na49view()
{
  //
  // This macro generates
  // a begin_html <a href="gif/na49canvas.gif">Canvas</a> end_html
  // with 2 views of the NA49 detector using the old obsolete geometry package.
  // Author: Rene Brun

  TCanvas *c1 = new TCanvas("c1", "The NA49 canvas", 200, 10, 700, 780);

  gBenchmark->Start("na49view");

  TPad *all = new TPad("all", "A Global view of NA49", 0.02, 0.02, 0.48, 0.82, 28);
  TPad *tof = new TPad("tof", "One Time Of Flight element", 0.52, 0.02, 0.98, 0.82, 28);
  all->Draw();
  tof->Draw();
  TPaveLabel *na49title = new TPaveLabel(0.04, 0.86, 0.96, 0.98, "Two views of the NA49 detector");
  na49title->SetFillColor(32);
  na49title->Draw();
  //
  TFile *nageom = new TFile("na49.root");
  if (!nageom || nageom->IsZombie()) return;
  TGeometry *n49 = (TGeometry *)gROOT->FindObject("na49");
  n49->SetBomb(1.2);
  n49->cd(); // Set current geometry
  all->cd(); // Set current pad
  n49->Draw();
  c1->Update();
  tof->cd();
  TNode *TOFR1 = n49->GetNode("TOFR1");
  TOFR1->Draw();
  c1->Update();

  gBenchmark->Show("na49view");

  // To have a better and dynamic view of any of these pads,
  // you can click with the middle button of your mouse to select it.
  // Then select "View with x3d" in the VIEW menu of the Canvas.
  // Once in x3d, you are in wireframe mode by default.
  // You can switch to:
  //   - Hidden Line mode by typing E
  //   - Solid mode by typing R
  //   - Wireframe mode by typing W
  //   - Stereo mode by clicking S (and you need special glasses)
  //   - To leave x3d type Q
}

void DrawGrid()
{

Int_t ncol = 100;
  Int_t nrow = 200;
  
  
   TPad *grid = new TPad("grid","",0,0,1,1);
   grid->Draw();
   grid->cd();
   grid->SetGrid();
   grid->SetFillStyle(4000);
   grid->SetFrameFillStyle(0);


   TH2 *hgrid = new TH2C("hgrid","", ncol+1, -5, ncol+5, nrow, -5, nrow-1+5);
   hgrid->Draw();
   hgrid->GetXaxis()->SetNdivisions(6,100);
   hgrid->GetYaxis()->SetNdivisions(200);
   hgrid->GetYaxis()->SetLabelOffset(999.);
   hgrid->GetXaxis()->SetLabelOffset(999.);
}

TCanvas* example_plot( int iPeriod, int iPos )
{ 
  //  if( iPos==0 ) relPosX = 0.12;

  int W = 800;
  int H = 600;

  // 
  // Simple example of macro: plot with CMS name and lumi text
  //  (this script does not pretend to work in all configurations)
  // iPeriod = 1*(0/1 7 TeV) + 2*(0/1 8 TeV)  + 4*(0/1 13 TeV) 
  // For instance: 
  //               iPeriod = 3 means: 7 TeV + 8 TeV
  //               iPeriod = 7 means: 7 TeV + 8 TeV + 13 TeV 
  // Initiated by: Gautier Hamel de Monchenault (Saclay)
  // Updated by:   Dinko Ferencek (Rutgers)
  //
  int H_ref = 600; 
  int W_ref = 800; 

  // references for T, B, L, R
  float T = 0.08*H_ref;
  float B = 0.12*H_ref; 
  float L = 0.12*W_ref;
  float R = 0.04*W_ref;

  TString canvName = "FigExample_";
  canvName += W;
  canvName += "-";
  canvName += H;
  canvName += "_";  
  canvName += iPeriod;
  if( writeExtraText ) canvName += "-prelim";
  if( iPos%10==0 ) canvName += "-out";
  else if( iPos%10==1 ) canvName += "-left";
  else if( iPos%10==2 )  canvName += "-center";
  else if( iPos%10==3 )  canvName += "-right";

  TCanvas* canv = new TCanvas(canvName,canvName,50,50,W,H);
  canv->SetFillColor(0);
  canv->SetBorderMode(0);
  canv->SetFrameFillStyle(0);
  canv->SetFrameBorderMode(0);
  canv->SetLeftMargin( L/W );
  canv->SetRightMargin( R/W );
  canv->SetTopMargin( T/H );
  canv->SetBottomMargin( B/H );
  canv->SetTickx(0);
  canv->SetTicky(0);

  TH1* h = new TH1F("h","h",40,70,110);
  h->GetXaxis()->SetNdivisions(6,5,0);
  h->GetXaxis()->SetTitle("m_{e^{+}e^{-}} (GeV)");  
  h->GetYaxis()->SetNdivisions(6,5,0);
  h->GetYaxis()->SetTitleOffset(1);
  h->GetYaxis()->SetTitle("Events / 0.5 GeV");  

  h->SetMaximum( 260 );
  if( iPos==1 ) h->SetMaximum( 300 );
  h->Draw();

  int histLineColor = kOrange+7;
  int histFillColor = kOrange-2;
  float markerSize  = 1.0;

  {
    TLatex latex;
				
    int n_ = 2;

    float x1_l = 0.92;
    float y1_l = 0.60;

    float dx_l = 0.30;
    float dy_l = 0.18;
    float x0_l = x1_l-dx_l;
    float y0_l = y1_l-dy_l;

    TPad* legend = new TPad("legend_0","legend_0",x0_l,y0_l,x1_l, y1_l );
    //    legend->SetFillColor( kGray );
    legend->Draw();
    legend->cd();
		
    float ar_l = dy_l/dx_l;
		
    float x_l[1];
    float ex_l[1];
    float y_l[1];
    float ey_l[1];
		
    //    float gap_ = 0.09/ar_l;
    float gap_ = 1./(n_+1);
		
    float bwx_ = 0.12;
    float bwy_ = gap_/1.5;
		
    x_l[0] = 1.2*bwx_;
    //    y_l[0] = 1-(1-0.10)/ar_l;
    y_l[0] = 1-gap_;
    ex_l[0] = 0;
//.........这里部分代码省略.........

TCanvas *drawRatioPlot(TH1D *prediction, TH1D *sr, TH1D *data, TString xTitle, TString filename, double ecaloCut){

  gStyle -> SetPadLeftMargin(0.20);

  data->SetMarkerStyle(20);
  data->SetMarkerColor(kGreen);
  data->SetLineColor(kGreen);


  
  TCanvas *c = new TCanvas("c"+filename,"c",0,0,500,500);

  float y = 0.3;

  TPad *pad1     = new TPad("pad1", "Control Plots 1", 0.01, y, 0.99, 0.99);
  TPad *padRatio = new TPad("rp1", "Ratio1", 0.01, 0.01, 0.99, y-0.01);

  pad1->SetNumber(100);
  pad1->SetTicks(0,1);
  cout<<"number pad1     = "<<pad1->GetNumber()<<endl;
  cout<<"number padRatio = "<<padRatio->GetNumber()<<endl;

  TH1D *ratio = 0;

  //ratio = (TH1D*) sr->Clone();
  //ratio->Divide(prediction);

  ratio = (TH1D*) prediction->Clone();
  ratio->Divide(data);



  for(int i=1; i<=ratio->GetNbinsX();i++){

    if(ratio->GetBinContent(i) != 0){
      cout<<"N in CR in "<<i<<". bin ="<<prediction->GetBinContent(i)<<endl;
      cout<<"N in SR in "<<i<<". bin ="<<sr->GetBinContent(i)<<endl;
      cout<<"Rel. difference in "<<i<<". bin ="<<(1./ratio->GetBinContent(i)-1.)*100<<"%"<<endl; 
    }
    else if(sr->GetBinContent(i) == 0 && prediction->GetBinContent(i) !=0)    cout<<"Scaling Factor in "<<i<<". bin <"<<prediction->GetBinContent(i)/1.15<<" +/- "<<ratio->GetBinError(i)<<endl;
    else if(sr->GetBinContent(i) != 0 && prediction->GetBinContent(i) ==0)    cout<<"Scaling Factor in "<<i<<". bin <"<<(sr->GetEntries()/prediction->GetEntries())/sr->GetBinContent(i)<<" +/- "<<ratio->GetBinError(i)<<endl;

  }

  ratio->GetYaxis()->SetTitle("#frac{CR (MC)}{CR (data)}");
  ratio->SetTitle("");
  ratio->SetLabelSize(0.1,"X");
  ratio->SetLabelSize(0.1,"Y");
  ratio->SetTitleOffset(0.5,"Y");
  ratio->SetTitleSize(0.15,"Y");


  padRatio->cd();
  //ratio->GetYaxis()->SetRangeUser(0,2);
  ratio->Draw("e");

  // Draw line at one!
  float xmin = ratio->GetXaxis()->GetXmin();
  float xmax = ratio->GetXaxis()->GetXmax();
  TLine *line = new TLine(xmin,1,xmax,1);
  line->SetLineWidth(2);
  line->Draw("same");
  padRatio->Modified();

  TLegend *leg = new TLegend(0.5,0.7,0.9,0.9);
  leg->AddEntry(sr,"SR (MC)","l"); 
  leg->AddEntry(prediction,"lepton CR (MC)","pel"); 
 
  pad1->cd();
  pad1->SetLogy();
  //  pad1->SetLogx();

  sr->SetLineColor(kRed);
  sr->SetMarkerColor(kRed);
  sr->SetMarkerStyle(20);
  sr->SetTitle("");
  prediction->SetMarkerStyle(20);
  prediction->SetTitle("");
  prediction->GetXaxis()->SetTitle(xTitle);
  sr->GetXaxis()->SetTitle(xTitle);

  
  

  prediction->SetTitleSize(0.07,"X");
  prediction->GetXaxis()->SetTitleOffset(0.7);
  sr->SetTitleSize(0.07,"X");
  sr->GetXaxis()->SetTitleOffset(0.7);

  double maximum = 0;
  double minimum = 1000000;
  if(sr->GetMinimum()!=0 && sr->GetMinimum()<minimum){
    minimum=sr->GetMinimum()*0.5;
  }
  if(prediction->GetMinimum()!=0 && prediction->GetMinimum()<minimum){
    minimum=prediction->GetMinimum()*0.5;
  }
  if(data->GetMinimum()!=0 && data->GetMinimum()<minimum){
    minimum=data->GetMinimum()*0.5;
  }
//.........这里部分代码省略.........

// draws the decay vertices of the Lambda_s
void genDisplay(std::string fullPath, int iGen, int iReco = -1, bool savePdf = false)
{
    const int fVerbose(0);
    setTDRStyle();
    //gStyle->SetOptStat(112211);
    gStyle->SetOptStat(0);
    gStyle->SetPalette(1);
    // Canvases
    crz = new TCanvas("crz","crz",1000,600);
    TPad* padrz = (TPad*)crz->cd(1);
    setPadMargins(padrz, 0.02, 0.05, 0.02, 0.05);
    crphi = new TCanvas("crphi","crphi",600,600);
    TPad* padrphi = (TPad*)crphi->cd(1);
    setPadMargins(padrphi, 0.02, 0.05, 0.02, 0.05);
    // Open file
    TFile *f = TFile::Open(fullPath.c_str());
    if (f==0)
    {
	cout << "File " << fullPath << " not found -- exiting" << endl;
	return;
    }
    if(fVerbose>0)
	cout << "Succesfully opened file " << fullPath << endl;
    // Get TTree with GenEvents
    TTree* tgen = (TTree*) f->Get("genevents");
    if(fVerbose>0) cout << "Got TTree with " << tgen->GetEntries() << " entries" << endl;
    // Do a cut, if needed
    tgen->Draw(">>lst","ptmu1>3&&ptmu2>3&&TMath::Abs(etamu1)<2.5&&TMath::Abs(etamu2)<2.5");
    TEventList *lst;
    lst = (TEventList*)gDirectory->Get("lst");
    tgen->SetEventList(lst);

    // Get TTree with iRecoEvents
    TTree* treco = (TTree*) f->Get("events");
    if(fVerbose>0) cout << "Got TTree with " << treco->GetEntries() << " entries" << endl;

    // Do plots
    const int nbins(30);
    const double rangeR(30), rangeZ(100); // standard
    //const double rangeR(60), rangeZ(200); // zoomed out
    //const double rangeR(100), rangeZ(200); // zoomed in
    //const double rangeR(12), rangeZ(30); // zoomed further in
    padrz->cd();
    TH2F *hrz = new TH2F("hdisp","", nbins,-rangeZ,rangeZ,nbins,-rangeR,rangeR);
    hrz->Draw();
    padrphi->cd();
    TH2F *hrphi = new TH2F("hdisp","", nbins,-rangeR,rangeR,nbins,-rangeR,rangeR);
    hrphi->Draw();

    // add tracker
    padrz->Modified();
    padrz->Update();
    drawTrackerRZ(padrz);
    padrphi->Modified();
    padrphi->Update();
    drawTrackerRPhi(padrphi);

    // set branch addresses for gentree
    double vrl0,vxl0,vyl0,vzl0,vrlb,vxlb,vylb,vzlb;
    double pmu1,phimu1,etamu1,pmu2,phimu2,etamu2;
    double ppr,ppi,phipr,phipi,etapr,etapi;
    tgen->SetBranchAddress("vrl0",&vrl0);
    tgen->SetBranchAddress("vxl0",&vxl0);
    tgen->SetBranchAddress("vyl0",&vyl0);
    tgen->SetBranchAddress("vzl0",&vzl0);
    tgen->SetBranchAddress("vrlb",&vrlb);
    tgen->SetBranchAddress("vxlb",&vxlb);
    tgen->SetBranchAddress("vylb",&vylb);
    tgen->SetBranchAddress("vzlb",&vzlb);
    tgen->SetBranchAddress("pmu1",&pmu1);
    tgen->SetBranchAddress("etamu1",&etamu1);
    tgen->SetBranchAddress("phimu1",&phimu1);
    tgen->SetBranchAddress("pmu2",&pmu2);
    tgen->SetBranchAddress("etamu2",&etamu2);
    tgen->SetBranchAddress("phimu2",&phimu2);
    tgen->SetBranchAddress("ppr",&ppr);
    tgen->SetBranchAddress("etapr",&etapr);
    tgen->SetBranchAddress("phipr",&phipr);
    tgen->SetBranchAddress("ppi",&ppi);
    tgen->SetBranchAddress("etapi",&etapi);
    tgen->SetBranchAddress("phipi",&phipi);
    int genrun, genls, genevt;
    tgen->SetBranchAddress("run",&genrun);
    tgen->SetBranchAddress("LS",&genls);
    tgen->SetBranchAddress("event",&genevt);

    // set branch addresses for recotree
    double rvrl0,rvxl0,rvyl0,rvzl0,rvrlb,rvxlb,rvylb,rvzlb;
    double rpmu1,retamu1,rphimu1,rpmu2,retamu2,rphimu2;
    double rppr,rppi,retapr,retapi,rphipr,rphipi;
    treco->SetBranchAddress("vrl0",&rvrl0);
    treco->SetBranchAddress("vxl0",&rvxl0);
    treco->SetBranchAddress("vyl0",&rvyl0);
    treco->SetBranchAddress("vzl0",&rvzl0);
    treco->SetBranchAddress("vrlb",&rvrlb);
    treco->SetBranchAddress("vxlb",&rvxlb);
    treco->SetBranchAddress("vylb",&rvylb);
    treco->SetBranchAddress("vzlb",&rvzlb);
    treco->SetBranchAddress("rpt1m",&rpmu1);
//.........这里部分代码省略.........

void
QA_Draw_Calorimeter_Sum_TrackProjEP(
    const char * qa_file_name_new =
        "/phenix/u/jinhuang/links/ePHENIX_work/sPHENIX_work/production_analysis_updates/spacal1d/fieldmap/G4Hits_sPHENIX_pi-_eta0.30_32GeV-0000.root_qa.root",
    const char * qa_file_name_ref =
        "/phenix/u/jinhuang/links/ePHENIX_work/sPHENIX_work/production_analysis_updates/spacal1d/fieldmap/G4Hits_sPHENIX_pi+_eta0.30_32GeV-0000.root_qa.root")
//QA_Draw_Calorimeter_Sum_TrackProjEP(const char * qa_file_name_new =
//    "data/G4sPHENIXCells_100e24GeV.root_qa.root",
//    const char * qa_file_name_ref =
//        "data/G4Hits_sPHENIX_e-_eta0_24GeV-0000.root_qa.root")
{

  SetOKStyle();
  gStyle->SetOptStat(0);
  gStyle->SetOptFit(1111);
  TVirtualFitter::SetDefaultFitter("Minuit2");

  TFile * qa_file_new = new TFile(qa_file_name_new);
  assert(qa_file_new->IsOpen());

  TFile * qa_file_ref = NULL;
  if (qa_file_name_ref)
    {
      qa_file_ref = new TFile(qa_file_name_ref);
      assert(qa_file_ref->IsOpen());
    }

  // obtain normalization
  double Nevent_new = 1;
  double Nevent_ref = 1;

  // obtain normalization
  double Ntrack_new = 0;
  double Ntrack_ref = 0;

  if (qa_file_new)
    {
      TH1D * h_norm = (TH1D *) qa_file_new->GetObjectChecked(
          TString("h_QAG4Sim_CalorimeterSum_Normalization"), "TH1D");
      assert(h_norm);

      Nevent_new = h_norm->GetBinContent(h_norm->GetXaxis()->FindBin("Event"));
      Ntrack_new = h_norm->GetBinContent(h_norm->GetXaxis()->FindBin("Track"));
    }
  if (qa_file_ref)
    {
      TH1D * h_norm = (TH1D *) qa_file_ref->GetObjectChecked(
          TString("h_QAG4Sim_CalorimeterSum_Normalization"), "TH1D");
      assert(h_norm);

      Nevent_ref = h_norm->GetBinContent(h_norm->GetXaxis()->FindBin("Event"));
      Ntrack_ref = h_norm->GetBinContent(h_norm->GetXaxis()->FindBin("Track"));
    }

  TCanvas *c1 = new TCanvas("QA_Draw_Calorimeter_Sum_TrackProjEP",
      "QA_Draw_Calorimeter_Sum_TrackProjEP", 1800, 600);
  c1->Divide(3, 1);
  int idx = 1;
  TPad * p;

  p = (TPad *) c1->cd(idx++);
  c1->Update();
  p->SetLogy();

  if (Ntrack_new>0)
    {
      TH1F * h_new = (TH1F *) qa_file_new->GetObjectChecked(
          "h_QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP", "TH1F");
      assert(h_new);

      h_new->Sumw2();
      h_new->Scale(1. / Ntrack_new);

      TH1F * h_ref = NULL;
      if (qa_file_ref)
        {
          TH1F * h_ref = (TH1F *) qa_file_ref->GetObjectChecked(
              "h_QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP", "TH1F");
          assert(h_ref);

          h_ref->Scale(1. / Ntrack_ref);
        }

      h_new->GetYaxis()->SetTitleOffset(1.5);
      h_new->GetYaxis()->SetTitle("Count / track / bin");
      //      h_new->GetXaxis()->SetRangeUser(-0, .1);

      DrawReference(h_new, h_ref);
    }


  p = (TPad *) c1->cd(idx++);
  c1->Update();
  p->SetLogy();

  if (Ntrack_new>0)
    {
      TH1F * h_new = (TH1F *) qa_file_new->GetObjectChecked(
          "h_QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP", "TH1F");
      assert(h_new);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    TFile *output = new TFile("20170721_KKbarAdded2_fixedwidth42_recon_pf100_scaled_error05.root", "RECREATE");

    TH1D *kstar0mass = new TH1D("kstar0mass", Form("Fit value of M*_{0} vs. p_{T} for %s", particles[PARTICLE_NUM].c_str()), NUM_PT_BINS, 0.0, 4.0);
    TH1D *kstar0width = new TH1D("kstar0width", Form("#Gamma_{tot}(M=M*_{0}) vs p_{T} for %s", particles[PARTICLE_NUM].c_str()), NUM_PT_BINS, 0.0, 4.0);
    TH1D *kstar0collWidth = new TH1D("kstar0collWidth", Form("Fit value of #Gamma_{coll} component vs. p_{T} for %s", particles[PARTICLE_NUM].c_str()), NUM_PT_BINS,0.0, 4.0);
    TH1D *kstar0decWidth = new TH1D("kstar0decWidth", Form("#Gamma_{dec}(M=M*_{0}) component vs. p_{T} for %s;p_{T} (GeV/c);Width (GeV/c^2)", particles[PARTICLE_NUM].c_str()), NUM_PT_BINS,0.0, 4.0);
   
    kstar0mass->GetXaxis()->SetTitle("p_{T} (GeV/c)");
    kstar0mass->GetYaxis()->SetTitle("Mass (GeV/c^{2})");
    kstar0width->GetXaxis()->SetTitle("p_{T} (GeV/c)");
    kstar0width->GetYaxis()->SetTitle("Width (GeV/c^2)");
    kstar0collWidth->GetXaxis()->SetTitle("p_{T} (GeV/c)");
    kstar0collWidth->GetYaxis()->SetTitle("Width (GeV/c^2)");

    kstar0mass->SetStats(kFALSE);
    kstar0width->SetStats(kFALSE);
    kstar0collWidth->SetStats(kFALSE);
    kstar0decWidth->SetStats(kFALSE);

    TF1 *massline = new TF1("massline", "[0]", 0.0, 4.0);
    massline->SetParameter(0, 0.892);
    massline->SetLineColor(2);
    massline->SetLineStyle(7);

    TF1 *widthline = new TF1("widthline", "[0]", 0.0, 4.0);
    widthline->SetParameter(0, 0.042);

    double mass = 0.0, width = 0.0, collWidth = 0.0, massBG=0.0;
    double massError = 0.0, widthError= 0.0, collWidthError = 0.0, massBGError=0.0;

    TCanvas *canvas[9];
    TCanvas *diffCanvas[9];
    TPaveStats *st;
    TPad *pad;

    //ofstream integrals;
    //integrals.open("kstarbar_integrals.txt");

    for(int nfile = 0; nfile < NUM_PT_BINS; nfile++){
        double meanPT = (double)(nfile*2+1)/10.0;
        string filename = folder+files[nfile];
        string ptLower = filename.substr(filename.find("[")+1, 3);
        string ptHigher = filename.substr(filename.find(",")+1, 3);   
        TH1D* histos[8];
        TH1D* newHistos[8];
        TH1D* diffHistos[8];
        TH1D* bg[8];
        for(int i=0; i<8; i++){
            if(nfile<5){
                histos[i] = new TH1D(Form("ptbin0%dparticle%d",nfile*2+1, i), Form("Invariant Mass for (%s), %s < p_{T} < %s",particles[i].c_str(), ptLower.c_str(), ptHigher.c_str()), NUM_MASS_BINS, MASS_LOW, MASS_HIGH);
            newHistos[i] = new TH1D(Form("newptbin0%dparticle%d",nfile*2+1, i), Form("Invariant Mass for (%s), %s < p_{T} < %s",particles[i].c_str(), ptLower.c_str(), ptHigher.c_str()), 250, MASS_LOW, MASS_HIGH);

            }else{
                histos[i] = new TH1D(Form("ptbin%dparticle%d",nfile*2+1, i), Form("Invariant Mass for (%s), %s < p_{T} < %s",particles[i].c_str(), ptLower.c_str(), ptHigher.c_str()), NUM_MASS_BINS, MASS_LOW, MASS_HIGH);
                newHistos[i] = new TH1D(Form("newptbin%dparticle%d",nfile*2+1, i), Form("Invariant Mass for (%s), %s < p_{T} < %s",particles[i].c_str(), ptLower.c_str(), ptHigher.c_str()), 250, MASS_LOW, MASS_HIGH);

            }
            histos[i]->GetXaxis()->SetTitle("Invariant Mass (GeV/c^{2})");
            histos[i]->GetYaxis()->SetTitle("Counts");
        }

        ifstream input;
        input.open(filename.c_str());
        string line = "";
        if(input.good()){
            getline(input, line);

//______________________________________________________________________________
void KVCanvas::HandleInput(EEventType event, Int_t px, Int_t py)
{
   // Handle Input Events.
   //
   //  Handle input events, like button up/down in current canvas.

   if (fFreezed) return;

   TPad*    pad;
   TPad*    prevSelPad = (TPad*) fSelectedPad;
   TObject* prevSelObj = fSelected;

   fPadSave = (TPad*)gPad;
   cd();        // make sure this canvas is the current canvas

   fEvent  = event;
   fEventX = px;
   fEventY = py;

   Int_t sign = 0;
   Bool_t sendOrder = true;

   if (fHasDisabledClasses && fSelected) {
      if (fDisabledClasses.Contains(fSelected->ClassName())) sendOrder = false;
   }

   if (fHasDisabledObject && fSelected) {
      if (fDisabledObjects.Contains(fSelected)) sendOrder = false;
   }

   switch (event) {

      case kMouseMotion:
         // highlight object tracked over
         pad = Pick(px, py, prevSelObj);
         if (!pad) return;

         EnterLeave(prevSelPad, prevSelObj);

         gPad = pad;   // don't use cd() we will use the current
         // canvas via the GetCanvas member and not via
         // gPad->GetCanvas

         if (sendOrder) fSelected->ExecuteEvent(event, px, py);

         RunAutoExec();

         if (fAgeOfEmpire && (fSelected->InheritsFrom("TH2"))) {
            TH2* TheHisto = (TH2*) FindHisto();//fSelected;

            Double_t size = 0.4 - 0.35 * fVenerMode;

            Int_t dX = 0;
            Int_t dY = 0;

            Double_t ppx = AbsPixeltoX(px);
            Double_t ppy = AbsPixeltoY(py);

            TAxis* ax = TheHisto->GetXaxis();
            Int_t X0 = ax->GetFirst();
            Int_t X1 = ax->GetLast();
            Int_t NbinsX = ax->GetNbins();
            px = ax->FindBin(ppx);

            Double_t ddX   = (X1 + X0) * 0.5 - px;
            Double_t distX = TMath::Abs(ddX) / (X1 - X0);
            if (distX >= 0.5) return;

            TAxis* ay = TheHisto->GetYaxis();
            Int_t Y0 = ay->GetFirst();
            Int_t Y1 = ay->GetLast();
            Int_t NbinsY = ay->GetNbins();
            py = ay->FindBin(ppy);

            Double_t ddY   = (Y1 + Y0) * 0.5 - py;
            Double_t distY = TMath::Abs(ddY) / (Y1 - Y0);
            if (distY >= 0.5) return;

            if ((distX <= size) && (distY <= size)) return;
            dX = TMath::Nint(ddX * (0.05 + 0.05 * fVenerMode));
            dY = TMath::Nint(ddY * (0.05 + 0.05 * fVenerMode));

            if (TMath::Abs(dX) < 1) dX = TMath::Sign(1., ddX);
            if (TMath::Abs(dY) < 1) dY = TMath::Sign(1., ddY);


            Bool_t up = false;
            if ((X0 - dX > 0) && (X1 - dX < NbinsX)) {
               ax->SetRange(X0 - dX, X1 - dX);
               up = true;
            }
            if ((Y0 - dY > 0) && (Y1 - dY < NbinsY)) {
               ay->SetRange(Y0 - dY, Y1 - dY);
               up = true;
            }
            if (up) {
               Modified();
               Update();
            }
//.........这里部分代码省略.........