TF1 * funGausQuad=NULL; TF1 * funSig=NULL; TF1 * funQuad=NULL; TF1 * funBkg=NULL; TF1 * funDGausPol0=NULL; TF1 * funPedestal = NULL; Float_t dGausPol0Fun(Double_t *x, Double_t *par){ Double_t xx = x[0]; return par[0] + par[1]*(par[2]*TMath::Gaus(xx, 0, par[3], kTRUE)+(1-par[2])*TMath::Gaus(xx, 0, par[4], kTRUE)); } Float_t gausFun(Double_t *x, Double_t *par){ Double_t xx = x[0]; return par[0]*TMath::Gaus(xx, 0, par[1], kTRUE); } Float_t pedestalFun(Double_t *x, Double_t *par){ Double_t xx = x[0]; return par[0]+par[1]*TMath::Gaus(xx, 0, par[2], kTRUE); } Float_t quadFun(Double_t *x, Double_t *par){ if (x[0] > par[3] && x[0] < par[4]){ TF1::RejectPoint(); return 0; } return par[0] + par[1]*x[0] + par[2]*x[0]*x[0] ; //a+bx+cx**2 } void fitHistDEta(TH1F* h1, Double_t mean, Double_t res, Double_t &yield, Double_t &yieldErr){ Double_t fitRangeLow = mean - 8.*res; Double_t fitRangeHigh = mean + 8.*res; funDGausPol0 = new TF1("funDGausPol0", dGausPol0Fun, fitRangeLow, fitRangeHigh, 5); funDGausPol0->SetParameter(0, h1->GetBinContent(5)); funDGausPol0->SetParameter(1, h1->Integral(h1->FindBin(mean-3.*res), h1->FindBin(mean+3.*res))*0.6); funDGausPol0->SetParameter(2, 0.9); funDGausPol0->SetParameter(3, res); funDGausPol0->SetParameter(4, 3.*res); h1->Fit(funDGausPol0, "ER+"); Double_t f = funDGausPol0->GetParameter(2); funSig = new TF1("funSig", gausFun, fitRangeLow, fitRangeHigh, 2); funSig->SetLineStyle(2); funPedestal = new TF1("funPedestal", pedestalFun, fitRangeLow, fitRangeHigh, 3); funPedestal->SetLineStyle(3); funPedestal->SetParameter(0, funDGausPol0->GetParameter(0)); funPedestal->SetParError(0, funDGausPol0->GetParError(0)); if(f > 0.5){ funSig->SetParameters(0, f*funDGausPol0->GetParameter(1)); funSig->SetParameters(1, funDGausPol0->GetParameter(3)); funPedestal->SetParameter(1, (1.-f)*funDGausPol0->GetParameter(1)); funPedestal->SetParameter(2, funDGausPol0->GetParameter(4)); funPedestal->SetParError(1, (1.-f)*funDGausPol0->GetParError(1)); funPedestal->SetParError(2, funDGausPol0->GetParError(4)); } else{ funSig->SetParameters(0, (1.-f)*funDGausPol0->GetParameter(1)); funSig->SetParameters(1, funDGausPol0->GetParameter(4)); funPedestal->SetParameter(1, f*funDGausPol0->GetParameter(1)); funPedestal->SetParameter(2, funDGausPol0->GetParameter(3)); funPedestal->SetParError(1, f*funDGausPol0->GetParError(1)); funPedestal->SetParError(2, funDGausPol0->GetParError(3)); } // funSig->Draw("same"); Double_t sigma = funSig->GetParameter(1); //Double_t C = funDGausPol0->GetParameter(0); Double_t bkg, sig, sigFit; Double_t bkgErr, sigErr, sigFitErr; Int_t binLow = h1->FindBin(-3.*sigma); Int_t binHigh = h1->FindBin(3.*sigma); Double_t binLowEdge = h1->GetXaxis()->GetBinLowEdge(binLow); Double_t binHighEdge = h1->GetXaxis()->GetBinUpEdge(binHigh); bkg = funPedestal->Integral(binLowEdge, binHighEdge)/h1->GetBinWidth(1); bkgErr = funPedestal->IntegralError(binLowEdge, binHighEdge)/h1->GetBinWidth(1); sig = h1->IntegralAndError(binLow, binHigh, sigErr)-bkg; sigFit = funDGausPol0->Integral(binLowEdge, binHighEdge)/h1->GetBinWidth(1)-bkg; Double_t sigFitTrue = funDGausPol0->GetParameter(1)/h1->GetBinWidth(1); cout<<"signal = "<SetLineWidth(1); funGausQuad->SetLineColor(kRed); funQuad = new TF1("quad", quadFun, fitRangeLow, fitRangeHigh, 5); funQuad->SetLineWidth(2); funQuad->SetLineColor(kBlue); funQuad->SetLineStyle(2); funBkg = new TF1("bkg", "[0]+[1]*x+[2]*x*x", fitRangeLow, fitRangeHigh); funBkg->SetLineWidth(2); funBkg->SetLineColor(kRed); funBkg->SetLineStyle(2); funBkg->SetFillStyle(3); funBkg->SetFillColor(kRed); funSig = new TF1("funSig", "gaus"); funQuad->SetParameter(0, 0/*h1->GetBinContent(1)*/); if(pol0) funQuad->SetParameter(0, h1->GetBinContent(5)); funQuad->SetParameter(1, 0); funQuad->SetParameter(2, 0); if(pol0){ funQuad->FixParameter(1, 0); funQuad->FixParameter(2, 0); } funQuad->FixParameter(3, massMean-3.*massRes); funQuad->FixParameter(4, massMean+3.*massRes); Int_t fitStatus0 = h1->Fit(funQuad, "WLRIN"); if(fitStatus0 != 0) fitStatus0 = h1->Fit(funQuad, "WLMRIN"); // h1->Fit(funQuad, "RIN"); /* Double_t paramA = (h1->Integral(h1->FindBin(massMean-3.*massRes), h1->FindBin(massMean+3.*massRes)) -funQuad->GetParameter(0)*6.*massRes/h1->GetBinWidth(1)) /2.5/massRes; */ funGausQuad->SetParameter(0, h1->GetEntries()/20.); if(pol0) funGausQuad->SetParameter(0, (h1->Integral(h1->FindBin(massMean-3.*massRes), h1->FindBin(massMean+3.*massRes)) -funQuad->GetParameter(0)*6.*massRes/h1->GetBinWidth(1))*.3); funGausQuad->SetParameter(1, massMean); funGausQuad->SetParameter(2, massRes); funGausQuad->SetParameter(3, funQuad->GetParameter(0)); funGausQuad->SetParameter(4, funQuad->GetParameter(1)); funGausQuad->SetParameter(5, funQuad->GetParameter(2)); if(pol0){ funGausQuad->FixParameter(4, funQuad->GetParameter(1)); funGausQuad->FixParameter(5, funQuad->GetParameter(2)); } Int_t fitStatus1 = h1->Fit(funGausQuad, "WLRI+"); if(fitStatus1 != 0) fitStatus1 = h1->Fit(funGausQuad, "WLMRI+"); //h1->Fit(funGausQuad, "RI+"); funBkg->SetParameter(0, funGausQuad->GetParameter(3)); funBkg->SetParameter(1, funGausQuad->GetParameter(4)); funBkg->SetParameter(2, funGausQuad->GetParameter(5)); if(pol0){ funBkg->FixParameter(1, funGausQuad->GetParameter(4)); funBkg->FixParameter(2, funGausQuad->GetParameter(5)); } funSig->SetParameter(0, funGausQuad->GetParameter(0)); funSig->SetParameter(1, funGausQuad->GetParameter(1)); funSig->SetParameter(2, funGausQuad->GetParameter(2)); Double_t mean = funGausQuad->GetParameter(1); Double_t sigma =TMath::Abs(funGausQuad->GetParameter(2)); cout<< "fitted mean = "<< mean <<", fitted sigma = " << sigma <FindBin(mean-3.*sigma); Int_t binHigh = h1->FindBin(mean+3.*sigma); Double_t binLowEdge = h1->GetXaxis()->GetBinLowEdge(binLow); Double_t binHighEdge = h1->GetXaxis()->GetBinUpEdge(binHigh); bkg = funBkg->Integral(binLowEdge, binHighEdge)/h1->GetBinWidth(1); sig = h1->IntegralAndError(binLow, binHigh, sigErr)-bkg; sigFit = funSig->Integral(binLowEdge, binHighEdge)/h1->GetBinWidth(1); Double_t sigFitTrue = funGausQuad->GetParameter(0)*TMath::Sqrt(2.*TMath::Pi())*sigma/h1->GetBinWidth(1); cout<<"signal = "<