C++ TH1::IntegralAndError方法代码示例

void fill_data(double e_min, double e_max, int mask)
{
	const char name_pattern[] = "danss_report_v4n-sect%d-calc.root";
	const char pos[3][20] = {"hUp_%d", "hMid_%d", "hDown_%d"};
	char fname[1024];
	char cpos[32];
	int i, j;
	TFile *f;
	TH1* h;
	double val, err;
	
	for (i=0; i<3; i++) {
		sprintf(fname, name_pattern, 3-i);
		f = new TFile(fname);
		if (!f->IsOpen()) return;
		for (j=0; j<3; j++) {
			sprintf(cpos, pos[j], mask);
			h = (TH1*) f->Get(cpos);
			if (!h) {
				printf("Something is wrong: %d %d %s.\n", i, j, cpos);
				return;
			}
			val = h->IntegralAndError(h->FindBin(e_min), h->FindBin(e_max), err);
			DataArray.cnt[3*i+j] = val;
			DataArray.ecnt[3*i+j] = err;
			printf("%d %d %f +- %f\n", i, j, val, err);
		}
		f->Close();
	}
}

void RDK2AnalysisPlotter::drawPlotDets(CoDet detType, PlotVarType plotVarType)
{
    gROOT->cd();

    int numDets;

    TString histTitleString;
    if(detType == DET_EPG)
    {
        numDets=12;
        histTitleString="Individual BGO Detectors;Detector Number;% Difference From Mean";
    }
    else
    {
        numDets=3;
        histTitleString="Individual BAPD Detectors;Detector Number;% Difference From Mean";
    }
    const int numHists=numExp+numMC;
    TH1* tempHist;
    double integral;
    double error;

    TH1** detsHists=new TH1*[numHists];
    TString histName;
    double max=-REALLY_BIG_DBL;
    double min=REALLY_BIG_DBL;

    for (int i = 0;i< numHists;i++)  //Set loop (either exp or mc)
    {
        histName="hist"+int2str(plotHists.size());
        detsHists[i]=new TH1D(histName,histTitleString,numDets,.5,numDets+.5);
        detsHists[i]->Sumw2();
        plotHists.push_back(detsHists[i]);

        for (int j = 0;j< numDets;j++)
        {
            if(i<numExp)
            {
                tempHist=getExpHist( i, detType, plotVarType,j+1);
            }
            else
            {
                tempHist=getMCHist( i-numExp, detType, plotVarType,j+1);
            }
            if(tempHist != nullptr)
            {
                integral = tempHist->IntegralAndError(1,tempHist->GetNbinsX(),error);
                detsHists[i]->SetBinContent(j+1,integral);
                detsHists[i]->SetBinError(j+1,error);
            }
        }

        if(i<numExp)
            setDisplayQuantitiesForHist(i,1,detsHists[i],detType,PLOT_DETS,plotVarType);
        else
            setDisplayQuantitiesForHist(-(i-numExp)-1,1,detsHists[i],detType,PLOT_DETS,plotVarType);

        double mean=detsHists[i]->Integral()/numDets; // Mean for the set
        detsHists[i]->Scale(100./mean);
        for (int j = 0;j< numDets;j++)
        {
            detsHists[i]->SetBinContent(j+1,detsHists[i]->GetBinContent(j+1)-100. );
        }
        if(detsHists[i]->GetMaximum() > max)
            max=detsHists[i]->GetMaximum()+detsHists[i]->GetBinError(detsHists[i]->GetMaximumBin());
        if(detsHists[i]->GetMinimum() < min)
            min=detsHists[i]->GetMinimum()-detsHists[i]->GetBinError(detsHists[i]->GetMinimumBin());
    }

    detsHists[0]->GetYaxis()->SetRangeUser(min-.05*abs(min),max+.05*abs(max));
    drawHist(numHists, detsHists, 0, nullptr);
    delete[] detsHists;

}

void EMuSigEff() {

    vector <int> masses;
    masses.push_back(40);
    masses.push_back(50);
    masses.push_back(60);
    masses.push_back(70);
    masses.push_back(80);
    masses.push_back(90);
    masses.push_back(100);
    masses.push_back(125);
    masses.push_back(150);
    masses.push_back(175);
    masses.push_back(200);
    masses.push_back(250);
    masses.push_back(300);
    masses.push_back(350);
    masses.push_back(400);
    masses.push_back(500);

    vector <TString> smasses;
    smasses.push_back("40");
    smasses.push_back("50");
    smasses.push_back("60");
    smasses.push_back("70");
    smasses.push_back("80");
    smasses.push_back("90");
    smasses.push_back("100");
    smasses.push_back("125");
    smasses.push_back("150");
    smasses.push_back("175");
    smasses.push_back("200");
    smasses.push_back("250");
    smasses.push_back("300");
    smasses.push_back("350");
    smasses.push_back("400");
    smasses.push_back("500");


    for(unsigned int i = 0 ; i < masses.size(); ++i) {

        TString im = smasses.at(i);

        TFile * file1 = new TFile(("/home/jalmond/HeavyNeutrino/Analysis/LQanalyzer/data/output/SSElectronMuon/HNEMu_SKHNmue" + im + "_nocut_5_3_14.root").Data());
        TFile * file2 = new TFile(("/home/jalmond/HeavyNeutrino/Analysis/LQanalyzer/data/output/SSElectronMuon/HNEMu_SKHNemu" + im + "_nocut_5_3_14.root").Data());

        TString cut = "SS_highmass_exc";// + im;
        //cut = "SS_highmass";
        //    if(i < 5) cut = "SS_lowmass_" + im;
        //if(i < 5) cut = "SS_lowmass_80_2";
        TString hist = (cut + "/h_Nelectrons_"+cut);

        TH1* hnsig =   (TH1F*)file1->Get(("CutFlow/_eventcutflow"));

        TH1* hnsig2 =   (TH1F*)file2->Get(("CutFlow/_eventcutflow"));

        float nsig = float(hnsig->GetBinContent(2));
        float nsig2= float(hnsig2->GetBinContent(2));

        TH1*  hpass = (TH1F*)file1->Get(hist);
        TH1*  hpass2 = (TH1F*)file2->Get(hist);

        cout  << " \n ------- " << endl;
        cout  << " Mass = " << masses.at(i) << endl;
        cout << "mu eacceptance = " << hpass->Integral()/nsig << endl;
        cout << "emu acceptance = " << hpass2->Integral()/nsig2 << endl;

        hpass->Add(hpass2,1.);
        double err ;
        hpass->IntegralAndError(1, hpass->GetNbinsX()+1, err    , "");
        cout << "Average = " << hpass->Integral()/(nsig + nsig2) << endl;
        cout << "%err = " << 100*err/(hpass->Integral())<< endl;;

    }
}