Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-06 08:14:30 

0001 #include <TNamed.h>
0002 #include <TH2.h>
0003 #include <TF1.h>
0004 #include <TFitResult.h>
0005 #include <TFitResultPtr.h>
0006 #include <TGraphErrors.h>
0007 #include <TStyle.h>
0008 #include <TProfile.h>
0009 #include <TString.h>
0010 #include <TList.h>
0011 #include <TDirectory.h>
0012 #include <TCanvas.h>
0013 #include <TLatex.h>
0014 #include <TBox.h>
0015 #include <iostream>
0016 #include <vector>
0017 #include "BinnedHistSet.h"
0018 
0019 void BinnedHistSet::printEdges(std::vector<double> edges) {
0020     int nbins = edges.size() - 1;
0021     std::cout << "Bin edges: (";
0022     for (int i=0; i<nbins; i++) std::cout << edges.at(i) << ", ";
0023     std::cout << edges.at(nbins) << ")\n";
0024 }
0025 
0026 std::vector<double> BinnedHistSet::getUniformBinEdges(int nbins, double lower, double upper) {
0027     double range = upper - lower;
0028     double delta = range/nbins;
0029     std::vector<double> edges;
0030     for (int i=0; i<nbins; i++) {
0031     double ledge = lower + i*delta;
0032     edges.push_back(ledge);
0033     }
0034     edges.push_back(upper);
0035     return edges;
0036 }
0037 
0038 std::vector<double> BinnedHistSet::getUniformBinEdges(int nbins, TH1* hist) {
0039     double lower = hist->GetXaxis()->GetXmin();
0040     double upper = hist->GetXaxis()->GetXmax();
0041     return getUniformBinEdges(nbins, lower, upper);
0042 }
0043 
0044 std::vector<double> BinnedHistSet::getEquallyPopulatedBinEdges(int nbins, TH1* hist) {
0045     std::vector<double> edges;
0046     edges.push_back(hist->GetXaxis()->GetXmin());
0047     double* quants = new double[nbins]; // positions in [0,1] at which to compute the quantiles
0048     for (int i=0; i<nbins; i++) {
0049     quants[i] = (double)(i+1)/nbins;
0050     }
0051     double* x_q = new double[nbins]; // array to hold the computed quantiles
0052     hist->GetQuantiles(nbins, x_q, quants);
0053     for (int i=0; i<nbins; i++) {
0054     edges.push_back(x_q[i]);
0055     }
0056     return edges;
0057 }
0058 
0059 std::vector<double> BinnedHistSet::getBinCenters(std::vector<double> edges) {
0060     std::vector<double> centers;
0061     int nbins = edges.size() - 1;
0062     for (int i=0; i<nbins; i++) {
0063     double center = (edges.at(i) + edges.at(i+1))/2.0;
0064     centers.push_back(center);
0065     }
0066     return centers;
0067 }
0068 
0069 double fitfunc(double* x, double* par) {
0070     double arg = (x[0] - par[1])/par[2];
0071     return (par[0]*TMath::Exp(-0.5*arg*arg) + par[3]);
0072 }
0073 
0074 std::vector<double> BinnedHistSet::doubleFitHist(TH1* hist, bool response) {
0075     std::vector<double> out;
0076     double mean;
0077     double stddev;
0078     /* TF1* func = new TF1("fit", fitfunc, -2.0, 2.0, 4); */
0079     /* func->SetParameters(100, hist->GetMean(), hist->GetRMS(), 10); */
0080     /* func->SetParLimits(2, 0, 999); */
0081     /* hist->Fit(func, "SQ"); */
0082     /* TF1* fit1 = hist->GetFunction("fit"); */
0083     if (response) hist->Fit("gaus", "SQ", "", 0.8, 1.3);
0084     else hist->Fit("gaus", "SQ");
0085     TF1* fit1 = hist->GetFunction("gaus");
0086     if (!fit1) {
0087     std::cout << "Greg info: bad fit #1 for hist ";
0088     hist->Print();
0089     out.push_back(9999);
0090     out.push_back(0);
0091     out.push_back(1);
0092     return out;
0093     }
0094     mean = fit1->GetParameter(1);
0095     stddev = fit1->GetParameter(2);
0096     hist->GetListOfFunctions()->Clear();
0097     /* hist->Fit(fit1, "SQ", "", (mean-1.5*stddev), (mean+1.5*stddev)); */
0098     /* TF1* fit2 = hist->GetFunction("fit"); */
0099     hist->Fit("gaus", "SQ", "", (mean-1.0*stddev), (mean+1.0*stddev));
0100     TF1* fit2 = hist->GetFunction("gaus");
0101     if (!fit2) {
0102     std::cout << "Greg info: bad fit #2 for hist ";
0103     hist->Print();
0104     out.push_back(9999);
0105     out.push_back(0);
0106     out.push_back(1);
0107     return out;
0108     }
0109     mean = fit2->GetParameter(1);
0110     stddev = fit2->GetParameter(2);
0111     hist->GetListOfFunctions()->Clear();
0112     hist->Fit("gaus", "SQ", "", (mean-1.0*stddev), (mean+1.0*stddev));
0113     TF1* fit3 = hist->GetFunction("gaus");
0114     if (!fit3) {
0115     std::cout << "Greg info: bad fit #3 for hist ";
0116     hist->Print();
0117     out.push_back(9999);
0118     out.push_back(0);
0119     out.push_back(1);
0120     return out;
0121     }
0122     mean = fit3->GetParameter(1);
0123     stddev = fit3->GetParameter(2);
0124     out.push_back(mean);
0125     out.push_back(stddev);
0126     out.push_back(hist->GetEntries());
0127     // next bit is only relevant for calculating error on the energy resolution
0128     double resolution_error = (fit3->GetParError(1))*(fit3->GetParError(1));
0129     resolution_error += (fit3->GetParError(2))*(fit3->GetParError(2));
0130     resolution_error = sqrt(resolution_error);
0131     out.push_back(resolution_error);
0132     return out;
0133 }
0134 
0135 BinnedHistSet::BinnedHistSet() {
0136     name_base = "EMPTY";
0137     title_base = "EMPTY";
0138     binned_quantity = "EMPTY";
0139 }
0140 
0141 /* BinnedHistSet::BinnedHistSet(std::string name, std::string title, int n_hist_bins, double lower, double upper, std::string binquantity, std::vector<double> edges): */ 
0142 /*     TNamed(name, title) { */
0143 BinnedHistSet::BinnedHistSet(std::string name, std::string title, int n_hist_bins, double lower, double upper, std::string binquantity, std::vector<double> edges) {
0144     name_base = name;
0145     title_base = title;
0146     hist_Nbins = n_hist_bins;
0147     hist_lower = lower;
0148     hist_upper = upper;
0149     binned_quantity = binquantity;
0150     bin_edges = edges;
0151     nbins = edges.size() - 1;
0152     /* std::cout << "In initializer: name_base=" << name_base << ", nbins=" << nbins << "\n"; */
0153 }
0154 
0155 /* BinnedHistSet::BinnedHistSet(TH1* h_prototype, std::string binquantity, std::vector<double> edges): */ 
0156     /* TNamed(Form("bhs_%s", h_prototype->GetName()), h_prototype->GetTitle()) { */
0157 BinnedHistSet::BinnedHistSet(TH1* h_prototype, std::string binquantity, std::vector<double> edges) {
0158     /* std::cout << "In initializer: h_prototype->GetName() = " << h_prototype->GetName() << "\n"; */
0159     name_base = std::string(h_prototype->GetName()) + "__" + binquantity;
0160     /* name_base.replace(name_base.find("h_"), std::string("h_").size(), "bhs_"); */
0161     /* std::cout << "In initializer: h_prototype->GetTitle() = " << h_prototype->GetTitle() << "\n"; */
0162     title_base = h_prototype->GetTitle();
0163     std::string xtitle = h_prototype->GetXaxis()->GetTitle();
0164     std::string ytitle = h_prototype->GetYaxis()->GetTitle();
0165     axislabels = ";" + xtitle + ";" + ytitle;
0166     hist_Nbins = h_prototype->GetNbinsX();
0167     hist_lower = h_prototype->GetXaxis()->GetXmin();
0168     hist_upper = h_prototype->GetXaxis()->GetXmax();;
0169     binned_quantity = binquantity;
0170     bin_edges = edges;
0171     nbins = edges.size() - 1;
0172 }
0173 
0174 BinnedHistSet::~BinnedHistSet() {
0175     std::cout << "Entering destructor for BinnedHitSet " << name_base << "\n";
0176 }
0177 
0178 void BinnedHistSet::GetAxisLabels() {
0179     // if we used a prototype histogram, axislabel should already be filled in
0180     if (axislabels != "") return;
0181     size_t pos = title_base.find(";");
0182     /* std::cout << "Getting axis labels\n"; */
0183     /* std::cout << "title_base = " << title_base << "\n"; */
0184     /* std::cout << "pos = " << pos << "\n"; */
0185     std::string axis_title = title_base.substr(pos);
0186     /* std::cout << "Got axis labels ( = " << axis_title << " )\n"; */
0187     axislabels = axis_title;
0188 }
0189 
0190 std::string BinnedHistSet::GetTitlePrefix() {
0191     size_t pos = title_base.find(";");
0192     /* std::cout << "Getting title prefix\n"; */
0193     std::string prefix = title_base.substr(0, pos);
0194     /* std::cout << "Got title prefix ( = " << prefix << " )\n"; */
0195     return prefix;
0196 }
0197 
0198 std::string BinnedHistSet::GetBinnedQuantityName() {
0199     size_t pos = binned_quantity.find_first_of("([");
0200     /* std::cout << "Getting binned quantity name\n"; */
0201     std::string name = binned_quantity.substr(0, pos);
0202     /* std::cout << "Got binned quantity name ( = " << name << " )\n"; */
0203     return name;
0204 }
0205 
0206 std::string BinnedHistSet::GetBinnedQuantityUnits() {
0207     size_t start_pos = binned_quantity.find_first_of("([");
0208     size_t end_pos = binned_quantity.find_first_of(")]");
0209     int length = end_pos - start_pos - 1;
0210     std::string units;
0211     /* std::cout << "Getting binned quantity units\n"; */
0212     if (start_pos == end_pos) units = "";
0213     else units = binned_quantity.substr(start_pos+1, length);
0214     /* std::cout << "Got binned quantity units ( = " << units << " )\n"; */
0215     return units;
0216 }
0217 
0218 void BinnedHistSet::MakeHists() {
0219     /* std::cout << "Entering MakeHists for " << name_base << "; nbins=" << nbins << "\n"; */
0220     std::string title_prefix = this->GetTitlePrefix(); // Hold the first part of the title string
0221     /* std::cout << "Got title_prefix = " << title_prefix << std::endl; */
0222     this->GetAxisLabels(); // Hold the axis label part of the title string
0223     /* std::cout << "Got axislabels = " << axislabels << std::endl; */
0224     char* title_range; // Hold the range of bin_quantity for the current bin */
0225     std::string binquant = this->GetBinnedQuantityName(); // Hold the name for bin_quantity
0226     /* std::cout << "Got binquant = " << binquant << std::endl; */
0227     std::string units = this->GetBinnedQuantityUnits(); // Hold the units for bin_quantity
0228     /* std::cout << "Got units = " << units << std::endl; */
0229 
0230     TH1* h;
0231     /* std::cout << "Making hists; hist0 has name " << Form("%s_0", name_base.c_str()) << " and title " << title_base.c_str() << "\n"; */
0232     h = new TH1D(Form("%s_0", name_base.c_str()), title_base.c_str(), hist_Nbins, hist_lower, hist_upper);
0233     hist_vec.push_back(h);
0234     for (int i=0; i<nbins; i++) {
0235     char* namestring = Form("%s_%d", name_base.c_str(), i+1);
0236     double lower = bin_edges.at(i);
0237     double upper = bin_edges.at(i+1);
0238     title_range = Form("%.2f%s #leq %s < %.2f%s", lower, units.c_str(), binquant.c_str(), upper, units.c_str());
0239     char* titlestring = Form("%s [%s]%s", title_prefix.c_str(), title_range, axislabels.c_str());
0240     /* std::cout << "Making hists; hist" << (i+1) << " has name " << namestring << " and title " << titlestring << "\n"; */
0241     h = new TH1D(namestring, titlestring, hist_Nbins, hist_lower, hist_upper);
0242     hist_vec.push_back(h);
0243     }
0244     /* hist_vec.at(0)->Print(); */
0245 }
0246 
0247 void BinnedHistSet::FillHists(double binValue, double fillValue) {
0248     if (hist_vec.size() < 1) { // don't try to fill anything if hist_arr isn't set
0249     std::cout << "Greg warning: trying to FillHists without setting hist_vec. Skipping this!\n";
0250     return;
0251     }
0252     // First fill the unbinned distribution
0253     hist_vec.at(0)->Fill(fillValue);
0254     /* std::cout << "Entering BinnedHistSet::fillHists\n"; */
0255     /* printEdges(edges); */
0256     /* int nbins = edges.size() - 1; */
0257     // Find which bin the event falls into
0258     int i; // index of the histogram we will fill
0259     for (i=0; i<nbins; i++) {
0260     double ledge = bin_edges.at(i);
0261     if ( (i==0) && (binValue<ledge) ) {
0262         i = -1;  // binValue is below the lowest bin edge; stop with i=-1
0263         break;
0264     }
0265     double redge = bin_edges.at(i+1);
0266     if ( (binValue>=ledge) && (binValue<redge) ) break;  // we have the correct value for i; stop here
0267     if ( (i==(nbins-1)) && (binValue>redge) ) {
0268         i = nbins;  // binValue is above the highest bin edge; stop with i=nbins
0269         break;
0270     }
0271     }
0272     /* std::cout << "Ended up with i=" << i << "\n"; */
0273     if ( (i==-1) || (i==nbins) ) return;  // binValue was outside the expected range; don't fill anything
0274     // Fill the corresponding histogram
0275     hist_vec.at(i+1)->Fill(fillValue);
0276     /* hist_vec.at(i+1)->Print(); */
0277     return;
0278 }
0279 
0280 void BinnedHistSet::WriteHists() {
0281     for (int i=0; i<nbins+1; i++) {
0282     /* std::cout << "hist_vec.at(" << i << ") = " << hist_vec.at(i) << "\nPrinting: "; */
0283     /* hist_vec.at(i)->Print(); */
0284     hist_vec.at(i)->Write();
0285     }
0286     std::cout << "Wrote hists for " << name_base << "\n";
0287 }
0288 
0289 void BinnedHistSet::GetHistsFromFile(std::string name) {
0290     /* std::cout << "In GetHistsFromFile for " << name_base << ": nbins=" << nbins << "\n"; */
0291     for (int i=0; i<nbins+1; i++) {
0292     char* namestring = Form("%s_%d", name.c_str(), i);
0293     /* std::cout << "In GetHistsFromFile for " << name_base << ": Looking for hist named " << namestring << "\n"; */
0294     TH1D* hist = (TH1D*)(gDirectory->Get(namestring));
0295     hist_vec.push_back(hist);
0296     /* std::cout << "Address is " << hist << "; title is " << hist->GetTitle() << "\n"; */
0297     }
0298 }
0299 
0300 TH1** BinnedHistSet::GetHistArr() {
0301     TH1** hist_arr = new TH1*[nbins+1];
0302     for (int i=0; i<nbins+1; i++) {
0303     hist_arr[i] = hist_vec.at(i);
0304     }
0305     return hist_arr;
0306 }
0307 
0308 void BinnedHistSet::PlotAllHistsWithFits(std::string outpdfname, bool response=false, std::string fitfunc="gaus") {
0309     /* std::cout << "In PlotAllHistsWithFits for " << name_base << "; nbins=" << nbins << "\n"; */
0310     TCanvas *c1 = new TCanvas("c1", "c1",0,50,1600,900);
0311     TH1* hist;
0312     for (int i=0; i<=nbins; i++) {
0313     hist = hist_vec.at(i);
0314     /* std::cout << "Plotting hist " << hist->GetName() << "\n"; */
0315     TLatex* fit_latex = new TLatex(0.75, 0.5, "blaahhh");
0316     fit_latex->SetTextSize(fit_latex->GetTextSize()/1.25);
0317     fit_latex->SetTextColor(kRed);
0318     fit_latex->SetTextAlign(21);
0319     double mean;
0320     double stddev;
0321     /* TF1* func = new TF1("fit", fitfunc, -2.0, 2.0, 4); */
0322     /* func->SetParameters(100, hist->GetMean(), hist->GetRMS(), 10); */
0323     /* func->SetParLimits(2, 0, 999); */
0324     c1->cd();
0325     hist->Draw("e1 x0");
0326     /* hist->Fit(func, "SQ"); */
0327     /* TF1* fit1 = hist->GetFunction("fit"); */
0328     if (fitfunc == "") {
0329         c1->Modified();
0330         c1->SaveAs(outpdfname.c_str());
0331     }
0332     if (fitfunc == "gaus") {
0333         if (response) hist->Fit("gaus", "SQ", "", 0.8, 1.3);
0334         else hist->Fit("gaus", "SQ");
0335         TF1* fit1 = hist->GetFunction("gaus");
0336         if (!fit1) {
0337         std::cout << "\n\nGreg info: bad fit #1 for hist ";
0338         hist->Print();
0339         fit_latex->DrawLatexNDC(0.50, 0.80, "First fit failed!");
0340         c1->Modified();
0341         c1->SaveAs(outpdfname.c_str());
0342         continue;
0343         }
0344         mean = fit1->GetParameter(1);
0345         stddev = fit1->GetParameter(2);
0346         hist->GetListOfFunctions()->Clear();
0347         hist->Fit("gaus", "SQ", "", (mean-stddev), (mean+stddev));
0348         /* hist->Fit(fit1, "SQ", "", (mean-stddev), (mean+stddev)); */
0349         /* TF1* fit2 = hist->GetFunction("fit"); */
0350         TF1* fit2 = hist->GetFunction("gaus");
0351         if (!fit2) {
0352         std::cout << "\n\nGreg info: bad fit #2 for hist ";
0353         hist->Print();
0354         fit_latex->DrawLatexNDC(0.50, 0.80, "Second fit failed! Showing first fit.");
0355         fit1->Draw("same");
0356         c1->Modified();
0357         c1->SaveAs(outpdfname.c_str());
0358         continue;
0359         }
0360         mean = fit2->GetParameter(1);
0361         stddev = fit2->GetParameter(2);
0362         hist->GetListOfFunctions()->Clear();
0363         hist->Fit("gaus", "SQ", "", (mean-stddev), (mean+stddev));
0364         TF1* fit3 = hist->GetFunction("gaus");
0365         if (!fit3) {
0366         std::cout << "\n\nGreg info: bad fit #3 for hist ";
0367         hist->Print();
0368         fit_latex->DrawLatexNDC(0.50, 0.80, "Third fit failed! Showing second fit.");
0369         fit2->Draw("same");
0370         c1->Modified();
0371         c1->SaveAs(outpdfname.c_str());
0372         continue;
0373         }
0374         c1->Modified();
0375         c1->SaveAs(outpdfname.c_str());
0376         delete fit_latex;
0377     }
0378     if (fitfunc == "mass") {
0379         TF1* fit = new TF1("fit", "[0] + [1]*sqrt(x) + [2]*x + [3]*x^2 + [4]*x^3 + [5]*exp(-0.5*((x-[6])/[7])^2)", 0.0, 1.0);
0380         fit->SetParLimits(5, hist->GetMaximum()/20.0, hist->GetMaximum()*1.2);
0381         fit->SetParLimits(6, 0.05, 0.35);
0382         fit->SetParLimits(7, 0.01, 0.10);
0383         fit->SetParameter(5, hist->GetMaximum()/2.0);
0384         fit->SetParameter(6, 0.145);
0385         fit->SetParameter(7, 0.025);
0386 
0387         hist->Fit(fit, "SQ", "E1", 0.05, 0.5);
0388         mean = fit->GetParameter(6);
0389         stddev = fit->GetParameter(7);
0390         fit->SetLineColor(kRed);
0391         TF1* signal = new TF1("signal", "[0]*exp(-0.5*((x-[1])/[2])^2)", 0.05, 0.5);
0392         TF1* background = new TF1("background", "[0] + [1]*sqrt(x) + [2]*x + [3]*x^2 + [4]*x^3", 0.05, 0.5);
0393         for (int i=0; i<5; i++) {
0394         background->SetParameter(i, fit->GetParameter(i));
0395         if (i<3) signal->SetParameter(i, fit->GetParameter(i+5));
0396         }
0397         signal->SetLineColor(kBlue);
0398         background->SetLineColor(kGreen);
0399         signal->Draw("same");
0400         background->SetLineStyle(kDashed);
0401         background->Draw("same");
0402         /* fit->Draw("same"); */
0403 
0404         TLatex* fit_latex = new TLatex(0.25, 0.8, "blaahhh");
0405         /* fit_latex->SetTextSize(fit_latex->GetTextSize()/1.25); */
0406         fit_latex->SetTextColor(kBlue);
0407         fit_latex->SetTextAlign(21);
0408         fit_latex->DrawLatexNDC(0.45, 0.70, Form("#splitline{#mu=%.5f#pm%.5f}{#sigma=%.5f#pm%.5f}", fit->GetParameter(6), fit->GetParError(6), fit->GetParameter(7), fit->GetParError(7)));
0409         c1->Modified();
0410         c1->SaveAs(outpdfname.c_str());
0411         hist->GetListOfFunctions()->Clear();
0412         delete fit_latex;
0413         delete fit;
0414     }
0415     if (fitfunc == "mass2") {
0416         double pi0range_low = 0.08;
0417         double pi0range_high = 0.22;
0418         double etarange_low = 0.45;
0419         double etarange_high = 0.75;
0420         double pi0bkgr_low_low = 0.02;
0421         double pi0bkgr_low_high = 0.07;
0422         double pi0bkgr_high_low = 0.23;
0423         double pi0bkgr_high_high = 0.28;
0424         double etabkgr_low_low = 0.33;
0425         double etabkgr_low_high = 0.43;
0426         double etabkgr_high_low = 0.77;
0427         double etabkgr_high_high = 0.87;
0428 
0429         // make a dummy hist to fit the eta background, excluding the signal region
0430         TH1* htemp = (TH1*)hist->Clone("htemp");
0431         for (int j=1; j<=htemp->GetNbinsX(); j++) {
0432         float bincenter = htemp->GetBinCenter(j);
0433         if (bincenter > etarange_low && bincenter < etarange_high) {
0434             htemp->SetBinError(j, 1.0e18);
0435         }
0436         }
0437         TF1* fit_etabkgr = new TF1("etabkgr", "[0] + [1]*x + [2]*x^2", 0.30, 1.0);
0438         TFitResultPtr fit_etabkgr_res = htemp->Fit(fit_etabkgr, "SQ", "E1", 0.30, 0.9);
0439         /* std::cout << "fit_etabkgr_res = " << fit_etabkgr_res << std::endl; */
0440 
0441         // make a dummy hist to fit the eta signal, after subtracting the background
0442         TH1* htemp2 = (TH1*)hist->Clone("htemp2");
0443         for (int j=1; j<=htemp2->GetNbinsX(); j++) {
0444         float bincenter = htemp2->GetBinCenter(j);
0445         if (bincenter > 0.40 && bincenter < 0.80) {
0446             float signal = htemp2->GetBinContent(j) - fit_etabkgr->Eval(bincenter);
0447             htemp2->SetBinContent(j, signal);
0448             htemp2->SetBinError(j, sqrt(signal));
0449         }
0450         }
0451         TF1* fit_etasignal = new TF1("fit_etasignal", "[0]*exp(-0.5*((x-[1])/[2])^2)", 0.30, 1.0);
0452         fit_etasignal->SetParLimits(1, 0.40, 0.80);
0453         fit_etasignal->SetParLimits(2, 0.01, 0.25);
0454         fit_etasignal->SetParameter(1, 0.60);
0455         fit_etasignal->SetParameter(2, 0.05);
0456         TFitResultPtr fit_etasignal_res = htemp2->Fit(fit_etasignal, "SQ", "E1", 0.40, 0.80);
0457         /* std::cout << "fit_etasignal_res = " << fit_etasignal_res << std::endl; */
0458         TF1* fit_eta = new TF1("fit_eta", "[0] + [1]*x + [2]*x^2 + [3]*exp(-0.5*((x-[4])/[5])^2)", 0.30, 1.0);
0459         fit_eta->SetParameter(0, fit_etabkgr->GetParameter(0));
0460         fit_eta->SetParameter(1, fit_etabkgr->GetParameter(1));
0461         fit_eta->SetParameter(2, fit_etabkgr->GetParameter(2));
0462         fit_eta->SetParameter(3, fit_etasignal->GetParameter(0));
0463         fit_eta->SetParameter(4, fit_etasignal->GetParameter(1));
0464         fit_eta->SetParameter(5, fit_etasignal->GetParameter(2));
0465 
0466         // make a dummy hist to fit the pi0 background, excluding the signal region
0467         htemp = (TH1*)hist->Clone("htemp");
0468         for (int j=1; j<=htemp->GetNbinsX(); j++) {
0469         float bincenter = htemp->GetBinCenter(j);
0470         if (bincenter > pi0range_low && bincenter < pi0range_high) {
0471             htemp->SetBinError(j, 1.0e18);
0472         }
0473         }
0474         TF1* fit_pi0bkgr = new TF1("pi0bkgr", "[0] + [1]*x + [2]*x^2 + [3]*x^3", 0.02, 0.30);
0475         TFitResultPtr fit_pi0bkgr_res = htemp->Fit(fit_pi0bkgr, "SQ", "E1", 0.02, 0.30);
0476         /* std::cout << "fit_pi0bkgr_res = " << fit_pi0bkgr_res << std::endl; */
0477 
0478         // make a dummy hist to fit the pi0 signal, after subtracting the background
0479         htemp2 = (TH1*)hist->Clone("htemp2");
0480         for (int j=1; j<=htemp2->GetNbinsX(); j++) {
0481         float bincenter = htemp2->GetBinCenter(j);
0482         if (bincenter > pi0range_low && bincenter < pi0range_high) {
0483             float signal = htemp2->GetBinContent(j) - fit_pi0bkgr->Eval(bincenter);
0484             htemp2->SetBinContent(j, signal);
0485             htemp2->SetBinError(j, sqrt(signal));
0486         }
0487         }
0488         TF1* fit_pi0signal = new TF1("fit_pi0signal", "[0]*exp(-0.5*((x-[1])/[2])^2)", 0.02, 0.30);
0489         fit_pi0signal->SetParLimits(1, 0.05, 0.25);
0490         fit_pi0signal->SetParLimits(2, 0.01, 0.25);
0491         fit_pi0signal->SetParameter(1, 0.15);
0492         fit_pi0signal->SetParameter(2, 0.05);
0493         TFitResultPtr fit_pi0signal_res = htemp2->Fit(fit_pi0signal, "SQ", "E1", pi0range_low, pi0range_high);
0494         /* std::cout << "fit_pi0signal_res = " << fit_pi0signal_res << std::endl; */
0495         TF1* fit_pi0 = new TF1("fit_pi0", "[0] + [1]*x + [2]*x^2 + [3]*x^3 + [4]*exp(-0.5*((x-[5])/[6])^2)", 0.02, 0.30);
0496         fit_pi0->SetParameter(0, fit_pi0bkgr->GetParameter(0));
0497         fit_pi0->SetParameter(1, fit_pi0bkgr->GetParameter(1));
0498         fit_pi0->SetParameter(2, fit_pi0bkgr->GetParameter(2));
0499         fit_pi0->SetParameter(3, fit_pi0bkgr->GetParameter(3));
0500         fit_pi0->SetParameter(4, fit_pi0signal->GetParameter(0));
0501         fit_pi0->SetParameter(5, fit_pi0signal->GetParameter(1));
0502         fit_pi0->SetParameter(6, fit_pi0signal->GetParameter(2));
0503 
0504         /* // do the pi0 fit last so hist is the most recent thing drawn */
0505         /* TF1* fit_pi0 = new TF1("fit_pi0", "[0] + [1]*x + [2]*x^2 + [3]*x^3 + [4]*exp(-0.5*((x-[5])/[6])^2)", 0.0, 1.0); */
0506         /* fit_pi0->SetParLimits(4, hist->GetMaximum()/20.0, hist->GetMaximum()*1.2); */
0507         /* fit_pi0->SetParLimits(5, 0.05, 0.35); */
0508         /* fit_pi0->SetParLimits(6, 0.01, 0.15); */
0509         /* fit_pi0->SetParameter(4, hist->GetMaximum()/2.0); */
0510         /* fit_pi0->SetParameter(5, 0.150); */
0511         /* fit_pi0->SetParameter(6, 0.025); */
0512 
0513         /* /1* hist->Fit(fit_pi0, "SQ", "E1", 0.05, 0.45); *1/ */
0514         /* TFitResultPtr fit_pi0_res = hist->Fit(fit_pi0, "SQ", "E1", 0.02, 0.29); */
0515         /* /1* gStyle->SetOptFit(1100); *1/ */
0516         /* gStyle->SetOptFit(0); */
0517         /* mean = fit_pi0->GetParameter(5); */
0518         /* stddev = fit_pi0->GetParameter(6); */
0519         /* fit_pi0->SetLineColor(kRed); */
0520         /* TF1* pi0signal = new TF1("pi0signal", "[0]*exp(-0.5*((x-[1])/[2])^2)", 0.02, 0.29); */
0521         /* TF1* pi0background = new TF1("pi0background", "[0] + [1]*x + [2]*x^2 + [3]*x^3", 0.02, 0.29); */
0522         /* for (int i=0; i<4; i++) { */
0523         /* pi0background->SetParameter(i, fit_pi0->GetParameter(i)); */
0524         /* pi0background->SetParError(i, fit_pi0->GetParError(i)); */
0525         /* if (i<3) { */
0526             /* pi0signal->SetParameter(i, fit_pi0->GetParameter(i+4)); */
0527             /* pi0signal->SetParError(i, fit_pi0->GetParError(i+4)); */
0528         /* } */
0529         /* } */
0530         /* pi0signal->SetLineColor(kBlue); */
0531         /* pi0background->SetLineColor(kGreen); */
0532         /* pi0background->SetLineStyle(kDashed); */
0533         hist->Draw("x0 e1");
0534         c1->Modified();
0535         c1->Update();
0536         double ymin = c1->GetUymin();
0537         double ymax = c1->GetUymax();
0538         hist->Print();
0539         std::cout << Form("ymin = %f, ymax = %f\n", ymin, ymax);
0540         fit_eta->SetLineColor(kRed);
0541         fit_etasignal->SetLineColor(kBlue);
0542         fit_etabkgr->SetLineColor(kGreen);
0543         fit_etabkgr->SetLineStyle(kDashed);
0544         fit_eta->Draw("same");
0545         fit_etasignal->Draw("same");
0546         fit_etabkgr->Draw("same");
0547         /* pi0signal->Draw("same"); */
0548         /* pi0background->Draw("same"); */
0549         /* fit->Draw("same"); */
0550         fit_pi0->SetLineColor(kRed);
0551         fit_pi0signal->SetLineColor(kBlue);
0552         fit_pi0bkgr->SetLineColor(kGreen);
0553         fit_pi0bkgr->SetLineStyle(kDashed);
0554         fit_pi0->Draw("same");
0555         fit_pi0signal->Draw("same");
0556         fit_pi0bkgr->Draw("same");
0557 
0558         TLatex* fit_pi0_latex = new TLatex(0.25, 0.8, "blaahhh");
0559         /* fit_pi0_latex->SetTextSize(fit_pi0_latex->GetTextSize()/1.25); */
0560         fit_pi0_latex->SetTextColor(kBlue);
0561         fit_pi0_latex->SetTextAlign(21);
0562         fit_pi0_latex->DrawLatexNDC(0.35, 0.80, Form("#splitline{#mu=%.5f#pm%.5f}{#sigma=%.5f#pm%.5f}", fit_pi0->GetParameter(5), fit_pi0->GetParError(5), fit_pi0->GetParameter(6), fit_pi0->GetParError(6)));
0563 
0564         TLatex* fit_eta_latex = new TLatex(0.25, 0.8, "blaahhh");
0565         /* fit_eta_latex->SetTextSize(fit_eta_latex->GetTextSize()/1.25); */
0566         fit_eta_latex->SetTextColor(kBlue);
0567         fit_eta_latex->SetTextAlign(21);
0568         fit_eta_latex->DrawLatexNDC(0.65, 0.50, Form("#splitline{#mu=%.4f#pm%.4f}{#sigma=%.4f#pm%.4f}", fit_etasignal->GetParameter(1), fit_etasignal->GetParError(1), fit_etasignal->GetParameter(2), fit_etasignal->GetParError(2)));
0569 
0570 
0571         double pi0signal_int = fit_pi0signal->Integral(pi0range_low, pi0range_high);
0572         double pi0signal_int_err = 0.0;
0573         if (fit_pi0signal_res != -1) pi0signal_int_err = fit_pi0signal->IntegralError(pi0range_low, pi0range_high, fit_pi0signal_res->GetParams(), fit_pi0signal_res->GetCovarianceMatrix().GetMatrixArray());
0574         double pi0background_int = fit_pi0bkgr->Integral(pi0range_low, pi0range_high);
0575         double pi0background_int_err = 0.0;
0576         if (fit_pi0bkgr_res != -1) pi0background_int_err = fit_pi0bkgr->IntegralError(pi0range_low, pi0range_high, fit_pi0bkgr_res->GetParams(), fit_pi0bkgr_res->GetCovarianceMatrix().GetMatrixArray());
0577         double etasignal_int = fit_etasignal->Integral(etarange_low, etarange_high);
0578         double etasignal_int_err = 0.0;
0579         if (fit_etasignal_res != -1) etasignal_int_err = fit_etasignal->IntegralError(etarange_low, etarange_high, fit_etasignal_res->GetParams(), fit_etasignal_res->GetCovarianceMatrix().GetMatrixArray());
0580         double etabackground_int = fit_etabkgr->Integral(etarange_low, etarange_high);
0581         double etabackground_int_err = 0.0;
0582         if (fit_etabkgr_res != -1) etabackground_int_err = fit_etabkgr->IntegralError(etarange_low, etarange_high, fit_etabkgr_res->GetParams(), fit_etabkgr_res->GetCovarianceMatrix().GetMatrixArray());
0583         double pi0_bkgr_fraction = pi0background_int / (pi0signal_int+pi0background_int);
0584         double pi0_sum = pi0signal_int + pi0background_int;
0585         double pi0_err1 = pi0background_int*pi0signal_int_err;
0586         double pi0_err2 = pi0background_int_err*pi0signal_int;
0587         double pi0_bkgr_fraction_err = sqrt((pi0_err1*pi0_err1) + (pi0_err2*pi0_err2))/(pi0_sum*pi0_sum);
0588         /* std::cout << Form("Greg info: pi0background_int = %f, pi0background_int_err = %f\n", pi0background_int, pi0background_int_err); */
0589         /* std::cout << Form("Greg info: pi0signal_int = %f, pi0signal_int_err = %f\n", pi0signal_int, pi0signal_int_err); */
0590         /* std::cout << Form("Greg info: dr = 1/(%f)^2 * sqrt( (%f)^2 + (%f)^2 )\n", pi0_sum, pi0_err1, pi0_err2); */
0591         double eta_bkgr_fraction = etabackground_int / (etasignal_int+etabackground_int);
0592         double eta_sum = etasignal_int + etabackground_int;
0593         double eta_err1 = etabackground_int*etasignal_int_err;
0594         double eta_err2 = etabackground_int_err*etasignal_int;
0595         double eta_bkgr_fraction_err = sqrt((eta_err1*eta_err1) + (eta_err2*eta_err2))/(eta_sum*eta_sum);
0596         /* std::cout << Form("Greg info: etabackground_int = %f, etabackground_int_err = %f\n", etabackground_int, etabackground_int_err); */
0597         /* std::cout << Form("Greg info: etasignal_int = %f, etasignal_int_err = %f\n", etasignal_int, etasignal_int_err); */
0598         /* double eta_bkgr_fraction = etabackground_int / (etasignal_int+etabackground_int); */
0599         
0600         fit_pi0_latex->DrawLatexNDC(0.35, 0.65, Form("r = %.5f#pm%0.5f", pi0_bkgr_fraction, pi0_bkgr_fraction_err));
0601         fit_eta_latex->DrawLatexNDC(0.65, 0.35, Form("r = %.5f#pm%0.5f", eta_bkgr_fraction, eta_bkgr_fraction_err));
0602         
0603         TBox* box_pi0 = new TBox(pi0range_low, ymin, pi0range_high, ymax);
0604         box_pi0->SetLineColorAlpha(2, 0.5);
0605         box_pi0->SetLineStyle(kDotted);
0606         box_pi0->SetFillColorAlpha(2, 0.35);
0607         box_pi0->Draw("same");
0608         TBox* box_eta = new TBox(etarange_low, ymin, etarange_high, ymax);
0609         box_eta->SetLineColorAlpha(4, 0.5);
0610         box_eta->SetLineStyle(kDotted);
0611         box_eta->SetFillColorAlpha(4, 0.35);
0612         box_eta->Draw("same");
0613 
0614         TBox* box_pi0bkgr_low = new TBox(pi0bkgr_low_low, ymin, pi0bkgr_low_high, ymax);
0615         box_pi0bkgr_low->SetLineColorAlpha(46, 0.5);
0616         box_pi0bkgr_low->SetLineStyle(kDotted);
0617         box_pi0bkgr_low->SetFillColorAlpha(46, 0.35);
0618         box_pi0bkgr_low->Draw("same");
0619         TBox* box_pi0bkgr_high = new TBox(pi0bkgr_high_low, ymin, pi0bkgr_high_high, ymax);
0620         box_pi0bkgr_high->SetLineColorAlpha(46, 0.5);
0621         box_pi0bkgr_high->SetLineStyle(kDotted);
0622         box_pi0bkgr_high->SetFillColorAlpha(46, 0.35);
0623         box_pi0bkgr_high->Draw("same");
0624         TBox* box_etabkgr_low = new TBox(etabkgr_low_low, ymin, etabkgr_low_high, ymax);
0625         box_etabkgr_low->SetLineColorAlpha(38, 0.5);
0626         box_etabkgr_low->SetLineStyle(kDotted);
0627         box_etabkgr_low->SetFillColorAlpha(38, 0.35);
0628         box_etabkgr_low->Draw("same");
0629         TBox* box_etabkgr_high = new TBox(etabkgr_high_low, ymin, etabkgr_high_high, ymax);
0630         box_etabkgr_high->SetLineColorAlpha(38, 0.5);
0631         box_etabkgr_high->SetLineStyle(kDotted);
0632         box_etabkgr_high->SetFillColorAlpha(38, 0.35);
0633         box_etabkgr_high->Draw("same");
0634 
0635         c1->Modified();
0636         c1->SaveAs(outpdfname.c_str());
0637         hist->GetListOfFunctions()->Clear();
0638         delete fit_pi0_latex; delete fit_eta_latex;
0639         delete fit_pi0; delete fit_pi0signal; delete fit_pi0bkgr;
0640         delete fit_eta; delete fit_etabkgr; delete fit_etasignal;
0641         delete htemp; delete htemp2;
0642     }
0643     }
0644     delete c1;
0645 }
0646 
0647 TGraphErrors* BinnedHistSet::GetEnergyResolution() {
0648     TGraphErrors* gr_res = new TGraphErrors(nbins);
0649     gr_res->SetName(Form("gr_res_%s", name_base.c_str()));
0650     gr_res->SetTitle(Form("Energy Resolution vs %s;%s;Resolution", binned_quantity.c_str(), binned_quantity.c_str()));
0651     std::vector<double> bc = getBinCenters(bin_edges);
0652     for (int i=0; i<nbins; i++) {
0653     std::vector<double> mean_std_n = doubleFitHist(hist_vec.at(i+1), true);
0654     double x = bc.at(i);
0655     double y = mean_std_n.at(1)/mean_std_n.at(0);
0656     gr_res->SetPoint(i, x, y);
0657     double ex = 0;
0658     double ey = 0;
0659         if (mean_std_n.size() == 4) ey = mean_std_n.at(3);
0660     gr_res->SetPointError(i, ex, ey);
0661     /* std::cout << Form("%f\t%f\t%f\n", x, y, ey); */
0662     }
0663     /* std::cout << "Greg info: done creating graph " << gr_res->GetName() << "\n"; */
0664     return gr_res;
0665 }
0666 
0667 void BinnedHistSet::PlotEnergyResolution(std::string outpdfname, double low, double high) {
0668     TCanvas *c1 = new TCanvas("c1", "c1",0,50,1600,900);
0669     TGraphErrors* gr = this->GetEnergyResolution();
0670     gr->SetMarkerColor(1);
0671     gr->SetMarkerStyle(20);
0672     gr->Draw("AP");
0673     gr->GetYaxis()->SetRangeUser(low, high);
0674     gr->GetXaxis()->CenterTitle(true);
0675     gr->GetYaxis()->CenterTitle(true);
0676     c1->Modified();
0677     c1->SaveAs(outpdfname.c_str());
0678     delete c1; delete gr;
0679 }
0680 
0681 TGraphErrors* BinnedHistSet::GetBinnedMeansWide() {
0682     TGraphErrors* gr_wide = new TGraphErrors(nbins);
0683     gr_wide->SetName(Form("gr_wide_%s", name_base.c_str()));
0684     /* std::string xlabel = binned_quantity + GetBinnedQuantityUnits(); */
0685     /* gr_wide->SetTitle(Form("%s vs %s;%s;%s", this->GetTitlePrefix().c_str(), binned_quantity.c_str(), xlabel.c_str(), this->GetTitlePrefix().c_str())); */
0686     gr_wide->SetTitle(Form("%s vs %s;%s;%s", this->GetTitlePrefix().c_str(), binned_quantity.c_str(), binned_quantity.c_str(), this->GetTitlePrefix().c_str()));
0687     std::vector<double> bc = getBinCenters(bin_edges);
0688     /* // to print bin centers and corresponding y, yerr */
0689     /* std::cout << "Greg info: in GetBinnedMeansWide for " << name_base << ", nbins=" << nbins << ", bin centers are\n"; */
0690     /* printEdges(bc); */
0691     /* std::cout << "x\ty\tey\n"; */
0692     for (int i=0; i<nbins; i++) {
0693     std::vector<double> mean_std_n = doubleFitHist(hist_vec.at(i+1));
0694     double x = bc.at(i);
0695     double y = mean_std_n.at(0);
0696     gr_wide->SetPoint(i, x, y);
0697     double ex = 0;
0698     double ey = mean_std_n.at(1);
0699     gr_wide->SetPointError(i, ex, ey);
0700     /* std::cout << Form("%f\t%f\t%f\n", x, y, ey); */
0701     }
0702     /* std::cout << "Greg info: done creating graph " << gr_wide->GetName() << "\n"; */
0703     return gr_wide;
0704 }
0705 
0706 TGraphErrors* BinnedHistSet::GetBinnedMeansNarrow() {
0707     TGraphErrors* gr_narrow = new TGraphErrors(nbins);
0708     gr_narrow->SetName(Form("gr_narrow_%s", name_base.c_str()));
0709     /* std::string xlabel = binned_quantity;// + GetBinnedQuantityUnits(); */
0710     /* gr_narrow->SetTitle(Form("%s vs %s;%s;%s", this->GetTitlePrefix().c_str(), binned_quantity.c_str(), xlabel.c_str(), this->GetTitlePrefix().c_str())); */
0711     gr_narrow->SetTitle(Form("%s vs %s;%s;%s", this->GetTitlePrefix().c_str(), binned_quantity.c_str(), binned_quantity.c_str(), this->GetTitlePrefix().c_str()));
0712     std::vector<double> bc = getBinCenters(bin_edges);
0713     for (int i=0; i<nbins; i++) {
0714     std::vector<double> mean_std_n = doubleFitHist(hist_vec.at(i+1));
0715     double x = bc.at(i);
0716     double y = mean_std_n.at(0);
0717     gr_narrow->SetPoint(i, x, y);
0718     double ex = 0;
0719     double ey;
0720         if (mean_std_n.at(2) == 0) ey = 0;
0721     else ey = mean_std_n.at(1)/sqrt(mean_std_n.at(2));
0722     gr_narrow->SetPointError(i, ex, ey);
0723     }
0724     /* std::cout << "Greg info: done creating graph " << gr_narrow->GetName() << "\n"; */
0725     return gr_narrow;
0726 }
0727 
0728 void BinnedHistSet::PlotBinnedMeans(std::string outpdfname, double low=0.85, double high=1.2) {
0729     /* std::cout << "Greg info: entering PlotBinnedMeans for bhs " << name_base << "\n"; */
0730     /* std::cout << "nbins = " << nbins << ", hist_vec size = " << hist_vec.size() << "\n"; */
0731     TCanvas *c1 = new TCanvas("c1", "c1",0,50,1600,900);
0732     TGraphErrors* gr_wide = this->GetBinnedMeansWide();
0733     TGraphErrors* gr_narrow = this->GetBinnedMeansNarrow();
0734     gr_wide->SetMarkerColor(1); gr_narrow->SetMarkerColor(1);
0735     gr_wide->SetMarkerStyle(20); gr_narrow->SetMarkerStyle(20);
0736     gr_wide->SetLineWidth(2);
0737     gr_wide->GetXaxis()->SetTitleSize(0.05);  gr_wide->GetYaxis()->SetTitleSize(0.05);
0738     gr_wide->GetXaxis()->SetTitleOffset(0.8); gr_wide->GetYaxis()->SetTitleOffset(0.8);
0739     /* gr_narrow->SetLineColor(kRed); */
0740     gr_wide->Draw("AP[]");
0741     gr_narrow->Draw("P");
0742     gr_wide->GetYaxis()->SetRangeUser(low, high);
0743     gr_wide->GetXaxis()->CenterTitle(true);
0744     gr_wide->GetYaxis()->CenterTitle(true);
0745     c1->Modified();
0746     c1->SaveAs(outpdfname.c_str());
0747     delete c1; delete gr_wide; delete gr_narrow;
0748 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!