GeoOffset/GetResult/FromdNdEta.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 void FindH1DUpperLowBound(TH1D * major_hist, TH1D * new_hist, bool IsForUpper)
0034 {
0035     for (int i = 1; i <= major_hist -> GetNbinsX(); i++)
0036     {
0037         if (IsForUpper)
0038         {
0039             if (new_hist -> GetBinContent(i) > major_hist -> GetBinContent(i))
0040             {
0041                 major_hist -> SetBinContent(i, new_hist -> GetBinContent(i));
0042             }
0043         }
0044         else
0045         {
0046             if (new_hist -> GetBinContent(i) < major_hist -> GetBinContent(i))
0047             {
0048                 major_hist -> SetBinContent(i, new_hist -> GetBinContent(i));
0049             }
0050         }
0051     }
0052 }
0053
0054 double vector_stddev (vector <double> input_vector){
0055
0056     double sum_subt = 0;
0057     double average  = accumulate( input_vector.begin(), input_vector.end(), 0.0 ) / double(input_vector.size());
0058
0059     // cout<<"average is : "<<average<<endl;
0060
0061     for (int i=0; i<input_vector.size(); i++){ sum_subt += pow((input_vector[i] - average),2); }
0062
0063     //cout<<"sum_subt : "<<sum_subt<<endl;
0064     // cout<<"print from the function, average : "<<average<<" std : "<<stddev<<endl;
0065
0066     return sqrt( sum_subt / double(input_vector.size()-1) );
0067 }
0068
0069 int FromdNdEta()
0070 {
0071
0072     std::string input_directory = "/sphenix/user/ChengWei/sPH_dNdeta/Run24AuAuMC/Sim_HIJING_MDC2_ana472_20250307/GeoOffset/completed";
0073     std::string input_foldername_NoIndex = "Run_00"; // note : Run_00XXX
0074     std::string filename_NoIndex = "MC_PreparedNdEtaEach_AllSensor_VtxZ10_Mbin70_test1"; // note : take out the _00001_dNdEta.root
0075
0076     std::string file_directory_ideal = "/sphenix/user/ChengWei/sPH_dNdeta/Run24AuAuMC/Sim_HIJING_MDC2_ana472_20250307/GeoOffset/completed/Run_NoGeoOffset_00000/MC_PreparedNdEtaEach_AllSensor_VtxZ10_Mbin70_test1_00000_dNdEta.root";
0077
0078     std::string label_text = "Simulation";
0079
0080     // Division : ------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0081
0082     input_foldername_NoIndex = (input_foldername_NoIndex.back() == '_') ? input_foldername_NoIndex.substr(0, input_foldername_NoIndex.size() - 1) : input_foldername_NoIndex;
0083     filename_NoIndex = (filename_NoIndex.back() == '_') ? filename_NoIndex.substr(0, filename_NoIndex.size() - 1) : filename_NoIndex;
0084
0085     std::string file_list_name = "file_list_dNdEta.txt";
0086     std::string output_directory = input_directory + "/" + "merged_result";
0087
0088     // Division : ------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0089
0090     // note : to generate the merged_result folder
0091     if(std::filesystem::exists(Form("%s",output_directory.c_str())) == false){
0092         system(Form("mkdir -p %s", output_directory.c_str()));
0093         cout<<"----------- check folder exists -----------"<<endl;
0094         system(Form("ls %s", output_directory.c_str()));
0095     }
0096     else
0097     {
0098         cout<<"----------- folder exists already -----------"<<endl;
0099         system(Form("ls %s", output_directory.c_str()));
0100     }
0101
0102     // note : to generate the file_list.txt
0103     // todo: the file_name
0104     if(true/*std::filesystem::exists(Form("%s/%s",input_directory.c_str(),file_list_name.c_str())) == false*/){
0105         system(Form("ls %s/%s*/%s_*_dNdEta.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()));
0106         cout<<"----------- file list generated -----------"<<endl;
0107         system(Form("cat %s/%s", input_directory.c_str(), file_list_name.c_str()));
0108         cout<<"----------- file list generated -----------"<<endl;
0109     }
0110     else
0111     {
0112         cout<<"----------- found the file list -----------"<<endl;
0113         system(Form("cat %s/%s", input_directory.c_str(), file_list_name.c_str()));
0114         cout<<"----------- found the file list -----------"<<endl;
0115     }
0116
0117     vector<string> file_list = read_list(input_directory, file_list_name);
0118
0119     // Division : ------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0120
0121     SetsPhenixStyle();
0122     TCanvas * c1 = new TCanvas("c1", "c1", 950, 800);
0123     TFile * file_in_ideal = TFile::Open(Form("%s", file_directory_ideal.c_str()));
0124     TH1D * h1D_BestPair_RecoTrackletEtaPerEvt_ideal = (TH1D*)file_in_ideal->Get("h1D_BestPair_RecoTrackletEtaPerEvt");
0125     h1D_BestPair_RecoTrackletEtaPerEvt_ideal -> SetLineColor(2);
0126     // h1D_BestPair_RecoTrackletEtaPerEvt_ideal -> Scale(1. / h1D_BestPair_RecoTrackletEtaPerEvt_ideal -> Integral(-1,-1));
0127
0128     TH1D * h1D_RotatedBkg_RecoTrackletEtaPerEvt_ideal = (TH1D*)file_in_ideal->Get("h1D_RotatedBkg_RecoTrackletEtaPerEvt");
0129     h1D_RotatedBkg_RecoTrackletEtaPerEvt_ideal -> SetLineColor(2);
0130     // h1D_RotatedBkg_RecoTrackletEtaPerEvt_ideal -> Scale(1. / h1D_RotatedBkg_RecoTrackletEtaPerEvt_ideal -> Integral(-1,-1));
0131
0132     std::map<int,std::vector<double>> EtaBin_Variation_vec; EtaBin_Variation_vec.clear();
0133     for (int i = 0; i < h1D_RotatedBkg_RecoTrackletEtaPerEvt_ideal -> GetNbinsX(); i++)
0134     {
0135         EtaBin_Variation_vec[i] = {};
0136     }
0137
0138     TLatex * ltx = new TLatex();
0139     ltx->SetNDC();
0140     ltx->SetTextSize(0.045);
0141     ltx->SetTextAlign(31);
0142
0143     // Division : ------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0144
0145     TH1D * h1D_RotatedBkg_RecoTrackletEtaPerEvt_LowerBound = (TH1D*) h1D_RotatedBkg_RecoTrackletEtaPerEvt_ideal -> Clone("h1D_RotatedBkg_RecoTrackletEtaPerEvt_LowerBound");
0146     TH1D * h1D_RotatedBkg_RecoTrackletEtaPerEvt_UpperBound = (TH1D*) h1D_RotatedBkg_RecoTrackletEtaPerEvt_ideal -> Clone("h1D_RotatedBkg_RecoTrackletEtaPerEvt_UpperBound");
0147
0148     // Division : ------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0149
0150     TFile * file_in;
0151     TH2D * h2D_BestPair_RecoTrackletEtaPerEvt;
0152     TH2D * h2D_RotatedBkg_RecoTrackletEtaPerEvt;
0153
0154     TH1D * h1D_BestPair_RecoTrackletEtaPerEvt;
0155     TH1D * h1D_RotatedBkg_RecoTrackletEtaPerEvt;
0156
0157     for (int i = 0; i < file_list.size(); i++)
0158     {
0159         file_in = TFile::Open(Form("%s", file_list[i].c_str()));
0160
0161         h1D_BestPair_RecoTrackletEtaPerEvt = (TH1D*) file_in -> Get("h1D_BestPair_RecoTrackletEtaPerEvt");
0162         h1D_RotatedBkg_RecoTrackletEtaPerEvt = (TH1D*) file_in -> Get("h1D_RotatedBkg_RecoTrackletEtaPerEvt");
0163
0164         // h1D_BestPair_RecoTrackletEtaPerEvt -> Scale(1. / h1D_BestPair_RecoTrackletEtaPerEvt -> Integral(-1, -1));
0165         // h1D_RotatedBkg_RecoTrackletEtaPerEvt -> Scale(1. / h1D_RotatedBkg_RecoTrackletEtaPerEvt -> Integral(-1, -1));
0166
0167         FindH1DUpperLowBound(h1D_RotatedBkg_RecoTrackletEtaPerEvt_LowerBound, h1D_RotatedBkg_RecoTrackletEtaPerEvt, false);
0168         FindH1DUpperLowBound(h1D_RotatedBkg_RecoTrackletEtaPerEvt_UpperBound, h1D_RotatedBkg_RecoTrackletEtaPerEvt, true);
0169
0170         for (int eta_bin = 0; eta_bin < h1D_RotatedBkg_RecoTrackletEtaPerEvt -> GetNbinsX(); eta_bin++)
0171         {
0172             EtaBin_Variation_vec[eta_bin].push_back(
0173                 (h1D_RotatedBkg_RecoTrackletEtaPerEvt -> GetBinContent(eta_bin + 1) - h1D_RotatedBkg_RecoTrackletEtaPerEvt_ideal -> GetBinContent(eta_bin + 1)) / h1D_RotatedBkg_RecoTrackletEtaPerEvt_ideal -> GetBinContent(eta_bin + 1)
0174             );
0175         }
0176
0177         if (i == 0)
0178         {
0179             h2D_BestPair_RecoTrackletEtaPerEvt = new TH2D(
0180                 "h2D_BestPair_RecoTrackletEtaPerEvt",
0181                 "h2D_BestPair_RecoTrackletEtaPerEvt;Tracklet #eta; N Reco. Tracklets per event",
0182                 h1D_BestPair_RecoTrackletEtaPerEvt -> GetNbinsX(), h1D_BestPair_RecoTrackletEtaPerEvt -> GetXaxis() -> GetXmin(), h1D_BestPair_RecoTrackletEtaPerEvt -> GetXaxis() -> GetXmax(),
0183                 200, h1D_BestPair_RecoTrackletEtaPerEvt -> GetMinimum(), h1D_BestPair_RecoTrackletEtaPerEvt -> GetBinContent(h1D_BestPair_RecoTrackletEtaPerEvt -> GetMaximumBin()) * 1.2
0184             );
0185
0186             h2D_RotatedBkg_RecoTrackletEtaPerEvt = new TH2D(
0187                 "h2D_RotatedBkg_RecoTrackletEtaPerEvt",
0188                 "h2D_RotatedBkg_RecoTrackletEtaPerEvt;Tracklet #eta; N Reco. Tracklets per event",
0189                 h1D_RotatedBkg_RecoTrackletEtaPerEvt -> GetNbinsX(), h1D_RotatedBkg_RecoTrackletEtaPerEvt -> GetXaxis() -> GetXmin(), h1D_RotatedBkg_RecoTrackletEtaPerEvt -> GetXaxis() -> GetXmax(),
0190                 200, h1D_RotatedBkg_RecoTrackletEtaPerEvt -> GetMinimum(), h1D_RotatedBkg_RecoTrackletEtaPerEvt -> GetBinContent(h1D_RotatedBkg_RecoTrackletEtaPerEvt -> GetMaximumBin()) * 1.2
0191             );
0192         }
0193
0194         for (int bin_i = 1; bin_i <= h1D_BestPair_RecoTrackletEtaPerEvt -> GetNbinsX(); bin_i++)
0195         {
0196             h2D_BestPair_RecoTrackletEtaPerEvt -> Fill(
0197                 h1D_BestPair_RecoTrackletEtaPerEvt -> GetBinCenter(bin_i),
0198                 h1D_BestPair_RecoTrackletEtaPerEvt -> GetBinContent(bin_i)
0199             );
0200         }
0201
0202         for (int bin_i = 1; bin_i <= h1D_RotatedBkg_RecoTrackletEtaPerEvt -> GetNbinsX(); bin_i++)
0203         {
0204             h2D_RotatedBkg_RecoTrackletEtaPerEvt -> Fill(
0205                 h1D_RotatedBkg_RecoTrackletEtaPerEvt -> GetBinCenter(bin_i),
0206                 h1D_RotatedBkg_RecoTrackletEtaPerEvt -> GetBinContent(bin_i)
0207             );
0208         }
0209     }
0210
0211     TH1D * h1D_RotatedBkg_RecoTrackletEtaPerEvt_VariationSigma = (TH1D*) h1D_RotatedBkg_RecoTrackletEtaPerEvt_ideal -> Clone("h1D_RotatedBkg_RecoTrackletEtaPerEvt_VariationSigma");
0212     h1D_RotatedBkg_RecoTrackletEtaPerEvt_VariationSigma -> Reset("ICESM");
0213     for (int i = 0; i < h1D_RotatedBkg_RecoTrackletEtaPerEvt_VariationSigma -> GetNbinsX(); i++)
0214     {
0215         double variation_stddev = vector_stddev(EtaBin_Variation_vec[i]);
0216         std::cout<<"EtaBin: "<<i+1<<" std: "<<variation_stddev<<std::endl;
0217
0218         if (variation_stddev != variation_stddev){continue;}// note : skip the nan value
0219
0220         h1D_RotatedBkg_RecoTrackletEtaPerEvt_VariationSigma -> SetBinContent(i + 1, 3 * vector_stddev(EtaBin_Variation_vec[i]));
0221     }
0222
0223
0224
0225     h1D_RotatedBkg_RecoTrackletEtaPerEvt_LowerBound -> Divide(h1D_RotatedBkg_RecoTrackletEtaPerEvt_ideal);
0226     h1D_RotatedBkg_RecoTrackletEtaPerEvt_UpperBound -> Divide(h1D_RotatedBkg_RecoTrackletEtaPerEvt_ideal);
0227     TH1D * h1D_RotatedBkg_RecoTrackletEtaPerEvt_VariationMax = (TH1D*) h1D_RotatedBkg_RecoTrackletEtaPerEvt_LowerBound -> Clone("h1D_RotatedBkg_RecoTrackletEtaPerEvt_VariationMax");
0228
0229     for (int i = 1; i <= h1D_RotatedBkg_RecoTrackletEtaPerEvt_VariationMax -> GetNbinsX(); i++)
0230     {
0231         double fabs_bin_content_lower = fabs(h1D_RotatedBkg_RecoTrackletEtaPerEvt_LowerBound -> GetBinContent(i) - 1);
0232         double fabs_bin_content_upper = fabs(h1D_RotatedBkg_RecoTrackletEtaPerEvt_UpperBound -> GetBinContent(i) - 1);
0233         double bin_content = (fabs_bin_content_lower > fabs_bin_content_upper) ? fabs_bin_content_lower : fabs_bin_content_upper;
0234
0235         h1D_RotatedBkg_RecoTrackletEtaPerEvt_VariationMax -> SetBinContent(i, bin_content);
0236     }
0237
0238     TFile * file_out = new TFile(Form("%s/FromdNdEta_New.root", output_directory.c_str()), "RECREATE");
0239     h2D_BestPair_RecoTrackletEtaPerEvt -> Write();
0240     h2D_RotatedBkg_RecoTrackletEtaPerEvt -> Write();
0241
0242     c1 -> cd();
0243     h2D_BestPair_RecoTrackletEtaPerEvt -> Draw("colz0");
0244     h1D_BestPair_RecoTrackletEtaPerEvt_ideal -> Draw("hist same");
0245     c1 -> Write("c1_BestPair_RecoTrackletEtaPerEvt");
0246     c1 -> Clear();
0247
0248     c1 -> cd();
0249     h2D_RotatedBkg_RecoTrackletEtaPerEvt -> Draw("colz0");
0250     h1D_RotatedBkg_RecoTrackletEtaPerEvt_ideal -> Draw("hist same");
0251     c1 -> Write("c1_RotatedBkg_RecoTrackletEtaPerEvt");
0252     c1 -> Clear();
0253
0254     h1D_RotatedBkg_RecoTrackletEtaPerEvt_LowerBound -> Write();
0255     h1D_RotatedBkg_RecoTrackletEtaPerEvt_UpperBound -> Write();
0256     h1D_RotatedBkg_RecoTrackletEtaPerEvt_VariationMax -> Write();
0257     h1D_RotatedBkg_RecoTrackletEtaPerEvt_VariationSigma -> Write();
0258
0259     file_out -> Close();
0260
0261     // Division: ---------------------------------------------------------------------------------------------------------------
0262     c1 -> cd();
0263     h2D_BestPair_RecoTrackletEtaPerEvt -> Draw("colz0");
0264     h1D_BestPair_RecoTrackletEtaPerEvt_ideal -> Draw("hist same");
0265     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", label_text.c_str()));
0266     c1 -> Print(Form("%s/c1_BestPair_RecoTrackletEtaPerEvt.pdf", output_directory.c_str()));
0267     c1 -> Clear();
0268
0269     c1 -> cd();
0270     h2D_RotatedBkg_RecoTrackletEtaPerEvt -> Draw("colz0");
0271     h1D_RotatedBkg_RecoTrackletEtaPerEvt_ideal -> Draw("hist same");
0272     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", label_text.c_str()));
0273     c1 -> Print(Form("%s/c1_RotatedBkg_RecoTrackletEtaPerEvt.pdf", output_directory.c_str()));
0274     c1 -> Clear();
0275
0276
0277     c1 -> cd();
0278     h1D_RotatedBkg_RecoTrackletEtaPerEvt_VariationMax -> GetXaxis() -> SetTitle("Tracklet #eta");
0279     h1D_RotatedBkg_RecoTrackletEtaPerEvt_VariationMax -> GetYaxis() -> SetTitle("Variation Max in Reco. Tracklets");
0280     h1D_RotatedBkg_RecoTrackletEtaPerEvt_VariationMax -> SetLineColor(4);
0281     h1D_RotatedBkg_RecoTrackletEtaPerEvt_VariationMax -> SetLineColorAlpha(4, 0);
0282     h1D_RotatedBkg_RecoTrackletEtaPerEvt_VariationMax -> SetMinimum(0);
0283     h1D_RotatedBkg_RecoTrackletEtaPerEvt_VariationMax -> SetMaximum(0.1);
0284     h1D_RotatedBkg_RecoTrackletEtaPerEvt_VariationMax -> Draw("p");
0285
0286     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", label_text.c_str()));
0287     c1 -> Print(Form("%s/h1D_RotatedBkg_RecoTrackletEtaPerEvt_VariationMax.pdf", output_directory.c_str()));
0288     c1 -> Clear();
0289
0290     c1 -> cd();
0291     h1D_RotatedBkg_RecoTrackletEtaPerEvt_VariationSigma -> GetXaxis() -> SetTitle("Tracklet #eta");
0292     h1D_RotatedBkg_RecoTrackletEtaPerEvt_VariationSigma -> GetYaxis() -> SetTitle("3#sigma of variation in Reco. Tracklets");
0293     h1D_RotatedBkg_RecoTrackletEtaPerEvt_VariationSigma -> SetLineColorAlpha(4,0);
0294     h1D_RotatedBkg_RecoTrackletEtaPerEvt_VariationSigma -> SetMinimum(0);
0295     h1D_RotatedBkg_RecoTrackletEtaPerEvt_VariationSigma -> SetMaximum(0.1);
0296     h1D_RotatedBkg_RecoTrackletEtaPerEvt_VariationSigma -> Draw("p");
0297
0298     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", label_text.c_str()));
0299     c1 -> Print(Form("%s/h1D_RotatedBkg_RecoTrackletEtaPerEvt_VariationSigma.pdf", output_directory.c_str()));
0300     c1 -> Clear();
0301
0302     return 0;
0303 }