Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2025-08-06 08:12:33

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.C"
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 MakeSummary(
0065 
0066     // note: MC
0067     int data_type = 1, // note : 1 = MC
0068     std::string input_folder = "/sphenix/user/ChengWei/sPH_dNdeta/Run24AuAuMC/Sim_AMPT_MDC2_ana472_20250310/Run7/AvgVtxXY/completed",
0069     std::string filename_NoIndex = "MC_AvgVtxXY",
0070 
0071     // // note : data
0072     // int data_type = 0, // note : 0 = data    
0073     // std::string input_folder = "/sphenix/tg/tg01/commissioning/INTT/work/cwshih/seflgendata/run_54280_HR_Feb102025/Run5/AvgVtxXY_nominal/completed",
0074     // std::string filename_NoIndex = "Data_AvgVtxXY_00054280",
0075 
0076     double frame_shift_forX = 0.03, // note : cm
0077     double frame_shift_forY = -0.01 // note : cm
0078 )
0079 {
0080     double X_range_l = (-0.1 + frame_shift_forX);
0081     double X_range_r = (0.05 + frame_shift_forX);
0082     double Y_range_l = (0.17 + frame_shift_forY);
0083     double Y_range_r = (0.32 + frame_shift_forY);
0084 
0085     filename_NoIndex = (filename_NoIndex.back() == '_') ? filename_NoIndex.substr(0, filename_NoIndex.size() - 1) : filename_NoIndex;
0086 
0087 
0088     string file_list_name = "file_list.txt";
0089     string output_directory = input_folder + "/" + "merged_result";
0090     string label_text = (data_type == 1)? "Simulation" : "Internal";
0091 
0092     // note : to generate the merged_result folder
0093     if(std::filesystem::exists(Form("%s",output_directory.c_str())) == false){
0094         system(Form("mkdir -p %s", output_directory.c_str()));
0095         cout<<"----------- check folder exists -----------"<<endl;
0096         system(Form("ls %s", output_directory.c_str()));
0097     }
0098     else 
0099     {
0100         cout<<"----------- folder exists already -----------"<<endl;
0101         system(Form("ls %s", output_directory.c_str()));
0102     }
0103 
0104     // note : to generate the file_list.txt
0105     // todo: the file_name
0106     if(true/*std::filesystem::exists(Form("%s/%s",input_folder.c_str(),file_list_name.c_str())) == false*/){
0107         system(Form("ls %s/%s_*.root > %s/%s", input_folder.c_str(), filename_NoIndex.c_str(), input_folder.c_str(), file_list_name.c_str()));
0108         cout<<"----------- file list generated -----------"<<endl;
0109         system(Form("cat %s/%s", input_folder.c_str(), file_list_name.c_str()));
0110         cout<<"----------- file list generated -----------"<<endl;
0111     }
0112     else 
0113     {
0114         cout<<"----------- found the file list -----------"<<endl;
0115         system(Form("cat %s/%s", input_folder.c_str(), file_list_name.c_str()));
0116         cout<<"----------- found the file list -----------"<<endl;
0117     }
0118     
0119     TChain *chain_in = new TChain("tree_vtxXY");
0120     vector<string> file_list = read_list(input_folder, file_list_name);
0121     for (string file_name : file_list) { chain_in -> Add((file_name).c_str()); }
0122     cout<<"chain_in -> GetEntries() : "<<chain_in -> GetEntries()<<endl;
0123 
0124 
0125     double in_vtxX_quadrant;
0126     double in_vtxXE_quadrant;
0127     double in_vtxY_quadrant;
0128     double in_vtxYE_quadrant;
0129     double in_vtxX_LineFill;
0130     double in_vtxXE_LineFill;
0131     double in_vtxXStdDev_LineFill;
0132     double in_vtxY_LineFill;
0133     double in_vtxYE_LineFill;
0134     double in_vtxYStdDev_LineFill;
0135     int in_run_nEvents;
0136     int in_job_index;
0137     int in_file_total_event;
0138 
0139 
0140     chain_in -> SetBranchAddress("vtxX_quadrant", &in_vtxX_quadrant);
0141     chain_in -> SetBranchAddress("vtxXE_quadrant", &in_vtxXE_quadrant);
0142     chain_in -> SetBranchAddress("vtxY_quadrant", &in_vtxY_quadrant);
0143     chain_in -> SetBranchAddress("vtxYE_quadrant", &in_vtxYE_quadrant);
0144     chain_in -> SetBranchAddress("vtxX_LineFill", &in_vtxX_LineFill);
0145     chain_in -> SetBranchAddress("vtxXE_LineFill", &in_vtxXE_LineFill);
0146     chain_in -> SetBranchAddress("vtxXStdDev_LineFill", &in_vtxXStdDev_LineFill);
0147     chain_in -> SetBranchAddress("vtxY_LineFill", &in_vtxY_LineFill);
0148     chain_in -> SetBranchAddress("vtxYE_LineFill", &in_vtxYE_LineFill);
0149     chain_in -> SetBranchAddress("vtxYStdDev_LineFill", &in_vtxYStdDev_LineFill);
0150     chain_in -> SetBranchAddress("run_nEvents", &in_run_nEvents);
0151     chain_in -> SetBranchAddress("job_index", &in_job_index);
0152     chain_in -> SetBranchAddress("file_total_event", &in_file_total_event);
0153 
0154 
0155     TH1F * quadrant_VtxX_1D = new TH1F("quadrant_VtxX",Form("quadrant_VtxX;X axis [cm];Entry"),100,X_range_l,X_range_r);
0156     TH1F * quadrant_VtxY_1D = new TH1F("quadrant_VtxY",Form("quadrant_VtxY;Y axis [cm];Entry"),100,Y_range_l,Y_range_r);
0157     TH1F * line_filled_mean_X_1D = new TH1F("line_filled_mean_X",Form("line_filled_mean_X;X axis [cm];Entry"),100,X_range_l,X_range_r);
0158     TH1F * line_filled_mean_Y_1D = new TH1F("line_filled_mean_Y",Form("line_filled_mean_Y;Y axis [cm];Entry"),100,Y_range_l,Y_range_r);
0159     
0160     TH1F * line_filled_stddev_X_1D = new TH1F("line_filled_stddev_X_1D",Form("line_filled_stddev_X_1D;StdDev, X axis [cm];Entry"),100,0,1 * 0.1);
0161     TH1F * line_filled_stddev_Y_1D = new TH1F("line_filled_stddev_Y_1D",Form("line_filled_stddev_Y_1D;StdDev, Y axis [cm];Entry"),100,0,1 * 0.1);
0162     
0163     // TH1F * line_filled_variance_X_1D = new TH1F("line_filled_variance_X_1D",Form("line_filled_variance_X_1D;Variance, X axis [cm];Entry"),100,0,1 * 0.1);    
0164     // TH1F * line_filled_variance_Y_1D = new TH1F("line_filled_variance_Y_1D",Form("line_filled_variance_Y_1D;Variance, Y axis [cm];Entry"),100,0,1 * 0.1);
0165     
0166     // TH1F * line_filled_covariance_1D = new TH1F("line_filled_covariance_1D",Form("line_filled_covariance_1D;Covariance [cm];Entry"),100,-0.06 * 0.1, 0.06 * 0.1);
0167 
0168 
0169     vector<double> evt_index_vec; evt_index_vec.clear();
0170     vector<double> evt_index_range_vec; evt_index_range_vec.clear();
0171 
0172     vector<double> quadrant_VtxX_vec; quadrant_VtxX_vec.clear();
0173     vector<double> quadrant_VtxY_vec; quadrant_VtxY_vec.clear();
0174     vector<double> quadrant_VtxXerror_vec; quadrant_VtxXerror_vec.clear();
0175     vector<double> quadrant_VtxYerror_vec; quadrant_VtxYerror_vec.clear();
0176 
0177     vector<double> line_filled_mean_X_vec; line_filled_mean_X_vec.clear();
0178     vector<double> line_filled_mean_Y_vec; line_filled_mean_Y_vec.clear();
0179     vector<double> line_filled_mean_Xerror_vec; line_filled_mean_Xerror_vec.clear();
0180     vector<double> line_filled_mean_Yerror_vec; line_filled_mean_Yerror_vec.clear();
0181 
0182     // note : for final determination, calculation or fitting
0183     vector<double> combined_index_vec; combined_index_vec.clear();
0184     vector<double> combined_index_range_vec; combined_index_range_vec.clear();
0185     vector<double> combined_X_vec; combined_X_vec.clear();
0186     vector<double> combined_Y_vec; combined_Y_vec.clear();
0187     vector<double> combined_XE_vec; combined_XE_vec.clear();
0188     vector<double> combined_YE_vec; combined_YE_vec.clear();
0189 
0190 
0191     vector<double> zero_vec; zero_vec.clear();
0192     vector<double> starteID_vec; starteID_vec.clear();
0193     vector<double> endeID_vec; endeID_vec.clear();
0194 
0195     SetsPhenixStyle();
0196     TCanvas * c1 = new TCanvas("c1","c1",950,800);
0197     c1 -> cd();
0198 
0199     TLatex * ltx = new TLatex();
0200     ltx->SetNDC();
0201     ltx->SetTextSize(0.045);
0202     ltx->SetTextAlign(31);
0203 
0204     TLatex * draw_text = new TLatex();
0205     draw_text -> SetNDC();
0206     draw_text -> SetTextSize(0.03);
0207 
0208     // int accumulate_runEvt = 0;
0209     int in_start_evt;
0210     int in_end_evt;
0211     int Previous_total_event;
0212 
0213     for (int i = 0; i < chain_in -> GetEntries(); i++)
0214     {
0215         chain_in -> GetEntry(i);
0216 
0217         // note : if changing from mm to cm is needed
0218 
0219         double quadrant_VtxX = in_vtxX_quadrant;
0220         double quadrant_VtxY = in_vtxY_quadrant;
0221         double quadrant_WindowX = in_vtxXE_quadrant;
0222         double quadrant_WindowY = in_vtxYE_quadrant;
0223         
0224         double line_filled_mean_X = in_vtxX_LineFill;
0225         double line_filled_mean_Y = in_vtxY_LineFill;
0226         double line_filled_meanE_X = in_vtxXE_LineFill;
0227         double line_filled_meanE_Y = in_vtxYE_LineFill;
0228         cout<<Form("unit : [cm]")<<", 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;
0229         cout<<Form("unit : [cm]")<<", 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;
0230 
0231         // todo : assuming all the files have the same number of events 
0232         // todo : and only the last file can have different number of events 
0233         // todo : the pros are, this one can handle if there is one file missing in the middle
0234         // todo : in addition, the number of events used in one file can have the freedom
0235         // todo : the cons are, it can not handle for the case if each file has different number of events
0236 
0237         in_start_evt = (i == chain_in -> GetEntries() - 1 && in_file_total_event > Previous_total_event) ? in_job_index * Previous_total_event : in_job_index * in_file_total_event;
0238         in_end_evt = in_start_evt + in_run_nEvents - 1;
0239 
0240         starteID_vec.push_back(in_start_evt);
0241         endeID_vec.push_back(in_end_evt);
0242 
0243         evt_index_vec.push_back((in_start_evt + in_end_evt)/2.);
0244         evt_index_range_vec.push_back((in_end_evt - in_start_evt)/2.);
0245         quadrant_VtxX_vec.push_back(quadrant_VtxX);
0246         quadrant_VtxY_vec.push_back(quadrant_VtxY);
0247         line_filled_mean_X_vec.push_back(line_filled_mean_X);
0248         line_filled_mean_Y_vec.push_back(line_filled_mean_Y);
0249         zero_vec.push_back(0);
0250 
0251         quadrant_VtxXerror_vec.push_back(quadrant_WindowX);
0252         quadrant_VtxYerror_vec.push_back(quadrant_WindowY);
0253         line_filled_mean_Xerror_vec.push_back(line_filled_meanE_X);
0254         line_filled_mean_Yerror_vec.push_back(line_filled_meanE_Y);
0255 
0256         combined_index_vec.push_back((in_start_evt + in_end_evt)/2.);
0257         combined_index_range_vec.push_back((in_end_evt - in_start_evt)/2.);
0258         combined_X_vec.push_back(quadrant_VtxX);
0259         combined_XE_vec.push_back(quadrant_WindowX);
0260         combined_Y_vec.push_back(quadrant_VtxY);
0261         combined_YE_vec.push_back(quadrant_WindowY);
0262 
0263         combined_index_vec.push_back((in_start_evt + in_end_evt)/2.);
0264         combined_index_range_vec.push_back((in_end_evt - in_start_evt)/2.);
0265         combined_X_vec.push_back(line_filled_mean_X);        
0266         combined_XE_vec.push_back(line_filled_meanE_X);
0267         combined_Y_vec.push_back(line_filled_mean_Y);
0268         combined_YE_vec.push_back(line_filled_meanE_Y);
0269 
0270 
0271 
0272 
0273         quadrant_VtxX_1D -> Fill(quadrant_VtxX);
0274         quadrant_VtxY_1D -> Fill(quadrant_VtxY);
0275         line_filled_mean_X_1D -> Fill(line_filled_mean_X);
0276         line_filled_mean_Y_1D -> Fill(line_filled_mean_Y);
0277 
0278         line_filled_stddev_X_1D -> Fill(in_vtxXStdDev_LineFill);
0279         line_filled_stddev_Y_1D -> Fill(in_vtxYStdDev_LineFill);
0280         // line_filled_variance_X_1D -> Fill(in_line_filled_variance_X);
0281         // line_filled_variance_Y_1D -> Fill(in_line_filled_variance_Y);
0282         // line_filled_covariance_1D -> Fill(in_line_filled_covariance);
0283 
0284         // accumulate_runEvt += in_run_nEvents;
0285 
0286         Previous_total_event = in_file_total_event;
0287     }
0288 
0289     for (int i = 0; i < starteID_vec.size(); i++)
0290     {
0291         cout<<"file : "<<i<<", start: "<<starteID_vec[i]<<", end: "<<endeID_vec[i]<<endl;
0292     }
0293 
0294     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]);
0295     g_quadrant_VtxX_index -> GetXaxis() -> SetNdivisions(505);
0296     g_quadrant_VtxX_index -> SetMarkerStyle(20);
0297     g_quadrant_VtxX_index -> SetMarkerSize(0.7);
0298     g_quadrant_VtxX_index -> SetMarkerColor(1);
0299     g_quadrant_VtxX_index -> GetXaxis() -> SetTitle("Event ID");
0300     g_quadrant_VtxX_index -> GetYaxis() -> SetTitle(Form("Average beam position in X axis [cm]"));
0301     g_quadrant_VtxX_index -> GetXaxis() -> SetLimits(-1000, g_quadrant_VtxX_index -> GetPointX(g_quadrant_VtxX_index -> GetN() - 1) + 10000);
0302     g_quadrant_VtxX_index -> GetYaxis() -> SetRangeUser(quadrant_VtxX_1D -> GetXaxis() -> GetXmin(),quadrant_VtxX_1D -> GetXaxis() -> GetXmax());
0303 
0304     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]);
0305     g_quadrant_VtxY_index -> GetXaxis() -> SetNdivisions(505);
0306     g_quadrant_VtxY_index -> SetMarkerStyle(20);
0307     g_quadrant_VtxY_index -> SetMarkerSize(0.7);
0308     g_quadrant_VtxY_index -> SetMarkerColor(1);
0309     g_quadrant_VtxY_index -> GetXaxis() -> SetTitle("Event ID");
0310     g_quadrant_VtxY_index -> GetYaxis() -> SetTitle(Form("Average beam position in Y axis [cm]"));
0311     g_quadrant_VtxY_index -> GetXaxis() -> SetLimits(-1000, g_quadrant_VtxY_index -> GetPointX(g_quadrant_VtxY_index -> GetN() - 1) + 10000);
0312     g_quadrant_VtxY_index -> GetYaxis() -> SetRangeUser(quadrant_VtxY_1D -> GetXaxis() -> GetXmin(),quadrant_VtxY_1D -> GetXaxis() -> GetXmax());
0313 
0314     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]);
0315     g_line_filled_mean_X_index -> GetXaxis() -> SetNdivisions(505);
0316     g_line_filled_mean_X_index -> SetMarkerStyle(20);
0317     g_line_filled_mean_X_index -> SetMarkerSize(0.7);
0318     g_line_filled_mean_X_index -> SetMarkerColor(4);
0319     g_line_filled_mean_X_index -> SetLineColor(4);
0320     g_line_filled_mean_X_index -> GetXaxis() -> SetTitle("Event ID");
0321     g_line_filled_mean_X_index -> GetYaxis() -> SetTitle(Form("Beam Position in X axis [cm]"));
0322 
0323     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]);
0324     g_line_filled_mean_Y_index -> GetXaxis() -> SetNdivisions(505);
0325     g_line_filled_mean_Y_index -> SetMarkerStyle(20);
0326     g_line_filled_mean_Y_index -> SetMarkerSize(0.7);
0327     g_line_filled_mean_Y_index -> SetMarkerColor(4);
0328     g_line_filled_mean_Y_index -> SetLineColor(4);
0329     g_line_filled_mean_Y_index -> GetXaxis() -> SetTitle("Event ID");
0330     g_line_filled_mean_Y_index -> GetYaxis() -> SetTitle(Form("Beam Position in Y axis [cm]"));
0331 
0332     
0333     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]);
0334     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);
0335     pol0_fit_X -> SetLineColor(2);
0336     pol0_fit_X -> SetLineStyle(9);
0337     pol0_fit_X -> SetLineWidth(1);
0338     combined_X_index_grE -> Fit(pol0_fit_X,"NQ");
0339 
0340     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]);
0341     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);
0342     pol0_fit_Y -> SetLineColor(2);
0343     pol0_fit_Y -> SetLineStyle(9);
0344     pol0_fit_Y -> SetLineWidth(1);
0345     combined_Y_index_grE -> Fit(pol0_fit_Y,"NQ");
0346 
0347     TLegend * leg = new TLegend(0.21,0.75,0.31,0.9);
0348     leg -> AddEntry(g_quadrant_VtxX_index, "Quadrant method", "lep");
0349     leg -> AddEntry(g_line_filled_mean_X_index, "2D tracklet fill method", "lep");
0350 
0351 
0352     c1 -> cd();
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("Numeric average vertex X: %.4f [cm]", ( vector_average(line_filled_mean_X_vec) + vector_average(quadrant_VtxX_vec) )/2. ));
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     c1 -> Print((output_directory + "/VtxX_index.pdf").c_str());
0359     c1 -> Clear();
0360 
0361     g_quadrant_VtxX_index -> Draw("AP");
0362     g_line_filled_mean_X_index -> Draw("p same");
0363     draw_text -> DrawLatex(0.21,0.71,Form("Fit vertex X: %.4f #pm%.4f [cm]", pol0_fit_X->GetParameter(0), pol0_fit_X->GetParError(0) ));
0364     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", label_text.c_str()));
0365     leg -> Draw("same");
0366     pol0_fit_X -> Draw("lsame");
0367     c1 -> Print((output_directory + "/VtxX_index_fit.pdf").c_str());
0368     c1 -> Clear();
0369 
0370     c1 -> cd();
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("Numeric average vertex Y: %.4f [cm]", ( vector_average(line_filled_mean_Y_vec) + vector_average(quadrant_VtxY_vec) )/2. ));
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     c1 -> Print((output_directory + "/VtxY_index.pdf").c_str());
0377     c1 -> Clear();
0378 
0379     g_quadrant_VtxY_index -> Draw("AP");
0380     g_line_filled_mean_Y_index -> Draw("p same");
0381     draw_text -> DrawLatex(0.21,0.71,Form("Fit vertex Y: %.4f #pm%.4f [cm]", pol0_fit_Y->GetParameter(0), pol0_fit_Y->GetParError(0) ));
0382     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", label_text.c_str()));
0383     leg -> Draw("same");
0384     pol0_fit_Y -> Draw("lsame");
0385     c1 -> Print((output_directory + "/VtxY_index_fit.pdf").c_str());
0386     c1 -> Clear();
0387 
0388     c1 -> cd();
0389     quadrant_VtxX_1D -> Draw("hist");
0390     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", label_text.c_str()));
0391     draw_text -> DrawLatex(0.2, 0.84, Form("Entry: %.0f", quadrant_VtxX_1D->GetEntries()));
0392     draw_text -> DrawLatex(0.2, 0.80, Form("Mean: %.4f", quadrant_VtxX_1D->GetMean()));
0393     draw_text -> DrawLatex(0.2, 0.76, Form("StdDev: %.4f", quadrant_VtxX_1D->GetStdDev()));
0394     c1 -> Print((output_directory + "/quadrant_VtxX_1D.pdf").c_str());
0395     c1 -> Clear();
0396 
0397     c1 -> cd();
0398     quadrant_VtxY_1D -> Draw("hist");
0399     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", label_text.c_str()));
0400     draw_text -> DrawLatex(0.2, 0.84, Form("Entry: %.0f", quadrant_VtxY_1D->GetEntries()));
0401     draw_text -> DrawLatex(0.2, 0.80, Form("Mean: %.4f", quadrant_VtxY_1D->GetMean()));
0402     draw_text -> DrawLatex(0.2, 0.76, Form("StdDev: %.4f", quadrant_VtxY_1D->GetStdDev()));
0403     c1 -> Print((output_directory + "/quadrant_VtxY_1D.pdf").c_str());
0404     c1 -> Clear();
0405 
0406     c1 -> cd();
0407     line_filled_mean_X_1D -> Draw("hist");
0408     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", label_text.c_str()));
0409     draw_text -> DrawLatex(0.2, 0.84, Form("Entry: %.0f", line_filled_mean_X_1D->GetEntries()));
0410     draw_text -> DrawLatex(0.2, 0.80, Form("Mean: %.4f", line_filled_mean_X_1D->GetMean()));
0411     draw_text -> DrawLatex(0.2, 0.76, Form("StdDev: %.4f", line_filled_mean_X_1D->GetStdDev()));
0412     c1 -> Print((output_directory + "/line_filled_mean_X_1D.pdf").c_str());
0413     c1 -> Clear();
0414 
0415     c1 -> cd();
0416     line_filled_mean_Y_1D -> Draw("hist");
0417     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", label_text.c_str()));
0418     draw_text -> DrawLatex(0.2, 0.84, Form("Entry: %.0f", line_filled_mean_Y_1D->GetEntries()));
0419     draw_text -> DrawLatex(0.2, 0.80, Form("Mean: %.4f", line_filled_mean_Y_1D->GetMean()));
0420     draw_text -> DrawLatex(0.2, 0.76, Form("StdDev: %.4f", line_filled_mean_Y_1D->GetStdDev()));
0421     c1 -> Print((output_directory + "/line_filled_mean_Y_1D.pdf").c_str());
0422     c1 -> Clear();
0423     
0424 
0425     c1 -> cd();
0426     line_filled_stddev_X_1D -> Draw("hist");
0427     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", label_text.c_str()));
0428     draw_text -> DrawLatex(0.2, 0.84, Form("Entry: %.0f", line_filled_stddev_X_1D->GetEntries()));
0429     draw_text -> DrawLatex(0.2, 0.80, Form("Mean: %.4f", line_filled_stddev_X_1D->GetMean()));
0430     draw_text -> DrawLatex(0.2, 0.76, Form("StdDev: %.4f", line_filled_stddev_X_1D->GetStdDev()));
0431     c1 -> Print((output_directory + "/line_filled_stddev_X_1D.pdf").c_str());
0432     c1 -> Clear();
0433 
0434     c1 -> cd();
0435     line_filled_stddev_Y_1D -> Draw("hist");
0436     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", label_text.c_str()));
0437     draw_text -> DrawLatex(0.2, 0.84, Form("Entry: %.0f", line_filled_stddev_Y_1D->GetEntries()));
0438     draw_text -> DrawLatex(0.2, 0.80, Form("Mean: %.4f", line_filled_stddev_Y_1D->GetMean()));
0439     draw_text -> DrawLatex(0.2, 0.76, Form("StdDev: %.4f", line_filled_stddev_Y_1D->GetStdDev()));
0440     c1 -> Print((output_directory + "/line_filled_stddev_Y_1D.pdf").c_str());
0441     c1 -> Clear();
0442 
0443     // c1 -> cd();
0444     // line_filled_variance_X_1D -> Draw("hist");
0445     // ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", label_text.c_str()));
0446     // draw_text -> DrawLatex(0.55, 0.84, Form("Entry: %.0f", line_filled_variance_X_1D->GetEntries()));
0447     // draw_text -> DrawLatex(0.55, 0.80, Form("Mean: %.4f", line_filled_variance_X_1D->GetMean()));
0448     // draw_text -> DrawLatex(0.55, 0.76, Form("StdDev: %.4f", line_filled_variance_X_1D->GetStdDev()));
0449     // c1 -> Print((output_directory + "/line_filled_variance_X_1D.pdf").c_str());
0450     // c1 -> Clear();
0451 
0452     // c1 -> cd();
0453     // line_filled_variance_Y_1D -> Draw("hist");
0454     // ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", label_text.c_str()));
0455     // draw_text -> DrawLatex(0.55, 0.84, Form("Entry: %.0f", line_filled_variance_Y_1D->GetEntries()));
0456     // draw_text -> DrawLatex(0.55, 0.80, Form("Mean: %.4f", line_filled_variance_Y_1D->GetMean()));
0457     // draw_text -> DrawLatex(0.55, 0.76, Form("StdDev: %.4f", line_filled_variance_Y_1D->GetStdDev()));
0458     // c1 -> Print((output_directory + "/line_filled_variance_Y_1D.pdf").c_str());
0459     // c1 -> Clear();
0460     
0461     // c1 -> cd();
0462     // line_filled_covariance_1D -> Draw("hist");
0463     // ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", label_text.c_str()));
0464     // draw_text -> DrawLatex(0.2, 0.84, Form("Entry: %.0f", line_filled_covariance_1D->GetEntries()));
0465     // draw_text -> DrawLatex(0.2, 0.80, Form("Mean: %.4f", line_filled_covariance_1D->GetMean()));
0466     // draw_text -> DrawLatex(0.2, 0.76, Form("StdDev: %.4f", line_filled_covariance_1D->GetStdDev()));
0467     // c1 -> Print((output_directory + "/line_filled_covariance_1D.pdf").c_str());
0468     // c1 -> Clear();
0469 
0470     double out_quadrant_avgX;
0471     double out_quadrant_avgXE;
0472     double out_quadrant_avgY;
0473     double out_quadrant_avgYE;
0474     double out_linefilled_avgX;
0475     double out_linefilled_avgXE;
0476     double out_linefilled_avgY;
0477     double out_linefilled_avgYE;
0478     double out_final_num_avgX;
0479     double out_final_num_avgY;
0480     double out_final_fit_avgX;
0481     double out_final_fit_avgY;
0482 
0483     TFile * file_out = new TFile (Form("%s/determined_avg_vtxXY.root",output_directory.c_str()),"recreate");
0484     TTree * tree_out = new TTree("tree","tree");
0485     tree_out -> Branch("quadrant_avgX",   &out_quadrant_avgX);
0486     tree_out -> Branch("quadrant_avgXE",    &out_quadrant_avgXE);
0487     tree_out -> Branch("quadrant_avgY",   &out_quadrant_avgY);
0488     tree_out -> Branch("quadrant_avgYE",    &out_quadrant_avgYE);
0489     tree_out -> Branch("linefilled_avgX",  &out_linefilled_avgX);
0490     tree_out -> Branch("linefilled_avgXE", &out_linefilled_avgXE);
0491     tree_out -> Branch("linefilled_avgY",  &out_linefilled_avgY);
0492     tree_out -> Branch("linefilled_avgYE", &out_linefilled_avgYE);
0493     tree_out -> Branch("final_num_avgX",   &out_final_num_avgX);
0494     tree_out -> Branch("final_num_avgY",   &out_final_num_avgY);
0495     tree_out -> Branch("final_fit_avgX",   &out_final_fit_avgX);
0496     tree_out -> Branch("final_fit_avgY",   &out_final_fit_avgY);
0497 
0498     out_quadrant_avgX = vector_average(quadrant_VtxX_vec); 
0499     out_quadrant_avgXE = vector_stddev(quadrant_VtxX_vec); 
0500     out_quadrant_avgY = vector_average(quadrant_VtxY_vec); 
0501     out_quadrant_avgYE = vector_stddev(quadrant_VtxY_vec); 
0502     
0503     out_linefilled_avgX = vector_average(line_filled_mean_X_vec);
0504     out_linefilled_avgXE = vector_stddev(line_filled_mean_X_vec);
0505     out_linefilled_avgY = vector_average(line_filled_mean_Y_vec);
0506     out_linefilled_avgYE = vector_stddev(line_filled_mean_Y_vec);
0507 
0508     out_final_num_avgX = ( vector_average(line_filled_mean_X_vec) + vector_average(quadrant_VtxX_vec) )/2.;
0509     out_final_num_avgY = ( vector_average(line_filled_mean_Y_vec) + vector_average(quadrant_VtxY_vec) )/2.;
0510     out_final_fit_avgX = pol0_fit_X->GetParameter(0);
0511     out_final_fit_avgY = pol0_fit_Y->GetParameter(0);
0512 
0513     tree_out -> Fill();
0514     file_out -> cd();
0515     // tree_out->SetDirectory(file_out);
0516     tree_out -> Write("", TObject::kOverwrite);
0517     file_out -> Close();
0518 
0519 
0520     cout<<"unit : [cm]"<<endl;
0521     // cout<<"line filled covariance info : "<<line_filled_covariance_1D->GetMean()<<" "<<line_filled_covariance_1D->GetStdDev()<<endl;
0522     // cout<<"line filled variance X info : "<<line_filled_variance_X_1D->GetMean()<<" "<<line_filled_variance_X_1D->GetStdDev()<<endl;
0523     // cout<<"line filled variance Y info : "<<line_filled_variance_Y_1D->GetMean()<<" "<<line_filled_variance_Y_1D->GetStdDev()<<endl;
0524     cout<<"line filled stdDev X   info : " <<line_filled_stddev_X_1D->GetMean()<<" "<<line_filled_stddev_X_1D->GetStdDev()<<endl;
0525     cout<<"line filled stdDev Y   info : " <<line_filled_stddev_Y_1D->GetMean()<<" "<<line_filled_stddev_Y_1D->GetStdDev()<<endl;
0526     cout<<endl;
0527     cout<<"unit : [cm]"<<endl;
0528     cout<<"final average vertex XY should be used : "<<endl;
0529     cout<<"line filled X : "<<vector_average(line_filled_mean_X_vec)<<" +/- "<<vector_stddev(line_filled_mean_X_vec)<<endl;
0530     cout<<"line filled Y : "<<vector_average(line_filled_mean_Y_vec)<<" +/- "<<vector_stddev(line_filled_mean_Y_vec)<<endl;
0531     cout<<"quadrant X : "<<vector_average(quadrant_VtxX_vec)<<" +/- "<<vector_stddev(quadrant_VtxX_vec)<<endl;
0532     cout<<"quadrant Y : "<<vector_average(quadrant_VtxY_vec)<<" +/- "<<vector_stddev(quadrant_VtxY_vec)<<endl;
0533     cout<<"avg: {"<< out_final_num_avgX <<Form(", ")<< out_final_num_avgY <<Form("} [cm]")<<endl;
0534     cout<<"Fit avg: {"<< out_final_fit_avgX <<Form(", ")<< out_final_fit_avgY <<Form("} [cm]")<<endl;
0535 
0536     return 0;
0537 }