LightFlavorRatios/charged_hadron_ratios/pion_ratio_from_Kshort.C

0001 #include <phool/PHRandomSeed.h>
0002 #include <gsl/gsl_rng.h>
0003 #include <algorithm>
0004
0005 using namespace std;
0006
0007 /*
0008  * This creates the output directory for the plots
0009  */
0010 std::string plotPath = "./Kshort_daughter_ratios";
0011 string makeDirectory = "mkdir -p " + plotPath;
0012 system(makeDirectory.c_str());
0013
0014 /*
0015  * This creates an nth-order polynomial for fitting the ratio to
0016  * Change the order for higher or lower order fits
0017  */
0018 const unsigned int polynomial_order = 5;
0019 TF1* ratio_function;
0020 Double_t fitfunc(Double_t *x, Double_t *par)
0021 {
0022   double value = 0;
0023   for (unsigned int i = 0; i <= polynomial_order; ++i) value += par[i]*pow(x[0], i);
0024
0025   return value;
0026 }
0027
0028 /*
0029  * Truncates numbers to better visualisation on plots
0030  */
0031 template <typename T>
0032 string to_string_with_precision(const T a_value, const int n = 1)
0033 {
0034     ostringstream out;
0035     out.precision(n);
0036     out << fixed << a_value;
0037     return out.str();
0038 }
0039
0040 /*
0041  * Sets all the canvas stuff and saves the plots with certain extensions
0042  * just add more extensions to the list
0043  */
0044 template <typename T>
0045 void savePlots(T myPlot, string plotName, bool logY = false)
0046 {
0047   TCanvas *c1  = new TCanvas("myCanvas", "myCanvas",800,800);
0048
0049   myPlot.SetMinimum(0);
0050   myPlot.SetMaximum(2.8);
0051
0052   if (strncmp(typeid(myPlot).name(), "4TH2F", 5) == 0)
0053   {
0054     myPlot.Draw("COLZ");
0055   }
0056   else
0057   {
0058     myPlot.Sumw2();
0059     if (logY) gPad->SetLogy();
0060     myPlot.Draw("PE1");
0061     ratio_function->Draw("SAME");
0062   }
0063
0064   TPaveText *pt;
0065   pt = new TPaveText(0.15,0.9,0.95,1., "NDC");
0066   pt->SetFillColor(0);
0067   pt->SetFillStyle(0);
0068   pt->SetTextFont(42);
0069   TText *pt_LaTex = pt->AddText("#it{#bf{sPHENIX}} Internal, #sqrt{s} = 200 GeV, pp");
0070   pt->SetBorderSize(0);
0071   pt->Draw();
0072   gPad->Modified();
0073
0074   string extensions[] = {".C", ".pdf", ".png"};
0075   for (auto extension : extensions)
0076   {
0077     string output = plotPath + "/" + plotName + extension;
0078     c1->SaveAs(output.c_str());
0079   }
0080 }
0081
0082 /*
0083  * Makes a historgram template with variable bin widths as defined below
0084  */
0085 float lower_bin_bounds[] = {0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.4, 1.8, 2.2, 3, 4};
0086 const unsigned int n_variable_bins = sizeof(lower_bin_bounds)/sizeof(lower_bin_bounds[0]) - 1;
0087
0088 TH1F makeHisto(int nBins, float* min, string xAxisTitle, string unit, int precision, string yAxisTitle = "Raw Ratio")
0089 {
0090   TH1F myHisto(xAxisTitle.c_str(), xAxisTitle.c_str(), nBins, min);
0091
0092   if (unit != "") xAxisTitle += " [" + unit +  "]";
0093   myHisto.GetXaxis()->SetTitle(xAxisTitle.c_str());
0094   myHisto.GetYaxis()->SetTitle(yAxisTitle.c_str());
0095
0096   return myHisto;
0097 }
0098
0099 /*
0100  * Makes a historgram template with fixed bin widths and ranges defines below
0101  */
0102 float min_pT = 0.2;
0103 float max_pT = 3.2;
0104 int nBins = 30;
0105
0106 TH1F makeHisto(int nBins, float min, float max, string xAxisTitle, string unit, int precision)
0107 {
0108   TH1F myHisto(xAxisTitle.c_str(), xAxisTitle.c_str(), nBins, min, max);
0109
0110   if (unit != "") xAxisTitle += " [" + unit +  "]";
0111   myHisto.GetXaxis()->SetTitle(xAxisTitle.c_str());
0112
0113   float binWidth = (float) (max - min)/nBins;
0114   string yAxisTitle;
0115   yAxisTitle = "#pi^{-}/#pi^{+} ratio from K_{S}^{0} / " + to_string_with_precision(binWidth) + " " + unit;
0116   myHisto.GetYaxis()->SetTitle(yAxisTitle.c_str());
0117
0118   return myHisto;
0119 }
0120
0121 void pion_ratio_from_Kshort()
0122 {
0123   //Input file should be public
0124   string fileName = "/gpfs/mnt/gpfs02/sphenix/user/cdean/scripts/fitters/files/KShort6RunCombined.root";
0125   TFile *file = new TFile(fileName.c_str());
0126   TTree* data = (TTree*)file->Get("DecayTree");
0127
0128   float piminus_pT; data->SetBranchAddress("track_1_pT", &piminus_pT);
0129   float piplus_pT; data->SetBranchAddress("track_2_pT", &piplus_pT);
0130   float K_S0_mass; data->SetBranchAddress("K_S0_mass", &K_S0_mass);
0131
0132   TH1F piminus_pT_histo = makeHisto(nBins, min_pT, max_pT, "p_{T}", "GeV", 1);
0133   TH1F piplus_pT_histo = makeHisto(nBins, min_pT, max_pT, "p_{T}", "GeV", 1);
0134
0135   int tmp = 0;
0136   int barWidth = 50;
0137   int num_entries = data->GetEntries();
0138   for (int  l = 0; l < num_entries; ++l)
0139   {
0140     if (tmp != (int)100*l/num_entries)
0141     {
0142       tmp = (int)100*l/num_entries;
0143       if ((tmp%1)  == 0)
0144       {
0145         cout << "[";
0146         int pos = barWidth * tmp/100;
0147         for (int i = 0; i < barWidth; ++i)
0148         {
0149           if (i < pos) cout << "=";
0150           else if (i == pos) cout << ">";
0151           else cout << " ";
0152         }
0153         cout << "] " << tmp << " %\r";
0154         cout.flush();
0155       }
0156     }
0157
0158     data->GetEntry(l);
0159
0160     if (0.47 <= K_S0_mass && K_S0_mass <= 0.52)
0161     {
0162       piminus_pT_histo.Fill(piminus_pT);
0163       piplus_pT_histo.Fill(piplus_pT);
0164     }
0165   }
0166
0167   cout << "[";
0168   int pos = barWidth * tmp;
0169   for (int i = 0; i < barWidth; ++i)
0170   {
0171     if (i < pos) cout << "=";
0172     else if (i == pos) cout << ">";
0173     else cout << " ";
0174   }
0175   cout << "] 100 %\r";
0176   cout.flush();
0177   cout<<endl;
0178
0179   piminus_pT_histo.Sumw2();
0180   piplus_pT_histo.Sumw2();
0181   piminus_pT_histo.Divide(&piplus_pT_histo);
0182
0183   ratio_function = new TF1("fit", fitfunc, min_pT, max_pT, polynomial_order);
0184   ratio_function->SetLineColor(kRed);
0185   TFitResultPtr r = piminus_pT_histo.Fit(ratio_function, "RS");
0186
0187   savePlots(piminus_pT_histo, "Kshort_data_piminus_pT_to_piplus_pT");
0188
0189   file->Close();
0190 }