parametric_fit/macros/fit_argusmodified_2.C

0001 #include "/sphenix/u/virgilemahaut/style/sPhenixStyle_Greg.C"
0002
0003 R__LOAD_LIBRARY(libCustomHistogramFit.so)
0004 #include <CustomHistogramFit/HistogramUtils.h>
0005 #include <CustomHistogramFit/FitFunctions.h>
0006 #include <CustomHistogramFit/ConstrainedFitter.h>
0007
0008 void fit_argusmodified_2(const std::string &production_date = "01312026",
0009                          const std::string &selection_label = "MBD_INCOMPLETE",
0010                          const std::string &inputfolder_mass = "",
0011                          const std::string &inputfolder_analysis = "",
0012                          const std::string &folder_inputs = "",
0013                          const std::string &outputfolder_ratios = "",
0014                          const std::string &plots_folder_template = "")
0015 {
0016   const int nPtBins = 9;
0017   const std::string pTBins[nPtBins + 1] = {"0", "1", "2", "3", "4", "5", "6", "7", "8", "9"};
0018   const std::string pTBounds[nPtBins][2] = {{"1","2"}, {"2","3"}, {"3", "4"}, {"4","5"}, {"5", "6"}, {"6","7"},{"7","8"},{"8","10"},{"10","20"}};
0019
0020   // Read mean pT values
0021   std::string input_csv_name = folder_inputs + "/input_pt_mean_" + production_date + "_" + selection_label + ".csv";
0022   float pTMeans[2][nPtBins] = {0};
0023   std::string line;
0024   std::stringstream linestream;
0025   std::string entry;
0026   std::ifstream input_csv;
0027   input_csv.open(input_csv_name);
0028   if (!input_csv) {
0029     std::cout << "File " << input_csv_name << " could not be opened." << std::endl;
0030     return;
0031   }
0032   std::getline(input_csv, line);
0033   std::stringstream line_pi0(line);
0034   for (int iPt = 0; iPt < nPtBins; iPt++) {
0035     std::getline(line_pi0, entry, ',');
0036     float pTMean = std::stof(entry);
0037     pTMeans[0][iPt] = pTMean;
0038   }
0039   std::getline(input_csv, line);
0040   std::stringstream line_eta(line);
0041   for (int iPt = 0; iPt < nPtBins; iPt++) {
0042     std::getline(line_eta, entry, ',');
0043     float pTMean = std::stof(entry);
0044     pTMeans[1][iPt] = pTMean;
0045   }
0046
0047   // Histogram to store the background ratios
0048
0049   TGraphErrors *pt_ratios_pi0 = new TGraphErrors();
0050   pt_ratios_pi0->SetName("pt_background_ratios_pi0");
0051   pt_ratios_pi0->SetTitle("#pi^{0} background ratio");
0052   pt_ratios_pi0->GetXaxis()->SetTitle("p_{T} [GeV]");
0053   pt_ratios_pi0->GetYaxis()->SetTitle("background fraction [%]");
0054   TGraphErrors *pt_ratios_eta = new TGraphErrors();
0055   pt_ratios_eta->SetTitle("#eta background ratio");
0056   pt_ratios_eta->GetXaxis()->SetTitle("p_{T} [GeV]");
0057   pt_ratios_eta->GetYaxis()->SetTitle("background ratio [%]");
0058   pt_ratios_eta->SetName("pt_background_ratios_eta");
0059
0060   std::string plots_folder = plots_folder_template + "/plots_mass_argusmodified_2/analysis_ana509/ana509_" + production_date + "_sigma_3/" + selection_label + "/";
0061
0062   gSystem->Exec(("mkdir -p " + plots_folder).c_str());
0063   double statSum = 0;
0064
0065   for (int iPt = 0; iPt < nPtBins + 1; iPt++)
0066   {
0067     std::string filename = inputfolder_analysis + "/analysis_complete_ana509_" + production_date + "_" + selection_label + ".root";
0068
0069     std::string hname = "h_pair_mass_pt_" + pTBins[iPt];
0070     if (iPt == nPtBins) hname = "h_pair_mass";
0071     gSystem->Exec(("mkdir -p " + plots_folder).c_str());
0072     SetsPhenixStyle();
0073
0074     // Open file
0075     TFile *f = TFile::Open(filename.c_str());
0076
0077     // Extract the unfitted histogram
0078     std::cout << "f = " << f << std::endl;
0079     TH1F *h_pair_mass = (TH1F*) f->Get(hname.c_str());
0080     std::cout << "hname = " << hname.c_str() << std::endl;
0081     std::cout << "h_pair_mass = " << h_pair_mass << std::endl;
0082     h_pair_mass->SetTitle("");
0083     h_pair_mass->GetXaxis()->SetTitle("M_{#gamma#gamma} [GeV]");
0084
0085     statSum += h_pair_mass->GetEntries();
0086     std::ofstream statList;
0087     statList.open((plots_folder + "/statList.txt"), ios::out | ios::app);
0088     statList << "bin : p_T in " << pTBins[iPt] << " GeV/c: nentries = " << h_pair_mass->GetEntries() << std::endl;
0089     statList.close();
0090
0091     // Rescale the histogram to make the fit easier
0092     double scaling_factor = h_pair_mass->Integral();
0093     h_pair_mass->Scale(1. / scaling_factor);
0094
0095
0096     // Determine when the histogram starts
0097     double thresh = HistogramUtils::GetMaxBinContentInRange(h_pair_mass, 0.005, 0.9) / 1000;
0098     double thresh_large = HistogramUtils::GetMaxBinContentInRange(h_pair_mass, 0.005, 0.9) / 10;
0099     std::cout << "thresholds (scaled): " << thresh << ", " << thresh_large << std::endl;
0100     const double offset = h_pair_mass->GetBinCenter(HistogramUtils::FindFirstOccupiedBin(h_pair_mass, thresh));
0101     const double offset_large = std::max(h_pair_mass->GetBinCenter(HistogramUtils::FindFirstOccupiedBin(h_pair_mass, thresh_large)), 0.10);
0102
0103     if (true) {
0104       std::cout << "offsets = (" << offset << ", " << offset_large << ")\n";
0105     }
0106
0107     // Define the limit between two fitting regions. There is some overlap
0108     double borderFitMin = 0.27;
0109     double borderFitMax = 0.35;
0110     double fitMax1 = 0.27;
0111     double fitMin2 = 0.35;
0112
0113     // Prior background fit
0114     ConstrainedFitter *fitMethodBkg1 = new ConstrainedFitter(h_pair_mass, {FitFunctions::FitModel::ARGUSMODIFIED});
0115     fitMethodBkg1->SetFitRange({{offset, offset_large}, {0.18, fitMax1}});
0116     if (iPt == 0) fitMethodBkg1->SetFitRange({{0.05, 0.11}, {0.20, 0.32}});
0117     fitMethodBkg1->SetInitialParams({0.002, offset, 0.15, 0.1, 0.1});
0118     fitMethodBkg1->SetParamBounds({{0.0001, 1.0}, {offset-0.001, offset+0.001}, {0.1, 0.3}, {0.001, 10}, {0.001, 10}});
0119     fitMethodBkg1->PerformFit();
0120     fitMethodBkg1->PrintResults();
0121     TF1 *fit_bkg_1 = fitMethodBkg1->GetFitFunction();
0122
0123     ConstrainedFitter *fitMethodBkg2 = new ConstrainedFitter(h_pair_mass, {FitFunctions::FitModel::POLY2});
0124     fitMethodBkg2->SetFitRange({{fitMin2, 0.45}, {0.75, 1.00}});
0125     fitMethodBkg2->SetInitialParams({0.5, 0.5, 0.5});
0126     fitMethodBkg2->SetParamBounds({{0.0, 10.0}, {0.0, 10.0}, {0.0, 10.0}});
0127     fitMethodBkg2->PerformFit();
0128     fitMethodBkg2->PrintResults();
0129     TF1 *fit_bkg_2 = fitMethodBkg2->GetFitFunction();
0130
0131     // pi0 fit, fixed background
0132     ConstrainedFitter *fitMethodPi0 = new ConstrainedFitter(h_pair_mass, {FitFunctions::FitModel::ARGUSMODIFIED, FitFunctions::FitModel::GAUSS});
0133     fitMethodPi0->SetFitRange({{offset, fitMax1}});
0134     fitMethodPi0->SetInitialParams({
0135         fit_bkg_1->GetParameter(0),
0136         offset,
0137         fit_bkg_1->GetParameter(2),
0138         fit_bkg_1->GetParameter(3),
0139         fit_bkg_1->GetParameter(4),
0140         0.01, 0.140, 0.02}); // Gaussian initial parameters
0141     fitMethodPi0->SetParamBounds({
0142         {fit_bkg_1->GetParameter(0) - 0.0001, fit_bkg_1->GetParameter(0) + 0.0001},
0143         {offset - 0.0001, offset + 0.0001},
0144         {fit_bkg_1->GetParameter(2) - 0.0001, fit_bkg_1->GetParameter(2) + 0.0001},
0145         {fit_bkg_1->GetParameter(3) - 0.0001, fit_bkg_1->GetParameter(2) + 0.0001},
0146         {fit_bkg_1->GetParameter(4) - 0.0001, fit_bkg_1->GetParameter(3) + 0.0001},
0147         {0.001, 1}, {0.135, 0.160}, {0.015, 0.025}}); // Gaussian parameter limits
0148     fitMethodPi0->PerformFit();
0149     fitMethodPi0->PrintResults();
0150     TF1 *fit_pi0 = fitMethodPi0->GetFitFunction();
0151
0152     // eta fit, fixed background
0153     ConstrainedFitter *fitMethodEta = new ConstrainedFitter(h_pair_mass, {FitFunctions::FitModel::POLY2, FitFunctions::FitModel::GAUSS});
0154     fitMethodEta->SetFitRange({{fitMin2, 0.75}});
0155     fitMethodEta->SetInitialParams({
0156         fit_bkg_2->GetParameter(0),
0157         fit_bkg_2->GetParameter(1),
0158         fit_bkg_2->GetParameter(2),
0159         0.01, 0.580, 0.05}); // Gaussian initial parameters
0160     fitMethodEta->SetParamBounds({
0161         {fit_bkg_2->GetParameter(0) - 0.0001, fit_bkg_2->GetParameter(0) + 0.0001},
0162         {fit_bkg_2->GetParameter(1) - 0.0001, fit_bkg_2->GetParameter(1) + 0.0001},
0163         {fit_bkg_2->GetParameter(2) - 0.0001, fit_bkg_2->GetParameter(2) + 0.0001},
0164         {0.0001, 1}, {0.560, 0.590}, {0.045, 0.055}}); // Gaussian parameter limits
0165     fitMethodEta->PerformFit();
0166     fitMethodEta->PrintResults();
0167     TF1 *fit_eta = fitMethodEta->GetFitFunction();
0168
0169
0170     // bkg, pi0 and eta fit
0171     ConstrainedFitter *fitMethodGlobal1 = new ConstrainedFitter(h_pair_mass, {FitFunctions::FitModel::ARGUSMODIFIED, FitFunctions::FitModel::GAUSS});
0172     fitMethodGlobal1->SetFitRange({{offset, fitMax1}});
0173     fitMethodGlobal1->SetInitialParams({
0174         fit_pi0->GetParameter(0),
0175         offset,
0176         fit_pi0->GetParameter(2),
0177         fit_pi0->GetParameter(3),
0178         fit_pi0->GetParameter(4),
0179         fit_pi0->GetParameter(5),
0180         fit_pi0->GetParameter(6),
0181         fit_pi0->GetParameter(7)
0182       });
0183     fitMethodGlobal1->SetParamBounds({
0184         {0.0001, 1.0}, {offset-0.001, offset+0.001}, {0.1, 0.3}, {0.001, 10}, {0.001, 10},
0185         {0.001, 1}, {0.135, 0.160}, {0.015, 0.025}, // Pi0 Gaussian parameter limits
0186       });
0187
0188     fitMethodGlobal1->PerformFit();
0189     fitMethodGlobal1->PrintResults();
0190     TF1 *fit_global_1 = fitMethodGlobal1->GetFitFunction();
0191
0192     ConstrainedFitter *fitMethodGlobal2 = new ConstrainedFitter(h_pair_mass, {FitFunctions::FitModel::POLY2, FitFunctions::FitModel::GAUSS});
0193     fitMethodGlobal2->SetFitRange({{fitMin2, 1.0}});
0194     fitMethodGlobal2->SetInitialParams({
0195         fit_eta->GetParameter(0),
0196         fit_eta->GetParameter(1),
0197         fit_eta->GetParameter(2),
0198         fit_eta->GetParameter(3),
0199         fit_eta->GetParameter(4),
0200         fit_eta->GetParameter(5)
0201       });
0202     fitMethodGlobal2->SetParamBounds({
0203         {0, 10}, {0, 10}, {0, 10},
0204         {0.0001, 1}, {0.560, 0.590}, {0.045, 0.055} // eta Gaussian parameter limits
0205       });
0206     fitMethodGlobal2->PerformFit();
0207     fitMethodGlobal2->PrintResults();
0208     TF1 *fit_global_2 = fitMethodGlobal2->GetFitFunction();
0209
0210     // Rescale
0211     fit_global_1->SetParameter(0, fit_global_1->GetParameter(0) * scaling_factor);
0212     fit_global_1->SetParameter(5, fit_global_1->GetParameter(5) * scaling_factor);
0213
0214     fit_global_2->SetParameter(0, fit_global_2->GetParameter(0) * scaling_factor);
0215     fit_global_2->SetParameter(1, fit_global_2->GetParameter(1) * scaling_factor);
0216     fit_global_2->SetParameter(2, fit_global_2->GetParameter(2) * scaling_factor);
0217     fit_global_2->SetParameter(3, fit_global_2->GetParameter(3) * scaling_factor);
0218
0219     // Extract the bkg, pi0 and eta contribution from the last fit
0220     ConstrainedFitter *fitMethodGlobalBkg1 = new ConstrainedFitter(h_pair_mass, {FitFunctions::FitModel::ARGUSMODIFIED});
0221     fitMethodGlobalBkg1->SetInitialParams({
0222         fit_global_1->GetParameter(0),
0223         fit_global_1->GetParameter(1),
0224         fit_global_1->GetParameter(2),
0225         fit_global_1->GetParameter(3),
0226         fit_global_1->GetParameter(4)
0227       });
0228     TF1 *fit_global_bkg_1 = fitMethodGlobalBkg1->GetFitFunction();
0229
0230     ConstrainedFitter *fitMethodGlobalBkg2 = new ConstrainedFitter(h_pair_mass, {FitFunctions::FitModel::POLY2});
0231     fitMethodGlobalBkg2->SetInitialParams({
0232         fit_global_2->GetParameter(0),
0233         fit_global_2->GetParameter(1),
0234         fit_global_2->GetParameter(2)
0235       });
0236     TF1 *fit_global_bkg_2 = fitMethodGlobalBkg2->GetFitFunction();
0237     fit_global_bkg_2->SetName("fit_global_bkg_2");
0238
0239
0240     //auto lambda_scaled_global_bkg_2 = [scaling_factor, fit_global_bkg_2] (double *x, double *p) { return scaling_factor * fit_global_bkg_2->Eval(x[0]); };
0241     //TF1 *fit_scaled_global_bkg_2 = new TF1("fit_scaled_global_bkg_2", lambda_scaled_global_bkg_2);
0242     /*TF1 *fit_scaled_global_bkg_2 = new TF1("fit_scaled_global_bkg_2", "[0] * fit_global_bkg_2");
0243       fit_scaled_global_bkg_2->SetParameter(0, scaling_factor);*/
0244
0245     ConstrainedFitter *fitMethodGlobalPi0 = new ConstrainedFitter(h_pair_mass, {FitFunctions::FitModel::GAUSS});
0246     fitMethodGlobalPi0->SetInitialParams({
0247         fit_global_1->GetParameter(5),
0248         fit_global_1->GetParameter(6),
0249         fit_global_1->GetParameter(7)
0250       });
0251     TF1 *fit_global_pi0 = fitMethodGlobalPi0->GetFitFunction();
0252
0253     ConstrainedFitter *fitMethodGlobalEta = new ConstrainedFitter(h_pair_mass, {FitFunctions::FitModel::GAUSS});
0254     fitMethodGlobalEta->SetInitialParams({
0255         fit_global_2->GetParameter(3),
0256         fit_global_2->GetParameter(4),
0257         fit_global_2->GetParameter(5)
0258       });
0259     TF1 *fit_global_eta = fitMethodGlobalEta->GetFitFunction();
0260
0261
0262     // Define transparent bands for the peak regions and the side bands
0263     const int nParticles = 2;
0264     const int nRegions = 3;
0265
0266     // 3 sigma
0267     float band_bounds[nParticles][nRegions][2] =
0268     {{{0.080,0.199}, {0.030, 0.070}, {0.209, 0.249}},
0269     {{0.399,0.739}, {0.257, 0.371}, {0.767, 0.880}}};
0270
0271     // 2 sigma
0272     /*float band_bounds[nParticles][nRegions][2] =
0273     {{{0.100,0.180}, {0.030, 0.070}, {0.209, 0.249}},
0274     {{0.456,0.683}, {0.257, 0.371}, {0.767, 0.880}}};*/
0275
0276     // 1.5 sigma
0277     /*float band_bounds[nParticles][nRegions][2] =
0278     {{{0.110,0.170}, {0.030, 0.070}, {0.209, 0.249}},
0279     {{0.484,0.654}, {0.257, 0.371}, {0.767, 0.880}}};*/
0280
0281     h_pair_mass->Scale(scaling_factor);
0282     // Compute ratios
0283     int firstbin, secondbin;
0284     firstbin = h_pair_mass->FindBin(band_bounds[0][0][0]);
0285     secondbin = h_pair_mass->FindBin(band_bounds[0][0][1]);
0286     float pi0_ratio = fit_global_bkg_1->Integral(band_bounds[0][0][0], band_bounds[0][0][1]) /
0287       fit_global_1->Integral(band_bounds[0][0][0], band_bounds[0][0][1]);
0288       // h_pair_mass->Integral(firstbin, secondbin) * (secondbin - firstbin) /
0289       // (band_bounds[0][0][1] - band_bounds[0][0][0]);
0290     float pi0_ratio_1 = fit_global_bkg_1->Integral(band_bounds[0][0][0] + 0.002, band_bounds[0][0][1] + 0.002) /
0291       fit_global_1->Integral(band_bounds[0][0][0] + 0.002, band_bounds[0][0][1] + 0.002);
0292       // h_pair_mass->Integral(firstbin+1, secondbin+1) * (secondbin - firstbin) /
0293       // (band_bounds[0][0][1] - band_bounds[0][0][0]);
0294     float pi0_ratio_2 = fit_global_bkg_1->Integral(band_bounds[0][0][0] - 0.002, band_bounds[0][0][1] - 0.002) /
0295       fit_global_1->Integral(band_bounds[0][0][0] - 0.002, band_bounds[0][0][1] - 0.002);
0296       // h_pair_mass->Integral(firstbin-1, secondbin-1) * (secondbin - firstbin) /
0297       // (band_bounds[0][0][1] - band_bounds[0][0][0]);
0298     float pi0_ratio_err = (std::abs(pi0_ratio_1 - pi0_ratio) + std::abs(pi0_ratio_2 - pi0_ratio)) / 2;
0299     //fit_global_1->Integral(band_bounds[0][0][0], band_bounds[0][0][1]);
0300     float pi0_signal_range = h_pair_mass->Integral(firstbin, secondbin);
0301
0302     firstbin = h_pair_mass->FindBin(band_bounds[1][0][0]);
0303     secondbin = h_pair_mass->FindBin(band_bounds[1][0][1]);
0304     float eta_ratio = fit_global_bkg_2->Integral(band_bounds[1][0][0], band_bounds[1][0][1]) /
0305       fit_global_2->Integral(band_bounds[1][0][0], band_bounds[1][0][1]);
0306       // h_pair_mass->Integral(firstbin, secondbin) * (secondbin - firstbin) /
0307       // (band_bounds[1][0][1] - band_bounds[1][0][0]);
0308     float eta_ratio_1 = fit_global_bkg_2->Integral(band_bounds[1][0][0] + 0.002, band_bounds[1][0][1] + 0.002) /
0309       fit_global_2->Integral(band_bounds[1][0][0] + 0.002, band_bounds[1][0][1] + 0.002);
0310       // h_pair_mass->Integral(firstbin + 1, secondbin + 1) * (secondbin - firstbin) /
0311       // (band_bounds[1][0][1] - band_bounds[1][0][0]);
0312     float eta_ratio_2 = fit_global_bkg_2->Integral(band_bounds[1][0][0] - 0.002, band_bounds[1][0][1] - 0.002) /
0313       fit_global_2->Integral(band_bounds[1][0][0] - 0.002, band_bounds[1][0][1] - 0.002);
0314       // h_pair_mass->Integral(firstbin + 1, secondbin + 1) * (secondbin - firstbin) /
0315       // (band_bounds[1][0][1] - band_bounds[1][0][0]);
0316     float eta_ratio_err = (std::abs(eta_ratio_1 - eta_ratio) + std::abs(eta_ratio_1 - eta_ratio)) / 2;
0317     //fit_global_2->Integral(band_bounds[1][0][0], band_bounds[1][0][1]);
0318     float eta_signal_range = h_pair_mass->Integral(firstbin, secondbin);
0319
0320
0321     // Access easily the total number of pi0 and eta
0322     if (iPt == 9) {
0323       float pi0_signal_stat = pi0_signal_range * (1 - pi0_ratio);
0324       float eta_signal_stat = eta_signal_range * (1 - eta_ratio);
0325       std::cout << "Total Number of entries = " << h_pair_mass->GetEntries() << std::endl;
0326       std::cout << "Total Number Integral = " << h_pair_mass->Integral() << std::endl;
0327       std::cout << "Total Number of pi0/eta mesons = " << pi0_signal_stat << "/" << eta_signal_stat << std::endl;
0328     }
0329
0330
0331     if (iPt < nPtBins) {
0332       //if (iPt > 1 && iPt < nPtBins - 1) {
0333       // Skip the first bins in pT (not included in the sPHENIX preliminary release)
0334       if (pi0_ratio > 0 && pi0_ratio < 1) {
0335         pt_ratios_pi0->SetPoint(iPt, pTMeans[0][iPt], pi0_ratio * 100);
0336         pt_ratios_pi0->SetPointError(iPt, 0, pi0_ratio_err * 100);
0337       }
0338       else {
0339         pt_ratios_pi0->SetPoint(iPt, pTMeans[0][iPt], 100);
0340         pt_ratios_pi0->SetPointError(iPt, 0, 0);
0341       }
0342       if (eta_ratio > 0 && eta_ratio < 1) {
0343         pt_ratios_eta->SetPoint(iPt, pTMeans[0][iPt], eta_ratio * 100);
0344         pt_ratios_eta->SetPointError(iPt, 0, eta_ratio_err * 100);
0345       }
0346       else {
0347         pt_ratios_eta->SetPoint(iPt, pTMeans[0][iPt], 100);
0348         pt_ratios_eta->SetPointError(iPt, 0, 0);
0349       }
0350     }
0351
0352     std::ofstream ratioList;
0353     ratioList.open((plots_folder + "/ratioList.txt"), ios::out | ios::app);
0354     ratioList << "bin : p_T in " << pTBins[iPt] << " GeV/c: r_{pi0} = " << pi0_ratio << ", r_{eta} = " << eta_ratio << std::endl;
0355     ratioList.close();
0356
0357     h_pair_mass->SetMinimum(0);
0358     h_pair_mass->GetYaxis()->SetTitle("counts / [2 MeV]");
0359     if (iPt == 8){
0360       double currentMax = h_pair_mass->GetMaximum();
0361       h_pair_mass->SetMaximum(3 * currentMax);
0362       h_pair_mass->SetMinimum(0);
0363     }
0364     if (iPt == 7){
0365       double currentMax = h_pair_mass->GetMaximum();
0366       h_pair_mass->SetMaximum(1.5 * currentMax);
0367       h_pair_mass->SetMinimum(0);
0368     }
0369
0370     //if (iPt == 2) h_pair_mass->SetMaximum(800000);
0371
0372     h_pair_mass->Draw();
0373     gPad->Update();
0374     double max_y = gPad->GetUymax();
0375
0376     TBox *pi0PeakBand = new TBox(band_bounds[0][0][0], 0, band_bounds[0][0][1], max_y);
0377     pi0PeakBand->SetFillColorAlpha(kRed, 0.4);
0378     TBox *pi0SideBandLeft = new TBox(band_bounds[0][1][0], 0, band_bounds[0][1][1], max_y);
0379     pi0SideBandLeft->SetFillColorAlpha(kRed, 0.2);
0380     TBox *pi0SideBandRight = new TBox(band_bounds[0][2][0], 0, band_bounds[0][2][1], max_y);
0381     pi0SideBandRight->SetFillColorAlpha(kRed, 0.2);
0382
0383     TBox *etaPeakBand = new TBox(band_bounds[1][0][0], 0, band_bounds[1][0][1], max_y);
0384     etaPeakBand->SetFillColorAlpha(kBlue, 0.4);
0385     TBox *etaSideBandLeft = new TBox(band_bounds[1][1][0], 0, band_bounds[1][1][1], max_y);
0386     etaSideBandLeft->SetFillColorAlpha(kBlue, 0.2);
0387     TBox *etaSideBandRight = new TBox(band_bounds[1][2][0], 0, band_bounds[1][2][1], max_y);
0388     etaSideBandRight->SetFillColorAlpha(kBlue, 0.2);
0389
0390     // Store the fit results
0391     std::string outfilename_root = plots_folder + "/parametric_fit_mass_pt_" + std::to_string(iPt) + ".root";
0392     TFile *outfile_root = new TFile(outfilename_root.c_str(), "RECREATE");
0393     outfile_root->cd();
0394
0395     h_pair_mass->SetName("h_pair_mass");
0396     h_pair_mass->Write();
0397
0398     fit_global_bkg_1->SetName("fit_global_bkg_pi0");
0399     fit_global_bkg_1->SetRange(offset, borderFitMin);
0400     fit_global_bkg_1->SetNpx(1000);
0401     fit_global_bkg_1->SetLineWidth(2);
0402     fit_global_bkg_1->SetLineColor(kGreen);
0403     fit_global_bkg_1->SetLineStyle(2);
0404     fit_global_bkg_1->Write();
0405
0406     fit_global_bkg_2->SetName("fit_global_bkg_eta");
0407     fit_global_bkg_2->SetRange(borderFitMax, 1.0);
0408     fit_global_bkg_2->SetNpx(1000);
0409     fit_global_bkg_2->SetLineWidth(2);
0410     fit_global_bkg_2->SetLineColor(kGreen);
0411     fit_global_bkg_2->SetLineStyle(2);
0412     fit_global_bkg_2->Write();
0413
0414     fit_global_pi0->SetName("fit_global_pi0");
0415     fit_global_pi0->SetRange(offset, 1.0);
0416     fit_global_pi0->SetNpx(1000);
0417     fit_global_pi0->SetLineWidth(2);
0418     fit_global_pi0->SetLineColor(kBlue);
0419     fit_global_pi0->SetLineStyle(kSolid);
0420     fit_global_pi0->Write();
0421
0422     fit_global_eta->SetName("fit_global_eta");
0423     fit_global_eta->SetRange(offset, 1.0);
0424     fit_global_eta->SetNpx(1000);
0425     fit_global_eta->SetLineWidth(2);
0426     fit_global_eta->SetLineColor(kBlue);
0427     fit_global_eta->SetLineStyle(kSolid);
0428     fit_global_eta->Write();
0429     outfile_root->Close();
0430     delete outfile_root;
0431
0432     // Print the histogram and the fits
0433     TFile *outfile = new TFile("fits.root", "RECREATE");
0434     std::string canvas_name = "canvas_pi0_eta_mass_pt_" + pTBins[iPt];
0435     TCanvas *canvas = new TCanvas(canvas_name.c_str(), "c", 1600, 900);
0436
0437     // Draw original histogram
0438     h_pair_mass->SetLineColor(kBlack);
0439     h_pair_mass->SetStats(0);
0440     h_pair_mass->Draw("E");
0441
0442     // Draw fits
0443     /*fit_bkg_1->SetRange(offset, borderFit);
0444     fit_bkg_1->SetNpx(1000);
0445     fit_bkg_1->SetLineWidth(3);
0446     fit_bkg_1->SetLineColor(kGreen);
0447     fit_bkg_1->SetLineStyle(kDashed);
0448     fit_bkg_1->Draw("same L");
0449
0450     fit_bkg_2->SetRange(borderFit,1);
0451     fit_bkg_2->SetNpx(1000);
0452     fit_bkg_2->SetLineWidth(3);
0453     fit_bkg_2->SetLineColor(kGreen);
0454     fit_bkg_2->SetLineStyle(kDashed);
0455     fit_bkg_2->Draw("same L");
0456
0457     fit_pi0->SetRange(offset, borderFit);
0458     fit_pi0->SetNpx(1000);
0459     fit_pi0->SetLineWidth(3);
0460     fit_pi0->SetLineColor(kBlue);
0461     fit_pi0->SetLineStyle(kDashed);
0462     fit_pi0->Draw("same L");
0463
0464     fit_eta->SetRange(borderFit, 1.0);
0465     fit_eta->SetNpx(1000);
0466     fit_eta->SetLineWidth(3);
0467     fit_eta->SetLineColor(kYellow);
0468     fit_eta->SetLineStyle(kSolid);
0469     fit_eta->Draw("same L");*/
0470
0471     fit_global_bkg_1->SetRange(offset, borderFitMin);
0472     fit_global_bkg_1->SetNpx(1000);
0473     fit_global_bkg_1->SetLineWidth(2);
0474     fit_global_bkg_1->SetLineColor(kGreen);
0475     fit_global_bkg_1->SetLineStyle(2);
0476     fit_global_bkg_1->Draw("same L");
0477
0478     fit_global_bkg_2->SetRange(borderFitMax, 1.0);
0479     fit_global_bkg_2->SetNpx(1000);
0480     fit_global_bkg_2->SetLineWidth(2);
0481     fit_global_bkg_2->SetLineColor(kGreen);
0482     fit_global_bkg_2->SetLineStyle(2);
0483     fit_global_bkg_2->Draw("same L");
0484
0485     fit_global_pi0->SetRange(offset, 1.0);
0486     fit_global_pi0->SetNpx(1000);
0487     fit_global_pi0->SetLineWidth(2);
0488     fit_global_pi0->SetLineColor(kBlue);
0489     fit_global_pi0->SetLineStyle(kSolid);
0490     fit_global_pi0->Draw("same L");
0491
0492     fit_global_eta->SetRange(offset, 1.0);
0493     fit_global_eta->SetNpx(1000);
0494     fit_global_eta->SetLineWidth(2);
0495     fit_global_eta->SetLineColor(kBlue);
0496     fit_global_eta->SetLineStyle(kSolid);
0497     fit_global_eta->Draw("same L");
0498
0499     fit_global_1->SetRange(offset, borderFitMin);
0500     fit_global_1->SetNpx(1000);
0501     fit_global_1->SetLineWidth(2);
0502     fit_global_1->SetLineColor(kRed);
0503     fit_global_1->SetLineStyle(kSolid);
0504     fit_global_1->Draw("same L");
0505
0506     fit_global_2->SetRange(borderFitMax, 1.0);
0507     fit_global_2->SetNpx(1000);
0508     fit_global_2->SetLineWidth(2);
0509     fit_global_2->SetLineColor(kRed);
0510     fit_global_2->SetLineStyle(kSolid);
0511     fit_global_2->Draw("same L");
0512
0513     // Draw peak regions
0514     pi0PeakBand->Draw("same");
0515     pi0SideBandLeft->Draw("same");
0516     pi0SideBandRight->Draw("same");
0517     etaPeakBand->Draw("same");
0518     etaSideBandLeft->Draw("same");
0519     etaSideBandRight->Draw("same");
0520
0521     // Show relevant numerical information
0522     TPad *p = new TPad("p","p",0.,0.,1.,1.); p->SetFillStyle(0); p->Draw(); p->cd();
0523     TBox *whiteBox = new TBox(0.67, 0.72, 0.88, 0.9);
0524     whiteBox->Draw();
0525     canvas->cd();
0526     whiteBox->SetFillColorAlpha(kWhite, 1);
0527     std::stringstream stream;
0528     TLatex latex;
0529     latex.SetNDC();
0530     latex.SetTextColor(kBlack);
0531     latex.DrawLatex(0.68, 0.85, "#font[72]{sPHENIX} Internal");
0532     latex.DrawLatex(0.68, 0.75, "p^{#uparrow}+p #sqrt{s} = 200 GeV");
0533     //std::cout << "text size: " << latex.GetTextSize() << std::endl;
0534     //latex.SetTextSize(18);
0535     if (iPt < nPtBins) {
0536       latex.DrawLatex(0.47,0.87, (pTBounds[iPt][0] + " < p_{T} [GeV] < " + pTBounds[iPt][1]).c_str());
0537     }
0538     /*latex.DrawLatex(0.5, 0.80,"#pi^{0}");
0539     stream << "mean: " << std::fixed << std::setprecision(5) << fit_global_pi0->GetParameter(1) << " #pm " << fit_global_1->GetParError(5);
0540     latex.DrawLatex(0.5, 0.75, stream.str().c_str());
0541     stream.str(""); stream << "#sigma: " << std::setprecision(5) << std::abs(fit_global_pi0->GetParameter(2)) << " #pm " << fit_global_1->GetParError(6);
0542     latex.DrawLatex(0.5, 0.70, stream.str().c_str());
0543     stream.str(""); stream << "r: " << std::setprecision(5) << pi0_ratio  << " #pm " << pi0_ratio_err;
0544     latex.DrawLatex(0.5, 0.65, stream.str().c_str());
0545     latex.DrawLatex(0.5, 0.55,"#eta");
0546     stream.str(""); stream << "mean: " << std::setprecision(5) << fit_global_eta->GetParameter(1) << " #pm " << fit_global_2->GetParError(5);
0547     latex.DrawLatex(0.5, 0.50, stream.str().c_str());
0548     stream.str(""); stream << "#sigma: " << std::setprecision(5) << std::abs(fit_global_eta->GetParameter(2)) << " #pm " << fit_global_2->GetParError(6);
0549     latex.DrawLatex(0.5, 0.45, stream.str().c_str());
0550     stream.str(""); stream << "r: " << std::setprecision(5) << eta_ratio << " #pm " << eta_ratio_err;
0551     latex.DrawLatex(0.5, 0.40, stream.str().c_str());*/
0552
0553     // Add the full legend
0554     TLegend *full_legend = new TLegend(0.55, 0.5, 0.85, 0.70);
0555     full_legend->SetTextSize(0.05);
0556     full_legend->AddEntry(h_pair_mass, "Data");
0557     full_legend->AddEntry(fit_global_1, "Full Fit", "l");
0558     full_legend->AddEntry(fit_global_bkg_1, "Fitted Background", "l");
0559     full_legend->AddEntry(fit_global_pi0, "Fitted Signal", "l");
0560     full_legend->Draw();
0561
0562     // Save histograms
0563     canvas->SaveAs((plots_folder + "/" + canvas_name + ".pdf").c_str());
0564     canvas->SaveAs((plots_folder + "/" + canvas_name + ".png").c_str());
0565
0566     // Close file
0567     outfile->cd();
0568     canvas->Write();
0569     outfile->Close();
0570     delete outfile;
0571
0572     delete fitMethodBkg1;
0573     delete fitMethodBkg2;
0574     delete fitMethodPi0;
0575     delete fitMethodEta;
0576     delete fitMethodGlobal1;
0577     delete fitMethodGlobal2;
0578     delete fitMethodGlobalBkg1;
0579     delete fitMethodGlobalBkg2;
0580     delete fitMethodGlobalPi0;
0581     delete fitMethodGlobalEta;
0582   }
0583
0584   // Draw background ratios
0585   std::string canvas_name = "canvas_pt_background_ratios";
0586   TCanvas *canvas = new TCanvas(canvas_name.c_str(), "c", 1600, 900);
0587   TMultiGraph *mg = new TMultiGraph();
0588   mg->GetXaxis()->SetTitle("p_{T} [GeV]");
0589   mg->GetYaxis()->SetTitle("Background Fraction [%]");
0590   pt_ratios_pi0->SetMarkerSize(2.5);
0591   pt_ratios_pi0->SetMarkerStyle(kFullCircle);
0592   pt_ratios_pi0->SetMarkerColor(kBlue);
0593   pt_ratios_eta->SetMarkerSize(2.5);
0594   pt_ratios_eta->SetMarkerStyle(kFullSquare);
0595   pt_ratios_eta->SetMarkerColor(kRed);
0596   mg->Add(pt_ratios_pi0);
0597   mg->Add(pt_ratios_eta);
0598   mg->SetMinimum(0);
0599   mg->SetMaximum(100);
0600   mg->GetXaxis()->SetRangeUser(2, 20);
0601   mg->Draw("AP");
0602
0603   TLatex latex;
0604   latex.SetNDC();
0605   latex.SetTextColor(kBlack);
0606   latex.DrawLatex(0.7, 0.75, "#font[72]{sPHENIX} Internal");
0607   latex.DrawLatex(0.7, 0.65, "p^{#uparrow}+p #sqrt{s} = 200 GeV");
0608
0609   TLegend *legend = new TLegend(0.2, 0.5, 0.45, 0.7);
0610   //TLegend *legend = new TLegend(0.6, 0.2, 0.9, 0.4);
0611   legend->AddEntry(pt_ratios_pi0, "p^{#uparrow}p #rightarrow #pi^{0} X");
0612   legend->AddEntry(pt_ratios_eta, "p^{#uparrow}p #rightarrow #eta X");
0613   legend->SetBorderSize(1);
0614   legend->Draw();
0615
0616   canvas->SaveAs((plots_folder + "/" + canvas_name + ".pdf").c_str());
0617   canvas->SaveAs((plots_folder + "/" + canvas_name + ".png").c_str());
0618
0619   canvas_name = "canvas_pt_background_ratios_pi0";
0620   canvas = new TCanvas("canvas_pt_background_ratios_pi0", "c", 1600, 900);
0621   canvas->cd();
0622   pt_ratios_pi0->SetMinimum(0);
0623   pt_ratios_pi0->SetMaximum(100);
0624   //pt_ratios_pi0->SetLineWidth(0);
0625   pt_ratios_pi0->SetMarkerSize(1.5);
0626   pt_ratios_pi0->SetMarkerStyle(kFullCircle);
0627   pt_ratios_pi0->SetMarkerColor(kBlue);
0628   pt_ratios_pi0->Draw("AP");
0629
0630   canvas->SaveAs((plots_folder + "/" + canvas_name + ".pdf").c_str());
0631   canvas->SaveAs((plots_folder + "/" + canvas_name + ".png").c_str());
0632   delete canvas;
0633
0634   canvas_name = "canvas_pt_background_ratios_eta";
0635   canvas = new TCanvas("canvas_pt_background_ratios_eta", "c", 1600, 900);
0636   canvas->cd();
0637   pt_ratios_eta->SetMinimum(70);
0638   pt_ratios_eta->SetMaximum(100);
0639   //pt_ratios_eta->SetLineWidth(0);
0640   pt_ratios_eta->SetMarkerSize(1.5);
0641   pt_ratios_eta->SetMarkerStyle(kFullCircle);
0642   pt_ratios_eta->SetMarkerColor(kRed);
0643   pt_ratios_eta->Draw();
0644
0645   std::ofstream ratio_csv;
0646   ratio_csv.open((outputfolder_ratios + "/" + production_date + "_" + selection_label + ".csv"));
0647   for (int iPt = 0; iPt < nPtBins - 1; iPt++) ratio_csv << pt_ratios_pi0->GetPointY(iPt) / 100. << ", ";
0648   ratio_csv << pt_ratios_pi0->GetPointY(nPtBins - 1) / 100. << "\n";
0649   for (int iPt = 0; iPt < nPtBins - 1; iPt++) ratio_csv << pt_ratios_eta->GetPointY(iPt) / 100. << ", ";
0650   ratio_csv << pt_ratios_eta->GetPointY(nPtBins - 1) / 100. << "\n";
0651   for (int iPt = 0; iPt < nPtBins - 1; iPt++) ratio_csv << 0  << ", ";
0652   ratio_csv << 0 << "\n";
0653   for (int iPt = 0; iPt < nPtBins - 1; iPt++) ratio_csv << 0 << ", ";
0654   ratio_csv << 0 << "\n";
0655   ratio_csv.close();
0656
0657   canvas->SaveAs((plots_folder + "/" + canvas_name + ".pdf").c_str());
0658   canvas->SaveAs((plots_folder + "/" + canvas_name + ".png").c_str());
0659   delete canvas;
0660
0661   gSystem->Exit(0);
0662 }