C++ TH1D::Add方法代码示例

TH1D* plot (std::string num, std::string den,std::vector<TFile*>& tfiles,std::vector<double>& weights) {

    double weightTot = 0.0;

    for ( int i = 0; i<weights.size(); i++) {
        weightTot+=weights[i];
    }

    std::string name = num;
    name+="NEW";
    TH1D* hNum = (TH1D*)tfiles[0]->Get(num.c_str());
    TH1D* HNum = (TH1D*)hNum->Clone(name.c_str());
    HNum->Sumw2();



    name=den;
    name+="NEW";
    TH1D* hDen = (TH1D*)tfiles[0]->Get(den.c_str());
    TH1D* HDen = (TH1D*)hDen->Clone(name.c_str());
    HDen->Sumw2();

    for (int i=1; i<tfiles.size(); i++) {
        TH1D* htempNum = (TH1D*)tfiles[i]->Get(num.c_str());
        TH1D* htempDen = (TH1D*)tfiles[i]->Get(den.c_str());
        HNum->Add(htempNum,weights[i]/weightTot);
        HDen->Add(htempDen,weights[i]/weightTot);
    }
    name = num;
    name+="Div";
    TH1D* hDiv = (TH1D*)HNum->Clone(name.c_str());
    //hDiv->Sumw2();
    hDiv->Divide(HDen);
    return hDiv;
}

TH1D* CutFlow::allMChisto(AllSamples samples, Variable variable){

	TH1D *allMC = (TH1D*)samples.ttbar->histo->Clone("all mc");

	allMC->Add(samples.qcd->histo);
	allMC->Add(samples.vjets->histo);
	allMC->Add(samples.single_t->histo);

	return allMC;
}

// Merge same modes
int dmMerge(TObject *a1, TObject *b1,TList *keys)
{
	TDecayMode *a=(TDecayMode*)a1;
	TDecayMode *b=(TDecayMode*)b1;
	TIter  nexthist(keys);
	TKey  *key_hist=0;
	while(key_hist=(TKey*)nexthist())
	{
		int cycle=-10;
		if(strcmp(key_hist->GetClassName(),"TH1D")==0)
		{
			TH1D  *h=0;
			int cycleh=key_hist->GetCycle();
			if(cycleh<cycle) { cout<<"Skipping..."<<endl; continue; }
			if(cycle<cycleh) cycle=cycleh;
			h=(TH1D*)key_hist->ReadObj();
			if(!h)
			{
				cout<<"Cannot read: "<<key_hist->GetName()<<endl;
				exit(-2);
			}
			TH1D *eh = (TH1D*) a->histograms->FindObject(h->GetName());
			if(!eh) continue;
			if(eh->GetNbinsX()!=h->GetNbinsX() ||
			   eh->GetXaxis()->GetXmax()!=h->GetXaxis()->GetXmax())
			   return -1;
			eh->Add(h);
		}
	}
        a->SetNEntries(a->GetNEntries()+b->GetNEntries());
        a->SetSumw(a->GetSumw()+b->GetSumw());
        a->SetSumw2(a->GetSumw2()+b->GetSumw2());
	return 0;
}

TH1D *CombineCFs(vector<TH1D*> cfs, vector<Double_t> counts)
{
  // Take in a vector of cfs and a vector of counts, and use them to
  // make an averaged correlation function
  UInt_t nCFs = cfs.size();

  if(nCFs == 0) {
    cout << "No CFs found in collection. Cannot combine" << endl;
    return NULL;
  } else if (nCFs == 1) {
    cout<<"Only one cf found in collection."<<endl
	<<"Cannot combine."<<endl
	<<"Returning the found CF."<<endl;;
    return cfs[0];
  }
  assert(nCFs == counts.size());

  TH1D *combinedCF = (TH1D*)cfs[0]->Clone();
  combinedCF->Scale(counts[0]);
  Double_t totalCounts = counts[0];

  // Add together weighted cfs
  for(UInt_t i = 1; i < nCFs; i++) {
    TH1D *copyCF = (TH1D*) cfs[i]->Clone();
    Double_t thisCount = counts[i];
    copyCF->Scale(thisCount);
    combinedCF->Add(copyCF);
    totalCounts += thisCount;
    delete copyCF; copyCF = NULL;
  }
  
  // Now get average
  combinedCF->Scale(1./totalCounts);
  return combinedCF;
}

void calculateBayesErrors(int BgID){
	TH1D *hist;
	for(int k=0; k<2; k++){
		//###sum over A B C D regions
		for(int i=0; i<2; i++){
			for(int j=0; j<2; j++){
				if(i==0 && j==0) hist = (TH1D*)histRaw[BgID-1][k][0][i][j]->Clone();
				else hist->Add(Form("Hist.hist%d%d",i,j),1);
			}
		}
		double scale_factor = L*pb2fb*X_bg[BgID-1];
		//###calculate Bayesian Errors
		for(int i=0; i<2; i++){
			for(int j=0; j<2; j++){
				for(int s=0; s<5; s++){//spectrum
					graph[BgID-1][k][s][i][j] = new TGraphAsymmErrors(Nbins);
					graph[BgID-1][k][s][i][j] -> BayesDivide(histRaw[BgID-1][k][s][i][j],hist);
					for(int bin=0; bin<graph[BgID-1][k][s][i][j]->GetN(); bin++){
						double x,y,err;
						graph[BgID-1][k][s][i][j]->GetPoint(bin,x,y);
						histYield[BgID-1][k][s][i][j]->SetBinContent(bin,y*scale_factor);
						err = sqrt( pow(Err_X_bg[BgID-1]/X_bgi[BgID-1],2) + pow(graph[BgID-1][k][s][i][j]->GetErrorYlow(bin)/y,2) );
						YieldErrors[BgID-1][k][s][i][j][bin][0] = err*(y*scale_factor);
						err = sqrt( pow(Err_X_bg[BgID-1]/X_bgi[BgID-1],2) + pow(graph[BgID-1][k][s][i][j]->GetErrorYhigh(bin)/y,2) );
						YieldErrors[BgID-1][k][s][i][j][bin][1] = err*(y*scale_factor);
					}
				}
			}
		}
	}//end of k
}

// Convert plots to paper format - no title, bigger fonts etc
void combineHists( TFile* fSig, TFile* fBg, TFile* fBg2, std::string histName, std::string plotOpt, std::string outputName, std::vector<double> scalingFactors, std::string label, std::string yTitle="DEFAULT"){
    TH1::SetDefaultSumw2();

    // for 3 hists - sig and 2 bg
    TCanvas c1;
    TH1D* hSig = fSig->Get(histName.c_str());
    hSig->SetLineColor(kRed);
    hSig->SetMarkerSize(0);
    doSignalHist(hSig);
    
    // Make combined BG hist
    // Need to rescale carefully
    std::vector<TFile*> files;
    files.push_back(fBg);
    files.push_back(fBg2);
    // TH1D* hBg = combine(files, histName, scalingFactors);
    TH1D* hBgA = fBg->Get(histName.c_str());
    TH1D* hBgB =  fBg2->Get(histName.c_str());
    TH1D* hBg = (TH1D*) hBgA->Clone();
    double total = scalingFactors[0]+scalingFactors[1];
    hBg->Scale(scalingFactors[0]/total);
    hBg->Add(hBgB, scalingFactors[1]/total);
    hBg->SetMarkerSize(0);
    doAltBGHist(hBg);
    
    THStack st("h","");
    st.Add(hSig);
    st.Add(hBg);
    st.Draw((plotOpt+"NOSTACK").c_str());
    st.GetXaxis()->SetTitle(hSig->GetXaxis()->GetTitle());
    if (yTitle == "DEFAULT") {
        st.GetYaxis()->SetTitle(hSig->GetYaxis()->GetTitle());
    } else {
        st.GetYaxis()->SetTitle(yTitle.c_str());
    }
    setAltTitleLabelSizes(&st.GetHistogram());
    st.SetTitle("");
    st.Draw((plotOpt+"NOSTACK").c_str());

    TLegend* l_all = new TLegend(0.65,0.6,0.89,0.89);
    l_all->AddEntry(hBg,"Gen. level QCD MC","lp");
    l_all->AddEntry((TObject*)0,"(b#bar{b} + q-g scatter,",""); //null pointers for blank entries
    l_all->AddEntry((TObject*)0,"q = b, #bar{b}, c, #bar{c})","");
    l_all->AddEntry(hSig, "Signal MC", "lp");
    l_all->AddEntry((TObject*)0,"m_{#phi} = 8 GeV", "");
    doStandardLegend(l_all);
    l_all->Draw();

    TPaveText t(0.15, 0.75, 0.5, 0.85, "NDC");
    t.AddText(label.c_str());
    doStandardText(&t);
    if (label != "") {
        t.Draw();
    }
    c1.SaveAs(outputName.c_str());

    if (!hSig) delete hSig;
    if (!hBg) delete hBg;
}

	// 1 file argument -> add e+m from same files e.g. mc
	// this is considered to be MC -> added to sum histo
	TH1D* get(const TString& File){ //
		TFile file(path+File);
		TH1D* mH = (TH1D*)file.Get(mName);
		TH1D* eH = (TH1D*)file.Get(eName);
		gROOT->cd();
		TH1D*  H = mH->Clone();// Sumw2() is copied as well
		H->Add(eH);
		return H;
	}

int main(){
	TFile * data = TFile::Open("TreesMu_Data_plots.root");
	TH1D * dataMtop = (TH1D*)data->Get("antiEtaFwDTrue_allW/antiEtaFwDTrue_allWcosTheta");
	
	TFile * tt = TFile::Open("TreesMu_TTBar_RW.root");
	
	TH1D * ttMtop = (TH1D*)tt->Get("MtopOutWindowTrue_allW/MtopOutWindowTrue_allWcosTheta");
	
	TH2D * ttMtop2D = (TH2D*)tt->Get("MtopOutWindowTrue_allW/MtopOutWindowTrue_allWcosTheta2D");
	
  std::pair<TF1, WeightFunctionCreator*> WeightFuncUD(WeightFunctionCreator::getWeightFunction("WeightFuncUDF", F0, FL));
  std::pair<TF1, WeightFunctionCreator*> WeightFuncDU(WeightFunctionCreator::getWeightFunction("WeightFuncDUF", F0, FL));
  cout<<F0 + F0Sys<<"\t"<<FL - FLSys<<endl;
	WeightFuncUD.first.SetParameters(F0 + F0Sys, FL - FLSys);
	WeightFuncDU.first.SetParameters(F0 - F0Sys, FL + FLSys);
	
	TH1D * ttUD = myReweightor(ttMtop2D,WeightFuncUD,"WeightFuncUD_");
	cout<<"ratio: "<<ttMtop2D->Integral()/ttUD->Integral()<<endl;
	ttUD->Scale(ttMtop2D->Integral()/ttUD->Integral());
	ttUD->Add(ttMtop);
	TH1D * wUD = (TH1D*)ttUD->Clone("wUD");
	wUD->Scale(-1.);
	wUD->Add(dataMtop);
	
	TH1D * ttDU = myReweightor(ttMtop2D,WeightFuncDU,"WeightFuncDU_");
	cout<<"ratio: "<<ttMtop2D->Integral()/ttDU->Integral()<<endl;
	ttDU->Scale(ttMtop2D->Integral()/ttDU->Integral());
	ttDU->Add(ttMtop);
	TH1D * wDU = (TH1D*)ttDU->Clone("wDU");
	wDU->Scale(-1.);
	wDU->Add(dataMtop);
	
	TFile * out = new TFile("file.root","recreate");
	out->cd();
	WeightFuncUD.first.Write();
	WeightFuncDU.first.Write();
	ttUD->Write();
	ttDU->Write();	
	wUD->Write();
	wDU->Write();	
	out->Close();
	return 1;
}

TH1D* ratioCalculator(TH1D* passTH1, TH1D* failTH1)
{
  passTH1->Sumw2();
  TH1D *sum = (TH1D*)passTH1->Clone();
  failTH1->Sumw2();
	
  sum->Add(failTH1);
  passTH1->Divide(passTH1,sum,1,1,"B");
  return passTH1;
}

/*
  Returns full reconstructed thetaC
  histogram where each MCP has been 
  shifted such that the mean of the
  gaussian aligns with the expected
  value
*/
TH1D* ThetaCorr( TTree *&tree,
				 bool prot = 1,
				 TString pidcut = "PID>1000",
				 double *&shifts,
				 TString corrtitle = "theta",
				 int bins = 120 )
{
	double angleP  = 0.8168; // assume 7 GeV for now
	double anglePi = 0.8249;
	double angle;
	if(prot) angle = angleP;
	else     angle = anglePi;
	
	const int nMCP = 15;
	TH1D *mcpHist[nMCP];
	TH1D *thetaCorr = new TH1D(corrtitle,corrtitle,bins,0.6,1);
	TF1 *mcpfit = new TF1("mcpfit","gaus");
	mcpfit->SetParameters(100,angle,0.007);

	double diffpeak = DiffPeak(tree,pidcut);
	pidcut += Form(" && abs(diff-%f)<1",diffpeak);

	// loop over mcps
	for(int mcpid = 0; mcpid < nMCP; mcpid++)
	{
        // get timing peak for MCP=mcpid
		TString cut = pidcut + Form(" && mcp==%d",mcpid);
		//cut += Form(" && mcp==%d",mcpid);
		//double diffpeak = DiffPeak(tree,cut);
		//cout << "MCP " << mcpid << " time shift\t" << diffpeak << endl;

		// project from tree
		// using PID cut and time cut
		//cut += Form(" && abs(diff-%f)<1",diffpeak);
		TString mcpname = Form("mcp%d",mcpid);
		mcpHist[mcpid] = new TH1D(mcpname,mcpname,bins,0.6,1);
		tree->Project(mcpname,"theta",cut);

		mcpHist[mcpid]->GetXaxis()->SetRangeUser(angle-0.04,angle+0.04);
		double max = mcpHist[mcpid]->GetXaxis()->GetBinCenter(mcpHist[mcpid]->GetMaximumBin());
		//cout << "max\t" << max << endl;
		mcpHist[mcpid]->Fit(mcpfit,"lq","",max-0.03,max+0.03);
		mcpHist[mcpid]->GetXaxis()->UnZoom();
		
		double shift =  angle - mcpfit->GetParameter(1);
		shifts[mcpid] = shift;
		//cout << "shift\t" << shift << endl;

		tree->Project(mcpname,Form("theta+%f",shift),cut);
		thetaCorr->Add(mcpHist[mcpid]);

	}

	return thetaCorr;
}

TH1D* addScale(std::vector<TH1D*> plots, std::vector<double> scalingFactors) {
    TH1::SetDefaultSumw2();

    TH1D* h = (TH1D*)plots[0]->Clone(plots[0]->GetName());
    h->Scale(scalingFactors[0]);
    for (unsigned i = 1; i < plots.size(); i++) {
        // TH1D* hTmp = (TH1D*)plots[i]->Rebin(nBinsX, plots[0]->GetTitle(), &massBins[0]);
        h->Add(plots[i], scalingFactors[i]);
    }
    return h;
}

TH1D* combine(std::vector<TFile*> files, std::string histName, std::vector<double> scalingFactors) {
    TH1::SetDefaultSumw2();
    // Get 1st hist in list
    TH1D* h = (TH1D*)files[0]->Get(histName.c_str())->Clone(files[0]->Get(histName.c_str())->GetName());
    h->Scale(scalingFactors[0]);
    for (unsigned i = 1; i < files.size(); i++) {
        TH1D* hTmp = (TH1D*) (files[i]->Get(histName.c_str()));
        h->Add(hTmp, scalingFactors[i]);
    }
    return h;
}

/* Get a fake estimate from the given histogram, i.e. return Data - MC.
 *
 * Subtract all background MCs from data but the one specified with
 * "notremove". Under- and overflows are ignored in computation.
 *
 * @param hname     Histogram name from which to compute estimate
 * @return          Histogram of fakes
 */
TH1D * get_fakes_1d(const char * hname)
{
  // get number of single fakes from data histogram
  plot(hname);
  legend();
  TH1D * hData = dataHisto();
  if (hData == 0) {
    THROW("get_fakes() needs a data histogram");
  }
  // data
  double N = hData->Integral();
  INFO("Data events: " << N);

  TH1D * hBack = 0;
  for (unsigned int i = 0; i < sizeof(removeNames)/sizeof(void *); i++) {
    TH1D * hSub = backgroundHisto(removeNames[i]);
    if (hSub == 0) {
      THROW(string("get_fakes() problem getting background ")+removeNames[i]);
    }
    // add backgrounds together
    if (hBack == 0) {
      hBack = hSub;
    }
    else {
      hBack->Add(hSub);
      delete hSub;
    }
  }
  N = hBack->Integral();
  INFO("Background events : " << N);

  // subtract
  hData->Add(hBack, -1.);
  // temporarily needed for function call
  delete hBack;
  return hData;
}

double getQCD(TString Variable, TString Obj, int rebinFact, double *errorRef) {
    TString dir = "rootFilesV4/central/";
    TFile* file = new TFile(dir + "SingleMu_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root");

    TFile* tt_file = new TFile(dir + "TTJet_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root");
    TFile* vjets_file = new TFile(dir + "VJets_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root");
    TFile* stop_file = new TFile(dir + "SingleTop_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root");

    cout << "TTbar_plus_X_analysis/MuPlusJets/QCD non iso mu+jets ge3j/"+Obj+"0btag" << endl;

    TString IsoFolder;
    if(Variable == "HT") {
        IsoFolder = "QCD mu+jets PFRelIso";
    } else {
        IsoFolder = "QCD non iso mu+jets ge3j";
    }


    TH1D* plot = (TH1D*) file->Get("TTbar_plus_X_analysis/MuPlusJets/"+IsoFolder+"/"+Obj+"0btag");
    plot->Sumw2();

    TH1D* tt_plot = (TH1D*) tt_file->Get("TTbar_plus_X_analysis/MuPlusJets/"+IsoFolder+"/"+Obj+"0btag");
    TH1D* vjets_plot = (TH1D*) vjets_file->Get("TTbar_plus_X_analysis/MuPlusJets/"+IsoFolder+"/"+Obj+"0btag");
    TH1D* stop_plot = (TH1D*) stop_file->Get("TTbar_plus_X_analysis/MuPlusJets/"+IsoFolder+"/"+Obj+"0btag");

    TH1D* allMC = (TH1D*)tt_plot->Clone("allMC");

    allMC->Add(vjets_plot);
    allMC->Add(stop_plot);

    plot->Add(allMC, -1);

    *errorRef = sqrt(plot->Integral() + pow(0.5*allMC->Integral(),2));

    return plot->Integral();

}

TH1D* combineRebin10bins(std::vector<TFile*> files, std::string histName, std::vector<double> scalingFactors) {
    TH1::SetDefaultSumw2();
    std::vector<double> massBins = generate10Bins();
    int nBinsX = massBins.size()-1;
    // Get 1st hist in list
    TH1D* h = (TH1D*)files[0]->Get(histName.c_str())->Clone(files[0]->Get(histName.c_str())->GetName());
    h->Scale(scalingFactors[0]);
    
    for (unsigned i = 1; i < files.size(); i++) {
        TH1D* hTmp = (TH1D*) (files[i]->Get(histName.c_str()));
        h->Add(hTmp, scalingFactors[i]);
    }
    // /files[0]->Get(histName.c_str())->GetTitle()
    TH1D* hNew = h->Rebin(nBinsX, files[0]->Get(histName.c_str())->GetTitle(), &massBins[0]);
    return hNew;
}