QM25_productions/D0/D0_QM25.C

0001 using namespace std;
0002
0003 string plotPath = "/sphenix/user/cdean/public/QM25";
0004
0005 string changeType[] = {"oldKFP_base_cut_ana475p017_optimized"};
0006 //string changeType[] = {"D0_eta", "D0_pT", "D0_decLenOverErr", "PV_z", "track_TPC_states", "track_MVTX_states", "track_pT", "track_eta", "track_INTT_hits", "D0_decLen", "D0_DIRA", "D0_chi2"};
0007
0008 class plotParameters
0009 {
0010   public:
0011     string branch;
0012     string title;
0013     string units;
0014     int bins;
0015     float min;
0016     float max;
0017     int precision;
0018
0019     plotParameters(string constructor_branch, string constructor_title, string constructor_units, int constructor_bins, float constructor_min, float constructor_max, int candidate_precision)
0020     {
0021       branch = constructor_branch;
0022       title = constructor_title;
0023       units = constructor_units;
0024       bins = constructor_bins;
0025       min = constructor_min;
0026       max = constructor_max;
0027       precision = candidate_precision;
0028     }
0029 };
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 template <typename T>
0041 void savePlots(T myPlot, string plotName, bool logY)
0042 {
0043   TCanvas *c1  = new TCanvas("myCanvas", "myCanvas",800,800);
0044
0045   myPlot.GetXaxis()->SetNdivisions(505);
0046   myPlot.GetYaxis()->SetNdivisions(505);
0047
0048   if (strncmp(typeid(myPlot).name(), "4TH2F", 5) == 0)
0049   {
0050     myPlot.Draw("COLZ");
0051   }
0052   else
0053   {
0054     myPlot.Sumw2();
0055     if (logY) gPad->SetLogy();
0056     myPlot.Draw("PE1");
0057   }
0058
0059   TPaveText *pt;
0060   pt = new TPaveText(0.25,0.94,0.99,1., "NDC");
0061   pt->SetFillColor(0);
0062   pt->SetFillStyle(0);
0063   pt->SetTextFont(42);
0064   TText *pt_LaTex = pt->AddText("#it{#bf{sPHENIX}} Internal, #sqrt{s} = 200 GeV, pp");
0065   pt->SetBorderSize(0);
0066   pt->Draw();
0067   gPad->Modified();
0068
0069   string extensions[] = {".png", ".C", ".pdf"};
0070   for (auto extension : extensions)
0071   {
0072     string output = plotPath + "/plots/" + plotName + extension;
0073     c1->SaveAs(output.c_str());
0074   }
0075 }
0076
0077 TH1F makeHisto(int nBins, float min, float max, string xAxisTitle, string unit, int precision)
0078 {
0079   TH1F myHisto(xAxisTitle.c_str(), xAxisTitle.c_str(), nBins, min, max);
0080
0081   if (unit != "") xAxisTitle += " [" + unit +  "]";
0082   myHisto.GetXaxis()->SetTitle(xAxisTitle.c_str());
0083
0084   float binWidth = (float) (max - min)/nBins;
0085   string yAxisTitle;
0086   if (unit != "") yAxisTitle = "Candidates / " + to_string_with_precision(binWidth, precision) + " " + unit;
0087   else yAxisTitle = "Candidates";
0088   myHisto.GetYaxis()->SetTitle(yAxisTitle.c_str());
0089
0090   return myHisto;
0091 }
0092
0093 TH2F makeHisto(int x_nBins, float x_min, float x_max, string xAxisTitle, string x_unit, int y_nBins, float y_min, float y_max, string yAxisTitle, string y_unit)
0094 {
0095   string combined_title = xAxisTitle + yAxisTitle;
0096   TH2F myHisto(combined_title.c_str(), combined_title.c_str(), x_nBins, x_min, x_max, y_nBins, y_min, y_max);
0097
0098   if (x_unit != "") xAxisTitle += " [" + x_unit +  "]";
0099   myHisto.GetXaxis()->SetTitle(xAxisTitle.c_str());
0100
0101   if (y_unit != "") yAxisTitle += " [" + y_unit +  "]";
0102   myHisto.GetYaxis()->SetTitle(yAxisTitle.c_str());
0103
0104   return myHisto;
0105 }
0106
0107 void D0_QM25()
0108 {
0109   int maxDigits = 3;
0110   TGaxis::SetMaxDigits(maxDigits);
0111   gStyle->SetOptStat(0);
0112   gStyle->SetPaintTextFormat(".3f");
0113
0114   TCanvas *c1  = new TCanvas("", "",800,800);
0115
0116   for (auto change : changeType)
0117   {
0118     string fileDir = "/sphenix/tg/tg01/hf/cdean/QM25_productions_ana475p017_20250401_2326/Kpi_reco/";
0119
0120     TChain *myChain = new TChain("DecayTree");
0121     TSystemDirectory dir("reco_files", fileDir.c_str());
0122     TList *files = dir.GetListOfFiles();
0123     TIter fileIterator(files);
0124     TSystemFile *file;
0125
0126     int maxFiles = -1;
0127     int currentLoadedFiles = 0;
0128     while ((file = (TSystemFile*)fileIterator()))
0129     {
0130       string fileName(file->GetName());
0131       if (fileName.find(".root") != string::npos)
0132       {
0133         if (fileName.find("53879") != string::npos) continue;
0134         ++currentLoadedFiles;
0135         if (currentLoadedFiles == maxFiles) break;
0136         string inputFile = fileDir + fileName;
0137         myChain->Add(inputFile.c_str());
0138       }
0139     }
0140     vector<TCut> individualCuts;
0141
0142     individualCuts.push_back("D0_pT > 0.75");
0143     individualCuts.push_back("D0_decayLength_xy > 0.008 && D0_decayLength_xy < 0.07");
0144     individualCuts.push_back("D0_decayLength_xy/D0_decayLengthErr_xy > 0.05");
0145     individualCuts.push_back("abs(D0_DIRA) > 0.985");
0146     individualCuts.push_back("abs(D0_chi2) < 14");
0147     individualCuts.push_back("D0_IP < 0.005");
0148     individualCuts.push_back("abs(primary_vertex_z) < 10.");
0149     individualCuts.push_back("min(track_1_TPC_nStates, track_2_TPC_nStates) > 25");
0150     individualCuts.push_back("min(track_1_INTT_nHits, track_2_INTT_nHits) > 1");
0151     individualCuts.push_back("min(track_1_MVTX_nStates, track_2_MVTX_nStates) == 3");
0152     individualCuts.push_back("0 <= track_1_bunch_crossing && track_1_bunch_crossing < 350");
0153     individualCuts.push_back("min(track_1_pT, track_2_pT) > 0.2");
0154     individualCuts.push_back("max(abs(track_1_IP_xy), abs(track_2_IP_xy)) < 0.02");
0155     individualCuts.push_back("track_1_track_2_DCA_xy < 0.002");
0156
0157     TCut total;
0158     for (auto cut : individualCuts)
0159     {
0160       total += cut;
0161     }
0162
0163     string outFileName = plotPath + "/output_D0_" + change + ".root";
0164     TFile *outFile = TFile::Open(outFileName.c_str(),"RECREATE");
0165     TTree* dataWithCut = myChain->CopyTree(total);
0166     outFile->Write();
0167
0168     int standard_bins = 50, standard_2D_bins = 150;
0169     vector<plotParameters> plots, twoD_plots;
0170
0171     twoD_plots.push_back(plotParameters("D0_mass", "m_{K^{+}#pi^{-}}", "GeV", standard_bins, 1.7, 2.0, 2));
0172     twoD_plots.push_back(plotParameters("D0_pseudorapidity", "D^{0} #eta", "", standard_bins, -1.2, 1.2, 2));
0173     twoD_plots.push_back(plotParameters("D0_pT", "D^{0} p_{T}", "GeV", standard_bins, 0.5, 4, 2));
0174     twoD_plots.push_back(plotParameters("D0_decayLength/D0_decayLengthErr", "D^{0}  Decay Length Over Error", "", standard_bins, 4, 10, 2));
0175     twoD_plots.push_back(plotParameters("D0_decayLength", "D^{0} decay length", "cm", standard_bins, 0, 0.06, 4));
0176     twoD_plots.push_back(plotParameters("D0_DIRA", "D^{0} DIRA", "", standard_bins, 0.9, 1, 2));
0177     twoD_plots.push_back(plotParameters("D0_chi2", "D^{0} SV #chi^{2}", "", standard_bins, 0, 25, 2));
0178     twoD_plots.push_back(plotParameters("D0_IP", "D^{0} IP", "cm", standard_bins, 0, 0.2, 4));
0179     twoD_plots.push_back(plotParameters("primary_vertex_z", "vtx_{z}", "cm", standard_bins, -13.5, 13.5, 2));
0180     twoD_plots.push_back(plotParameters("min(track_1_TPC_nStates, track_2_TPC_nStates)", "min. track TPC states", "", 31, 20, 50, 2));
0181     twoD_plots.push_back(plotParameters("min(track_1_INTT_nHits, track_2_INTT_nHits)", "min. track INTT clusters", "", 4, 0, 3, 2));
0182     twoD_plots.push_back(plotParameters("min(track_1_MVTX_nStates, track_2_MVTX_nStates)", "min. track MVTX states", "", 6, 0, 5, 2));
0183     twoD_plots.push_back(plotParameters("min(track_1_pT, track_2_pT)", "min. track p_{T}", "GeV", standard_bins, 0, 2, 2));
0184     twoD_plots.push_back(plotParameters("track_1_pseudorapidity", "track 1 #eta", "", standard_bins, -1.2, 1.2, 2));
0185     twoD_plots.push_back(plotParameters("track_2_pseudorapidity", "track 2 #eta", "", standard_bins, -1.2, 1.2, 2));
0186     twoD_plots.push_back(plotParameters("max(abs(track_1_IP_xy), abs(track_2_IP_xy))", "max. |track IP_{xy}|", "cm", standard_bins, 0, 0.15, 4));
0187     twoD_plots.push_back(plotParameters("track_1_track_2_DCA_xy", "track-to-track DCA", "cm", standard_bins, 0, 0.005, 4));
0188
0189     for (unsigned int i = 0; i < twoD_plots.size() - 1; ++i)
0190     {
0191       for (unsigned int j = i+1; j < twoD_plots.size(); ++j)
0192       {
0193         TH2F twoD_plot = makeHisto(standard_2D_bins, twoD_plots[i].min, twoD_plots[i].max, twoD_plots[i].title, twoD_plots[i].units,
0194                                    standard_2D_bins, twoD_plots[j].min, twoD_plots[j].max, twoD_plots[j].title, twoD_plots[j].units);
0195         string send_to_plot = twoD_plots[j].branch + ":" + twoD_plots[i].branch +  " >> " + twoD_plots[i].title + twoD_plots[j].title;
0196         dataWithCut->Draw(send_to_plot.c_str(), total);
0197         string saveName_first = twoD_plots[i].branch == "D0_decayLength/D0_decayLengthErr" ? "D0_decayLength_over_D0_decayLengthErr" : twoD_plots[i].branch;
0198         string saveName_second = twoD_plots[j].branch == "D0_decayLength/D0_decayLengthErr" ? "D0_decayLength_over_D0_decayLengthErr" : twoD_plots[j].branch;
0199         string saveName = saveName_first + "_" + saveName_second + "_altered_" + change;
0200         savePlots(twoD_plot, saveName, false);
0201       }
0202     }
0203
0204     for (auto& plot : twoD_plots)
0205     {
0206       TH1F histo = makeHisto(plot.bins, plot.min, plot.max, plot.title, plot.units, plot.precision);
0207       string send_to_plot =  plot.branch + " >> " + plot.title;
0208       dataWithCut->Draw(send_to_plot.c_str(), total);
0209       if (plot.branch == "D0_decayLength/D0_decayLengthErr") savePlots(histo, "D0_decayLength_over_D0_decayLengthErr", false);
0210       else savePlots(histo, plot.branch, false);
0211     }
0212   }
0213 }