analysis_INTT/plot/INTTVtxZ_Sim.C

0001 #include <algorithm>
0002 #include <cmath>
0003 #include <iostream>
0004 #include <type_traits>
0005 #include <vector>
0006
0007 // Template function to calculate percentile bounds in descending order
0008 template <typename T>
0009 std::vector<std::pair<T, T>> calculatePercentileBounds(const std::vector<T>& data, int interval) {
0010     // Ensure T is a numeric type
0011     static_assert(std::is_arithmetic<T>::value, "T must be a numeric type.");
0012
0013     if (data.empty()) {
0014         throw std::invalid_argument("Data vector is empty.");
0015     }
0016
0017     if (interval <= 0 || interval > 100) {
0018         throw std::invalid_argument("Interval must be between 1 and 100.");
0019     }
0020
0021     // Create a copy of the input vector and sort it in descending order
0022     std::vector<T> sortedData = data;
0023     std::sort(sortedData.begin(), sortedData.end(), std::greater<T>());
0024
0025     std::vector<std::pair<T, T>> percentileBounds;
0026     size_t n = sortedData.size();
0027
0028     // Calculate percentile bounds based on the interval
0029     for (int p = 0; p < 100; p += interval) {
0030         double lowerRank = p / 100.0 * (n - 1);             // Lower bound rank
0031         double upperRank = (p + interval) / 100.0 * (n - 1); // Upper bound rank
0032
0033         size_t lowerIdx = static_cast<size_t>(std::floor(lowerRank)); // Floor index for lower bound
0034         size_t upperIdx = static_cast<size_t>(std::ceil(upperRank));  // Ceil index for upper bound
0035
0036         // Interpolate lower bound
0037         T lowerValue = sortedData[lowerIdx];
0038         if (lowerIdx < n - 1 && lowerIdx != upperRank) {
0039             double lowerWeight = lowerRank - lowerIdx;
0040             lowerValue = static_cast<T>(
0041                 sortedData[lowerIdx] * (1 - lowerWeight) + sortedData[lowerIdx + 1] * lowerWeight
0042             );
0043         }
0044
0045         // Interpolate upper bound
0046         T upperValue;
0047         if (p + interval < 100) {
0048             upperValue = sortedData[upperIdx];
0049             if (upperIdx > 0 && upperIdx != upperRank) {
0050                 double upperWeight = upperRank - upperIdx;
0051                 upperValue = static_cast<T>(
0052                     sortedData[upperIdx - 1] * (1 - upperWeight) + sortedData[upperIdx] * upperWeight
0053                 );
0054             }
0055         } else {
0056             // For the last interval, set the upper bound to exceed the max value
0057             upperValue = sortedData[n - 1] - 1; // Subtract 1 or adjust as needed for descending
0058         }
0059
0060         percentileBounds.emplace_back(lowerValue, upperValue);
0061     }
0062
0063     return percentileBounds;
0064 }
0065
0066 // make a function to draw 2D histograms
0067 void draw2DHist(TH2 *hist, bool logz, const std::string exttext, const std::string &filename)
0068 {
0069     TCanvas *c = new TCanvas("canvas", "canvas", 800, 700);
0070     c->cd();
0071     c->SetLogz(logz);
0072     gPad->SetRightMargin(0.2);
0073     hist->GetZaxis()->SetTitleOffset(1.5);
0074     hist->Draw("colz");
0075     // Add exttext to the plot
0076     TLatex text;
0077     text.SetNDC();
0078     text.SetTextFont(43);
0079     text.SetTextSize(25);
0080     text.DrawLatex(0.2, 0.85, exttext.c_str());
0081     // Save the plot
0082     c->SaveAs(Form("%s.pdf", filename.c_str()));
0083     c->SaveAs(Form("%s.png", filename.c_str()));
0084
0085     //properly delete the canvas and memory
0086     delete c;
0087 }
0088
0089 void INTTVtxZ_Sim()
0090 {
0091     int Npercentile = 10;
0092
0093     std::string fileprefix = "Sim_HIJING_ananew_20250131";
0094     std::string plotdir = "./RecoPV_ana/" + fileprefix;
0095     system(Form("mkdir -p %s", plotdir.c_str()));
0096
0097     ROOT::EnableImplicitMT();
0098     ROOT::RDataFrame df("minitree", Form("/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/minitree/VtxEvtMap_%s/minitree_00*.root", fileprefix.c_str()));
0099     // get NClusLayer1 in each event in std::vector<int>
0100     auto NClusLayer1 = df.Take<int>("NClusLayer1");
0101     // print out the min, max of NClusLayer1
0102     std::cout << "NClusLayer1: min = " << *std::min_element(NClusLayer1.begin(), NClusLayer1.end()) << ", max = " << *std::max_element(NClusLayer1.begin(), NClusLayer1.end()) << std::endl;
0103     // get the percentile intervals
0104     auto percentiles = calculatePercentileBounds(*NClusLayer1, Npercentile);
0105     std::cout << "Percentiles (NClusLayer1): ";
0106     for (const auto &p : percentiles)
0107     {
0108         std::cout << "[" << p.first << ", " << p.second << "] ";
0109     }
0110     std::cout << std::endl;
0111
0112     // Add a new column to the dataframe: the difference between the reconstructed and truth vertex Z PV_z - TruthPV_z
0113     auto df_with_diff = df.Define("RecoTruthDiff", "PV_z - TruthPV_z");
0114     // Add a new column to label the events based on the percentile intervals
0115     auto df_with_percentile = df_with_diff.Define("NClusLayer1Percentile", [percentiles](int NClusLayer1) {
0116         for (size_t i = 0; i < percentiles.size(); i++)
0117         {
0118             if (NClusLayer1 <= percentiles[i].first && NClusLayer1 > percentiles[i].second)
0119             {
0120                 return static_cast<int>(i);
0121             }
0122         }
0123         return -1;
0124     }, {"NClusLayer1"});
0125
0126     // histogram for reconstructed and truth vertex Z
0127     auto hM_TruthVtxZ_RecoTruthDiff = df_with_percentile.Filter("is_min_bias").Histo2D({"hM_TruthVtxZ_RecoTruthDiff", ";vtx_{Z}^{Truth} [cm];vtx_{Z}^{Reco}-vtx_{Z}^{Truth} [cm];Entries", 200, -25, 25, 200, -5, 5}, "TruthPV_z", "RecoTruthDiff");
0128     draw2DHist(hM_TruthVtxZ_RecoTruthDiff.GetPtr(), true, "", Form("%s/TruthVtxZ_RecoTruthDiff_Sim", plotdir.c_str()));
0129
0130     // create histograms for each percentile interval and draw them
0131     for (size_t i = 0; i < percentiles.size(); i++)
0132     {
0133         auto hM_TruthVtxZ_RecoTruthDiff_Percentile = df_with_percentile.Filter("is_min_bias").Filter([i](int NClusLayer1Percentile) { return NClusLayer1Percentile == i; }, {"NClusLayer1Percentile"}).Histo2D({"hM_TruthVtxZ_RecoTruthDiff_Percentile", ";vtx_{Z}^{Truth} [cm];vtx_{Z}^{Reco}-vtx_{Z}^{Truth} [cm];Entries", 200, -25, 25, 200, -5, 5}, "TruthPV_z", "RecoTruthDiff");
0134         draw2DHist(hM_TruthVtxZ_RecoTruthDiff_Percentile.GetPtr(), true, Form("Cluster multiplicity percentile %d-%d%%", static_cast<int>(i*100/Npercentile), static_cast<int>((i+1)*100/Npercentile)), Form("%s/TruthVtxZ_RecoTruthDiff_Percentile%d_Sim", plotdir.c_str(), static_cast<int>(i)));
0135     }
0136 }