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

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

TH1 *PullHisto(const TList *list, const char *name, Int_t min, Int_t max,
               Double_t &mean)
{
  THnSparse *hs = list->FindObject(name);
  if (!hs) return 0;
  TAxis *atmp = hs->GetAxis(1);
  atmp->SetRange(min, max);
  // !!!!!!!!!!!!!!!!!!!!
  hs->GetAxis(2)->SetRangeUser(-0.5, 0.5);
  TH1 *hfin = hs->Projection(0);
  hfin->SetTitle(Form("p_{t} #in (%4.2f, %4.2f) GeV/c",
                      atmp->GetBinLowEdge(min),
                      atmp->GetBinLowEdge(max) + atmp->GetBinWidth(max)));
  mean = atmp->GetBinLowEdge(min) +
    (atmp->GetBinLowEdge(max) + atmp->GetBinWidth(max) -
     atmp->GetBinLowEdge(min))/2.0;
  // !!!!!!!!!!!!!!!!!!!!
  return hfin;//->Rebin();
}

/** 
 * Make a GraphSysErr object 
 * 
 * @param o         Output stream
 * @param h         Histogram
 * @param c1        Lower centrality bound
 * @param c2        Upper centrality bound 
 * @param reweigh   True if reweighed 
 * @param fac       Scaling factor 
 */
void MakeGSE(std::ostream& o, const TH1* h,
	     Double_t c1, Double_t c2, Bool_t reweigh, Double_t fac=1)
{
  // These are correlated
  Double_t bgMid   = CSysEval(c2, 100*0.02,  100*0.002);
  Double_t cSys    = CSysEval(c2, 100*0.005, 100*0.075);
  Double_t strMid  = 100*0.005;
  Double_t strFwd  = 100*0.05;
  Double_t pidMid  = 100*0;
  Double_t pidFwd  = 100*0.01;
  // Double_t bgSys = (1-c1/100)*(2-0.2) + 0.2;
  // Double_t cSys  = TMath::Power(c1/100,2)*(7.5-0.5) + 0.5;
  
  o << "*dataset:\n"
    << "*dscomment: The pseudo-rapidity density of charged particle\n"
    << "*reackey: PB PB --> CHARGED X\n"
    << "*obskey: DN/DETARAP\n"
    << "*qual: CENTRALITY IN PCT : " << c1 << " TO " << c2 << "\n"
    << "*qual: SQRT(S)/NUCLEON IN GEV : 5023\n";
  if (!reweigh) {
    o << "*dserror: " << strMid << " PCT : Weak decays\n"
      << "*dserror: " << pidMid << " PCT : Particle composition\n"
      << "*dserror: " << bgMid  << " PCT : Background subtraction\n";
  }
  o << "*dserror: 1   PCT : pT extrapolation\n"
    << "*dserror: " << cSys << " PCT : Centrality\n"
    << "*xheader: ETARAP\n"
    << "*yheader: DN/DETARAP\n"
    << "*data: x : y" << std::endl;
  // Define points
  TAxis* xa = h->GetXaxis();
  for (Int_t i = 0; i < h->GetNbinsX(); i++) {
    Int_t    j  = i+1;
    Double_t cc = h->GetBinContent(j)*fac;
    Double_t x  = h->GetXaxis()->GetBinCenter(j);
    Double_t ex = h->GetXaxis()->GetBinWidth(j)/2;
    Double_t xo = TMath::Abs(x)+ex;
    if (cc < 1e-8) continue;
    o << " " << xa->GetBinLowEdge(j) << " TO " << xa->GetBinUpEdge(j)
      << "; " << cc << " +-" << h->GetBinError(j)
      << " (DSYS=" << 0.01+0.01*TMath::Power(xo/2,2) << " PCT:Acceptance";
    if (reweigh) {
      Double_t bg     = EtaSysEval(xo, bgMid, 3*bgMid);
      Double_t pid    = EtaSysEval(xo, pidMid, pidFwd);
      Double_t str    = EtaSysEval(xo, strMid, strFwd);
      o << ",DSYS=" << bg  << " PCT:Background subtraction"
	<< ",DSYS=" << pid << " PCT:Particle composition"
	<< ",DSYS=" << str << " PCT:Weak decay";
      
    }
    o << ");"
      << std::endl;
  }
  o << "*dataend:\n" << std::endl;
}

double integralInRange( TH1F * var_h, float lower, float upper ){
  TAxis *axis = var_h->GetXaxis();
  int bmin = axis->FindBin(lower);
  int bmax = axis->FindBin(upper);
  double integr = var_h->Integral(bmin,bmax);
  integr -= ( (var_h->GetBinContent(bmin)) * (lower - (axis->GetBinLowEdge(bmin))) )/
    axis->GetBinWidth(bmin);
  integr -= ( (var_h->GetBinContent(bmax)) * ( (axis->GetBinUpEdge(bmax)) - upper) )/
    axis->GetBinWidth(bmax);
  return integr;
}

void GetInt(TH1F *hzj,float xmin,float xmax,double& integralfull, double& integral)
{
  TAxis *axis = hzj->GetXaxis();
  int bmin = axis->FindBin(xmin); 
  int bmax = axis->FindBin(xmax); 
  integralfull = hzj->Integral();
  integral = hzj->Integral(bmin,bmax);
  integral -= hzj->GetBinContent(bmin)*(xmin-axis->GetBinLowEdge(bmin))/axis->GetBinWidth(bmin);
  integral -= hzj->GetBinContent(bmax)*(axis->GetBinUpEdge(bmax)-xmax)/axis->GetBinWidth(bmax);
  fout<<"Full Integral: "<<integralfull<<endl;
}

double integrate(TH1F* h, double xmin, double xmax){
  TAxis *axis = h->GetXaxis();
  int bmin = axis->FindBin(xmin); //in your case xmin=-1.5
  int bmax = axis->FindBin(xmax); //in your case xmax=0.8
  double integral = h->Integral(bmin,bmax);
  integral -= h->GetBinContent(bmin)*(xmin-axis->GetBinLowEdge(bmin))/axis->GetBinWidth(bmin);
  integral -= h->GetBinContent(bmax)*(axis->GetBinUpEdge(bmax)-xmax)/axis->GetBinWidth(bmax);

  return integral;

}

//______________________________________________________________________________
TEveCaloDataVec* MakeVecData(Int_t ncells=0)
{
   // Example how to fill data when bins can be iregular.
   // If ncells = 0 (default) whole histogram is taken,
   // otherwise just ncells cells around the maximum.

   TFile::SetCacheFileDir(".");
   TFile* hf = TFile::Open(histFile, "CACHEREAD");
   TH2F* h1 = (TH2F*)hf->Get("ecalLego");
   TH2F* h2 = (TH2F*)hf->Get("hcalLego");

   TEveCaloDataVec* data = new TEveCaloDataVec(2);
   data->RefSliceInfo(0).Setup("ECAL", 0.3, kRed);
   data->RefSliceInfo(1).Setup("HCAL", 0.1, kBlue);

   TAxis *ax =  h1->GetXaxis();
   TAxis *ay =  h1->GetYaxis();

   Int_t xm = 1, xM = ax->GetNbins();
   Int_t ym = 1, yM = ay->GetNbins();
   if (ncells != 0)
   {
      Int_t cx, cy, cz;
      h1->GetMaximumBin(cx, cy, cz);
      xm = TMath::Max(xm, cx-ncells);
      xM = TMath::Min(xM, cx+ncells);
      ym = TMath::Max(ym, cy-ncells);
      yM = TMath::Min(yM, cy+ncells);
   }

   // Take every second cell and set a random size.
   for(Int_t i=xm; i<=xM; i+=2)
   {
      for(Int_t j=ym; j<=yM; j+=2)
      {
         if ( (i+j) % 3)
         {
            data->AddTower(ax->GetBinLowEdge(i), ax->GetBinUpEdge(i),
                           ay->GetBinLowEdge(j), ay->GetBinUpEdge(j));
            data->FillSlice(0, h1->GetBinContent(i, j));
            data->FillSlice(1, h2->GetBinContent(i, j));
         }
         else
         {
            data->AddTower(ax->GetBinLowEdge(i),
                           2 * ax->GetBinWidth(i) + ax->GetBinLowEdge(i),
                           ay->GetBinLowEdge(j),
                           2 * ay->GetBinWidth(j) + ay->GetBinLowEdge(j));
            data->FillSlice(0, h2->GetBinContent(i, j));
            data->FillSlice(1, h2->GetBinContent(i, j));
         }
      }
   }

   data->SetEtaBins(ax);
   data->SetPhiBins(ay);
   data->DataChanged();
   return data;
}

void GetIntg(TH1F *hzj,float xmin,float xmax,double& count1, double& count2)
{
 
  TAxis *axis = hzj->GetXaxis();
  int bmin = axis->FindBin(xmin); 
  int bmax = axis->FindBin(xmax); 
  double integralfull = hzj->Integral();
  double integral = hzj->Integral(bmin,bmax);
  integral -= hzj->GetBinContent(bmin)*(xmin-axis->GetBinLowEdge(bmin))/axis->GetBinWidth(bmin);
  integral -= hzj->GetBinContent(bmax)*(axis->GetBinUpEdge(bmax)-xmax)/axis->GetBinWidth(bmax);
  count1=count1+integralfull;
  count2=count2+integral;
  // cout<<count1<<'\t'<<count2<<endl;

}

void Bin2DTree::constrainSplit(int axis, double& cut, bool& veto)
/*****************************************************************/
{
    if(m_gridConstraint && !m_vetoSplitXY[axis])
    {
        TAxis* gridAxis = NULL;
        if(axis==0)
        {
            gridAxis = m_gridConstraint->GetXaxis();
        }
        else
        {
            gridAxis = m_gridConstraint->GetYaxis();
        }
        // Find the closest grid constraint for the cut
        // And modify the cut according to this constraint
        int b   = gridAxis->FindBin(cut);
        double low = gridAxis->GetBinLowEdge(b);
        double up  = gridAxis->GetBinUpEdge(b);
        if(fabs(up-cut)<fabs(cut-low))
        {
            cut = up;
            // If the constrained cut is outside the bin boundaries, try the other grid constraint
            if(cut>=getBinBoundaries()[axis].second)
            {
                cut = low;
            }
        }
        else
        {
            cut = low;
            // If the constrained cut is outside the bin boundaries, try the other grid constraint
            if(cut<=getBinBoundaries()[axis].first)
            {
                cut = up;
            }
        }
        //  If the constrained cut is still outside the bin boundaries, veto this bin and axis
        if(cut<=getBinBoundaries()[axis].first || cut>=getBinBoundaries()[axis].second)
        {
            m_vetoSplitXY[axis] = true;
        }
    }
    veto = m_vetoSplitXY[axis];
}

void rebin() {
   //create a fix bin histogram
   TH1F *h = new TH1F("h","test rebin",100,-3,3);
   Int_t nentries = 1000;
   h->FillRandom("gaus",nentries);
   Double_t xbins[1001];
   Int_t k=0;
   TAxis *axis = h->GetXaxis();
   for (Int_t i=1;i<=100;i++) {
      Int_t y = (Int_t)h->GetBinContent(i);
      if (y <=0) continue;
      Double_t dx = axis->GetBinWidth(i)/y;
      Double_t xmin = axis->GetBinLowEdge(i);
      for (Int_t j=0;j<y;j++) {
         xbins[k] = xmin +j*dx;
         k++;
      }
   }
   xbins[k] = axis->GetXmax();
   //create a variable bin-width histogram out of fix bin histogram
   //new rebinned histogram should have about 10 entries per bin
   TH1F *hnew = new TH1F("hnew","rebinned",k,xbins);
   hnew->FillRandom("gaus",10*nentries);

   //rebin hnew keeping only 50% of the bins
   Double_t xbins2[501];
   Int_t kk=0;
   for (Int_t j=0;j<k;j+=2) {
      xbins2[kk] = xbins[j];
      kk++;
   }
   xbins2[kk] = xbins[k];
   TH1F *hnew2 = (TH1F*)hnew->Rebin(kk,"hnew2",xbins2);

   //draw the 3 histograms
   TCanvas *c1 = new TCanvas("c1","c1",800,1000);
   c1->Divide(1,3);
   c1->cd(1);
   h->Draw();
   c1->cd(2);
   hnew->Draw();
   c1->cd(3);
   hnew2->Draw();
}

TH1D* TH2toTH1(TH2D* h2, Int_t massbin)
{
	TString massbinname = (TString)_s(massbin);
	TString origname = (TString)h2->GetName();
	// origname.ReplaceAll("hg4Mass_template_","");
	TString name = origname+"_mll_"+massbinname;
	//------------------------------------------------------------
	const Int_t nbins = h2->GetNbinsY();
	Double_t bins[nbins+1];
	TAxis* yaxis = (TAxis*)h2->GetYaxis();
	for(int i=0 ; i<nbins ; i++) bins[i] = yaxis->GetBinLowEdge(i+1);
	bins[nbins] = yaxis->GetBinUpEdge(nbins);
	//------------------------------------------------------------
	TH1D* h1 = new TH1D(name,"#frac{dN}{dg^{4}};g^{4};dN/dg^{4}",nbins,bins);
	for(Int_t bin=1 ; bin<=nbins ; bin++)
	{
		h1->SetBinContent(bin, h2->GetBinContent(massbin+bins2chop,bin));
	}
	return h1;
}

void th22mama(TH2* m, const char* filename, const char* comment="none")
{
    char tmp[128];
    time_t now = time(0);
    strftime(tmp, sizeof(tmp), "%Y-%m-%d %T", localtime(&now)); // not the 23-Mar-07 16:02:34 format
    
    TAxis *xax = m->GetXaxis(), *yax = m->GetYaxis();
    
    const int nx = std::min(xax->GetNbins(), 2048);
    const int ny = yax->GetNbins();
    std::cout << "matrix is " << nx << 'x' << ny << "; comment='" << comment << "'" << std::endl;
    
    std::ofstream mama(filename);
    mama << "!FILE=Disk \n"
    "!KIND=Matrix \n"
    "!LABORATORY=Oslo Cyclotron Laboratory (OCL) \n"
    "!EXPERIMENT=siri2root \n"
    "!COMMENT=" << comment << "\n"
    "!TIME=" << tmp << "\n"
    "!CALIBRATION EkeV=6";
    const float cal[6] = { xax->GetBinLowEdge(1), xax->GetBinWidth(1), 0,
        yax->GetBinLowEdge(1), yax->GetBinWidth(1), 0 };
    for(int i=0; i<6; ++i) {
        snprintf(tmp, sizeof(tmp), ",%13.6E", cal[i]);
        mama << tmp;
    }
    mama << "\n!PRECISION=16\n";
    snprintf(tmp, sizeof(tmp), "!DIMENSION=2,0:%4d,0:%4d\n", nx-1, ny-1);
    mama << tmp;
    snprintf(tmp, sizeof(tmp), "!CHANNEL=(0:%4d,0:%4d)\n", nx-1, ny-1);
    mama << tmp;
    
    for(int iy=1; iy<=ny; ++iy) {
        for(int ix=1; ix<=nx; ++ix)
            mama << m->GetBinContent(ix, iy) << ' ';
        mama << "\n";
    }
    mama << "!IDEND=" << endl << endl;
    
    mama.close();
}

void compare()
{
	msglvl[DBG] = SILENT;
	msglvl[INF] = VISUAL;
	msglvl[WRN] = VISUAL;
	msglvl[ERR] = VISUAL;
	msglvl[FAT] = VISUAL;

	TDirectory* oldDir = gDirectory;

	style();

	TFile* f2ds = NULL;
	TFile* f2d  = NULL;
	TH2D* h2ds  = NULL;
	TH2D* h2d   = NULL;
	TH1D* h1s   = NULL;
	TH1D* h1    = NULL;
	TH1D* h1dy  = NULL;

	Int_t g4bin = 50;
	TString model = "ZP";

	int imassmin = 330;
	int imassmax = 4930;
	int dimass   = 100;
	for(int imass=imassmin ; imass<=imassmax ; imass+=dimass)
	{
		TString mass = (TString)_s(imass);
		_INFO("enter "+(string)mass);

		TString modeltmp = (model=="ZP") ? "P" : "KK";

		oldDir->cd();
		f2ds = new TFile("~yiftahsi/KK_analysis_t301/2dtemplates_FEB1/templates2d_Xmass"+mass+".root","READ");
		h2ds = (TH2D*)f2ds->Get("h2"+model+"/Z"+modeltmp+"_mXp"+mass)->Clone();

		f2d  = new TFile("~yiftahsi/KK_analysis_t301/2dtemplates_FEB2/templates2d_Xmass"+mass+".root","READ");
		h2d = (TH2D*)f2d->Get("h2"+model+"/Z"+modeltmp+"_mXp"+mass)->Clone();
                oldDir->cd();

		//------------------------------------------------------------
		const Int_t nbins = h2d->GetNbinsX();
		Double_t bins[nbins+1];
		TAxis* xaxis = (TAxis*)h2d->GetXaxis();
		for(int i=0 ; i<nbins ; i++) bins[i] = xaxis->GetBinLowEdge(i+1);
		bins[nbins] = xaxis->GetBinUpEdge(nbins);
		//------------------------------------------------------------
		
		Double_t g4 = h2d->GetYaxis()->GetBinLowEdge(g4bin);
		TString g = (TString)_s(sqrt(sqrt(g4)),2);
	
		TString modelname = (model=="ZP") ? "Z'" : "KK";
	
		h1s  = new TH1D(mass+"_subtracted","g="+g+", m_{"+modelname+"} = "+mass+" GeV;m_{ee} TeV;Events",nbins,bins);
		h1   = new TH1D(mass+"_inamplitude","g="+g+", m_{"+modelname+"} = "+mass+" GeV;m_{ee} TeV;Events",nbins,bins);
		h1dy = new TH1D(mass+"_dy","g="+g+", m_{"+modelname+"} = "+mass+" GeV;m_{ee} TeV;Events",nbins,bins);
		for(Int_t bin=0 ; bin<=nbins+1 ; bin++)
		{
			h1s->SetBinContent(bin, h2ds->GetBinContent(bin,g4bin));
			h1->SetBinContent(bin, h2d->GetBinContent(bin,g4bin));
			h1dy->SetBinContent(bin, h2d->GetBinContent(bin,1));
		}

		TLegend* leg = new TLegend(0.55,0.72,0.85,0.85,NULL,"brNDC");
		leg->SetFillStyle(4000); //will be transparent
		leg->SetFillColor(0);
		leg->SetTextFont(42);
		
		TCanvas* c = new TCanvas("c"+mass,"c"+mass,600,400);
		c->Draw();
		c->cd();
		c->SetTicks(1,1);
		c->SetLogy();
		c->SetLogx();
		c->SetGridy();	
		c->cd();

		double ymin = getYmin(h1s);
		ymin = (getYmin(h1dy)<ymin) ? getYmin(h1dy): ymin;
		h1dy->SetMinimum( ymin*0.1);
		TLine* l = new TLine(0.12805,h1dy->GetMinimum(),0.12805,h1dy->GetMaximum());
		l->SetLineColor(kGreen+2);
		
		h1dy->SetLineColor(kBlue);
		h1dy->SetLineStyle(1);
		h1dy->SetLineWidth(2);
		leg->AddEntry(h1dy,"DY, in-amplitude k_{F}","l");
		h1dy->DrawCopy();

		h1s->SetLineColor(kBlack);
		h1s->SetLineStyle(1);
		h1s->SetLineWidth(2);
		leg->AddEntry(h1s,modelname+", subtracted k_{F}","l");
		h1s->DrawCopy("SAMES");

		h1->SetLineColor(kRed);
		h1->SetLineStyle(3);
		h1->SetLineWidth(2);
		leg->AddEntry(h1,modelname+", in-amplitude k_{F}","l");
//.........这里部分代码省略.........

TauAnalysisSelector::TauAnalysisSelector(const std::string& puWeight, bool isEmbedded,
                                         const std::string& embeddingWTauMuFile, const std::string& embeddingWTauMuPath,
                                         const std::string& mcTauMultiplicity, const std::string& mcTauMode,
                                         const std::string& embeddingNormalizationMode,
                                         TH1 *embeddingMuonWeights):
  BaseSelector(),
  //fMuons("Emb"),
  fGenTaus(isEmbedded ? "gentausOriginal" : "gentaus", true),
  fGenTausEmbedded("gentausEmbedded", true),
  fPuWeightName(puWeight),
  fIsEmbedded(isEmbedded),
  fEmbeddingWTauMuWeights(0),
  fMCTauMultiplicity(parseMCTauMultiplicity(mcTauMultiplicity)),
  fMCTauMode(parseMCTauMode(mcTauMode)),
  fEmbeddingNormalizationMode(parseEmbeddingNormalizationMode(embeddingNormalizationMode)),
  fEmbeddingMuonWeights(embeddingMuonWeights),
  cAll(fEventCounter.addCounter("All events")),
  //cElectronVeto(fEventCounter.addCounter("Electron veto")),
  cGeneratorWeight(fEventCounter.addCounter("Vis. pt weight")),
  cTauBRWeight(fEventCounter.addCounter("Tau BR weighting")),
  cGenTauFound(fEventCounter.addCounter("Gen tau found")),
  cTauMCSelection(fEventCounter.addCounter("Tau MC requirement")),
  cOnlyWMu(fEventCounter.addCounter("Only W->mu")),
  cWTauMuWeight(fEventCounter.addCounter("W->tau->mu weighting")),
  cTriggerEffWeight(fEventCounter.addCounter("Muon trigger eff weighting")),
  cIdEffWeight(fEventCounter.addCounter("Muon ID eff weighting")),
  cMuonWeight(fEventCounter.addCounter("Muon corr weighting (from argument)")),
  cJetSelection(fEventCounter.addCounter("Jet selection")),
  cPrimaryVertex(fEventCounter.addCounter("Primary vertex")),
  cAllTauCandidates(fEventCounter.addCounter(">= 1 tau candidate")),
  cPrePtCut(fEventCounter.addCounter("Pre Pt cut")),
  cDecayModeFinding(fEventCounter.addCounter("Decay mode finding")),
  cEtaCut(fEventCounter.addCounter("Eta cut")),
  cPtCut(fEventCounter.addCounter("Pt cut")),
  cLeadingTrackPtCut(fEventCounter.addCounter("Leading track pt")),
  cEcalCracks(fEventCounter.addCounter("ECAL fiducial: cracks")),
  cEcalGap(fEventCounter.addCounter("ECAL fiducial: gap")),
  cAgainstElectron(fEventCounter.addCounter("Against electron")),
  cAgainstMuon(fEventCounter.addCounter("Against muon")),
  cIsolation(fEventCounter.addCounter("Isolation")),
  cOneProng(fEventCounter.addCounter("One prong")),
  cRtau(fEventCounter.addCounter("Rtau")),
  cMuTrigger(fEventCounter.addCounter("Mu trigger"))
{
  if(isEmbedded && !embeddingWTauMuFile.empty()) {
    TFile *file = TFile::Open(embeddingWTauMuFile.c_str());
    TEfficiency *eff = dynamic_cast<TEfficiency *>(file->Get(embeddingWTauMuPath.c_str()));
    TAxis *xaxis = eff->GetPassedHistogram()->GetXaxis();
    fEmbeddingWTauMuWeights = new TH1F("weights", "weights", xaxis->GetNbins(), xaxis->GetBinLowEdge(xaxis->GetFirst()), xaxis->GetBinUpEdge(xaxis->GetLast()));
    fEmbeddingWTauMuWeights->SetDirectory(0);
    for(int bin=1; bin <= xaxis->GetNbins(); ++bin) {
      fEmbeddingWTauMuWeights->SetBinContent(bin, eff->GetEfficiency(bin));
      std::cout << "Bin " << bin << " low edge " << fEmbeddingWTauMuWeights->GetXaxis()->GetBinLowEdge(bin) << " value " << fEmbeddingWTauMuWeights->GetBinContent(bin) << std::endl;
    }
    fEmbeddingWTauMuWeights->SetBinContent(xaxis->GetNbins()+1, fEmbeddingWTauMuWeights->GetBinContent(xaxis->GetNbins()));
    file->Close();
  }

  if(isEmbedded) {
    fMuons.setIdEfficiencyName("efficiency_Run2011AB");
    fMuons.setTriggerEfficiencyName("efficiency_trigger");
  }
}

void make( string charge = "p", bool lines = false, int iCen = 0, string cs = "Pi", string hFile = "histograms.root" ){

	Bichsel bgen;
	gStyle->SetPadLeftMargin(0.16);
	gStyle->SetPadBottomMargin(0.12);
	gStyle->SetPadRightMargin(0.12);

	set_plot_style();

	RooPlotLib rpl;
	string rpName = "rp_" + charge + "_dedx_beta_lines.pdf";
	if ( false == lines )
		rpName = "rp_" + charge + "_dedx_beta.pdf";

	string chargeName = "Positive Tracks : ";
	if ( "n" == charge )
		chargeName = "Negative Tracks : ";

	Reporter rp( rpName, 600, 800 );

	TFile * f = new TFile( hFile.c_str(), "READ" );

	TH1D * nlBeta = (TH1D*)f->Get( "nlBeta" );


	for ( int iPt = 0; iPt < 100; iPt ++ ){
		if ( iPt != 16 ) continue;


		TAxis * x = nlBeta->GetXaxis();
		double lpT = x->GetBinLowEdge( iPt + 1 );
		double hpT = x->GetBinLowEdge( iPt + 1 ) + x->GetBinWidth( iPt + 1 );
		double avgP = (lpT + hpT) / 2.0;

		string name = "dedx_tof/dedx_tof_" + cs + "_" + charge + "_" + ts(iCen) + "_" + ts( iPt);
		TH2D * h2 = (TH2D*)f->Get( name.c_str() );

		if ( !h2 )
			break;
		if ( iPt == 5 ) // what went wrong?
			continue;


		double x1 = h2->GetXaxis()->GetBinLowEdge( 1 );
		double x2 = h2->GetXaxis()->GetBinLowEdge( h2->GetNbinsX() ) + h2->GetXaxis()->GetBinWidth( h2->GetNbinsX() );
		double y1 = h2->GetYaxis()->GetBinLowEdge( 1 );
		double y2 = h2->GetYaxis()->GetBinLowEdge( h2->GetNbinsY() ) + h2->GetYaxis()->GetBinWidth( h2->GetNbinsY() );

		  

		rpl.style( h2 ).set( "draw", "col" )
			.set( "optstat", 0 ).set( "logz", 1 )
			.set( "title", " ; ln(dE/dx) - ln(dE/dx)_{#pi}  ; #beta^{-1} - #beta^{-1}_{#pi} " )
			.set( "ts", 0.16 )
			.set( "xts", 0.06 )
			.set( "xto", 0.75 )
			.set( "yts", 0.07 )
			.set( "yto", 0.85 )
			.set( "yoffset", 0.01 ).draw(); 



		if ( 16 == iPt ){
			gStyle->SetTitleFontSize( 0.45 );
			rpl.style(h2)
				.set( "xr", -0.3, 0.95 )
				.set( "yr", -0.15, 0.6 ).draw();
		}
		TLatex *text = new TLatex( -0.5, 0.64, (chargeName + dts(lpT) + " < p_{T} [GeV/c] < " + dts(hpT)).c_str() );
		text->SetTextSize(0.055);
		text->Draw("");


		if ( lines  ){
			double bKaon = one_beta( mK, avgP ) - one_beta( mPi, avgP );
			double bProton = one_beta( mP, avgP ) - one_beta( mPi, avgP );

			double dKaon = bgen.meanLog( avgP, mK, -1, 1000 ) - bgen.meanLog( avgP, mPi, -1, 1000 );
			double dProton = bgen.meanLog( avgP, mP, -1, 1000 ) - bgen.meanLog( avgP, mPi, -1, 1000 );

			drawTofLines( 0, kRed+1, x1, x2 );		
			drawTpcLines( 0, kRed + 1, y1, y2);
			// around Kaons
			drawTofLines( bKaon, kOrange+1, x1, x2 );	
			drawTpcLines( dKaon, kOrange + 1, y1, y2);
			// around Protons
			drawTofLines( bProton, kBlack, x1, x2 );	
			drawTpcLines( dProton, kBlack, y1, y2);
		}

		// double dElec = bgen.meanLog( avgP, mE, -1, 1000 ) - bgen.meanLog( avgP, mPi, -1, 1000 );
		// double bElec = one_beta( mE, avgP ) - one_beta( mPi, avgP );
		// //drawTpcLines( dMerged, kBlack, y1, y2);
		// TEllipse * ell = new TEllipse( dElec,bElec, sigTpc * 3, sigTof * 3 );
		// ell->SetFillColorAlpha( kBlack, 0 );
		// ell->Draw();
		// 
		
		gPad->SetRightMargin( 0.01 );
		gPad->SetBottomMargin( 0.10 );
//.........这里部分代码省略.........