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

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

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

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

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

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;

}

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

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

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

TH1F *QCD_MatrixMethod_qcdhist(TH1F *looseH, TH1F *tightH, TH1F*epsqcd_hist, Double_t EpsSig, Double_t EpsSigErr,Bool_t makeGlobal, Double_t *NQCD_glob,  Double_t *NQCD_glob_err, Bool_t witherror)
{
  printf("we use the histo method for qcd determination!\n");
  Double_t EpsQCD=0.0;
  Double_t EpsQCDErr=0.0;
  Int_t binhelp=0;
  TH1F *QCDH(0);
  char NewName[500], NewTitle[500], NewHist[1000];
  
  sprintf(NewName,"%s_QCD",looseH->GetName());

  TAxis *xaxis = looseH ->GetXaxis();
  TAxis *yaxis = looseH->GetYaxis();
  
  Int_t Nbins=looseH->GetNbinsX();
  Double_t firstbinx= xaxis->GetBinLowEdge(1);
  Double_t lastbin= xaxis->GetBinUpEdge(Nbins);

  const char* xTitle=xaxis->GetTitle();
  const char* yTitle=yaxis->GetTitle();

  sprintf(NewTitle,"%s - QCD-Histogramm",xTitle);
  sprintf(NewHist,"%s;%s;%s",NewTitle, xTitle ,yTitle);
 
// for global: number of QCD-Background alltogether
    Double_t NQCD_all=0.0, NQCD_all_Error=0.0 ;
  
    // The same idea as in QCD_MatrixMethod but this time assume Loose,
    // Tight and EpsilonQCD are histograms with several bins (instead of
    // scalar values). Thus calculate NQCD_all as a sum of corresponding
    // values per bin.
    for(Int_t i=0;i<Nbins;i++)
      {

	if(i<12) //just to have not too crappy values
	  {
	    binhelp=epsqcd_hist->FindBin(looseH->GetBinCenter(i+1));
	    EpsQCD=epsqcd_hist->GetBinContent(binhelp);
	    EpsQCDErr=epsqcd_hist->GetBinError(binhelp);
	  }

	Double_t NQCD_all_help=(EpsSig*looseH->GetBinContent(i+1)- tightH->GetBinContent(i+1))/(EpsSig -  EpsQCD);
	//	if(NQCD_all_help>0.0)
	  NQCD_all+=EpsQCD*NQCD_all_help;
	Double_t N1=looseH->GetBinContent(i+1) - tightH->GetBinContent(i+1);
	//          Double_t N1Error=TMath::Sqrt(N1);
	Double_t N1Error=TMath::Sqrt(pow(looseH->GetBinError(i+1),2.0) - pow(tightH->GetBinError(i+1),2.0));
	Double_t N2=tightH->GetBinContent(i+1);
	//          Double_t N2Error=TMath::Sqrt(N2);
	Double_t N2Error=tightH->GetBinError(i+1);
	Double_t NtQEpsSig=EpsQCD*(N2 - EpsQCD* (N1 + N2))/pow((EpsSig - EpsQCD),2);
	Double_t NtQEpsQCD=EpsSig*(EpsSig*(N1+N2) - N2 )/pow((EpsSig - EpsQCD),2);
	Double_t NtQN1=EpsQCD*EpsSig/(EpsSig - EpsQCD);
	Double_t NtQN2=EpsQCD*(EpsSig-1)/(EpsSig - EpsQCD);
	Double_t binContentErrorSquare=pow(NtQEpsSig,2)*pow(EpsSigErr,2) + pow(NtQEpsQCD,2)*pow(EpsQCDErr,2) + pow(NtQN1,2)*pow(N1Error,2) + pow(NtQN2,2)* pow(N2Error,2);
	Double_t binContentError=TMath::Sqrt(binContentErrorSquare);
	
	//  	if(NQCD_all_help>0.0)
	  NQCD_all_Error+=binContentErrorSquare;
	
	//   printf("Anzahl QCD gesamt:%f, Fehler dieser Anzahl:%f\n\n",NQCD_gesamt,NQCD_gesamt_Error);
      }
      // End loop over bins

    NQCD_all_Error=TMath::Sqrt(NQCD_all_Error);
    // End for global
    
    if(makeGlobal && NQCD_glob)
      {
	*NQCD_glob=NQCD_all;
	if(NQCD_glob_err)
	  *NQCD_glob_err=NQCD_all_Error;
	return QCDH;
      }

    QCDH=new TH1F(NewName,NewHist,Nbins,firstbinx, lastbin);
    QCDH->Sumw2();


    for(Int_t i=0;i<Nbins;i++)
      {
	if(i<12) //just to have not too crappy values
	  {
	    binhelp=epsqcd_hist->FindBin(looseH->GetBinCenter(i+1));
	    EpsQCD=epsqcd_hist->GetBinContent(binhelp);
	    EpsQCDErr=epsqcd_hist->GetBinError(binhelp);
	  }
        Double_t binContentloose=looseH->GetBinContent(i+1);
        Double_t binContenttight=tightH->GetBinContent(i+1);
        Double_t NQCD=(EpsSig*binContentloose - binContenttight)/(EpsSig -  EpsQCD);
        Double_t binContent=EpsQCD*NQCD;
	//  QCDEffizienz+= binContent;
        if(binContent <0.0)
	  binContent=0.0;
        QCDH->SetBinContent(i+1,binContent);
	// printf("QCDEffizienz:%f\n",QCDEffizienz);
	
	
	
        if(witherror)
//.........这里部分代码省略.........

TH1F *QCD_MatrixMethod(TH1F *looseH, TH1F *tightH, Double_t EpsQCD, Double_t EpsSig,  Double_t EpsQCDErr, Double_t EpsSigErr,Bool_t makeGlobal, Double_t *NQCD_glob,  Double_t *NQCD_glob_err, Bool_t witherror,Double_t *NW_glob, Double_t *NW_glob_err)
{

  TH1F *QCDH(0);
  char NewName[500], NewTitle[500], NewHist[1000];
  
  sprintf(NewName,"%s_QCD",looseH->GetName());

  TAxis *xaxis = looseH ->GetXaxis();
  TAxis *yaxis = looseH->GetYaxis();
  
  Int_t Nbins=looseH->GetNbinsX();
  Double_t firstbinx= xaxis->GetBinLowEdge(1);
  Double_t lastbin= xaxis->GetBinUpEdge(Nbins);
  const char* xTitle=xaxis->GetTitle();
  const char* yTitle=yaxis->GetTitle();
  sprintf(NewTitle,"%s - QCD-Histogramm",xTitle);
  sprintf(NewHist,"%s;%s;%s",NewTitle, xTitle ,yTitle);
 
// for global: number of QCD-Background alltogether
    Double_t NQCD_all=0.0, NQCD_all_Error=0.0 ;
    Double_t NW_all=0.0, NW_all_Error=0.0 ;
    // Functions with one bin only that would count how many loose or tight
    // events are there.
    TH1F *qcd_all_loose=new TH1F(hn(),"total number of QCD events",1,0.0,1.0);
    TH1F *qcd_all_tight=new TH1F(hn(),"total number of QCD events",1,0.0,1.0);
    qcd_all_loose->Sumw2();
    qcd_all_tight->Sumw2();

    // Count events
    for(Int_t i=0;i<=Nbins;i++)
      {
        Double_t binContentloose=looseH->GetBinContent(i+1);
        Double_t binContenttight=tightH->GetBinContent(i+1);
        qcd_all_loose->Fill(0.5,binContentloose);
        qcd_all_tight->Fill(0.5,binContenttight);
      }

    for(Int_t i=0;i<1;i++)
      {
	NQCD_all=(EpsSig*qcd_all_loose->GetBinContent(i+1)- qcd_all_tight->GetBinContent(i+1))/(EpsSig -  EpsQCD);
	NQCD_all=EpsQCD*NQCD_all;
	Double_t N1=qcd_all_loose->GetBinContent(i+1) - qcd_all_tight->GetBinContent(i+1);
	//          Double_t N1Error=TMath::Sqrt(N1);
	//	Double_t N1Error=TMath::Sqrt(pow(qcd_all_loose->GetBinError(i+1),2.0) - pow(qcd_all_tight->GetBinError(i+1),2.0));
	Double_t N1Error=TMath::Sqrt(N1);

	Double_t N2=qcd_all_tight->GetBinContent(i+1);
	//          Double_t N2Error=TMath::Sqrt(N2);
	//	Double_t N2Error=qcd_all_tight->GetBinError(i+1);
	Double_t N2Error=TMath::Sqrt(N2);

	Double_t NtQEpsSig=EpsQCD*(N2 - EpsQCD* (N1 + N2))/pow((EpsSig - EpsQCD),2);
	Double_t NtQEpsQCD=EpsSig*(EpsSig*(N1+N2) - N2 )/pow((EpsSig - EpsQCD),2);

	Double_t NtQN1=EpsQCD*EpsSig/(EpsSig - EpsQCD);
	Double_t NtQN2=EpsQCD*(EpsSig-1)/(EpsSig - EpsQCD);

	Double_t binContentErrorSquare=pow(NtQEpsSig,2)*pow(EpsSigErr,2) + pow(NtQEpsQCD,2)*pow(EpsQCDErr,2) + pow(NtQN1,2)*pow(N1Error,2) + pow(NtQN2,2)* pow(N2Error,2);
	Double_t binContentError=TMath::Sqrt(binContentErrorSquare);
	
	NQCD_all_Error=binContentError;
	//	printf("N1: %f N2: %f\n", N1, N2);
	//	printf("NtQEpsSig: %f NtQEpsQCD:%f %f %f %f\n", NtQEpsSig, NtQEpsQCD, N2 - EpsQCD* (N1 + N2), EpsSig*(N1+N2) - N2 , pow((EpsSig - EpsQCD),2));
	//	printf("the bin error: %f where value is %f\n",  binContentError, NQCD_all);
	//	printf("the single contributions: %f %f %f %f\n", pow(NtQEpsSig,2)*pow(EpsSigErr,2),pow(NtQEpsQCD,2)*pow(EpsQCDErr,2), pow(NtQN1,2)*pow(N1Error,2),  pow(NtQN2,2)* pow(N2Error,2));
	//   printf("Anzahl QCD gesamt:%f, Fehler dieser Anzahl:%f\n\n",NQCD_gesamt,NQCD_gesamt_Error);
      }
    // End for global
    
    if(makeGlobal && NQCD_glob)
      {
	//set the NQCD value to 0 if it's smaller
	if(NQCD_all < 0.0) NQCD_all=0.0;
	*NQCD_glob=NQCD_all;
	if(NQCD_glob_err)
	  *NQCD_glob_err=NQCD_all_Error;

      }

    for(Int_t i=0;i<1;i++)
      {
	NW_all=( qcd_all_tight->GetBinContent(i+1) - EpsQCD * qcd_all_loose->GetBinContent(i+1))/(EpsSig -  EpsQCD);
	NW_all=EpsSig*NW_all;

	Double_t N1=qcd_all_loose->GetBinContent(i+1) - qcd_all_tight->GetBinContent(i+1);
	//          Double_t N1Error=TMath::Sqrt(N1);
	//	Double_t N1Error=TMath::Sqrt(pow(qcd_all_loose->GetBinError(i+1),2.0) - pow(qcd_all_tight->GetBinError(i+1),2.0));
	Double_t N1Error=TMath::Sqrt(N1);

	Double_t N2=qcd_all_tight->GetBinContent(i+1);
	//          Double_t N2Error=TMath::Sqrt(N2);
	//	Double_t N2Error=qcd_all_tight->GetBinError(i+1);
	Double_t N2Error=TMath::Sqrt(N2);

	Double_t NtWEpsSig=EpsQCD*(N2 - EpsQCD* (N1 + N2))/pow((EpsSig - EpsQCD),2);
	Double_t NtWEpsQCD=EpsSig*(EpsSig*(N1+N2) - N2 )/pow((EpsSig - EpsQCD),2);

	Double_t NtWN1=EpsQCD*EpsSig/(EpsSig - EpsQCD);
	Double_t NtWN2=EpsSig*(1-EpsQCD)/(EpsSig - EpsQCD);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
   TCanvas *c1 = new TCanvas("c1","c1",1200,600);
   c1->SetFrameFillColor(kYellow-3);
   c1->SetGridx();
   c1->SetGridy();
   h->Draw();
   //create a TPaveText to show the summary results
   TPaveText *pvt = new TPaveText(.5,.9,.75,.99,"brNDC");
   pvt->Draw();
   //create a TPaveLabel to show the time
   TPaveLabel *ptime = new TPaveLabel(.905,.7,.995,.76,"time","brNDC");
   ptime->SetFillColor(kYellow-3);
   ptime->Draw();
   //draw producer identifier
   TNamed *named = (TNamed*)T->GetUserInfo()->FindObject("SysInfo");
   TText tmachine;
   tmachine.SetTextSize(0.02);
   tmachine.SetNDC();
   if (named) tmachine.DrawText(0.01,0.01,named->GetTitle());

   //start loop on selected rows
   Int_t bin,nb=0,j;
   Long64_t ipos;
   Double_t dbin,rest,time;
   Double_t updateLast = 0;
   Int_t nleaks = 0;
   Int_t i;
   for (i=0;i<nsel;i++) {
      pos    = V1[i];
      ipos = (Long64_t)(pos-ivmin);
      nbytes = (Int_t)V2[i];
      time = 0.0001*V3[i];
      bin = axis->FindBin(pos);
      if (bin<1 || bin>nbins) continue;
      dbin = axis->GetBinUpEdge(bin)-pos;
      if (nbytes > 0) {
         halloc->Fill(pos);
         if (dbin > nbytes) dbin = nbytes;
         //fill bytes in the first page
         h->AddBinContent(bin,100*dbin/dv);
	 //fill bytes in full following pages
         nb = Int_t((nbytes-dbin)/dv);
	 if (bin+nb >nbins) nb = nbins-bin;
         for (j=1;j<=nb;j++) h->AddBinContent(bin+j,100);
	 //fill the bytes remaining in last page
         rest = nbytes-nb*dv-dbin;
	 if (rest > 0) h->AddBinContent(bin+nb+1,100*rest/dv);
	 //we save nbytes at pos. This info will be used when we free this slot
         if (nbold[ipos] > 0) printf("reallocating %d bytes (was %lld) at %lld, entry=%d\n",nbytes,nbold[ipos],ipos,i);
         if (nbold[ipos] == 0) {
            nleaks++;            
            //save the Tree entry number where we made this allocation
            ientry[ipos] = i;
         }
         nbold[ipos] = nbytes;
      } else {
         hfree->Fill(pos);
         nbytes = nbold[ipos];
	 if (bin+nb >nbins) nb = nbins-bin;
	 nbold[ipos] = 0; nleaks--;
	 if (nbytes <= 0) continue;
         //fill bytes free in the first page
         if (dbin > nbytes) dbin = nbytes;
	 h->AddBinContent(bin,-100*dbin/dv);
	 //fill bytes free in full following pages
	 nb = Int_t((nbytes-dbin)/dv);
	 if (bin+nb >nbins) nb = nbins-bin;

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

void writeSignalHistosForModel(std::vector<TH1D *>& vsd,
			       const TString& sigmodel,
			       TFile *allHistFile)
{
  for (int ichan=0; ichan<NUMCHAN; ichan++) {
    TH1D * sdh = vsd[ichan];

    // Find limit window from gaussian fit to signal peak.
    //
    double wid  = sdh->GetRMS();
    double mean = sdh->GetMean();

    //TCanvas *c1 = new TCanvas(s,s,300,300);
    TFitResultPtr r = sdh->Fit("gaus","QNS","",mean-2.5*wid,mean+2.5*wid);

    cout<<" mean= "<<mean<<", RMS= "<<wid<<", Fit sigma= "<<r->Parameter(2)<<endl;
    //cout<<r->Parameter(0)<<" "<<r->Parameter(1)<<" "<<r->Parameter(2)<<endl;

    TAxis *xax = sdh->GetXaxis();

#if 0
    int lobin = xax->FindFixBin(r->Parameter(1)-2*r->Parameter(2));
    int hibin = xax->FindFixBin(r->Parameter(1)+2*r->Parameter(2));
    sd.sumwinmin = xax->GetBinLowEdge(lobin);
    sd.sumwinmax = xax->GetBinUpEdge(hibin);
#elif 0
    int lobin = xax->FindFixBin(140);   sd.sumwinmin=140; // 1 bin  left,
    int hibin = xax->FindFixBin(170)-1; sd.sumwinmax=170; // 2 bins right
#else
    int lobin = xax->FindFixBin(sumwinmin);
    int hibin = xax->FindFixBin(sumwinmax)-1;
#endif
    int nbins = hibin-lobin+1;

    // for variable binning - all histos must have the same binning per channel
    TVectorD xbins   = TVectorD(sdh->GetNbinsX(),sdh->GetXaxis()->GetXbins()->GetArray());
    TVectorD xwindow = xbins.GetSub(lobin-1,hibin);

    xax->SetRange(lobin,hibin);

    // Copy contents to window-restricted signal histogram
    // and write to output file.
    //
    TString name = Form("Signal%s_%s",sigmodel.Data(),channames[ichan]);
    printf("Booking TH1D(%s,%s,%d,xwindowarray)\n",
	   name.Data(),sdh->GetTitle(),nbins);
    TH1D *signm = new TH1D(name.Data(),
			   sdh->GetTitle(),
			   nbins,
			   xwindow.GetMatrixArray());

    // make copies of the histograms that are restricted to the
    // bin range lobin-hibin
    //
    for (int ibin=lobin; ibin<=hibin; ibin++) 
      signm->SetBinContent((ibin-lobin+1),
			   sdh->GetBinContent(ibin)
			   *sdh->GetBinWidth(ibin)
			   );

    if (!sigmodel.CompareTo("wh"))
      signm->Scale(whsigscaleto);
    else
      signm->Scale(sigscaleto);

    allHistFile->WriteTObject(signm);

  } // channel loop
}                                                     // writeSignalHistosForModel