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

double Prediction::GetLimitedRatio(TH1D* h1, TH1D* h2, float min, float max, bool binomial, float &relerror){
	if(h1->GetNbinsX() != h2->GetNbinsX()){cout << "GetIntegral ERROR! " << endl; exit(-1);}
	TH1D* h1cp = h1->Clone("h1cp");
	TH1D* h2cp = h2->Clone("h2cp");
	for(int i=1; i<=h1cp->GetNbinsX(); ++i){
		if(h1cp->GetBinLowEdge(i) <min || h1cp->GetBinLowEdge(i)+h1cp->GetBinWidth(i) >max){
			h1cp->SetBinContent(i, 0);
			h1cp->SetBinError(i, 0);
			h2cp->SetBinContent(i, 0);
			h2cp->SetBinError(i, 0);
		}
	}
	TH1D* h1cp_2 = GetScaledHisto(h1cp, 1, 0, h1cp->GetNbinsX());
	TH1D* h2cp_2 = GetScaledHisto(h2cp, 1, 0, h1cp->GetNbinsX());

	TH1D* r      = h1cp_2->Clone("r");
	r->Divide(h1cp_2, h2cp_2, 1, 1, binomial? "B":"");

	*fLogStream << "Getting Limited Ratio for " << h1->GetName()  << " " << h2->GetName() << endl;
	*fLogStream << ">>> ratio = " << r->GetBinContent(1)  << " pm " << r->GetBinError(1) << endl;
	if(fDoRatioStatErr){
		relerror = r->GetBinError(1)/r->GetBinContent(1);
	}
	double  ratio = r->GetBinContent(1);
	delete h1cp;
	delete h2cp;
	delete h1cp_2;
	delete h2cp_2;
	delete r;
	return ratio;
}

void MCefficiency(TString inputmc, TString outputfile, TString tfend, TString selmcgen, TString cut, TString weight, Float_t centmin, Float_t centmax)
{
  Float_t hiBinMin,hiBinMax;
  hiBinMin = centmin*2;
  hiBinMax = centmax*2;

  selmcgen = Form("%s&&hiBin>=%f&&hiBin<%f",selmcgen.Data(),hiBinMin,hiBinMax);
  cut = Form("%s&&hiBin>=%f&&hiBin<%f",cut.Data(),hiBinMin,hiBinMax);

  TFile* infMC = new TFile(inputmc.Data());
  TTree* ntMC = (TTree*)infMC->Get("ntDkpi");
  ntMC->AddFriend("ntHi");
  ntMC->AddFriend("ntSkim");
  TTree* ntGen = (TTree*)infMC->Get("ntGen");
  ntGen->AddFriend("ntHi");
  ntGen->AddFriend("ntSkim");

  TH1D* hPtMC = new TH1D("hPtMC","",nPtBins,ptBins);
  TH1D* hPtGen = new TH1D("hPtGen","",nPtBins,ptBins);

  ntMC->Project("hPtMC","Dpt",TCut(weight)*(TCut(cut.Data())&&"(Dgen==23333)"));
  divideBinWidth(hPtMC);
  ntGen->Project("hPtGen","Gpt",TCut(weight)*(TCut(selmcgen.Data())));
  divideBinWidth(hPtGen);
  TH1D* hEff = (TH1D*)hPtMC->Clone("hEff");
  hEff->Divide(hPtGen);

  TFile* fout=new TFile(Form("%s_cent_%.0f_%.0f_%s.root",outputfile.Data(),centmin,centmax,tfend.Data()),"recreate");
  fout->cd();
  hEff->Write();
  fout->Close();
}

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

void plotter::draw_purity(TH1D* numerator_, TH1D* denominator_, TString file_name){
  TH1D* numerator = (TH1D*) numerator_->Clone("numerator");
  TH1D* denominator = (TH1D*) denominator_->Clone("denominator");

  TH1D* purity = numerator; // just to set correct binning
  purity->Divide(numerator, denominator, 1., 1., "B");

  purity->SetTitle(" ");
  purity->GetXaxis()->SetTitle("m_{gen}");
  purity->GetYaxis()->SetTitle("purity");
  purity->GetYaxis()->SetRangeUser(0,1);

  purity->GetXaxis()->SetTitleSize(0.05);
  purity->GetYaxis()->SetTitleSize(0.05);
  purity->GetXaxis()->SetTitleOffset(0.9);
  purity->GetYaxis()->SetTitleOffset(0.8);
  purity->GetXaxis()->SetNdivisions(505);
  purity->GetYaxis()->SetNdivisions(505);

  purity->SetMarkerStyle(20);
  purity->SetMarkerSize(0.8);
  purity->SetLineColor(1);

  TCanvas *c= new TCanvas("Purity","",600,600);
  gPad->SetLeftMargin(0.15);
  TGaxis::SetMaxDigits(3);
  purity->Draw("E1");
  gPad->RedrawAxis();
  c->SaveAs(directory + file_name + ".pdf");
  delete c;
}

TH1* processFile(TString fname, double scaleFactor){

  static int hcounter=0;

  cout<<"Doing fname="<<fname<<endl;
  TFile *_file2 = TFile::Open(fname);
  PS->cd();
  TString plot =  "lLV[0].lpfIso:lLV[0].fP.Eta()";
  TString hname = Form("h%i",hcounter);
  jets2p->Draw(plot +" >> " + hname+"(24,-3,3,100,0,2)","","qn");
  TH2D *aux = (TH2D*)gDirectory->Get(hname);
  //aux->Scale(scaleFactor);
  TH1D * ss = aux->ProjectionX(hname+"ss",1,10); // signal region 
  TH1D * ni = aux->ProjectionX(hname+"ni",16,100);// nonIso region (use 11
  //cout<<"h="<<hcounter<<" ss="<<ss->GetEntries()<<" ni="<<ni->GetEntries()<<endl;

  TString statis = Form("( %i / %i )",(int)ss->GetEntries(),(int)ni->GetEntries());
  ss->Divide(ni);
  ss->GetYaxis()->SetTitle(TString("N_{signal}/N_{nonIso}")+statis);
  ss->GetXaxis()->SetTitle("#eta");
  
  hcounter++;

  return ss;

}

TH1D* Prediction::GetRatio(TH1D* h1, TH1D* h2, int mergenbins, bool binomial){
	// compute Zll/Photon ratio
	*fLogStream << "------------------------------------------------------------------------" << endl;
	*fLogStream << "Compute ratio between " << h1->GetName() << " and " << h2->GetName()      << endl;           

	
	if(h1==0 || h2==0 ){
		cout << "GetMCZnunuToPhotonRatio: received 0 pointer!" << endl;
		exit(-1);
	}

	// GetScaled histos with only one bin and probagated errors: stat error and error on scale factor
	TH1D *currh1   = GetScaledHisto(h1, 1, 0, mergenbins); 
	TH1D *currh2   = GetScaledHisto(h2, 1, 0, mergenbins); 

	TString name = h1->GetName();
	name         += "_";
	name         += h2->GetName();
	name         += "_Ratio";
	TH1D *ratio = (TH1D*) currh1->Clone(name);
	ratio->Divide(currh1, currh2, 1,1, binomial?"B":""); // binomial errors
	
	ratio->SetLineColor(kBlack);
	ratio->SetMarkerColor(kBlack);
	ratio->SetMarkerStyle(4);

	if(fWriteToFile){
		ratio ->Write();
	}
	delete currh1;
	delete currh2;
	return ratio;
}

void CheckTaggEffButVsH() {
	Char_t Name[256];
	Char_t *VarName;

	//Connect to all existing histograms
	RootFileBut = new TFile(Str_RootFilesBut);
	TH2D *hBut2D = (TH2D*)gROOT->FindObject("hTaggerTime");
	RootFileH = new TFile(Str_RootFilesH);
	TH2D *hH2D = (TH2D*)gROOT->FindObject("hTaggerTime");

	gROOT->cd();

	TH1D *hButBg = hBut2D->ProjectionY("TempBut1D", 800, 950);
	TH1D *hHBg = hH2D->ProjectionY("TempH1D", 800, 950);


	TCanvas *c1 = new TCanvas();
	gStyle->SetPalette(1);
	c1->Divide(2,2);
	c1->cd(1);
	hButBg->Draw("");
	c1->cd(2);
	hHBg ->Draw("");
	c1->cd(3);
	delete gROOT->FindObject("Diff");
	TH1D *hDiff = (TH1D*) hButBg->Clone("Diff");
	hDiff->Divide(hHBg);
	hDiff->Draw();
	

	RootFileBut->Close();
	RootFileH->Close();

	printf("Finished.\n");
}

void merge(){
	TVectorD Nevent;	Nevent.ResizeTo(nbin);  Nevent.Zero();
     //   TVectorD totmultall;	totmultall.ResizeTo(nbin);      totmultall.Zero();
     //   TVectorD avgmultall;	avgmultall.ResizeTo(nbin);      avgmultall.Zero();
        TVectorD tottrk;	tottrk.ResizeTo(nbin);      tottrk.Zero();
//	TVectorD totptall;      totptall.ResizeTo(nbin);    totptall.Zero();
//	TVectorD totetaall;      totetaall.ResizeTo(nbin);    totetaall.Zero();
        TVectorD avgtrk;	avgtrk.ResizeTo(nbin);      avgtrk.Zero();
//	TVectorD avgmult;       avgmult.ResizeTo(nbin);
            
        TH2F* s[nbin];
        TH2F* b[nbin];
        for(int ibin=0;ibin<nbin;ibin++){
            s[ibin] = new TH2F(Form("s_%d",ibin),Form("signal",ibin),detastep,detamin,detamax,dphistep,dphimin,dphimax);
            s[ibin]->Sumw2();
            b[ibin]  = new TH2F(Form("b_%d",ibin), "background",detastep,detamin,detamax,dphistep,dphimin,dphimax);
            b[ibin]->Sumw2();
        }
        TFile *fout = new TFile(Form("Anav3_merged.root"),"Recreate");
        TFile *f[nFileAll];
        for(int ifile=; ifile<63; ifile++){
                f[ifile] = TFile::Open(Form("%s/Anav3_%d.root",outdir.Data(),ifile));
                if(!f[ifile]) continue;
		TVectorD* Nevent_t =  (TVectorD*)f[ifile]->Get(Form("Nevent"));
	//	TVectorD* totmultall_t =  (TVectorD*)f[ifile]->Get(Form("totmultall"));
		TVectorD* tottrk_t =  (TVectorD*)f[ifile]->Get(Form("tottrk"));
	//	TVectorD* totptall_t =  (TVectorD*)f[ifile]->Get(Form("totptall"));
//		TVectorD* totetaall_t =  (TVectorD*)f[ifile]->Get(Form("totetaall"));
		for(int ibin=0;ibin<nbin;ibin++){
			//	totptall[ibin] += (*totptall_t)[ibin];
			//	totetaall[ibin] += (*totetaall_t)[ibin];
			        Nevent[ibin] += (*Nevent_t)[ibin];
			  //      totmultall[ibin] += (*totmultall_t)[ibin];	
			        tottrk[ibin] += (*tottrk_t)[ibin];	
                                TH1D* s_t = (TH1D*)f[ifile]->Get(Form("signal_%d",ibin));
                                TH1D* b_t = (TH1D*)f[ifile]->Get(Form("background_%d",ibin));
                                s[ibin]->Add(s_t);
                                b[ibin]->Add(b_t);
		}
		f[ifile]->Close();
        }
    for(int ibin=0;ibin<nbin;ibin++){
    avgtrk[ibin] = tottrk[ibin]/Nevent[ibin];
    TH1D* hr = (TH1D*)s[ibin]->Clone("hr");
    hr->Divide(b[ibin]);
    hr->Scale(b[ibin]->GetBinContent(b[ibin]->FindBin(0,0)));

    fout->cd();
    Nevent.Write("Nevent");
    avgtrk.Write("avgtrk");
    TDirectory *dir0 = (TDirectory*)fout->mkdir(Form("D_%d",ibin));
    dir0->cd();
    s[ibin]->Write("s");
    b[ibin]->Write("b");
    hr->Write();
    }
}

TH1D* GetRawCorrFunc1D_ratio(int itrg, int jass)
{
                    TH1D* hsignalphi = GetRawSignal1D(itrg,jass);
                    TH1D* hcorrphi = (TH1D*)hsignalphi->Clone(Form("corrphi_trg%d_ass%d",itrg,jass));
		    hcorrphi->SetYTitle("C(#Delta#phi)");  
		    TH1D* hbackphi = GetRawBackground1D(itrg,jass);
		    hcorrphi->Divide(hbackphi); 
                    float max = hcorrphi->GetBinContent(hcorrphi->GetMaximumBin());
//		    hcorrphi->SetAxisRange(ymin,max*1.3,"Y");
                    return hcorrphi;
}

void compareOldandNew()
{
  TCanvas * c1 = new TCanvas("c1","c1",800,600);
  TLatex lat;  

  int centBounds[6] = {60,80,100,120,140,200};
  for(int i = 0 ; i<5; i++){
  const char * histName = "pPbPbp_FF";
  TFile * fold = TFile::Open("FragmentationFunctionsUE2_nominal.root","read");
  TFile * fnew = TFile::Open("FragmentationFunctionsUE2.root","read");
  TH1D * num = (TH1D*)fnew->Get(Form("%s_%d_%d",histName,centBounds[i],centBounds[i+1]));
  TH1D * den = (TH1D*)fold->Get(Form("%s_%d_%d",histName,centBounds[i],centBounds[i+1]));
  //den->Add(num);
  num->Divide(den);
  //num->Scale(2);
  num->GetYaxis()->SetTitle("(no L2 residual)/(Nominal)");
  num->GetYaxis()->SetRangeUser(0.7,1.3);
  num->GetXaxis()->SetTitle("p_{T}");
  num->GetXaxis()->SetRangeUser(0.5,200);
  num->Print("All");
  num->Draw();
  
  lat.DrawLatex(1,1.2,Form("%d < p_{T}^{jet} < %d",centBounds[i],centBounds[i+1]));

  c1->SetLogx();
  c1->SaveAs(Form("diffPlots/noL2Residual_%d_%d.png",centBounds[i],centBounds[i+1]));
  c1->SaveAs(Form("diffPlots/noL2Residual_%d_%d.pdf",centBounds[i],centBounds[i+1]));
  fnew->Close();
  fold->Close();
  }  

  /*const char * histName = "PbPbTrackSpectrum_0_5";
  TFile * fnew = TFile::Open("Spectra_Jun9_noChi2Cut.root","read");
  TFile * fold = TFile::Open("Spectra_Jun9_withChi2Cut.root","read");
  TH1D * num = (TH1D*)fnew->Get(histName);
  TH1D * den = (TH1D*)fold->Get(histName);
  //den->Add(num);
  num->Divide(den);
  //num->Scale(2);
  num->GetYaxis()->SetTitle("(No Chi2)/(With Chi2)");
  num->GetYaxis()->SetRangeUser(0.7,1.5);
  num->GetXaxis()->SetRangeUser(0.7,350);
  num->Print("All");
  num->Draw();
  c1->SetLogx();
  c1->SaveAs("plots/comparisonPlots/Chi2CutTest_PbPb_0_5.png");
  c1->SaveAs("plots/comparisonPlots/Chi2CutTest_PbPb_0_5.pdf");
  */

  //c1->SaveAs("plots/comparisonPlots/ppChargeFraction2.C");
}

TH1D* makeRatioHist(TH1D& Hist1, TH1D& Hist2, char* name="CaloJetPtRatio2over1", double maxX=100.0) {
  
  TH1D* hist = Hist2.Clone(name);
  TH1D* temphist = Hist1.Clone("temphist");
  //temphist->Add( &Hist2, -1);
  hist->Divide(temphist);
  hist->GetXaxis()->SetTitle("jet p_{T} (GeV/c)");
  hist->SetMinimum(0.0);
  if(maxX > 110.0) hist->SetMaximum(0.35);
  else hist->SetMaximum(0.3);
  hist->GetXaxis()->SetRangeUser(0,maxX);
  delete temphist;
  return hist;
}

void SaveClosure(TH1D* prediction, TH1D* expectation, TDirectory* Folder) // prediction durch expectation
{
  TH1D* Closure = (TH1D*) prediction->Clone();
  Closure->Divide(prediction,expectation,1,1,"B");
  TString title = prediction->GetTitle();
  title +="_closure";
  // 	title = "#mu & e Control-Sample Ratio in Search Bins; Search bins; #mu / e CS";
  Closure->SetTitle(title);
  title = prediction->GetName();
  title+="_closure";
  Closure->SetName(title);
  Folder->cd();
  Closure->Write();
}

void mytestreweight(){


  TFile*finput=new TFile("fworking.root");
  TH1D*hBPtFONLL=(TH1D*)finput->Get("hBPt");
  TH1D*hDPtFONLL=(TH1D*)finput->Get("hDPt");
  TH2D*hD=(TH2D*)finput->Get("hD");
  
   hD->SetXTitle("B p_{T} (GeV/c)");
   hD->SetYTitle("D^{0} p_{T} (GeV/c)");
   hBPtFONLL->SetXTitle("B p_{T} (GeV/c)");
   hDPtFONLL->SetXTitle("D p_{T} (GeV/c)");
  
  TH2D *hDreweight = (TH2D*)hD->Clone("hDreweight");
  TH1D *hBPtPythia = hD->ProjectionX("hBPtPythia");
  TH1D *hDPtPythia = hD->ProjectionY("hDPtPythia");

   for (int x=1;x<=hDreweight->GetNbinsX()+1;x++){  //loop over the B pt bins
      if (hBPtPythia->GetBinContent(x)==0) continue;
      double ratio = hBPtFONLL->GetBinContent(x)/hBPtPythia->GetBinContent(x);  // in each pt bin we calculate the ratio of pythiaB/FONLLB 
       
      for (int y=1;y<=hDreweight->GetNbinsY()+1;y++){   //loop over the D pt bins
         double ratio2 = ratio;                                              //copy the ratio for each B bin
         double val = hDreweight->GetBinContent(x,y)*ratio2;
         double valError = hDreweight->GetBinError(x,y)*ratio2;
	 hDreweight->SetBinContent(x,y,val);
	 hDreweight->SetBinError(x,y,valError);
      }   
   }

  TH1D *hBPtPythiaReweight = hDreweight->ProjectionX("hBPtPythiaReweight");
  TH1D *hDPtPythiaReweight = hDreweight->ProjectionY("hDPtPythiaReweight");
  
  hBPtPythiaReweight->Divide(hBPtFONLL);
   
  TCanvas*canvas=new TCanvas("canvas","canvas",1000,500);  
  canvas->SetLogy();
  canvas->Divide(3,1);
  canvas->cd(1);
  hBPtFONLL->Draw();
  canvas->cd(2);
  hBPtPythia->Draw();
  canvas->cd(3);
  hBPtPythiaReweight->Draw();

}

//========================
void PlotRatio1D(Int_t iEff, Int_t iVar){

  Char_t name[100];
  sprintf(name, "c1_%sEff_%s", effName[iEff], varName[iVar]);
  TCanvas *c1 = new TCanvas(name, name);

  sprintf(name, "hRatio_%sEff_%s", effName[iEff], varName[iVar]);
  TH1D *hRatio = (TH1D *) hEffTP[iEff][iVar]->Clone(name);
  hRatio->Divide(hEffMCTruth1D[iEff][iVar]);
  hRatio->SetMarkerStyle(20);
  hRatio->SetMinimum(0.);
  hRatio->SetMaximum(1.5);
  hRatio->Draw("p");

  sprintf(name, "Figures/ratio_%sEff_%s.pdf", effName[iEff], varName[iVar]);
  c1->Print(name);
}

void quickDphi()
{
  TH1::SetDefaultSumw2();
  
  TFile *pAFile = TFile::Open("gammaJets_pA_Data.root");
  TTree *jets = (TTree*)pAFile->Get("jetTree");
  jets->AddFriend("photonTree","gammaJets_pA_Data.root");

  TH1D *hfp20 = new TH1D("hfp20",";#Delta #phi_{J#gamma}", 50, 0, TMath::Pi());
  TH1D *hfl20 = (TH1D*)hfp20->Clone("hfl20");

  TCut anaCut = "photonTree.hadronicOverEm<0.1 && abs(photonTree.eta)<1.44 && abs(jetTree.eta)<1.6 && jetTree.pt>30 && photonTree.ecalRecHitSumEtConeDR04 <4.2  &&  photonTree.hcalTowerSumEtConeDR04 < 2.2  &&  photonTree.trkSumPtHollowConeDR04 < 2";

  // jets->Project("hfl20","dPhi",anaCut&&"(photonTree.HFplusEta4+photonTree.HFminusEta4 < 20)");
  // jets->Project("hfp20","dPhi",anaCut&&"(photonTree.HFplusEta4+photonTree.HFminusEta4 > 20)");
  jets->Project("hfl20","dPhi",anaCut&&"(photonTree.pt > 40)");
  jets->Project("hfp20","dPhi",anaCut&&"(photonTree.pt > 60)");
  
  hfp20->Scale(1./hfp20->GetEntries());
  hfl20->Scale(1./hfl20->GetEntries());

  TCanvas *c1 = new TCanvas();
  c1->SetLogy();

  // hfp20->Draw();
  // hfp20->GetXaxis()->CenterTitle();
  // hfl20->SetMarkerColor(kBlue);
  // hfl20->SetLineColor(kBlue);
  // hfl20->Draw("same");

  // TLegend *leg = new TLegend(0.15,0.5,0.5,0.8);
  // leg->SetFillColor(0);
  // // leg->AddEntry(hfp20,"E^{HF |#eta|>4} > 20","l");
  // // leg->AddEntry(hfl20,"E^{HF |#eta|>4} < 20","l");
  // leg->AddEntry(hfp20,"#gamma p_{T} > 60GeV","l");
  // leg->AddEntry(hfl20,"#gamma p_{T} > 40GeV","l");
  // leg->Draw();

  TH1D *ratio = (TH1D*)hfp20->Clone("ratio");
  //ratio->SetYTitle("(E^{HF |#eta|>4} > 20)/(E^{HF |#eta|>4} < 20)");
  ratio->Divide(hfp20,hfl20,1,1);
  ratio->Draw();
}