analysis_INTT_CW/DST_MC/INTTXYvtx.h

0001 #ifndef INTTXYvtx_h
0002 #define INTTXYvtx_h
0003
0004 #include <filesystem>
0005 #include "../private_cluster_gen/InttConversion_new.h"
0006 #include "../private_cluster_gen/InttClustering.h"
0007 #include "../sigmaEff.h"
0008 #include "../sPhenixStyle.C"
0009 #include "gaus_func.h"
0010 #include "ana_map_folder/ana_map_v1.h"
0011 namespace ana_map_version = ANA_MAP_V3;
0012
0013 // note : this class mainly focus on two things
0014 // note : 1. find a single vertex for the whole run
0015 // note :     a. the functions prepared for this purpose are : ProcessEvt(), MacroVTXxyCorrection(), MacroVTXxyCorrection_new(), GetFinalVTXxy(), MacrovtxSquare()
0016 // note : 2. find the vertex for each event
0017 // note :     a. the functions prepared for this purpose are : ProcessEvt()
0018
0019 struct type_pos{
0020     double x;
0021     double y;
0022 };
0023
0024 double cos_func(double *x, double *par)
0025 {
0026     return -1 * par[0] * cos(par[1] * (x[0] - par[2])) + par[3];
0027 }
0028
0029 // double gaus_func(double *x, double *par)
0030 // {
0031 //     // note : par[0] : size
0032 //     // note : par[1] : mean
0033 //     // note : par[2] : width
0034 //     // note : par[3] : offset
0035 //     return par[0] * TMath::Gaus(x[0],par[1],par[2]) + par[3];
0036 // }
0037
0038 class INTTXYvtx {
0039     public:
0040         INTTXYvtx(string run_type, string out_folder_directory, pair<double,double> beam_origin, int geo_mode_id, double phi_diff_cut = 0.11, pair<double, double> DCA_cut = {-1,1}, int N_clu_cutl = 20, int N_clu_cut = 10000, bool draw_event_display = true, double peek = 3.32405, double angle_diff_new_l = 0.0, double angle_diff_new_r = 3, bool print_message_opt = true);
0041         INTTXYvtx(string run_type, string out_folder_directory);
0042         void ProcessEvt(
0043             int event_i,
0044             vector<clu_info> temp_sPH_inner_nocolumn_vec,
0045             vector<clu_info> temp_sPH_outer_nocolumn_vec,
0046             vector<vector<double>> temp_sPH_nocolumn_vec,
0047             vector<vector<double>> temp_sPH_nocolumn_rz_vec,
0048             int NvtxMC,
0049             double TrigZvtxMC,
0050             bool PhiCheckTag,
0051             Long64_t bco_full,
0052             int is_min_bias, // note : for data and MC
0053             int is_min_bias_wozdc // note : for data only
0054         );
0055         virtual void ClearEvt();
0056         virtual void PrintPlots();
0057         virtual void EndRun();
0058
0059         unsigned long GetVecNele();
0060         void MacroVTXxyCorrection(bool phi_correction, bool calculate_XY, int N_trial, double fit_range_l, double fit_range_r);
0061         void MacroVTXxyCorrection_new(double fit_range_l, double fit_range_r, int N_trial, pair<double,double> vertex_in=make_pair(-999,-999));
0062         virtual vector<pair<double,double>> MacroVTXSquare(double length, int N_trial, bool draw_plot_opt = true);
0063         pair<double,double> GetFinalVTXxy();
0064         pair<vector<TH2F *>, vector<TH1F*>> GetHistFinal();
0065         void TH1F_FakeClone(TH1F*hist_in, TH1F*hist_out);
0066         void TH2F_FakeClone(TH2F*hist_in, TH2F*hist_out);
0067         void TH2F_FakeRebin(TH2F*hist_in, TH2F*hist_out);
0068         vector<pair<double,double>> FillLine_FindVertex(pair<double,double> window_center, double segmentation = 0.005, double window_width = 5.0, int N_bins = 100, bool draw_plot = true);
0069         vector<double> LineFill_bkgrm_Get_covariance() {return Get_covariance_TH2(xy_hist_bkgrm);};
0070
0071         // note : for the delta phi study
0072         void RunDeltaPhiStudy();
0073         TProfile * Get_final_angle_diff_inner_phi_degree_coarseX_peak_profile() {return final_angle_diff_inner_phi_degree_coarseX_peak_profile;};
0074         TH1F * Get_final_angle_diff_inner_phi_degree_coarseX_peak_profile_hist() {
0075
0076             TH1F * angle_diff_correction_degree_hist = new TH1F(
0077                 "",
0078                 "",
0079                 final_angle_diff_inner_phi_degree_coarseX_peak_profile -> GetNbinsX(),
0080                 final_angle_diff_inner_phi_degree_coarseX_peak_profile -> GetXaxis() -> GetXmin(),
0081                 final_angle_diff_inner_phi_degree_coarseX_peak_profile -> GetXaxis() -> GetXmax()
0082             );
0083
0084             for (int i = 0; i < final_angle_diff_inner_phi_degree_coarseX_peak_profile -> GetNbinsX(); i++) {
0085                 angle_diff_correction_degree_hist -> SetBinContent(i+1, final_angle_diff_inner_phi_degree_coarseX_peak_profile -> GetBinContent(i+1));
0086             }
0087
0088             return angle_diff_correction_degree_hist;
0089         }
0090         TH1F * Get_final_DCA_mm_inner_phi_degree_coarseX_peak_profile_hist() {
0091
0092             TH1F * DCA_mm_correction_degree_hist = new TH1F(
0093                 "",
0094                 "",
0095                 final_DCA_mm_inner_phi_degree_coarseX_peak_profile -> GetNbinsX(),
0096                 final_DCA_mm_inner_phi_degree_coarseX_peak_profile -> GetXaxis() -> GetXmin(),
0097                 final_DCA_mm_inner_phi_degree_coarseX_peak_profile -> GetXaxis() -> GetXmax()
0098             );
0099
0100             for (int i = 0; i < final_DCA_mm_inner_phi_degree_coarseX_peak_profile -> GetNbinsX(); i++) {
0101                 DCA_mm_correction_degree_hist -> SetBinContent(i+1, final_DCA_mm_inner_phi_degree_coarseX_peak_profile -> GetBinContent(i+1));
0102             }
0103
0104             return DCA_mm_correction_degree_hist;
0105         }
0106
0107     protected:
0108         TCanvas * c1;
0109         TLatex * ltx;
0110         TLatex * draw_text;
0111         TF1 * cos_fit;
0112         TF1 * gaus_fit;
0113         TF1 * gaus_fit_angle_diff;
0114         TF1 * gaus_fit_MC;
0115         TF1 * correlation_Line;
0116         TF1 * bkg_fit_pol2;
0117         TF1 * draw_pol2_line;
0118         TF1 * correlation_fit_MC;
0119         TF1 * gaus_pol2_fit;
0120
0121         TF1 * horizontal_fit_inner;
0122         TF1 * horizontal_fit_angle_diff_inner;
0123         TF1 * horizontal_fit_outer;
0124         TF1 * horizontal_fit_angle_diff_outer;
0125         TF1 * fit_rz;
0126         TF1 * gaus_pol1_fit;
0127         TF1 * draw_gaus_line;
0128         TF1 * draw_pol1_line;
0129         TF1 * d_gaus_pol1_fit;
0130         TF1 * draw_d_gaus;
0131         // TF1 * track_pol1_fit;
0132
0133         TLine * track_line;
0134         TLine * coord_line;
0135         TLine * ladder_line;
0136
0137         TH2F * angle_correlation;
0138         TH2F * angle_diff_DCA_dist;
0139         TH1F * angle_diff;
0140         TH1F * angle_diff_new;
0141         TH1F * angle_diff_new_bkg_remove;
0142         TH1F * angle_diff_new_bkg_remove_final;
0143         TH1F * angle_diff_radian;
0144         TH2F * inner_pos_xy;
0145         TH2F * outer_pos_xy;
0146         TH2F * inner_outer_pos_xy;
0147         TH2F * DCA_point;
0148         TH1F * DCA_distance;
0149         TH2F * N_cluster_correlation;
0150         TH2F * N_cluster_correlation_close;
0151
0152         TH2F * DCA_distance_inner_X;
0153         TH2F * DCA_distance_inner_Y;
0154         TH2F * DCA_distance_outer_X;
0155         TH2F * DCA_distance_outer_Y;
0156
0157         TH2F *     DCA_distance_inner_phi;
0158         TH2F *     DCA_distance_inner_phi_peak;
0159         TProfile * DCA_distance_inner_phi_peak_profile;
0160         TH2F *     DCA_distance_outer_phi;
0161         TH2F *     DCA_distance_outer_phi_peak;
0162         TProfile * DCA_distance_outer_phi_peak_profile;
0163
0164         TH2F *     angle_diff_inner_phi;
0165         TH2F *     angle_diff_inner_phi_peak;
0166         TProfile * angle_diff_inner_phi_peak_profile;
0167         TH2F *     angle_diff_outer_phi;
0168         TH2F *     angle_diff_outer_phi_peak;
0169         TProfile * angle_diff_outer_phi_peak_profile;
0170
0171         // note : those are for shifting the peak to the center
0172         TH2F *     final_DCA_mm_inner_phi_degree_coarseX;
0173         TH2F *     final_DCA_mm_inner_phi_degree_coarseX_peak;
0174         TProfile * final_DCA_mm_inner_phi_degree_coarseX_peak_profile;
0175
0176         TH2F *     final_DCA_cm_inner_phi_radian_coarseX;
0177         TH2F *     final_DCA_cm_inner_phi_radian_coarseX_peak;
0178         TProfile * final_DCA_cm_inner_phi_radian_coarseX_peak_profile;
0179
0180         TH2F *     final_angle_diff_inner_phi_radian_coarseX;
0181         TH2F *     final_angle_diff_inner_phi_radian_coarseX_peak;
0182         TProfile * final_angle_diff_inner_phi_radian_coarseX_peak_profile;
0183         TH2F *     final_angle_diff_inner_phi_degree_coarseX;
0184         TH2F *     final_angle_diff_inner_phi_degree_coarseX_peak;
0185         TProfile * final_angle_diff_inner_phi_degree_coarseX_peak_profile;
0186         TH1F *     final_angle_diff_radian_before;
0187         TH1F *     final_angle_diff_radian_post;
0188         TH2F *     final_angle_diff_radian_DCA_2D_before;
0189         TH2F *     final_angle_diff_radian_DCA_2D_post;
0190         TH1F *     final_DCA_cm_before;
0191         TH1F *     final_DCA_cm_post;
0192
0193         TGraph * angle_diff_inner_phi_peak_profile_graph;
0194         TGraph * angle_diff_outer_phi_peak_profile_graph;
0195         TGraph * DCA_distance_inner_phi_peak_profile_graph;
0196         TGraph * DCA_distance_outer_phi_peak_profile_graph;
0197
0198         // note : it's for the geometry correction
0199         TH2F * DCA_distance_inner_phi_peak_final;
0200         TH2F * angle_diff_inner_phi_peak_final;
0201         TH2F * DCA_distance_outer_phi_peak_final;
0202         TH2F * angle_diff_outer_phi_peak_final;
0203
0204         TH2F * xy_hist;
0205         TH2F * xy_hist_bkgrm;
0206
0207         // note : to keep the cluster pair information
0208         // note : this is the vector for the whole run, not event by event
0209         vector<pair<type_pos,type_pos>> cluster_pair_vec;
0210         double Clus_InnerPhi_Offset;
0211         double Clus_OuterPhi_Offset;
0212         double Clus_InnerPhi_Offset_radian;
0213         double Clus_OuterPhi_Offset_radian;
0214         double vtxXcorrection;
0215         double vtxYcorrection;
0216         double current_vtxX;
0217         double current_vtxY;
0218         int    MacroVTXxyCorrection_new_function_call_count;
0219
0220         // note : not used in the INTTXYvtx class
0221         // pair<double,double> zvtx_QA_width;
0222         // double zvtx_QA_ratio;
0223         int zvtx_cal_require = 15;
0224         map<pair<int,int>,int> axis_quadrant_map = {
0225             {{1,1},0},
0226             {{-1,1},1},
0227             {{-1,-1},2},
0228             {{1,-1},3}
0229         };
0230
0231         double radian_correction = 180./M_PI;
0232
0233         pair<double,double> beam_origin;
0234         pair<double,double> DCA_cut;
0235         string out_folder_directory;
0236         string plot_text;
0237         string run_type;
0238         double phi_diff_cut;
0239         double peek;
0240         double angle_diff_new_l;
0241         double angle_diff_new_r;
0242         bool draw_event_display;
0243         bool print_message_opt;
0244         long total_NClus;
0245         int geo_mode_id;
0246         int N_clu_cutl;
0247         int N_clu_cut;
0248
0249         // TFile * file_out;
0250         // TTree * tree_out;
0251         // double out_quadrant_corner_X;
0252         // double out_quadrant_corner_Y;
0253         // double out_quadrant_center_X;
0254         // double out_quadrant_center_X;
0255         // double out_line_filled_mean_X;
0256         // double out_line_filled_mean_Y;
0257         // double out_line_filled_stddev_X;
0258         // double out_line_filled_stddev_Y;
0259
0260         void Init();
0261         void InitHist();
0262         void InitCanvas();
0263         virtual void InitTreeOut();
0264         void InitRest();
0265         void InitGraph();
0266
0267         pair<double,double> GetVTXxyCorrection(bool phi_correction, bool calculate_XY, int trial_index, double fit_range_l, double fit_range_r);
0268         vector<double> GetVTXxyCorrection_new(int trial_index);
0269         void PrintPlotsVTXxy(string sub_out_folder_name, int print_option = 0);
0270         void ClearHist(int print_option = 0);
0271
0272         double get_dist_offset(TH1F * hist_in, int check_N_bin); // note : check_N_bin 1 to N bins of hist
0273         double get_radius(double x, double y);
0274         void Characterize_Pad(TPad *pad, float left = 0.15, float right = 0.1, float top = 0.1, float bottom = 0.12, bool set_logY = false, int setgrid_bool = 0);
0275         std::vector<double> calculateDistanceAndClosestPoint(double x1, double y1, double x2, double y2, double target_x, double target_y);
0276         double calculateAngleBetweenVectors(double x1, double y1, double x2, double y2, double targetX, double targetY);
0277         void TH2F_threshold(TH2F * hist, double threshold);
0278         void TH2F_threshold_advanced(TH2F * hist, double threshold); // note : for each column
0279         void TH2F_threshold_advanced_row(TH2F * hist, double threshold); // note : for each row
0280         void TH2F_threshold_advanced_2(TH2F * hist, double threshold);
0281         pair<type_pos,type_pos> GetTangentPointsAtCircle(double CenterX, double CenterY, double R, double XX, double YY);
0282         type_pos PolarToCartesian(double radius, double angleDegrees);
0283         pair<type_pos,type_pos> findIntersectionPoints(double A, double mu, double sigma, double B, double C);
0284         vector<double> find_Ngroup(TH1F * hist_in);
0285         void Hist_1D_bkg_remove(TH1F * hist_in, double factor);
0286         void DrawTGraphErrors(vector<double> x_vec, vector<double> y_vec, vector<double> xE_vec, vector<double> yE_vec, string output_directory, vector<string> plot_name);
0287         void Draw2TGraph(vector<double> x1_vec, vector<double> y1_vec, vector<double> x2_vec, vector<double> y2_vec, string output_directory, vector<string> plot_name);
0288         void DrawTH2F(TH2F * hist_in, string output_directory, vector<string> plot_name);
0289         pair<double,double> GetCircleAngle(double fit_range_l, double fit_range_r);
0290         vector<double> SumTH2FColumnContent(TH2F * hist_in);
0291         vector<double> SumTH2FColumnContent_row(TH2F * hist_in);
0292         vector<pair<double,double>> Get4vtx(pair<double,double> origin, double length);
0293         vector<pair<double,double>> Get4vtxAxis(pair<double,double> origin, double length);
0294         vector<double> subMacroVTXxyCorrection(int test_index, int trial_index, bool draw_plot_opt);
0295         virtual void subMacroPlotWorking(bool phi_correction, double cos_fit_rangel, double cos_fit_ranger, double guas_fit_range);
0296         pair<double,double> Get_possible_zvtx(double rvtx, vector<double> p0, vector<double> p1);
0297         int find_quadrant(pair<double,double> Origin, pair<double,double> check_point);
0298         vector<double> Get_covariance_TH2(TH2F * hist_in);
0299
0300         TH1F * PrintPlots_bkgrm_helper(TH1F * hist_in, double signal_region);
0301
0302         // note : from the INTTXYvtxEvt.h
0303         void TH2FSampleLineFill(TH2F * hist_in, double segmentation, std::pair<double,double> inner_clu, std::pair<double,double> outer_clu);
0304 };
0305
0306 INTTXYvtx::INTTXYvtx(string run_type, string out_folder_directory, pair<double,double> beam_origin, int geo_mode_id, double phi_diff_cut, pair<double,double> DCA_cut, int N_clu_cutl, int N_clu_cut, bool draw_event_display, double peek, double angle_diff_new_l, double angle_diff_new_r, bool print_message_opt)
0307 :run_type(run_type), out_folder_directory(out_folder_directory), beam_origin(beam_origin), geo_mode_id(geo_mode_id), peek(peek), N_clu_cut(N_clu_cut), N_clu_cutl(N_clu_cutl), phi_diff_cut(phi_diff_cut), DCA_cut(DCA_cut), draw_event_display(draw_event_display), angle_diff_new_l(angle_diff_new_l), angle_diff_new_r(angle_diff_new_r), print_message_opt(print_message_opt)
0308 {
0309     SetsPhenixStyle();
0310     if (filesystem::exists(out_folder_directory.c_str()) == false) {system(Form("mkdir %s",out_folder_directory.c_str()));}
0311     gErrorIgnoreLevel = kWarning; // note : To not print the "print plot info."
0312
0313     Init();
0314     plot_text = (run_type == "MC") ? "Simulation" : "Work-in-progress";
0315
0316     cluster_pair_vec.clear();
0317     vtxXcorrection = 0;
0318     vtxYcorrection = 0;
0319
0320     current_vtxX = beam_origin.first;
0321     current_vtxY = beam_origin.second;
0322
0323     MacroVTXxyCorrection_new_function_call_count = 0;
0324 }
0325
0326 INTTXYvtx::INTTXYvtx(string run_type, string out_folder_directory)
0327 :run_type(run_type), out_folder_directory(out_folder_directory)
0328 {
0329     SetsPhenixStyle();
0330     if (filesystem::exists(out_folder_directory.c_str()) == false) {system(Form("mkdir %s",out_folder_directory.c_str()));}
0331     gErrorIgnoreLevel = kWarning; // note : To not print the "print plot info."
0332
0333     InitCanvas();
0334     InitRest();
0335     plot_text = (run_type == "MC") ? "Simulation" : "Work-in-progress";
0336 }
0337
0338 void INTTXYvtx::Init()
0339 {
0340     InitHist();
0341     InitCanvas();
0342     InitTreeOut();
0343     InitRest();
0344     InitGraph();
0345 }
0346
0347 void INTTXYvtx::InitTreeOut()
0348 {
0349     // file_out = new TFile(Form("%s/run_XY_tree.root",out_folder_directory.c_str()),"RECREATE");
0350     // file_out -> cd();
0351
0352     // tree_out = new TTree("tree", "tree avg VtxXY");
0353     // tree_out -> Branch("quadrant_corner_X",&out_quadrant_corner_X);
0354     // tree_out -> Branch("quadrant_corner_Y",&out_quadrant_corner_Y);
0355     // tree_out -> Branch("quadrant_center_X",&out_quadrant_center_X);
0356     // tree_out -> Branch("quadrant_center_X",&out_quadrant_center_X);
0357     // tree_out -> Branch("line_filled_mean_X",&out_line_filled_mean_X);
0358     // tree_out -> Branch("line_filled_mean_Y",&out_line_filled_mean_Y);
0359     // tree_out -> Branch("line_filled_stddev_X",&out_line_filled_stddev_X);
0360     // tree_out -> Branch("line_filled_stddev_Y",&out_line_filled_stddev_Y);
0361
0362     return;
0363 }
0364
0365 void INTTXYvtx::InitGraph()
0366 {
0367     angle_diff_inner_phi_peak_profile_graph = new TGraph();
0368     angle_diff_outer_phi_peak_profile_graph = new TGraph();
0369     DCA_distance_inner_phi_peak_profile_graph = new TGraph();
0370     DCA_distance_outer_phi_peak_profile_graph = new TGraph();
0371 }
0372
0373 void INTTXYvtx::InitRest()
0374 {
0375     ltx = new TLatex();
0376     ltx->SetNDC();
0377     ltx->SetTextSize(0.045);
0378     ltx->SetTextAlign(31);
0379
0380     draw_text = new TLatex();
0381     draw_text -> SetNDC();
0382     draw_text -> SetTextSize(0.03);
0383
0384     cos_fit = new TF1("cos_fit",cos_func,0,360,4);
0385     cos_fit -> SetParNames("[A]", "[B]", "[C]", "[D]");
0386     cos_fit -> SetLineColor(2);
0387
0388     gaus_fit = new TF1("gaus_fit",gaus_func,0,360,4);
0389     gaus_fit -> SetLineColor(4);
0390     gaus_fit -> SetLineWidth(1);
0391     gaus_fit -> SetParNames("size", "mean", "width", "offset");
0392     gaus_fit -> SetNpx(1000);
0393
0394     gaus_pol2_fit = new TF1("gaus_pol2_fit", gaus_pol2_func, -5, 5, 7);
0395     gaus_pol2_fit -> SetLineWidth(2);
0396     gaus_pol2_fit -> SetLineColor(2);
0397     gaus_pol2_fit -> SetNpx(1000);
0398
0399     gaus_fit_MC = new TF1("gaus_fit_MC",gaus_func,-10,10,4);
0400     gaus_fit_MC -> SetLineColor(2);
0401     gaus_fit_MC -> SetLineWidth(2);
0402     gaus_fit_MC -> SetNpx(1000);
0403
0404     correlation_fit_MC = new TF1("correlation_fit_MC","pol1",-10,10);
0405     correlation_fit_MC -> SetLineColor(2);
0406     correlation_fit_MC -> SetLineWidth(2);
0407
0408     correlation_Line = new TF1("correlation_Line","pol1",-10000,10000);
0409     correlation_Line -> SetLineColor(2);
0410     correlation_Line -> SetLineWidth(2);
0411     correlation_Line -> SetParameters(0,1);
0412     // gaus_fit_MC -> SetNpx(1000);
0413
0414
0415     gaus_fit_angle_diff = new TF1("gaus_fit_angle_diff",gaus_func,0,360,4);
0416     gaus_fit_angle_diff -> SetLineColor(4);
0417     gaus_fit_angle_diff -> SetLineWidth(2);
0418     gaus_fit_angle_diff -> SetParNames("size", "mean", "width", "offset");
0419     gaus_fit_angle_diff -> SetNpx(1000);
0420
0421     horizontal_fit_inner = new TF1("horizontal_fit_inner","pol0",-360,360);
0422     horizontal_fit_inner -> SetLineWidth(2);
0423     horizontal_fit_inner -> SetLineColor(2);
0424
0425     horizontal_fit_angle_diff_inner = new TF1("horizontal_fit_angle_diff_inner","pol0",-360,360);
0426     horizontal_fit_angle_diff_inner -> SetLineWidth(2);
0427     horizontal_fit_angle_diff_inner -> SetLineColor(2);
0428
0429     horizontal_fit_outer = new TF1("horizontal_fit_outer","pol0",-360,360);
0430     horizontal_fit_outer -> SetLineWidth(2);
0431     horizontal_fit_outer -> SetLineColor(2);
0432
0433     horizontal_fit_angle_diff_outer = new TF1("horizontal_fit_angle_diff_outer","pol0",-360,360);
0434     horizontal_fit_angle_diff_outer -> SetLineWidth(2);
0435     horizontal_fit_angle_diff_outer -> SetLineColor(2);
0436
0437     coord_line = new TLine();
0438     coord_line -> SetLineWidth(1);
0439     coord_line -> SetLineColor(16);
0440     coord_line -> SetLineStyle(2);
0441
0442     track_line = new TLine();
0443     track_line -> SetLineWidth(1);
0444     track_line -> SetLineColor(38);
0445     // track_line -> SetLineColorAlpha(38,0.5);
0446
0447     ladder_line = new TLine();
0448     ladder_line -> SetLineWidth(1);
0449
0450     fit_rz = new TF1("fit_rz","pol1",-500,500);
0451
0452     gaus_pol1_fit = new TF1("gaus_pol1_fit",gaus_pol1_func,-3,3,5);
0453     gaus_pol1_fit -> SetLineWidth(2);
0454     gaus_pol1_fit -> SetLineColor(2);
0455     gaus_pol1_fit -> SetNpx(1000);
0456
0457     draw_gaus_line = new TF1("draw_gaus_line",gaus_func,-3,3,4);
0458     draw_gaus_line -> SetLineWidth(2);
0459     draw_gaus_line -> SetLineStyle(9);
0460     draw_gaus_line -> SetLineColor(6);
0461     draw_gaus_line -> SetNpx(1000);
0462
0463     draw_pol1_line = new TF1("draw_pol1_line","pol1",-3,3);
0464     draw_pol1_line -> SetLineStyle(9);
0465     draw_pol1_line -> SetLineWidth(2);
0466     draw_pol1_line -> SetLineColor(3);
0467
0468     bkg_fit_pol2 = new TF1("bkg_fit_pol2",bkg_pol2_func,-5,5,5);
0469     bkg_fit_pol2 -> SetLineWidth(2);
0470     bkg_fit_pol2 -> SetLineColor(2);
0471     bkg_fit_pol2 -> SetNpx(1000);
0472
0473     draw_pol2_line = new TF1("draw_pol2_line",full_pol2_func,-5,5,4);
0474     draw_pol2_line -> SetLineWidth(2);
0475     draw_pol2_line -> SetLineColor(2);
0476     draw_pol2_line -> SetNpx(1000);
0477
0478
0479
0480
0481     d_gaus_pol1_fit  = new TF1("d_gaus_pol1_fit",d_gaus_pol1_func,-3,3,7);
0482     d_gaus_pol1_fit -> SetLineWidth(2);
0483     d_gaus_pol1_fit -> SetLineColor(2);
0484     d_gaus_pol1_fit -> SetNpx(1000);
0485
0486     draw_d_gaus = new TF1("draw_d_gaus", d_gaus_func,-3,3,5);
0487     draw_d_gaus -> SetLineStyle(9);
0488     draw_d_gaus -> SetLineWidth(2);
0489     draw_d_gaus -> SetLineColor(3);
0490     draw_d_gaus -> SetNpx(1000);
0491
0492
0493     // track_pol1_fit = new TF1("track_pol1_fit","pol1",-500,500);
0494 }
0495
0496 void INTTXYvtx::InitCanvas()
0497 {
0498     c1 = new TCanvas("","",950,800);
0499     c1 -> cd();
0500 }
0501
0502 void INTTXYvtx::InitHist()
0503 {
0504     angle_correlation = new TH2F("","angle_correlation",361,0,361,361,0,361);
0505     angle_correlation -> SetStats(0);
0506     angle_correlation -> GetXaxis() -> SetTitle("Inner cluster #Phi [degree]");
0507     angle_correlation -> GetYaxis() -> SetTitle("Outer cluster #Phi [degree]");
0508     angle_correlation -> GetXaxis() -> SetNdivisions(505);
0509
0510     angle_diff_DCA_dist = new TH2F("","angle_diff_DCA_dist",100,-1.5,1.5,100,-3.5,3.5);
0511     angle_diff_DCA_dist -> SetStats(0);
0512     angle_diff_DCA_dist -> GetXaxis() -> SetTitle("Inner - Outer [degree]");
0513     angle_diff_DCA_dist -> GetYaxis() -> SetTitle("DCA distance [mm]");
0514     angle_diff_DCA_dist -> GetXaxis() -> SetNdivisions(505);
0515
0516     angle_diff = new TH1F("angle_diff","angle_diff",100,0,5);
0517     angle_diff -> SetStats(0);
0518     angle_diff -> GetXaxis() -> SetTitle("|Inner - Outer| [degree]");
0519     angle_diff -> GetYaxis() -> SetTitle("Entry");
0520     angle_diff -> GetXaxis() -> SetNdivisions(505);
0521
0522     angle_diff_radian = new TH1F("","angle_diff_radian;#Delta#phi (Inner - Outer) [radian];Entry",100,-0.052359878,0.052359878);
0523     angle_diff_radian -> SetStats(0);
0524     angle_diff_radian -> GetXaxis() -> SetNdivisions(505);
0525
0526     angle_diff_new = new TH1F("angle_diff_new","angle_diff_new",100, angle_diff_new_l, angle_diff_new_r);
0527     angle_diff_new -> SetStats(0);
0528     angle_diff_new -> GetXaxis() -> SetTitle("|Inner - Outer| [degree]");
0529     angle_diff_new -> GetYaxis() -> SetTitle("Entry");
0530     angle_diff_new -> GetXaxis() -> SetNdivisions(505);
0531
0532     angle_diff_new_bkg_remove_final = new TH1F("","angle_diff_new_bkg_remove_final",100, angle_diff_new_l, angle_diff_new_r);
0533     angle_diff_new_bkg_remove_final -> SetStats(0);
0534     angle_diff_new_bkg_remove_final -> GetXaxis() -> SetTitle("|Inner - Outer| [degree]");
0535     angle_diff_new_bkg_remove_final -> GetYaxis() -> SetTitle("Entry");
0536     angle_diff_new_bkg_remove_final -> GetXaxis() -> SetNdivisions(505);
0537
0538     inner_pos_xy = new TH2F("","inner_pos_xy",360,-100,100,360,-100,100);
0539     inner_pos_xy -> SetStats(0);
0540     inner_pos_xy -> GetXaxis() -> SetTitle("X axis [mm]");
0541     inner_pos_xy -> GetYaxis() -> SetTitle("Y axis [mm]");
0542     inner_pos_xy -> GetXaxis() -> SetNdivisions(505);
0543
0544     outer_pos_xy = new TH2F("","outer_pos_xy",360,-150,150,360,-150,150);
0545     outer_pos_xy -> SetStats(0);
0546     outer_pos_xy -> GetXaxis() -> SetTitle("X axis [mm]");
0547     outer_pos_xy -> GetYaxis() -> SetTitle("Y axis [mm]");
0548     outer_pos_xy -> GetXaxis() -> SetNdivisions(505);
0549
0550     inner_outer_pos_xy = new TH2F("","inner_outer_pos_xy",360,-150,150,360,-150,150);
0551     inner_outer_pos_xy -> SetStats(0);
0552     inner_outer_pos_xy -> GetXaxis() -> SetTitle("X axis [mm]");
0553     inner_outer_pos_xy -> GetYaxis() -> SetTitle("Y axis [mm]");
0554     inner_outer_pos_xy -> GetXaxis() -> SetNdivisions(505);
0555
0556     DCA_point = new TH2F("","DCA_point",100,-10,10,100,-10,10);
0557     DCA_point -> SetStats(0);
0558     DCA_point -> GetXaxis() -> SetTitle("X pos [mm]");
0559     DCA_point -> GetYaxis() -> SetTitle("Y pos [mm]");
0560     DCA_point -> GetXaxis() -> SetNdivisions(505);
0561
0562     DCA_distance = new TH1F("","DCA_distance",100,-10,10);
0563     DCA_distance -> SetStats(0);
0564     DCA_distance -> GetXaxis() -> SetTitle("DCA [mm]");
0565     DCA_distance -> GetYaxis() -> SetTitle("Entry");
0566     DCA_distance -> GetXaxis() -> SetNdivisions(505);
0567
0568     N_cluster_correlation = new TH2F("","N_cluster_correlation",100,0,4000,100,0,4000);
0569     N_cluster_correlation -> SetStats(0);
0570     N_cluster_correlation -> GetXaxis() -> SetTitle("inner N_cluster");
0571     N_cluster_correlation -> GetYaxis() -> SetTitle("Outer N_cluster");
0572     N_cluster_correlation -> GetXaxis() -> SetNdivisions(505);
0573
0574     N_cluster_correlation_close = new TH2F("","N_cluster_correlation_close",100,0,500,100,0,500);
0575     N_cluster_correlation_close -> SetStats(0);
0576     N_cluster_correlation_close -> GetXaxis() -> SetTitle("inner N_cluster");
0577     N_cluster_correlation_close -> GetYaxis() -> SetTitle("Outer N_cluster");
0578     N_cluster_correlation_close -> GetXaxis() -> SetNdivisions(505);
0579
0580     DCA_distance_inner_X = new TH2F("","DCA_distance_inner_X",100,-100,100,100,-10,10);
0581     DCA_distance_inner_X -> SetStats(0);
0582     DCA_distance_inner_X -> GetXaxis() -> SetTitle("Inner cluster X [mm]");
0583     DCA_distance_inner_X -> GetYaxis() -> SetTitle("DCA [mm]");
0584     DCA_distance_inner_X -> GetXaxis() -> SetNdivisions(505);
0585
0586     DCA_distance_inner_Y = new TH2F("","DCA_distance_inner_Y",100,-100,100,100,-10,10);
0587     DCA_distance_inner_Y -> SetStats(0);
0588     DCA_distance_inner_Y -> GetXaxis() -> SetTitle("Inner cluster Y [mm]");
0589     DCA_distance_inner_Y -> GetYaxis() -> SetTitle("DCA [mm]");
0590     DCA_distance_inner_Y -> GetXaxis() -> SetNdivisions(505);
0591
0592     DCA_distance_outer_X = new TH2F("","DCA_distance_outer_X",100,-130,130,100,-10,10);
0593     DCA_distance_outer_X -> SetStats(0);
0594     DCA_distance_outer_X -> GetXaxis() -> SetTitle("Outer cluster X [mm]");
0595     DCA_distance_outer_X -> GetYaxis() -> SetTitle("DCA [mm]");
0596     DCA_distance_outer_X -> GetXaxis() -> SetNdivisions(505);
0597
0598     DCA_distance_outer_Y = new TH2F("","DCA_distance_outer_Y",100,-130,130,100,-10,10);
0599     DCA_distance_outer_Y -> SetStats(0);
0600     DCA_distance_outer_Y -> GetXaxis() -> SetTitle("Outer cluster Y [mm]");
0601     DCA_distance_outer_Y -> GetYaxis() -> SetTitle("DCA [mm]");
0602     DCA_distance_outer_Y -> GetXaxis() -> SetNdivisions(505);
0603
0604     DCA_distance_inner_phi = new TH2F("","DCA_distance_inner_phi;Inner cluster #Phi [radian];DCA [cm]",100,-M_PI,M_PI,100,-1,1);
0605     DCA_distance_inner_phi -> SetStats(0);
0606     DCA_distance_inner_phi -> GetXaxis() -> SetNdivisions(505);
0607
0608     DCA_distance_outer_phi = new TH2F("","DCA_distance_outer_phi;Outer cluster #Phi [radian];DCA [cm]",100,-M_PI,M_PI,100,-1,1);
0609     DCA_distance_outer_phi -> SetStats(0);
0610     DCA_distance_outer_phi -> GetXaxis() -> SetNdivisions(505);
0611
0612     DCA_distance_inner_phi_peak_final = new TH2F("","DCA_distance_inner_phi_peak_final;Inner cluster #Phi [radian];DCA [cm]",100,-M_PI,M_PI,100,-1,1);
0613     DCA_distance_inner_phi_peak_final -> SetStats(0);
0614     DCA_distance_inner_phi_peak_final -> GetXaxis() -> SetNdivisions(505);
0615
0616     DCA_distance_outer_phi_peak_final = new TH2F("","DCA_distance_outer_phi_peak_final;Outer cluster #Phi [radian];DCA [cm]",100,-M_PI,M_PI,100,-1,1);
0617     DCA_distance_outer_phi_peak_final -> SetStats(0);
0618     DCA_distance_outer_phi_peak_final -> GetXaxis() -> SetNdivisions(505);
0619
0620     angle_diff_inner_phi = new TH2F("","angle_diff_inner_phi;Inner cluster #Phi [radian];#Delta#Phi (Inner - Outer) [radian]",100,-M_PI,M_PI,100,-0.03,0.03);
0621     angle_diff_inner_phi -> SetStats(0);
0622     angle_diff_inner_phi -> GetXaxis() -> SetNdivisions(505);
0623
0624     angle_diff_outer_phi = new TH2F("","angle_diff_outer_phi;Outer cluster #Phi [radian];#Delta#Phi (Inner - Outer) [radian]",100,-M_PI,M_PI,100,-0.03,0.03);
0625     angle_diff_outer_phi -> SetStats(0);
0626     angle_diff_outer_phi -> GetXaxis() -> SetNdivisions(505);
0627
0628     angle_diff_inner_phi_peak_final = new TH2F("","angle_diff_inner_phi_peak_final;Inner cluster #Phi [radian];#Delta#Phi (Inner - Outer) [radian]",100,-M_PI,M_PI,100,-0.03,0.03);
0629     angle_diff_inner_phi_peak_final -> SetStats(0);
0630     angle_diff_inner_phi_peak_final -> GetXaxis() -> SetNdivisions(505);
0631
0632     angle_diff_outer_phi_peak_final = new TH2F("","angle_diff_outer_phi_peak_final;Outer cluster #Phi [radian];#Delta#Phi (Inner - Outer) [radian]",100,-M_PI,M_PI,100,-0.03,0.03);
0633     angle_diff_outer_phi_peak_final -> SetStats(0);
0634     angle_diff_outer_phi_peak_final -> GetXaxis() -> SetNdivisions(505);
0635
0636
0637
0638
0639     final_angle_diff_inner_phi_radian_coarseX = new TH2F(
0640         "",
0641         "final_angle_diff_inner_phi_radian_coarseX;Inner cluster #Phi [radian];#Delta#Phi (Inner - Outer) [radian]",
0642         ana_map_version::final_angle_diff_inner_phi_radian_coarseX_XNbins, ana_map_version::final_angle_diff_inner_phi_radian_coarseX_Xmin, ana_map_version::final_angle_diff_inner_phi_radian_coarseX_Xmax,
0643         ana_map_version::final_angle_diff_inner_phi_radian_coarseX_YNbins, ana_map_version::final_angle_diff_inner_phi_radian_coarseX_Ymin, ana_map_version::final_angle_diff_inner_phi_radian_coarseX_Ymax
0644     );
0645     final_angle_diff_inner_phi_radian_coarseX -> SetStats(0);
0646     final_angle_diff_inner_phi_radian_coarseX -> GetXaxis() -> SetNdivisions(505);
0647
0648     final_angle_diff_inner_phi_degree_coarseX = new TH2F(
0649         "",
0650         "final_angle_diff_inner_phi_degree_coarseX;Inner cluster #Phi [degree];#Delta#Phi (Inner - Outer) [degree]",
0651         final_angle_diff_inner_phi_radian_coarseX -> GetNbinsX(), 0, 360,
0652         final_angle_diff_inner_phi_radian_coarseX -> GetNbinsY(), final_angle_diff_inner_phi_radian_coarseX -> GetYaxis() -> GetXmin() * radian_correction, final_angle_diff_inner_phi_radian_coarseX -> GetYaxis() -> GetXmax() * radian_correction
0653     );
0654     final_angle_diff_inner_phi_degree_coarseX -> SetStats(0);
0655     final_angle_diff_inner_phi_degree_coarseX -> GetXaxis() -> SetNdivisions(505);
0656
0657     final_angle_diff_radian_before = new TH1F("","final_angle_diff_radian_before;#Delta#phi (Inner - Outer) [radian];Entry",200,-0.045,0.045);
0658     final_angle_diff_radian_before -> SetStats(0);
0659     final_angle_diff_radian_before -> GetXaxis() -> SetNdivisions(505);
0660
0661     final_angle_diff_radian_post = new TH1F(
0662         "",
0663         "final_angle_diff_radian_post;#Delta#phi (Inner - Outer) [radian];Entry",
0664         final_angle_diff_radian_before -> GetNbinsX(),
0665         final_angle_diff_radian_before -> GetXaxis() -> GetXmin(),
0666         final_angle_diff_radian_before -> GetXaxis() -> GetXmax()
0667     );
0668     final_angle_diff_radian_post -> SetStats(0);
0669     final_angle_diff_radian_post -> GetXaxis() -> SetNdivisions(505);
0670
0671
0672     final_angle_diff_radian_DCA_2D_before = new TH2F("","final_angle_diff_radian_DCA_2D_before;Inner - Outer [radian];DCA distance [cm]",100,-1.5/radian_correction,1.5/radian_correction,100,-0.35,0.35);
0673     final_angle_diff_radian_DCA_2D_before -> SetStats(0);
0674     final_angle_diff_radian_DCA_2D_before -> GetXaxis() -> SetNdivisions(505);
0675
0676     final_angle_diff_radian_DCA_2D_post = new TH2F(
0677         "",
0678         "final_angle_diff_radian_DCA_2D_post;Inner - Outer [radian];DCA distance [cm]",
0679         final_angle_diff_radian_DCA_2D_before -> GetNbinsX(),
0680         final_angle_diff_radian_DCA_2D_before -> GetXaxis() -> GetXmin(),
0681         final_angle_diff_radian_DCA_2D_before -> GetXaxis() -> GetXmax(),
0682         final_angle_diff_radian_DCA_2D_before -> GetNbinsY(),
0683         final_angle_diff_radian_DCA_2D_before -> GetYaxis() -> GetXmin(),
0684         final_angle_diff_radian_DCA_2D_before -> GetYaxis() -> GetXmax()
0685     );
0686     final_angle_diff_radian_DCA_2D_post -> GetXaxis() -> SetNdivisions(505);
0687
0688     final_DCA_cm_inner_phi_radian_coarseX = new TH2F(
0689         "",
0690         "final_DCA_cm_inner_phi_radian_coarseX;Inner cluster #Phi [radian];DCA [cm]",
0691         ana_map_version::final_angle_diff_inner_phi_radian_coarseX_XNbins, ana_map_version::final_angle_diff_inner_phi_radian_coarseX_Xmin, ana_map_version::final_angle_diff_inner_phi_radian_coarseX_Xmax,
0692         100, -1., 1. // note : unit cm
0693     );
0694     final_DCA_cm_inner_phi_radian_coarseX -> GetXaxis() -> SetNdivisions(505);
0695
0696     final_DCA_mm_inner_phi_degree_coarseX = new TH2F(
0697         "",
0698         "final_DCA_mm_inner_phi_degree_coarseX;Inner cluster #Phi [degree];DCA [mm]",
0699         final_DCA_cm_inner_phi_radian_coarseX -> GetNbinsX(), 0, 360,
0700         final_DCA_cm_inner_phi_radian_coarseX -> GetNbinsY(), final_DCA_cm_inner_phi_radian_coarseX -> GetYaxis() -> GetXmin() * 10, final_DCA_cm_inner_phi_radian_coarseX -> GetYaxis() -> GetXmax() * 10 // note : unit mm
0701     );
0702     final_DCA_mm_inner_phi_degree_coarseX -> GetXaxis() -> SetNdivisions(505);
0703
0704     final_DCA_cm_before = new TH1F(
0705         "",
0706         "final_DCA_cm_before;DCA [cm];Entry",
0707         final_DCA_cm_inner_phi_radian_coarseX -> GetNbinsY(),
0708         final_DCA_cm_inner_phi_radian_coarseX -> GetYaxis() -> GetXmin(),
0709         final_DCA_cm_inner_phi_radian_coarseX -> GetYaxis() -> GetXmax()
0710     );
0711     final_DCA_cm_before -> GetXaxis() -> SetNdivisions(505);
0712
0713     final_DCA_cm_post = new TH1F(
0714         "",
0715         "final_DCA_cm_post;DCA [cm];Entry",
0716         final_DCA_cm_inner_phi_radian_coarseX -> GetNbinsY(),
0717         final_DCA_cm_inner_phi_radian_coarseX -> GetYaxis() -> GetXmin(),
0718         final_DCA_cm_inner_phi_radian_coarseX -> GetYaxis() -> GetXmax()
0719     );
0720     final_DCA_cm_post -> GetXaxis() -> SetNdivisions(505);
0721
0722
0723 }
0724
0725 // note : this function only prepare the pairs for the vertex XY calculation, it's like a general vertex for the whole run
0726 void INTTXYvtx::ProcessEvt(
0727     int event_i,
0728     vector<clu_info> temp_sPH_inner_nocolumn_vec,
0729     vector<clu_info> temp_sPH_outer_nocolumn_vec,
0730     vector<vector<double>> temp_sPH_nocolumn_vec,
0731     vector<vector<double>> temp_sPH_nocolumn_rz_vec,
0732     int NvtxMC,
0733     double TrigZvtxMC,
0734     bool PhiCheckTag,
0735     Long64_t bco_full,
0736     int is_min_bias,
0737     int is_min_bias_wozdc
0738 )
0739 {
0740     if (event_i%10000 == 0) {cout<<"In INTTXYvtx class, running event : "<<event_i<<endl;}
0741
0742     total_NClus = temp_sPH_inner_nocolumn_vec.size() + temp_sPH_outer_nocolumn_vec.size();
0743
0744     // note : the Move these two in the beginning of the function, in order to avoid those event-reject cuts
0745     N_cluster_correlation -> Fill(temp_sPH_inner_nocolumn_vec.size(),temp_sPH_outer_nocolumn_vec.size());
0746     N_cluster_correlation_close -> Fill(temp_sPH_inner_nocolumn_vec.size(),temp_sPH_outer_nocolumn_vec.size());
0747
0748     // todo : for the data, the is_min_bias_wozdc is used
0749     if (run_type == "MC" && is_min_bias_wozdc == 0) {return; cout<<"In INTTXYvtx class, MC, not min bias event"<<endl;}
0750     if (run_type != "MC" && is_min_bias_wozdc == 0) {return; cout<<"In INTTXYvtx class, data, not min bias event"<<endl;}
0751
0752     if (total_NClus < zvtx_cal_require) {return; cout<<"return confirmation"<<endl;}
0753
0754     if (run_type == "MC" && NvtxMC != 1) { return; cout<<"In INTTXYvtx class, event : "<<event_i<<" Nvtx : "<<NvtxMC<<" Nvtx more than one "<<endl;}
0755     if (PhiCheckTag == false)            { return; cout<<"In INTTXYvtx class, event : "<<event_i<<" Nvtx : "<<NvtxMC<<" Not full phi has hits "<<endl;}
0756
0757     if (temp_sPH_inner_nocolumn_vec.size() < 10 || temp_sPH_outer_nocolumn_vec.size() < 10 || total_NClus > N_clu_cut || total_NClus < N_clu_cutl)
0758     {
0759         return;
0760         printf("In INTTXYvtx class, event : %i, low clu continue, NClus : %lu \n", event_i, total_NClus);
0761     }
0762
0763     for ( int inner_i = 0; inner_i < temp_sPH_inner_nocolumn_vec.size(); inner_i++ )
0764     {
0765         for ( int outer_i = 0; outer_i < temp_sPH_outer_nocolumn_vec.size(); outer_i++ )
0766         {
0767             // note : try to ease the analysis and also make it quick.
0768             if (fabs(temp_sPH_inner_nocolumn_vec[inner_i].phi - temp_sPH_outer_nocolumn_vec[outer_i].phi) < 7) // todo : the pre phi cut is here, can be optimized
0769             {
0770                 cluster_pair_vec.push_back({{temp_sPH_inner_nocolumn_vec[inner_i].x,temp_sPH_inner_nocolumn_vec[inner_i].y},{temp_sPH_outer_nocolumn_vec[outer_i].x,temp_sPH_outer_nocolumn_vec[outer_i].y}});
0771             }
0772
0773         }
0774     }
0775
0776     for ( int inner_i = 0; inner_i < temp_sPH_inner_nocolumn_vec.size(); inner_i++ )
0777     {
0778         inner_pos_xy -> Fill(temp_sPH_inner_nocolumn_vec[inner_i].x,temp_sPH_inner_nocolumn_vec[inner_i].y);
0779         inner_outer_pos_xy -> Fill(temp_sPH_inner_nocolumn_vec[inner_i].x,temp_sPH_inner_nocolumn_vec[inner_i].y);
0780     }
0781
0782     for ( int outer_i = 0; outer_i < temp_sPH_outer_nocolumn_vec.size(); outer_i++ )
0783     {
0784         outer_pos_xy -> Fill(temp_sPH_outer_nocolumn_vec[outer_i].x,temp_sPH_outer_nocolumn_vec[outer_i].y);
0785         inner_outer_pos_xy -> Fill(temp_sPH_outer_nocolumn_vec[outer_i].x,temp_sPH_outer_nocolumn_vec[outer_i].y);
0786     }
0787 }
0788
0789 unsigned long INTTXYvtx::GetVecNele()
0790 {
0791     return cluster_pair_vec.size();
0792 }
0793
0794 pair<double,double> INTTXYvtx::GetFinalVTXxy()
0795 {
0796     return {current_vtxX, current_vtxY};
0797 }
0798
0799 void INTTXYvtx::ClearEvt()
0800 {
0801     return;
0802 }
0803
0804 pair<vector<TH2F *>, vector<TH1F*>> INTTXYvtx::GetHistFinal()
0805 {
0806     return {
0807         {DCA_distance_inner_phi_peak_final, angle_diff_inner_phi_peak_final, DCA_distance_outer_phi_peak_final, angle_diff_outer_phi_peak_final, xy_hist, xy_hist_bkgrm},
0808         {angle_diff_new_bkg_remove_final}
0809     };
0810 }
0811
0812 void INTTXYvtx::PrintPlots()
0813 {
0814     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
0815     inner_outer_pos_xy -> Draw("colz0");
0816     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
0817     c1 -> Print(Form("%s/inner_outer_pos_xy.pdf",out_folder_directory.c_str()));
0818     c1 -> Clear();
0819
0820     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
0821     inner_pos_xy -> Draw("colz0");
0822     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
0823     c1 -> Print(Form("%s/inner_pos_xy.pdf",out_folder_directory.c_str()));
0824     c1 -> Clear();
0825
0826     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
0827     outer_pos_xy -> Draw("colz0");
0828     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
0829     c1 -> Print(Form("%s/outer_pos_xy.pdf",out_folder_directory.c_str()));
0830     c1 -> Clear();
0831
0832     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
0833     N_cluster_correlation -> Draw("colz0");
0834     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
0835     c1 -> Print(Form("%s/N_cluster_correlation.pdf",out_folder_directory.c_str()));
0836     c1 -> Clear();
0837
0838     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
0839     N_cluster_correlation_close -> Draw("colz0");
0840     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
0841     c1 -> Print(Form("%s/N_cluster_correlation_close.pdf",out_folder_directory.c_str()));
0842     c1 -> Clear();
0843 }
0844
0845 // note : the pair here would be {radius of the circle, possible correction angle}
0846 pair<double,double> INTTXYvtx::GetCircleAngle(double fit_range_l, double fit_range_r)
0847 {
0848     cout<<"Get Circle of correction and possible correction angle ---------------------------- ---------------------------- ----------------------------"<<endl;
0849
0850     cos_fit -> SetParameters(4, 1.49143e-02,  185, 0.3); // todo : we may have to apply more constaints on the fitting
0851     cos_fit -> SetParLimits(2,0,360); // note : the peak location has to be positive
0852     // cos_fit -> SetParLimits(0,0,1000); // note : the amplitude has to be positive
0853
0854     // note : here is the test with a gaus fitting to find the peak
0855     gaus_fit -> SetParameters(-4.5, cos_fit->GetParameter(2), 50, 0);
0856     gaus_fit -> SetParLimits(0,-100,0); // note : the gaus distribution points down
0857
0858     subMacroPlotWorking(0,fit_range_l, fit_range_r, 40);
0859
0860     cout<<"circle radius by gaussian : "<<abs(gaus_fit->GetParameter(0) + gaus_fit->GetParameter(3))<<" possible correction angle : "<<gaus_fit->GetParameter(1)<<endl;
0861     cout<<"circle radius by cosine : "<<abs(cos_fit->GetParameter(0) - cos_fit->GetParameter(3))<<" possible correction angle : "<<cos_fit->GetParameter(2)<<endl;
0862
0863     // return {abs(gaus_fit->GetParameter(0) + gaus_fit->GetParameter(3)), gaus_fit->GetParameter(1)}; // note : gaussian based
0864     return {abs(cos_fit->GetParameter(0) - cos_fit->GetParameter(3)), cos_fit->GetParameter(2)}; // note : cosine wave based
0865 }
0866
0867 vector<pair<double,double>> INTTXYvtx::MacroVTXSquare(double length, int N_trial, bool draw_plot_opt)
0868 {
0869     double original_length = length;
0870     pair<double,double> origin = {0,0};
0871     vector<pair<double,double>> vtx_vec = Get4vtx(origin,length); // vtx_vec.push_back(origin);
0872     int small_index;
0873     vector<double> small_info_vec;
0874     vector<double> grr_x; grr_x.clear();
0875     vector<double> grr_E; grr_E.clear();
0876     vector<double> grr_y; grr_y.clear();
0877
0878     vector<double> All_FitError_DCA_Y; All_FitError_DCA_Y.clear();
0879     vector<double> All_FitError_DCA_X; All_FitError_DCA_X.clear();
0880     vector<double> All_FitError_angle_Y; All_FitError_angle_Y.clear();
0881     vector<double> All_FitError_angle_X; All_FitError_angle_X.clear();
0882
0883     vector<double> Winner_FitError_DCA_Y; Winner_FitError_DCA_Y.clear();
0884     vector<double> Winner_FitError_DCA_X; Winner_FitError_DCA_X.clear();
0885     vector<double> Winner_FitError_angle_Y; Winner_FitError_angle_Y.clear();
0886     vector<double> Winner_FitError_angle_X; Winner_FitError_angle_X.clear();
0887
0888     cout<<"In INTTXYvtx::MacroVTXSquare, N pairs : "<<cluster_pair_vec.size()<<endl;
0889
0890     if (print_message_opt == true) {cout<<N_trial<<" runs, smart. which gives you the resolution down to "<<length/pow(2,N_trial)<<" mm"<<endl;}
0891
0892     for (int i = 0; i < N_trial; i++)
0893     {
0894         if (print_message_opt == true) {cout<<"~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~~~~"<<" step "<<i<<" ~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~~~~"<<endl;}
0895         for (int i1 = 0; i1 < vtx_vec.size(); i1++)
0896         {
0897             if (print_message_opt == true) {cout<<"tested vertex : "<<vtx_vec[i1].first<<" "<<vtx_vec[i1].second<<endl;}
0898             current_vtxX = vtx_vec[i1].first;
0899             current_vtxY = vtx_vec[i1].second;
0900             vector<double> info_vec = subMacroVTXxyCorrection(i,i1, draw_plot_opt);
0901             cout<<"trial : "<<i<<" vertex : "<<i1<<" DCA fit error : "<<info_vec[3]<<" angle diff fit error : "<<info_vec[5]<<endl;
0902             All_FitError_DCA_Y.push_back(info_vec[3]);
0903             All_FitError_DCA_X.push_back(i);
0904             All_FitError_angle_Y.push_back(info_vec[5]);
0905             All_FitError_angle_X.push_back(i);
0906
0907             if (i1 == 0)
0908             {
0909                 small_info_vec = info_vec;
0910                 small_index = i1;
0911
0912                 TH2F_FakeClone(DCA_distance_inner_phi_peak,DCA_distance_inner_phi_peak_final);
0913                 TH2F_FakeClone(angle_diff_inner_phi_peak,angle_diff_inner_phi_peak_final);
0914                 TH2F_FakeClone(DCA_distance_outer_phi_peak,DCA_distance_outer_phi_peak_final);
0915                 TH2F_FakeClone(angle_diff_outer_phi_peak,angle_diff_outer_phi_peak_final);
0916                 TH1F_FakeClone(angle_diff_new_bkg_remove, angle_diff_new_bkg_remove_final);
0917             }
0918             else
0919             {
0920                 if (info_vec[3] < small_info_vec[3] && info_vec[5] < small_info_vec[5]) // note : the fit error of the pol0 fit
0921                 {
0922                     small_info_vec = info_vec;
0923                     small_index = i1;
0924
0925                     TH2F_FakeClone(DCA_distance_inner_phi_peak,DCA_distance_inner_phi_peak_final);
0926                     TH2F_FakeClone(angle_diff_inner_phi_peak,angle_diff_inner_phi_peak_final);
0927                     TH2F_FakeClone(DCA_distance_outer_phi_peak,DCA_distance_outer_phi_peak_final);
0928                     TH2F_FakeClone(angle_diff_outer_phi_peak,angle_diff_outer_phi_peak_final);
0929                     TH1F_FakeClone(angle_diff_new_bkg_remove, angle_diff_new_bkg_remove_final);
0930                 }
0931             }
0932             if (print_message_opt == true){cout<<" "<<endl;}
0933
0934             ClearHist(1);
0935         }
0936
0937         if (print_message_opt == true) {cout<<"the Quadrant "<<small_index<<" won the competition"<<endl;}
0938
0939         Winner_FitError_DCA_Y.push_back(small_info_vec[3]);
0940         Winner_FitError_DCA_X.push_back(i);
0941         Winner_FitError_angle_Y.push_back(small_info_vec[5]);
0942         Winner_FitError_angle_X.push_back(i);
0943
0944         grr_x.push_back(i);
0945         grr_y.push_back(small_index);
0946         grr_E.push_back(0);
0947
0948         // note : generating the new 4 vertex for the next comparison
0949         // note : start to shrink the square
0950         if (i != N_trial - 1)
0951         {
0952             origin = {(vtx_vec[small_index].first + origin.first)/2., (vtx_vec[small_index].second + origin.second)/2.};
0953             // cout<<"test : "<<origin.first<<" "<<origin.second<<" length: "<<length<<endl;
0954             // if (small_index == 4) {length /= 1.5;}
0955             // else {length /= 2.;}
0956             length /= 2.;
0957             vtx_vec = Get4vtx(origin,length); // vtx_vec.push_back(origin);
0958         }
0959     }
0960
0961     if (draw_plot_opt == true) {DrawTGraphErrors(grr_x, grr_y, grr_E, grr_E, out_folder_directory, {Form("Square_scan_history_%.1fmm_%iTrials", original_length, N_trial), "nth scan", "Winner index", "APL"});}
0962     if (draw_plot_opt == true) {Draw2TGraph(All_FitError_angle_X, All_FitError_angle_Y, Winner_FitError_angle_X, Winner_FitError_angle_Y, out_folder_directory, {Form("Angle_diff_fit_error_%iTrials", N_trial), "n iteration", "#Delta#phi fit error [radian]"});}
0963     if (draw_plot_opt == true) {Draw2TGraph(All_FitError_DCA_X, All_FitError_DCA_Y, Winner_FitError_DCA_X, Winner_FitError_DCA_Y, out_folder_directory, {Form("DCA_fit_error_%iTrials", N_trial), "n iteration", "DCA fit error [cm]"});}
0964
0965     return {
0966         vtx_vec[small_index], // note : the best vertex
0967         origin,               // note : the origin in that trial
0968         {small_info_vec[3],small_info_vec[5]},   // note : horizontal_fit_inner -> GetParError(0),  horizontal_fit_angle_diff_inner -> GetParError(0)
0969         {small_info_vec[10],small_info_vec[11]}, // note : horizontal_fit_inner -> GetParameter(0), horizontal_fit_angle_diff_inner -> GetParameter(0)
0970         {small_info_vec[7],small_info_vec[9]},   // note : horizontal_fit_outer -> GetParError(0),  horizontal_fit_angle_diff_outer -> GetParError(0)
0971         {small_info_vec[12],small_info_vec[13]}, // note : horizontal_fit_outer -> GetParameter(0), horizontal_fit_angle_diff_outer -> GetParameter(0)
0972         {small_info_vec[14],small_info_vec[15]},  // note : the mean and stddev of angle_diff
0973         {small_info_vec[16],small_info_vec[17]},  // note : the mean and stddev of DCA_distance
0974         {small_info_vec[0],small_info_vec[1]},    // note : the mean and stddev of angle_diff, but with the background removed
0975     };
0976 }
0977
0978 vector<double> INTTXYvtx::subMacroVTXxyCorrection(int test_index, int trial_index, bool draw_plot_opt)
0979 {
0980     int true_trial_index = test_index * 4 + trial_index;
0981     string sub_out_folder_name;
0982     if (draw_plot_opt == true){
0983         sub_out_folder_name = Form("%s/New_trial_square_%i_%i",out_folder_directory.c_str(), test_index, trial_index);
0984         if (filesystem::exists(sub_out_folder_name.c_str()) == false) {system(Form("mkdir %s",sub_out_folder_name.c_str()));}
0985     }
0986
0987     vector<double> out_vec = GetVTXxyCorrection_new(true_trial_index);
0988     if (draw_plot_opt == true){PrintPlotsVTXxy(sub_out_folder_name, 1);}
0989     // ClearHist(1);
0990
0991     return out_vec;
0992 }
0993
0994 // note : {circle radius, possible correction angle, the chi2/NDF of pol0 fit}
0995 vector<double> INTTXYvtx::GetVTXxyCorrection_new(int trial_index)
0996 {
0997     if (print_message_opt == true) {
0998         cout<<"Trial : "<<trial_index<<"---------------------------- ---------------------------- ----------------------------"<<endl;
0999         cout<<"Given vertex: "<<current_vtxX <<" "<<current_vtxY<<endl;
1000     }
1001
1002
1003     cos_fit -> SetParameters(4, 1.49143e-02,  185, 0.3); // todo : we may have to apply more constaints on the fitting
1004     // cos_fit -> SetParLimits(0,0,1000); // note : the amplitude has to be positive
1005     cos_fit -> SetParLimits(2,0,360); // note : the peak location has to be positive
1006
1007     // note : here is the test with a gaus fitting to find the peak
1008     gaus_fit -> SetParameters(-4.5, 197, 50, 0);
1009     gaus_fit -> SetParLimits(0,-100,0); // note : the gaus distribution points down
1010     // DCA_distance_inner_phi_peak_profile -> Fit(gaus_fit, "N","",100, 260);
1011     // cout<<"test, gaus fit range : "<<gaus_fit->GetParameter(1) - 25<<" "<<gaus_fit->GetParameter(1) + 25<<endl;
1012
1013     subMacroPlotWorking(1,100,260,25);
1014
1015     return {
1016         angle_diff_new_bkg_remove -> GetMean(), angle_diff_new_bkg_remove -> GetStdDev(), // note : angle diff stddev and error (1D histogram)
1017         horizontal_fit_inner -> GetChisquare() / double(horizontal_fit_inner -> GetNDF()), horizontal_fit_inner->GetParError(0),  // note : inner DCA, pol0
1018         horizontal_fit_angle_diff_inner -> GetChisquare() / double(horizontal_fit_angle_diff_inner -> GetNDF()), horizontal_fit_angle_diff_inner->GetParError(0), // note : inner angle diff, pol0
1019         horizontal_fit_outer -> GetChisquare() / double(horizontal_fit_outer -> GetNDF()), horizontal_fit_outer->GetParError(0), // note : outer DCA, pol0
1020         horizontal_fit_angle_diff_outer -> GetChisquare() / double(horizontal_fit_angle_diff_outer -> GetNDF()), horizontal_fit_angle_diff_outer->GetParError(0),  // note : outer angle diff, pol0
1021         horizontal_fit_inner -> GetParameter(0), horizontal_fit_angle_diff_inner -> GetParameter(0), // note : 10, 11
1022         horizontal_fit_outer -> GetParameter(0), horizontal_fit_angle_diff_outer -> GetParameter(0), // note : 12, 13
1023         angle_diff -> GetMean(), angle_diff -> GetStdDev(), // note : 14, 15
1024         DCA_distance -> GetMean(), DCA_distance -> GetStdDev() // note : 16, 17
1025     };
1026 }
1027
1028 void INTTXYvtx::subMacroPlotWorking(bool phi_correction, double cos_fit_rangel, double cos_fit_ranger, double guas_fit_range)
1029 {
1030
1031     for (int i = 0; i < cluster_pair_vec.size(); i++)
1032     {
1033         vector<double> DCA_info_vec = calculateDistanceAndClosestPoint(
1034             cluster_pair_vec[i].first.x, cluster_pair_vec[i].first.y,
1035             cluster_pair_vec[i].second.x, cluster_pair_vec[i].second.y,
1036             current_vtxX, current_vtxY
1037         );
1038
1039         double DCA_sign = calculateAngleBetweenVectors(
1040             cluster_pair_vec[i].second.x, cluster_pair_vec[i].second.y,
1041             cluster_pair_vec[i].first.x, cluster_pair_vec[i].first.y,
1042             current_vtxX, current_vtxY
1043         );
1044
1045         if (phi_correction == true)
1046         {
1047             // cout<<"option selected "<<endl;
1048             Clus_InnerPhi_Offset = (cluster_pair_vec[i].first.y - current_vtxY < 0) ? atan2(cluster_pair_vec[i].first.y - current_vtxY, cluster_pair_vec[i].first.x - current_vtxX) * (180./TMath::Pi()) + 360 : atan2(cluster_pair_vec[i].first.y - current_vtxY, cluster_pair_vec[i].first.x - current_vtxX) * (180./TMath::Pi());
1049             Clus_OuterPhi_Offset = (cluster_pair_vec[i].second.y - current_vtxY < 0) ? atan2(cluster_pair_vec[i].second.y - current_vtxY, cluster_pair_vec[i].second.x - current_vtxX) * (180./TMath::Pi()) + 360 : atan2(cluster_pair_vec[i].second.y - current_vtxY, cluster_pair_vec[i].second.x - current_vtxX) * (180./TMath::Pi());
1050
1051             // note : to have the radian
1052             Clus_InnerPhi_Offset_radian = atan2(cluster_pair_vec[i].first.y  - current_vtxY, cluster_pair_vec[i].first.x  - current_vtxX);
1053             Clus_OuterPhi_Offset_radian = atan2(cluster_pair_vec[i].second.y - current_vtxY, cluster_pair_vec[i].second.x - current_vtxX);
1054
1055         }
1056         else // note : phi_correction == false
1057         {
1058             Clus_InnerPhi_Offset = (cluster_pair_vec[i].first.y < 0) ? atan2(cluster_pair_vec[i].first.y, cluster_pair_vec[i].first.x) * (180./TMath::Pi()) + 360 : atan2(cluster_pair_vec[i].first.y, cluster_pair_vec[i].first.x) * (180./TMath::Pi());
1059             Clus_OuterPhi_Offset = (cluster_pair_vec[i].second.y < 0) ? atan2(cluster_pair_vec[i].second.y, cluster_pair_vec[i].second.x) * (180./TMath::Pi()) + 360 : atan2(cluster_pair_vec[i].second.y, cluster_pair_vec[i].second.x) * (180./TMath::Pi());
1060
1061             // note : to have the radian
1062             Clus_InnerPhi_Offset_radian = atan2(cluster_pair_vec[i].first.y,  cluster_pair_vec[i].first.x);
1063             Clus_OuterPhi_Offset_radian = atan2(cluster_pair_vec[i].second.y, cluster_pair_vec[i].second.x);
1064
1065         }
1066
1067         // double phi_test = (cluster_pair_vec[i].first.y < 0) ? atan2(cluster_pair_vec[i].first.y, cluster_pair_vec[i].first.x) * (180./TMath::Pi()) + 360 : atan2(cluster_pair_vec[i].first.y, cluster_pair_vec[i].first.x) * (180./TMath::Pi());
1068         // double phi_test_offset = (cluster_pair_vec[i].first.y - current_vtxY < 0) ? atan2(cluster_pair_vec[i].first.y - current_vtxY, cluster_pair_vec[i].first.x - current_vtxX) * (180./TMath::Pi()) + 360 : atan2(cluster_pair_vec[i].first.y - current_vtxY, cluster_pair_vec[i].first.x - current_vtxX) * (180./TMath::Pi());
1069         // cout<<"angle offset test : "<<phi_test<<", "<<phi_test_offset<<endl;
1070
1071         angle_correlation -> Fill(Clus_InnerPhi_Offset, Clus_OuterPhi_Offset);
1072         angle_diff_DCA_dist -> Fill(Clus_InnerPhi_Offset - Clus_OuterPhi_Offset, DCA_sign);
1073         angle_diff -> Fill(abs(Clus_InnerPhi_Offset - Clus_OuterPhi_Offset));
1074         DCA_point -> Fill(DCA_info_vec[1], DCA_info_vec[2]);
1075         DCA_distance -> Fill(DCA_sign);
1076         DCA_distance_inner_X -> Fill(cluster_pair_vec[i].first.x, DCA_sign);
1077         DCA_distance_inner_Y -> Fill(cluster_pair_vec[i].first.y, DCA_sign);
1078         DCA_distance_outer_X -> Fill(cluster_pair_vec[i].second.x, DCA_sign);
1079         DCA_distance_outer_Y -> Fill(cluster_pair_vec[i].second.y, DCA_sign);
1080
1081         angle_diff_new -> Fill(abs(Clus_InnerPhi_Offset - Clus_OuterPhi_Offset));
1082         angle_diff_radian -> Fill(Clus_InnerPhi_Offset_radian - Clus_OuterPhi_Offset_radian);
1083
1084         // note : for the special treatment
1085         DCA_distance_inner_phi -> Fill(Clus_InnerPhi_Offset_radian, DCA_sign / 10.);
1086         DCA_distance_outer_phi -> Fill(Clus_OuterPhi_Offset_radian, DCA_sign / 10.);
1087         angle_diff_inner_phi   -> Fill(Clus_InnerPhi_Offset_radian, Clus_InnerPhi_Offset_radian - Clus_OuterPhi_Offset_radian);
1088         angle_diff_outer_phi   -> Fill(Clus_OuterPhi_Offset_radian, Clus_InnerPhi_Offset_radian - Clus_OuterPhi_Offset_radian);
1089
1090     } // note : end of the loop for the cluster pair
1091
1092     // note : ----------- ----------- ----------- ---------- ----------- ----------- ---------- ----------- ----------- ----------- -----------
1093     DCA_distance_inner_phi_peak = (TH2F*)DCA_distance_inner_phi -> Clone();
1094     TH2F_threshold_advanced_2(DCA_distance_inner_phi_peak, 0.5); // todo : the background cut can be modified, the ratio 0.5
1095     DCA_distance_inner_phi_peak_profile =  DCA_distance_inner_phi_peak->ProfileX("DCA_distance_inner_phi_peak_profile");
1096     // TGraph * DCA_distance_inner_phi_peak_profile_graph = new TGraph();
1097     double point_index = 0;
1098     vector<double> hist_column_content = SumTH2FColumnContent(DCA_distance_inner_phi_peak);
1099     for (int i = 0; i < DCA_distance_inner_phi_peak_profile->GetNbinsX(); i++){
1100         if (hist_column_content[i] < 5){continue;} // note : in order to remove some remaining background
1101
1102         DCA_distance_inner_phi_peak_profile_graph -> SetPoint(point_index, DCA_distance_inner_phi_peak_profile->GetBinCenter(i+1), DCA_distance_inner_phi_peak_profile->GetBinContent(i+1));
1103         // cout<<"("<<DCA_distance_inner_phi_peak_profile->GetBinCenter(i+1)<<", "<< DCA_distance_inner_phi_peak_profile->GetBinContent(i+1)<<")"<<endl;
1104         point_index += 1;
1105     }
1106
1107     // cos_fit -> SetParameters(4, 1.49143e-02,  185, 0.3); // todo : we may have to apply more constaints on the fitting
1108     // cos_fit -> SetParLimits(0,0,1000); // note : the amplitude has to be positive
1109     // cos_fit -> SetParLimits(2,0,360); // note : the peak location has to be positive
1110
1111     // note : here is the test with a gaus fitting to find the peak
1112     // gaus_fit -> SetParameters(-4.5, cos_fit->GetParameter(2), 50, 0);
1113     // gaus_fit -> SetParLimits(0,-100,0); // note : the gaus distribution points down
1114     // DCA_distance_inner_phi_peak_profile -> Fit(gaus_fit, "N","",100, 260);
1115     // cout<<"test, gaus fit range : "<<gaus_fit->GetParameter(1) - 25<<" "<<gaus_fit->GetParameter(1) + 25<<endl;
1116
1117     horizontal_fit_inner -> SetParameter(0,0);
1118
1119     // todo : the fit range of the gaussian fit can be modified here
1120     DCA_distance_inner_phi_peak_profile_graph -> Fit(horizontal_fit_inner,"NQ","",-M_PI, M_PI);
1121     // DCA_distance_inner_phi_peak_profile_graph -> Fit(gaus_fit, "NQ","",cos_fit->GetParameter(2) - guas_fit_range, cos_fit->GetParameter(2) + guas_fit_range); // note : what we want and need is the peak position, so we fit the peak again
1122     // DCA_distance_inner_phi_peak_profile_graph -> Fit(cos_fit, "NQ","",cos_fit_rangel, cos_fit_ranger);
1123
1124     // note : ----------- ----------- ----------- ---------- ----------- ----------- ---------- ----------- ----------- ----------- -----------
1125     angle_diff_new_bkg_remove = (TH1F*)angle_diff_new -> Clone();
1126     angle_diff_new_bkg_remove -> SetLineColor(2);
1127     Hist_1D_bkg_remove(angle_diff_new_bkg_remove, 1.5);
1128
1129     // note : ----------- ----------- ----------- ---------- ----------- ----------- ---------- ----------- ----------- ----------- -----------
1130     angle_diff_inner_phi_peak = (TH2F*)angle_diff_inner_phi -> Clone();
1131     // TH2F_threshold_advanced(angle_diff_inner_phi_peak, 0.5);
1132     TH2F_threshold_advanced_2(angle_diff_inner_phi_peak, 0.5); // todo : threshold ratio can be modified here
1133     hist_column_content = SumTH2FColumnContent(angle_diff_inner_phi_peak);
1134     angle_diff_inner_phi_peak_profile =  angle_diff_inner_phi_peak->ProfileX("angle_diff_inner_phi_peak_profile");
1135     // TGraph * angle_diff_inner_phi_peak_profile_graph = new TGraph();
1136     point_index = 0;
1137     for (int i = 0; i < angle_diff_inner_phi_peak_profile->GetNbinsX(); i++){
1138         if (hist_column_content[i] < 5){continue;} // note : in order to remove some remaining background
1139
1140         angle_diff_inner_phi_peak_profile_graph -> SetPoint(point_index, angle_diff_inner_phi_peak_profile->GetBinCenter(i+1), angle_diff_inner_phi_peak_profile->GetBinContent(i+1));
1141         // cout<<"("<<angle_diff_inner_phi_peak_profile->GetBinCenter(i+1)<<", "<< angle_diff_inner_phi_peak_profile->GetBinContent(i+1)<<")"<<endl;
1142         point_index += 1;
1143     }
1144
1145     horizontal_fit_angle_diff_inner -> SetParameter(0,0);
1146     angle_diff_inner_phi_peak_profile_graph -> Fit(horizontal_fit_angle_diff_inner,"NQ","",-M_PI, M_PI);
1147
1148     // note : ----------- ----------- ----------- ---------- ----------- ----------- ---------- ----------- ----------- ----------- -----------
1149     DCA_distance_outer_phi_peak = (TH2F*)DCA_distance_outer_phi -> Clone();
1150     TH2F_threshold_advanced_2(DCA_distance_outer_phi_peak, 0.5); // todo : the background cut can be modified, the ratio 0.5
1151     DCA_distance_outer_phi_peak_profile =  DCA_distance_outer_phi_peak->ProfileX("DCA_distance_outer_phi_peak_profile");
1152     // TGraph * DCA_distance_outer_phi_peak_profile_graph = new TGraph();
1153     point_index = 0;
1154     hist_column_content = SumTH2FColumnContent(DCA_distance_outer_phi_peak);
1155     for (int i = 0; i < DCA_distance_outer_phi_peak_profile->GetNbinsX(); i++){
1156         if (hist_column_content[i] < 5){continue;} // note : in order to remove some remaining background
1157
1158         DCA_distance_outer_phi_peak_profile_graph -> SetPoint(point_index, DCA_distance_outer_phi_peak_profile->GetBinCenter(i+1), DCA_distance_outer_phi_peak_profile->GetBinContent(i+1));
1159         // cout<<"("<<DCA_distance_outer_phi_peak_profile->GetBinCenter(i+1)<<", "<< DCA_distance_outer_phi_peak_profile->GetBinContent(i+1)<<")"<<endl;
1160         point_index += 1;
1161     }
1162
1163     horizontal_fit_outer -> SetParameter(0,0);
1164     // todo : the fit range of the gaussian fit can be modified here
1165     DCA_distance_outer_phi_peak_profile_graph -> Fit(horizontal_fit_outer,"NQ","",-M_PI, M_PI);
1166
1167     // note : ----------- ----------- ----------- ---------- ----------- ----------- ---------- ----------- ----------- ----------- -----------
1168     angle_diff_outer_phi_peak = (TH2F*)angle_diff_outer_phi -> Clone();
1169     // TH2F_threshold_advanced(angle_diff_outer_phi_peak, 0.5);
1170     TH2F_threshold_advanced_2(angle_diff_outer_phi_peak, 0.5); // todo : threshold ratio can be modified here
1171     hist_column_content = SumTH2FColumnContent(angle_diff_outer_phi_peak);
1172     angle_diff_outer_phi_peak_profile =  angle_diff_outer_phi_peak->ProfileX("angle_diff_outer_phi_peak_profile");
1173     // TGraph * angle_diff_outer_phi_peak_profile_graph = new TGraph();
1174     point_index = 0;
1175     for (int i = 0; i < angle_diff_outer_phi_peak_profile->GetNbinsX(); i++){
1176         if (hist_column_content[i] < 5){continue;} // note : in order to remove some remaining background
1177
1178         angle_diff_outer_phi_peak_profile_graph -> SetPoint(point_index, angle_diff_outer_phi_peak_profile->GetBinCenter(i+1), angle_diff_outer_phi_peak_profile->GetBinContent(i+1));
1179         // cout<<"("<<angle_diff_outer_phi_peak_profile->GetBinCenter(i+1)<<", "<< angle_diff_outer_phi_peak_profile->GetBinContent(i+1)<<")"<<endl;
1180         point_index += 1;
1181     }
1182
1183     horizontal_fit_angle_diff_outer -> SetParameter(0,0);
1184     angle_diff_outer_phi_peak_profile_graph -> Fit(horizontal_fit_angle_diff_outer,"NQ","",-M_PI, M_PI);
1185     // note : ----------- ----------- ----------- ---------- ----------- ----------- ---------- ----------- ----------- ----------- -----------
1186
1187
1188
1189     angle_diff_inner_phi_peak_profile_graph -> Set(0);
1190     angle_diff_outer_phi_peak_profile_graph -> Set(0);
1191     DCA_distance_inner_phi_peak_profile_graph -> Set(0);
1192     DCA_distance_outer_phi_peak_profile_graph -> Set(0);
1193
1194     if (print_message_opt == true) {cout<<"circle radius : "<<abs(gaus_fit->GetParameter(0) + gaus_fit->GetParameter(3))<<" possible correction angle : "<<gaus_fit->GetParameter(1)<<endl;}
1195     // return {abs(gaus_fit->GetParameter(0) + gaus_fit->GetParameter(3)), gaus_fit->GetParameter(1)}; // note : gaussian based
1196
1197     // return {
1198     //     abs(cos_fit->GetParameter(0) - cos_fit->GetParameter(3)), cos_fit->GetParameter(2), // note : cosine wave based
1199     //     angle_diff_new_bkg_remove -> GetStdDev(), angle_diff_new_bkg_remove -> GetStdDevError() // note : the std dev of the angle diff 1D distribution
1200     // };
1201 }
1202
1203
1204
1205 void INTTXYvtx::PrintPlotsVTXxy(string sub_out_folder_name, int print_option)
1206 {
1207     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1208     angle_correlation -> Draw("colz0");
1209     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1210     c1 -> Print(Form("%s/angle_correlation.pdf", sub_out_folder_name.c_str()));
1211     c1 -> Clear();
1212
1213     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1214     angle_diff_DCA_dist -> Draw("colz0");
1215     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1216     c1 -> Print(Form("%s/angle_diff_DCA_dist.pdf", sub_out_folder_name.c_str()));
1217     c1 -> Clear();
1218
1219     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1220     angle_diff -> SetMinimum(0);
1221     angle_diff -> Draw("hist");
1222     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1223     c1 -> Print(Form("%s/angle_diff.pdf", sub_out_folder_name.c_str()));
1224     c1 -> Clear();
1225
1226     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1227     angle_diff_inner_phi -> Draw("colz0");
1228     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1229     c1 -> Print(Form("%s/angle_diff_inner_phi.pdf", sub_out_folder_name.c_str()));
1230     c1 -> Clear();
1231
1232     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1233     angle_diff_outer_phi -> Draw("colz0");
1234     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1235     c1 -> Print(Form("%s/angle_diff_outer_phi.pdf", sub_out_folder_name.c_str()));
1236     c1 -> Clear();
1237
1238     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1239     DCA_point -> Draw("colz0");
1240     // ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s, peak : %f", plot_text.c_str(), peek));
1241     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1242     c1 -> Print(Form("%s/DCA_point.pdf", sub_out_folder_name.c_str()));
1243     c1 -> Clear();
1244
1245     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1246     DCA_distance_inner_phi -> Draw("colz0");
1247     // ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s, peak : %f", plot_text.c_str(), peek));
1248     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1249     c1 -> Print(Form("%s/DCA_distance_inner_phi.pdf", sub_out_folder_name.c_str()));
1250     c1 -> Clear();
1251
1252     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1253     DCA_distance_outer_phi -> Draw("colz0");
1254     // ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s, peak : %f", plot_text.c_str(), peek));
1255     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1256     c1 -> Print(Form("%s/DCA_distance_outer_phi.pdf", sub_out_folder_name.c_str()));
1257     c1 -> Clear();
1258
1259     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1260     DCA_distance -> Draw("hist");
1261     // ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s, peak : %f", plot_text.c_str(), peek));
1262     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1263     c1 -> Print(Form("%s/DCA_distance.pdf", sub_out_folder_name.c_str()));
1264     c1 -> Clear();
1265
1266     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1267     DCA_distance_inner_phi_peak -> SetStats(0);
1268     DCA_distance_inner_phi_peak -> GetXaxis() -> SetTitle( DCA_distance_inner_phi -> GetXaxis() -> GetTitle() );
1269     DCA_distance_inner_phi_peak -> GetYaxis() -> SetTitle( DCA_distance_inner_phi -> GetYaxis() -> GetTitle() );
1270     DCA_distance_inner_phi_peak -> Draw("colz0");
1271     DCA_distance_inner_phi_peak_profile -> Draw("same");
1272     // cos_fit -> Draw("l same");
1273     // gaus_fit -> Draw("l same");
1274     horizontal_fit_inner -> Draw("l same");
1275     // ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s, peak : %f", plot_text.c_str(), peek));
1276     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1277     draw_text -> DrawLatex(0.25, 0.84, Form("#color[2]{Assumed vertex: %.4f cm, %.4f cm}", current_vtxX/10., current_vtxY/10.));
1278     // draw_text -> DrawLatex(0.25, 0.80, Form("#color[2]{Fit: -%.2f cos(%.2f(x - %.2f)) + %.2f}", cos_fit -> GetParameter(0), cos_fit -> GetParameter(1), cos_fit -> GetParameter(2), cos_fit -> GetParameter(3)));
1279     draw_text -> DrawLatex(0.25, 0.80, Form("#color[2]{Pol0 fit chi2/NDF: %.3f, fit error: %.4f}", horizontal_fit_inner -> GetChisquare() / double(horizontal_fit_inner -> GetNDF()), horizontal_fit_inner->GetParError(0)));
1280     c1 -> Print(Form("%s/DCA_distance_inner_phi_peak.pdf", sub_out_folder_name.c_str()));
1281     c1 -> Clear();
1282
1283     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1284
1285     DCA_distance_outer_phi_peak -> SetStats(0);
1286     DCA_distance_outer_phi_peak -> GetXaxis() -> SetTitle( DCA_distance_outer_phi -> GetXaxis() -> GetTitle() );
1287     DCA_distance_outer_phi_peak -> GetYaxis() -> SetTitle( DCA_distance_outer_phi -> GetYaxis() -> GetTitle() );
1288     DCA_distance_outer_phi_peak -> Draw("colz0");
1289     DCA_distance_outer_phi_peak_profile -> Draw("same");
1290     horizontal_fit_outer -> Draw("l same");
1291     // ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s, peak : %f", plot_text.c_str(), peek));
1292     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1293     draw_text -> DrawLatex(0.25, 0.84, Form("#color[2]{Assumed vertex: %.4f cm, %.4f cm}", current_vtxX/10., current_vtxY/10.));
1294     draw_text -> DrawLatex(0.25, 0.80, Form("#color[2]{Pol0 fit chi2/NDF: %.3f, fit error: %.4f}", horizontal_fit_outer -> GetChisquare() / double(horizontal_fit_outer -> GetNDF()), horizontal_fit_outer->GetParError(0)));
1295     c1 -> Print(Form("%s/DCA_distance_outer_phi_peak.pdf", sub_out_folder_name.c_str()));
1296     c1 -> Clear();
1297
1298
1299     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1300     DCA_distance_inner_X -> Draw("colz0");
1301     // ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s, peak : %f", plot_text.c_str(), peek));
1302     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1303     c1 -> Print(Form("%s/DCA_distance_inner_X.pdf", sub_out_folder_name.c_str()));
1304     c1 -> Clear();
1305
1306     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1307     DCA_distance_inner_Y -> Draw("colz0");
1308     // ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s, peak : %f", plot_text.c_str(), peek));
1309     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1310     c1 -> Print(Form("%s/DCA_distance_inner_Y.pdf", sub_out_folder_name.c_str()));
1311     c1 -> Clear();
1312
1313     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1314     DCA_distance_outer_X -> Draw("colz0");
1315     // ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s, peak : %f", plot_text.c_str(), peek));
1316     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1317     c1 -> Print(Form("%s/DCA_distance_outer_X.pdf", sub_out_folder_name.c_str()));
1318     c1 -> Clear();
1319
1320     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1321     DCA_distance_outer_Y -> Draw("colz0");
1322     // ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s, peak : %f", plot_text.c_str(), peek));
1323     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1324     c1 -> Print(Form("%s/DCA_distance_outer_Y.pdf", sub_out_folder_name.c_str()));
1325     c1 -> Clear();
1326
1327     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1328     angle_diff_new -> SetMinimum(0);
1329     angle_diff_new -> Draw("hist");
1330     angle_diff_new_bkg_remove -> Draw("hist same");
1331     // ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s, peak : %f", plot_text.c_str(), peek));
1332     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1333     draw_text -> DrawLatex(0.4, 0.80, Form("#color[2]{Dist. StdDev: %.4f}", angle_diff_new_bkg_remove->GetStdDev()));
1334     c1 -> Print(Form("%s/angle_diff_new.pdf", sub_out_folder_name.c_str()));
1335     c1 -> Clear();
1336
1337     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1338     angle_diff_radian -> SetMinimum(0);
1339     angle_diff_radian -> Draw("hist");
1340     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1341     // gaus_fit_MC -> SetParameters(angle_diff_radian -> GetBinContent(angle_diff_radian -> GetMaximumBin()), angle_diff_radian -> GetMean(), angle_diff_radian -> GetStdDev() * 0.2, 0);
1342
1343     // angle_diff_radian -> Fit(gaus_fit_MC,"NQ","", -1. * (angle_diff_radian -> GetStdDev() * 0.4), (angle_diff_radian -> GetStdDev() * 0.4 ) );
1344     // gaus_fit_MC -> SetRange( (gaus_fit_MC -> GetParameter(1) - gaus_fit_MC -> GetParameter(2)), (gaus_fit_MC -> GetParameter(1) + gaus_fit_MC -> GetParameter(2)) );
1345
1346     // draw_text -> DrawLatex(0.2, 0.84, Form("Fit mean: %.4f", gaus_fit_MC -> GetParameter(1)));
1347     // draw_text -> DrawLatex(0.2, 0.80, Form("Fit width: %.4f", gaus_fit_MC -> GetParameter(2)));
1348
1349     double angle_diff_radian_bkg_level = get_dist_offset(angle_diff_radian, 15);
1350     gaus_pol2_fit -> SetParameters(
1351         angle_diff_radian -> GetBinContent(angle_diff_radian -> GetMaximumBin()) - angle_diff_radian_bkg_level, angle_diff_radian -> GetMean(), angle_diff_radian -> GetStdDev() * 0.2,
1352         angle_diff_radian_bkg_level, 0, -0.003490, 0
1353     );
1354     angle_diff_radian -> Fit(gaus_pol2_fit,"NQ");
1355     gaus_pol2_fit -> Draw("l same");
1356
1357     // gaus_fit_MC -> Draw("l same");
1358     c1 -> Print(Form("%s/angle_diff_radian.pdf", sub_out_folder_name.c_str()));
1359     c1 -> Clear();
1360
1361     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1362     angle_diff_inner_phi_peak -> SetStats(0);
1363     angle_diff_inner_phi_peak -> GetXaxis() -> SetTitle( angle_diff_inner_phi -> GetXaxis() -> GetTitle() );
1364     angle_diff_inner_phi_peak -> GetYaxis() -> SetTitle( angle_diff_inner_phi -> GetYaxis() -> GetTitle() );
1365     angle_diff_inner_phi_peak -> Draw("colz0");
1366     angle_diff_inner_phi_peak_profile -> Draw("same");
1367     horizontal_fit_angle_diff_inner -> Draw("l same");
1368     // gaus_fit_angle_diff -> Draw("lsame");
1369     // ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s, peak : %f", plot_text.c_str(), peek));
1370     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1371     draw_text -> DrawLatex(0.25, 0.84, Form("#color[2]{Assumed vertex: %.4f cm, %.4f cm}", current_vtxX/10., current_vtxY/10.));
1372     draw_text -> DrawLatex(0.25, 0.80, Form("#color[2]{Pol0 fit chi2/NDF: %.3f, fit error: %.4f}", horizontal_fit_angle_diff_inner -> GetChisquare() / double(horizontal_fit_angle_diff_inner -> GetNDF()), horizontal_fit_angle_diff_inner->GetParError(0)));
1373     c1 -> Print(Form("%s/angle_diff_inner_phi_peak.pdf", sub_out_folder_name.c_str()));
1374     c1 -> Clear();
1375
1376     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1377     angle_diff_outer_phi_peak -> SetStats(0);
1378     angle_diff_outer_phi_peak -> GetXaxis() -> SetTitle( angle_diff_outer_phi -> GetXaxis() -> GetTitle() );
1379     angle_diff_outer_phi_peak -> GetYaxis() -> SetTitle( angle_diff_outer_phi -> GetYaxis() -> GetTitle() );
1380     angle_diff_outer_phi_peak -> Draw("colz0");
1381     angle_diff_outer_phi_peak_profile -> Draw("same");
1382     horizontal_fit_angle_diff_outer -> Draw("l same");
1383     // ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s, peak : %f", plot_text.c_str(), peek));
1384     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1385     draw_text -> DrawLatex(0.25, 0.84, Form("#color[2]{Assumed vertex: %.4f cm, %.4f cm}", current_vtxX/10., current_vtxY/10.));
1386     draw_text -> DrawLatex(0.25, 0.80, Form("#color[2]{Pol0 fit chi2/NDF: %.3f, fit error: %.4f}", horizontal_fit_angle_diff_outer -> GetChisquare() / double(horizontal_fit_angle_diff_outer -> GetNDF()), horizontal_fit_angle_diff_outer->GetParError(0)));
1387     c1 -> Print(Form("%s/angle_diff_outer_phi_peak.pdf", sub_out_folder_name.c_str()));
1388     c1 -> Clear();
1389 }
1390
1391 TH1F * INTTXYvtx::PrintPlots_bkgrm_helper(TH1F * hist_in, double signal_region)
1392 {
1393     c1 -> cd();
1394     hist_in -> SetMinimum(0);
1395     hist_in -> SetMaximum( hist_in -> GetBinContent( hist_in -> GetMaximumBin() ) * 1.8);
1396     hist_in -> Draw("hist");
1397     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1398
1399     double hist_in_bkg_level = get_dist_offset(hist_in, 15);
1400
1401     bkg_fit_pol2 -> SetParameters(hist_in_bkg_level, 0, -0.0034906585, 0, signal_region);
1402     bkg_fit_pol2 -> FixParameter(4, signal_region);
1403     bkg_fit_pol2 -> SetParLimits(2, -100, 0);
1404     hist_in -> Fit(bkg_fit_pol2,"NQ");
1405
1406     draw_pol2_line -> SetParameters(bkg_fit_pol2 -> GetParameter(0), bkg_fit_pol2 -> GetParameter(1), bkg_fit_pol2 -> GetParameter(2), bkg_fit_pol2 -> GetParameter(3));
1407
1408     TH1F * hist_in_bkgrm = (TH1F*)hist_in -> Clone();
1409     hist_in_bkgrm -> SetLineColor(8);
1410     for (int i = 0; i < hist_in_bkgrm -> GetNbinsX(); i++)
1411     {
1412         double rest_bincontent = hist_in -> GetBinContent(i+1) - draw_pol2_line -> Eval(hist_in -> GetBinCenter(i+1));
1413         rest_bincontent = (rest_bincontent < 0) ? 0 : rest_bincontent;
1414         hist_in_bkgrm -> SetBinContent(i+1, rest_bincontent);
1415     }
1416
1417     gaus_fit_MC -> SetParameters(hist_in_bkgrm -> GetBinContent(hist_in_bkgrm -> GetMaximumBin()), hist_in_bkgrm -> GetMean(), hist_in_bkgrm -> GetStdDev() * 0.2, 0);
1418     hist_in_bkgrm -> Fit(gaus_fit_MC,"NQ","", -1. * (hist_in_bkgrm -> GetStdDev() * 0.7), (hist_in_bkgrm -> GetStdDev() * 0.7 ));
1419     gaus_fit_MC -> SetRange(gaus_fit_MC -> GetParameter(1) - gaus_fit_MC -> GetParameter(2) * 2., gaus_fit_MC -> GetParameter(1) + gaus_fit_MC -> GetParameter(2) * 2.);
1420
1421     draw_text -> DrawLatex(0.21, 0.86, Form("Fit mean: %.6f", gaus_fit_MC -> GetParameter(1)));
1422     draw_text -> DrawLatex(0.21, 0.82, Form("Fit width: %.6f", gaus_fit_MC -> GetParameter(2)));
1423     draw_text -> DrawLatex(0.21, 0.78, Form("Hist StdDev: %.6f", hist_in_bkgrm -> GetStdDev()));
1424
1425     draw_pol2_line -> Draw("l same");
1426     hist_in_bkgrm -> Draw("hist same");
1427     gaus_fit_MC -> Draw("l same");
1428
1429     c1 -> Print(Form("%s/%s.pdf", out_folder_directory.c_str(), hist_in -> GetName()));
1430     c1 -> Clear();
1431
1432     return hist_in_bkgrm;
1433 }
1434
1435 void INTTXYvtx::RunDeltaPhiStudy()
1436 {
1437     for (int i = 0; i < cluster_pair_vec.size(); i++)
1438     {
1439         if ( i % 100000 == 0 ) {cout<<" In RunDeltaPhiStudy, i : "<<i<<endl;}
1440
1441         Clus_InnerPhi_Offset = (cluster_pair_vec[i].first.y - beam_origin.second < 0) ? atan2(cluster_pair_vec[i].first.y - beam_origin.second, cluster_pair_vec[i].first.x - beam_origin.first) * (180./TMath::Pi()) + 360 : atan2(cluster_pair_vec[i].first.y - beam_origin.second, cluster_pair_vec[i].first.x - beam_origin.first) * (180./TMath::Pi());
1442         Clus_OuterPhi_Offset = (cluster_pair_vec[i].second.y - beam_origin.second < 0) ? atan2(cluster_pair_vec[i].second.y - beam_origin.second, cluster_pair_vec[i].second.x - beam_origin.first) * (180./TMath::Pi()) + 360 : atan2(cluster_pair_vec[i].second.y - beam_origin.second, cluster_pair_vec[i].second.x - beam_origin.first) * (180./TMath::Pi());
1443
1444         // note : to have the radian
1445         Clus_InnerPhi_Offset_radian = atan2(cluster_pair_vec[i].first.y  - beam_origin.second, cluster_pair_vec[i].first.x  - beam_origin.first);
1446         Clus_OuterPhi_Offset_radian = atan2(cluster_pair_vec[i].second.y - beam_origin.second, cluster_pair_vec[i].second.x - beam_origin.first);
1447
1448         double DCA_sign = calculateAngleBetweenVectors(
1449             cluster_pair_vec[i].second.x, cluster_pair_vec[i].second.y,
1450             cluster_pair_vec[i].first.x, cluster_pair_vec[i].first.y,
1451             beam_origin.first, beam_origin.second
1452         );
1453
1454         final_angle_diff_inner_phi_degree_coarseX -> Fill(Clus_InnerPhi_Offset, Clus_InnerPhi_Offset - Clus_OuterPhi_Offset);
1455         final_angle_diff_inner_phi_radian_coarseX -> Fill(Clus_InnerPhi_Offset_radian, Clus_InnerPhi_Offset_radian - Clus_OuterPhi_Offset_radian);
1456         final_angle_diff_radian_before            -> Fill(Clus_InnerPhi_Offset_radian - Clus_OuterPhi_Offset_radian);
1457
1458         final_angle_diff_radian_DCA_2D_before -> Fill(Clus_InnerPhi_Offset_radian - Clus_OuterPhi_Offset_radian, DCA_sign * 0.1);
1459
1460         final_DCA_cm_before -> Fill(DCA_sign * 0.1);
1461         final_DCA_cm_inner_phi_radian_coarseX -> Fill(Clus_InnerPhi_Offset_radian, DCA_sign * 0.1);
1462         final_DCA_mm_inner_phi_degree_coarseX -> Fill(Clus_InnerPhi_Offset, DCA_sign);
1463
1464     } // note : end of the loop for the cluster pair
1465
1466     final_angle_diff_inner_phi_degree_coarseX_peak = (TH2F*)final_angle_diff_inner_phi_degree_coarseX -> Clone();
1467     TH2F_threshold_advanced(final_angle_diff_inner_phi_degree_coarseX_peak, 0.5); // todo : threshold ratio can be modified here
1468     final_angle_diff_inner_phi_degree_coarseX_peak_profile =  final_angle_diff_inner_phi_degree_coarseX_peak->ProfileX("final_angle_diff_inner_phi_degree_coarseX_peak_profile");
1469
1470     final_angle_diff_inner_phi_radian_coarseX_peak = (TH2F*)final_angle_diff_inner_phi_radian_coarseX -> Clone();
1471     TH2F_threshold_advanced(final_angle_diff_inner_phi_radian_coarseX_peak, 0.5); // todo : threshold ratio can be modified here
1472     final_angle_diff_inner_phi_radian_coarseX_peak_profile =  final_angle_diff_inner_phi_radian_coarseX_peak->ProfileX("final_angle_diff_inner_phi_radian_coarseX_peak_profile");
1473
1474     final_DCA_cm_inner_phi_radian_coarseX_peak = (TH2F*)final_DCA_cm_inner_phi_radian_coarseX -> Clone();
1475     TH2F_threshold_advanced(final_DCA_cm_inner_phi_radian_coarseX_peak, 0.5); // todo : threshold ratio can be modified here
1476     final_DCA_cm_inner_phi_radian_coarseX_peak_profile =  final_DCA_cm_inner_phi_radian_coarseX_peak->ProfileX("final_DCA_cm_inner_phi_radian_coarseX_peak_profile");
1477
1478     final_DCA_mm_inner_phi_degree_coarseX_peak = (TH2F*)final_DCA_mm_inner_phi_degree_coarseX -> Clone();
1479     TH2F_threshold_advanced(final_DCA_mm_inner_phi_degree_coarseX_peak, 0.5); // todo : threshold ratio can be modified here
1480     final_DCA_mm_inner_phi_degree_coarseX_peak_profile =  final_DCA_mm_inner_phi_degree_coarseX_peak->ProfileX("final_DCA_mm_inner_phi_degree_coarseX_peak_profile");
1481
1482     cout<<"// for the radian case"<<endl;
1483     cout<<"vector<double> angle_diff_correction_radian = {";
1484     for (int i = 0; i < final_angle_diff_inner_phi_radian_coarseX_peak_profile -> GetNbinsX(); i++)
1485     {
1486         if (i == final_angle_diff_inner_phi_radian_coarseX_peak_profile -> GetNbinsX() - 1)
1487         {
1488             cout<<final_angle_diff_inner_phi_radian_coarseX_peak_profile -> GetBinContent(i+1)<<"};"<<endl;
1489             break;
1490         }
1491         cout<<final_angle_diff_inner_phi_radian_coarseX_peak_profile -> GetBinContent(i+1)<<", ";
1492     }
1493
1494     cout<<"// for the degree case"<<endl;
1495     cout<<"vector<double> angle_diff_correction_degree = {";
1496     for (int i = 0; i < final_angle_diff_inner_phi_degree_coarseX_peak_profile -> GetNbinsX(); i++)
1497     {
1498         if (i == final_angle_diff_inner_phi_degree_coarseX_peak_profile -> GetNbinsX() - 1)
1499         {
1500             cout<<final_angle_diff_inner_phi_degree_coarseX_peak_profile -> GetBinContent(i+1)<<"};"<<endl;
1501             break;
1502         }
1503         cout<<final_angle_diff_inner_phi_degree_coarseX_peak_profile -> GetBinContent(i+1)<<", ";
1504     }
1505
1506     // note : after the correction
1507     for (int i = 0; i < cluster_pair_vec.size(); i++)
1508     {
1509         // note : to have the radian
1510         Clus_InnerPhi_Offset_radian = atan2(cluster_pair_vec[i].first.y  - beam_origin.second, cluster_pair_vec[i].first.x  - beam_origin.first);
1511         Clus_OuterPhi_Offset_radian = atan2(cluster_pair_vec[i].second.y - beam_origin.second, cluster_pair_vec[i].second.x - beam_origin.first);
1512
1513         double angle_correction_radian = final_angle_diff_inner_phi_radian_coarseX_peak_profile -> GetBinContent( final_angle_diff_inner_phi_radian_coarseX_peak_profile -> FindBin(Clus_InnerPhi_Offset_radian) );
1514         double eventual_angle_diff_radian = (Clus_InnerPhi_Offset_radian - Clus_OuterPhi_Offset_radian) - angle_correction_radian;
1515
1516         // note : unit : mm
1517         double DCA_sign = calculateAngleBetweenVectors(
1518             cluster_pair_vec[i].second.x, cluster_pair_vec[i].second.y,
1519             cluster_pair_vec[i].first.x, cluster_pair_vec[i].first.y,
1520             beam_origin.first, beam_origin.second
1521         );
1522
1523         double DCA_sign_cm_correction = DCA_sign * 0.1 - final_DCA_cm_inner_phi_radian_coarseX_peak_profile -> GetBinContent( final_DCA_cm_inner_phi_radian_coarseX_peak_profile -> FindBin(Clus_InnerPhi_Offset_radian) );
1524
1525         final_angle_diff_radian_post -> Fill( eventual_angle_diff_radian );
1526
1527         final_angle_diff_radian_DCA_2D_post -> Fill(eventual_angle_diff_radian, DCA_sign_cm_correction);
1528
1529         final_DCA_cm_post -> Fill(DCA_sign_cm_correction);
1530
1531     } // note : end of the loop for the cluster pair
1532
1533     // note : print the plots
1534
1535     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1536     final_DCA_cm_inner_phi_radian_coarseX -> Draw("colz0");
1537     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1538     c1 -> Print(Form("%s/final_DCA_cm_inner_phi_radian_coarseX.pdf", out_folder_directory.c_str()));
1539     c1 -> Clear();
1540
1541     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1542     final_DCA_cm_inner_phi_radian_coarseX_peak -> Draw("colz0");
1543     final_DCA_cm_inner_phi_radian_coarseX_peak_profile -> Draw("same");
1544     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1545     c1 -> Print(Form("%s/final_DCA_cm_inner_phi_radian_coarseX_peak.pdf", out_folder_directory.c_str()));
1546     c1 -> Clear();
1547
1548
1549
1550     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1551     final_DCA_mm_inner_phi_degree_coarseX -> Draw("colz0");
1552     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1553     c1 -> Print(Form("%s/final_DCA_mm_inner_phi_degree_coarseX.pdf", out_folder_directory.c_str()));
1554     c1 -> Clear();
1555
1556     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1557     final_DCA_mm_inner_phi_degree_coarseX_peak -> Draw("colz0");
1558     final_DCA_mm_inner_phi_degree_coarseX_peak_profile -> Draw("same");
1559     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1560     c1 -> Print(Form("%s/final_DCA_mm_inner_phi_degree_coarseX_peak.pdf", out_folder_directory.c_str()));
1561     c1 -> Clear();
1562
1563
1564     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1565     final_angle_diff_inner_phi_degree_coarseX -> Draw("colz0");
1566     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1567     c1 -> Print(Form("%s/final_angle_diff_inner_phi_degree_coarseX.pdf", out_folder_directory.c_str()));
1568     c1 -> Clear();
1569
1570     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1571     final_angle_diff_inner_phi_degree_coarseX_peak -> Draw("colz0");
1572     final_angle_diff_inner_phi_degree_coarseX_peak_profile -> Draw("same");
1573     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1574     c1 -> Print(Form("%s/final_angle_diff_inner_phi_degree_coarseX_peak.pdf", out_folder_directory.c_str()));
1575     c1 -> Clear();
1576
1577     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1578     final_angle_diff_inner_phi_radian_coarseX -> Draw("colz0");
1579     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1580     c1 -> Print(Form("%s/final_angle_diff_inner_phi_radian_coarseX.pdf", out_folder_directory.c_str()));
1581     c1 -> Clear();
1582
1583     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1584     final_angle_diff_inner_phi_radian_coarseX_peak -> Draw("colz0");
1585     final_angle_diff_inner_phi_radian_coarseX_peak_profile -> Draw("same");
1586     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1587     c1 -> Print(Form("%s/final_angle_diff_inner_phi_radian_coarseX_peak.pdf", out_folder_directory.c_str()));
1588     c1 -> Clear();
1589
1590     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1591     final_angle_diff_radian_DCA_2D_before -> Draw("colz0");
1592     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1593     c1 -> Print(Form("%s/final_angle_diff_radian_DCA_2D_before.pdf", out_folder_directory.c_str()));
1594     c1 -> Clear();
1595
1596     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1597     final_angle_diff_radian_DCA_2D_post -> Draw("colz0");
1598     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1599     c1 -> Print(Form("%s/final_angle_diff_radian_DCA_2D_post.pdf", out_folder_directory.c_str()));
1600     c1 -> Clear();
1601
1602     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1603     TH1F * final_angle_diff_radian_before_bkgrm = PrintPlots_bkgrm_helper(final_angle_diff_radian_before, ana_map_version::signal_region / radian_correction);
1604
1605     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1606     TH1F * final_angle_diff_radian_post_bkgrm = PrintPlots_bkgrm_helper(final_angle_diff_radian_post, ana_map_version::signal_region / radian_correction);
1607
1608     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1609     TLegend * legend = new TLegend(0.2,0.75,0.45,0.9);
1610     // legend -> SetMargin(0);
1611     legend->SetTextSize(0.03);
1612
1613     final_angle_diff_radian_before_bkgrm -> SetLineColor(4);
1614     final_angle_diff_radian_before_bkgrm -> Draw("hist");
1615     final_angle_diff_radian_post_bkgrm -> SetLineColor(2);
1616     final_angle_diff_radian_post_bkgrm -> Draw("hist same");
1617
1618     legend -> AddEntry(final_angle_diff_radian_before_bkgrm, "Original", "f");
1619     legend -> AddEntry(final_angle_diff_radian_post_bkgrm, "#Delta#phi mean shift", "f");
1620
1621     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1622     legend -> Draw("same");
1623
1624     c1 -> Print(Form("%s/final_angle_diff_radian_comp.pdf", out_folder_directory.c_str()));
1625     c1 -> Clear();
1626     legend -> Clear();
1627
1628
1629     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1630     TH1F * final_DCA_cm_before_bkgrm = PrintPlots_bkgrm_helper(final_DCA_cm_before, 0.03);
1631
1632     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1633     TH1F * final_DCA_cm_post_bkgrm = PrintPlots_bkgrm_helper(final_DCA_cm_post, 0.03);
1634
1635     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1636     final_DCA_cm_before_bkgrm -> SetLineColor(4);
1637     final_DCA_cm_before_bkgrm -> Draw("hist");
1638     final_DCA_cm_post_bkgrm -> SetLineColor(2);
1639     final_DCA_cm_post_bkgrm -> Draw("hist same");
1640
1641     legend -> AddEntry(final_DCA_cm_before_bkgrm, "Original", "f");
1642     legend -> AddEntry(final_DCA_cm_post_bkgrm, "DCA mean shift", "f");
1643
1644     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1645     legend -> Draw("same");
1646
1647     c1 -> Print(Form("%s/final_DCA_cm_comp.pdf", out_folder_directory.c_str()));
1648     c1 -> Clear();
1649     legend -> Clear();
1650
1651     cout<<"End of RunDeltaPhiStudy"<<endl;
1652 }
1653
1654 void INTTXYvtx::ClearHist(int print_option)
1655 {
1656     angle_correlation -> Reset("ICESM");
1657     angle_diff_DCA_dist -> Reset("ICESM");
1658     angle_diff -> Reset("ICESM");
1659     DCA_point -> Reset("ICESM");
1660     DCA_distance -> Reset("ICESM");
1661     DCA_distance_inner_X -> Reset("ICESM");
1662     DCA_distance_inner_Y -> Reset("ICESM");
1663     DCA_distance_outer_X -> Reset("ICESM");
1664     DCA_distance_outer_Y -> Reset("ICESM");
1665
1666     DCA_distance_inner_phi -> Reset("ICESM");
1667     // DCA_distance_inner_phi_peak -> Reset("ICESM");
1668     DCA_distance_outer_phi -> Reset("ICESM");
1669     // DCA_distance_outer_phi_peak -> Reset("ICESM");
1670
1671     angle_diff_inner_phi -> Reset("ICESM");
1672     angle_diff_outer_phi -> Reset("ICESM");
1673
1674     // angle_diff_inner_phi_peak -> Reset("ICESM");
1675     // angle_diff_outer_phi_peak -> Reset("ICESM");
1676
1677     angle_diff_new -> Reset("ICESM");
1678     angle_diff_radian -> Reset("ICESM");
1679     // angle_diff_new_bkg_remove -> Reset("ICESM");
1680
1681     delete angle_diff_new_bkg_remove;
1682     delete DCA_distance_inner_phi_peak;
1683     delete DCA_distance_outer_phi_peak;
1684     delete angle_diff_inner_phi_peak;
1685     delete angle_diff_outer_phi_peak;
1686 }
1687
1688 void INTTXYvtx::EndRun()
1689 {
1690     inner_pos_xy -> Reset("ICESM");
1691     outer_pos_xy -> Reset("ICESM");
1692     inner_outer_pos_xy -> Reset("ICESM");
1693     N_cluster_correlation -> Reset("ICESM");
1694     N_cluster_correlation_close -> Reset("ICESM");
1695
1696     // file_out -> cd();
1697     // tree_out -> SetDirectory(file_out);
1698     // tree_out -> Write();
1699     // file_out -> Close();
1700
1701     // c1 -> Delete();
1702     // ltx -> Delete();
1703     // draw_text -> Delete();
1704     // cos_fit -> Delete();
1705     // gaus_fit -> Delete();
1706     // gaus_fit_angle_diff -> Delete();
1707
1708     // horizontal_fit_inner -> Delete();
1709     // horizontal_fit_angle_diff_inner -> Delete();
1710     // horizontal_fit_outer -> Delete();
1711     // horizontal_fit_angle_diff_outer -> Delete();
1712
1713     delete gROOT->FindObject("angle_diff");
1714     delete gROOT->FindObject("angle_diff_new");
1715
1716     // angle_correlation -> Delete();
1717     // angle_diff_DCA_dist -> Delete();
1718     // angle_diff -> Delete();
1719     // angle_diff_new -> Delete();
1720     // angle_diff_new_bkg_remove -> Delete();
1721     // inner_pos_xy -> Delete();
1722     // outer_pos_xy -> Delete();
1723     // inner_outer_pos_xy -> Delete();
1724     // DCA_point -> Delete();
1725     // DCA_distance -> Delete();
1726     // N_cluster_correlation -> Delete();
1727     // N_cluster_correlation_close -> Delete();
1728
1729     // DCA_distance_inner_X -> Delete();
1730     // DCA_distance_inner_Y -> Delete();
1731     // DCA_distance_outer_X -> Delete();
1732     // DCA_distance_outer_Y -> Delete();
1733
1734     // DCA_distance_inner_phi -> Delete();
1735     // DCA_distance_inner_phi_peak -> Delete();
1736     // DCA_distance_inner_phi_peak_profile -> Delete();
1737     // DCA_distance_outer_phi -> Delete();
1738     // DCA_distance_outer_phi_peak -> Delete();
1739     // DCA_distance_outer_phi_peak_profile -> Delete();
1740
1741     // angle_diff_inner_phi -> Delete();
1742     // angle_diff_inner_phi_peak -> Delete();
1743     // angle_diff_inner_phi_peak_profile -> Delete();
1744     // angle_diff_outer_phi -> Delete();
1745     // angle_diff_outer_phi_peak -> Delete();
1746     // angle_diff_outer_phi_peak_profile -> Delete();
1747
1748     // DCA_distance_inner_phi_peak_final -> Delete();
1749     // angle_diff_inner_phi_peak_final -> Delete();
1750
1751     return;
1752 }
1753
1754 void INTTXYvtx::TH2F_threshold(TH2F * hist, double threshold)
1755 {
1756     double max_cut = hist -> GetMaximum() * threshold;
1757
1758     for (int xi = 0; xi < hist -> GetNbinsX(); xi++){
1759         for(int yi = 0; yi < hist -> GetNbinsY(); yi++){
1760             if (hist -> GetBinContent(xi + 1, yi +1) < max_cut){ hist -> SetBinContent(xi + 1, yi +1, 0); }
1761         }
1762     }
1763 }
1764
1765 void INTTXYvtx::TH2F_threshold_advanced(TH2F * hist, double threshold)
1766 {
1767     double big_num, max_cut;
1768
1769     // note : this function removes the background of the 2D histogram by the following method
1770     // note : 1. find the maximum bin content of each column and then time with the threshold ratio
1771     // note : so the background removal is performed column by column
1772
1773     for (int xi = 0; xi < hist -> GetNbinsX(); xi++){
1774
1775         for(int yi = 0; yi < hist -> GetNbinsY(); yi++){
1776             if (yi == 0) {big_num = hist -> GetBinContent(xi + 1, yi +1);}
1777             else
1778             {
1779                 if (hist -> GetBinContent(xi + 1, yi +1) > big_num) {big_num = hist -> GetBinContent(xi + 1, yi +1);}
1780             }
1781         }
1782
1783         max_cut = big_num * threshold;
1784
1785         for(int yi = 0; yi < hist -> GetNbinsY(); yi++){
1786             if (hist -> GetBinContent(xi + 1, yi +1) < max_cut){ hist -> SetBinContent(xi + 1, yi +1, 0); }
1787         }
1788     }
1789 }
1790
1791 void INTTXYvtx::TH2F_threshold_advanced_row(TH2F * hist, double threshold)
1792 {
1793     double big_num, max_cut;
1794
1795     // note : this function removes the background of the 2D histogram by the following method
1796     // note : 1. find the maximum bin content of each row and then time with the threshold ratio
1797     // note : so the background removal is performed row by row
1798
1799     for (int yi = 0; yi < hist -> GetNbinsY(); yi++)
1800     {
1801         for(int xi = 0; xi < hist -> GetNbinsX(); xi++)
1802         {
1803             if (xi == 0) {big_num = hist -> GetBinContent(xi + 1, yi +1);}
1804             else
1805             {
1806                 if (hist -> GetBinContent(xi + 1, yi +1) > big_num) {big_num = hist -> GetBinContent(xi + 1, yi +1);}
1807             }
1808         }
1809
1810         max_cut = big_num * threshold;
1811
1812         for(int xi = 0; xi < hist -> GetNbinsX(); xi++)
1813         {
1814             if (hist -> GetBinContent(xi + 1, yi +1) < max_cut){ hist -> SetBinContent(xi + 1, yi +1, 0); }
1815         }
1816     }
1817 }
1818
1819 void INTTXYvtx::TH2F_threshold_advanced_2(TH2F * hist, double threshold)
1820 {
1821     // note : this function is to remove the background of the 2D histogram
1822     // note : but the threshold is given by average of the contents of the top "chosen_bin" bins and timing the threshold
1823     double max_cut = 0;
1824     int chosen_bin = 7;
1825
1826     vector<float> all_bin_content_vec; all_bin_content_vec.clear();
1827     for (int xi = 0; xi < hist -> GetNbinsX(); xi++){
1828         for(int yi = 0; yi < hist -> GetNbinsY(); yi++){
1829             all_bin_content_vec.push_back(hist -> GetBinContent(xi + 1, yi +1));
1830         }
1831     }
1832     vector<unsigned long> ind(all_bin_content_vec.size(),0);
1833     TMath::Sort(all_bin_content_vec.size(), &all_bin_content_vec[0], &ind[0]);
1834     for (int i = 0; i < chosen_bin; i++) {max_cut += all_bin_content_vec[ind[i]]; /*cout<<"test : "<<all_bin_content_vec[ind[i]]<<endl;*/}
1835
1836     max_cut = (max_cut / double(chosen_bin)) * threshold;
1837
1838     for (int xi = 0; xi < hist -> GetNbinsX(); xi++){
1839         for(int yi = 0; yi < hist -> GetNbinsY(); yi++){
1840             if (hist -> GetBinContent(xi + 1, yi +1) < max_cut){ hist -> SetBinContent(xi + 1, yi +1, 0); }
1841         }
1842     }
1843 }
1844
1845 double INTTXYvtx::get_radius(double x, double y)
1846 {
1847     return sqrt(pow(x,2)+pow(y,2));
1848 }
1849
1850 type_pos INTTXYvtx::PolarToCartesian(double radius, double angleDegrees) {
1851     // Convert angle from degrees to radians
1852     double angleRadians = angleDegrees * M_PI / 180.0;
1853
1854     // Calculate Cartesian coordinates
1855     double x = radius * cos(angleRadians);
1856     double y = radius * sin(angleRadians);
1857
1858     return {x, y};
1859 }
1860
1861 void INTTXYvtx::Characterize_Pad(TPad *pad, float left = 0.15, float right = 0.1, float top = 0.1, float bottom = 0.12, bool set_logY = false, int setgrid_bool = 0)
1862 {
1863     if (setgrid_bool == true) {pad -> SetGrid (1, 1);}
1864     pad -> SetLeftMargin   (left);
1865     pad -> SetRightMargin  (right);
1866     pad -> SetTopMargin    (top);
1867     pad -> SetBottomMargin (bottom);
1868     pad -> SetTicks(1,1);
1869     if (set_logY == true)
1870     {
1871         pad -> SetLogy (1);
1872     }
1873
1874 }
1875
1876 std::vector<double> INTTXYvtx::calculateDistanceAndClosestPoint(double x1, double y1, double x2, double y2, double target_x, double target_y)
1877 {
1878
1879     if (x1 != x2)
1880     {
1881         // Calculate the slope and intercept of the line passing through (x1, y1) and (x2, y2)
1882         double a = (y2 - y1) / (x2 - x1);
1883         double b = y1 - a * x1;
1884
1885         // cout<<"slope : y="<<a<<"x+"<<b<<endl;
1886
1887         // Calculate the closest distance from (target_x, target_y) to the line y = ax + b
1888         double closest_distance = std::abs(a * target_x - target_y + b) / std::sqrt(a * a + 1);
1889
1890         // Calculate the coordinates of the closest point (Xc, Yc) on the line y = ax + b
1891         double Xc = (target_x + a * target_y - a * b) / (a * a + 1);
1892         double Yc = a * Xc + b;
1893
1894         return { closest_distance, Xc, Yc };
1895     }
1896     else
1897     {
1898         double closest_distance = std::abs(x1 - target_x);
1899         double Xc = x1;
1900         double Yc = target_y;
1901
1902         return { closest_distance, Xc, Yc };
1903     }
1904
1905
1906 }
1907
1908 // note : Function to calculate the angle between two vectors in degrees using the cross product
1909 double INTTXYvtx::calculateAngleBetweenVectors(double x1, double y1, double x2, double y2, double targetX, double targetY) {
1910     // Calculate the vectors vector_1 (point_1 to point_2) and vector_2 (point_1 to target)
1911     double vector1X = x2 - x1;
1912     double vector1Y = y2 - y1;
1913
1914     double vector2X = targetX - x1;
1915     double vector2Y = targetY - y1;
1916
1917     // Calculate the cross product of vector_1 and vector_2 (z-component)
1918     double crossProduct = vector1X * vector2Y - vector1Y * vector2X;
1919
1920     // cout<<" crossProduct : "<<crossProduct<<endl;
1921
1922     // Calculate the magnitudes of vector_1 and vector_2
1923     double magnitude1 = std::sqrt(vector1X * vector1X + vector1Y * vector1Y);
1924     double magnitude2 = std::sqrt(vector2X * vector2X + vector2Y * vector2Y);
1925
1926     // Calculate the angle in radians using the inverse tangent of the cross product and dot product
1927     double dotProduct = vector1X * vector2X + vector1Y * vector2Y;
1928
1929     double angleInRadians = std::atan2(std::abs(crossProduct), dotProduct);
1930     // Convert the angle from radians to degrees and return it
1931     double angleInDegrees = angleInRadians * 180.0 / M_PI;
1932
1933     double angleInRadians_new = std::asin( crossProduct/(magnitude1*magnitude2) );
1934     double angleInDegrees_new = angleInRadians_new * 180.0 / M_PI;
1935
1936     // cout<<"angle : "<<angleInDegrees_new<<endl;
1937
1938     double DCA_distance = sin(angleInRadians_new) * magnitude2;
1939
1940     return DCA_distance;
1941 }
1942
1943
1944 pair<type_pos,type_pos> INTTXYvtx::GetTangentPointsAtCircle(double CenterX, double CenterY, double R, double XX, double YY) {
1945
1946     double XT0;
1947     double YT0;
1948     double XT1;
1949     double YT1;
1950
1951     if (R == 0) {
1952         // This behavior can be modified
1953         cout<<"In INTTXYvtx, the input vertex is zero"<<endl;
1954         return {{-999.,-999.},{-999.,-999.}};
1955     }
1956
1957     double nx = (XX - CenterX) / R; // Shift and scale
1958     double ny = (YY - CenterY) / R;
1959     double xy = nx * nx + ny * ny;
1960
1961     if (abs(xy - 1.0) < 1e-10) { // Point lies at circumference, one tangent
1962         XT0 = XX;
1963         YT0 = YY;
1964         return {{XT0,YT0},{XT0,YT0}};
1965     }
1966
1967     if (xy < 1.0){ // Point lies inside the circle, no tangents
1968         cout<<"IN INTTXYvtx, the point is inside the circle, no tangnet points"<<endl;
1969         return {{-999.,-999.},{-999.,-999.}};
1970     }
1971
1972     // Common case, two tangents
1973     int Result = 2;
1974     double D = ny * sqrt(xy - 1);
1975     double tx0 = (nx - D) / xy;
1976     double tx1 = (nx + D) / xy;
1977
1978     if (abs(ny) > 1e-10) { // Common case
1979         YT0 = CenterY + R * (1 - tx0 * nx) / ny;
1980         YT1 = CenterY + R * (1 - tx1 * nx) / ny;
1981     } else { // Point at the center horizontal, Y=0
1982         D = R * sqrt(1 - tx0 * tx0);
1983         YT0 = CenterY + D;
1984         YT1 = CenterY - D;
1985     }
1986
1987     // Restore scale and position
1988     XT0 = CenterX + R * tx0;
1989     XT1 = CenterX + R * tx1;
1990
1991     return {{XT0,YT0},{XT1,YT1}};
1992 }
1993
1994 pair<type_pos,type_pos> INTTXYvtx::findIntersectionPoints(double A, double mu, double sigma, double B, double C) {
1995     // Solve for x in the equation A * e^(-((x - mu)^2) / (2 * sigma^2)) + B = C
1996     // This involves solving a quadratic equation, and the solutions will be the x values of intersection points
1997
1998     double delta = mu * mu - 2 * sigma * sigma * log((C - B) / A);
1999
2000     // Check if there are real solutions
2001     if (delta >= 0) {
2002         double x1 = mu + sqrt(delta);
2003         double x2 = mu - sqrt(delta);
2004
2005         // Calculate corresponding y values
2006         double y1 = A * exp(-pow(x1 - mu, 2) / (2 * sigma * sigma)) + B;
2007         double y2 = A * exp(-pow(x2 - mu, 2) / (2 * sigma * sigma)) + B;
2008
2009         return {{x1, y1}, {x2, y2}};
2010
2011     } else {
2012         cout << "No real intersection points." << endl;
2013         return {{-999.,-999.},{-999.,-999.}};
2014     }
2015 }
2016
2017
2018 //  note : N group, group size, group tag, group width ?  // note : {group size, group entry, group tag, group widthL, group widthR}
2019 // note : {N_group, ratio (if two), peak widthL, peak widthR}
2020 vector<double> INTTXYvtx::find_Ngroup(TH1F * hist_in)
2021 {
2022     double Highest_bin_Content  = hist_in -> GetBinContent(hist_in -> GetMaximumBin());
2023     double Highest_bin_Center   = hist_in -> GetBinCenter(hist_in -> GetMaximumBin());
2024
2025     int group_Nbin = 0;
2026     int peak_group_ID;
2027     double group_entry = 0;
2028     double peak_group_ratio;
2029     vector<int> group_Nbin_vec; group_Nbin_vec.clear();
2030     vector<double> group_entry_vec; group_entry_vec.clear();
2031     vector<double> group_widthL_vec; group_widthL_vec.clear();
2032     vector<double> group_widthR_vec; group_widthR_vec.clear();
2033
2034     for (int i = 0; i < hist_in -> GetNbinsX(); i++){
2035         // todo : the background rejection is here : Highest_bin_Content/2. for the time being
2036         double bin_content = ( hist_in -> GetBinContent(i+1) <= Highest_bin_Content/2.) ? 0. : ( hist_in -> GetBinContent(i+1) - Highest_bin_Content/2. );
2037
2038         if (bin_content != 0){
2039
2040             if (group_Nbin == 0) {
2041                 group_widthL_vec.push_back(hist_in -> GetBinCenter(i+1) - (hist_in -> GetBinWidth(i+1)/2.));
2042             }
2043
2044             group_Nbin += 1;
2045             group_entry += bin_content;
2046         }
2047         else if (bin_content == 0 && group_Nbin != 0){
2048             group_widthR_vec.push_back(hist_in -> GetBinCenter(i+1) - (hist_in -> GetBinWidth(i+1)/2.));
2049             group_Nbin_vec.push_back(group_Nbin);
2050             group_entry_vec.push_back(group_entry);
2051             group_Nbin = 0;
2052             group_entry = 0;
2053         }
2054     }
2055     if (group_Nbin != 0) {
2056         group_Nbin_vec.push_back(group_Nbin);
2057         group_entry_vec.push_back(group_entry);
2058         group_widthR_vec.push_back(hist_in -> GetXaxis()->GetXmax());
2059     } // note : the last group at the edge
2060
2061     // note : find the peak group
2062     for (int i = 0; i < group_Nbin_vec.size(); i++){
2063         if (group_widthL_vec[i] < Highest_bin_Center && Highest_bin_Center < group_widthR_vec[i]){
2064             peak_group_ID = i;
2065             break;
2066         }
2067     }
2068
2069
2070     peak_group_ratio = group_entry_vec[peak_group_ID] / (accumulate( group_entry_vec.begin(), group_entry_vec.end(), 0.0 ));
2071
2072     // for (int i = 0; i < group_Nbin_vec.size(); i++)
2073     // {
2074     //     cout<<" "<<endl;
2075     //     cout<<"group size : "<<group_Nbin_vec[i]<<endl;
2076     //     cout<<"group entry : "<<group_entry_vec[i]<<endl;
2077     //     cout<<group_widthL_vec[i]<<" "<<group_widthR_vec[i]<<endl;
2078     // }
2079
2080     // cout<<" "<<endl;
2081     // cout<<"N group : "<<group_Nbin_vec.size()<<endl;
2082     // cout<<"Peak group ID : "<<peak_group_ID<<endl;
2083     // cout<<"peak group width : "<<group_widthL_vec[peak_group_ID]<<" "<<group_widthR_vec[peak_group_ID]<<endl;
2084     // cout<<"ratio : "<<peak_group_ratio<<endl;
2085
2086     // note : {N_group, ratio (if two), peak widthL, peak widthR}
2087     return {double(group_Nbin_vec.size()), peak_group_ratio, group_widthL_vec[peak_group_ID], group_widthR_vec[peak_group_ID]};
2088 }
2089
2090 void INTTXYvtx::Hist_1D_bkg_remove(TH1F * hist_in, double factor)
2091 {
2092     // todo : N bins considered to be used in the background quantification
2093     vector<double> Nbin_content_vec; Nbin_content_vec.clear();
2094     for (int i = hist_in -> GetNbinsX() - 5; i < hist_in -> GetNbinsX(); i++){Nbin_content_vec.push_back(hist_in -> GetBinContent(i+1));}
2095     double bkg_level = accumulate( Nbin_content_vec.begin(), Nbin_content_vec.end(), 0.0 ) / Nbin_content_vec.size();
2096     // cout<<"test, bkg cut : "<<bkg_level * factor<<endl;
2097
2098     for (int i = 0; i < hist_in -> GetNbinsX(); i++){
2099         // note : the background rejection is here : bkg_level * 1.5 for the time being
2100         double bin_content = ( hist_in -> GetBinContent(i+1) <= bkg_level * factor) ? 0. : ( hist_in -> GetBinContent(i+1) - bkg_level * factor );
2101         hist_in -> SetBinContent(i+1, bin_content);
2102     }
2103 }
2104
2105 void INTTXYvtx::DrawTGraphErrors(vector<double> x_vec, vector<double> y_vec, vector<double> xE_vec, vector<double> yE_vec, string output_directory, vector<string> plot_name)
2106 {
2107     c1 -> cd();
2108
2109     TGraphErrors * g = new TGraphErrors(x_vec.size(), &x_vec[0], &y_vec[0], &xE_vec[0], &yE_vec[0]);
2110     g -> SetMarkerStyle(20);
2111     g -> SetMarkerSize(1.5);
2112     g -> SetMarkerColor(1);
2113     g -> SetLineColor(1);
2114     g -> GetXaxis() -> SetTitle(plot_name[1].c_str());
2115     g -> GetYaxis() -> SetTitle(plot_name[2].c_str());
2116     if (plot_name.size() == 4){g -> Draw(plot_name[3].c_str());}
2117     else {g -> Draw("AP");}
2118
2119     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
2120     c1 -> Print(Form("%s/%s.pdf", output_directory.c_str(), plot_name[0].c_str()));
2121     c1 -> Clear();
2122     g -> Delete();
2123 }
2124
2125 void INTTXYvtx::Draw2TGraph(vector<double> x1_vec, vector<double> y1_vec, vector<double> x2_vec, vector<double> y2_vec, string output_directory, vector<string> plot_name)
2126 {
2127     c1 -> cd();
2128
2129     c1 -> SetLogy(1);
2130
2131     TGraph * g1 = new TGraph(x1_vec.size(), &x1_vec[0], &y1_vec[0]);
2132     g1 -> SetMarkerStyle(5);
2133     g1 -> SetMarkerSize(1);
2134     g1 -> SetMarkerColor(1);
2135     g1 -> SetLineColor(1);
2136     g1 -> GetXaxis() -> SetTitle(plot_name[1].c_str());
2137     g1 -> GetYaxis() -> SetTitle(plot_name[2].c_str());
2138     g1 -> GetXaxis() -> SetNdivisions(505);
2139     g1 -> GetXaxis() -> SetLimits(-1, x1_vec[x1_vec.size()-1] + 1);
2140     g1 -> Draw("AP");
2141
2142     TGraph * g2 = new TGraph(x2_vec.size(), &x2_vec[0], &y2_vec[0]);
2143     g2 -> SetMarkerStyle(5);
2144     g2 -> SetMarkerSize(1);
2145     g2 -> SetMarkerColor(2);
2146     g2 -> SetLineColor(2);
2147     g2 -> Draw("PL same");
2148
2149     TLegend * legend = new TLegend(0.4,0.75,0.65,0.9);
2150     // legend -> SetMargin(0);
2151     legend->SetTextSize(0.03);
2152     legend -> AddEntry(g1, "Tested vertex candidates", "p");
2153     legend -> AddEntry(g2, "Better performed candidates", "p");
2154     legend -> Draw("same");
2155     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
2156     c1 -> Print(Form("%s/%s.pdf", output_directory.c_str(), plot_name[0].c_str()));
2157     c1 -> Clear();
2158     g1 -> Delete();
2159     g2 -> Delete();
2160
2161     c1 -> SetLogy(0);
2162
2163     return;
2164 }
2165
2166 void INTTXYvtx::DrawTH2F(TH2F * hist_in, string output_directory, vector<string> plot_name)
2167 {
2168     c1 -> cd();
2169
2170     hist_in -> SetStats(0);
2171     hist_in -> GetXaxis() -> SetTitle(plot_name[1].c_str());
2172     hist_in -> GetYaxis() -> SetTitle(plot_name[2].c_str());
2173     hist_in -> GetZaxis() -> SetTitle(plot_name[3].c_str());
2174     hist_in -> GetXaxis() -> SetNdivisions(505);
2175     hist_in -> Draw("colz0");
2176
2177
2178     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
2179     c1 -> Print(Form("%s/%s.pdf", output_directory.c_str(), plot_name[0].c_str()));
2180     c1 -> Clear();
2181 }
2182
2183 vector<double> INTTXYvtx::SumTH2FColumnContent(TH2F * hist_in)
2184 {
2185     vector<double> sum_vec; sum_vec.clear();
2186     for (int i = 0; i < hist_in -> GetNbinsX(); i++){
2187         double sum = 0;
2188         for (int j = 0; j < hist_in -> GetNbinsY(); j++){
2189             sum += hist_in -> GetBinContent(i+1, j+1);
2190         }
2191         sum_vec.push_back(sum);
2192     }
2193     return sum_vec;
2194 }
2195
2196 vector<double> INTTXYvtx::SumTH2FColumnContent_row(TH2F * hist_in)
2197 {
2198     vector<double> sum_vec; sum_vec.clear();
2199     for (int i = 0; i < hist_in -> GetNbinsY(); i++){
2200         double sum = 0;
2201         for (int j = 0; j < hist_in -> GetNbinsX(); j++){
2202             sum += hist_in -> GetBinContent(j+1, i+1);
2203         }
2204         sum_vec.push_back(sum);
2205     }
2206     return sum_vec;
2207 }
2208
2209
2210 vector<pair<double,double>> INTTXYvtx::Get4vtx(pair<double,double> origin, double length)
2211 {
2212     vector<pair<double,double>> unit_vtx = {{1,1},{-1,1},{-1,-1},{1,-1}};
2213     vector<pair<double,double>> vec_out; vec_out.clear();
2214
2215     for (pair i1 : unit_vtx)
2216     {
2217         vec_out.push_back({i1.first * length + origin.first, i1.second * length + origin.second});
2218     }
2219
2220     return vec_out;
2221 }
2222
2223 vector<pair<double,double>> INTTXYvtx::Get4vtxAxis(pair<double,double> origin, double length)
2224 {
2225     vector<pair<double,double>> unit_vtx = {{1,0},{-1,0},{0,1},{0,-1}}; //note : X axis right and left, Y axis up and down
2226     vector<pair<double,double>> vec_out; vec_out.clear();
2227
2228     for (pair i1 : unit_vtx)
2229     {
2230         vec_out.push_back({i1.first * length + origin.first, i1.second * length + origin.second});
2231     }
2232
2233     return vec_out;
2234 }
2235
2236 void INTTXYvtx::TH2F_FakeClone(TH2F*hist_in, TH2F*hist_out)
2237 {
2238     if (hist_in -> GetNbinsX() != hist_out -> GetNbinsX() || hist_in -> GetNbinsY() != hist_out -> GetNbinsY())
2239     {
2240         cout<<"In INTTXYvtx::TH2F_FakeClone, the input and output histogram have different binning!"<<endl;
2241         return;
2242     }
2243
2244     for (int i = 0; i < hist_in -> GetNbinsX(); i++){
2245         for (int j = 0; j < hist_in -> GetNbinsY(); j++){
2246             hist_out -> SetBinContent(i+1, j+1, hist_in -> GetBinContent(i+1, j+1));
2247         }
2248     }
2249 }
2250
2251 void INTTXYvtx::TH1F_FakeClone(TH1F*hist_in, TH1F*hist_out)
2252 {
2253     if (hist_in -> GetNbinsX() != hist_out -> GetNbinsX())
2254     {
2255         cout<<"In INTTXYvtx::TH1F_FakeClone, the input and output histogram have different binning!"<<endl;
2256         return;
2257     }
2258
2259     for (int i = 0; i < hist_in -> GetNbinsX(); i++){
2260         hist_out -> SetBinContent(i+1, hist_in -> GetBinContent(i+1));
2261     }
2262 }
2263
2264 // void INTTXYvtx::TH2F_FakeRebin(TH2F*hist_in, TH2F*hist_out)
2265 // {
2266
2267 // }
2268
2269 int INTTXYvtx::find_quadrant(pair<double,double> origin, pair<double,double> check_point)
2270 {
2271     int check_x, check_y;
2272
2273     if (origin.first < check_point.first){ check_x = 1; }
2274     else if (origin.first == check_point.first) { check_x = 0; }
2275     else { check_x = -1; }
2276
2277     if (origin.second < check_point.second){ check_y = 1; }
2278     else if (origin.second == check_point.second){ check_y = 0; }
2279     else { check_y = -1; }
2280
2281     if (check_x != 0 && check_y != 0) // note : normal 4 quadrants
2282     {
2283         return axis_quadrant_map[{check_x,check_y}];
2284     }
2285     else if (check_x == 0 && check_y != 0) // note : on the Y axis
2286     {
2287         return (check_y > 0) ? 4 : 5; // note : 4 is up, 5 is down, along the Y axis
2288     }
2289     else if (check_x != 0 && check_y == 0) // note : on the X axis
2290     {
2291         return (check_x > 0) ? 6 : 7; // note : 6 is right, 7 is left, along the X axis
2292     }
2293     else // note : the two position are identical!
2294     {
2295         return 8;
2296     }
2297 }
2298
2299 void INTTXYvtx::TH2FSampleLineFill(TH2F * hist_in, double segmentation, std::pair<double,double> inner_clu, std::pair<double,double> outer_clu)
2300 {
2301     double x_min = hist_in -> GetXaxis() -> GetXmin();
2302     double x_max = hist_in -> GetXaxis() -> GetXmax();
2303     double y_min = hist_in -> GetYaxis() -> GetXmin();
2304     double y_max = hist_in -> GetYaxis() -> GetXmax();
2305
2306     double seg_x, seg_y;
2307     double angle;
2308     int n_seg = 0;
2309
2310     while (true)
2311     {
2312         angle = atan2(inner_clu.second-outer_clu.second, inner_clu.first-outer_clu.first);
2313         seg_x = (n_seg * segmentation) * cos(angle) + outer_clu.first; // note : atan2(y,x), point.first is the radius
2314         seg_y = (n_seg * segmentation) * sin(angle) + outer_clu.second;
2315
2316         if ( (seg_x > x_min && seg_x < x_max && seg_y > y_min && seg_y < y_max) != true ) {break;}
2317         hist_in -> Fill(seg_x, seg_y);
2318         n_seg += 1;
2319     }
2320
2321     n_seg = 1;
2322     while (true)
2323     {
2324         angle = atan2(inner_clu.second-outer_clu.second, inner_clu.first-outer_clu.first);
2325         seg_x = (-1 * n_seg * segmentation) * cos(angle) + outer_clu.first; // note : atan2(y,x), point.first is the radius
2326         seg_y = (-1 * n_seg * segmentation) * sin(angle) + outer_clu.second;
2327
2328         if ( (seg_x > x_min && seg_x < x_max && seg_y > y_min && seg_y < y_max) != true ) {break;}
2329         hist_in -> Fill(seg_x, seg_y);
2330         n_seg += 1;
2331     }
2332 }
2333
2334 // note : the input is still mm
2335 vector<pair<double,double>> INTTXYvtx::FillLine_FindVertex(pair<double,double> window_center, double segmentation, double window_width, int N_bins, bool draw_plot)
2336 {
2337     // note : set the histogram to be cm
2338     xy_hist = new TH2F(
2339         "",
2340         "xy_hist",
2341         N_bins,
2342         (-1 * window_width / 2. + window_center.first) * 0.1,
2343         (window_width / 2. + window_center.first) * 0.1,
2344         N_bins,
2345         (-1 * window_width / 2. + window_center.second) * 0.1,
2346         (window_width / 2. + window_center.second) * 0.1
2347     );
2348     xy_hist -> GetXaxis() -> SetTitle("X axis [cm]");
2349     xy_hist -> GetYaxis() -> SetTitle("Y axis [cm]");
2350
2351     xy_hist_bkgrm = new TH2F(
2352         "",
2353         "xy_hist_bkgrm",
2354         N_bins,
2355         (-1 * window_width / 2. + window_center.first) * 0.1,
2356         (window_width / 2. + window_center.first) * 0.1,
2357         N_bins,
2358         (-1 * window_width / 2. + window_center.second) * 0.1,
2359         (window_width / 2. + window_center.second) * 0.1
2360     );
2361     xy_hist_bkgrm -> GetXaxis() -> SetTitle("X axis [cm]");
2362     xy_hist_bkgrm -> GetYaxis() -> SetTitle("Y axis [cm]");
2363
2364     // cout<<"test test size and bin of the hist xy_hist : "<<xy_hist -> GetNbinsX()<<" "<<xy_hist -> GetNbinsY()<<endl;
2365     // cout<<"test test bin width of the hist xy_hist : "<<xy_hist -> GetXaxis() -> GetBinWidth(1)<<" "<<xy_hist -> GetYaxis() -> GetBinWidth(1)<<endl;
2366     // cout<<"draw_plot status : "<<draw_plot<<endl;
2367
2368
2369     for (int i = 0; i < cluster_pair_vec.size(); i++)
2370     {
2371         vector<double> DCA_info_vec = calculateDistanceAndClosestPoint(
2372             cluster_pair_vec[i].first.x, cluster_pair_vec[i].first.y,
2373             cluster_pair_vec[i].second.x, cluster_pair_vec[i].second.y,
2374             window_center.first, window_center.second
2375         );
2376
2377         double DCA_sign = calculateAngleBetweenVectors(
2378             cluster_pair_vec[i].second.x, cluster_pair_vec[i].second.y,
2379             cluster_pair_vec[i].first.x, cluster_pair_vec[i].first.y,
2380             window_center.first, window_center.second
2381         );
2382
2383         if (DCA_info_vec[0] != fabs(DCA_sign) && fabs( DCA_info_vec[0] - fabs(DCA_sign) ) > 0.1){
2384             cout<<"different DCA : "<<DCA_info_vec[0]<<" "<<DCA_sign<<" diff : "<<DCA_info_vec[0] - fabs(DCA_sign)<<endl;
2385         }
2386
2387         Clus_InnerPhi_Offset = (cluster_pair_vec[i].first.y - window_center.second < 0) ? atan2(cluster_pair_vec[i].first.y - window_center.second, cluster_pair_vec[i].first.x - window_center.first) * (180./TMath::Pi()) + 360 : atan2(cluster_pair_vec[i].first.y - window_center.second, cluster_pair_vec[i].first.x - window_center.first) * (180./TMath::Pi());
2388         Clus_OuterPhi_Offset = (cluster_pair_vec[i].second.y - window_center.second < 0) ? atan2(cluster_pair_vec[i].second.y - window_center.second, cluster_pair_vec[i].second.x - window_center.first) * (180./TMath::Pi()) + 360 : atan2(cluster_pair_vec[i].second.y - window_center.second, cluster_pair_vec[i].second.x - window_center.first) * (180./TMath::Pi());
2389
2390         if (fabs(Clus_InnerPhi_Offset - Clus_OuterPhi_Offset) < 5)
2391         {
2392             TH2FSampleLineFill(xy_hist, segmentation * 0.1, {(cluster_pair_vec[i].first.x) * 0.1, (cluster_pair_vec[i].first.y) * 0.1}, {(DCA_info_vec[1]) * 0.1, (DCA_info_vec[2]) * 0.1});
2393             // note : the DCA cut may be biased since the DCA was calculated with respect to the vertex calculated by the quadrant method
2394             // if (DCA_cut.first < DCA_sign && DCA_sign < DCA_cut.second)
2395             // {
2396             //     TH2FSampleLineFill(xy_hist, segmentation, {cluster_pair_vec[i].first.x, cluster_pair_vec[i].first.y}, {DCA_info_vec[1], DCA_info_vec[2]});
2397             // }
2398         }
2399
2400     }
2401
2402     // note : try to remove the background
2403     TH2F_FakeClone(xy_hist,xy_hist_bkgrm);
2404     TH2F_threshold_advanced_2(xy_hist_bkgrm, 0.7);
2405
2406     // note : the unit is cm here
2407     double reco_vtx_x     = xy_hist_bkgrm->GetMean(1); // + xy_hist_bkgrm -> GetXaxis() -> GetBinWidth(1) / 2.; // note : the TH2F calculate the GetMean based on the bin center, no need to apply additional offset
2408     double reco_vtx_y     = xy_hist_bkgrm->GetMean(2); // + xy_hist_bkgrm -> GetYaxis() -> GetBinWidth(1) / 2.; // note : the TH2F calculate the GetMean based on the bin center, no need to apply additional offset
2409     double reco_vtx_x_err = xy_hist_bkgrm->GetMeanError(1);
2410     double reco_vtx_y_err = xy_hist_bkgrm->GetMeanError(2);
2411     double reco_vtx_x_std = xy_hist_bkgrm->GetStdDev(1);
2412     double reco_vtx_y_std = xy_hist_bkgrm->GetStdDev(2);
2413
2414     // cout<<"test : in the line filled, the process is almost done"<<endl;
2415
2416     if (draw_plot)
2417     {
2418         TGraph * reco_vertex_gr = new TGraph();
2419         reco_vertex_gr -> SetMarkerStyle(50);
2420         reco_vertex_gr -> SetMarkerColor(2);
2421         reco_vertex_gr -> SetMarkerSize(1);
2422         reco_vertex_gr -> SetPoint(reco_vertex_gr -> GetN(), reco_vtx_x, reco_vtx_y); // note : change the unit to cm
2423
2424         // note : change the unit to cm, for the presentation
2425         // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
2426         xy_hist -> SetStats(0);
2427         // xy_hist -> GetXaxis() -> SetLimits( xy_hist -> GetXaxis() -> GetBinLowEdge(1) * 0.1, (xy_hist -> GetXaxis() -> GetBinLowEdge( xy_hist -> GetNbinsX() ) + xy_hist -> GetXaxis() -> GetBinWidth( xy_hist -> GetNbinsX() )) * 0.1 );
2428         // xy_hist -> GetYaxis() -> SetLimits( xy_hist -> GetYaxis() -> GetBinLowEdge(1) * 0.1, (xy_hist -> GetYaxis() -> GetBinLowEdge( xy_hist -> GetNbinsY() ) + xy_hist -> GetYaxis() -> GetBinWidth( xy_hist -> GetNbinsY() )) * 0.1 );
2429         xy_hist -> GetXaxis() -> SetNdivisions(505);
2430         xy_hist -> Draw("colz0");
2431         // draw_text -> DrawLatex(0.21, 0.71+0.13, Form("Vertex of the Run: %.3f mm, %.3f mm", reco_vtx_x, reco_vtx_y));
2432         // draw_text -> DrawLatex(0.21, 0.67+0.13, Form("Vertex error: %.3f mm, %.3f mm", xy_hist_bkgrm->GetMeanError(1), xy_hist_bkgrm->GetMeanError(2)));
2433         // reco_vertex_gr -> Draw("p same");
2434         ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
2435         c1 -> Print(Form("%s/Run_xy_hist.pdf",out_folder_directory.c_str()));
2436         c1 -> Clear();
2437
2438         // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
2439         xy_hist_bkgrm -> SetStats(0);
2440         // xy_hist_bkgrm -> GetYaxis() -> SetLimits( xy_hist_bkgrm -> GetYaxis() -> GetBinLowEdge(1) * 0.1, (xy_hist_bkgrm -> GetYaxis() -> GetBinLowEdge( xy_hist_bkgrm -> GetNbinsY() ) + xy_hist_bkgrm -> GetYaxis() -> GetBinWidth( xy_hist_bkgrm -> GetNbinsY() )) * 0.1 );
2441         // xy_hist_bkgrm -> GetXaxis() -> SetLimits( xy_hist_bkgrm -> GetXaxis() -> GetBinLowEdge(1) * 0.1, (xy_hist_bkgrm -> GetXaxis() -> GetBinLowEdge( xy_hist_bkgrm -> GetNbinsX() ) + xy_hist_bkgrm -> GetXaxis() -> GetBinWidth( xy_hist_bkgrm -> GetNbinsX() )) * 0.1 );
2442         xy_hist_bkgrm -> GetXaxis() -> SetNdivisions(505);
2443         xy_hist_bkgrm -> Draw("colz0");
2444         draw_text -> DrawLatex(0.21, 0.71+0.13, Form("Vertex of the Run: %.5f cm, %.5f cm", reco_vtx_x, reco_vtx_y));
2445         draw_text -> DrawLatex(0.21, 0.67+0.13, Form("Vertex error: %.5f cm, %.5f cm", reco_vtx_x_err, reco_vtx_y_err));
2446         reco_vertex_gr -> Draw("p same");
2447         ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
2448         c1 -> Print(Form("%s/Run_xy_hist_bkgrm.pdf",out_folder_directory.c_str()));
2449         c1 -> Clear();
2450
2451         // cout<<"test : hello, can you see me ?"<<endl;
2452     }
2453
2454     // note: change the unit from cm to mm
2455     return {{reco_vtx_x * 10.,reco_vtx_y * 10.}, {reco_vtx_x_err * 10.,reco_vtx_y_err * 10.}, {reco_vtx_x_std * 10.,reco_vtx_y_std * 10.}};
2456 }
2457
2458 // todo : the input histogram should be in the unit of cm
2459 vector<double> INTTXYvtx::Get_covariance_TH2(TH2F * hist_in)
2460 {
2461     double X_mean = hist_in -> GetMean(1);
2462     double Y_mean = hist_in -> GetMean(2);
2463
2464     double denominator = 0;
2465     double variance_x  = 0;
2466     double variance_y  = 0;
2467     double covariance = 0;
2468
2469     for (int xi = 0; xi < hist_in -> GetNbinsX(); xi++){
2470         for (int yi = 0; yi < hist_in -> GetNbinsY(); yi++){
2471             double cell_weight = hist_in -> GetBinContent(xi+1, yi+1);
2472             double cell_x = hist_in -> GetXaxis() -> GetBinCenter(xi+1);
2473             double cell_y = hist_in -> GetYaxis() -> GetBinCenter(yi+1);
2474
2475             denominator += pow(cell_weight, 2);
2476             variance_x  += pow(cell_x - X_mean, 2) * pow(cell_weight,2);
2477             variance_y  += pow(cell_y - Y_mean, 2) * pow(cell_weight,2);
2478             covariance  += (cell_x - X_mean) * (cell_y - Y_mean) * pow(cell_weight,2);
2479         }
2480     }
2481
2482     // note : the histogram is in the unit of cm, we change it back to mm
2483     return {X_mean * 10., Y_mean * 10., (variance_x/denominator) * 10., (variance_y/denominator) * 10., (covariance/denominator) * 10.};
2484 }
2485
2486 //note : accumulate the number of entries from both sides of the histogram
2487 double INTTXYvtx::get_dist_offset(TH1F * hist_in, int check_N_bin) // note : check_N_bin 1 to N bins of hist
2488 {
2489     if (check_N_bin < 0 || check_N_bin > hist_in -> GetNbinsX()) {cout<<" wrong check_N_bin "<<endl; exit(1);}
2490     double total_entry = 0;
2491     for (int i = 0; i < check_N_bin; i++)
2492     {
2493         total_entry += hist_in -> GetBinContent(i+1); // note : 1, 2, 3.....
2494         total_entry += hist_in -> GetBinContent(hist_in -> GetNbinsX() - i);
2495     }
2496
2497     return total_entry/double(2. * check_N_bin);
2498 }
2499
2500 #endif
2501
2502
2503     // //note : N_trial has to be odd -> even_number + 0 + event_number
2504     // for (int i = 0; i < N_trial; i++)
2505     // {
2506
2507     //     sub_out_folder_name = Form("%s/New_trial_%i",out_folder_directory.c_str(), i);
2508     //     system(Form("mkdir %s",sub_out_folder_name.c_str()));
2509     //     tested_angle = correction_circle.second + (i - ((N_trial - 1)/2) ) * moving_unit - 90.;
2510
2511     //     // cout<<"test : "<<correction_circle.first<<" "<<beam_origin.first<<" "<<beam_origin.second<<endl;
2512     //     // cout<<"test : "<<correction.x<<" "<<correction.y<<endl;
2513
2514     //     correction = PolarToCartesian(correction_circle.first, tested_angle);
2515     //     current_vtxX = beam_origin.first + correction.x;
2516     //     current_vtxY = beam_origin.second + correction.y;
2517
2518     //     info_vec = GetVTXxyCorrection_new(i);
2519     //     cout<<"given angle: "<<tested_angle<<endl;
2520     //     cout<<"Pol1 fit, chi2: "<<horizontal_fit_inner->GetChisquare()<<" NDF :"<<horizontal_fit_inner->GetNDF()<<" fit error: #pm"<<horizontal_fit_inner->GetParError(0)<<endl;
2521
2522     //     Xpos_vec.push_back(tested_angle);
2523     //     XE_vec.push_back(0);
2524     //     Ypos_vec_stddev.push_back(info_vec[0]);
2525     //     YE_vec_stddev.push_back(info_vec[1]);
2526     //     Ypos_vec_pol0_chi2.push_back(info_vec[2]);
2527     //     Ypos_vec_pol0_fitE.push_back(info_vec[3]);
2528     //     YE_vec_pol0.push_back(0);
2529
2530     //     PrintPlotsVTXxy(sub_out_folder_name, 1);
2531     //     ClearHist(1);
2532     // }
2533
2534     // DrawTGraphErrors(Xpos_vec, Ypos_vec_stddev, XE_vec, YE_vec_stddev, out_folder_directory, {Form("angle_diff_pos_relation"), "Correction angle [degree]", "Angle diff StdDev [degree]"});
2535     // DrawTGraphErrors(Xpos_vec, Ypos_vec_pol0_chi2, XE_vec, YE_vec_pol0, out_folder_directory, {Form("Pol0_chi2_pos_relation"), "Correction angle [degree]", "Pol0 fit reduced #chi^{2}"});
2536     // DrawTGraphErrors(Xpos_vec, Ypos_vec_pol0_fitE, XE_vec, YE_vec_pol0, out_folder_directory, {Form("Pol0_fitE_pos_relation"), "Correction angle [degree]", "Pol0 fit fit error [mm]"});