C++ TAxis::GetBinCenter方法代码示例

    static bool CheckAxisCompatibility(const TAxis& a, const TAxis& b)
    {
        static const double delta_max = 0.001;
        static const double rel_delta_max = 0.01;

        if(a.GetNbins() != b.GetNbins()) return false;
        for(Int_t n = 1; n <= a.GetNbins(); ++n) {
            const double c_a = a.GetBinCenter(n);
            const double c_b = b.GetBinCenter(n);
            const double delta = c_a - c_b;
            if( ( c_a && std::abs(delta/c_a) > rel_delta_max ) || ( !c_a && std::abs(delta) > delta_max ) )
                return false;
        }
        return true;
    }

//////////////////////////////////////////////////
// addXAxis - add X axis information
unsigned int addXAxis(SEXP data, SEXP dataNames, unsigned int j, TH1* hist)
{
  int n = hist->GetNbinsX();
  TAxis* axis = hist->GetXaxis();
  
  // Determine breaks--
  // Add to list
  SEXP breaks = addNumericVector(data, dataNames, j++, n+1, "breaks");

  // Get information
  for ( unsigned int i=0; i<n; ++i ) {
    NUMERIC_POINTER(breaks)[i] = axis->GetBinLowEdge(i+1);
  }

  // Add the high edge
  NUMERIC_POINTER(breaks)[n] = axis->GetBinUpEdge(n);

  // Determine mids--
  SEXP mids = addNumericVector(data, dataNames, j++, n, "mids");

  // Get information
  for ( unsigned int i=0; i<n; ++i ) {
    NUMERIC_POINTER(mids)[i] = axis->GetBinCenter(i+1);
  }

  // Get name of axis
  SEXP xname = addCharVector(data, dataNames, j++, 1, "xname");
  SET_STRING_ELT( xname, 0, mkChar( axis->GetTitle() ) );
  
  // Done
  return j;
}

TH2F * histo_rebin(TH2F *h, const int nx, const Double_t *xbins, const int ny, const Double_t *ybins){

  // NB: the errors of the rebinned histogram are wrong!

  TH2F * h_new = new TH2F("h_new",h->GetTitle(),nx,xbins,ny,ybins);
  h_new->Sumw2();

  TAxis *xaxis = h->GetXaxis();
  TAxis *yaxis = h->GetYaxis();
  
  for (int jbin=1; jbin<=yaxis->GetNbins(); ++jbin){
    for (int ibin=1; ibin<=xaxis->GetNbins(); ++ibin){
      h_new->Fill(xaxis->GetBinCenter(ibin),yaxis->GetBinCenter(jbin),h->GetBinContent(ibin,jbin));
    }
  }
  
  return (TH2F*) h_new->Clone();
  
}

void dumpHistogram(TH1* histogram)
{
  std::cout << "<dumpHistogram>:" << std::endl;
  std::cout << " histogram: name = " << histogram->GetName() << ", title = " << histogram->GetTitle() << std::endl;
  std::cout << "  fillColor = " << histogram->GetFillColor() << ", fillStyle = " << histogram->GetFillStyle() << ","
	    << " lineColor = " << histogram->GetLineColor() << ", lineStyle = " << histogram->GetLineStyle() << ", lineWidth = " << histogram->GetLineWidth() << "," 
	    << " markerColor = " << histogram->GetMarkerColor() << ", markerStyle = " << histogram->GetMarkerStyle() << ", markerSize = " << histogram->GetMarkerSize() << std::endl;
  TAxis* xAxis = histogram->GetXaxis();
  int numBins = xAxis->GetNbins();
  for ( int iBin = 1; iBin <= numBins; ++iBin ) {
    std::cout << "bin #" << iBin << " (x = " << xAxis->GetBinCenter(iBin) << "): " << histogram->GetBinContent(iBin) << " +/- " << histogram->GetBinError(iBin) << std::endl;
  }
  std::cout << "integral = " << compIntegral(histogram, true, true) << std::endl;
}

//TH1D *makeDiffHistWptLog(TH1D* hData, TH1D* hFit, const TString name)
TH1D *makeDiffHistWptLog(TH1D* hData, TH1D* hFit, TH1D* hFitError, const TString name)
{
  //double WptBins[14]={1,7.5,12.5,17.5,24,30,40,50,70,110,150,190,250,600};
  double WptBins[14]={0.68,7.5,12.5,17.5,24,30,40,50,70,110,150,190,250,600};
  TH1D *hDiff = new TH1D(name,"",13,WptBins);
  TAxis *xaxis = hData->GetXaxis();
  for(int ibin=1; ibin<=hData->GetNbinsX(); ibin++) {
    //cout << "Wpt Data: " << hData->GetBinContent(ibin) << "\t sqrt(data): " << sqrt(hData->GetBinContent(ibin)) << "\t Wpt Fit: " << hFit->GetBinContent(ibin) << "\t fit error : " << hFitError->GetBinError(ibin) << endl;
    Double_t diff = (hData->GetBinContent(ibin)-hFit->GetBinContent(ibin)); // data-mc
    //Double_t err = sqrt(hData->GetBinContent(ibin));
    Double_t err = hFitError->GetBinError(ibin); // error from fit
    if(err==0) err= sqrt(hFit->GetBinContent(ibin));
    if(err>0) hDiff->Fill(xaxis->GetBinCenter(ibin),diff/err);
    else      hDiff->Fill(xaxis->GetBinCenter(ibin),0);
    cout << 
      "Wpt Data:\t " <<   hData->GetBinContent(ibin) << 
      "\t sqrt(data):\t " << sqrt(hData->GetBinContent(ibin)) << 
      "\t Wpt Fit:\t " << hFit->GetBinContent(ibin) << 
      "\t error:\t " << hFitError->GetBinError(ibin) <<
      "\t Wpt Data - Wpt Fit:\t " <<   hData->GetBinContent(ibin)-hFit->GetBinContent(ibin) << 
      "\tdiff/err: " << diff/err  <<
      endl;
    //cout << ibin <<  "\tdiff:   " << diff  <<  "\tfab(diff):   " << fab(diff)  << endl;
    //cout << ibin<< "\tdiff/err: " << diff/err  << endl;
  }

  //hDiff->GetYaxis()->SetTitleOffset(0.42);
  hDiff->GetYaxis()->SetTitleOffset(0.55);
  //hDiff->GetYaxis()->SetTitleSize(0.13);
  hDiff->GetYaxis()->SetTitleSize(0.10);
  hDiff->GetYaxis()->SetLabelSize(0.10);
  hDiff->GetYaxis()->SetNdivisions(104);
  hDiff->GetYaxis()->CenterTitle();
  hDiff->GetXaxis()->SetTitleOffset(1.2);
  hDiff->GetXaxis()->SetTitleSize(0.13);
  hDiff->GetXaxis()->SetLabelSize(0.12);
  //hDiff->GetXaxis()->CenterTitle();

  return hDiff;
}

 void fillHistograms(double mvaOutput, double pt, double eta, double Nvtx, double evtWeight)
 {    
   histogramPt_denominator_->Fill(pt, evtWeight);
   histogramEta_denominator_->Fill(eta, evtWeight);
   histogramNvtx_denominator_->Fill(Nvtx, evtWeight);
   bool passesCuts = (mvaOutput > mvaCut_);
   //std::cout << "passesCuts = " << passesCuts << std::endl;
   if ( passesCuts ) {
     histogramPt_numerator_->Fill(pt, evtWeight);
     histogramEta_numerator_->Fill(eta, evtWeight);
     histogramNvtx_numerator_->Fill(Nvtx, evtWeight);
   }
   double y = mvaOutput;
   TAxis* yAxis = histogramMVAoutput_vs_Pt_->GetYaxis();
   int binY = yAxis->FindBin(y);
   int numBinsY = yAxis->GetNbins();
   if ( binY <  1       ) binY = 1;
   if ( binY > numBinsY ) binY = numBinsY;
   double yWithinRange = yAxis->GetBinCenter(binY);
   histogramMVAoutput_vs_Pt_->Fill(pt, y, evtWeight);
   histogramPt_->Fill(pt, evtWeight);
 }

void DrawMeasurement( TH1 *h, Int_t nbin, Double_t val, Double_t err )
{
   TAxis *axis   = h->GetYaxis();
   
   Double_t dy;
   Double_t x1,y1,x2,y2, ymin;
   
   Int_t i = nbin;
   
   ymin = axis->GetBinCenter(i)+1;
   dy   = axis->GetBinWidth(i);
   x1   = val - err;
   y1   = ymin - 0.05*dy;
   x2   = val + err;
   y2   = ymin + 0.05*dy;
  
   TLine* line = new TLine;
   TLine* lbar = new TLine;
   TLine* rbar = new TLine;
   TMarker* m2 = new TMarker( val, ymin, 20 );
   Int_t color = 1;
   line->SetLineColor( color );
   line->SetLineWidth( 2 );
   lbar->SetLineColor( color );
   lbar->SetLineWidth( 2 );
   rbar->SetLineColor( color );
   rbar->SetLineWidth( 2 );
   m2->SetMarkerColor( color );
   m2->SetMarkerSize( 1.2 );
   m2->SetMarkerStyle( 21 );
   
   line->DrawLine( x1, ymin, x2, ymin );
   lbar->DrawLine( x1, y1, x1, y2 );
   rbar->DrawLine( x2, y1, x2, y2 );
   m2->Draw();
}

void RebinProfile(){

	f = new TFile("../../../rootfiles/MC/EPOS_PbPb_Ntrk_0.root");

	TH2D* h2;
	h2 = (TH2D*)f->Get("demo/scatterHist_effcorr");



	TCanvas* C = new TCanvas("fitscatterplot","EPOS Hydjet Scatterplot",1800,1200);

	C->Divide(3,2);
	C->cd(1);

	double newbins[8] = {-0.3,-0.075,-0.045,-0.015,0.015,0.045,0.075,0.3};
	//double newbins[6] = {-0.3,-0.045,-0.015,0.015,0.045,0.3};
	TH2D* hnew = new TH2D("rebin","rebin",7,newbins,1000,-0.3,0.3); 

	TAxis *xaxis = h2->GetXaxis();
	TAxis *yaxis = h2->GetYaxis();
	for (int k=1; k<=yaxis->GetNbins();k++) {
		for (int l=1; l<=xaxis->GetNbins();l++) {
			hnew->Fill(xaxis->GetBinCenter(l),yaxis->GetBinCenter(k),h2->GetBinContent(l,k));
		}
	}

	TProfile *prof = hnew->ProfileX();

	TF1* fit1 = new TF1("Linear fitting case 1", "[0]+x*[1]", -0.2, 0.2);
	fit1->SetLineColor(kBlue);
	fit1->SetLineStyle(1);
	fit1->SetLineWidth(3);
	prof->Fit(fit1,"RN0");

	TLatex* text2 = makeLatex(Form("slope : %.3f #pm %.3f",fit1->GetParameter(1),fit1->GetParError(1)),0.55,0.25) ;
	prof->Draw();

	text2->Draw("same");
	fit1->Draw("same");


	


/*
	for (Int_t i=0;i<2;i++) {
		for (Int_t j=0;j<3;j++) {


			h = (TH2D*)f->Get("demo/scatterHist_noeffcorr");
			h -> new TH2D(" ");


			corrTable[oc]->SetBinContent(eta+1, pt+1, c);

			

			double c1 = testload.getWeight(ptbins[pt], etabins[eta], ocbin[oc], "eff");
			double c2 = testload.getWeight(ptbins[pt], etabins[eta], ocbin[oc], "fak");
			double c = c1/(1-c2);

			if( oc == 0 && eta == 15 && pt == 20 ) cout << "1: " << c << endl;
			if( oc == 4 && eta == 15 && pt == 20 ) cout << "4: " << c << endl;
				//cout << "corr factor: " << c << endl; 
			corrTable[oc]->SetBinContent(eta+1, pt+1, c);


			c = c1/(1-c2);
			c2 = c3/(1-c1);
			c4 = c2/(1-c3);
			c5 = c6/(1-c7);

*/
		}

//void pi0_mfitpeak(char *FileName, char *HistName, Int_t etapi0flag) 
void pi0_mfitpeak(TH1F *mh1, Int_t etapi0flag, float xmin, float xmax, int npol,float res[],int posFlag, const char *dirName, const char *histName, float text_x, float text_y, const char *texName) 
  
{
  TGaxis::SetMaxDigits(3);


  // TVirtualFitter::SetDefaultFitter("Minuit");
  


  // This script attempts to fit any pi0 peak so that the freaking fit function would converge
  // currently background is fitted to a pol4 function and the peak by a gaussian;
  // results are not very nice
  // usage  .x pi0_mfitpeak.C++ ("pi0calib.root","minv_spb")
  // or eg.  .x pi0_mfitpeak.C ("../pi0anal/pi0ana_punorm.root","minv_spb",0)

  gROOT->Reset();
  //  gStyle->SetOptFit();
  //  gStyle->SetOptFit(0);
  //  gStyle->SetOptStat(0);
  //  gStyle->SetOptTitle(0);
Bool_t NOTE=1;
  if(NOTE) gStyle->SetCanvasBorderMode(0);

  gStyle->SetPadTopMargin(0.08);
  gStyle->SetPadBottomMargin(0.12);
  gStyle->SetPadLeftMargin(0.15);
  gStyle->SetPadRightMargin(0.08);

  mh1->GetXaxis()->SetRangeUser(xmin,xmax);
  
  Int_t highx=500;
  TCanvas *c2 = new TCanvas("c2","",200,10,highx,500);


    // cout<<FileName<<" "<<HistName<<endl;

  //     TFile f(FileName);
  //   TH1F *mh1 = (TH1F*) f.Get(HistName);
      mh1->SetMarkerStyle(20);
      mh1->SetMarkerSize(1.);
      mh1->SetStats(0); // 1/0 to set the stat box

   mh1->GetXaxis()->SetTitle("Invariant Mass of Photon Pairs (GeV/c^{2})");


   float binwidth = mh1->GetBinWidth(1);
   
   char *ytitle = new char[100];

   sprintf(ytitle,"Photon Pairs / %4.3f GeV/c^{2}",binwidth);
   

   mh1->GetYaxis()->SetTitle(ytitle);


   mh1->GetXaxis()->SetTitleSize(0.055);
   mh1->GetYaxis()->SetTitleSize(0.055);
   mh1->GetXaxis()->SetLabelSize(0.045);
   mh1->GetYaxis()->SetLabelSize(0.045);
   mh1->GetXaxis()->SetTitleOffset(0.90);
   mh1->GetXaxis()->CenterTitle();
   mh1->GetYaxis()->SetTitleOffset(1.32);
   

   // First work with the histogram and find the peak and fit ranges
   TAxis *xaxis = mh1->GetXaxis();
   Float_t binsiz= xaxis->GetBinCenter(3) - xaxis->GetBinCenter(2);
   Int_t nbins = xaxis->GetNbins(); 
   Float_t nevtperbin0[10000];
   Float_t errorbin0[10000];
   Float_t nevttot;
   Float_t maxbin=0; Int_t nmaxbin=0, nminbord=0, nmaxbord=nbins;
   
   for (Int_t nn=1; nn <= nbins; nn++)
     {
       nevtperbin0[nn] = mh1->GetBinContent(nn); 
       if(nevtperbin0[nn] > maxbin) { maxbin=nevtperbin0[nn]; nmaxbin=nn; }
       errorbin0[nn] = mh1->GetBinError(nn); 
       nevttot+=nevtperbin0[nn];
       if(nevtperbin0[nn] > 0 && nminbord == 0) nminbord=nn; 
       if(nevtperbin0[nn] == 0 && (nn > nminbord +10) && nmaxbord==0 && nminbord > 0) nmaxbord=nn; 
     }
   cout<<"Minbordl "<<nminbord<<" with events: "<<nevtperbin0[nminbord]<<endl;
   cout<<"Maxbordl "<<nmaxbord<<" with events: "<<nevtperbin0[nmaxbord]<<endl;
   nminbord+=0;
   nmaxbord-=0;
   Int_t nmin0=nminbord;
   while(nevtperbin0[nminbord] < nevtperbin0[nmaxbin]*0.025) nminbord++;
   while(nevtperbin0[nmaxbord] < nevtperbin0[nmaxbin]*0.025) nmaxbord--;
   // the above was just to get the info and low/high bins

   // Set the fit range ! This is for total fit !	 
   Float_t fitl=xmin;
     float fith=xmax;
   //     Float_t fitl=0.07, fith=0.2;// this works better for pileup
   //         Float_t fitl=0.08, fith=0.18;// this works even better for pileup
     //  if(etapi0flag == 1)
     // {
//.........这里部分代码省略.........

//.........这里部分代码省略.........
      for(unsigned int i=0; i<fr->J1_PseudoPoints.size(); i++)
	{
	  for(unsigned int j=0; j<fr->J3_PseudoPoints.size(); j++)
	    {
	      for(unsigned int k=0; k<fr->J4_PseudoPoints.size(); k++)
		{
		  double Rad  = FindR(fr->J1_PseudoPoints[i],fr->J3_PseudoPoints[j],fr->J4_PseudoPoints[k]);
		  double Xcen = FindX(fr->J1_PseudoPoints[i].X(),fr->J1_PseudoPoints[i].Y(),fr->J3_PseudoPoints[j].X(),fr->J3_PseudoPoints[j].Y(),fr->J4_PseudoPoints[k].X(),fr->J4_PseudoPoints[k].Y());
		  double Ycen = FindY(fr->J1_PseudoPoints[i].X(),fr->J1_PseudoPoints[i].Y(),fr->J3_PseudoPoints[j].X(),fr->J3_PseudoPoints[j].Y(),fr->J4_PseudoPoints[k].X(),fr->J4_PseudoPoints[k].Y());
		  Hough->Fill(Xcen,Ycen,Rad);
		  Houghxy->Fill(Xcen,Ycen);
		  Houghyr->Fill(Ycen,Rad);
		  Houghxr->Fill(Xcen,Rad);
		}
	    }
	}
    }
  
  //J2,J3,J4;
  if(fr->J2_PseudoPoints.size() >= 1 && fr->J3_PseudoPoints.size() >= 1 && fr->J4_PseudoPoints.size() >= 1)
    {
      for(unsigned int i=0; i<fr->J2_PseudoPoints.size(); i++)
	{
	  for(unsigned int j=0; j<fr->J3_PseudoPoints.size(); j++)
	    {
	      for(unsigned int k=0; k<fr->J4_PseudoPoints.size(); k++)
		{
		  double Rad  = FindR(fr->J2_PseudoPoints[i],fr->J3_PseudoPoints[j],fr->J4_PseudoPoints[k]);
		  double Xcen = FindX(fr->J2_PseudoPoints[i].X(),fr->J2_PseudoPoints[i].Y(),fr->J3_PseudoPoints[j].X(),fr->J3_PseudoPoints[j].Y(),fr->J4_PseudoPoints[k].X(),fr->J4_PseudoPoints[k].Y());
		  double Ycen = FindY(fr->J2_PseudoPoints[i].X(),fr->J2_PseudoPoints[i].Y(),fr->J3_PseudoPoints[j].X(),fr->J3_PseudoPoints[j].Y(),fr->J4_PseudoPoints[k].X(),fr->J4_PseudoPoints[k].Y());
		  Hough->Fill(Xcen,Ycen,Rad);
		  Houghxy->Fill(Xcen,Ycen);
		  Houghyr->Fill(Ycen,Rad);
		  Houghxr->Fill(Xcen,Rad);
		}
	    }
	}
    }

  /*
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!DO NOT ERASE!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  //Forming PseudoPoints:
  for(unsigned int ps=0; ps<fr->AllX.size(); ps++)
    {
      for(unsigned int s=0; s<fr->AllY.size(); s++)
	{
	   fr->PseudoPoints.push_back(TVector3(fr->AllX[ps].Cor(),fr->AllY[s].Cor(),0));
	}
    }
    
  //Full Combinatorial:
   for(unsigned int i=0; i<fr->PseudoPoints.size()-2; i++)
     {
       for(unsigned int j=i+1; j<fr->PseudoPoints.size()-1; j++)
 	{
 	  for(unsigned int k=j+1; k<fr->PseudoPoints.size(); k++)
 	    {
 	      double Rad = FindR(fr->PseudoPoints[i],fr->PseudoPoints[j],fr->PseudoPoints[k]);
 	      double Xcen = FindX(fr->PseudoPoints[i].X(),fr->PseudoPoints[i].Y(),fr->PseudoPoints[j].X(),fr->PseudoPoints[j].Y(),fr->PseudoPoints[k].X(),fr->PseudoPoints[k].Y());
 	      double Ycen = FindY(fr->PseudoPoints[i].X(),fr->PseudoPoints[i].Y(),fr->PseudoPoints[j].X(),fr->PseudoPoints[j].Y(),fr->PseudoPoints[k].X(),fr->PseudoPoints[k].Y());
 	      Hough->Fill(Xcen,Ycen,Rad);
 	      Houghxy->Fill(Xcen,Ycen);
 	      Houghyr->Fill(Ycen,Rad);
 	      Houghxr->Fill(Xcen,Rad);
	      }
 	}
     }
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  */      
  int xbin, ybin, zbin;
  xbin=0; ybin=0; zbin=0;
  Hough->GetMaximumBin(xbin, ybin, zbin);
  //cout << "xbin= " << xbin << " ybin= " << ybin << " zbin= " << zbin << endl;

  TAxis *xaxis = Hough->GetXaxis();  
  double Xrec=xaxis->GetBinCenter(xbin);
  cout<<"Xrec= "<< Xrec <<endl;

  TAxis *yaxis = Hough->GetYaxis();  
  double Yrec=yaxis->GetBinCenter(ybin);
  cout<<"Yrec= "<< Yrec <<endl;

  TAxis *zaxis = Hough->GetZaxis();  
  double Rrec=zaxis->GetBinCenter(zbin);
  cout<<"Rrec= "<< Rrec <<endl;

  //Drawing the founded ring and its center:
  TEllipse ring (Xrec,Yrec,Rrec,Rrec,1,360,0);
  ring.SetLineColor(2);
  ring.SetFillStyle(0);
  ring.DrawEllipse(Xrec,Yrec,Rrec,Rrec,0,360,0);
  HBD->Update();

  TMarker circumcenter (Xrec,Yrec,8);
  circumcenter.SetMarkerColor(2);
  circumcenter.SetMarkerSize(1);
  circumcenter.DrawMarker(Xrec,Yrec);
  HBD->Update();  
  
}

//.........这里部分代码省略.........
   leg->AddEntry(eMeanInter,"not reconstructed tracks 2 < p_{T}^{#pi^{#pm}} < 10 GeV","P");
   leg->Draw();

   c3->cd(2);
   TAxis* xaxis = gleak->GetXaxis();
   gleak->GetXaxis()->SetTitle("#eta");
   gleak->GetYaxis()->SetTitle("ratio = <E_{T for lost tracks}^{raw reco} / E_{T for found tracks}^{raw reco}>");
   xaxis->SetLimits(0.0,2.4);
   gleak->SetMarkerStyle(21);
   gleak->SetMaximum(1.0);
   //  leak->SetMinimum(0.7);
   gleak->SetMinimum(0.5);
   gleak->Draw("AP");
   TLatex *t = new TLatex();
   t->SetTextSize(0.08);
   //   t->DrawLatex(0.1,0.75,"<E_{T for lost tracks}^{raw reco} / E_{T for found tracks}^{raw reco}> for p_{T}^{#pi^{#pm}} >2 GeV");
   c3->SaveAs("eMean_vs_eta_pt2-10.gif");
   c3->SaveAs("eMean_vs_eta_pt2-10.eps");

   //  original distribtions
   setTDRStyle(0,0);
   gStyle->SetOptFit(0);
   TCanvas* c4 = new TCanvas("X","Y",1);
   hEnTrkNotInter[0]->GetYaxis()->SetTitle("");
   hEnTrkNotInter[0]->GetXaxis()->SetTitle("calo E_{T}^{raw reco} in cone 0.5/p_{T}^{MC}");
   Double_t scale = 1./hEnTrkNotInter[0]->Integral();
   hEnTrkNotInter[0]->Scale(scale);
   hEnTrkNotInter[0]->SetLineWidth(4);
   hEnTrkNotInter[0]->SetMaximum(0.14);
   //   hEnTrkNotInter[0]->SetMinimum(0.55);
   // Fitting
   Int_t binMax = hEnTrkNotInter[0]->GetMaximumBin();
   TAxis* xaxis = hEnTrkNotInter[0]->GetXaxis();
   Double_t binCenter = xaxis->GetBinCenter(binMax);
   Double_t rms = hEnTrkNotInter[0]->GetRMS();
   Double_t rFitMin = binCenter - 1.5 * rms; 
   Double_t rFitMax = binCenter + 1.5 * rms;
   hEnTrkNotInter[0]->Fit("gaus","","",rFitMin,rFitMax);
   TF1 *fit = hEnTrkNotInter[0]->GetFunction("gaus"); 
   fit->SetLineWidth(4);
   fit->SetLineStyle(2);
   fit->Draw("same");
   scale = 1./hEnTrkInter[0]->Integral();
   hEnTrkInter[0]->Scale(scale);
   hEnTrkInter[0]->SetLineWidth(2);
  // Fitting
   Int_t binMax = hEnTrkInter[0]->GetMaximumBin();
   TAxis* xaxis = hEnTrkInter[0]->GetXaxis();
   Double_t binCenter = xaxis->GetBinCenter(binMax);
   Double_t rms = hEnTrkNotInter[0]->GetRMS();
   Double_t rFitMin = binCenter - 1.5 * rms; 
   Double_t rFitMax = binCenter + 1.5 * rms;
   hEnTrkInter[0]->Fit("gaus","","same",rFitMin,rFitMax);
   TF1 *fit = hEnTrkInter[0]->GetFunction("gaus"); 
   fit->SetLineWidth(2);
   fit->SetLineStyle(2);
   fit->Draw("same");
   TLatex *t = new TLatex();
   t->SetTextSize(0.08);
   t->DrawLatex(0.2,0.9,"CMSSW_1_6_9");
   TLegend *leg = new TLegend(0.15,0.7,0.9,0.9,NULL,"brNDC");
   leg->SetFillColor(10);
   leg->AddEntry(hEnTrkNotInter[0],"recontructed tracks 2< p_{T}^{#pi^{#pm}} < 10 GeV, 0<#eta<0.2","L");
   leg->AddEntry(hEnTrkInter[0],"not reconstructed tracks 2< p_{T}^{#pi^{#pm}} < 10 GeV, 0<#eta<0.2","L");
   leg->Draw();

//.........这里部分代码省略.........
   TLatex * s1400 = new TLatex( 790, 580, "#tilde{q} (1400 GeV)" );
  s1400->SetTextAlign( 11 );
  s1400->SetTextAngle(-60);
  s1400->SetTextSize( 0.025 );
  s1400->SetTextColor( 16 );
  s1400->Draw();   

  /*TLatex * g400 = new TLatex( 1100, 140, "#tilde{g} (400 GeV)" );
  g400->SetTextAlign( 11 );
  g400->SetTextSize( 0.025 );
  g400->SetTextColor( 203 );
  g400->Draw();*/
  /*TLatex * g500 = new TLatex( 1000, 185, "#tilde{g} (500 GeV)" );
  g500->SetTextAlign( 11 );
  g500->SetTextSize( 0.025 );
  g500->SetTextColor( TColor::GetColor("#dddddd") );
  g500->Draw();*/
  TLatex * g600 = new TLatex( 1100, 225, "#tilde{g} (600 GeV)" );
  g600->SetTextAlign( 11 );
  g600->SetTextAngle(-4);
  g600->SetTextSize( 0.025 );
  g600->SetTextColor( 16 );
  g600->Draw();
  /*TLatex * g900 = new TLatex( 550, 380, "#tilde{g} (900 GeV)" );
  g900->SetTextAlign( 11 );
  g900->SetTextSize( 0.025 );
  g900->SetTextColor( TColor::GetColor("#dddddd") );
  g900->Draw();*/
  TLatex * g800 = new TLatex( 690, 330, "#tilde{g} (800 GeV)" );
  g800->SetTextAlign( 11 );
  g800->SetTextSize( 0.025 );
  g800->SetTextColor( 16 );
  //g800->Draw();
  TLatex * g1000 = new TLatex( 1400, 399, "#tilde{g} (1000 GeV)" );
  g1000->SetTextAlign( 11 );
  g1000->SetTextAngle(-5); 
  g1000->SetTextSize( 0.025 );
  g1000->SetTextColor( 16 );
  g1000->Draw();
  TLatex * g1200 = new TLatex( 1550, 489, "#tilde{g} (1200 GeV)" );
  g1200->SetTextAlign( 11 );
  g1200->SetTextAngle(-6); 
  g1200->SetTextSize( 0.025 );
  g1200->SetTextColor( 16 );
  //g1200->Draw();  
  TLatex * g1400 = new TLatex( 1650, 582, "#tilde{g} (1400 GeV)" );
  g1400->SetTextAlign( 11 );
  g1400->SetTextAngle(-6); 
  g1400->SetTextSize( 0.025 );
  g1400->SetTextColor( 16 );
  g1400->Draw();
  
  // island hacks
  if (true && channel==4) { // muon fixes 
    cout << "removing islands in muon channel ..." << endl;
    // contour line is drawn for values at 1.64485
    TAxis* ax = contour_obs->GetXaxis();
    TAxis* ay = contour_obs->GetYaxis();

    TH2F* contour_fix = contour_em1s;

    for (int xbin = 1; xbin <= contour_fix->GetNbinsX(); xbin++) {
      for (int ybin = 1; ybin <= contour_fix->GetNbinsY(); ybin++) {
	// island 1
	if ( ax->GetBinCenter( xbin) > 1350.  && ax->GetBinCenter( xbin) < 1500. && ay->GetBinCenter( ybin) < 130. && ay->GetBinCenter( ybin) > 89. ) {
	  cout << "Found spot here: " << xbin << " (" << ax->GetBinCenter( xbin)  << "), "
	       << ybin << " (" << ay->GetBinCenter( ybin) << "), "
	       << " value: " << contour_fix->GetBinContent(xbin,ybin) <<   endl;
	  cout << "   HACK : Setting above point by hand to 1.65 (!)" << endl;
	  if (contour_fix->GetBinContent(xbin,ybin)<1.65) contour_fix->SetBinContent(xbin, ybin, 1.66);
	}

      }
    }

  } 
  if (false && channel==1) { // electron

    cout << "removing islands in electron channel ..." << endl;
    // contour line is drawn for values at 1.64485
    TAxis* ax = contour_obs->GetXaxis();
    TAxis* ay = contour_obs->GetYaxis();

    contour_em1s

    for (int xbin = 1; xbin <= contour_obs->GetNbinsX(); xbin++) {
      for (int ybin = 1; ybin <= contour_obs->GetNbinsY(); ybin++) {
	// island 2
	if ( ax->GetBinCenter( xbin) > 420. && ax->GetBinCenter( xbin) < 480. &&
	     ay->GetBinCenter( ybin) > 140. && ay->GetBinCenter( ybin) < 160. ) {
	  cout << "Found spot here: " << xbin << " (" << ax->GetBinCenter( xbin)  << "), "
	       << ybin << " (" << ay->GetBinCenter( ybin) << "), "
	       << " value: " << contour->GetBinContent(xbin,ybin) <<   endl;
	  cout << "   HACK : Setting above point by hand to 1.50 (!)" << endl;
	  contour->SetBinContent(xbin, ybin, 1.50);
	}
      }
    }

  }

//.........这里部分代码省略.........
  }

/*
  ///////////////////
  // Fit pT-function:
  Char_t pt_fit[] = "[0]+[1]*x+[2]*x^2";
  Char_t fitopt[] = "QMR";
  Double_t xmin = 150;//45;
  Double_t xmax = 1050;//275;
  Double_t pt[2][3];
  TF1 *ptfit = new TF1("ptfit",pt_fit,xmin,xmax);

  p0->Fit(ptfit,fitopt,"sames",xmin,xmax);
  pt[0][0] = ptfit->GetParameter(0);
  pt[0][1] = ptfit->GetParameter(1);
  pt[0][2] = ptfit->GetParameter(2);

  p1->Fit(ptfit,fitopt,"sames",xmin,xmax);
  pt[1][0] = ptfit->GetParameter(0);
  pt[1][1] = ptfit->GetParameter(1);
  pt[1][2] = ptfit->GetParameter(2);
*/

  ////////////////
  // Closure test:
  // ratio = 1 + (pt[0][0]+pt[0][1]*pt+pt[0][2]*pt^2)*KOW + (pt[1][0]+pt[1][1]*pt+[1][2]*pt^2)*KOW^2
  TH1* test = new TH1D("test","Test aQGC Photon ET",pv.NBINS,pv.MINRange,pv.MAXRange);
  TH1* test2 = new TH1D("test2","Test aQGC Photon ET",pv.NBINS,pv.MINRange,pv.MAXRange);

  for(Int_t j=1;j<=pv.NBINS;j++){

     Double_t value,bincent;
     TAxis* xaxis = test->GetXaxis();
     bincent = xaxis->GetBinCenter(j);

     value = p[0][j-1]+p[1][j-1]*(8E-11)+p[2][j-1]*pow(8E-11,2);
     test->SetBinContent(j,value);

     test2->SetBinContent(j,p[0][j-1]+p[1][j-1]*(-8E-11)+p[2][j-1]*pow(-8E-11,2));

  }

//  test->SetBinContent(8,1+p[0][7]*(8E-11)+p[1][7]*pow(8E-11,2));
//  test2->SetBinContent(8,1+p[0][7]*(-8E-11)+p[1][7]*pow(-8E-11,2));

  TH1D* test_r = (TH1D*)test->Clone("test_r");
  TH1D* test2_r = (TH1D*)test2->Clone("test2_r");

  test_r->Divide(signal_lt0_80);
  test2_r->Divide(signal_lt0_m80);

  test->Multiply(th1wwa);
  test2->Multiply(th1wwa);
  test->Add(th1wwa,-1);
  test2->Add(th1wwa,-1);


  /////////////////////
  // Simulate new aQGC:
  TH1* signal_lt0_18 = new TH1D("signal_lt0_18","signal_lt0_1p8",pv.NBINS,pv.MINRange,pv.MAXRange);
  TH1* signal_lt0_20 = new TH1D("signal_lt0_20","signal_lt0_2p0",pv.NBINS,pv.MINRange,pv.MAXRange);
  TH1* signal_lt0_22 = new TH1D("signal_lt0_22","signal_lt0_2p2",pv.NBINS,pv.MINRange,pv.MAXRange);
  TH1* signal_lt0_24 = new TH1D("signal_lt0_24","signal_lt0_2p4",pv.NBINS,pv.MINRange,pv.MAXRange);
  TH1* signal_lt0_26 = new TH1D("signal_lt0_26","signal_lt0_2p6",pv.NBINS,pv.MINRange,pv.MAXRange);
  TH1* signal_lt0_28 = new TH1D("signal_lt0_28","signal_lt0_2p8",pv.NBINS,pv.MINRange,pv.MAXRange);
  TH1* signal_lt0_32 = new TH1D("signal_lt0_32","signal_lt0_3p2",pv.NBINS,pv.MINRange,pv.MAXRange);

//.........这里部分代码省略.........
  CPlot* plotWMptDiffLog;

//W plus pt distribution
  TH1D* hWptMC_p = (TH1D*)hDYToTauTauP->Clone("hWptMC_p");
  hWptMC_p->Add(hTTJetsP);
  hWptMC_p->Add(hWToTauNuP);
  hWptMC_p->Add(hDYToMuMuP);
  hWptMC_p->Add(hQCDWPpt);
  hWptMC_p->Add(hSigWPpt);
 
  TH1D* hWPptDiffLog;
  TH1D* hWMptDiffLog;

  //double WptBinsLog[14]={1,7.5,12.5,17.5,24,30,40,50,70,110,150,190,250,600};
  double WptBinsLog[14]={0.68,7.5,12.5,17.5,24,30,40,50,70,110,150,190,250,600};
  double x1,x2,x3,x4,x5,x6,x7,err;
  TH1D *hDYToTauTauLogP = new TH1D("hDYToTauTauLogP","",13,WptBinsLog);
  TH1D *hTTJetsLogP     = new TH1D("hTTJetsLogP","",13,WptBinsLog);
  TH1D *hWToTauNuLogP   = new TH1D("hWToTauNuLogP","",13,WptBinsLog);
  TH1D *hDYToMuMuLogP   = new TH1D("hDYToMuMuLogP","",13,WptBinsLog);
  TH1D *hQCDWptLogP     = new TH1D("hQCDWptLogP","",13,WptBinsLog);
  TH1D *hSigWptLogP     = new TH1D("hSigWptLogP","",13,WptBinsLog);
  TH1D *hdataWptLogP    = new TH1D("hdataWptLogP","",13,WptBinsLog);
  TAxis *xaxis = hSigWPpt->GetXaxis();
  for (int i=1; i<=hSigWPpt->GetNbinsX(); i++){
    x1=hDYToTauTauP->GetBinContent(i);
    x2=hTTJetsP->GetBinContent(i);
    x3=hWToTauNuP->GetBinContent(i);
    x4=hDYToMuMuP->GetBinContent(i);
    x5=hQCDWPpt->GetBinContent(i);
    x6=hSigWPpt->GetBinContent(i);
    x7=hdataWPpt->GetBinContent(i);
    err = hdataWPpt->GetBinError(i);
    hDYToTauTauLogP->Fill(xaxis->GetBinCenter(i),x1);
    hTTJetsLogP->Fill(xaxis->GetBinCenter(i),x2);
    hWToTauNuLogP->Fill(xaxis->GetBinCenter(i),x3);
    hDYToMuMuLogP->Fill(xaxis->GetBinCenter(i),x4);
    hQCDWptLogP->Fill(xaxis->GetBinCenter(i),x5);
    hSigWptLogP->Fill(xaxis->GetBinCenter(i),x6);
    hdataWptLogP->Fill(xaxis->GetBinCenter(i),x7);
    hdataWptLogP->SetBinError(i,err);
  }
  
  TH1D* hEwkWPptLog = (TH1D*)hDYToTauTauLogP->Clone("hEwkWPptLog");
  hEwkWPptLog->Add(hTTJetsLogP);
  hEwkWPptLog->Add(hWToTauNuLogP);
  hEwkWPptLog->Add(hDYToMuMuLogP);
  
  TString plotName = "FitWDistribution_MuonPLog";
  if (filetype == "Electron")
    plotName = "FitWDistribution_ElePLog";
  plotWPptLog=new CPlot(plotName,"","p_{T}^{W} [GeV]","Events");
  plotWPptLog->setOutDir(CPlot::sOutDir);
  plotWPptLog->AddHist1D(hdataWptLogP,"#font[42]{Data}","E");
  plotWPptLog->AddToStack(hEwkWPptLog,"#font[42]{EW+t#bar{t}}",fillcolorEWK,linecolorEWK);
  plotWPptLog->AddToStack(hQCDWptLogP,"#font[42]{QCD}",fillcolorQCD,linecolorQCD);
  if(filetype == "Muon")
    plotWPptLog->AddToStack(hSigWptLogP,"#font[42]{W^{+}#rightarrow #mu^{+}#nu}",fillcolorW,linecolorW);
  if(filetype == "Electron")
    plotWPptLog->AddToStack(hSigWptLogP,"#font[42]{W^{+}#rightarrow e^{+}#nu}",fillcolorW,linecolorW);
  plotWPptLog->SetLegend(0.67,0.70,.98,0.88);
  //plotWPptLog->SetYRange(0.25,1.4*(hWptMC_p->GetMaximum()));
  plotWPptLog->SetYRange(1,50*(hWptMC_p->GetMaximum()));
  plotWPptLog->SetLogx();
  plotWPptLog->SetLogy();
  plotWPptLog->AddTextBox(CMStext,0.14,0.90,0.24,0.97,0);

//.........这里部分代码省略.........
  //J1,J3,J4:
  if(fr->J1_PseudoPoints.size() >= 1 && fr->J3_PseudoPoints.size() >= 1 && fr->J4_PseudoPoints.size() >= 1)
    {
      for(unsigned int i=0; i<fr->J1_PseudoPoints.size(); i++)
	{
	  for(unsigned int j=0; j<fr->J3_PseudoPoints.size(); j++)
	    {
	      for(unsigned int k=0; k<fr->J4_PseudoPoints.size(); k++)
		{
		  double Rad  = RFindR(fr->J1_PseudoPoints[i],fr->J3_PseudoPoints[j],fr->J4_PseudoPoints[k]);
		  double Xcen = RFindX(fr->J1_PseudoPoints[i].X(),fr->J1_PseudoPoints[i].Y(),fr->J3_PseudoPoints[j].X(),fr->J3_PseudoPoints[j].Y(),fr->J4_PseudoPoints[k].X(),fr->J4_PseudoPoints[k].Y());
		  double Ycen = RFindY(fr->J1_PseudoPoints[i].X(),fr->J1_PseudoPoints[i].Y(),fr->J3_PseudoPoints[j].X(),fr->J3_PseudoPoints[j].Y(),fr->J4_PseudoPoints[k].X(),fr->J4_PseudoPoints[k].Y());
		  HoughF->Fill(Xcen,Ycen,Rad);
		}
	    }
	}
    }
  
  //J2,J3,J4;
  if(fr->J2_PseudoPoints.size() >= 1 && fr->J3_PseudoPoints.size() >= 1 && fr->J4_PseudoPoints.size() >= 1)
    {
      for(unsigned int i=0; i<fr->J2_PseudoPoints.size(); i++)
	{
	  for(unsigned int j=0; j<fr->J3_PseudoPoints.size(); j++)
	    {
	      for(unsigned int k=0; k<fr->J4_PseudoPoints.size(); k++)
		{
		  double Rad  = RFindR(fr->J2_PseudoPoints[i],fr->J3_PseudoPoints[j],fr->J4_PseudoPoints[k]);
		  double Xcen = RFindX(fr->J2_PseudoPoints[i].X(),fr->J2_PseudoPoints[i].Y(),fr->J3_PseudoPoints[j].X(),fr->J3_PseudoPoints[j].Y(),fr->J4_PseudoPoints[k].X(),fr->J4_PseudoPoints[k].Y());
		  double Ycen = RFindY(fr->J2_PseudoPoints[i].X(),fr->J2_PseudoPoints[i].Y(),fr->J3_PseudoPoints[j].X(),fr->J3_PseudoPoints[j].Y(),fr->J4_PseudoPoints[k].X(),fr->J4_PseudoPoints[k].Y());
		  HoughF->Fill(Xcen,Ycen,Rad);
		}
	    }
	}
    }

  /*
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!DO NOT ERASE!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  //Forming PseudoPoints:
  for(unsigned int ps=0; ps<fr->AllX.size(); ps++)
    {
      for(unsigned int s=0; s<fr->AllY.size(); s++)
	{
	   fr->PseudoPoints.push_back(TVector3(fr->AllX[ps].Cor(),fr->AllY[s].Cor(),0));
	}
    }
    
  //Full Combinatorial:
   for(unsigned int i=0; i<fr->PseudoPoints.size()-2; i++)
     {
       for(unsigned int j=i+1; j<fr->PseudoPoints.size()-1; j++)
 	{
 	  for(unsigned int k=j+1; k<fr->PseudoPoints.size(); k++)
 	    {
 	      double Rad = RFindR(fr->PseudoPoints[i],fr->PseudoPoints[j],fr->PseudoPoints[k]);
 	      double Xcen = RFindX(fr->PseudoPoints[i].X(),fr->PseudoPoints[i].Y(),fr->PseudoPoints[j].X(),fr->PseudoPoints[j].Y(),fr->PseudoPoints[k].X(),fr->PseudoPoints[k].Y());
 	      double Ycen = RFindY(fr->PseudoPoints[i].X(),fr->PseudoPoints[i].Y(),fr->PseudoPoints[j].X(),fr->PseudoPoints[j].Y(),fr->PseudoPoints[k].X(),fr->PseudoPoints[k].Y());
 	      HoughF->Fill(Xcen,Ycen,Rad);
 	      HoughFxy->Fill(Xcen,Ycen);
 	      HoughFyr->Fill(Ycen,Rad);
 	      HoughFxr->Fill(Xcen,Rad);
	      }
 	}
     }
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  */      
  int xbin, ybin, zbin;
  xbin=0; ybin=0; zbin=0;
  HoughF->GetMaximumBin(xbin, ybin, zbin);
  // cout << "xbin= " << xbin << " ybin= " << ybin << " zbin= " << zbin << endl;

  TAxis *xaxis = HoughF->GetXaxis();  
  double Xrec=xaxis->GetBinCenter(xbin);
  //cout<<"Xrec= "<< Xrec <<endl;

  TAxis *yaxis = HoughF->GetYaxis();  
  double Yrec=yaxis->GetBinCenter(ybin);
  //cout<<"Yrec= "<< Yrec <<endl;

  TAxis *zaxis = HoughF->GetZaxis();  
  double Rrec=zaxis->GetBinCenter(zbin);
  //cout<<"Rrec= "<< Rrec <<endl;

   double Xerr = fr->CX[0].Cor()-Xrec;  
   double Yerr = fr->CY[0].Cor()-Yrec;  
   double Rerr = fr->Ring[0].R()-Rrec;  

//  if(fabs(Xerr) < 3*fr->Ring[0].SigmaX() && fabs(Yerr) < 3*fr->Ring[0].SigmaY() && fabs(Rerr) < 3*fr->Ring[0].SigmaR()) fr->Nsuc.push_back(1);
       
//    double Xdif = fabs(fr->Ring[0].Xmean())-fabs(Xerr); //cout<<"Xdif= " <<Xdif<<endl;
//    double Ydif = fabs(fr->Ring[0].Ymean())-fabs(Yerr); //cout<<"Ydif= " <<Ydif<<endl;
//    double Rdif = fabs(fr->Ring[0].Rmean())-fabs(Rerr); //cout<<"Rdif= " <<Rdif<<endl;

//    if(fabs(Xdif) < 3*fr->Ring[0].SigmaX() && fabs(Ydif) < 3*fr->Ring[0].SigmaY() && fabs(Rdif) < 3*fr->Ring[0].SigmaR()) fr->Nsuc_off.push_back(1);

   xerr->Fill(Xerr);
   yerr->Fill(Yerr);
   rerr->Fill(Rerr);
   
}