DST_MC/runXY_stability/Check_stability.cpp

0001 #include <iostream>
0002 #include <fstream>
0003 #include <filesystem>
0004 #include <cstdlib>
0005 #include <numeric>
0006 using namespace std;
0007
0008 #include <TFile.h>
0009 #include <TChain.h>
0010 #include <TTree.h>
0011 #include <TGraph.h>
0012 #include <TGraphErrors.h>
0013 #include <TCanvas.h>
0014 #include <TH1.h>
0015 #include <TH2.h>
0016 #include <TF1.h>
0017 #include <TLatex.h>
0018 #include <TLegend.h>
0019
0020 #include "sPhenixStyle.h"
0021
0022 double  vector_average (vector <double> input_vector) {
0023     return accumulate( input_vector.begin(), input_vector.end(), 0.0 ) / double(input_vector.size());
0024 }
0025
0026 double vector_stddev (vector <double> input_vector){
0027
0028     double sum_subt = 0;
0029     double average  = accumulate( input_vector.begin(), input_vector.end(), 0.0 ) / double(input_vector.size());
0030
0031     // cout<<"average is : "<<average<<endl;
0032
0033     for (int i=0; i<input_vector.size(); i++){ sum_subt += pow((input_vector[i] - average),2); }
0034
0035     //cout<<"sum_subt : "<<sum_subt<<endl;
0036     // cout<<"print from the function, average : "<<average<<" std : "<<stddev<<endl;
0037
0038     return sqrt( sum_subt / double(input_vector.size()-1) );
0039 }
0040
0041 vector<string> read_list(string folder_direction, string MC_list_name)
0042 {
0043     vector<string> file_list;
0044     string list_buffer;
0045     ifstream data_list;
0046     data_list.open((folder_direction + "/" + MC_list_name).c_str());
0047
0048    file_list.clear();
0049
0050     while (1)
0051     {
0052         data_list >> list_buffer;
0053         if (!data_list.good())
0054         {
0055             break;
0056         }
0057         file_list.push_back(list_buffer);
0058     }
0059     cout<<"size in the" <<MC_list_name<<": "<<file_list.size()<<endl;
0060
0061     return file_list;
0062 }
0063
0064 int main(int argc, char* argv[])
0065 {
0066     int data_type = std::stod(argv[1]); // note : 1 = MC, 0 = data
0067     int unit_tag = 1; // note : 1 = cm, 0 = mm
0068     double unit_correction = (unit_tag == 1) ? 0.1 : 1.;
0069     double frame_shift_forX = (data_type == 1) ? 0. : 0.3;
0070     double frame_shift_forY = (data_type == 1) ? 0. : 0.6;
0071
0072     // todo : modify the followings if the vertex is expected to be in different region
0073     // note : still using mm, here
0074     double X_range_l = (-1. + frame_shift_forX) * unit_correction;
0075     double X_range_r = (0.5 + frame_shift_forX) * unit_correction;
0076     double Y_range_l = (1.7 + frame_shift_forY) * unit_correction;
0077     double Y_range_r = (3.2 + frame_shift_forY) * unit_correction;
0078
0079     string input_folder = argv[2];
0080     // string input_folder = "/sphenix/user/ChengWei/sPH_dNdeta/Sim_Ntuple_HIJING_new_20240424_HR/avg_vtxXY";
0081     // string input_folder = "/sphenix/user/ChengWei/INTT/INTT_commissioning/ZeroField/F4A_20869/2024_05_07/folder_Data_CombinedNtuple_Run20869_HotDead_BCO_ADC_Survey/avg_vtxXY";
0082     string file_list_name = "file_list.txt";
0083     string output_directory = input_folder + "/" + "merged_result";
0084
0085
0086
0087     string label_text = (data_type == 1)? "Simulation" : "Work-in-progress";
0088     string unit_text = (unit_tag == 1) ? "cm" : "mm";
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 %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     if(std::filesystem::exists(Form("%s/%s",input_folder.c_str(),file_list_name.c_str())) == false){
0104         system(Form("ls %s/runXY_*/run_XY_tree.root > %s/%s", input_folder.c_str(), input_folder.c_str(), file_list_name.c_str()));
0105         cout<<"----------- file list generated -----------"<<endl;
0106         system(Form("cat %s/%s", input_folder.c_str(), file_list_name.c_str()));
0107         cout<<"----------- file list generated -----------"<<endl;
0108     }
0109     else
0110     {
0111         cout<<"----------- found the file list -----------"<<endl;
0112         system(Form("cat %s/%s", input_folder.c_str(), file_list_name.c_str()));
0113         cout<<"----------- found the file list -----------"<<endl;
0114     }
0115
0116     TChain *chain_in = new TChain("tree");
0117     vector<string> file_list = read_list(input_folder, file_list_name);
0118     for (string file_name : file_list) { chain_in -> Add((file_name).c_str()); }
0119     cout<<"chain_in -> GetEntries() : "<<chain_in -> GetEntries()<<endl;
0120
0121     long in_start_evt;
0122     long in_end_evt;
0123     double in_quadrant_corner_X;
0124     double in_quadrant_corner_Y;
0125     double in_quadrant_center_X;
0126     double in_quadrant_center_Y;
0127     double in_line_filled_mean_X;
0128     double in_line_filled_mean_Y;
0129     double in_line_filled_meanE_X;
0130     double in_line_filled_meanE_Y;
0131     double in_line_filled_stddev_X;
0132     double in_line_filled_stddev_Y;
0133     double in_line_filled_variance_X;
0134     double in_line_filled_variance_Y;
0135     double in_line_filled_covariance;
0136
0137     chain_in -> SetBranchAddress("start_evt",&in_start_evt);
0138     chain_in -> SetBranchAddress("end_evt",&in_end_evt);
0139     chain_in -> SetBranchAddress("quadrant_corner_X",&in_quadrant_corner_X);
0140     chain_in -> SetBranchAddress("quadrant_corner_Y",&in_quadrant_corner_Y);
0141     chain_in -> SetBranchAddress("quadrant_center_X",&in_quadrant_center_X);
0142     chain_in -> SetBranchAddress("quadrant_center_Y",&in_quadrant_center_Y);
0143     chain_in -> SetBranchAddress("line_filled_mean_X",&in_line_filled_mean_X);
0144     chain_in -> SetBranchAddress("line_filled_mean_Y",&in_line_filled_mean_Y);
0145     chain_in -> SetBranchAddress("line_filled_meanE_X",&in_line_filled_meanE_X);
0146     chain_in -> SetBranchAddress("line_filled_meanE_Y",&in_line_filled_meanE_Y);
0147     chain_in -> SetBranchAddress("line_filled_stddev_X",&in_line_filled_stddev_X);
0148     chain_in -> SetBranchAddress("line_filled_stddev_Y",&in_line_filled_stddev_Y);
0149     chain_in -> SetBranchAddress("line_filled_variance_X",&in_line_filled_variance_X);
0150     chain_in -> SetBranchAddress("line_filled_variance_Y",&in_line_filled_variance_Y);
0151     chain_in -> SetBranchAddress("line_filled_covariance",&in_line_filled_covariance);
0152
0153     TH1F * quadrant_VtxX_1D = new TH1F("quadrant_VtxX",Form("quadrant_VtxX;X axis [%s];Entry",unit_text.c_str()),100,X_range_l,X_range_r);
0154     TH1F * quadrant_VtxY_1D = new TH1F("quadrant_VtxY",Form("quadrant_VtxY;Y axis [%s];Entry",unit_text.c_str()),100,Y_range_l,Y_range_r);
0155     TH1F * line_filled_mean_X_1D = new TH1F("line_filled_mean_X",Form("line_filled_mean_X;X axis [%s];Entry",unit_text.c_str()),100,X_range_l,X_range_r);
0156     TH1F * line_filled_mean_Y_1D = new TH1F("line_filled_mean_Y",Form("line_filled_mean_Y;Y axis [%s];Entry",unit_text.c_str()),100,Y_range_l,Y_range_r);
0157
0158     TH1F * line_filled_stddev_X_1D = new TH1F("line_filled_stddev_X_1D",Form("line_filled_stddev_X_1D;StdDev, X axis [%s];Entry",unit_text.c_str()),100,0,1 * unit_correction);
0159     TH1F * line_filled_stddev_Y_1D = new TH1F("line_filled_stddev_Y_1D",Form("line_filled_stddev_Y_1D;StdDev, Y axis [%s];Entry",unit_text.c_str()),100,0,1 * unit_correction);
0160
0161     TH1F * line_filled_variance_X_1D = new TH1F("line_filled_variance_X_1D",Form("line_filled_variance_X_1D;Variance, X axis [%s];Entry",unit_text.c_str()),100,0,1 * unit_correction);
0162     TH1F * line_filled_variance_Y_1D = new TH1F("line_filled_variance_Y_1D",Form("line_filled_variance_Y_1D;Variance, Y axis [%s];Entry",unit_text.c_str()),100,0,1 * unit_correction);
0163
0164     TH1F * line_filled_covariance_1D = new TH1F("line_filled_covariance_1D",Form("line_filled_covariance_1D;Covariance [%s];Entry",unit_text.c_str()),100,-0.06 * unit_correction,0.06 * unit_correction);
0165
0166
0167     vector<double> evt_index_vec; evt_index_vec.clear();
0168     vector<double> evt_index_range_vec; evt_index_range_vec.clear();
0169
0170     vector<double> quadrant_VtxX_vec; quadrant_VtxX_vec.clear();
0171     vector<double> quadrant_VtxY_vec; quadrant_VtxY_vec.clear();
0172     vector<double> quadrant_VtxXerror_vec; quadrant_VtxXerror_vec.clear();
0173     vector<double> quadrant_VtxYerror_vec; quadrant_VtxYerror_vec.clear();
0174
0175     vector<double> line_filled_mean_X_vec; line_filled_mean_X_vec.clear();
0176     vector<double> line_filled_mean_Y_vec; line_filled_mean_Y_vec.clear();
0177     vector<double> line_filled_mean_Xerror_vec; line_filled_mean_Xerror_vec.clear();
0178     vector<double> line_filled_mean_Yerror_vec; line_filled_mean_Yerror_vec.clear();
0179
0180     vector<double> combined_index_vec; combined_index_vec.clear();
0181     vector<double> combined_index_range_vec; combined_index_range_vec.clear();
0182     vector<double> combined_X_vec; combined_X_vec.clear();
0183     vector<double> combined_Y_vec; combined_Y_vec.clear();
0184     vector<double> combined_XE_vec; combined_XE_vec.clear();
0185     vector<double> combined_YE_vec; combined_YE_vec.clear();
0186
0187
0188     vector<double> zero_vec; zero_vec.clear();
0189     vector<double> starteID_vec; starteID_vec.clear();
0190     vector<double> endeID_vec; endeID_vec.clear();
0191
0192     SetsPhenixStyle();
0193     TCanvas * c1 = new TCanvas("c1","c1",950,800);
0194     c1 -> cd();
0195
0196     TLatex * ltx = new TLatex();
0197     ltx->SetNDC();
0198     ltx->SetTextSize(0.045);
0199     ltx->SetTextAlign(31);
0200
0201     TLatex * draw_text = new TLatex();
0202     draw_text -> SetNDC();
0203     draw_text -> SetTextSize(0.03);
0204
0205     for (int i = 0; i < chain_in -> GetEntries(); i++)
0206     {
0207         chain_in -> GetEntry(i);
0208
0209         // note : if changing from mm to cm is needed
0210         in_quadrant_corner_X = in_quadrant_corner_X * unit_correction;
0211         in_quadrant_corner_Y = in_quadrant_corner_Y * unit_correction;
0212         in_quadrant_center_X = in_quadrant_center_X * unit_correction;
0213         in_quadrant_center_Y = in_quadrant_center_Y * unit_correction;
0214         in_line_filled_mean_X = in_line_filled_mean_X * unit_correction;
0215         in_line_filled_mean_Y = in_line_filled_mean_Y * unit_correction;
0216         in_line_filled_meanE_X = in_line_filled_meanE_X * unit_correction;
0217         in_line_filled_meanE_Y = in_line_filled_meanE_Y * unit_correction;
0218         in_line_filled_stddev_X = in_line_filled_stddev_X * unit_correction;
0219         in_line_filled_stddev_Y = in_line_filled_stddev_Y * unit_correction;
0220         in_line_filled_variance_X = in_line_filled_variance_X * unit_correction;
0221         in_line_filled_variance_Y = in_line_filled_variance_Y * unit_correction;
0222         in_line_filled_covariance = in_line_filled_covariance * unit_correction;
0223
0224         double quadrant_VtxX = (in_quadrant_corner_X + in_quadrant_center_X)/2.;
0225         double quadrant_VtxY = (in_quadrant_corner_Y + in_quadrant_center_Y)/2.;
0226         double quadrant_WindowX = (fabs(in_quadrant_corner_X - in_quadrant_center_X))/2.;
0227         double quadrant_WindowY = (fabs(in_quadrant_corner_Y - in_quadrant_center_Y))/2.;
0228
0229         double line_filled_mean_X = in_line_filled_mean_X;
0230         double line_filled_mean_Y = in_line_filled_mean_Y;
0231         double line_filled_meanE_X = in_line_filled_meanE_X;
0232         double line_filled_meanE_Y = in_line_filled_meanE_Y;
0233         cout<<Form("unit : %s",unit_text.c_str())<<", index : "<<i<<" quadrant_VtxX : "<<quadrant_VtxX<<" quadrant_VtxY : "<<quadrant_VtxY<<" line_filled_mean_X : "<<line_filled_mean_X<<" line_filled_mean_Y : "<<line_filled_mean_Y<<endl;
0234         cout<<Form("unit : %s",unit_text.c_str())<<", index : "<<i<<" quadrant_VtxXE : "<<quadrant_WindowX<<" quadrant_VtxYE : "<<quadrant_WindowY<<" line_filled_meanE_X : "<<line_filled_meanE_X<<" line_filled_meaE_Y : "<<line_filled_meanE_Y<<endl;
0235
0236         starteID_vec.push_back(in_start_evt);
0237         endeID_vec.push_back(in_end_evt);
0238
0239         evt_index_vec.push_back((in_start_evt + in_end_evt)/2.);
0240         evt_index_range_vec.push_back((in_end_evt - in_start_evt)/2.);
0241         quadrant_VtxX_vec.push_back(quadrant_VtxX);
0242         quadrant_VtxY_vec.push_back(quadrant_VtxY);
0243         line_filled_mean_X_vec.push_back(line_filled_mean_X);
0244         line_filled_mean_Y_vec.push_back(line_filled_mean_Y);
0245         zero_vec.push_back(0);
0246
0247         quadrant_VtxXerror_vec.push_back(quadrant_WindowX);
0248         quadrant_VtxYerror_vec.push_back(quadrant_WindowY);
0249         line_filled_mean_Xerror_vec.push_back(line_filled_meanE_X);
0250         line_filled_mean_Yerror_vec.push_back(line_filled_meanE_Y);
0251
0252         combined_index_vec.push_back((in_start_evt + in_end_evt)/2.);
0253         combined_index_range_vec.push_back((in_end_evt - in_start_evt)/2.);
0254         combined_X_vec.push_back(quadrant_VtxX);
0255         combined_XE_vec.push_back(quadrant_WindowX);
0256         combined_Y_vec.push_back(quadrant_VtxY);
0257         combined_YE_vec.push_back(quadrant_WindowY);
0258
0259         combined_index_vec.push_back((in_start_evt + in_end_evt)/2.);
0260         combined_index_range_vec.push_back((in_end_evt - in_start_evt)/2.);
0261         combined_X_vec.push_back(line_filled_mean_X);
0262         combined_XE_vec.push_back(line_filled_meanE_X);
0263         combined_Y_vec.push_back(line_filled_mean_Y);
0264         combined_YE_vec.push_back(line_filled_meanE_Y);
0265
0266
0267
0268
0269         quadrant_VtxX_1D -> Fill(quadrant_VtxX);
0270         quadrant_VtxY_1D -> Fill(quadrant_VtxY);
0271         line_filled_mean_X_1D -> Fill(line_filled_mean_X);
0272         line_filled_mean_Y_1D -> Fill(line_filled_mean_Y);
0273
0274         line_filled_stddev_X_1D -> Fill(in_line_filled_stddev_X);
0275         line_filled_stddev_Y_1D -> Fill(in_line_filled_stddev_Y);
0276         line_filled_variance_X_1D -> Fill(in_line_filled_variance_X);
0277         line_filled_variance_Y_1D -> Fill(in_line_filled_variance_Y);
0278         line_filled_covariance_1D -> Fill(in_line_filled_covariance);
0279     }
0280
0281     for (int i = 0; i < starteID_vec.size(); i++)
0282     {
0283         cout<<"file : "<<i<<", start: "<<starteID_vec[i]<<", end: "<<endeID_vec[i]<<endl;
0284     }
0285
0286     TGraphErrors * g_quadrant_VtxX_index = new TGraphErrors(chain_in -> GetEntries(), &evt_index_vec[0], &quadrant_VtxX_vec[0], &evt_index_range_vec[0], &quadrant_VtxXerror_vec[0]);
0287     g_quadrant_VtxX_index -> GetXaxis() -> SetNdivisions(505);
0288     g_quadrant_VtxX_index -> SetMarkerStyle(20);
0289     g_quadrant_VtxX_index -> SetMarkerSize(0.7);
0290     g_quadrant_VtxX_index -> SetMarkerColor(1);
0291     g_quadrant_VtxX_index -> GetXaxis() -> SetTitle("Event ID");
0292     g_quadrant_VtxX_index -> GetYaxis() -> SetTitle(Form("Average beam position in X axis [%s]",unit_text.c_str()));
0293     g_quadrant_VtxX_index -> GetXaxis() -> SetLimits(-1000, g_quadrant_VtxX_index -> GetPointX(g_quadrant_VtxX_index -> GetN() - 1) + 10000);
0294     g_quadrant_VtxX_index -> GetYaxis() -> SetRangeUser(quadrant_VtxX_1D -> GetXaxis() -> GetXmin(),quadrant_VtxX_1D -> GetXaxis() -> GetXmax());
0295
0296     TGraphErrors * g_quadrant_VtxY_index = new TGraphErrors(chain_in -> GetEntries(), &evt_index_vec[0], &quadrant_VtxY_vec[0], &evt_index_range_vec[0], &quadrant_VtxYerror_vec[0]);
0297     g_quadrant_VtxY_index -> GetXaxis() -> SetNdivisions(505);
0298     g_quadrant_VtxY_index -> SetMarkerStyle(20);
0299     g_quadrant_VtxY_index -> SetMarkerSize(0.7);
0300     g_quadrant_VtxY_index -> SetMarkerColor(1);
0301     g_quadrant_VtxY_index -> GetXaxis() -> SetTitle("Event ID");
0302     g_quadrant_VtxY_index -> GetYaxis() -> SetTitle(Form("Average beam position in Y axis [%s]",unit_text.c_str()));
0303     g_quadrant_VtxY_index -> GetXaxis() -> SetLimits(-1000, g_quadrant_VtxY_index -> GetPointX(g_quadrant_VtxY_index -> GetN() - 1) + 10000);
0304     g_quadrant_VtxY_index -> GetYaxis() -> SetRangeUser(quadrant_VtxY_1D -> GetXaxis() -> GetXmin(),quadrant_VtxY_1D -> GetXaxis() -> GetXmax());
0305
0306     TGraphErrors * g_line_filled_mean_X_index = new TGraphErrors(chain_in -> GetEntries(), &evt_index_vec[0], &line_filled_mean_X_vec[0], &evt_index_range_vec[0], &line_filled_mean_Xerror_vec[0]);
0307     g_line_filled_mean_X_index -> GetXaxis() -> SetNdivisions(505);
0308     g_line_filled_mean_X_index -> SetMarkerStyle(20);
0309     g_line_filled_mean_X_index -> SetMarkerSize(0.7);
0310     g_line_filled_mean_X_index -> SetMarkerColor(4);
0311     g_line_filled_mean_X_index -> SetLineColor(4);
0312     g_line_filled_mean_X_index -> GetXaxis() -> SetTitle("Event ID");
0313     g_line_filled_mean_X_index -> GetYaxis() -> SetTitle(Form("Beam Position in X axis [%s]",unit_text.c_str()));
0314
0315     TGraphErrors * g_line_filled_mean_Y_index = new TGraphErrors(chain_in -> GetEntries(), &evt_index_vec[0], &line_filled_mean_Y_vec[0], &evt_index_range_vec[0], &line_filled_mean_Yerror_vec[0]);
0316     g_line_filled_mean_Y_index -> GetXaxis() -> SetNdivisions(505);
0317     g_line_filled_mean_Y_index -> SetMarkerStyle(20);
0318     g_line_filled_mean_Y_index -> SetMarkerSize(0.7);
0319     g_line_filled_mean_Y_index -> SetMarkerColor(4);
0320     g_line_filled_mean_Y_index -> SetLineColor(4);
0321     g_line_filled_mean_Y_index -> GetXaxis() -> SetTitle("Event ID");
0322     g_line_filled_mean_Y_index -> GetYaxis() -> SetTitle(Form("Beam Position in Y axis [%s]",unit_text.c_str()));
0323
0324
0325     TGraphErrors * combined_X_index_grE = new TGraphErrors(combined_index_vec.size(), &combined_index_vec[0], &combined_X_vec[0], &combined_index_range_vec[0], &combined_XE_vec[0]);
0326     TF1 * pol0_fit_X = new TF1("pol0_fit_X", "pol0", 0, g_line_filled_mean_Y_index->GetPointX(g_line_filled_mean_Y_index -> GetN() - 1) * 1.5);
0327     pol0_fit_X -> SetLineColor(2);
0328     pol0_fit_X -> SetLineStyle(9);
0329     pol0_fit_X -> SetLineWidth(1);
0330     combined_X_index_grE -> Fit(pol0_fit_X,"NQ");
0331
0332     TGraphErrors * combined_Y_index_grE = new TGraphErrors(combined_index_vec.size(), &combined_index_vec[0], &combined_Y_vec[0], &combined_index_range_vec[0], &combined_YE_vec[0]);
0333     TF1 * pol0_fit_Y = new TF1("pol0_fit_Y", "pol0", 0, g_line_filled_mean_Y_index->GetPointX(g_line_filled_mean_Y_index -> GetN() - 1) * 1.5);
0334     pol0_fit_Y -> SetLineColor(2);
0335     pol0_fit_Y -> SetLineStyle(9);
0336     pol0_fit_Y -> SetLineWidth(1);
0337     combined_Y_index_grE -> Fit(pol0_fit_Y,"NQ");
0338
0339     TLegend * leg = new TLegend(0.21,0.75,0.31,0.9);
0340     leg -> AddEntry(g_quadrant_VtxX_index, "Quadrant method", "lep");
0341     leg -> AddEntry(g_line_filled_mean_X_index, "Line-filled method", "lep");
0342
0343
0344     c1 -> cd();
0345     g_quadrant_VtxX_index -> Draw("AP");
0346     g_line_filled_mean_X_index -> Draw("p same");
0347     draw_text -> DrawLatex(0.21,0.71,Form("Numeric average vertex X: %.4f %s", ( vector_average(line_filled_mean_X_vec) + vector_average(quadrant_VtxX_vec) )/2., unit_text.c_str() ));
0348     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", label_text.c_str()));
0349     leg -> Draw("same");
0350     c1 -> Print((output_directory + "/VtxX_index.pdf").c_str());
0351     c1 -> Clear();
0352
0353     g_quadrant_VtxX_index -> Draw("AP");
0354     g_line_filled_mean_X_index -> Draw("p same");
0355     draw_text -> DrawLatex(0.21,0.71,Form("Fit vertex X: %.4f #pm%.4f %s", pol0_fit_X->GetParameter(0), pol0_fit_X->GetParError(0), unit_text.c_str() ));
0356     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", label_text.c_str()));
0357     leg -> Draw("same");
0358     pol0_fit_X -> Draw("lsame");
0359     c1 -> Print((output_directory + "/VtxX_index_fit.pdf").c_str());
0360     c1 -> Clear();
0361
0362     c1 -> cd();
0363     g_quadrant_VtxY_index -> Draw("AP");
0364     g_line_filled_mean_Y_index -> Draw("p same");
0365     draw_text -> DrawLatex(0.21,0.71,Form("Numeric average vertex Y: %.4f %s", ( vector_average(line_filled_mean_Y_vec) + vector_average(quadrant_VtxY_vec) )/2., unit_text.c_str() ));
0366     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", label_text.c_str()));
0367     leg -> Draw("same");
0368     c1 -> Print((output_directory + "/VtxY_index.pdf").c_str());
0369     c1 -> Clear();
0370
0371     g_quadrant_VtxY_index -> Draw("AP");
0372     g_line_filled_mean_Y_index -> Draw("p same");
0373     draw_text -> DrawLatex(0.21,0.71,Form("Fit vertex Y: %.4f #pm%.4f %s", pol0_fit_Y->GetParameter(0), pol0_fit_Y->GetParError(0), unit_text.c_str() ));
0374     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", label_text.c_str()));
0375     leg -> Draw("same");
0376     pol0_fit_Y -> Draw("lsame");
0377     c1 -> Print((output_directory + "/VtxY_index_fit.pdf").c_str());
0378     c1 -> Clear();
0379
0380     c1 -> cd();
0381     quadrant_VtxX_1D -> Draw("hist");
0382     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", label_text.c_str()));
0383     draw_text -> DrawLatex(0.2, 0.84, Form("Entry: %.0f", quadrant_VtxX_1D->GetEntries()));
0384     draw_text -> DrawLatex(0.2, 0.80, Form("Mean: %.4f", quadrant_VtxX_1D->GetMean()));
0385     draw_text -> DrawLatex(0.2, 0.76, Form("StdDev: %.4f", quadrant_VtxX_1D->GetStdDev()));
0386     c1 -> Print((output_directory + "/quadrant_VtxX_1D.pdf").c_str());
0387     c1 -> Clear();
0388
0389     c1 -> cd();
0390     quadrant_VtxY_1D -> Draw("hist");
0391     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", label_text.c_str()));
0392     draw_text -> DrawLatex(0.2, 0.84, Form("Entry: %.0f", quadrant_VtxY_1D->GetEntries()));
0393     draw_text -> DrawLatex(0.2, 0.80, Form("Mean: %.4f", quadrant_VtxY_1D->GetMean()));
0394     draw_text -> DrawLatex(0.2, 0.76, Form("StdDev: %.4f", quadrant_VtxY_1D->GetStdDev()));
0395     c1 -> Print((output_directory + "/quadrant_VtxY_1D.pdf").c_str());
0396     c1 -> Clear();
0397
0398     c1 -> cd();
0399     line_filled_mean_X_1D -> Draw("hist");
0400     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", label_text.c_str()));
0401     draw_text -> DrawLatex(0.2, 0.84, Form("Entry: %.0f", line_filled_mean_X_1D->GetEntries()));
0402     draw_text -> DrawLatex(0.2, 0.80, Form("Mean: %.4f", line_filled_mean_X_1D->GetMean()));
0403     draw_text -> DrawLatex(0.2, 0.76, Form("StdDev: %.4f", line_filled_mean_X_1D->GetStdDev()));
0404     c1 -> Print((output_directory + "/line_filled_mean_X_1D.pdf").c_str());
0405     c1 -> Clear();
0406
0407     c1 -> cd();
0408     line_filled_mean_Y_1D -> Draw("hist");
0409     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", label_text.c_str()));
0410     draw_text -> DrawLatex(0.2, 0.84, Form("Entry: %.0f", line_filled_mean_Y_1D->GetEntries()));
0411     draw_text -> DrawLatex(0.2, 0.80, Form("Mean: %.4f", line_filled_mean_Y_1D->GetMean()));
0412     draw_text -> DrawLatex(0.2, 0.76, Form("StdDev: %.4f", line_filled_mean_Y_1D->GetStdDev()));
0413     c1 -> Print((output_directory + "/line_filled_mean_Y_1D.pdf").c_str());
0414     c1 -> Clear();
0415
0416
0417     c1 -> cd();
0418     line_filled_stddev_X_1D -> Draw("hist");
0419     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", label_text.c_str()));
0420     draw_text -> DrawLatex(0.2, 0.84, Form("Entry: %.0f", line_filled_stddev_X_1D->GetEntries()));
0421     draw_text -> DrawLatex(0.2, 0.80, Form("Mean: %.4f", line_filled_stddev_X_1D->GetMean()));
0422     draw_text -> DrawLatex(0.2, 0.76, Form("StdDev: %.4f", line_filled_stddev_X_1D->GetStdDev()));
0423     c1 -> Print((output_directory + "/line_filled_stddev_X_1D.pdf").c_str());
0424     c1 -> Clear();
0425
0426     c1 -> cd();
0427     line_filled_stddev_Y_1D -> Draw("hist");
0428     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", label_text.c_str()));
0429     draw_text -> DrawLatex(0.2, 0.84, Form("Entry: %.0f", line_filled_stddev_Y_1D->GetEntries()));
0430     draw_text -> DrawLatex(0.2, 0.80, Form("Mean: %.4f", line_filled_stddev_Y_1D->GetMean()));
0431     draw_text -> DrawLatex(0.2, 0.76, Form("StdDev: %.4f", line_filled_stddev_Y_1D->GetStdDev()));
0432     c1 -> Print((output_directory + "/line_filled_stddev_Y_1D.pdf").c_str());
0433     c1 -> Clear();
0434
0435     c1 -> cd();
0436     line_filled_variance_X_1D -> Draw("hist");
0437     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", label_text.c_str()));
0438     draw_text -> DrawLatex(0.55, 0.84, Form("Entry: %.0f", line_filled_variance_X_1D->GetEntries()));
0439     draw_text -> DrawLatex(0.55, 0.80, Form("Mean: %.4f", line_filled_variance_X_1D->GetMean()));
0440     draw_text -> DrawLatex(0.55, 0.76, Form("StdDev: %.4f", line_filled_variance_X_1D->GetStdDev()));
0441     c1 -> Print((output_directory + "/line_filled_variance_X_1D.pdf").c_str());
0442     c1 -> Clear();
0443
0444     c1 -> cd();
0445     line_filled_variance_Y_1D -> Draw("hist");
0446     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", label_text.c_str()));
0447     draw_text -> DrawLatex(0.55, 0.84, Form("Entry: %.0f", line_filled_variance_Y_1D->GetEntries()));
0448     draw_text -> DrawLatex(0.55, 0.80, Form("Mean: %.4f", line_filled_variance_Y_1D->GetMean()));
0449     draw_text -> DrawLatex(0.55, 0.76, Form("StdDev: %.4f", line_filled_variance_Y_1D->GetStdDev()));
0450     c1 -> Print((output_directory + "/line_filled_variance_Y_1D.pdf").c_str());
0451     c1 -> Clear();
0452
0453     c1 -> cd();
0454     line_filled_covariance_1D -> Draw("hist");
0455     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", label_text.c_str()));
0456     draw_text -> DrawLatex(0.2, 0.84, Form("Entry: %.0f", line_filled_covariance_1D->GetEntries()));
0457     draw_text -> DrawLatex(0.2, 0.80, Form("Mean: %.4f", line_filled_covariance_1D->GetMean()));
0458     draw_text -> DrawLatex(0.2, 0.76, Form("StdDev: %.4f", line_filled_covariance_1D->GetStdDev()));
0459     c1 -> Print((output_directory + "/line_filled_covariance_1D.pdf").c_str());
0460     c1 -> Clear();
0461
0462     double out_quadrant_avgX;
0463     double out_quadrant_avgXE;
0464     double out_quadrant_avgY;
0465     double out_quadrant_avgYE;
0466     double out_linefilled_avgX;
0467     double out_linefilled_avgXE;
0468     double out_linefilled_avgY;
0469     double out_linefilled_avgYE;
0470     double out_final_num_avgX;
0471     double out_final_num_avgY;
0472     double out_final_fit_avgX;
0473     double out_final_fit_avgY;
0474
0475     TFile * file_out = new TFile (Form("%s/determined_avg_vtxXY.root",output_directory.c_str()),"recreate");
0476     TTree * tree_out = new TTree("tree","tree");
0477     tree_out -> Branch("quadrant_avgX",   &out_quadrant_avgX);
0478     tree_out -> Branch("quadrant_avgXE",    &out_quadrant_avgXE);
0479     tree_out -> Branch("quadrant_avgY",   &out_quadrant_avgY);
0480     tree_out -> Branch("quadrant_avgYE",    &out_quadrant_avgYE);
0481     tree_out -> Branch("linefilled_avgX",  &out_linefilled_avgX);
0482     tree_out -> Branch("linefilled_avgXE", &out_linefilled_avgXE);
0483     tree_out -> Branch("linefilled_avgY",  &out_linefilled_avgY);
0484     tree_out -> Branch("linefilled_avgYE", &out_linefilled_avgYE);
0485     tree_out -> Branch("final_num_avgX",   &out_final_num_avgX);
0486     tree_out -> Branch("final_num_avgY",   &out_final_num_avgY);
0487     tree_out -> Branch("final_fit_avgX",   &out_final_fit_avgX);
0488     tree_out -> Branch("final_fit_avgY",   &out_final_fit_avgY);
0489
0490     out_quadrant_avgX = vector_average(quadrant_VtxX_vec);
0491     out_quadrant_avgXE = vector_stddev(quadrant_VtxX_vec);
0492     out_quadrant_avgY = vector_average(quadrant_VtxY_vec);
0493     out_quadrant_avgYE = vector_stddev(quadrant_VtxY_vec);
0494
0495     out_linefilled_avgX = vector_average(line_filled_mean_X_vec);
0496     out_linefilled_avgXE = vector_stddev(line_filled_mean_X_vec);
0497     out_linefilled_avgY = vector_average(line_filled_mean_Y_vec);
0498     out_linefilled_avgYE = vector_stddev(line_filled_mean_Y_vec);
0499
0500     out_final_num_avgX = ( vector_average(line_filled_mean_X_vec) + vector_average(quadrant_VtxX_vec) )/2.;
0501     out_final_num_avgY = ( vector_average(line_filled_mean_Y_vec) + vector_average(quadrant_VtxY_vec) )/2.;
0502     out_final_fit_avgX = pol0_fit_X->GetParameter(0);
0503     out_final_fit_avgY = pol0_fit_Y->GetParameter(0);
0504
0505     tree_out -> Fill();
0506     file_out -> cd();
0507     // tree_out->SetDirectory(file_out);
0508     tree_out -> Write("", TObject::kOverwrite);
0509     file_out -> Close();
0510
0511
0512     cout<<"unit : "<<unit_text<<endl;
0513     cout<<"line filled covariance info : "<<line_filled_covariance_1D->GetMean()<<" "<<line_filled_covariance_1D->GetStdDev()<<endl;
0514     cout<<"line filled vairance X info : "<<line_filled_variance_X_1D->GetMean()<<" "<<line_filled_variance_X_1D->GetStdDev()<<endl;
0515     cout<<"line filled vairance Y info : "<<line_filled_variance_Y_1D->GetMean()<<" "<<line_filled_variance_Y_1D->GetStdDev()<<endl;
0516     cout<<"line filled stdDev X   info : " <<line_filled_stddev_X_1D->GetMean()<<" "<<line_filled_stddev_X_1D->GetStdDev()<<endl;
0517     cout<<"line filled stdDev Y   info : " <<line_filled_stddev_Y_1D->GetMean()<<" "<<line_filled_stddev_Y_1D->GetStdDev()<<endl;
0518     cout<<endl;
0519     cout<<"unit : "<<unit_text<<endl;
0520     cout<<"final average vertex XY should be used : "<<endl;
0521     cout<<"line filled X : "<<vector_average(line_filled_mean_X_vec)<<" +/- "<<vector_stddev(line_filled_mean_X_vec)<<endl;
0522     cout<<"line filled Y : "<<vector_average(line_filled_mean_Y_vec)<<" +/- "<<vector_stddev(line_filled_mean_Y_vec)<<endl;
0523     cout<<"quadrant X : "<<vector_average(quadrant_VtxX_vec)<<" +/- "<<vector_stddev(quadrant_VtxX_vec)<<endl;
0524     cout<<"quadrant Y : "<<vector_average(quadrant_VtxY_vec)<<" +/- "<<vector_stddev(quadrant_VtxY_vec)<<endl;
0525     cout<<"avg: {"<< out_final_num_avgX <<Form(" * %s, ",unit_text.c_str())<< out_final_num_avgY <<Form(" * %s}",unit_text.c_str())<<endl;
0526     cout<<"Fit avg: {"<< out_final_fit_avgX <<Form(" * %s, ",unit_text.c_str())<< out_final_fit_avgY <<Form(" * %s}",unit_text.c_str())<<endl;
0527     return 0;
0528 }