当前位置: 首页>>代码示例>>C++>>正文

C++ TH1F::Chi2Test方法代码示例

本文整理汇总了C++中TH1F::Chi2Test方法的典型用法代码示例。如果您正苦于以下问题：C++ TH1F::Chi2Test方法的具体用法？C++ TH1F::Chi2Test怎么用？C++ TH1F::Chi2Test使用的例子？那么恭喜您, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在类TH1F的用法示例。


在下文中一共展示了TH1F::Chi2Test方法的6个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于系统推荐出更棒的C++代码示例。

示例1: ComputeMin

void Analyzer::ComputeMin(){
	g2=new TGraph2D(); //TODO
	
	alpha=1.0;beta=0;
	Loop(t_data,1);
	if(varName=="QGLMLP")
		Loop(t_mc,4);
	//scan
	alpha=1.0;beta=0;
	for(float ai=0.7; ai<=1.1; ai+=0.02)
		{
		Reset(h_mc);	
		alpha=ai;
		Loop(t_mc,2);
		h_mc->Scale(h_data->Integral()/h_mc->Integral());
		g2->SetPoint(g2->GetN(),alpha,beta, h_data->Chi2Test(h_mc,opt.c_str())  );	
		}
	alpha=1.0;beta=0;
	for(float bi=-0.5; bi<=0.5; bi+=0.01)
		{
		Reset(h_mc);	
		beta=bi;
		Loop(t_mc,2);
		h_mc->Scale(h_data->Integral()/h_mc->Integral());
		g2->SetPoint(g2->GetN(),alpha,beta, h_data->Chi2Test(h_mc,opt.c_str())  );	
		}
	
	//Find min0;min1
	float min0=1,min1=0;
	pair<float,float> R=MinG(g2);
	min0=R.first;min1=R.second;

	for(int i=-nstep;i<=nstep;i++)
	for(int j=-nstep;j<=nstep;j++)
        	{
        	alpha=min0+i*stp0;
        	beta=min1+j*stp1;
		Loop(t_mc,2);
		h_mc->Scale(h_data->Integral()/h_mc->Integral());
		g2->SetPoint(g2->GetN(),alpha,beta, h_data->Chi2Test(h_mc,opt.c_str())  );
		}
	//double min0,min1;	
	R=MinG(g2);
	printf("a=%.3f;b=%.3f;lmin=%.3f;lmax=%.3f;break;\n",R.first,R.second,lmin,lmax);
	return;
	
}
开发者ID:amarini，项目名称:UserCode，代码行数:47，代码来源:ComputeDoubleMinDiJet.C

示例2: 

pair<float,float> Analyzer::SmearDoubleMinFast(float a0_q,float b0_q , float a0_g,float b0_g,int type,int WriteOut){ //type = 0 Q, 1 G

	TGraph2D *g2_q=new TGraph2D(); g2_q->SetName("g2_q");
	TGraph2D *g2_g=new TGraph2D(); g2_g->SetName("g2_g");
	
	//scan
	a_q=a0_q;b_q=b0_q;
	a_g=a0_g;b_g=b0_g;

	alpha=1.0;beta=0;
	for(float ai=0.7; ai<=1.1; ai+=0.02)
		{
		Reset(h_mc);	
		if(type==0)a_q=ai;
		if(type==1)a_g=ai;
		LoopFast();	
		h_mc->Scale(h_data->Integral()/h_mc->Integral());

		if(type==0)g2_q->SetPoint(g2_q->GetN(),a_q,b_q, h_data->Chi2Test(h_mc,opt.c_str())  );	
		if(type==1)g2_g->SetPoint(g2_g->GetN(),a_g,b_g, h_data->Chi2Test(h_mc,opt.c_str())  );	
		}
	alpha=1.0;beta=0;
	a_q=a0_q;b_q=b0_q;
	a_g=a0_g;b_g=b0_g;

	for(float bi=-0.5; bi<=0.5; bi+=0.01)
		{
		Reset(h_mc);	
		if(type==0)b_q=bi;
		if(type==1)b_g=bi;
		LoopFast();
		h_mc->Scale(h_data->Integral()/h_mc->Integral());
		if(type==0)g2_q->SetPoint(g2_q->GetN(),a_q,b_q, h_data->Chi2Test(h_mc,opt.c_str())  );	
		if(type==1)g2_g->SetPoint(g2_g->GetN(),a_g,b_g, h_data->Chi2Test(h_mc,opt.c_str())  );	
		}
	
	//Find min0;min1
	float min0=1,min1=0;
	pair<float,float> R;
		if(type==0)R=MinG(g2_q);
		if(type==1)R=MinG(g2_g);
	min0=R.first;min1=R.second;

	for(int i=-nstep;i<=nstep;i++)
	for(int j=-nstep;j<=nstep;j++)
        	{
        	if(type==0)a_q=min0+i*stp0;
        	if(type==1)a_g=min0+i*stp0;
        	if(type==0)b_q=min1+j*stp1;
        	if(type==1)b_g=min1+j*stp1;
		LoopFast();
		h_mc->Scale(h_data->Integral()/h_mc->Integral());
		if(type==0)g2_q->SetPoint(g2_q->GetN(),a_q,b_q, h_data->Chi2Test(h_mc,opt.c_str())  );	
		if(type==1)g2_g->SetPoint(g2_g->GetN(),a_g,b_g, h_data->Chi2Test(h_mc,opt.c_str())  );	
		}
	
		if(type==0)R=MinG(g2_q);
		if(type==1)R=MinG(g2_g);
	//SAME ON G,& REDO
	//printf("a=%.3f;b=%.3f;lmin=%.3f;lmax=%.3f;break;\n",R.first,R.second,lmin,lmax);
	
	if(WriteOut){	
	string name=Form("Results/outputZJet2_%s_pt%.0f_%.0f_rho%.0f_%.0f_eta%.0f_%.0f",varName.c_str(),PtMin,PtMax,RhoMin,RhoMax,EtaMin,EtaMax);
	if(type==0)g2_q->SaveAs((name+"g2_q.root").c_str());
	if(type==1)g2_g->SaveAs((name+"g2_g.root").c_str());
	}
	
	return R;
}
开发者ID:amarini，项目名称:UserCode，代码行数:69，代码来源:ComputeDoubleMinDiJet.C

示例3: ComputeMinFast

void Analyzer::ComputeMinFast(){
	g2=new TGraph2D(); 
	
	alpha=1.0;beta=0;
	Loop(t_data,1);
	for(int j=0;j<=h_data->GetNbinsX()+1;j++)if(h_data->GetBinError(j)==0)h_data->SetBinError(j,1);
	if(varName=="QGLMLP")
		Loop(t_mc,4);
	//scan
	//reset Fast
	ResetFast();
	
	alpha=1.0;beta=0;
	for(float ai=aMin; ai<=aMax; ai+=0.02)
		{
		Reset(h_mc);	
		alphaFast.push_back(ai);	
		betaFast.push_back(0);	
		h_mcFast.push_back(new TH1F( Form("hmc_%d",int(h_mcFast.size())),"hmc",nBins,xMin,xMax) );
		}
	alpha=1.0;beta=0;
	for(float bi=bMin; bi<=bMax; bi+=0.01)
		{
		Reset(h_mc);	
		alphaFast.push_back(1.0);	
		betaFast.push_back(bi);	
		h_mcFast.push_back(new TH1F( Form("hmc_%d",int(h_mcFast.size())),"hmc",nBins,xMin,xMax) );
		}

	for(int j=0;j< int(h_mcFast.size());j++) h_mcFast[j]->Sumw2();	
	Loop(t_mc,16);
		for(int i=0 ;i<int(alphaFast.size());i++)
			{
			for(int j=0;j<=h_mcFast[i]->GetNbinsX()+1;j++)if(h_mcFast[i]->GetBinError(j)==0)h_mcFast[i]->SetBinError(j,1);
			h_mcFast[i]->Scale(h_data->Integral()/h_mcFast[i]->Integral());
			g2->SetPoint(g2->GetN(),alphaFast[i],betaFast[i], h_data->Chi2Test(h_mcFast[i],opt.c_str())  );
			}
	//Find min0;min1
	float min0=1,min1=0;
	pair<float,float> R=MinG(g2);
	min0=R.first;min1=R.second;
		
	ResetFast();

	delete g2;
	g2=new TGraph2D();

	for(int i=-nstep;i<=nstep;i++)
	for(int j=-nstep;j<=nstep;j++)
        	{
        	alphaFast.push_back(min0+i*stp0 );
        	betaFast.push_back(min1+j*stp1 );
		h_mcFast.push_back(new TH1F( Form("hmc_%d",int(h_mcFast.size())),"hmc",nBins,xMin,xMax) );
		//g2->SetPoint(g2->GetN(),alpha,beta, h_data->Chi2Test(h_mc,opt.c_str())  );
		}
	alphaFast.push_back(min0);
	betaFast.push_back(0);
	h_mcFast.push_back(new TH1F( Form("hmc_%d",int(h_mcFast.size())),"hmc",nBins,xMin,xMax) );

	alphaFast.push_back(min1);
	betaFast.push_back(bMax);
	h_mcFast.push_back(new TH1F( Form("hmc_%d",int(h_mcFast.size())),"hmc",nBins,xMin,xMax) );
	
	alphaFast.push_back(1);
	betaFast.push_back(0);
	h_mcFast.push_back(new TH1F( Form("hmc_%d",int(h_mcFast.size())),"hmc",nBins,xMin,xMax) );
	
	for(int j=0;j<int(h_mcFast.size());j++) h_mcFast[j]->Sumw2();	
	Loop(t_mc,16);
		for(int i=0 ;i<int(alphaFast.size());i++)
			{
			for(int j=0;j<=h_mcFast[i]->GetNbinsX()+1;j++)if(h_mcFast[i]->GetBinError(j)==0)h_mcFast[i]->SetBinError(j,1);
			h_mcFast[i]->Scale(h_data->Integral()/h_mcFast[i]->Integral());
			g2->SetPoint(g2->GetN(),alphaFast[i],betaFast[i], h_data->Chi2Test(h_mcFast[i],opt.c_str())  );
			}
	double m0,m1;
	R=MinG(g2,&m0,&m1);
	printf("a=%.3f;b=%.3f;lmin=%.3f;lmax=%.3f;break;//chi2=%.3lf; chi2_0=%.3lf\n",R.first,R.second,lmin,lmax,m0,m1);
	{
	TFile *out=TFile::Open("output.root","UPDATE");out->cd();
	for(int i=0;i<int(h_mcFast.size());i++)
		{
		h_mcFast[i]->SetName(Form("%s_alpha%.2f_beta%.2f_lmin%.3f_lmax%.3f_pt%0f_%.0f_rho%.0f_%.0f_eta%.0f_%.0f",varName.c_str(),alphaFast[i],betaFast[i],lmin,lmax,PtMin,PtMax,RhoMin,RhoMax,EtaMin,EtaMax));
		h_mcFast[i]->Write();
		}
	h_data->SetName(Form("%s_data_pt%0f_%.0f_rho%.0f_%.0f_eta%.0f_%.0f",varName.c_str(),PtMin,PtMax,RhoMin,RhoMax,EtaMin,EtaMax));
	h_data->Write();
	}
	return;
}
开发者ID:amarini，项目名称:UserCode，代码行数:90，代码来源:ComputeDoubleMinDiJet.C

示例4: main


//.........这里部分代码省略.........
    double NewMax = max;

    int rebinning = Nbins; 

    if (RMS>0 && hd->GetType() ) 
      {
	dNdX = 100. / ( 10 * RMS);
	NewMin = Mean - 10 * RMS;
	NewMax = Mean + 10 * RMS;
      }
    
    if ((dX * dNdX)>0  && hd->GetType() ) 
      rebinning = (int)(double(Nbins) / (dX * dNdX));
    
    if ( rebinning > 1 && hd->GetType() ) 
      {
	href->Rebin(rebinning);
	hnew->Rebin(rebinning);
      }

    if ( hd->GetType() == 1 )
      { 
	href->GetXaxis()->SetRangeUser(0.0, NewMax);
	hnew->GetXaxis()->SetRangeUser(0.0, NewMax);
      }
    else if ( hd->GetType() == 2 )
      {
	href->GetXaxis()->SetRangeUser(NewMin, NewMax);
	hnew->GetXaxis()->SetRangeUser(NewMin, NewMax);
      }

    // perform statistical tests
    double ks_score = hnew->KolmogorovTest(href,"D");
    double chi2_score = hnew->Chi2Test(href, "p");
    //double result = KS_TEST ? ks_score : chi2_score;
    double result = (ks_score>chi2_score) ? ks_score : chi2_score;
    
      href->SetNormFactor(Nevents_new);     	
      hnew->SetNormFactor(Nevents_new);
    //hnew->SetNormFactor(1);

    // ensure that the peaks of both histograms will be shown by making a dummy histogram
    float Nentries_ref = href->GetEntries();
    float Nentries_new = hnew->GetEntries();
    float XaxisMin_ref = 0, XaxisMax_ref = 0, YaxisMin_ref = 0, YaxisMax_ref = 0;
    float XaxisMin_new = 0, XaxisMax_new = 0, YaxisMin_new = 0, YaxisMax_new = 0;
    if (Nentries_ref>0) YaxisMax_ref = (href->GetMaximum()+TMath::Sqrt(href->GetMaximum()))*(Nentries_new/Nentries_ref);
    if (Nentries_new>0) YaxisMax_new = (hnew->GetMaximum()+TMath::Sqrt(hnew->GetMaximum()));

    XaxisMin_ref = href->GetXaxis()->GetXmin()>NewMin  ? href->GetXaxis()->GetXmin() : NewMin;
    XaxisMax_ref = href->GetXaxis()->GetXmax()<=NewMax ? href->GetXaxis()->GetXmax() : NewMax;
    YaxisMax_ref = (YaxisMax_ref>=YaxisMax_new) ? YaxisMax_ref : YaxisMax_new;

    if (TMath::Abs(XaxisMin_ref - XaxisMax_ref)<1E-6)
      {
	XaxisMin_ref = 0;
	XaxisMax_ref = 1;
      }
    
    TH1F *hdumb = new TH1F("hdumb","", rebinning, XaxisMin_ref, XaxisMax_ref);
    hdumb->SetMinimum(1E-1); //--For Rick
    hdumb->SetMaximum(1.05*YaxisMax_ref);
    //    if (href->GetMaximum() < hnew->GetMaximum())
    //  href->SetAxisRange(0, 1.1 * hnew->GetMaximum(), "Y");
        
    // set drawing options on the reference histogram
开发者ID:aashaqshah，项目名称:cmssw-1，代码行数:67，代码来源:htMetplotCompare.cpp

示例5: TCanvas


//.........这里部分代码省略.........

  //d->SaveAs("tmp/d_hist.root");
  s->SaveAs("tmp/s_hist.root");
  s2->SaveAs("tmp/s2_hist.root");


  yLabel.Form("Events /(%.2f %s)", s->GetBinWidth(2), yLabelUnit.Data());
  float max = 1.1 * std::max(
			     d->GetMaximum(),///d->Integral(),
			     s->GetMaximum() ///s->Integral()
			     );
  max=1.1*d->GetMaximum();
  std::cout << "max = " << max << std::endl;
  std::cout << "nEvents data: " << d->Integral() << "\t" << d->GetEntries() << std::endl;
  std::cout << "nEvents signal: " << s->Integral() << "\t" << s->GetEntries() << std::endl;
  std::cout << "nEvents signal2: " << s2->Integral() << "\t" << s2->GetEntries() << std::endl;
  if(logy){
    max*=10;
    d->GetYaxis()->SetRangeUser(0.1,max);
    s->GetYaxis()->SetRangeUser(0.1,max);
    s2->GetYaxis()->SetRangeUser(0.1,max);
    c->SetLogy();
  } else {
    d->GetYaxis()->SetRangeUser(0,max);
    s->GetYaxis()->SetRangeUser(0,max);
    s2->GetYaxis()->SetRangeUser(0,max);
  }
  s->GetYaxis()->SetTitle(yLabel);
  s->GetXaxis()->SetTitle(xLabel);
  s2->GetYaxis()->SetTitle(yLabel);
  s2->GetXaxis()->SetTitle(xLabel);
  d->GetYaxis()->SetTitle(yLabel);
  d->GetXaxis()->SetTitle(xLabel);


  d->SetMarkerStyle(20);
  d->SetMarkerSize(1);
  if(d != s){
    s->SetMarkerStyle(20);
    s->SetMarkerSize(1);
    s->SetFillStyle(3001);
    s->SetFillColor(kRed);
  }
  if(s2 != s){
    //s2->SetMarkerStyle(1);
    //s2->SetMarkerSize(0);
    //s->SetFillStyle(0);
    //s->SetFillColor(kB);
    s2->SetLineWidth(3);
    s2->SetLineColor(kBlack);
  }


  TH1F* s_norm = (TH1F *) (s->DrawNormalized("hist", d->Integral()));
  TH1F* s2_norm = (TH1F *) (s2->DrawNormalized("hist same", d->Integral()));
  //TH1F* d_norm = s_norm;
  //if(d!=s) d_norm = (TH1F *) (d->DrawNormalized("p same", d->Integral()));
  if(d!=s) d->Draw("p same");

  if(logy){
    //d_norm->GetYaxis()->SetRangeUser(0.1,max);
    s_norm->GetYaxis()->SetRangeUser(0.1,max);
    c->SetLogy();
  } else {
    //d_norm->GetYaxis()->SetRangeUser(0,max);  
    s_norm->GetYaxis()->SetRangeUser(0,max);  
  }
  std::cout << "Variable  & Data & Simulation & Simulation2 \\" << std::endl;
  std::cout << "Mean      & " << d->GetMean() << " " << d->GetMeanError() 
	    << " & " << s_norm->GetMean() <<  " " << s_norm->GetMeanError() 
	    << " & " << s2_norm->GetMean() <<  " " << s2_norm->GetMeanError() 
	    << " \\" << std::endl;
  std::cout << "Std. dev. & " << d->GetRMS() << " " << d->GetRMSError() 
	    << " & " << s_norm->GetRMS() << " " << s_norm->GetRMSError() 
	    << " & " << s2_norm->GetRMS() << " " << s2_norm->GetRMSError() 
	    << " \\" << std::endl;
  std::cout << "\\hline" << std::endl;
  std::cout << "$\\Chi^2$ " <<  d->Chi2Test(s_norm, "UW CHI2/NDF NORM") << std::endl;
  

  TLegend *leg = new TLegend(0.6,0.8,1,1);
  if(dataLabel !="") leg->AddEntry(d,dataLabel,"p");
  if(mcLabel   !="") leg->AddEntry(s,mcLabel, "lf");
  if(mc2Label   !="") leg->AddEntry(s2,mc2Label, "l");
  leg->SetBorderSize(1);
  leg->SetFillColor(0);
  leg->SetTextSize(0.04);
  if(dataLabel !="" && mcLabel !="") leg->Draw();
  //c->GetListOfPrimitives()->Add(leg,"");

  TPaveText *pv = new TPaveText(0.23,0.95,0.6,1,"NDC");
  pv->AddText("CMS Preliminary 2016");
  pv->SetFillColor(0);
  pv->SetBorderSize(0);
  pv->Draw();


  return c;

}
开发者ID:GiuseppeFasanella，项目名称:ECALELF，代码行数:101，代码来源:PlotDataMC.C

示例6: ReactorNuAnalysis

void ReactorNuAnalysis()
{
	// -------- VARIABLE DEFINITIONS and SETUP --------
	// Style settings
	gStyle->SetOptStat("");
	
	// Constants
	Double_t NuSpectrumMinE = 0.0;  // [MeV]
	Double_t NuSpectrumMaxE = 10.0; // [MeV]
	Int_t seed = 43534;
	const Double_t Gfermi = 1.16637e-11;																			 // [MeV^-2]
	const Double_t Sin2ThetaW = 0.2387;
	const Double_t InverseMeVtoCm = 1.97e-11;
	const Double_t elementaryCharge = 1.602176565e-19;
	
	// Spectrum files (these are just samples for now)
	const char ReactorNeutronBackgroundFile[] = "SampleData.txt";
	
	// Define constants relevant for this run
	Double_t OnOffTimeRatio = 1.0/1.0;
	Double_t time = 100 * 24.0*3600.0;																				 // [sec]
	Double_t detMass = 5000.0;																								 // [g]
	Double_t distance = 200.0;																								 // [cm]
	Double_t activity = 6.0/200.0/elementaryCharge;														 // [sec^-1]
	const Int_t nNeutrons = 14;
	const Int_t nProtons = 14;
	const Double_t Qweak = nNeutrons - (1 - 4*Sin2ThetaW) * nProtons;
	const Double_t NucleonMass = (nNeutrons * 939.565) + (nProtons * 938.272);			 // [MeV]
	
	// Define the histograms
	TH1F* RecoilEvtHistogram = new TH1F("RecoilEvtHistogram","^{28}Si recoil energy spectrum;Energy [eV];Events / 10 eV",50,0.0,500.0);
	TH1F* NeutrinoEvtHistogram = new TH1F("NeutrinoEvtHistogram","#bar{#nu}_{e} energy spectrum;Energy [MeV];Events / 0.5 MeV",50,0.0,20);
	TH1F* TheoryRecoilEvtHistogram = new TH1F("TheoryRecoilEvtHistogram","High-statistics (\"theoretical\") Coherent Recoil Spectrum;Energy [eV];Events / 10 eV",50,0.0,500.0);
	TH1F* TheoryNeutrinoEvtHistogram = new TH1F("TheoryNeutrinoEvtHistogram","High-statistics (\"theoretical\") Neutrino Energy Spectrum;Energy [MeV] / 0.5 MeV;Events",50,0.0,20);
	TH1F* TheoryRecoilEvtHistogramFit = new TH1F("TheoryRecoilEvtHistogram","MC Fit;Events",50,0.0,0.002);
	TH1F* EMNoLukeBackground = new TH1F("EMNoLukeBackground","Background from EM recoils (before Luke Effect amplification;Energy [MeV];Events)",50,0.0,0.002);
	TH1F* ReactorNeutronBackground = new TH1F("ReactorNeutronBackground","Background from reactor neutrons;Energy [MeV];Events",50,0.0,0.002);
	TH1F* ReactorOnCosmoNeutronBackground = new TH1F("ReactorOnCosmoNeutronBackground","Reactor-on background from muon-induced neutrons;Energy [MeV];Events)",50,0.0,0.002);
	TH1F* ReactorOffCosmoNeutronBackground = new TH1F("ReactorOffCosmoNeutronBackground","Reactor-off background from muon-induced neutrons;Energy [MeV];Events)",50,0.0,0.002);	
	TH1F* ReactorOnHisto = new TH1F("ReactorOnHisto","Recoil Spectrum for Reactor-On Data;Energy [MeV];Events",50,0.0,0.002);
	TH1F* ReactorOffHisto = new TH1F("ReactorOffHisto","Recoil Spectrum for Reactor-Off Data;Energy [MeV];Events",50,0.0,0.002);
	TH1F* BackgroundSubtractedSignal = new TH1F("BackgroundSubtractedSignal","Recoil Spectrum for Reactor-Off Data;Energy [MeV];Events",50,0.0,0.002);
	
	// Energy spectra
	// (Spectral parameterizations from arXiv:1101.2663v3)
	TF1* NeutrinoEnergySpectrum = new TF1("NeutrinoSpectrum", "TMath::Exp([0] + [1]*x + [2]*TMath::Power(x,2) + [3]*TMath::Power(x,3) + [4]*TMath::Power(x,4) + [5]*TMath::Power(x,5))", NuSpectrumMinE, NuSpectrumMaxE);
	NeutrinoEnergySpectrum->SetParameters(3.217, -3.111, 1.395, -0.369, 0.04445, -0.002053);
	TF1* IntegratedRecoilSpectrum = new TF1("IntegratedRecoilSpectrum", "TMath::Power([0]*[1]*[3],2)/(4*TMath::Pi()) * [2] * (x/(1 + [2]/(2*x))) * (1 - ([2]/(4*x*x)) * (x/(1 + [2]/(2*x))))", NuSpectrumMinE, NuSpectrumMaxE);
	IntegratedRecoilSpectrum->SetParameters(Gfermi, Qweak, NucleonMass, InverseMeVtoCm);
	TF1* RecoilSpectrum = new TF1("RecoilSpectrum","IntegratedRecoilSpectrum * NeutrinoSpectrum", 0.0, 10.0);
	TF1* DiffRecoilSpectrumAtConstE = new TF1("DiffRecoilSpectrumAtConstE", "(x<[4])*TMath::Power([0]*[1],2)/(4*TMath::Pi()) * [2] * (1 - ([2] * x)/(2 * TMath::Power([3],2))) + (x>[4])*0", NuSpectrumMinE, 0.01);
	DiffRecoilSpectrumAtConstE->SetParameters(Gfermi, Qweak, NucleonMass);
	
	
	
	
	// -------- HISTOGRAM FILLING --------
	// Fill "experimental" histograms
	Int_t nEvt = GenerateNumOfNuRecoils(time, detMass, distance, activity, nNeutrons, nProtons, RecoilSpectrum, NuSpectrumMinE, NuSpectrumMaxE, seed);
	FillNuRecoilSpectrum(RecoilEvtHistogram, NeutrinoEvtHistogram, nEvt, nNeutrons, nProtons, RecoilSpectrum, DiffRecoilSpectrumAtConstE, NuSpectrumMinE, NuSpectrumMaxE, seed+4);
	FillRecoilSpectrumFromFile(ReactorNeutronBackground, 100.0, 10.0, ReactorNeutronBackgroundFile, seed+1);      // Testing purposes only.  CHANGE ME!!
	FillRecoilSpectrumFromFile(ReactorOnCosmoNeutronBackground, 50.0, 10.0, ReactorNeutronBackgroundFile, seed+2);		// Testing purposes only.  CHANGE ME!!
	FillRecoilSpectrumFromFile(ReactorOffCosmoNeutronBackground, 50.0, 10.0, ReactorNeutronBackgroundFile, seed+3);		// Testing purposes only.  CHANGE ME!!
	
	cout << "nEvt: " << nEvt << endl;
	
	// Fill high-statistics "theoretical" histograms
	FillNuRecoilSpectrum(TheoryRecoilEvtHistogram, TheoryNeutrinoEvtHistogram, 10000, nNeutrons, nProtons, RecoilSpectrum, DiffRecoilSpectrumAtConstE, NuSpectrumMinE, NuSpectrumMaxE, seed+5);
	TheoryNeutrinoEvtHistogram->Scale(nEvt/TheoryNeutrinoEvtHistogram->GetEntries());
	TheoryRecoilEvtHistogram->Scale(nEvt/TheoryNeutrinoEvtHistogram->GetEntries());
	
	// Combine the histograms into total
	ReactorOnHisto->Add(RecoilEvtHistogram);
	ReactorOnHisto->Add(ReactorNeutronBackground);
	ReactorOnHisto->Add(ReactorOnCosmoNeutronBackground);
	
	ReactorOffHisto->Add(ReactorOffCosmoNeutronBackground);
	
	
	
	
	// -------- HYPOTHESIS TESTING --------
	cout << "p-value between Reactor-On and Reactor-Off Data: " << ReactorOnHisto->Chi2Test(ReactorOffHisto) << endl;
	cout << "p-value between the simulated data and the Monte Carlo histogram: " << RecoilEvtHistogram->Chi2Test(TheoryRecoilEvtHistogram) << endl;
	
	
	
	
	// -------- BACKGROUND SUBTRACTION and FITTING --------
	// Normalize reactor-off data to reactor-on data by exposure
	BackgroundSubtractedSignal->Add(ReactorOnHisto, ReactorOffHisto, 1.0, -1.0*OnOffTimeRatio);
	
	// Use TFractionFitter to do fitting
	TObjArray *FractionFitData = new TObjArray(2);
	FractionFitData->Add(TheoryRecoilEvtHistogram);
	FractionFitData->Add(ReactorOffHisto);
	TFractionFitter* ffit = new TFractionFitter(ReactorOnHisto, FractionFitData);
	ffit->Constrain(0,0.1,10.0);
	ffit->Constrain(1,1.0,1.0);
	Int_t status = ffit->Fit();
//.........这里部分代码省略.........
开发者ID:spitzj，项目名称:RicochetMC，代码行数:101，代码来源:ReactorNuAnalysis.C


注：本文中的TH1F::Chi2Test方法示例由纯净天空整理自Github/MSDocs等开源代码及文档管理平台，相关代码片段筛选自各路编程大神贡献的开源项目，源码版权归原作者所有，传播和使用请参考对应项目的License；未经允许，请勿转载。