NewCode2024/KshortMomentum/Kshort_sample.C

0001 #include "KshortMomentum.C"
0002
0003 // construct a Exponentially modified Gaussian distribution function
0004 double EMG(double *x, double *par)
0005 {
0006     double N = par[0];      // Normalization
0007     double mu = par[1];     // Mean of the Gaussian
0008     double sigma = par[2];  // Width of the Gaussian
0009     double lambda = par[3]; // Decay constant of the exponential
0010
0011     double t = x[0];
0012     double z = (mu + lambda * sigma * sigma - t) / (sqrt(2) * sigma);
0013
0014     double prefactor = lambda / 2.0;
0015     double exp_part = exp((lambda / 2.0) * (2.0 * mu + lambda * sigma * sigma - 2.0 * t));
0016     double erfc_part = TMath::Erfc(z);
0017
0018     return N * prefactor * exp_part * erfc_part;
0019 }
0020
0021 // Double Gaussian function: simplfied model, fraction is fixed to 0.5, gaussian widths are the same for both Gaussians
0022 double DBGaussian(double *x, double *par)
0023 {
0024     double N = par[0];     // Normalization
0025     double mu1 = par[1];   // Mean of the first Gaussian
0026     double mu2 = par[2];   // Mean of the second Gaussian
0027     double sigma = par[3]; // Width of the Gaussian
0028
0029     return N * (TMath::Gaus(x[0], mu1, sigma) + TMath::Gaus(x[0], mu2, sigma));
0030 }
0031
0032 void Kshort_sample()
0033 {
0034     ROOT::EnableImplicitMT();
0035     ROOT::RDataFrame df("HFTrackEfficiency", "/sphenix/user/cdean/public/forHaoRen/outputHFEffFile_Kshort.root");
0036
0037     auto h_motherP = df.Histo1D({"h_motherP", ";K_{S}^{0} p [GeV]", 100, 0, 10}, "true_mother_p");
0038     auto h_motherPt = df.Histo1D({"h_motherPt", ";K_{S}^{0} p_{T} [GeV]", 100, 0, 10}, "true_mother_pT");
0039     auto h_motherPz = df.Histo1D({"h_motherPz", ";K_{S}^{0} p_{z} [GeV]", 50, 0, 5}, "true_mother_pz");
0040     auto h_motherEta = df.Histo1D({"h_motherEta", ";K_{S}^{0} #eta", 40, -2, 2}, "true_mother_eta");
0041     auto h_mother_Pt_Eta = df.Histo2D({"h_mother_Pt_Eta", ";Truth K_{S}^{0} p_{T} [GeV];Truth K_{S}^{0} #eta", 100, 0, 4, 44, -1.1, 1.1}, "true_mother_pT", "true_mother_eta");
0042
0043     // TF1 for the Exponentially modified Gaussian for pt and Double Gaussian for eta
0044     TF1 *fpt = new TF1("fpt", EMG, -10, 10, 4);
0045     fpt->SetParameters(h_motherPt->GetMaximum(), 0.5, 1, 1);
0046     fpt->SetParName(0, "N");
0047     fpt->SetParName(1, "#mu");
0048     fpt->SetParName(2, "#sigma");
0049     fpt->SetParName(3, "#lambda");
0050     fpt->SetParLimits(0, 0, h_motherPt->GetMaximum() * 2);
0051     fpt->SetParLimits(1, 0, 2);
0052     fpt->SetParLimits(2, -100, 100);
0053     fpt->SetParLimits(3, -100, 100);
0054     fpt->SetLineColor(kRed);
0055     fpt->SetLineStyle(2);
0056     // Fit the histogram
0057     h_motherPt->Scale(1.0 / h_motherPt->Integral(-1, -1));
0058     h_motherPt->Fit(fpt, "R");
0059
0060     TF1 *feta = new TF1("feta", DBGaussian, -1, 1, 4);
0061     feta->SetParameters(1, -0.5, 0.5, 0.5);
0062     feta->SetParName(0, "N");
0063     feta->SetParName(1, "#mu_{1}");
0064     feta->SetParName(2, "#mu_{2}");
0065     feta->SetParName(3, "#sigma");
0066     feta->SetParLimits(0, 0, 1);
0067     feta->SetParLimits(1, -1, 0);
0068     feta->SetParLimits(2, 0, 1);
0069     feta->SetParLimits(3, 0, 10);
0070     // Fit the histogram
0071     h_motherEta->Scale(1.0 / h_motherEta->Integral(-1, -1));
0072     h_motherEta->Fit(feta, "R");
0073
0074     // Fit and draw the histograms
0075     TCanvas *c = new TCanvas("c", "c", 800, 700);
0076     c->cd();
0077     h_motherPt->GetYaxis()->SetTitle("Normalized counts");
0078     h_motherPt->GetYaxis()->SetTitleOffset(1.5);
0079     h_motherPt->Draw("hist");
0080     fpt->Draw("same");
0081     TLatex *pt = new TLatex();
0082     pt->SetNDC();
0083     pt->SetTextFont(42);
0084     pt->SetTextSize(0.04);
0085     pt->DrawLatex(0.4, 0.85, "Exponentially-modified Gaussian");
0086     pt->DrawLatex(0.4, 0.8, Form("#mu = %.3e#pm%.3e", fpt->GetParameter(1), fpt->GetParError(1)));
0087     pt->DrawLatex(0.4, 0.75, Form("#sigma = %.3e#pm%.3e", fpt->GetParameter(2), fpt->GetParError(2)));
0088     pt->DrawLatex(0.4, 0.7, Form("#lambda = %.3e#pm%.3e", fpt->GetParameter(3), fpt->GetParError(3)));
0089     pt->Draw();
0090     c->Draw();
0091     c->SaveAs("Kshort_motherPt_fit.png");
0092     c->SaveAs("Kshort_motherPt_fit.pdf");
0093
0094     c->Clear();
0095     c->cd();
0096     h_motherEta->GetYaxis()->SetRangeUser(0, h_motherEta->GetMaximum() * 1.65);
0097     h_motherEta->GetYaxis()->SetTitle("Normalized counts");
0098     h_motherEta->GetYaxis()->SetTitleOffset(1.5);
0099     h_motherEta->Draw("hist");
0100     feta->SetLineColor(kRed);
0101     feta->SetLineStyle(2);
0102     feta->Draw("same");
0103     pt->DrawLatex(0.55, 0.85, "Double Gaussian");
0104     pt->DrawLatex(0.55, 0.8, Form("#mu_{1} = %.3e#pm%.3e", feta->GetParameter(1), feta->GetParError(1)));
0105     pt->DrawLatex(0.55, 0.75, Form("#mu_{2} = %.3e#pm%.3e", feta->GetParameter(2), feta->GetParError(2)));
0106     pt->DrawLatex(0.55, 0.7, Form("#sigma = %.3e#pm%.3e", feta->GetParameter(3), feta->GetParError(3)));
0107     pt->Draw();
0108     c->Draw();
0109     c->SaveAs("Kshort_sampled_eta_fit.png");
0110     c->SaveAs("Kshort_sampled_eta_fit.pdf");
0111
0112     // sample pt and eta from the fitted functions
0113     // and calculate the corresponding kinematics
0114     int Nsample = h_mother_Pt_Eta->Integral(-1, -1, -1, -1);
0115     double M_kshort = 0.497611; // GeV
0116     TH1F *h_sampled_p = new TH1F("h_sampled_p", ";Sampled K_{S}^{0} p [GeV]", 100, 0, 10);
0117     TH1F *h_sampled_pz = new TH1F("h_sampled_pz", ";Sampled K_{S}^{0} p_{z} [GeV]", 50, 0, 5);
0118     TH2F *h_sampled_pt_eta = new TH2F("h_sampled_pt_eta", ";Sampled K_{S}^{0} p_{T} [GeV];Sampled K_{S}^{0} #eta", 100, 0, 4, 44, -1.1, 1.1);
0119     for (int i = 0; i < Nsample; i++)
0120     {
0121         double pt = fpt->GetRandom();
0122         double eta = feta->GetRandom();
0123         double phi = gRandom->Uniform(-TMath::Pi(), TMath::Pi());
0124         TLorentzVector lvec_rand;
0125         lvec_rand.SetPtEtaPhiM(pt, eta, phi, M_kshort);
0126         h_sampled_pt_eta->Fill(pt, eta);
0127         h_sampled_p->Fill(lvec_rand.P());
0128         h_sampled_pz->Fill(lvec_rand.Pz());
0129     }
0130     c->Clear();
0131     gPad->SetRightMargin(0.13);
0132     c->cd();
0133     h_sampled_pt_eta->Draw("colz");
0134     c->Draw();
0135     c->SaveAs("Kshort_sampled_pt_eta.png");
0136     c->SaveAs("Kshort_sampled_pt_eta.pdf");
0137
0138     c->Clear();
0139     gPad->SetRightMargin(0.13);
0140     c->cd();
0141     h_mother_Pt_Eta->Draw("colz");
0142     c->Draw();
0143     c->SaveAs("Kshort_truth_pt_eta.png");
0144     c->SaveAs("Kshort_truth_pt_eta.pdf");
0145
0146     // Draw the sampled histograms
0147     Draw1D(h_motherP.GetPtr(), h_sampled_p, "Kshort_motherP_TF1sampled", "Kshort p", "K_{S}^{0} p [GeV]", "Normalized counts");
0148     Draw1D(h_motherPz.GetPtr(), h_sampled_pz, "Kshort_motherPz_TF1sampled", "Kshort pz", "K_{S}^{0} p_{z} [GeV]", "Normalized counts");
0149 }