GeoOffset/GetResult/FromAvgVtxXY.C

0001 #include <iostream>
0002 #include <fstream>
0003 #include <filesystem>
0004 #include <cstdlib>
0005 #include <numeric>
0006 using namespace std;
0007
0008 #include "sPhenixStyle.C"
0009
0010 vector<string> read_list(string folder_direction, string MC_list_name)
0011 {
0012     vector<string> file_list;
0013     string list_buffer;
0014     ifstream data_list;
0015     data_list.open((folder_direction + "/" + MC_list_name).c_str());
0016
0017    file_list.clear();
0018
0019     while (1)
0020     {
0021         data_list >> list_buffer;
0022         if (!data_list.good())
0023         {
0024             break;
0025         }
0026         file_list.push_back(list_buffer);
0027     }
0028     cout<<"size in the" <<MC_list_name<<": "<<file_list.size()<<endl;
0029
0030     return file_list;
0031 }
0032
0033 int FromAvgVtxXY()
0034 {
0035
0036     std::string input_directory = "/sphenix/user/ChengWei/sPH_dNdeta/Run24AuAuMC/Sim_HIJING_MDC2_ana472_20250307/GeoOffset/completed";
0037     std::string input_foldername_NoIndex = "Run_00"; // note : Run_00XXX
0038     std::string filename_NoIndex = "MC_AvgVtxXY_GeoOffset_test1";
0039     std::string label_text = "Simulation";
0040
0041     // std::pair<double,double> ideal_recoXY = {-0.02239, 0.22385}; // -0.02239 [cm], 0.22385 [cm] //note: deltaphi = 0.04
0042     std::pair<double,double> ideal_recoXY = {-0.0223617, 0.223504}; // -0.0223617 0.223504 //note: deltaphi = 0.026
0043
0044     double frame_shift_forX = ideal_recoXY.first; // note : cm
0045     double frame_shift_forY = ideal_recoXY.second; // note : cm
0046
0047     // Division : ------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0048
0049     input_foldername_NoIndex = (input_foldername_NoIndex.back() == '_') ? input_foldername_NoIndex.substr(0, input_foldername_NoIndex.size() - 1) : input_foldername_NoIndex;
0050     filename_NoIndex = (filename_NoIndex.back() == '_') ? filename_NoIndex.substr(0, filename_NoIndex.size() - 1) : filename_NoIndex;
0051
0052     std::string file_list_name = "file_list_AvgVtxXY.txt";
0053     std::string output_directory = input_directory + "/" + "merged_result";
0054
0055     double X_range_l = (frame_shift_forX - 0.15 );
0056     double X_range_r = (frame_shift_forX + 0.15 );
0057     double Y_range_l = (frame_shift_forY - 0.15 );
0058     double Y_range_r = (frame_shift_forY + 0.15 );
0059
0060     // Division : ------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0061
0062     // note : to generate the merged_result folder
0063     if(std::filesystem::exists(Form("%s",output_directory.c_str())) == false){
0064         system(Form("mkdir -p %s", output_directory.c_str()));
0065         cout<<"----------- check folder exists -----------"<<endl;
0066         system(Form("ls %s", output_directory.c_str()));
0067     }
0068     else
0069     {
0070         cout<<"----------- folder exists already -----------"<<endl;
0071         system(Form("ls %s", output_directory.c_str()));
0072     }
0073
0074     // note : to generate the file_list.txt
0075     // todo: the file_name
0076     if(true/*std::filesystem::exists(Form("%s/%s",input_directory.c_str(),file_list_name.c_str())) == false*/){
0077         system(Form("ls %s/%s*/%s_*.root > %s/%s", input_directory.c_str(), input_foldername_NoIndex.c_str(), filename_NoIndex.c_str(), input_directory.c_str(), file_list_name.c_str()));
0078         cout<<"----------- file list generated -----------"<<endl;
0079         system(Form("cat %s/%s", input_directory.c_str(), file_list_name.c_str()));
0080         cout<<"----------- file list generated -----------"<<endl;
0081     }
0082     else
0083     {
0084         cout<<"----------- found the file list -----------"<<endl;
0085         system(Form("cat %s/%s", input_directory.c_str(), file_list_name.c_str()));
0086         cout<<"----------- found the file list -----------"<<endl;
0087     }
0088
0089     vector<string> file_list = read_list(input_directory, file_list_name);
0090
0091     // Division : ------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0092     SetsPhenixStyle();
0093     TCanvas * c1 = new TCanvas("c1", "c1", 950, 800);
0094     TGraph * g_vtxXY_ideal = new TGraph();
0095     g_vtxXY_ideal -> SetPoint(0, ideal_recoXY.first, ideal_recoXY.second);
0096     g_vtxXY_ideal -> SetMarkerStyle(28);
0097     g_vtxXY_ideal -> SetMarkerSize(1.5);
0098     g_vtxXY_ideal -> SetLineWidth(2);
0099     g_vtxXY_ideal -> SetMarkerColor(kRed);
0100
0101     TLatex * ltx = new TLatex();
0102     ltx->SetNDC();
0103     ltx->SetTextSize(0.045);
0104     ltx->SetTextAlign(31);
0105
0106     // Division : ------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0107     TH1D * h1D_vtxX = new TH1D("h1D_vtxX", "h1D_vtxX;Reco. Vtx X [cm];Entries", 100, X_range_l, X_range_r);
0108     TH1D * h1D_vtxY = new TH1D("h1D_vtxY", "h1D_vtxY;Reco. Vtx Y [cm];Entries", 100, Y_range_l, Y_range_r);
0109     TH2D * h2D_vtxXY = new TH2D("h2D_vtxXY", "h2D_vtxXY;Reco. Vtx X [cm];Reco. Vtx Y [cm]", 100, X_range_l, X_range_r, 100, Y_range_l, Y_range_r);
0110
0111     TH1D * h1D_GeoVaryX = new TH1D("h1D_GeoVaryX", "h1D_GeoVaryX;Ladder Variation in X [cm];Entries", 100, -0.03, 0.03);
0112     TH1D * h1D_GeoVaryY = new TH1D("h1D_GeoVaryY", "h1D_GeoVaryY;Ladder Variation in Y [cm];Entries", 100, -0.03, 0.03);
0113     TH1D * h1D_GeoVaryZ = new TH1D("h1D_GeoVaryZ", "h1D_GeoVaryZ;Ladder Variation in Z [cm];Entries", 100, -0.03, 0.03);
0114
0115
0116     // Division : ------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0117
0118     TChain *chain_vtxXY_in = new TChain("tree_vtxXY");
0119     for (string file_name : file_list) { chain_vtxXY_in -> Add((file_name).c_str()); }
0120     cout<<"chain_vtxXY_in -> GetEntries() : "<<chain_vtxXY_in -> GetEntries()<<endl;
0121
0122     double in_vtxX_quadrant;
0123     double in_vtxXE_quadrant;
0124     double in_vtxY_quadrant;
0125     double in_vtxYE_quadrant;
0126     double in_vtxX_LineFill;
0127     double in_vtxXE_LineFill;
0128     double in_vtxXStdDev_LineFill;
0129     double in_vtxY_LineFill;
0130     double in_vtxYE_LineFill;
0131     double in_vtxYStdDev_LineFill;
0132
0133     chain_vtxXY_in -> SetBranchAddress("vtxX_quadrant", &in_vtxX_quadrant);
0134     chain_vtxXY_in -> SetBranchAddress("vtxXE_quadrant", &in_vtxXE_quadrant);
0135     chain_vtxXY_in -> SetBranchAddress("vtxY_quadrant", &in_vtxY_quadrant);
0136     chain_vtxXY_in -> SetBranchAddress("vtxYE_quadrant", &in_vtxYE_quadrant);
0137     chain_vtxXY_in -> SetBranchAddress("vtxX_LineFill", &in_vtxX_LineFill);
0138     chain_vtxXY_in -> SetBranchAddress("vtxXE_LineFill", &in_vtxXE_LineFill);
0139     chain_vtxXY_in -> SetBranchAddress("vtxXStdDev_LineFill", &in_vtxXStdDev_LineFill);
0140     chain_vtxXY_in -> SetBranchAddress("vtxY_LineFill", &in_vtxY_LineFill);
0141     chain_vtxXY_in -> SetBranchAddress("vtxYE_LineFill", &in_vtxYE_LineFill);
0142     chain_vtxXY_in -> SetBranchAddress("vtxYStdDev_LineFill", &in_vtxYStdDev_LineFill);
0143
0144     for (int i = 0; i < chain_vtxXY_in -> GetEntries(); i++)
0145     {
0146         chain_vtxXY_in -> GetEntry(i);
0147
0148         double final_X = (in_vtxX_quadrant + in_vtxX_LineFill) / 2.;
0149         double final_Y = (in_vtxY_quadrant + in_vtxY_LineFill) / 2.;
0150
0151         h1D_vtxX -> Fill(final_X);
0152         h1D_vtxY -> Fill(final_Y);
0153         h2D_vtxXY -> Fill(final_X, final_Y);
0154     }
0155
0156     TChain *chain_GeoOffset_in = new TChain("tree_geooffset");
0157     for (string file_name : file_list) { chain_GeoOffset_in -> Add((file_name).c_str()); }
0158     cout<<"chain_GeoOffset_in -> GetEntries() : "<<chain_GeoOffset_in -> GetEntries()<<endl;
0159
0160     double offset_x;
0161     double offset_y;
0162     double offset_z;
0163
0164     chain_GeoOffset_in -> SetBranchAddress("offset_x", &offset_x);
0165     chain_GeoOffset_in -> SetBranchAddress("offset_y", &offset_y);
0166     chain_GeoOffset_in -> SetBranchAddress("offset_z", &offset_z);
0167
0168
0169     for (int geo_i = 0; geo_i < chain_GeoOffset_in -> GetEntries(); geo_i++)
0170     {
0171         chain_GeoOffset_in -> GetEntry(geo_i);
0172
0173         h1D_GeoVaryX -> Fill(offset_x);
0174         h1D_GeoVaryY -> Fill(offset_y);
0175         h1D_GeoVaryZ -> Fill(offset_z);
0176     }
0177
0178     TFile * file_out = new TFile(Form("%s/AvgXY_GeoOffset.root", output_directory.c_str()), "RECREATE");
0179     h1D_vtxX -> Write();
0180     h1D_vtxY -> Write();
0181     h2D_vtxXY -> Write();
0182
0183     h1D_GeoVaryX -> Write();
0184     h1D_GeoVaryY -> Write();
0185     h1D_GeoVaryZ -> Write();
0186
0187     c1 -> cd();
0188     h2D_vtxXY -> Draw("colz0");
0189     g_vtxXY_ideal -> Draw("P same");
0190     c1 -> Write("c1_h2D_vtxXY");
0191     c1 -> Clear();
0192
0193     c1 -> cd();
0194     h2D_vtxXY -> Draw("colz0");
0195     g_vtxXY_ideal -> Draw("P same");
0196     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", label_text.c_str()));
0197     c1 -> Print(Form("%s/h2D_AvgXY_GeoOffset.pdf", output_directory.c_str()));
0198     c1 -> Clear();
0199
0200     c1 -> cd();
0201     h1D_vtxX -> Draw("hist");
0202     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", label_text.c_str()));
0203     c1 -> Print(Form("%s/h1D_AvgX_GeoOffset.pdf", output_directory.c_str()));
0204     c1 -> Clear();
0205
0206
0207     c1 -> cd();
0208     h1D_vtxY -> Draw("hist");
0209     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", label_text.c_str()));
0210     c1 -> Print(Form("%s/h1D_AvgY_GeoOffset.pdf", output_directory.c_str()));
0211     c1 -> Clear();
0212
0213
0214     c1 -> cd();
0215     h1D_GeoVaryX -> Draw("hist");
0216     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", label_text.c_str()));
0217     c1 -> Print(Form("%s/h1D_GeoVaryX.pdf", output_directory.c_str()));
0218     c1 -> Clear();
0219
0220     h1D_GeoVaryY -> Draw("hist");
0221     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", label_text.c_str()));
0222     c1 -> Print(Form("%s/h1D_GeoVaryY.pdf", output_directory.c_str()));
0223     c1 -> Clear();
0224
0225     h1D_GeoVaryZ -> Draw("hist");
0226     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", label_text.c_str()));
0227     c1 -> Print(Form("%s/h1D_GeoVaryZ.pdf", output_directory.c_str()));
0228     c1 -> Clear();
0229
0230
0231
0232     file_out -> Close();
0233     return 0;
0234 }