C++ TFitResultPtr类代码示例

Double_t dExpFunc(const TF1 *fcn, const Double_t x, const TFitResultPtr fs) {
  Double_t df[3];
  Double_t a = fcn->GetParameter(0);
  Double_t b = fcn->GetParameter(1);
  Double_t c = fcn->GetParameter(2);

  TF1 *deriv_par0 = new TF1("dfdp0",df_dParExp,130,1500,4);
  deriv_par0->SetParameters(0,a,b,c); // This will set the derivative for the first parameter.
  TF1 *deriv_par1 = new TF1("dfdp1",df_dParExp,130,1500,4);
  deriv_par1->SetParameters(1,a,b,c); // This will set the derivative for the second parameter
  TF1 *deriv_par2 = new TF1("dfdp2",df_dParExp,130,1500,4);
  deriv_par2->SetParameters(2,a,b,c); // This will set the derivative for the 3rd parameter
  df[0] = deriv_par0->Eval(x);
  df[1] = deriv_par1->Eval(x);
  df[2] = deriv_par2->Eval(x);
  Double_t err2=0;
  for(Int_t i=0; i<3; i++) {
    err2 += df[i]*df[i]*(fs->GetCovarianceMatrix()[i][i]);
    for(Int_t j=i+1; j<3; j++) {
      err2 += 2.0*df[i]*df[j]*(fs->GetCovarianceMatrix()[i][j]);
    }
  }
  assert(err2>=0);
  return sqrt(err2);
}

std::pair<double,double> CalibTree::fitMean(TH1D* hist, int mode) {
  std::pair<double,double> results = std::pair<double,double>(1,0);
  if (hist != 0) {
    double mean = hist->GetMean(), rms = hist->GetRMS();
    double LowEdge(0.7), HighEdge(1.3);
    char   option[20];
    if (mode == 1) {
      LowEdge  = mean - 1.5*rms;
      HighEdge = mean + 1.5*rms;
      int nbin = hist->GetNbinsX();
      if (LowEdge < hist->GetBinLowEdge(1)) LowEdge = hist->GetBinLowEdge(1);
      if (HighEdge > hist->GetBinLowEdge(nbin)+hist->GetBinWidth(nbin))
	HighEdge = hist->GetBinLowEdge(nbin)+hist->GetBinWidth(nbin);
    }
    if (hist->GetEntries() > 100) sprintf (option, "+QRLS");
    else                          sprintf (option, "+QRWLS");
    double value(mean);
    double error = rms/sqrt(hist->GetEntries());
    if (hist->GetEntries() > 20) {
      TFitResultPtr Fit = hist->Fit("gaus",option,"",LowEdge,HighEdge);
      value    = Fit->Value(1);
      error    = Fit->FitResult::Error(1); 
      /*
      LowEdge  = value - 1.5*error;
      HighEdge = value + 1.5*error;
      value    = Fit->Value(1);
      error    = Fit->FitResult::Error(1); 
      Fit = hist->Fit("gaus",option,"",LowEdge,HighEdge);
      */
    }
    results =  std::pair<double,double>(value,error);
  }
  return results;
}

Double_t dpol1Func(const TF1 *fcn, const Double_t x, const TFitResultPtr fs) {
  Double_t df[2];
  Double_t a = fcn->GetParameter(0);
  Double_t b = fcn->GetParameter(1);
  
  df[0] = 1;
  df[1] = x;
  Double_t err2 = df[0]*df[0]*(fs->GetCovarianceMatrix()[0][0])
                  + df[1]*df[1]*(fs->GetCovarianceMatrix()[1][1])
		  + 2.0*df[0]*df[1]*(fs->GetCovarianceMatrix()[0][1]);
  assert(err2>=0);
  return sqrt(err2);
}

void CalibTree::fitPol0(TH1D* hist, bool debug) {

  hist->GetXaxis()->SetTitle("i#eta");
  hist->GetYaxis()->SetTitle("<E_{HCAL}/(p-E_{ECAL})>");
  hist->GetYaxis()->SetRangeUser(0.4,1.6);
  TFitResultPtr Fit = hist->Fit("pol0","+QRWLS");
  if (debug) {
    std::cout << "Fit to Pol0 to " << hist->GetTitle() << ": " 
	      << Fit->Value(0) << " +- " << Fit->FitResult::Error(0)
	      << std::endl;
  }
  hist->Write();
}

TGraphErrors* ExtractCalCurve(vector<calData>& FitPositions, vector<Float_t>& Energies)
//void ExtractCalCurve(vector<Float_t>& FitPositions, vector<Float_t>& Energies)
{
    //TF1* calFit = new TF1("calFit","[0]*x + [1]",0,800E3);
    TGraphErrors* calCurve = new TGraphErrors(Energies.size());
    
    //Get highest energy peak and estimate [0]
    std::vector<calData>::iterator pit;
    std::vector<Float_t>::iterator eit;
    
    pit = max_element(FitPositions.begin(),FitPositions.end(),CompareByadc);
    eit = max_element(Energies.begin(),Energies.end());
    
    //Extrapolate coefficient and place first point
    int currPoint = 1;
    calCurve->SetPoint(currPoint,(*pit).fadc,*eit);
    calCurve->SetPointError(currPoint,(*pit).efadc,0.0);
    TFitResultPtr r = calCurve->Fit("pol1","SQ");   
    Float_t a = r->Parameter(1);

    Float_t CurrentPeak = 0.0;
    Float_t CurrentEnergy = 0.0;
    Float_t CurrentEnergyEst = 0.0;
    // Loop through found peaks and locate closest estimated energy
    // Assume fitted peaks are already ordered from lowest to highest
    for(std::vector<calData>::iterator i = --(FitPositions.end()); i!=FitPositions.begin(); --i)
    {
        currPoint++;
        CurrentPeak = (*i).fadc;
        CurrentEnergyEst = CurrentPeak*a;
        Float_t CurrentDelta = 800E3;
        for(std::vector<Float_t>::iterator j = Energies.begin(); j!=Energies.end(); j++)
	{
		if( abs(*j - CurrentEnergyEst) < CurrentDelta)
		{
			CurrentDelta = abs(*j - CurrentEnergyEst);
                        CurrentEnergy = *j;
		}
	} 
	
	calCurve->SetPoint(currPoint,CurrentPeak,CurrentEnergy);
        calCurve->SetPointError(currPoint,(*i).efadc,CurrentDelta);
        r = calCurve->Fit("pol1","SQ");
        a = r->Parameter(1);
    }
    r->Print("V");
    return calCurve; 
}

Double_t dsqrtFunc(const TF1 *fcn, const Double_t x, const TFitResultPtr fs) {
  Double_t df[2];
  Double_t a = fcn->GetParameter(0);
  Double_t b = fcn->GetParameter(1);

  TF1 *deriv_par0 = new TF1("dfdp0",df_dParsqrt,130,1620,3);
  deriv_par0->SetParameters(0,a,b); // This will set the derivative for the first parameter.
  TF1 *deriv_par1 = new TF1("dfdp1",df_dParsqrt,130,1620,3);
  deriv_par1->SetParameters(1,a,b); // This will set the derivative for the second parameter
  df[0] = deriv_par0->Eval(x);
  df[1] = deriv_par1->Eval(x);
  Double_t err2 = df[0]*df[0]*(fs->GetCovarianceMatrix()[0][0])
                  + df[1]*df[1]*(fs->GetCovarianceMatrix()[1][1])
		  + 2.0*df[0]*df[1]*(fs->GetCovarianceMatrix()[0][1]);
  assert(err2>=0);
  return sqrt(err2);
}

void PlotERes(string FilterSet){


  stringstream name;
  name<<FilterSet<<"_vsEnergy";

  
  TH2F* thePlot = (TH2F*)gDirectory->Get(name.str().c_str());

  if (thePlot==NULL){
    cout<<"No plot"<<endl;
    return;
  }

  int numYBins = thePlot->GetNbinsY();

  TF1 * myFunc = new TF1("myFunc","gaus",-0.3,0.3);
  

  double * x = (double*)malloc(numYBins*sizeof(double));
  double * y = (double*)malloc(numYBins*sizeof(double));
  double * ey = (double*)malloc(numYBins*sizeof(double));
  
  for (int i=0;i<numYBins;i++){
    


    TH1D* proj=thePlot->ProjectionX("_px",i,numYBins-1);
    if (proj->GetEntries()>40){
      TFitResultPtr result = proj->Fit("myFunc","RQS");
      Int_t status = result;
      if (status==0){
	cout<<result->Value(2)*2.35*4<<endl;
	x[i]=i;
	y[i]=result->Value(2)*2.35*4;
	ey[i]=result->UpperError(2)*2.35*4;
      }
    }
  }
  
  
  TGraphErrors * graph = new TGraphErrors(numYBins,x,y,0,ey);
  graph->GetHistogram()->SetMarkerStyle(5);
  graph->Draw("AP");
}

hw()
{
    Float_t data;
    TFile *f=new TFile("data.root");
    TTree* tree = f->Get("mytree");
    TBranch* branch = tree->GetBranch("data");
    branch->SetAddress(&data);

    TH1F* h1 = new TH1F("h1","h1",200,0.0,5.0);
    for(Int_t i= 0 ; i< tree->GetEntries() ; i++) {
	tree->GetEntry(i);
        h1->Fill( data) ;
    }
       
    TF1 *f2 = new TF1("f2", "[0]*x+[1]+gaus(2)", 0.0, 5.0);
    f2->SetParameter(3,3.7);

    TVirtualFitter *min = TVirtualFitter::Fitter(0,2);
    TVirtualFitter::SetDefaultFitter("Minuit");

    TFitResultPtr fitresult = h1->Fit("f2", "MES");
    cout<< "Parameter 0" << fitresult->Parameter(0) << "+-" << fitresult->ParError(0) << endl;
    cout<<  "Parameter 1" <<fitresult->Parameter(1) << "+-" << fitresult->ParError(1) << endl;
    cout << "Chi2 " << fitresult->Chi2() << endl;
    cout << "NDF " << fitresult->Ndf() << endl;
    cout << "Chi2 probability " << fitresult->Prob() << endl;
    TMatrixD covmat(fitresult->GetCovarianceMatrix());
    cout << "Covariance matrix" << endl;
    covmat.Print();
    TCanvas *c1 = new TCanvas("c1", "binned line fit", 600, 600);
    c1->Draw();
    h1->Draw();


}

TH1F*	makehist( string id, int ndf )
{
	int	hits_count = ndf + 2;
	TString	hist_name, cut, title;
	hist_name += id;
	hist_name += "_chisq_";
	hist_name += hits_count;
	hist_name += "hits";
	cut += id;
	cut += "_hits_count == ";
	cut += hits_count;
	float	min_x = 0;
	float	max_x = (ndf == 2) ? 2.0 : 0.5;
	TH1F	*hist = new TH1F(hist_name.Data(), cut, 1000, min_x, max_x);
	events->Draw(Form("%s_chisq >> %s", id.c_str(), hist_name.Data()), cut);
	if (hits_count > 3)
	{
		TF1 *func = new TF1("func", "[0]*exp([1]*x)+[2]");
		TFitResultPtr fit = hist->Fit(func, "SQ", "", 0.0, 0.5); // S - return fit result, Q - quiet
		std::cout << hist_name << "\tslope: " << fit->Parameter(1) << "\tadd const: " << fit->Parameter(2) << std::endl;
	}
	if (id[1] == '3')
	{
		title += "Left ";
	}
	else
	{
		title += "Right ";
	}
	title += id[2];
	title += " (";
	title += cut;
	title += ")";
	hist->SetTitle(title);
	hist->SetLabelSize(0.04, "X");
	hist->SetLabelSize(0.04, "Y");
	hist->SetTitleSize(0.05, "X");
	hist->SetTitleSize(0.05, "Y");
	hist->SetTitleOffset(0.9, "X");
	hist->SetTitleOffset(1.1, "Y");
	hist->GetXaxis()->SetTitle("#chi^2, [mm]");
	hist->GetYaxis()->SetTitle("N");

	return hist;
}

	double FitConfidence::choleskyUncertainty( double xx, TFitResultPtr fitResult, TF1 * f, int nSamples ){
		int nP = f->GetNpar();

		TMatrixDSym cov = fitResult->GetCovarianceMatrix();
		double *covArray = new double[ nP * nP ]; // number of parameters x number of parameters
		covArray = cov.GetMatrixArray();

		return choleskyUncertainty( xx, covArray, f, nSamples );
	}

void estimateEfficiencyAndPurity(TH1* fraction,double cut,double& efficiency, double& purity, double& efferror, double& purerror, TFitResultPtr fitRes) {
  // efficiency loss is defined as the estimated signal below the cut over the estimated total signal
  double rangeLow, rangeHigh;
  TF1* expo = fraction->GetFunction("expo");
  expo->GetRange(rangeLow, rangeHigh);

  double signalLoss         = expo->Integral(0,cut)/fraction->GetBinWidth(1);
  double signalLossError = expo->IntegralError(0,cut,fitRes->GetParams(),fitRes->GetCovarianceMatrix().GetMatrixArray())/fraction->GetBinWidth(1);
  
  double signal_expopart         = expo->Integral(cut,fraction->GetBinLowEdge(fraction->FindBin(rangeHigh)+1))/fraction->GetBinWidth(1);
  double signal_expopartError = expo->IntegralError(cut, fraction->GetBinLowEdge(fraction->FindBin(rangeHigh)+1),fitRes->GetParams(),fitRes->GetCovarianceMatrix().GetMatrixArray())/fraction->GetBinWidth(1);

  double signal_toppoart = fraction->Integral(fraction->FindBin(rangeHigh)+1,fraction->FindBin(1.)+1);
  
  efficiency = (signal_expopart+signal_toppoart)/(signal_expopart+signal_toppoart+signalLoss);
  efferror = TMath::Sqrt( pow(signalLoss,2)*pow(signal_expopartError,2) + pow(signal_toppoart+signal_expopart,2)*pow(signalLossError,2) )/pow(signal_expopart+signal_toppoart+signalLoss,2);
  
  // purity is defined as the signal above the cut (signal_expopart+signal_toppoart) 
  // over the total data above that cut
  
  double data_kept = fraction->Integral(fraction->FindBin(cut),fraction->FindBin(1.)+1);
  double allBeforeRangeHigh = fraction->Integral(fraction->FindBin(cut),fraction->FindBin(rangeHigh));
  purity = (signal_expopart+signal_toppoart)/data_kept;
  purerror = signal_expopartError/data_kept; 
 
  if(purity > 1.){
    std::cout << "#################################################"<<std::endl;
    std::cout << "cut: "<<cut<<" purity:" << purity <<std::endl;
    std::cout << "signalLoss: " << signalLoss << " signalExtrapolated: " << signal_expopart << " signalTop: " << signal_toppoart << std::endl;
    std::cout << "all b.r.h.: " << allBeforeRangeHigh << std::endl;
    std::cout << "total signal with expo: "<< (signal_expopart+signal_toppoart)<< std::endl;
    std::cout << "total signal integrating: "<< (allBeforeRangeHigh+signal_toppoart)<< std::endl;
    std::cout << "data kept:  " << data_kept  << std::endl;
    std::cout << "ratio:      " << allBeforeRangeHigh/signal_expopart << std::endl;
    std::cout << "#################################################"<<std::endl;
  }
 
  //   std::cout << "estimateEfficiencyAndPurity for cut=" << cut << std::endl;
  //   std::cout << "signal: " << signalLoss << " " << signal_expopart << " " << signal_toppoart << std::endl;
  //   std::cout << "data kept: " << data_kept << std::endl;
  //   std::cout << "eff: " << efficiency << "+/- " << efferror << " pur: " << purity << " +/- " << purerror << std::endl;
}

Double_t dpol2Func(const TF1 *fcn, const Double_t x, const TFitResultPtr fs) {
  Double_t df[3];
  Double_t a = fcn->GetParameter(0);
  Double_t b = fcn->GetParameter(1);
  Double_t c = fcn->GetParameter(2);
  
  df[0] = 1;
  df[1] = x;
  df[2] = x*x;
  
  Double_t err2=0;
  for(Int_t i=0; i<3; i++) {
    err2 += df[i]*df[i]*(fs->GetCovarianceMatrix()[i][i]);
    for(Int_t j=i+1; j<3; j++) {
      err2 += 2.0*df[i]*df[j]*(fs->GetCovarianceMatrix()[i][j]);
    }
  }

  assert(err2>=0);
  return sqrt(err2);
}

TH1F*	make_drift_efficiency_hist(char arm, char axis, TTree *events)
{
	if ((arm != 'l') && (arm != 'r'))
	{
		throw "make_drift_efficiency_hist: Invalid arm value";
	}
	if ((axis != 'X') && (axis != 'Y'))
	{
		throw "make_drift_efficiency_hist: Invalid axis value";
	}
	char arm_chamber = (arm == 'l') ? '3' : '4';
	int canvas_cell = 1 + (int)(arm == 'l') + 2 * (int)(axis == 'X');

	TH2F *four_hit_theta_x = new TH2F(
		Form("%c%c_four_hit_theta_x", arm, axis), "",
		ANGLE_BINS, 0, ANGLE_MAX,
		X_BINS, X_MIN, X_MAX
		);
	four_hit_theta_x->SetOption("zcol");
	TH2F *any_theta_x = new TH2F(
		Form("%c%c_any_theta_x", arm, axis), "",
		ANGLE_BINS, 0, ANGLE_MAX,
		X_BINS, X_MIN, X_MAX
		);
	any_theta_x->SetOption("zcol");

	TH1F *four_hit_theta = new TH1F(
		Form("%c%c_four_hit_theta", arm, axis), "",
		ANGLE_BINS, 0, ANGLE_MAX
		);
	TH1F *any_theta = new TH1F(
		Form("%c%c_any_theta", arm, axis), "",
		ANGLE_BINS, 0, ANGLE_MAX
		);

	c1.cd();
	TH2F	*beam_profile = new TH2F("beam_profile", ";Z, [mm];Y, [mm]", 100, -20, 20, 100, -20, 20);
	beam_profile->SetOption("zcol");
	events->Draw("RL_y:RL_z >> beam_profile", "", "ZCOL");
	TF2	*xygaus = new TF2("xygaus", "xygaus");
	TFitResultPtr fit = beam_profile->Fit(xygaus, "S"); // S - return fit result 
	const double *params = fit->GetParams();
	if (!params)
	{
		throw "Error determining beam profile parameters.";
	}
	const char	*cut = Form(
		"(abs(RL_z - (%f)) < %f) && "
		"(abs(RL_y - (%f)) < %f) && "
		"(RL_x < %f) && "
		"(RL_x > %f) && "
		"(abs(LR_z - (%f)) < %f) && "
		"(abs(LR_y - (%f)) < %f) && "
		"(LR_x < %f) && "
		"(LR_x > %f)",
		params[1], params[2] * 2,
		params[3], params[4] * 2,
		X_MAX, X_MIN,
		params[1], params[2] * 2,
		params[3], params[4] * 2,
		X_MAX, X_MIN
		);

	c2.cd(canvas_cell);
	events->Draw(
		Form("RL_x:theta_%c >> %s", arm, any_theta_x->GetName()),
		cut
		);
	events->Draw(
		Form("RL_x:theta_%c >> %s", arm, four_hit_theta_x->GetName()),
		Form("(t%c%c_hits_count[0] == 4) && %s", arm_chamber, axis, cut)
		);
	TH2F *u_tx = new TH2F(
		Form("%c%c_efficiency_theta_x", arm, axis), "",
		ANGLE_BINS, 0, ANGLE_MAX,
		X_BINS, X_MIN, X_MAX
		);
	u_tx->Divide(four_hit_theta_x, any_theta_x);
	u_tx->GetXaxis()->SetTitle("\\Theta, rad");
	u_tx->GetYaxis()->SetTitle("X coordinate on target, mm");
	u_tx->SetOption("zcol");
	u_tx->Draw("zcol");

	c3.cd(canvas_cell);
	events->Draw(
		Form("theta_%c >> %s", arm, any_theta->GetName()),
		cut
		);
	events->Draw(
		Form("theta_%c >> %s", arm, four_hit_theta->GetName()),
		Form("(t%c%c_hits_count[0] == 4) && %s", arm_chamber, axis, cut)
		);
	TH1F *u_t = new TH1F(
		Form("%c%c_efficiency_theta", arm, axis), "",
		ANGLE_BINS, 0, ANGLE_MAX
		);
	u_t->Divide(four_hit_theta, any_theta);
	u_t->GetXaxis()->SetTitle("\\Theta, rad");
	u_t->GetYaxis()->SetTitle("Efficiency, 1");
	u_t->GetYaxis()->SetRangeUser(0.0, 1.0);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
            c[iCat] -> cd();
            g[iCat] -> Draw("P,same");
            c_all -> cd(iCat+1);
            g[iCat] -> Draw("P,same");
	  }
        }
        
        char funcName[150];
        sprintf(funcName,"f_prefit_%s_%d",methods.at(iMeth).c_str(),iCat);
        TF1* f_prefit;
        if( fitMethod == "pol0")
        {
          f_prefit = new TF1(funcName,"[0]",55., 230.);
          f_prefit -> SetParameter(0,1.);         
        }
        if( fitMethod == "pol1")
        {
          f_prefit = new TF1(funcName,"[0]+[1]*(x-90.)",55., 210.);
          f_prefit -> SetParameters(1.,0.00002);         
        }
        if( fitMethod == "exp3par")
        {
          f_prefit = new TF1(funcName,"1.+[0]*(1.-exp(-1.*[1]*(x-90.)))+[2]",55., 230.);
          f_prefit -> SetParameters(0.005,0.02,0.);
        }
        if( fitMethod == "exp")
        {
	  f_prefit = new TF1(funcName,"1.+[0]*(1.-exp(-1.*[1]*(x-90.)) )",55., 230.);
          f_prefit -> SetParameters(0.005,0.02);
          f_prefit -> SetParLimits(0,0.,1.);
          f_prefit -> SetParLimits(1,0.,0.05);
        }
        
        TFitResultPtr fitResult1;
        int fitStatus1 = -1;
	gClone[iCat] -> Fit(funcName,"QRHNS");
        //g[iCat] -> Fit(funcName,"QRHNS");
        fitStatus1 = fitResult1;
        
        //std::cout << " fitStatus1 = " << fitStatus1 << std::endl;
        //std::cout << " f_prefit->GetParameter(0) = " << f_prefit->GetParameter(0) << std::endl;
        //std::cout << " f_prefit->GetParameter(1) = " << f_prefit->GetParameter(1) << std::endl;
        std::cout << " f_prefit->GetChisquare()/f_prefit->GetNDF() = " << f_prefit->GetChisquare()/f_prefit->GetNDF() << std::endl;
        
        if( drawFitFunc == true && rescaleErrors == true )
        {
          for(int point = 0; point < g[iCat]->GetN(); ++point)
          {
  	    double ey = g[iCat] -> GetErrorY(point);
	    if(addSyst == false){
	      g[iCat] -> SetPointEYhigh(point,ey*sqrt(f_prefit->GetChisquare()/f_prefit->GetNDF()));
	      g[iCat] -> SetPointEYlow (point,ey*sqrt(f_prefit->GetChisquare()/f_prefit->GetNDF()));
	  
	      if(point == 4 || point == 5 || point == 6 || point == 7){
		g[iCat] -> SetPointEYhigh(point, sqrt(pow(g[iCat] -> GetErrorY(point),2)+pow(0.0005,2)) );
		g[iCat] -> SetPointEYlow(point, sqrt(pow(g[iCat] -> GetErrorY(point),2)+pow(0.0005,2)) );
		}
	  
	    }
	    //std::cout << " >>>  prima ey = " << ey  << std::endl;

	    if(SysError[iCat][point] != 0. && addSyst == true){
	      double statE2 = (pow(ey,2) - pow(0.01*SysError[iCat][point], 2));
	      //double statE2 = pow(ey,2);
	      double sysE2 = pow(SysError[iCat][point], 2);
	      //double sysE2 = pow(SysError[iCat][point], 2);

//.........这里部分代码省略.........
	
	TTree* ic5t = (TTree*)inf->Get("icPu5JetAnalyzer/t");
	TTree* t = (TTree*)inf->Get("hltanalysis/HltTree");
	
	ic5t->SetBranchAddress("jtpt",jtpt_ic5);
	ic5t->SetBranchAddress("jteta",jteta_ic5);
	ic5t->SetBranchAddress("jtphi",jtphi_ic5);
	ic5t->SetBranchAddress("nref",&njets_ic5);
	
	TTree* ak3t = (TTree*)inf->Get("akPu3PFJetAnalyzer/t");
	
	ak3t->SetBranchAddress("jtpt",jtpt_ak3);
	ak3t->SetBranchAddress("jteta",jteta_ak3);
	ak3t->SetBranchAddress("jtphi",jtphi_ak3);
	ak3t->SetBranchAddress("nref",&njets_ak3);
	
	t->SetBranchAddress("Run",&run);
	t->SetBranchAddress("Event",&evt);
	t->SetBranchAddress("hiBin",&bin);
	
	TString outname = "JetSkimTree.root";
			
	int Nevents = t->GetEntries();
	for(int iev = 0; iev < Nevents; ++iev){
		t->GetEntry(iev);
		ic5t->GetEntry(iev);
		ak3t->GetEntry(iev);
		
		for(int i = 0; i < njets_ak3; ++i){
			//if(jtpt_ak3[i] < cutPtJet) continue;
			if(fabs(jteta_ak3[i]) > cutEtaJet) continue;
			hak3pujet->Fill(jtpt_ak3[i]);
							}//1st for loop
		
		int id=0;
		for(int i = 0; i < njets_ak3; ++i){
			if(fabs(jteta_ic5[i]) < cutEtaJet)
				id =1;
			
		}// 2nd for loop
		if(id == 1){
				for(int i = 0; i < njets_ak3; ++i){
					if(jtpt_ic5[i] < cutPtJet) continue;
					//if(jtpt_ak3[i] < cutPtJet) continue;
					if(fabs(jteta_ak3[i]) > cutEtaJet) continue;
					hak3pujet_cut->Fill(jtpt_ak3[i]);
				}
			
	}// if loop
	}//evt loop
	
	
	TCanvas *c1 = new TCanvas("c1", "",46,477,700,509);
	c1->Range(-125,-11.1492,1125,3.24528);
	c1->SetBorderSize(0);
	c1->SetBorderMode(0);
	c1->SetLeftMargin(0.16);
	c1->SetTopMargin(0.06);
	c1->SetBottomMargin(0.17);
	c1->SetTicks(1,1);
	
	
	
	gROOT->LoadMacro("erf.C");
		
	TGraphAsymmErrors *dividedthingy= new TGraphAsymmErrors;
	dividedthingy->BayesDivide(hak3pujet_cut,hak3pujet);
	
	
	float xstart=20.0;
	float xend=300.0;
	
	
	TF1 *fitfcn = new TF1("fitfcn",erf,xstart,xend,2); //<============
	fitfcn->SetParameter(0,40); //<============= 
	TFitResultPtr results =dividedthingy->Fit("fitfcn","VSR");
	
	cout << "Chi2: " << results->Chi2() << " / " << results->Ndf() << " NDF" << endl;
	
	float turnon=fitfcn->GetX(0.99,10,300);
	cout << "99% at: " << turnon << endl;
	char turnontext[50]="";
	sprintf(turnontext," ak5pu at %4.1f",turnon);
	
	
	
	dividedthingy->GetFunction("fitfcn")->SetLineColor(3);
	dividedthingy->SetLineColor(1);
	dividedthingy->Draw("Ap");
	
	
	
	TText *t2 = new TText(95,0.08,turnontext);
	t2->SetTextSize(0.055);
	t2->SetTextColor(1);
	t2->Draw();
		
		
	
}