analysis_INTT_CW/DST_MC/INTTZvtx.h

0001 #ifndef INTTZvtx_h
0002 #define INTTZvtx_h
0003
0004 #include "../private_cluster_gen/InttConversion_new.h"
0005 #include "../private_cluster_gen/InttClustering.h"
0006 #include "../sigmaEff.h"
0007 #include "gaus_func.h"
0008 #include "MegaTrackFinder.h"
0009 #include "../vector_stddev.h"
0010 // #include "../sPhenixStyle.C"
0011
0012 // note : the centrality is still float here.
0013
0014 double double_gaus_func(double *x, double *par)
0015 {
0016     // note : par[0] : size of gaus1
0017     // note : par[1] : mean of gaus1
0018     // note : par[2] : width of gaus1
0019
0020     // note : par[3] : size of gaus2
0021     // note : par[4] : mean of gaus2
0022     // note : par[5] : width of gaus2
0023
0024     // note : par[6] : systematic ffset
0025     double gaus1 = par[0] * TMath::Gaus(x[0],par[1],par[2]);
0026     double gaus2 = par[3] * TMath::Gaus(x[0],par[4],par[5]);
0027
0028     return gaus1 + gaus2 + par[6];
0029 }
0030
0031 class INTTZvtx
0032 {
0033     public :
0034         INTTZvtx(
0035             string run_type,
0036             string out_folder_directory,
0037             pair<double,double> beam_origin,
0038             int geo_mode_id,
0039             double phi_diff_cut = 0.11,
0040             pair<double, double> DCA_cut = {-1,1},
0041             int N_clu_cutl = 20,
0042             int N_clu_cut = 10000,
0043             int zvtx_cal_require = 15,
0044             pair<double,double> zvtx_QA_width = {39.62, 65.36},
0045             double zvtx_QA_ratio = 0.00001,
0046             bool draw_event_display = true,
0047             double peek = 3.32405,
0048             bool print_message_opt = true,
0049             TH1F * delta_phi_inner_correction_hist_in = nullptr // note : special_tag
0050         );
0051         void ProcessEvt(
0052             int event_i,
0053             vector<clu_info> temp_sPH_inner_nocolumn_vec,
0054             vector<clu_info> temp_sPH_outer_nocolumn_vec,
0055             vector<vector<double>> temp_sPH_nocolumn_vec,
0056             vector<vector<double>> temp_sPH_nocolumn_rz_vec,
0057             int NvtxMC,
0058             double TrigZvtxMC,
0059             bool PhiCheckTag,
0060             Long64_t bco_full,
0061             float centrality_float,
0062             double MBD_reco_z,
0063             int is_min_bias,
0064             int is_min_bias_wozdc,
0065             double MBD_north_charge_sum,
0066             double MBD_south_charge_sum
0067         );
0068         void ClearEvt();
0069         void PrintPlots();
0070         void EndRun();
0071
0072         double GetZdiffPeakMC();
0073         double GetZdiffWidthMC();
0074
0075         vector<double> GetEvtZPeak();
0076
0077
0078     private :
0079         TCanvas * c2;
0080         TCanvas * c1;
0081         TPad * pad_xy;
0082         TPad * pad_rz;
0083         TPad * pad_z;
0084         TPad * pad_z_hist;
0085         TPad * pad_z_line;
0086         TPad * pad_phi_diff;
0087         TPad * pad_track_phi;
0088         TPad * pad_inner_outer_phi;
0089         TPad * pad_phi_diff_1D;
0090         TPad * pad_EvtZDist;
0091         TPad * pad_ZoomIn_EvtZDist;
0092
0093         TH1F * avg_event_zvtx;
0094         TH1F * zvtx_evt_fitError;
0095         TH2F * zvtx_evt_fitError_corre;
0096         TH2F * zvtx_evt_width_corre;
0097         TH2F * zvtx_evt_nclu_corre;
0098         TH1F * width_density;         // note : N good hits / width
0099         TH1F * ES_width;
0100         TH1F * ES_width_ratio;
0101         TH2F * Z_resolution_Nclu;
0102         TH2F * Z_resolution_pos;
0103         TH2F * Z_resolution_pos_cut;
0104         TH1F * Z_resolution;
0105         TH1F * evt_possible_z;
0106         TH1F * line_breakdown_hist_cm;            // note : same thing, but with the unit cm
0107         TH1F * line_breakdown_hist_cm_zoomin;     // note : same thing, but with the unit cm
0108         TH1F * line_breakdown_hist;               // note : try to fill the line into the histogram
0109         TH1F * line_breakdown_hist_zoomin;        // note : try to fill the line into the histogram
0110         TH1F * line_breakdown_gaus_ratio_hist;    // note : the distribution of the entry/width of gaus fit
0111         TH1F * line_breakdown_gaus_width_hist;    // note : the distribution of the gaus fit width
0112         TH2F * gaus_width_Nclu;
0113         TH2F * gaus_rchi2_Nclu;
0114         TH2F * final_fit_width;
0115         TH2F * N_track_candidate_Nclu; // note : Number of tracklet candidate (In xy plane) vs number of clusters
0116         TH1F * temp_event_zvtx;
0117         TH1F * peak_group_width_hist;
0118         TH1F * peak_group_ratio_hist;
0119         TH1F * N_group_hist;
0120         TH1F * peak_group_detail_width_hist;
0121         TH1F * peak_group_detail_ratio_hist;
0122         TH1F * N_group_detail_hist;
0123         TH1F * evt_select_track_phi;
0124         TH1F * evt_phi_diff_1D;
0125         TH2F * evt_phi_diff_inner_phi;
0126         TH2F * evt_inner_outer_phi;
0127
0128
0129
0130         InttConversion * ch_pos_DB;
0131         TLatex * ltx;
0132         TLatex * draw_text;
0133         TLine  * ladder_line;
0134         TLine  * eff_sig_range_line;
0135         TLine  * final_fit_range_line;
0136         TLine  * track_line;
0137         TLine  * coord_line;
0138         TFile  * out_file;
0139         TTree  * tree_out;
0140         vector<TF1 *> gaus_fit_vec; // note : mm
0141         vector<TF1 *> gaus_fit_cm_vec; // note : cm
0142         TF1    * gaus_fit;
0143         TF1    * gaus_fit_cm;
0144         TF1    * gaus_fit_width_cm;
0145         TF1    * gaus_fit_2;
0146         TF1    * double_gaus_fit;
0147         TF1    * zvtx_finder;
0148
0149         pair<double,double> evt_possible_z_range = {-700, 700};
0150         vector<string> conversion_mode_BD = {"ideal", "survey_1_XYAlpha_Peek", "full_survey_3.32"};
0151         double Integrate_portion_final = 0.68;
0152         double Integrate_portion = 0.35; // todo : Width selection per event (the range finder for fitting)
0153         double width_density_cut = 0.00055; // Todo : change the width density here
0154         double high_multi_line = 1000; // todo : the cut to classify the high-low multiplicity, which are fit with different method
0155         double zvtx_hist_l = -500;
0156         double zvtx_hist_r = 500;
0157         int print_rate = 50; // todo : the print rate is here
0158
0159         pair<double, double> beam_origin;
0160         pair<double, double> DCA_cut;       // note : if (< DCA_cut)          -> pass      unit mm
0161         pair<double, double> zvtx_QA_width; // note : for the zvtx range Quality check, check the width
0162         string out_folder_directory;
0163         string run_type;
0164         double peek;
0165         double phi_diff_cut;          // note : if (< phi_diff_cut)      -> pass      unit degree
0166         double zvtx_QA_ratio;         // note : for the zvtx range QUality check, check the Norm. size / width
0167         double width_density_par;
0168         double Clus_InnerPhi_Offset; // note : the vertex in XY is not at zero, so the "offset" moves the offset back to the orign which is (0,0)
0169         double Clus_OuterPhi_Offset; // note : the vertex in XY is not at zero, so the "offset" moves the offset back to the orign which is (0,0)
0170         bool draw_event_display;
0171         bool print_message_opt;
0172         int zvtx_cal_require;         // note : if (> zvtx_cal_require)  -> pass
0173         int geo_mode_id;
0174         int N_clu_cut;                // note : if (> N_clu_cut)         -> continue  unit number
0175         int N_clu_cutl;               // note : if (< N_clu_cutl)        -> continue  unit number
0176
0177         TH1F * delta_phi_inner_correction_hist; // note : special_tag
0178
0179         vector<vector<pair<bool,clu_info>>> inner_clu_phi_map; // note: phi
0180         vector<vector<pair<bool,clu_info>>> outer_clu_phi_map; // note: phi
0181
0182         // note : for tree_out
0183         double out_ES_zvtx, out_ES_zvtxE, out_ES_rangeL, out_ES_rangeR, out_ES_width_density, MC_true_zvtx;
0184         double out_LB_Gaus_Mean_mean, out_LB_Gaus_Mean_meanE, out_LB_Gaus_Mean_width, out_LB_Gaus_Mean_chi2;
0185         double out_LB_Gaus_Width_width, out_LB_Gaus_Width_size_width, out_LB_Gaus_Width_offset, out_LB_geo_mean;
0186         double out_LB_refit, out_LB_refitE; // note : not used in the out tree
0187         double out_mid_cut_peak_width, out_mid_cut_peak_ratio, out_LB_cut_peak_width, out_LB_cut_peak_ratio;
0188         double out_MBD_reco_z;
0189         int    out_is_min_bias, out_is_min_bias_wozdc;
0190         double out_MBD_south_charge_sum, out_MBD_north_charge_sum;
0191         float out_centrality_float;
0192         bool out_good_zvtx_tag;
0193         int out_eID, N_cluster_outer_out, N_cluster_inner_out, out_ES_N_good, out_mid_cut_Ngroup, out_LB_cut_Ngroup;
0194         int out_N_cluster_north, out_N_cluster_south;
0195         vector<double> out_evtZ_hist_cut_content_vec; // note : this is for the training test
0196         Long64_t bco_full_out;
0197
0198         // note : for out parameters
0199         double MC_z_diff_peak, MC_z_diff_width; // note : the comparison between Reco - true in MC. Values are from fitting foucsing on the peak region.
0200
0201         void Init();
0202         void InitHist();
0203         void InitCanvas();
0204         void InitTreeOut();
0205         void InitRest();
0206
0207         vector<vector<double>> good_track_xy_vec; // note : for event drawing, not used now
0208         vector<vector<double>> good_track_rz_vec; // note : for event drawing, not used now
0209         vector<vector<double>> good_comb_rz_vec;  // note : not used
0210         vector<vector<double>> good_comb_xy_vec;  // note : not used
0211         vector<vector<double>> good_comb_phi_vec; // note : not used
0212         vector<vector<double>> good_tracklet_rz;  // note : not used
0213         vector<float> temp_event_zvtx_vec;        // note : not used
0214         vector<float> temp_event_zvtx_info;
0215         vector<float> avg_event_zvtx_vec;
0216         vector<float> Z_resolution_vec;
0217         vector<float> line_breakdown_vec;
0218         vector<double> N_group_info; // note : the information of groups remaining in the histogram after the strong background suppression
0219         vector<double> N_group_info_detail; // note : detail
0220         string plot_text;
0221         double final_selection_widthD;
0222         double final_selection_widthU;
0223         double gaus_ratio;
0224         double gaus_fit_offset; // note : when doing the integral, removing the offset is more preferable, so keep the offset here
0225         double final_zvtx;
0226         double tight_offset_width;
0227         double tight_offset_peak;
0228         double loose_offset_peak;
0229         double loose_offset_peakE;
0230         bool good_zvtx_tag;
0231         double good_zvtx_tag_int;
0232         long total_NClus;
0233         int fit_winner_id;
0234         int good_comb_id;
0235
0236
0237         // note : for the trapezoidal method for the vertex determination
0238         TH1F * combination_zvtx_range_shape;
0239         TRandom3 * rand_gen;
0240         TH2F * best_fitwidth_NClus;
0241         vector<double> fit_mean_mean_vec;
0242         vector<double> fit_mean_reducedChi2_vec;
0243         vector<double> fit_mean_width_vec;
0244         std::vector<std::string> color_code = {
0245             "#9e0142",
0246             "#d53e4f",
0247             "#f46d43",
0248             "#fdae61",
0249             "#fee08b",
0250             "#e6f598",
0251             "#abdda4",
0252             "#66c2a5",
0253             "#3288bd",
0254             "#5e4fa2"
0255         };
0256
0257         // note : the class that handles the 3/4-cluster tracklet finder
0258         MegaTrackFinder * mega_track_finder;
0259
0260         TGraphErrors * z_range_gr;
0261         TGraphErrors * z_range_gr_draw;
0262         TGraph * temp_event_xy;
0263         TGraph * temp_event_rz;
0264         TGraph * bkg;
0265
0266
0267         // note : in the event process
0268         vector<float> N_comb; vector<float> N_comb_e; vector<float> z_mid; vector<float> z_range; vector<double> N_comb_phi;
0269         vector<double> eff_N_comb; vector<double> eff_z_mid; vector<double> eff_N_comb_e; vector<double> eff_z_range; // note : eff_sig
0270         vector<double> LB_N_comb; vector<double> LB_z_mid; vector<double> LB_N_comb_e; vector<double> LB_z_range; // note : LB gaus
0271
0272         std::vector<double> calculateDistanceAndClosestPoint(double x1, double y1, double x2, double y2, double target_x, double target_y);
0273         pair<double,double> Get_possible_zvtx(double rvtx, vector<double> p0, vector<double> p1);
0274         vector<double> find_Ngroup(TH1F * hist_in);
0275         vector<double> get_half_hist_vec(TH1F * hist_in);
0276         double get_radius(double x, double y);
0277         double calculateAngleBetweenVectors(double x1, double y1, double x2, double y2, double targetX, double targetY);
0278         double Get_extrapolation(double given_y, double p0x, double p0y, double p1x, double p1y);
0279         double LB_geo_mean(TH1F * hist_in, pair<double, double> search_range, int event_i);
0280         void temp_bkg(TPad * c1, double peek, pair<double,double> beam_origin, InttConversion * ch_pos_DB);
0281         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);
0282         void line_breakdown(TH1F* hist_in, pair<double,double> line_range);
0283         void trapezoidal_line_breakdown(TH1F * hist_in, double inner_r, double inner_z, double outer_r, double outer_z);
0284         int min_width_ID(vector<double> input_width_vec);
0285         double get_delta_phi(double angle_1, double angle_2);
0286         double get_track_phi(double inner_clu_phi_in, double delta_phi_in);
0287         void FakeClone1D(TH1F * hist_in, TH1F * hist_out);
0288
0289
0290 };
0291
0292
0293
0294 INTTZvtx::INTTZvtx(
0295     string run_type,
0296     string out_folder_directory,
0297     pair<double,double> beam_origin,
0298     int geo_mode_id,
0299     double phi_diff_cut,
0300     pair<double,double> DCA_cut,
0301     int N_clu_cutl,
0302     int N_clu_cut,
0303     int zvtx_cal_require,
0304     pair<double,double> zvtx_QA_width,
0305     double zvtx_QA_ratio,
0306     bool draw_event_display,
0307     double peek,
0308     bool print_message_opt,
0309     TH1F * delta_phi_inner_correction_hist_in // note : special_tag
0310 )
0311 :run_type(run_type),
0312 out_folder_directory(out_folder_directory),
0313 beam_origin(beam_origin),
0314 geo_mode_id(geo_mode_id),
0315 peek(peek),
0316 N_clu_cut(N_clu_cut),
0317 N_clu_cutl(N_clu_cutl),
0318 phi_diff_cut(phi_diff_cut),
0319 DCA_cut(DCA_cut),
0320 zvtx_cal_require(zvtx_cal_require),
0321 zvtx_QA_width(zvtx_QA_width),
0322 zvtx_QA_ratio(zvtx_QA_ratio),
0323 draw_event_display(draw_event_display),
0324 print_message_opt(print_message_opt),
0325 delta_phi_inner_correction_hist(delta_phi_inner_correction_hist_in) // note : special_tag
0326 {
0327     SetsPhenixStyle();
0328     system(Form("mkdir %s",out_folder_directory.c_str()));
0329     gErrorIgnoreLevel = kWarning; // note : To not print the "print plot info."
0330
0331     fit_mean_mean_vec.clear();
0332     fit_mean_reducedChi2_vec.clear();
0333     fit_mean_width_vec.clear();
0334
0335     out_evtZ_hist_cut_content_vec.clear();
0336
0337     temp_event_zvtx_vec.clear();
0338     temp_event_zvtx_info.clear();
0339     good_track_xy_vec.clear();
0340     good_track_rz_vec.clear();
0341     good_comb_rz_vec.clear();
0342     good_comb_xy_vec.clear();
0343     good_comb_phi_vec.clear();
0344     good_tracklet_rz.clear();
0345     avg_event_zvtx_vec.clear();
0346     Z_resolution_vec.clear();
0347     line_breakdown_vec.clear();
0348     good_comb_id = 0;
0349     MC_z_diff_peak = -777.;
0350     MC_z_diff_width = -777.;
0351
0352     out_N_cluster_south = 0;
0353     out_N_cluster_north = 0;
0354
0355     N_group_info_detail = {-1.,-1.,-1.,-1.};
0356
0357     N_comb.clear(); N_comb_e.clear(); z_mid.clear(); z_range.clear(); N_comb_phi.clear();
0358     eff_N_comb.clear(); eff_z_mid.clear(); eff_N_comb_e.clear(); eff_z_range.clear(); // note : eff_sig
0359     LB_N_comb.clear(); LB_z_mid.clear(); LB_N_comb_e.clear(); LB_z_range.clear(); // note : LB gaus
0360
0361     inner_clu_phi_map.clear();
0362     outer_clu_phi_map.clear();
0363     inner_clu_phi_map = vector<vector<pair<bool,clu_info>>>(360);
0364     outer_clu_phi_map = vector<vector<pair<bool,clu_info>>>(360);
0365
0366     Init();
0367
0368     plot_text = (run_type == "MC") ? "Simulation" : "Work-in-progress";
0369
0370     // mega_track_finder = new MegaTrackFinder(run_type, out_folder_directory, centrality_map.size(), beam_origin);
0371
0372     if (draw_event_display)
0373     {
0374         c2 -> cd();
0375         c2 -> Print(Form("%s/temp_event_display.pdf(",out_folder_directory.c_str()));
0376         c1 -> cd();
0377         c1 -> Print(Form("%s/evt_linebreak_display.pdf(",out_folder_directory.c_str()));
0378     }
0379 }
0380
0381 void INTTZvtx::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)
0382 {
0383     if (setgrid_bool == true) {pad -> SetGrid (1, 1);}
0384     pad -> SetLeftMargin   (left);
0385     pad -> SetRightMargin  (right);
0386     pad -> SetTopMargin    (top);
0387     pad -> SetBottomMargin (bottom);
0388     pad -> SetTicks(1,1);
0389     if (set_logY == true)
0390     {
0391         pad -> SetLogy (1);
0392     }
0393
0394 }
0395
0396 void INTTZvtx::Init()
0397 {
0398     InitCanvas();
0399     InitHist();
0400     InitTreeOut();
0401     InitRest();
0402 }
0403
0404 void INTTZvtx::InitHist()
0405 {
0406     avg_event_zvtx = new TH1F("","avg_event_zvtx",100,zvtx_hist_l,zvtx_hist_r);
0407     // avg_event_zvtx -> SetMarkerStyle(20);
0408     // avg_event_zvtx -> SetMarkerSize(0.8);
0409     // avg_event_zvtx -> SetMarkerColor(TColor::GetColor("#1A3947"));
0410     avg_event_zvtx -> SetLineColor(TColor::GetColor("#1A3947"));
0411     avg_event_zvtx -> SetLineWidth(2);
0412     avg_event_zvtx -> GetYaxis() -> SetTitle("Entry");
0413     avg_event_zvtx -> GetXaxis() -> SetTitle("Z vertex position [mm]");
0414     avg_event_zvtx -> GetYaxis() -> SetRangeUser(0,50);
0415     avg_event_zvtx -> SetTitleSize(0.06, "X");
0416     avg_event_zvtx -> SetTitleSize(0.06, "Y");
0417     avg_event_zvtx -> GetXaxis() -> SetTitleOffset(0.82);
0418     avg_event_zvtx -> GetYaxis() -> SetTitleOffset(1.1);
0419     avg_event_zvtx -> GetXaxis() -> CenterTitle(true);
0420     avg_event_zvtx -> GetYaxis() -> CenterTitle(true);
0421     avg_event_zvtx -> GetXaxis() -> SetNdivisions(505);
0422
0423     zvtx_evt_fitError = new TH1F("","zvtx_evt_fitError",150,0,50);
0424     zvtx_evt_fitError -> GetXaxis() -> SetTitle(" mm ");
0425     zvtx_evt_fitError -> GetYaxis() -> SetTitle("entry");
0426     zvtx_evt_fitError -> GetXaxis() -> SetNdivisions(505);
0427
0428     zvtx_evt_fitError_corre = new TH2F("","zvtx_evt_fitError_corre",200,0,10000,200,0,20);
0429     zvtx_evt_fitError_corre -> GetXaxis() -> SetTitle(" # of clusters ");
0430     zvtx_evt_fitError_corre -> GetYaxis() -> SetTitle(" #pm mm ");
0431     zvtx_evt_fitError_corre -> GetXaxis() -> SetNdivisions(505);
0432
0433     zvtx_evt_width_corre = new TH2F("","zvtx_evt_width_corre",200,0,10000,200,0,300);
0434     zvtx_evt_width_corre -> GetXaxis() -> SetTitle(" # of clusters ");
0435     zvtx_evt_width_corre -> GetYaxis() -> SetTitle(" mm ");
0436     zvtx_evt_width_corre -> GetXaxis() -> SetNdivisions(505);
0437
0438     zvtx_evt_nclu_corre = new TH2F("","zvtx_evt_nclu_corre",200,0,10000,200,-1000,1000);
0439     zvtx_evt_nclu_corre -> GetXaxis() -> SetTitle(" # of clusters ");
0440     zvtx_evt_nclu_corre -> GetYaxis() -> SetTitle(" zvtx [mm] ");
0441     zvtx_evt_nclu_corre -> GetXaxis() -> SetNdivisions(505);
0442
0443     width_density = new TH1F("","width_density",200,0,0.006); // note : N good hits / width
0444     width_density -> GetXaxis() -> SetTitle(" (N good track / width) / NClus ");
0445     width_density -> GetYaxis() -> SetTitle(" Entry ");
0446     width_density -> GetXaxis() -> SetNdivisions(505);
0447
0448     ES_width = new TH1F("","ES_width",200,0,150);
0449     ES_width -> GetXaxis() -> SetTitle(" Width [mm] ");
0450     ES_width -> GetYaxis() -> SetTitle(" Entry ");
0451     ES_width -> GetXaxis() -> SetNdivisions(505);
0452
0453     ES_width_ratio = new TH1F("","ES_width_ratio",200,0,60);
0454     ES_width_ratio -> GetXaxis() -> SetTitle(" NClus / Width ");
0455     ES_width_ratio -> GetYaxis() -> SetTitle(" Entry ");
0456     ES_width_ratio -> GetXaxis() -> SetNdivisions(505);
0457
0458     Z_resolution_Nclu = new TH2F("","Z_resolution_Nclu",200,0,10000,200,-100,100);
0459     Z_resolution_Nclu -> GetXaxis() -> SetTitle(" # of clusters ");
0460     Z_resolution_Nclu -> GetYaxis() -> SetTitle("#Delta Z (Reco - True) [mm]");
0461     Z_resolution_Nclu -> GetXaxis() -> SetNdivisions(505);
0462
0463     Z_resolution_pos = new TH2F("","Z_resolution_pos",200,-350,350,200,-50,50);
0464     Z_resolution_pos -> GetXaxis() -> SetTitle("True Zvtx [mm]");
0465     Z_resolution_pos -> GetYaxis() -> SetTitle("#Delta Z (Reco - True) [mm]");
0466     Z_resolution_pos -> GetXaxis() -> SetNdivisions(505);
0467
0468     Z_resolution_pos_cut = new TH2F("","Z_resolution_pos_cut",200,-350,350,200,-50,50);
0469     Z_resolution_pos_cut -> GetXaxis() -> SetTitle("True Zvtx [mm]");
0470     Z_resolution_pos_cut -> GetYaxis() -> SetTitle("#Delta Z (Reco - True) [mm]");
0471     Z_resolution_pos_cut -> GetXaxis() -> SetNdivisions(505);
0472
0473     Z_resolution = new TH1F("","Z_resolution",200,-30,30);
0474     Z_resolution -> GetXaxis() -> SetTitle("#Delta Z (Reco - True) [mm]");
0475     Z_resolution -> GetYaxis() -> SetTitle(" Entry ");
0476     Z_resolution -> GetXaxis() -> SetNdivisions(505);
0477
0478     evt_possible_z = new TH1F("","evt_possible_z",50, evt_possible_z_range.first, evt_possible_z_range.second);
0479     evt_possible_z -> SetLineWidth(1);
0480     evt_possible_z -> GetXaxis() -> SetTitle("Z [mm]");
0481     evt_possible_z -> GetYaxis() -> SetTitle("Entry");
0482
0483     int N = 1200;        // note : N bins for each side, regardless the bin at zero
0484     double width = 0.5;  // note : bin width with the unit [mm]
0485     line_breakdown_hist = new TH1F("", "line_breakdown_hist", 2*N+1, -1*(width*N + width/2.), width*N + width/2.);
0486     line_breakdown_hist -> SetLineWidth(1);
0487     line_breakdown_hist -> GetXaxis() -> SetTitle("Z [mm]");
0488     line_breakdown_hist -> GetYaxis() -> SetTitle("Entries [A.U.]");
0489     line_breakdown_hist -> GetXaxis() -> SetNdivisions(705);
0490
0491     line_breakdown_hist_zoomin = new TH1F(
0492         "",
0493         "line_breakdown_hist_zoomin;Z [mm];Entries [A.U.]",
0494         line_breakdown_hist -> GetNbinsX(),
0495         line_breakdown_hist -> GetXaxis() -> GetXmin(),
0496         line_breakdown_hist -> GetXaxis() -> GetXmax()
0497     );
0498     line_breakdown_hist_zoomin -> GetXaxis() -> SetNdivisions(705);
0499
0500     if (print_message_opt == true) {cout<<"class INTTZvtx, Line brakdown hist, range : "<<line_breakdown_hist->GetXaxis()->GetXmin()<<" "<<line_breakdown_hist->GetXaxis()->GetXmax()<<" "<<line_breakdown_hist->GetBinWidth(1)<<endl;}
0501
0502     line_breakdown_hist_cm = new TH1F(
0503         "",
0504         "line_breakdown_hist_cm;Z [cm];Entries [A.U.]",
0505         line_breakdown_hist -> GetNbinsX(),
0506         line_breakdown_hist -> GetXaxis() -> GetXmin() * 0.1,
0507         line_breakdown_hist -> GetXaxis() -> GetXmax() * 0.1
0508     );
0509     line_breakdown_hist_cm -> SetLineWidth(1);
0510     line_breakdown_hist_cm -> GetXaxis() -> SetNdivisions(705);
0511
0512     line_breakdown_hist_cm_zoomin = new TH1F(
0513         "",
0514         "line_breakdown_hist_cm_zoomin;Z [cm];Entries [A.U.]",
0515         line_breakdown_hist -> GetNbinsX(),
0516         line_breakdown_hist -> GetXaxis() -> GetXmin() * 0.1,
0517         line_breakdown_hist -> GetXaxis() -> GetXmax() * 0.1
0518     );
0519     line_breakdown_hist_cm_zoomin -> GetXaxis() -> SetNdivisions(705);
0520
0521     // note : unit : mm
0522     combination_zvtx_range_shape = new TH1F(
0523         "",
0524         "",
0525         line_breakdown_hist -> GetNbinsX(),
0526         line_breakdown_hist -> GetXaxis() -> GetXmin(),
0527         line_breakdown_hist -> GetXaxis() -> GetXmax()
0528     );
0529
0530     best_fitwidth_NClus = new TH2F(
0531         "",
0532         "best_fitwidth_NClus;N proto-tracklet;Best fit width [mm]",
0533         200,0,10000,
0534         100,0,200
0535     );
0536
0537     line_breakdown_gaus_ratio_hist = new TH1F("","line_breakdown_gaus_ratio_hist",200,0,0.0005);
0538     line_breakdown_gaus_ratio_hist -> GetXaxis() -> SetTitle("(Norm. size) / width");
0539     line_breakdown_gaus_ratio_hist -> GetYaxis() -> SetTitle(" Entry ");
0540     line_breakdown_gaus_ratio_hist -> GetXaxis() -> SetNdivisions(505);
0541
0542     line_breakdown_gaus_width_hist = new TH1F("","line_breakdown_gaus_width_hist",200,0,100);
0543     line_breakdown_gaus_width_hist -> GetXaxis() -> SetTitle("Width [mm]");
0544     line_breakdown_gaus_width_hist -> GetYaxis() -> SetTitle(" Entry ");
0545     line_breakdown_gaus_width_hist -> GetXaxis() -> SetNdivisions(505);
0546
0547     gaus_width_Nclu = new TH2F("","gaus_width_Nclu",200,0,10000,200,0,100);
0548     gaus_width_Nclu -> GetXaxis() -> SetTitle(" # of clusters ");
0549     gaus_width_Nclu -> GetYaxis() -> SetTitle("Gaus fit width [mm]");
0550     gaus_width_Nclu -> GetXaxis() -> SetNdivisions(505);
0551
0552     gaus_rchi2_Nclu = new TH2F("","gaus_rchi2_Nclu",200,0,10000,200,0,1);
0553     gaus_rchi2_Nclu -> GetXaxis() -> SetTitle(" # of clusters ");
0554     gaus_rchi2_Nclu -> GetYaxis() -> SetTitle("Gaus fit #chi2^{2}/NDF");
0555     gaus_rchi2_Nclu -> GetXaxis() -> SetNdivisions(505);
0556
0557     final_fit_width = new TH2F("","final_fit_width",200,0,10000,200,0,400);
0558     final_fit_width -> GetXaxis() -> SetTitle(" # of clusters ");
0559     final_fit_width -> GetYaxis() -> SetTitle("Fit width [mm]");
0560     final_fit_width -> GetXaxis() -> SetNdivisions(505);
0561
0562     N_track_candidate_Nclu = new TH2F("","N_track_candidate_Nclu",200,0,10000,200,0,13000);
0563     N_track_candidate_Nclu -> GetXaxis() -> SetTitle(" # of clusters ");
0564     N_track_candidate_Nclu -> GetYaxis() -> SetTitle("N tracklet candidates");
0565     N_track_candidate_Nclu -> GetXaxis() -> SetNdivisions(505);
0566
0567     temp_event_zvtx = new TH1F("","Z vertex dist",125,zvtx_hist_l,zvtx_hist_r);
0568     temp_event_zvtx -> GetXaxis() -> SetTitle("Z vertex position (mm)");
0569     temp_event_zvtx -> GetYaxis() -> SetTitle("Entry");
0570
0571     peak_group_width_hist = new TH1F("","peak_group_width_hist",100,0,300);
0572     peak_group_width_hist -> GetXaxis() -> SetTitle("Width [mm]");
0573     peak_group_width_hist -> GetYaxis() -> SetTitle("Entry");
0574     peak_group_width_hist -> GetXaxis() -> SetNdivisions(505);
0575
0576     peak_group_ratio_hist = new TH1F("","peak_group_ratio_hist",110,0,1.1);
0577     peak_group_ratio_hist -> GetXaxis() -> SetTitle("Peak group ratio");
0578     peak_group_ratio_hist -> GetYaxis() -> SetTitle("Entry");
0579     peak_group_ratio_hist -> GetXaxis() -> SetNdivisions(505);
0580
0581     N_group_hist  = new TH1F("","N_group_hist",30,0,30);
0582     N_group_hist -> GetXaxis() -> SetTitle("N group post cut");
0583     N_group_hist -> GetYaxis() -> SetTitle("Entry");
0584     N_group_hist -> GetXaxis() -> SetNdivisions(505);
0585
0586     peak_group_detail_width_hist = new TH1F("","peak_group_detail_width_hist",200,0,300);
0587     peak_group_detail_width_hist -> GetXaxis() -> SetTitle("Width [mm]");
0588     peak_group_detail_width_hist -> GetYaxis() -> SetTitle("Entry");
0589     peak_group_detail_width_hist -> GetXaxis() -> SetNdivisions(505);
0590
0591     peak_group_detail_ratio_hist  = new TH1F("","peak_group_detail_ratio_hist",110,0,1.1);
0592     peak_group_detail_ratio_hist -> GetXaxis() -> SetTitle("Peak group ratio");
0593     peak_group_detail_ratio_hist -> GetYaxis() -> SetTitle("Entry");
0594     peak_group_detail_ratio_hist -> GetXaxis() -> SetNdivisions(505);
0595
0596     N_group_detail_hist = new TH1F("","N_group_detail_hist",30,0,30);
0597     N_group_detail_hist -> GetXaxis() -> SetTitle("N group post cut");
0598     N_group_detail_hist -> GetYaxis() -> SetTitle("Entry");
0599     N_group_detail_hist -> GetXaxis() -> SetNdivisions(505);
0600
0601     evt_select_track_phi = new TH1F("","evt_select_track_phi",361,0,361);
0602     evt_select_track_phi -> GetXaxis() -> SetTitle("Track phi [degree]");
0603     evt_select_track_phi -> GetYaxis() -> SetTitle("Entry");
0604     evt_select_track_phi -> GetXaxis() -> SetNdivisions(505);
0605
0606     evt_phi_diff_inner_phi = new TH2F("","evt_phi_diff_inner_phi",361,0,361,100,-1.5,1.5);
0607     evt_phi_diff_inner_phi -> GetXaxis() -> SetTitle("Inner phi [degree]");
0608     evt_phi_diff_inner_phi -> GetYaxis() -> SetTitle("Inner - Outer [degree]");
0609     evt_phi_diff_inner_phi -> GetXaxis() -> SetNdivisions(505);
0610
0611     evt_inner_outer_phi = new TH2F("","evt_inner_outer_phi",120,0,360,120,0,360);
0612     evt_inner_outer_phi -> GetXaxis() -> SetTitle("Inner phi [degree]");
0613     evt_inner_outer_phi -> GetYaxis() -> SetTitle("Outer phi [degree]");
0614     evt_inner_outer_phi -> GetXaxis() -> SetNdivisions(505);
0615
0616     evt_phi_diff_1D = new TH1F("","evt_phi_diff_1D",100,-10,10);
0617     evt_phi_diff_1D -> GetXaxis() -> SetTitle("Inner - Outer [degree]");
0618     evt_phi_diff_1D -> GetYaxis() -> SetTitle("Entry");
0619     evt_phi_diff_1D -> GetXaxis() -> SetNdivisions(505);
0620 }
0621
0622 void INTTZvtx::InitCanvas()
0623 {
0624     c2 = new TCanvas("","",4000,1600);
0625     c2 -> cd();
0626     pad_xy = new TPad(Form("pad_xy"), "", 0.0, 0.5, 0.2, 1.0);
0627     Characterize_Pad(pad_xy, 0.15, 0.1, 0.1, 0.2, 0, 0);
0628     pad_xy -> Draw();
0629
0630     pad_rz = new TPad(Form("pad_rz"), "", 0.2, 0.5, 0.40, 1.0);
0631     Characterize_Pad(pad_rz, 0.15, 0.1, 0.1, 0.2, 0, 0);
0632     pad_rz -> Draw();
0633
0634     pad_z = new TPad(Form("pad_z"), "", 0.40, 0.5, 0.6, 1.0);
0635     Characterize_Pad(pad_z, 0.15, 0.1, 0.1, 0.2, 0, 0);
0636     pad_z -> Draw();
0637
0638     pad_z_hist = new TPad(Form("pad_z_hist"), "", 0.6, 0.5, 0.8, 1.0);
0639     Characterize_Pad(pad_z_hist, 0.15, 0.1, 0.1, 0.2, 0, 0);
0640     pad_z_hist -> Draw();
0641
0642     pad_z_line = new TPad(Form("pad_z_line"), "", 0.8, 0.5, 1, 1.0);
0643     Characterize_Pad(pad_z_line, 0.15, 0.1, 0.1, 0.2, 0, 0);
0644     pad_z_line -> Draw();
0645
0646     pad_phi_diff = new TPad(Form("pad_phi_diff"), "", 0.0, 0.0, 0.2, 0.5);
0647     Characterize_Pad(pad_phi_diff, 0.15, 0.1, 0.1, 0.2, 0, 0);
0648     pad_phi_diff -> Draw();
0649
0650     pad_track_phi = new TPad(Form("pad_track_phi"), "", 0.2, 0.0, 0.40, 0.5);
0651     Characterize_Pad(pad_track_phi, 0.15, 0.1, 0.1, 0.2, 0, 0);
0652     pad_track_phi -> Draw();
0653
0654     pad_inner_outer_phi = new TPad(Form("pad_inner_outer_phi"), "", 0.4, 0.0, 0.60, 0.5);
0655     Characterize_Pad(pad_inner_outer_phi, 0.15, 0.1, 0.1, 0.2, 0, 0);
0656     pad_inner_outer_phi -> SetLogz(1);
0657     pad_inner_outer_phi -> Draw();
0658
0659     pad_phi_diff_1D = new TPad(Form("pad_phi_diff_1D"), "", 0.6, 0.0, 0.80, 0.5);
0660     Characterize_Pad(pad_phi_diff_1D, 0.15, 0.1, 0.1, 0.2, 0, 0);
0661     // pad_phi_diff_1D -> SetLogz(1);
0662     pad_phi_diff_1D -> Draw();
0663
0664
0665
0666     c1 = new TCanvas("","",950,800);
0667     c1 -> cd();
0668
0669     pad_EvtZDist = new TPad("pad_EvtZDist", "pad_EvtZDist", 0.0, 0.0, 1.0, 1.0);
0670     // Characterize_Pad(pad_EvtZDist, 0.15, 0.1, 0.1, 0.2, 0, 0);
0671     pad_EvtZDist -> Draw();
0672
0673     pad_ZoomIn_EvtZDist = new TPad("pad_ZoomIn_EvtZDist", "pad_ZoomIn_EvtZDist", 0.52, 0.15+0.1, 0.82, 0.5+0.1);
0674     Characterize_Pad(pad_ZoomIn_EvtZDist, 0.15, 0.1, 0.1, 0.2, 0, 0);
0675     pad_ZoomIn_EvtZDist -> SetFillColor(0); // Set background color of the pad to white
0676     pad_ZoomIn_EvtZDist -> Draw();
0677 }
0678
0679 void INTTZvtx::InitTreeOut()
0680 {
0681     out_file = new TFile(Form("%s/INTT_zvtx.root",out_folder_directory.c_str()),"RECREATE");
0682     tree_out =  new TTree ("tree_Z", "INTT Z info.");
0683
0684     tree_out -> Branch("eID",&out_eID);
0685     tree_out -> Branch("bco_full",&bco_full_out);
0686     tree_out -> Branch("nclu_inner",&N_cluster_inner_out);
0687     tree_out -> Branch("nclu_outer",&N_cluster_outer_out);
0688     tree_out -> Branch("nclu_south",&out_N_cluster_south);
0689     tree_out -> Branch("nclu_north",&out_N_cluster_north);
0690     tree_out -> Branch("ES_zvtx",&out_ES_zvtx);                   // note : effective sigma, pol0 fit z-vertex
0691     tree_out -> Branch("ES_zvtxE",&out_ES_zvtxE);                 // note : effective sigma, pol0 fit z-vertex error
0692     tree_out -> Branch("ES_rangeL",&out_ES_rangeL);               // note : effective sigma, selected range left
0693     tree_out -> Branch("ES_rangeR",&out_ES_rangeR);               // note : effective sigma, selected range right
0694     tree_out -> Branch("ES_N_good",&out_ES_N_good);               // note : effective sigma, number of z-vertex candidates in the range
0695     tree_out -> Branch("ES_Width_density",&out_ES_width_density); // note : effective sigma, N good z-vertex candidates divided by width
0696     tree_out -> Branch("LB_Gaus_Mean_mean",&out_LB_Gaus_Mean_mean);              // note : Line break loose offset - gaus mean
0697     tree_out -> Branch("LB_Gaus_Mean_meanE",&out_LB_Gaus_Mean_meanE);            // note : Line break loose offset - gaus mean error
0698     tree_out -> Branch("LB_Gaus_Mean_chi2", &out_LB_Gaus_Mean_chi2);             // note : Line break loose offset - reduce chi2
0699     tree_out -> Branch("LB_Gaus_Mean_width",&out_LB_Gaus_Mean_width);            // note : Line break loose offset - width
0700     tree_out -> Branch("LB_Gaus_Mean_mean_vec", &fit_mean_mean_vec);                // note : Line break loose offset - the set of mean from the fitting with different fit ranges
0701     tree_out -> Branch("LB_Gaus_Mean_width_vec", &fit_mean_width_vec);              // note : Line break loose offset - the set of width from the fitting with different fit ranges
0702     tree_out -> Branch("LB_Gaus_Mean_chi2_vec", &fit_mean_reducedChi2_vec);         // note : Line break loose offset - the set of reduced chi2 from the fitting with different fit ranges
0703     tree_out -> Branch("LB_Gaus_Width_width",&out_LB_Gaus_Width_width);                 // note : Line break tight offset - gaus width
0704     tree_out -> Branch("LB_Gaus_Width_offset", &out_LB_Gaus_Width_offset);              // note : Line break tight offset - offset
0705     tree_out -> Branch("LB_Gaus_Width_size_width", &out_LB_Gaus_Width_size_width);      // note : Line break tight offset - norm. height / width
0706     tree_out -> Branch("LB_geo_mean", &out_LB_geo_mean);          // note : Line break peak position directly from the distribution (with the bin width 0.5 mm)
0707     tree_out -> Branch("good_zvtx_tag", &out_good_zvtx_tag);
0708     tree_out -> Branch("mid_cut_Ngroup",     &out_mid_cut_Ngroup);            // note : mid cut Ngroup
0709     tree_out -> Branch("mid_cut_peak_width", &out_mid_cut_peak_width);        // note : mid cut peak width
0710     tree_out -> Branch("mid_cut_peak_ratio", &out_mid_cut_peak_ratio);        // note : mid cut peak ratio
0711     tree_out -> Branch("LB_cut_Ngroup",     &out_LB_cut_Ngroup);         // note : LB cut Ngroup
0712     tree_out -> Branch("LB_cut_peak_width", &out_LB_cut_peak_width);     // note : LB cut peak width
0713     tree_out -> Branch("LB_cut_peak_ratio", &out_LB_cut_peak_ratio);     // note : LB cut peak ratio
0714     tree_out -> Branch("MC_true_zvtx",&MC_true_zvtx);
0715     tree_out -> Branch("Centrality_float",&out_centrality_float);
0716     tree_out -> Branch("MBD_reco_z", &out_MBD_reco_z);
0717     tree_out -> Branch("is_min_bias", &out_is_min_bias);
0718     tree_out -> Branch("is_min_bias_wozdc", &out_is_min_bias_wozdc);
0719     tree_out -> Branch("MBD_north_charge_sum", &out_MBD_north_charge_sum);
0720     tree_out -> Branch("MBD_south_charge_sum", &out_MBD_south_charge_sum);
0721     tree_out -> Branch("evtZ_hist_cut_content_vec", &out_evtZ_hist_cut_content_vec); // note : this is for the training test of the z-vertex finder
0722 }
0723
0724 void INTTZvtx::InitRest()
0725 {
0726     gaus_fit_vec.clear();
0727     gaus_fit_cm_vec.clear();
0728     for (int i = 0; i < 7; i++)
0729     {
0730         gaus_fit_vec.push_back(new TF1("gaus_fit_vec",gaus_func,evt_possible_z_range.first,evt_possible_z_range.second,4));
0731         gaus_fit_vec[i] -> SetLineColor(TColor::GetColor(color_code[i].c_str()));
0732         gaus_fit_vec[i] -> SetLineWidth(1);
0733         gaus_fit_vec[i] -> SetNpx(1000);
0734
0735         gaus_fit_cm_vec.push_back(new TF1("gaus_fit_cm_vec",gaus_func,evt_possible_z_range.first,evt_possible_z_range.second,4));
0736         gaus_fit_cm_vec[i] -> SetLineColor(TColor::GetColor(color_code[i].c_str()));
0737         gaus_fit_cm_vec[i] -> SetLineWidth(1);
0738         gaus_fit_cm_vec[i] -> SetNpx(1000);
0739
0740     }
0741
0742     gaus_fit = new TF1("gaus_fit",gaus_func,evt_possible_z_range.first,evt_possible_z_range.second,4);
0743     gaus_fit -> SetLineColor(2);
0744     gaus_fit -> SetLineWidth(1);
0745     gaus_fit -> SetNpx(1000);
0746
0747     gaus_fit_cm = new TF1("gaus_fit_cm",gaus_func,evt_possible_z_range.first,evt_possible_z_range.second,4);
0748     gaus_fit_cm -> SetLineColor(2);
0749     gaus_fit_cm -> SetLineWidth(1);
0750     gaus_fit_cm -> SetNpx(1000);
0751
0752     gaus_fit_width_cm = new TF1("gaus_fit_width_cm",gaus_func,evt_possible_z_range.first,evt_possible_z_range.second,4);
0753     gaus_fit_width_cm -> SetLineColor(3);
0754     gaus_fit_width_cm -> SetLineStyle(9);
0755     gaus_fit_width_cm -> SetLineWidth(1);
0756     gaus_fit_width_cm -> SetNpx(1000);
0757
0758     gaus_fit_2 = new TF1("gaus_fit_2",gaus_func,evt_possible_z_range.first,evt_possible_z_range.second,4);
0759     gaus_fit_2 -> SetLineColor(2);
0760     gaus_fit_2 -> SetLineWidth(2);
0761     gaus_fit_2 -> SetNpx(1000);
0762
0763     double_gaus_fit = new TF1("double_gaus_fit",double_gaus_func,-600,600,7);
0764     double_gaus_fit -> SetLineColor(3);
0765     double_gaus_fit -> SetLineWidth(1);
0766     double_gaus_fit -> SetNpx(1000);
0767
0768     eff_sig_range_line = new TLine();
0769     eff_sig_range_line -> SetLineWidth(2);
0770     eff_sig_range_line -> SetLineColor(TColor::GetColor("#A08144"));
0771     eff_sig_range_line -> SetLineStyle(2);
0772
0773     final_fit_range_line = new TLine();
0774     final_fit_range_line -> SetLineWidth(1);
0775     final_fit_range_line -> SetLineColor(TColor::GetColor("#44a04d"));
0776     // final_fit_range_line -> SetLineStyle(2);
0777
0778     zvtx_finder = new TF1("zvtx_finder","pol0",-1,20000);
0779     zvtx_finder -> SetLineColor(2);
0780     zvtx_finder -> SetLineWidth(1);
0781
0782     ltx = new TLatex();
0783     ltx->SetNDC();
0784     ltx->SetTextSize(0.045);
0785     ltx->SetTextAlign(31);
0786
0787     draw_text = new TLatex();
0788     draw_text -> SetNDC();
0789     draw_text -> SetTextSize(0.03);
0790
0791     ladder_line = new TLine();
0792     ladder_line -> SetLineWidth(1);
0793
0794     track_line = new TLine();
0795     track_line -> SetLineWidth(1);
0796     track_line -> SetLineColorAlpha(38,0.5);
0797
0798     coord_line = new TLine();
0799     coord_line -> SetLineWidth(1);
0800     coord_line -> SetLineColor(16);
0801     coord_line -> SetLineStyle(2);
0802
0803     rand_gen = new TRandom3(0);
0804
0805     ch_pos_DB = new InttConversion(conversion_mode_BD[geo_mode_id], peek);
0806 }
0807
0808 void INTTZvtx::ProcessEvt(
0809     int event_i,
0810     vector<clu_info> temp_sPH_inner_nocolumn_vec,
0811     vector<clu_info> temp_sPH_outer_nocolumn_vec,
0812     vector<vector<double>> temp_sPH_nocolumn_vec,
0813     vector<vector<double>> temp_sPH_nocolumn_rz_vec,
0814     int NvtxMC,
0815     double TrigZvtxMC,
0816     bool PhiCheckTag,
0817     Long64_t bco_full,
0818     float centrality_float,
0819     double MBD_reco_z,
0820     int is_min_bias,
0821     int is_min_bias_wozdc,
0822     double MBD_north_charge_sum,
0823     double MBD_south_charge_sum
0824 )
0825 {
0826     out_eID = event_i;
0827     N_cluster_inner_out = -1;
0828     N_cluster_outer_out = -1;
0829     bco_full_out = bco_full;
0830
0831     out_ES_zvtx = -1;
0832     out_ES_zvtxE = -1;
0833     out_ES_rangeL = -1;
0834     out_ES_rangeR = -1;
0835     out_ES_N_good = -1;
0836     out_ES_width_density = -1;
0837
0838     out_LB_Gaus_Mean_mean = -1;
0839     out_LB_Gaus_Mean_meanE = -1;
0840     out_LB_Gaus_Mean_chi2 = -1;
0841     out_LB_Gaus_Mean_width = -1;
0842
0843     out_LB_Gaus_Width_width = -1;
0844     out_LB_Gaus_Width_size_width = -1;
0845     out_LB_Gaus_Width_offset = -1;
0846
0847     out_mid_cut_Ngroup = -1;
0848     out_mid_cut_peak_width = -1;
0849     out_mid_cut_peak_ratio = -1;
0850
0851     out_LB_cut_Ngroup = -1;
0852     out_LB_cut_peak_width = -1;
0853     out_LB_cut_peak_ratio = -1;
0854
0855     out_LB_geo_mean = -1;
0856     out_good_zvtx_tag = 0;
0857     MC_true_zvtx = TrigZvtxMC * 10.;
0858
0859     out_LB_refit = -1;
0860     out_LB_refitE = -1;
0861
0862     good_zvtx_tag = false;
0863     good_zvtx_tag_int = 0;
0864     loose_offset_peak = -999; // note : unit [mm]
0865     loose_offset_peakE = -999; // note : unit [mm]
0866
0867     out_centrality_float = centrality_float;
0868     out_MBD_reco_z = MBD_reco_z;
0869     out_is_min_bias = is_min_bias;
0870     out_is_min_bias_wozdc = is_min_bias_wozdc;
0871     out_MBD_north_charge_sum = MBD_north_charge_sum;
0872     out_MBD_south_charge_sum = MBD_south_charge_sum;
0873
0874     out_N_cluster_north = 0;
0875     out_N_cluster_south = 0;
0876
0877     out_evtZ_hist_cut_content_vec.clear();
0878
0879
0880     if (event_i%10 == 0 && print_message_opt == true) {cout<<"In INTTZvtx class, running event : "<<event_i<<endl;}
0881
0882     total_NClus = temp_sPH_inner_nocolumn_vec.size() + temp_sPH_outer_nocolumn_vec.size();
0883
0884     if (total_NClus < zvtx_cal_require)  {N_cluster_inner_out = -1; tree_out -> Fill(); return; cout<<"return confirmation"<<endl;}
0885
0886     if (run_type == "MC" && NvtxMC != 1) {N_cluster_inner_out = -2; tree_out -> Fill(); return; cout<<"In INTTZvtx class, event : "<<event_i<<" Nvtx : "<<NvtxMC<<" Nvtx more than one "<<endl; }
0887     if (PhiCheckTag == false)            {N_cluster_inner_out = -3; tree_out -> Fill(); return; cout<<"In INTTZvtx class, event : "<<event_i<<" Nvtx : "<<NvtxMC<<" Not full phi has hits "<<endl;}
0888
0889     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)
0890     {
0891         N_cluster_inner_out = -4;
0892         tree_out -> Fill();
0893         return;
0894         printf("In INTTZvtx class, event : %i, low clu continue, NClus : %lu \n", event_i, total_NClus);
0895     }
0896
0897     // note : put the cluster into the phi map, the first bool is for the cluster usage.
0898     // note : false means the cluster is not used
0899     for (int inner_i = 0; inner_i < temp_sPH_inner_nocolumn_vec.size(); inner_i++) {
0900         Clus_InnerPhi_Offset = (temp_sPH_inner_nocolumn_vec[inner_i].y - beam_origin.second < 0) ? atan2(temp_sPH_inner_nocolumn_vec[inner_i].y - beam_origin.second, temp_sPH_inner_nocolumn_vec[inner_i].x - beam_origin.first) * (180./TMath::Pi()) + 360 : atan2(temp_sPH_inner_nocolumn_vec[inner_i].y - beam_origin.second, temp_sPH_inner_nocolumn_vec[inner_i].x - beam_origin.first) * (180./TMath::Pi());
0901         // cout<<"inner clu phi : "<<Clus_InnerPhi_Offset<<" origin: "<< temp_sPH_inner_nocolumn_vec[inner_i].phi <<endl;
0902         // cout<<" ("<<Clus_InnerPhi_Offset<<", "<< temp_sPH_inner_nocolumn_vec[inner_i].phi<<")" <<endl;
0903         inner_clu_phi_map[ int(Clus_InnerPhi_Offset) ].push_back({false,temp_sPH_inner_nocolumn_vec[inner_i]});
0904
0905         if (temp_sPH_inner_nocolumn_vec[inner_i].z > 0) {out_N_cluster_north += 1;}
0906         else {out_N_cluster_south += 1;}
0907     }
0908     for (int outer_i = 0; outer_i < temp_sPH_outer_nocolumn_vec.size(); outer_i++) {
0909         Clus_OuterPhi_Offset = (temp_sPH_outer_nocolumn_vec[outer_i].y - beam_origin.second < 0) ? atan2(temp_sPH_outer_nocolumn_vec[outer_i].y - beam_origin.second, temp_sPH_outer_nocolumn_vec[outer_i].x - beam_origin.first) * (180./TMath::Pi()) + 360 : atan2(temp_sPH_outer_nocolumn_vec[outer_i].y - beam_origin.second, temp_sPH_outer_nocolumn_vec[outer_i].x - beam_origin.first) * (180./TMath::Pi());
0910         outer_clu_phi_map[ int(Clus_OuterPhi_Offset) ].push_back({false,temp_sPH_outer_nocolumn_vec[outer_i]});
0911
0912         if (temp_sPH_outer_nocolumn_vec[outer_i].z > 0) {out_N_cluster_north += 1;}
0913         else {out_N_cluster_south += 1;}
0914     }
0915
0916     // note : for the mega cluster finder, the 3/4-cluster tracklet that happens at the overlap region
0917     // mega_track_finder -> FindMegaTracks(inner_clu_phi_map, outer_clu_phi_map, {MC_true_zvtx,0.5}, centrality_map[centrality_float]);
0918
0919     double good_pair_count = 0;
0920
0921     for (int inner_phi_i = 0; inner_phi_i < 360; inner_phi_i++) // note : each phi cell (1 degree)
0922     {
0923         // note : N cluster in this phi cell
0924         for (int inner_phi_clu_i = 0; inner_phi_clu_i < inner_clu_phi_map[inner_phi_i].size(); inner_phi_clu_i++)
0925         {
0926             if (inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].first == true) {continue;}
0927
0928             Clus_InnerPhi_Offset = (inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.y - beam_origin.second < 0) ? atan2(inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.y - beam_origin.second, inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.x - beam_origin.first) * (180./TMath::Pi()) + 360 : atan2(inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.y - beam_origin.second, inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.x - beam_origin.first) * (180./TMath::Pi());
0929
0930             // todo: change the outer phi scan range
0931             // note : the outer phi index, -1, 0, 1
0932             // note : the outer phi index, -2, -1, 0, 1, 2, for the scan test
0933             for (int scan_i = -2; scan_i < 3; scan_i++)
0934             {
0935                 int true_scan_i = ((inner_phi_i + scan_i) < 0) ? 360 + (inner_phi_i + scan_i) : ((inner_phi_i + scan_i) > 359) ? (inner_phi_i + scan_i)-360 : inner_phi_i + scan_i;
0936
0937                 // note : N clusters in that outer phi cell
0938                 for (int outer_phi_clu_i = 0; outer_phi_clu_i < outer_clu_phi_map[true_scan_i].size(); outer_phi_clu_i++)
0939                 {
0940                     if (outer_clu_phi_map[true_scan_i][outer_phi_clu_i].first == true) {continue;}
0941
0942                     Clus_OuterPhi_Offset = (outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.y - beam_origin.second < 0) ? atan2(outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.y - beam_origin.second, outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.x - beam_origin.first) * (180./TMath::Pi()) + 360 : atan2(outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.y - beam_origin.second, outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.x - beam_origin.first) * (180./TMath::Pi());
0943                     double delta_phi = get_delta_phi(Clus_InnerPhi_Offset, Clus_OuterPhi_Offset);
0944
0945                     double delta_phi_correct;
0946                     if (delta_phi_inner_correction_hist != nullptr) {
0947                         delta_phi_correct = delta_phi - delta_phi_inner_correction_hist->GetBinContent(delta_phi_inner_correction_hist->GetXaxis()->FindBin(Clus_InnerPhi_Offset)); //note: special_tag
0948                     }
0949                     else {
0950                         delta_phi_correct = delta_phi;
0951                     }
0952
0953                     evt_phi_diff_inner_phi -> Fill(Clus_InnerPhi_Offset, delta_phi_correct);
0954                     evt_inner_outer_phi -> Fill(Clus_InnerPhi_Offset, Clus_OuterPhi_Offset);
0955                     evt_phi_diff_1D -> Fill(delta_phi_correct);
0956
0957                     if (fabs(delta_phi_correct) < phi_diff_cut)
0958                     {
0959                         double DCA_sign = calculateAngleBetweenVectors(
0960                             outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.x, outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.y,
0961                             inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.x, inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.y,
0962                             beam_origin.first, beam_origin.second
0963                         );
0964
0965                         if (DCA_cut.first < DCA_sign && DCA_sign < DCA_cut.second){
0966                             good_pair_count += 1;
0967                             // used_outer_check[outer_i] = 1; //note : this outer cluster was already used!
0968
0969                             // note : we basically transform the coordinate from cartesian to cylinder
0970                             // note : we should set the offset first, otherwise it provides the bias
0971                             // todo : which point should be used, DCA point or vertex xy ? Has to be studied
0972                             pair<double,double> z_range_info = Get_possible_zvtx(
0973                                 0., // get_radius(beam_origin.first,beam_origin.second),
0974                                 {get_radius(inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.x - beam_origin.first, inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.y - beam_origin.second), inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.z}, // note : unsign radius
0975                                 {get_radius(outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.x - beam_origin.first, outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.y - beam_origin.second), outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.z}  // note : unsign radius
0976                             );
0977
0978                             // note : try to remove some crazy background candidates. Can be a todo
0979                             if (evt_possible_z_range.first < z_range_info.first && z_range_info.first < evt_possible_z_range.second)
0980                             {
0981                                 N_comb.push_back(good_comb_id);
0982                                 N_comb_e.push_back(0);
0983                                 // N_comb_phi.push_back( (inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.phi + outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.phi)/2. );
0984                                 N_comb_phi.push_back( get_track_phi(Clus_InnerPhi_Offset, delta_phi) );
0985                                 z_mid.push_back(z_range_info.first);
0986                                 z_range.push_back(z_range_info.second);
0987
0988                                 evt_possible_z -> Fill(z_range_info.first);
0989
0990                                 // note : fill the line_breakdown histogram as well as a vector for the width determination
0991                                 // line_breakdown(line_breakdown_hist,{z_range_info.first - z_range_info.second, z_range_info.first + z_range_info.second}); // note: special_tag
0992                                 trapezoidal_line_breakdown( // note: special_tag
0993                                     line_breakdown_hist, // note: special_tag
0994                                     // note : inner_r and inner_z
0995                                     get_radius(inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.x - beam_origin.first, inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.y - beam_origin.second), // note: special_tag
0996                                     inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.z,// note: special_tag
0997                                     // note : outer_r and outer_z
0998                                     get_radius(outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.x - beam_origin.first, outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.y - beam_origin.second), // note: special_tag
0999                                     outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.z// note: special_tag
1000                                 );
1001
1002                                 good_comb_id += 1;
1003                             }
1004
1005                         }
1006                     }
1007
1008                 }
1009             } // note : end of outer clu loop
1010         }
1011
1012     } // note : end of inner clu loop
1013
1014     // cout<<"good pair count : "<<good_pair_count<<endl;
1015
1016     // if (event_i == 906) {
1017     //     for (int hist_i = 0; hist_i < line_breakdown_hist->GetNbinsX(); hist_i++){
1018     //         cout<<line_breakdown_hist->GetBinContent(hist_i+1)<<","<<endl;
1019     //     }
1020     // }
1021     N_track_candidate_Nclu -> Fill(total_NClus, N_comb.size());
1022
1023     if (N_comb.size() > zvtx_cal_require)
1024     {
1025         // note : this is for data, basically there is some background candidates that are not from the z-vertex
1026         // todo : it's  forcaibly written here. Will need to be revisited
1027         for (int bin_i = 0; bin_i < line_breakdown_hist->GetNbinsX(); bin_i++){ // note : special_tag
1028             if (line_breakdown_hist -> GetBinCenter(bin_i+1) < -500 || line_breakdown_hist -> GetBinCenter(bin_i+1) > 500){ // note : special_tag
1029                 line_breakdown_hist -> SetBinContent(bin_i+1,0); // note : special_tag
1030             } // note : special_tag
1031         } // note : special_tag
1032
1033         N_group_info = find_Ngroup(evt_possible_z);
1034         N_group_info_detail = find_Ngroup(line_breakdown_hist);
1035         out_evtZ_hist_cut_content_vec = get_half_hist_vec(line_breakdown_hist);
1036
1037         double Half_N_group_half_width = fabs(N_group_info_detail[3] - N_group_info_detail[2]) / 2.;
1038
1039         FakeClone1D(line_breakdown_hist, line_breakdown_hist_cm);
1040         FakeClone1D(line_breakdown_hist, line_breakdown_hist_zoomin);
1041         FakeClone1D(line_breakdown_hist, line_breakdown_hist_cm_zoomin);
1042
1043         double line_breakdown_hist_max_content = line_breakdown_hist -> GetBinContent( line_breakdown_hist -> GetMaximumBin() );
1044         double line_breakdown_hist_max_center = line_breakdown_hist -> GetBinCenter( line_breakdown_hist -> GetMaximumBin() );
1045
1046         // note : first fit is for the width, so apply the constraints on the Gaussian offset
1047         gaus_fit -> SetParameters(line_breakdown_hist_max_content, line_breakdown_hist_max_center, Half_N_group_half_width, 0);
1048         gaus_fit -> SetParLimits(0,0,100000);  // note : size
1049         gaus_fit -> SetParLimits(2,5,10000);   // note : Width
1050         gaus_fit -> SetParLimits(3,0,10000);   // note : offset
1051         double width_fit_range_l = line_breakdown_hist_max_center - Half_N_group_half_width * 0.6;
1052         double width_fit_range_r = line_breakdown_hist_max_center + Half_N_group_half_width * 0.6;
1053         // todo : try to use single gaus to fit the distribution, and try to only fit the peak region (peak - 100 mm + peak + 100 mm)
1054         line_breakdown_hist -> Fit(gaus_fit, "NQ", "", width_fit_range_l, width_fit_range_r);
1055         gaus_fit_width_cm -> SetParameters(gaus_fit -> GetParameter(0), gaus_fit -> GetParameter(1) * 0.1, gaus_fit -> GetParameter(2) * 0.1, gaus_fit -> GetParameter(3));
1056         gaus_fit_width_cm -> SetRange( width_fit_range_l * 0.1, width_fit_range_r * 0.1);
1057         tight_offset_peak = gaus_fit -> GetParameter(1);
1058         tight_offset_width = fabs(gaus_fit -> GetParameter(2));
1059         final_selection_widthU = (tight_offset_peak + tight_offset_width);
1060         final_selection_widthD = (tight_offset_peak - tight_offset_width);
1061         gaus_fit_offset = gaus_fit -> GetParameter(3); gaus_fit -> SetParameter(3,0); // note : in order to calculate the integration
1062         gaus_ratio = (fabs(gaus_fit -> GetParameter(0))/gaus_fit->Integral(-600,600)) / fabs(gaus_fit -> GetParameter(2));
1063         gaus_fit -> SetParameter(3, gaus_fit_offset); // note : put the offset back to the function
1064
1065         for (int fit_i = 0; fit_i < gaus_fit_vec.size(); fit_i++)
1066         {
1067             gaus_fit_vec[fit_i] -> SetParameters(line_breakdown_hist_max_content, line_breakdown_hist_max_center, Half_N_group_half_width, 0);
1068             gaus_fit_vec[fit_i] -> SetParLimits(0,0,100000);  // note : size
1069             gaus_fit_vec[fit_i] -> SetParLimits(2,5,10000);   // note : Width
1070             gaus_fit_vec[fit_i] -> SetParLimits(3,-200,10000);   // note : offset
1071             double N_width_ratio = 0.2 + 0.15 * fit_i;
1072             double mean_fit_range_l = line_breakdown_hist_max_center - Half_N_group_half_width * N_width_ratio;
1073             double mean_fit_range_r = line_breakdown_hist_max_center + Half_N_group_half_width * N_width_ratio;
1074             // todo : try to use single gaus to fit the distribution, and try to only fit the peak region (peak - 100 mm + peak + 100 mm)
1075             line_breakdown_hist -> Fit(gaus_fit_vec[fit_i], "NQ", "", mean_fit_range_l, mean_fit_range_r);
1076             gaus_fit_vec[fit_i] -> SetRange(mean_fit_range_l, mean_fit_range_r);
1077
1078             fit_mean_mean_vec.push_back(gaus_fit_vec[fit_i] -> GetParameter(1));
1079             fit_mean_reducedChi2_vec.push_back(gaus_fit_vec[fit_i] -> GetChisquare() / double(gaus_fit_vec[fit_i] -> GetNDF()));
1080             fit_mean_width_vec.push_back(fabs(gaus_fit_vec[fit_i] -> GetParameter(2)));
1081
1082             gaus_fit_cm_vec[fit_i] -> SetParameters(
1083                 gaus_fit_vec[fit_i] -> GetParameter(0),
1084                 gaus_fit_vec[fit_i] -> GetParameter(1) * 0.1,
1085                 gaus_fit_vec[fit_i] -> GetParameter(2) * 0.1,
1086                 gaus_fit_vec[fit_i] -> GetParameter(3)
1087             );
1088             gaus_fit_cm_vec[fit_i] -> SetRange(mean_fit_range_l * 0.1, mean_fit_range_r * 0.1);
1089         }
1090
1091
1092         // note : second fit is for the peak position, therefore, loose the constraints on the Gaussian offset
1093         // gaus_fit -> SetParameters(line_breakdown_hist_max_content, line_breakdown_hist_max_center, Half_N_group_half_width, 0);
1094         // gaus_fit -> SetParLimits(0,0,100000);  // note : size
1095         // gaus_fit -> SetParLimits(2,5,10000);   // note : Width
1096         // gaus_fit -> SetParLimits(3,-200,10000);   // note : offset
1097         // double mean_fit_range_l = line_breakdown_hist_max_center - best_width; // note : special_tag // todo
1098         // double mean_fit_range_r = line_breakdown_hist_max_center + best_width; // note : special_tag // todo change back
1099         // // todo : try to use single gaus to fit the distribution, and try to only fit the peak region (peak - 100 mm + peak + 100 mm)
1100         // line_breakdown_hist -> Fit(gaus_fit, "NQ", "", mean_fit_range_l, mean_fit_range_r);
1101         // gaus_fit -> SetRange(gaus_fit -> GetParameter(1) - gaus_fit -> GetParameter(2), gaus_fit -> GetParameter(1) + gaus_fit -> GetParameter(2));
1102         // // line_breakdown_hist -> Fit(gaus_fit, "NQ", "", N_group_info_detail[2]-10, N_group_info_detail[3]+10);
1103         // loose_offset_peak = gaus_fit -> GetParameter(1);
1104         // loose_offset_peakE = gaus_fit -> GetParError(1);
1105
1106         loose_offset_peak = vector_average(fit_mean_mean_vec);
1107         loose_offset_peakE = vector_stddev(fit_mean_mean_vec);
1108
1109         line_breakdown_hist_zoomin -> GetXaxis() -> SetRangeUser(
1110             line_breakdown_hist_zoomin -> GetBinCenter(line_breakdown_hist_zoomin -> GetMaximumBin()) - Half_N_group_half_width,
1111             line_breakdown_hist_zoomin -> GetBinCenter(line_breakdown_hist_zoomin -> GetMaximumBin()) + Half_N_group_half_width
1112         );
1113         line_breakdown_hist_zoomin -> SetMinimum(line_breakdown_hist_zoomin -> GetBinContent(line_breakdown_hist_zoomin -> GetMaximumBin()) * 0.5);
1114         line_breakdown_hist_zoomin -> SetMaximum(line_breakdown_hist_zoomin -> GetBinContent(line_breakdown_hist_zoomin -> GetMaximumBin()) * 1.5);
1115
1116         // note : ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
1117         // note : first fit is for the width, so apply the constraints on the Gaussian offset
1118         // gaus_fit_cm -> SetParameters(gaus_fit -> GetParameter(0), gaus_fit -> GetParameter(1) * 0.1, gaus_fit -> GetParameter(2) * 0.1, gaus_fit -> GetParameter(3));
1119         // gaus_fit_cm -> SetRange( mean_fit_range_l * 0.1, mean_fit_range_r * 0.1);
1120         // note : ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
1121
1122         line_breakdown_hist_cm_zoomin -> GetXaxis() -> SetRangeUser(
1123             gaus_fit_cm_vec.front() -> GetParameter(1) - Half_N_group_half_width * 0.1,
1124             gaus_fit_cm_vec.front() -> GetParameter(1) + Half_N_group_half_width * 0.1
1125         );
1126         line_breakdown_hist_cm_zoomin -> SetMinimum( gaus_fit_cm_vec.front() -> GetParameter(0) * 0.5);
1127         line_breakdown_hist_cm_zoomin -> SetMaximum( gaus_fit_cm_vec.front() -> GetParameter(0) * 1.5);
1128
1129
1130         // note : eff sigma method, relatively sensitive to the background
1131         // note : use z-mid to do the effi_sig, because that line_breakdown takes too long time
1132         temp_event_zvtx_info = sigmaEff_avg(z_mid,Integrate_portion);
1133         for (int track_i = 0; track_i < N_comb.size(); track_i++) {
1134
1135             if ( N_group_info[2] <= z_mid[track_i] && z_mid[track_i] <= N_group_info[3] ){
1136                 eff_N_comb.push_back(N_comb[track_i]);
1137                 eff_z_mid.push_back(z_mid[track_i]);
1138                 eff_N_comb_e.push_back(N_comb_e[track_i]);
1139                 eff_z_range.push_back(z_range[track_i]);
1140             }
1141
1142             if (final_selection_widthD <= z_mid[track_i] && z_mid[track_i] <= final_selection_widthU){
1143                 // note : for monitoring the the phi distribution that is used for the z vertex determination.
1144                 // note : in principle, I expect it should be something uniform.
1145                 evt_select_track_phi -> Fill(N_comb_phi[track_i]);
1146             }
1147
1148             // if ((z_mid[track_i] - z_range[track_i]) > final_selection_widthU || (z_mid[track_i] + z_range[track_i]) < final_selection_widthD) {continue;}
1149             // else {
1150             //     // eff_N_comb.push_back(N_comb[track_i]);
1151             //     // eff_z_mid.push_back(z_mid[track_i]);
1152             //     // eff_N_comb_e.push_back(N_comb_e[track_i]);
1153             //     // eff_z_range.push_back(z_range[track_i]);
1154             // }
1155         }
1156
1157         z_range_gr = new TGraphErrors(eff_N_comb.size(),&eff_N_comb[0],&eff_z_mid[0],&eff_N_comb_e[0],&eff_z_range[0]);
1158         z_range_gr -> Fit(zvtx_finder,"NQ","",0,N_comb[N_comb.size() - 1]);
1159         width_density_par = (double(eff_N_comb.size()) / fabs(temp_event_zvtx_info[2] - temp_event_zvtx_info[1]));
1160
1161
1162
1163         // if ( zvtx_QA_width.first < tight_offset_width && tight_offset_width < zvtx_QA_width.second && 100 < fabs(N_group_info_detail[3] - N_group_info_detail[2]) && fabs(N_group_info_detail[3] - N_group_info_detail[2]) < 190){
1164         //     if (N_group_info[0] < 4 && N_group_info[1] >= 0.6 && N_group_info_detail[0] < 7 && N_group_info_detail[1] > 0.9) {good_zvtx_tag = true;}
1165         //     else {good_zvtx_tag = false;}
1166         // }
1167         // else {good_zvtx_tag = false;}
1168
1169         if ( zvtx_QA_width.first < tight_offset_width && tight_offset_width < zvtx_QA_width.second && 100 < fabs(N_group_info_detail[3] - N_group_info_detail[2]) && fabs(N_group_info_detail[3] - N_group_info_detail[2]) < 190){
1170             good_zvtx_tag = true;
1171         }
1172         else {good_zvtx_tag = false;}
1173
1174         good_zvtx_tag_int = (good_zvtx_tag == true) ? 1 : 0; // note : so stupid way to convert bool to int...
1175
1176         // if (good_zvtx_tag == false) {
1177         //     good_zvtx_tag = (zvtx_QA_width.first < fabs(gaus_fit -> GetParameter(2)) && fabs(gaus_fit -> GetParameter(2)) < zvtx_QA_width.second && gaus_ratio > zvtx_QA_ratio && width_density_par > width_density_cut) ? true : false;
1178         // } // note : kinda double check
1179
1180         // note : for the group information post the background cut
1181         peak_group_width_hist -> Fill(fabs(N_group_info[3] - N_group_info[2]));
1182         peak_group_ratio_hist -> Fill(N_group_info[1]);
1183         N_group_hist -> Fill(N_group_info[0]);
1184         peak_group_detail_width_hist -> Fill(fabs(N_group_info_detail[3] - N_group_info_detail[2]));
1185         peak_group_detail_ratio_hist -> Fill(N_group_info_detail[1]);
1186         N_group_detail_hist -> Fill(N_group_info_detail[0]);
1187
1188         // note : final
1189         final_fit_width -> Fill(total_NClus, fabs(final_selection_widthU - final_selection_widthD)); // note : from LB gaus for the moment
1190         final_zvtx =  loose_offset_peak; // zvtx_finder -> GetParameter(0);
1191
1192         // note : gaus fit on the linebreak
1193         gaus_width_Nclu -> Fill(total_NClus, tight_offset_width);
1194         gaus_rchi2_Nclu -> Fill(total_NClus, vector_average(fit_mean_reducedChi2_vec));
1195         line_breakdown_gaus_ratio_hist -> Fill(gaus_ratio);
1196         line_breakdown_gaus_width_hist -> Fill(tight_offset_width);
1197
1198         // note : the effi sig on the z_mid vector
1199         zvtx_evt_width_corre -> Fill(total_NClus, fabs(temp_event_zvtx_info[2] - temp_event_zvtx_info[1]));
1200         width_density -> Fill( width_density_par );
1201         ES_width -> Fill(fabs(temp_event_zvtx_info[2] - temp_event_zvtx_info[1]));
1202         ES_width_ratio -> Fill(total_NClus / fabs(temp_event_zvtx_info[2] - temp_event_zvtx_info[1]));
1203
1204         // note : regarding to the zvtx determination, pol0 fit with error bar considered, fit range given by eff sigma method
1205         zvtx_evt_fitError -> Fill(fabs(zvtx_finder -> GetParError(0)));
1206         zvtx_evt_fitError_corre -> Fill(total_NClus, fabs(zvtx_finder -> GetParError(0)));
1207
1208         if ( good_zvtx_tag ){
1209
1210             zvtx_evt_nclu_corre -> Fill(total_NClus, final_zvtx);
1211             avg_event_zvtx -> Fill(final_zvtx);
1212             avg_event_zvtx_vec.push_back(final_zvtx);
1213             if (run_type == "MC") {
1214                 Z_resolution -> Fill( final_zvtx - (TrigZvtxMC * 10.) );
1215                 Z_resolution_vec.push_back( final_zvtx - (TrigZvtxMC * 10.) );
1216                 Z_resolution_Nclu -> Fill( total_NClus , final_zvtx - (TrigZvtxMC * 10.) );
1217                 Z_resolution_pos -> Fill(TrigZvtxMC * 10., final_zvtx - (TrigZvtxMC * 10.));
1218                 if (total_NClus > high_multi_line) {Z_resolution_pos_cut -> Fill(TrigZvtxMC * 10., final_zvtx - (TrigZvtxMC * 10.));}
1219             }
1220
1221
1222         }
1223
1224         // note : output the root tree
1225         out_eID = event_i;
1226         N_cluster_inner_out = temp_sPH_inner_nocolumn_vec.size();
1227         N_cluster_outer_out = temp_sPH_outer_nocolumn_vec.size();
1228
1229         out_ES_zvtx = zvtx_finder -> GetParameter(0);
1230         out_ES_zvtxE = zvtx_finder -> GetParError(0);
1231         out_ES_rangeL = temp_event_zvtx_info[1];
1232         out_ES_rangeR = temp_event_zvtx_info[2];
1233         out_ES_N_good = eff_N_comb.size();
1234         out_ES_width_density = width_density_par;
1235
1236         out_LB_Gaus_Mean_mean = loose_offset_peak;
1237         out_LB_Gaus_Mean_meanE = loose_offset_peakE; // note : -> loose one
1238         out_LB_Gaus_Mean_chi2 = vector_average(fit_mean_reducedChi2_vec); // note : -> loose one
1239         out_LB_Gaus_Mean_width = vector_average(fit_mean_width_vec); // note : -> loose one
1240
1241         out_LB_Gaus_Width_width = tight_offset_width;
1242         out_LB_Gaus_Width_offset = gaus_fit_offset;
1243         out_LB_Gaus_Width_size_width = gaus_ratio;
1244
1245         out_mid_cut_Ngroup = N_group_info[0];
1246         out_mid_cut_peak_ratio = N_group_info[1];
1247         out_mid_cut_peak_width = fabs(N_group_info[3] - N_group_info[2]) / 2.;
1248
1249         out_LB_cut_Ngroup = N_group_info_detail[0];
1250         out_LB_cut_peak_ratio = N_group_info_detail[1];
1251         out_LB_cut_peak_width = fabs(N_group_info_detail[3] - N_group_info_detail[2]) / 2.;
1252
1253         out_centrality_float = centrality_float;
1254         out_MBD_reco_z = MBD_reco_z;
1255         out_is_min_bias = is_min_bias;
1256         out_is_min_bias_wozdc = is_min_bias_wozdc;
1257         out_MBD_north_charge_sum = MBD_north_charge_sum;
1258         out_MBD_south_charge_sum = MBD_south_charge_sum;
1259
1260         out_LB_geo_mean = LB_geo_mean(line_breakdown_hist, { (tight_offset_peak - tight_offset_width), (tight_offset_peak + tight_offset_width) }, event_i);
1261         out_good_zvtx_tag = good_zvtx_tag;
1262         bco_full_out = bco_full;
1263         MC_true_zvtx = TrigZvtxMC * 10.;
1264         tree_out -> Fill();
1265
1266         out_LB_refit = -1;  // note : not in the out tree
1267         out_LB_refitE = -1; // note : not in the out tree
1268
1269     } // note : if N good tracks in xy found > certain value
1270     else {N_cluster_inner_out = -5; tree_out -> Fill();}
1271
1272     if (N_comb.size() > zvtx_cal_require && draw_event_display == true)
1273     {
1274         c2 -> cd();
1275         temp_event_xy = new TGraph(temp_sPH_nocolumn_vec[0].size(),&temp_sPH_nocolumn_vec[0][0],&temp_sPH_nocolumn_vec[1][0]);
1276         temp_event_xy -> SetTitle("INTT event display X-Y plane");
1277         temp_event_xy -> GetXaxis() -> SetLimits(-150,150);
1278         temp_event_xy -> GetYaxis() -> SetRangeUser(-150,150);
1279         temp_event_xy -> GetXaxis() -> SetTitle("X [mm]");
1280         temp_event_xy -> GetYaxis() -> SetTitle("Y [mm]");
1281         temp_event_xy -> SetMarkerStyle(20);
1282         temp_event_xy -> SetMarkerColor(2);
1283         temp_event_xy -> SetMarkerSize(1);
1284         temp_event_rz = new TGraph(temp_sPH_nocolumn_rz_vec[0].size(),&temp_sPH_nocolumn_rz_vec[0][0],&temp_sPH_nocolumn_rz_vec[1][0]);
1285         temp_event_rz -> SetTitle("INTT event display r-Z plane");
1286         temp_event_rz -> GetXaxis() -> SetLimits(-500,500);
1287         temp_event_rz -> GetYaxis() -> SetRangeUser(-150,150);
1288         temp_event_rz -> GetXaxis() -> SetTitle("Z [mm]");
1289         temp_event_rz -> GetYaxis() -> SetTitle("Radius [mm]");
1290         temp_event_rz -> SetMarkerStyle(20);
1291         temp_event_rz -> SetMarkerColor(2);
1292         temp_event_rz -> SetMarkerSize(1);
1293
1294         // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1295
1296         pad_xy -> cd();
1297         temp_bkg(pad_xy, peek, beam_origin, ch_pos_DB);
1298         temp_event_xy -> Draw("p same");
1299         for (int track_i = 0; track_i < good_track_xy_vec.size(); track_i++){
1300             track_line -> DrawLine(good_track_xy_vec[track_i][0],good_track_xy_vec[track_i][1],good_track_xy_vec[track_i][2],good_track_xy_vec[track_i][3]);
1301         }
1302         track_line -> Draw("l same");
1303         draw_text -> DrawLatex(0.2, 0.85, Form("eID : %i, inner Ncluster : %zu, outer Ncluster : %zu",event_i,temp_sPH_inner_nocolumn_vec.size(),temp_sPH_outer_nocolumn_vec.size()));
1304
1305         // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1306
1307         pad_rz -> cd();
1308         temp_event_rz -> Draw("ap");
1309         // eff_sig_range_line -> DrawLine(temp_event_zvtx_info[0],-150,temp_event_zvtx_info[0],150);
1310         coord_line -> DrawLine(0,-150,0,150);
1311         coord_line -> DrawLine(-500,0,500,0);
1312         for (int track_i = 0; track_i < good_track_rz_vec.size(); track_i++){
1313             track_line -> DrawLine(good_track_rz_vec[track_i][0],good_track_rz_vec[track_i][1],good_track_rz_vec[track_i][2],good_track_rz_vec[track_i][3]);
1314         }
1315         draw_text -> DrawLatex(0.2, 0.85, Form("Negative radius : Clu_{outer} > 180^{0}"));
1316         // draw_text -> DrawLatex(0.2, 0.81, Form("EffSig avg : %.2f mm",temp_event_zvtx_info[0]));
1317
1318         // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1319         // cout<<"test tag 2-5"<<endl;
1320         pad_z -> cd();
1321         z_range_gr_draw = new TGraphErrors(N_comb.size(),&N_comb[0],&z_mid[0],&N_comb_e[0],&z_range[0]);
1322         z_range_gr_draw -> GetYaxis() -> SetRangeUser(-650,650);
1323         z_range_gr_draw -> GetXaxis() -> SetTitle("Index");
1324         z_range_gr_draw -> GetYaxis() -> SetTitle("Z [mm]");
1325         z_range_gr_draw -> SetMarkerStyle(20);
1326         z_range_gr_draw -> Draw("ap");
1327
1328         // eff_sig_range_line -> DrawLine(z_range_gr_draw->GetXaxis()->GetXmin(),temp_event_zvtx_info[1],z_range_gr_draw->GetXaxis()->GetXmax(),temp_event_zvtx_info[1]);
1329         // eff_sig_range_line -> DrawLine(z_range_gr_draw->GetXaxis()->GetXmin(),temp_event_zvtx_info[2],z_range_gr_draw->GetXaxis()->GetXmax(),temp_event_zvtx_info[2]);
1330         eff_sig_range_line -> DrawLine(z_range_gr_draw->GetXaxis()->GetXmin(),N_group_info[2],z_range_gr_draw->GetXaxis()->GetXmax(),N_group_info[2]);
1331         eff_sig_range_line -> DrawLine(z_range_gr_draw->GetXaxis()->GetXmin(),N_group_info[3],z_range_gr_draw->GetXaxis()->GetXmax(),N_group_info[3]);
1332         draw_text -> DrawLatex(0.2, 0.82, Form("#color[2]{Event Zvtx %.2f mm, error : #pm%.2f}", zvtx_finder -> GetParameter(0), zvtx_finder -> GetParError(0)));
1333         draw_text -> DrawLatex(0.2, 0.78, Form("#color[2]{Width density : %.6f}", ( width_density_par )));
1334         draw_text -> DrawLatex(0.2, 0.74, Form("#color[2]{Width (%.0f%%) : %.2f to %.2f mm #rightarrow %.2f mm}", Integrate_portion*100., temp_event_zvtx_info[2] , temp_event_zvtx_info[1], fabs(temp_event_zvtx_info[2] - temp_event_zvtx_info[1])/2.));
1335
1336         // todo : the final range finder is here
1337         // final_fit_range_line -> DrawLine(z_range_gr_draw->GetXaxis()->GetXmin(),temp_event_zvtx_info[1],z_range_gr_draw->GetXaxis()->GetXmax(),temp_event_zvtx_info[1]);
1338         // final_fit_range_line -> DrawLine(z_range_gr_draw->GetXaxis()->GetXmin(),temp_event_zvtx_info[2],z_range_gr_draw->GetXaxis()->GetXmax(),temp_event_zvtx_info[2]);
1339
1340         // z_range_gr_draw -> Draw("p same");
1341         zvtx_finder -> Draw("lsame");
1342
1343         // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1344         pad_z_hist -> cd();
1345         evt_possible_z -> SetMaximum(evt_possible_z -> GetBinContent( evt_possible_z->GetMaximumBin() ) * 1.4);
1346         evt_possible_z -> Draw("hist");
1347         eff_sig_range_line -> DrawLine(evt_possible_z->GetXaxis()->GetXmin(), evt_possible_z -> GetBinContent(evt_possible_z->GetMaximumBin())/2., evt_possible_z->GetXaxis()->GetXmax(), evt_possible_z -> GetBinContent(evt_possible_z->GetMaximumBin())/2.);
1348         draw_text -> DrawLatex(0.2, 0.82, Form("N group : %.0f", N_group_info[0]));
1349         draw_text -> DrawLatex(0.2, 0.78, Form("Main peak ratio : %.2f", N_group_info[1]));
1350         draw_text -> DrawLatex(0.2, 0.74, Form("good z tag : %i", good_zvtx_tag));
1351
1352         // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1353         pad_z_line -> cd();
1354         line_breakdown_hist -> SetMinimum(0);
1355         line_breakdown_hist -> SetMaximum(line_breakdown_hist -> GetBinContent( line_breakdown_hist->GetMaximumBin() ) * 2);
1356         line_breakdown_hist -> Draw("hist");
1357         gaus_fit_vec.back() -> Draw("lsame");
1358         if (!good_zvtx_tag) {final_fit_range_line -> DrawLine(line_breakdown_hist->GetXaxis()->GetXmin(), line_breakdown_hist -> GetMaximum(), line_breakdown_hist->GetXaxis()->GetXmax(), line_breakdown_hist -> GetMinimum());}
1359         final_fit_range_line -> DrawLine(final_selection_widthD, line_breakdown_hist -> GetMinimum(),final_selection_widthD, line_breakdown_hist -> GetMaximum());
1360         final_fit_range_line -> DrawLine(final_selection_widthU, line_breakdown_hist -> GetMinimum(),final_selection_widthU, line_breakdown_hist -> GetMaximum());
1361         draw_text -> DrawLatex(0.2, 0.82, Form("Gaus mean %.2f mm", loose_offset_peak));
1362         draw_text -> DrawLatex(0.2, 0.78, Form("Width : %.2f mm", tight_offset_width));
1363         draw_text -> DrawLatex(0.2, 0.74, Form("Reduced #chi2 : %.3f", gaus_fit_vec.back() -> GetChisquare() / double(gaus_fit_vec.back() -> GetNDF())));
1364         draw_text -> DrawLatex(0.2, 0.70, Form("Norm. entry / Width : %.6f mm",  gaus_ratio  ));
1365         draw_text -> DrawLatex(0.2, 0.66, Form("LB Geo mean : %.3f mm",  out_LB_geo_mean  ));
1366
1367         if (N_group_info_detail[0] != -1){
1368             eff_sig_range_line -> DrawLine(line_breakdown_hist->GetXaxis()->GetXmin(), line_breakdown_hist -> GetBinContent(line_breakdown_hist->GetMaximumBin())/2., line_breakdown_hist->GetXaxis()->GetXmax(), line_breakdown_hist -> GetBinContent(line_breakdown_hist->GetMaximumBin())/2.);
1369             draw_text -> DrawLatex(0.2, 0.62, Form("N group : %.0f", N_group_info_detail[0]));
1370             draw_text -> DrawLatex(0.2, 0.58, Form("Main peak ratio : %.2f", N_group_info_detail[1]));
1371         }
1372
1373         if (run_type == "MC") {draw_text -> DrawLatex(0.2, 0.54, Form("True MCz : %.3f mm",  TrigZvtxMC * 10.));}
1374         draw_text -> DrawLatex(0.2, 0.5, Form("INTTz: %.3f, MBDz: %.3f, diff: %.3f",final_zvtx, MBD_reco_z, final_zvtx - MBD_reco_z));
1375
1376         // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1377         pad_phi_diff -> cd();
1378         evt_phi_diff_inner_phi -> Draw("colz0");
1379
1380         // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1381         pad_track_phi -> cd();
1382         evt_select_track_phi -> Draw("hist");
1383
1384         // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1385         pad_inner_outer_phi -> cd();
1386         evt_inner_outer_phi -> Draw("colz0");
1387
1388         // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1389         pad_phi_diff_1D -> cd();
1390         evt_phi_diff_1D -> Draw("hist");
1391
1392
1393         // temp_event_zvtx -> Draw("hist");
1394         // // zvtx_fitting -> Draw("lsame");
1395         // eff_sig_range_line -> DrawLine(temp_event_zvtx_info[1],0,temp_event_zvtx_info[1],temp_event_zvtx -> GetMaximum()*1.05);
1396         // eff_sig_range_line -> DrawLine(temp_event_zvtx_info[2],0,temp_event_zvtx_info[2],temp_event_zvtx -> GetMaximum()*1.05);
1397         // draw_text -> DrawLatex(0.2, 0.85, Form("eID : %i, Total event hit : %i, innter Ncluster : %i, outer Ncluster : %i",event_i,N_hits,temp_sPH_inner_nocolumn_vec.size(),temp_sPH_outer_nocolumn_vec.size()));
1398         // // draw_text -> DrawLatex(0.2, 0.84, Form("Gaus fit mean : %.3f mm",zvtx_fitting -> GetParameter(1)));
1399         // draw_text -> DrawLatex(0.2, 0.82, Form("EffSig min : %.2f mm, max : %.2f mm",temp_event_zvtx_info[1],temp_event_zvtx_info[2]));
1400         // draw_text -> DrawLatex(0.2, 0.79, Form("EffSig avg : %.2f mm",temp_event_zvtx_info[0]));
1401
1402         if(draw_event_display && (event_i % print_rate) == 0 /*&& good_zvtx_tag == true*/){c2 -> Print(Form("%s/temp_event_display.pdf",out_folder_directory.c_str()));}
1403         // else if(good_zvtx_tag == false){c2 -> Print(Form("%s/temp_event_display.pdf",out_folder_directory.c_str()));}
1404         else if(draw_event_display && (final_zvtx > 0 || final_zvtx < -450)) {cout<<"In INTTZvtx class, event :"<<event_i<<" weird zvtx "<<endl; c2 -> Print(Form("%s/temp_event_display.pdf",out_folder_directory.c_str()));}
1405         // else if(total_NClus > 5000) {cout<<"In INTTZvtx class, event :"<<event_i<<" high Nclus "<<endl; c2 -> Print(Form("%s/temp_event_display.pdf",out_folder_directory.c_str()));}
1406         else if ( draw_event_display && ( (final_zvtx - MBD_reco_z) < -45 || (final_zvtx - MBD_reco_z) > 30 ) && total_NClus > 1500 ) { cout<<"In INTTZvtx class, event :"<<event_i<<" High NClus, poor INTTrecoZ with MBDrecoZ, diff : "<< final_zvtx - MBD_reco_z <<endl; c2 -> Print(Form("%s/temp_event_display.pdf",out_folder_directory.c_str())); }
1407         else if ( draw_event_display && run_type == "MC" && fabs(final_zvtx - (TrigZvtxMC * 10.)) > 2. && total_NClus > 3000) { cout<<"In INTTZvtx class, event :"<<event_i<<" High NClus, poor Z : "<< fabs(final_zvtx - (TrigZvtxMC * 10.)) <<endl; c2 -> Print(Form("%s/temp_event_display.pdf",out_folder_directory.c_str())); }
1408         else if ( draw_event_display && run_type == "MC" && fabs(final_zvtx - (TrigZvtxMC * 10.)) > 2.) { cout<<"In INTTZvtx class, event :"<<event_i<<" low NClus, poor Z : "<< fabs(final_zvtx - (TrigZvtxMC * 10.)) <<endl; c2 -> Print(Form("%s/temp_event_display.pdf",out_folder_directory.c_str())); }
1409
1410         pad_xy -> Clear();
1411         pad_rz -> Clear();
1412         pad_z  -> Clear();
1413         pad_z_hist -> Clear();
1414         pad_z_line -> Clear();
1415         pad_phi_diff -> Clear();
1416         pad_track_phi -> Clear();
1417         pad_inner_outer_phi -> Clear();
1418         pad_phi_diff_1D -> Clear();
1419
1420         c1 -> cd();
1421         pad_EvtZDist -> cd();
1422
1423         line_breakdown_hist_cm -> SetMinimum(0);
1424         line_breakdown_hist_cm -> SetMaximum(line_breakdown_hist_cm -> GetBinContent( line_breakdown_hist_cm->GetMaximumBin() ) * 1.65);
1425         line_breakdown_hist_cm -> Draw("hist");
1426         // gaus_fit_width_cm -> Draw("l same");
1427         for (int fit_i = 0; fit_i < gaus_fit_cm_vec.size(); fit_i++){
1428             gaus_fit_cm_vec[gaus_fit_cm_vec.size() - fit_i - 1] -> Draw("l same");
1429         }
1430         ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1431
1432         draw_text -> DrawLatex(0.2, 0.9, Form("Event ID: %i", event_i));
1433         draw_text -> DrawLatex(0.2, 0.86, Form("NClus: %li", total_NClus));
1434         draw_text -> DrawLatex(0.2, 0.82, Form("Reco. vertex Z: %.3f cm, StdDev: %.3f cm", loose_offset_peak * 0.1, loose_offset_peakE * 0.1));
1435         if (run_type == "MC") {
1436             draw_text -> DrawLatex(0.2, 0.78, Form("True vertex Z: %.3f cm",TrigZvtxMC));
1437             draw_text -> DrawLatex(0.2, 0.74, Form("Difference: %.3f cm", loose_offset_peak * 0.1 - TrigZvtxMC));
1438         }
1439
1440         pad_ZoomIn_EvtZDist -> cd();
1441         line_breakdown_hist_cm_zoomin -> Draw("hist");
1442         for (int fit_i = 0; fit_i < gaus_fit_cm_vec.size(); fit_i++){
1443             gaus_fit_cm_vec[gaus_fit_cm_vec.size() - fit_i - 1] -> Draw("l same");
1444         }
1445
1446         // todo : the print rate is here
1447         if(draw_event_display && (event_i % print_rate) == 0 /*&& good_zvtx_tag == true*/){c1 -> Print(Form("%s/evt_linebreak_display.pdf",out_folder_directory.c_str()));}
1448         else if(draw_event_display && (final_zvtx > 0 || final_zvtx < -450)) {cout<<"In INTTZvtx class, event :"<<event_i<<" weird zvtx "<<endl; c1 -> Print(Form("%s/evt_linebreak_display.pdf",out_folder_directory.c_str()));}
1449         else if ( draw_event_display && ( (final_zvtx - MBD_reco_z) < -45 || (final_zvtx - MBD_reco_z) > 30 ) && total_NClus > 1500 ) { cout<<"In INTTZvtx class, event :"<<event_i<<" High NClus, poor INTTrecoZ with MBDrecoZ, diff : "<< final_zvtx - MBD_reco_z <<endl; c1 -> Print(Form("%s/evt_linebreak_display.pdf",out_folder_directory.c_str())); }
1450         else if ( draw_event_display && run_type == "MC" && fabs(final_zvtx - (TrigZvtxMC * 10.)) > 2. && total_NClus > 3000) { cout<<"In INTTZvtx class, event :"<<event_i<<" High NClus, poor Z : "<< fabs(final_zvtx - (TrigZvtxMC * 10.)) <<endl; c1 -> Print(Form("%s/evt_linebreak_display.pdf",out_folder_directory.c_str())); }
1451         else if ( draw_event_display && run_type == "MC" && fabs(final_zvtx - (TrigZvtxMC * 10.)) > 2.) { cout<<"In INTTZvtx class, event :"<<event_i<<" low NClus, poor Z : "<< fabs(final_zvtx - (TrigZvtxMC * 10.)) <<endl; c1 -> Print(Form("%s/evt_linebreak_display.pdf",out_folder_directory.c_str())); }
1452
1453         pad_EvtZDist -> Clear();
1454         pad_ZoomIn_EvtZDist -> Clear();
1455
1456         // if (pad_ZoomIn_EvtZDist) {
1457         //     pad_ZoomIn_EvtZDist -> Delete();
1458         // }
1459
1460         // temp_event_xy -> Delete();
1461         // temp_event_rz -> Delete();
1462         // z_range_gr_draw -> Delete();
1463
1464     }
1465
1466 }
1467
1468 void INTTZvtx::ClearEvt()
1469 {
1470     good_comb_id = 0;
1471     out_N_cluster_north = 0;
1472     out_N_cluster_south = 0;
1473
1474     // note : ultra-stupid way to avoid the errors
1475     z_range_gr_draw = new TGraphErrors(); z_range_gr_draw -> Delete();
1476     z_range_gr = new TGraphErrors(); z_range_gr -> Delete();
1477     temp_event_xy = new TGraph(); temp_event_xy -> Delete();
1478     temp_event_rz = new TGraph(); temp_event_rz -> Delete();
1479     bkg = new TGraph(); bkg -> Delete();
1480
1481     N_comb_phi.clear();
1482     N_comb.clear();
1483     N_comb_e.clear();
1484     z_mid.clear();
1485     z_range.clear();
1486
1487     eff_N_comb.clear();
1488     eff_z_mid.clear();
1489     eff_N_comb_e.clear();
1490     eff_z_range.clear();
1491     line_breakdown_vec.clear();
1492     temp_event_zvtx_info = {0,-1000,-999.99};
1493
1494     LB_N_comb.clear();
1495     LB_z_mid.clear();
1496     LB_N_comb_e.clear();
1497     LB_z_range.clear();
1498     N_group_info.clear();
1499     N_group_info_detail = {-1., -1., -1., -1.};
1500
1501     fit_mean_mean_vec.clear();
1502     fit_mean_reducedChi2_vec.clear();
1503     fit_mean_width_vec.clear();
1504
1505     evt_possible_z -> Reset("ICESM");
1506     line_breakdown_hist -> Reset("ICESM");
1507     line_breakdown_hist_zoomin -> Reset("ICESM");
1508     line_breakdown_hist_cm -> Reset("ICESM");
1509     line_breakdown_hist_cm_zoomin -> Reset("ICESM");
1510     evt_phi_diff_inner_phi -> Reset("ICESM");
1511     evt_select_track_phi -> Reset("ICESM");
1512     evt_inner_outer_phi -> Reset("ICESM");
1513     evt_phi_diff_1D -> Reset("ICESM");
1514
1515     inner_clu_phi_map.clear();
1516     outer_clu_phi_map.clear();
1517     inner_clu_phi_map = vector<vector<pair<bool,clu_info>>>(360);
1518     outer_clu_phi_map = vector<vector<pair<bool,clu_info>>>(360);
1519
1520     out_evtZ_hist_cut_content_vec.clear();
1521
1522     // mega_track_finder -> ClearEvt();
1523
1524     // note : this is the distribution for full run
1525     // line_breakdown_gaus_ratio_hist -> Reset("ICESM");
1526 }
1527
1528 void INTTZvtx::PrintPlots()
1529 {
1530     if (draw_event_display) {c2 -> Print(Form("%s/temp_event_display.pdf)",out_folder_directory.c_str()));}
1531     c2 -> Clear();
1532
1533     if (draw_event_display) {c1 -> Print(Form("%s/evt_linebreak_display.pdf)",out_folder_directory.c_str()));}
1534     c1 -> Clear();
1535
1536     cout<<"avg_event_zvtx_vec size : "<<avg_event_zvtx_vec.size()<<endl;
1537
1538     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1539
1540     c1 -> cd();
1541     vector<float> avg_event_zvtx_info = {0,0,0};
1542     if (avg_event_zvtx_vec.size() > 10) {avg_event_zvtx_info = sigmaEff_avg(avg_event_zvtx_vec,Integrate_portion_final);}
1543
1544     avg_event_zvtx -> SetMinimum( 0 );  avg_event_zvtx -> SetMaximum( avg_event_zvtx->GetBinContent(avg_event_zvtx->GetMaximumBin()) * 1.5 );
1545     avg_event_zvtx -> Draw("hist");
1546
1547     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1548     // ltx->DrawLatex(0.54, 0.86, Form("Run %s",run_ID.c_str()));
1549     // ltx->DrawLatex(0.54, 0.86, "Au+Au #sqrt{s_{NN}} = 200 GeV");
1550
1551     eff_sig_range_line -> DrawLine(avg_event_zvtx_info[1],0,avg_event_zvtx_info[1],avg_event_zvtx -> GetMaximum());
1552     eff_sig_range_line -> DrawLine(avg_event_zvtx_info[2],0,avg_event_zvtx_info[2],avg_event_zvtx -> GetMaximum());
1553     draw_text -> DrawLatex(0.21, 0.87, Form("EffSig min : %.2f mm",avg_event_zvtx_info[1]));
1554     draw_text -> DrawLatex(0.21, 0.83, Form("EffSig max : %.2f mm",avg_event_zvtx_info[2]));
1555     draw_text -> DrawLatex(0.21, 0.79, Form("EffSig avg : %.2f mm",avg_event_zvtx_info[0]));
1556     c1 -> Print(Form("%s/avg_event_zvtx.pdf",out_folder_directory.c_str()));
1557     c1 -> Clear();
1558
1559     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1560
1561     width_density -> Draw("hist");
1562     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1563     c1 -> Print(Form("%s/width_density.pdf",out_folder_directory.c_str()));
1564     c1 -> Clear();
1565
1566     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1567
1568     ES_width -> Draw("hist");
1569     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1570     c1 -> Print(Form("%s/ES_width.pdf",out_folder_directory.c_str()));
1571     c1 -> Clear();
1572
1573     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1574
1575     ES_width_ratio -> Draw("hist");
1576     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1577     c1 -> Print(Form("%s/ES_width_ratio.pdf",out_folder_directory.c_str()));
1578     c1 -> Clear();
1579
1580     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1581
1582     zvtx_evt_fitError -> Draw("hist");
1583     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1584     c1 -> Print(Form("%s/zvtx_evt_fitError.pdf",out_folder_directory.c_str()));
1585     c1 -> Clear();
1586
1587     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1588     if (run_type == "MC")
1589     {
1590         vector<float> Z_resolution_vec_info = sigmaEff_avg(Z_resolution_vec,Integrate_portion_final);
1591
1592         gaus_fit_2 -> SetParameters(Z_resolution -> GetBinContent( Z_resolution -> GetMaximumBin() ), Z_resolution -> GetBinCenter( Z_resolution -> GetMaximumBin() ), 3, 0);
1593         gaus_fit_2 -> SetParLimits(0,0,100000);  // note : size
1594         gaus_fit_2 -> SetParLimits(2,0,10000);   // note : Width
1595         gaus_fit_2 -> SetParLimits(3,0,10000);   // note : offset
1596         Z_resolution -> Fit(gaus_fit_2, "NQ", "", Z_resolution -> GetBinCenter( Z_resolution -> GetMaximumBin() ) - (2 * Z_resolution -> GetStdDev() ), Z_resolution -> GetBinCenter( Z_resolution -> GetMaximumBin() ) + (2 * Z_resolution -> GetStdDev() ) );
1597
1598         Z_resolution -> Draw("hist");
1599         gaus_fit_2 -> SetRange( gaus_fit_2->GetParameter(1) - gaus_fit_2->GetParameter(2) * 2.5, gaus_fit_2->GetParameter(1) + gaus_fit_2->GetParameter(2) * 2.5 );
1600         gaus_fit_2 -> Draw("lsame");
1601         eff_sig_range_line -> DrawLine(Z_resolution_vec_info[1],0,Z_resolution_vec_info[1],Z_resolution -> GetMaximum());
1602         eff_sig_range_line -> DrawLine(Z_resolution_vec_info[2],0,Z_resolution_vec_info[2],Z_resolution -> GetMaximum());
1603         draw_text -> DrawLatex(0.21, 0.87, Form("EffSig min : %.2f mm",Z_resolution_vec_info[1]));
1604         draw_text -> DrawLatex(0.21, 0.83, Form("EffSig max : %.2f mm",Z_resolution_vec_info[2]));
1605         draw_text -> DrawLatex(0.21, 0.79, Form("EffSig avg : %.2f mm",Z_resolution_vec_info[0]));
1606         draw_text -> DrawLatex(0.21, 0.71, Form("Gaus mean  : %.2f mm",gaus_fit_2 -> GetParameter(1)));
1607         draw_text -> DrawLatex(0.21, 0.67, Form("Gaus width : %.2f mm",fabs(gaus_fit_2 -> GetParameter(2))));
1608         ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1609         c1 -> Print(Form("%s/Z_resolution.pdf",out_folder_directory.c_str()));
1610         c1 -> Clear();
1611
1612         MC_z_diff_peak = gaus_fit_2 -> GetParameter(1);
1613         MC_z_diff_width = fabs(gaus_fit_2 -> GetParameter(2));
1614     }
1615
1616
1617     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1618
1619     if (run_type == "MC")
1620     {
1621         Z_resolution_Nclu -> Draw("colz0");
1622         ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1623         c1 -> Print(Form("%s/Z_resolution_Nclu.pdf",out_folder_directory.c_str()));
1624         c1 -> Clear();
1625     }
1626     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1627
1628     if (run_type == "MC")
1629     {
1630         Z_resolution_pos -> Draw("colz0");
1631         ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1632         c1 -> Print(Form("%s/Z_resolution_pos.pdf",out_folder_directory.c_str()));
1633         c1 -> Clear();
1634     }
1635     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1636
1637     if (run_type == "MC")
1638     {
1639         Z_resolution_pos_cut -> Draw("colz0");
1640         ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1641         c1 -> Print(Form("%s/Z_resolution_pos_cut.pdf",out_folder_directory.c_str()));
1642         c1 -> Clear();
1643     }
1644     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1645
1646     zvtx_evt_fitError_corre -> Draw("colz0");
1647     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1648     c1 -> Print(Form("%s/zvtx_evt_fitError_corre.pdf",out_folder_directory.c_str()));
1649     c1 -> Clear();
1650
1651     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1652
1653     zvtx_evt_nclu_corre -> Draw("colz0");
1654     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1655     c1 -> Print(Form("%s/zvtx_evt_nclu_corre.pdf",out_folder_directory.c_str()));
1656     c1 -> Clear();
1657
1658     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1659
1660     zvtx_evt_width_corre -> Draw("colz0");
1661     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1662     c1 -> Print(Form("%s/zvtx_evt_width_corre.pdf",out_folder_directory.c_str()));
1663     c1 -> Clear();
1664
1665     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1666
1667     gaus_width_Nclu -> Draw("colz0");
1668     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1669     c1 -> Print(Form("%s/gaus_width_Nclu.pdf",out_folder_directory.c_str()));
1670     c1 -> Clear();
1671
1672     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1673
1674     gaus_rchi2_Nclu -> Draw("colz0");
1675     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1676     c1 -> Print(Form("%s/gaus_rchi2_Nclu.pdf",out_folder_directory.c_str()));
1677     c1 -> Clear();
1678
1679     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1680
1681     final_fit_width -> Draw("colz0");
1682     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1683     c1 -> Print(Form("%s/final_fit_width.pdf",out_folder_directory.c_str()));
1684     c1 -> Clear();
1685
1686     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1687
1688     N_track_candidate_Nclu -> Draw("colz0");
1689     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1690     c1 -> Print(Form("%s/N_track_candidate_Nclu.pdf",out_folder_directory.c_str()));
1691     c1 -> Clear();
1692
1693     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1694     peak_group_width_hist -> Draw("hist");
1695     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1696     c1 -> Print(Form("%s/peak_group_width_hist.pdf",out_folder_directory.c_str()));
1697     c1 -> Clear();
1698
1699     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1700     peak_group_ratio_hist -> Draw("hist");
1701     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1702     c1 -> Print(Form("%s/peak_group_ratio_hist.pdf",out_folder_directory.c_str()));
1703     c1 -> Clear();
1704
1705     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1706     N_group_hist -> Draw("hist");
1707     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1708     c1 -> Print(Form("%s/N_group_hist.pdf",out_folder_directory.c_str()));
1709     c1 -> Clear();
1710
1711     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1712     peak_group_detail_width_hist -> Draw("hist");
1713     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1714     c1 -> Print(Form("%s/peak_group_detail_width_hist.pdf",out_folder_directory.c_str()));
1715     c1 -> Clear();
1716
1717     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1718     peak_group_detail_ratio_hist -> Draw("hist");
1719     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1720     c1 -> Print(Form("%s/peak_group_detail_ratio_hist.pdf",out_folder_directory.c_str()));
1721     c1 -> Clear();
1722
1723     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1724     N_group_detail_hist -> Draw("hist");
1725     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1726     c1 -> Print(Form("%s/N_group_detail_hist.pdf",out_folder_directory.c_str()));
1727     c1 -> Clear();
1728
1729     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1730     // gaus_fit -> SetParameters(line_breakdown_gaus_ratio_hist -> GetBinContent( line_breakdown_gaus_ratio_hist -> GetMaximumBin() ), line_breakdown_gaus_ratio_hist -> GetBinCenter( line_breakdown_gaus_ratio_hist -> GetMaximumBin() ), 0.00005, 0);
1731     // gaus_fit -> SetParLimits(0,0,100000);  // note : size
1732     // gaus_fit -> SetParLimits(2,0,10000);   // note : Width
1733     // gaus_fit -> SetParLimits(3,0,10000);   // note : offset
1734
1735     line_breakdown_gaus_ratio_hist -> Draw("hist");
1736     // line_breakdown_gaus_ratio_hist -> Fit(gaus_fit, "NQ", "", line_breakdown_gaus_ratio_hist -> GetBinCenter( line_breakdown_gaus_ratio_hist -> GetMaximumBin() ) - 1, line_breakdown_gaus_ratio_hist -> GetBinCenter( line_breakdown_gaus_ratio_hist -> GetMaximumBin() ) + 1);
1737     // gaus_fit -> Draw("lsame");
1738
1739     // draw_text -> DrawLatex(0.45, 0.82, Form("Gaus mean %.2f mm", gaus_fit -> GetParameter(1)));
1740     // draw_text -> DrawLatex(0.45, 0.78, Form("Width : %.2f mm", fabs(gaus_fit -> GetParameter(2))));
1741     // draw_text -> DrawLatex(0.45, 0.74, Form("Reduced #chi2 : %.3f", gaus_fit -> GetChisquare() / double(gaus_fit -> GetNDF())));
1742     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1743
1744     c1 -> Print(Form("%s/line_breakdown_gaus_ratio_hist.pdf",out_folder_directory.c_str()));
1745     c1 -> Clear();
1746
1747     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1748     gaus_fit -> SetParameters(line_breakdown_gaus_width_hist -> GetBinContent( line_breakdown_gaus_width_hist -> GetMaximumBin() ), line_breakdown_gaus_width_hist -> GetBinCenter( line_breakdown_gaus_width_hist -> GetMaximumBin() ), 8, 0);
1749     gaus_fit -> SetParLimits(0,0,100000);  // note : size
1750     gaus_fit -> SetParLimits(2,1,10000);   // note : Width
1751     gaus_fit -> SetParLimits(3,0,10000);   // note : offset
1752
1753     // todo : try to use single gaus to fit the distribution, and try to only fit the peak region (peak - 100 mm + peak + 100 mm)
1754     line_breakdown_gaus_width_hist -> Fit(gaus_fit, "NQ", "", line_breakdown_gaus_width_hist -> GetBinCenter( line_breakdown_gaus_width_hist -> GetMaximumBin() ) - 100, line_breakdown_gaus_width_hist -> GetBinCenter( line_breakdown_gaus_width_hist -> GetMaximumBin() ) + 100);
1755     line_breakdown_gaus_width_hist -> Draw("hist");
1756     gaus_fit -> SetRange(gaus_fit -> GetParameter(1) - gaus_fit -> GetParameter(2) * 1.5, gaus_fit -> GetParameter(1) + gaus_fit -> GetParameter(2) * 1.5);
1757     gaus_fit -> Draw("lsame");
1758
1759     draw_text -> DrawLatex(0.2, 0.82, Form("Gaus mean %.2f mm", gaus_fit -> GetParameter(1)));
1760     draw_text -> DrawLatex(0.2, 0.78, Form("Width : %.2f mm", fabs(gaus_fit -> GetParameter(2))));
1761     draw_text -> DrawLatex(0.2, 0.74, Form("Reduced #chi2 : %.3f", gaus_fit -> GetChisquare() / double(gaus_fit -> GetNDF())));
1762     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1763
1764     c1 -> Print(Form("%s/line_breakdown_gaus_width_hist.pdf",out_folder_directory.c_str()));
1765     c1 -> Clear();
1766 }
1767
1768 void INTTZvtx::EndRun()
1769 {
1770
1771     if (run_type == "MC")
1772     {
1773         gaus_fit_2 -> SetParameters(Z_resolution -> GetBinContent( Z_resolution -> GetMaximumBin() ), Z_resolution -> GetBinCenter( Z_resolution -> GetMaximumBin() ), 3, 0);
1774         gaus_fit_2 -> SetParLimits(0,0,100000);  // note : size
1775         gaus_fit_2 -> SetParLimits(2,0,10000);   // note : Width
1776         gaus_fit_2 -> SetParLimits(3,0,10000);   // note : offset
1777         Z_resolution -> Fit(gaus_fit_2, "NQ", "", Z_resolution -> GetBinCenter( Z_resolution -> GetMaximumBin() ) - (2 * Z_resolution -> GetStdDev() ), Z_resolution -> GetBinCenter( Z_resolution -> GetMaximumBin() ) + (2 * Z_resolution -> GetStdDev() ) );
1778
1779         MC_z_diff_peak = gaus_fit_2 -> GetParameter(1);
1780         MC_z_diff_width = fabs(gaus_fit_2 -> GetParameter(2));
1781     }
1782
1783     out_file -> cd();
1784     tree_out -> SetDirectory(out_file);
1785     best_fitwidth_NClus -> Write("best_fitwidth_NClus", TObject::kOverwrite);
1786     tree_out -> Write("", TObject::kOverwrite);
1787
1788     out_file -> Close();
1789 }
1790
1791 double INTTZvtx::GetZdiffPeakMC()
1792 {
1793     if (run_type == "MC" || run_type == "MC_geo_test")
1794     {
1795         if (MC_z_diff_peak != -777.)
1796         {
1797             return MC_z_diff_peak;
1798         }
1799         else
1800         {
1801             cout<<"In INTTZvtx. Are you playing with data? The MC_z_diff_peak wasn't assigned, the value is still -777. Pleak check1"<<endl;
1802             return -777.;
1803         }
1804     }
1805     else
1806     {
1807         cout<<"In INTTZvtx. Are you playing with data? The MC_z_diff_peak wasn't assigned, the value is still -777. Pleak check2"<<endl;
1808         return -777.;
1809     }
1810 }
1811
1812 double INTTZvtx::GetZdiffWidthMC()
1813 {
1814     if (run_type == "MC" || run_type == "MC_geo_test")
1815     {
1816         if (MC_z_diff_width != -777.)
1817         {
1818             return MC_z_diff_width;
1819         }
1820         else
1821         {
1822             cout<<"In INTTZvtx. Are you playing with data? The MC_z_diff_width wasn't assigned, the value is still -777. Pleak check1"<<endl;
1823             return -777.;
1824         }
1825     }
1826     else
1827     {
1828         cout<<"In INTTZvtx. Are you playing with data? The MC_z_diff_width wasn't assigned, the value is still -777. Pleak check2"<<endl;
1829         return -777.;
1830     }
1831 }
1832
1833 vector<double> INTTZvtx::GetEvtZPeak()
1834 {
1835     return{good_zvtx_tag_int, loose_offset_peak, loose_offset_peakE};
1836 }
1837
1838 double INTTZvtx::get_radius(double x, double y)
1839 {
1840     return sqrt(pow(x,2)+pow(y,2));
1841 }
1842
1843 std::vector<double> INTTZvtx::calculateDistanceAndClosestPoint(double x1, double y1, double x2, double y2, double target_x, double target_y)
1844 {
1845
1846     if (x1 != x2)
1847     {
1848         // Calculate the slope and intercept of the line passing through (x1, y1) and (x2, y2)
1849         double a = (y2 - y1) / (x2 - x1);
1850         double b = y1 - a * x1;
1851
1852         // cout<<"slope : y="<<a<<"x+"<<b<<endl;
1853
1854         // Calculate the closest distance from (target_x, target_y) to the line y = ax + b
1855         double closest_distance = std::abs(a * target_x - target_y + b) / std::sqrt(a * a + 1);
1856
1857         // Calculate the coordinates of the closest point (Xc, Yc) on the line y = ax + b
1858         double Xc = (target_x + a * target_y - a * b) / (a * a + 1);
1859         double Yc = a * Xc + b;
1860
1861         return { closest_distance, Xc, Yc };
1862     }
1863     else
1864     {
1865         double closest_distance = std::abs(x1 - target_x);
1866         double Xc = x1;
1867         double Yc = target_y;
1868
1869         return { closest_distance, Xc, Yc };
1870     }
1871
1872
1873 }
1874
1875 // note : Function to calculate the angle between two vectors in degrees using the cross product
1876 double INTTZvtx::calculateAngleBetweenVectors(double x1, double y1, double x2, double y2, double targetX, double targetY) {
1877     // Calculate the vectors vector_1 (point_1 to point_2) and vector_2 (point_1 to target)
1878     double vector1X = x2 - x1;
1879     double vector1Y = y2 - y1;
1880
1881     double vector2X = targetX - x1;
1882     double vector2Y = targetY - y1;
1883
1884     // Calculate the cross product of vector_1 and vector_2 (z-component)
1885     double crossProduct = vector1X * vector2Y - vector1Y * vector2X;
1886
1887     // cout<<" crossProduct : "<<crossProduct<<endl;
1888
1889     // Calculate the magnitudes of vector_1 and vector_2
1890     double magnitude1 = std::sqrt(vector1X * vector1X + vector1Y * vector1Y);
1891     double magnitude2 = std::sqrt(vector2X * vector2X + vector2Y * vector2Y);
1892
1893     // Calculate the angle in radians using the inverse tangent of the cross product and dot product
1894     double dotProduct = vector1X * vector2X + vector1Y * vector2Y;
1895
1896     double angleInRadians = std::atan2(std::abs(crossProduct), dotProduct);
1897     // Convert the angle from radians to degrees and return it
1898     double angleInDegrees = angleInRadians * 180.0 / M_PI;
1899
1900     double angleInRadians_new = std::asin( crossProduct/(magnitude1*magnitude2) );
1901     double angleInDegrees_new = angleInRadians_new * 180.0 / M_PI;
1902
1903     // cout<<"angle : "<<angleInDegrees_new<<endl;
1904
1905     double DCA_distance = sin(angleInRadians_new) * magnitude2;
1906
1907     return DCA_distance;
1908 }
1909
1910 void INTTZvtx::temp_bkg(TPad * c1, double peek, pair<double,double> beam_origin, InttConversion * ch_pos_DB)
1911 {
1912     c1 -> cd();
1913
1914     int N_ladder[4] = {12, 12, 16, 16};
1915     string ladder_index_string[16] = {"00","01","02","03","04","05","06","07","08","09","10","11","12","13","14","15"};
1916
1917     vector<double> x_vec; x_vec.clear();
1918     vector<double> y_vec; y_vec.clear();
1919
1920     vector<double> x_vec_2; x_vec_2.clear();
1921     vector<double> y_vec_2; y_vec_2.clear();
1922
1923     bkg = new TGraph();
1924     bkg -> SetTitle("INTT event display X-Y plane");
1925     bkg -> SetMarkerStyle(20);
1926     bkg -> SetMarkerSize(0.1);
1927     bkg -> SetPoint(0,0,0);
1928     bkg -> SetPoint(1,beam_origin.first,beam_origin.second);
1929     bkg -> GetXaxis() -> SetLimits(-150,150);
1930     bkg -> GetYaxis() -> SetRangeUser(-150,150);
1931     bkg -> GetXaxis() -> SetTitle("X [mm]");
1932     bkg -> GetYaxis() -> SetTitle("Y [mm]");
1933
1934     bkg -> Draw("ap");
1935
1936
1937
1938     for (int server_i = 0; server_i < 4; server_i++)
1939     {
1940         for (int FC_i = 0; FC_i < 14; FC_i++)
1941         {
1942             ladder_line -> DrawLine(
1943                 ch_pos_DB -> Get_XY_all(Form("intt%i",server_i),FC_i,14,0).x, ch_pos_DB -> Get_XY_all(Form("intt%i",server_i),FC_i,14,0).y,
1944                 ch_pos_DB -> Get_XY_all(Form("intt%i",server_i),FC_i,1,0).x, ch_pos_DB -> Get_XY_all(Form("intt%i",server_i),FC_i,1,0).y
1945             );
1946         }
1947     }
1948
1949     ladder_line -> Draw("l same");
1950
1951 }
1952
1953 double INTTZvtx::Get_extrapolation(double given_y, double p0x, double p0y, double p1x, double p1y) // note : x : z, y : r
1954 {
1955     if ( fabs(p0x - p1x) < 0.00001 ){ // note : the line is vertical (if z is along the x axis)
1956         return p0x;
1957     }
1958     else {
1959         double slope = (p1y - p0y) / (p1x - p0x);
1960         double yIntercept = p0y - slope * p0x;
1961         double xCoordinate = (given_y - yIntercept) / slope;
1962         return xCoordinate;
1963     }
1964 }
1965
1966 pair<double,double> INTTZvtx::Get_possible_zvtx(double rvtx, vector<double> p0, vector<double> p1) // note : inner p0, outer p1, vector {r,z}, -> {y,x}
1967 {
1968     vector<double> p0_z_edge = { ( fabs( p0[1] ) < 130 ) ? p0[1] - 8. : p0[1] - 10., ( fabs( p0[1] ) < 130 ) ? p0[1] + 8. : p0[1] + 10.}; // note : vector {left edge, right edge}
1969     vector<double> p1_z_edge = { ( fabs( p1[1] ) < 130 ) ? p1[1] - 8. : p1[1] - 10., ( fabs( p1[1] ) < 130 ) ? p1[1] + 8. : p1[1] + 10.}; // note : vector {left edge, right edge}
1970
1971     double edge_first  = Get_extrapolation(rvtx,p0_z_edge[0],p0[0],p1_z_edge[1],p1[0]);
1972     double edge_second = Get_extrapolation(rvtx,p0_z_edge[1],p0[0],p1_z_edge[0],p1[0]);
1973
1974     double mid_point = (edge_first + edge_second) / 2.;
1975     double possible_width = fabs(edge_first - edge_second) / 2.;
1976
1977     return {mid_point, possible_width}; // note : first : mid point, second : width
1978
1979 }
1980
1981 int INTTZvtx::min_width_ID(vector<double> input_width_vec)
1982 {
1983     double big_num;
1984     int big_num_id;
1985     for (int i = 0; i < input_width_vec.size(); i++)
1986     {
1987         if (i == 0) {
1988             big_num = input_width_vec[i];
1989             big_num_id = i;
1990         }
1991         else if ( input_width_vec[i] < big_num )
1992         {
1993             big_num = input_width_vec[i];
1994             big_num_id = i;
1995         }
1996     }
1997
1998     return big_num_id;
1999 }
2000
2001 void INTTZvtx::line_breakdown(TH1F* hist_in, pair<double,double> line_range)
2002 {
2003     int first_bin = int((line_range.first  - hist_in->GetXaxis()->GetXmin()) /  hist_in->GetBinWidth(1)) + 1;
2004     int last_bin  = int((line_range.second - hist_in->GetXaxis()->GetXmin()) /  hist_in->GetBinWidth(1)) + 1;
2005
2006     if (first_bin < 1) {first_bin = 0;}
2007     else if (first_bin > hist_in -> GetNbinsX()) {first_bin = hist_in -> GetNbinsX() + 1;}
2008
2009     if (last_bin < 1) {last_bin = 0;}
2010     else if (last_bin > hist_in -> GetNbinsX()) {last_bin = hist_in -> GetNbinsX() + 1;}
2011
2012     // double fill_weight = 1./fabs(line_range.second - line_range.first); // note : the entry is weitghted by the width of the line, by Akiba san's suggestion // special_tag cancel
2013     double fill_weight = 1.; // note : the weight is 1, it's for testing the trapezoidal method // special_tag
2014
2015     // cout<<"Digitize the bin : "<<first_bin<<" "<<last_bin<<endl;
2016
2017     // note : if first:last = (0:0) or (N+1:N+1) -> the subtraction of them euqals to zero.
2018     for (int i = 0; i < (last_bin - first_bin) + 1; i++){
2019         hist_in -> SetBinContent(first_bin + i,  hist_in -> GetBinContent(first_bin + i) + fill_weight ); // note : special_tag
2020         // line_breakdown_vec.push_back(hist_in -> GetBinCenter(first_bin + i));
2021     }
2022
2023 }
2024
2025 // note : for each combination
2026 // note : corrected inner_r and outer_r
2027 void INTTZvtx::trapezoidal_line_breakdown(TH1F * hist_in, double inner_r, double inner_z, double outer_r, double outer_z) // note : special_tag
2028 {
2029     combination_zvtx_range_shape -> Reset("ICESM");
2030
2031     vector<double> inner_edge = { ( fabs( inner_z ) < 130 ) ? inner_z - 8. : inner_z - 10., ( fabs( inner_z ) < 130 ) ? inner_z + 8. : inner_z + 10.}; // note : vector {left edge, right edge}
2032     vector<double> outer_edge = { ( fabs( outer_z ) < 130 ) ? outer_z - 8. : outer_z - 10., ( fabs( outer_z ) < 130 ) ? outer_z + 8. : outer_z + 10.}; // note : vector {left edge, right edge}
2033
2034     for (int possible_i = 0; possible_i < 2001; possible_i++)
2035     {
2036         // double random_inner_z = rand_gen -> Uniform(inner_edge[0], inner_edge[1]);
2037         double random_inner_z = inner_edge[0] + ((inner_edge[1] - inner_edge[0]) / 2000.) * possible_i;
2038         double edge_first  = Get_extrapolation(0, random_inner_z, inner_r, outer_edge[1], outer_r); // note : inner p0, outer p1, vector {r,z}, -> {y,x}
2039         double edge_second = Get_extrapolation(0, random_inner_z, inner_r, outer_edge[0], outer_r); // note : inner p0, outer p1, vector {r,z}, -> {y,x}
2040
2041         double mid_point = (edge_first + edge_second) / 2.;
2042         double possible_width = fabs(edge_first - edge_second) / 2.;
2043
2044         line_breakdown(combination_zvtx_range_shape, {mid_point - possible_width, mid_point + possible_width});
2045     }
2046
2047     combination_zvtx_range_shape -> Scale(1./ combination_zvtx_range_shape -> Integral(-1,-1));
2048
2049     for (int bin_i = 1; bin_i <= combination_zvtx_range_shape -> GetNbinsX(); bin_i++)
2050     {
2051         hist_in -> SetBinContent(bin_i, hist_in -> GetBinContent(bin_i) + combination_zvtx_range_shape -> GetBinContent(bin_i));
2052     }
2053 }
2054
2055 // note : search_range : should be the gaus fit range
2056 double INTTZvtx::LB_geo_mean(TH1F * hist_in, pair<double, double> search_range, int event_i)
2057 {
2058     int    Highest_bin_index   = hist_in -> GetMaximumBin();
2059     double Highest_bin_center  = hist_in -> GetBinCenter(Highest_bin_index);
2060     double Highest_bin_content = hist_in -> GetBinContent(Highest_bin_index);
2061     if (Highest_bin_center < search_range.first || search_range.second < Highest_bin_center )
2062     {
2063         // cout<<"In INTTZvtx class, event : "<<event_i<<", interesting event, different group was fit"<<endl;
2064         return -999.;
2065     }
2066
2067     vector<float> same_height_bin; same_height_bin.clear(); same_height_bin.push_back(Highest_bin_center);
2068
2069     int search_index = 1;
2070     while( (hist_in -> GetBinCenter(Highest_bin_index + search_index)) < search_range.second ){
2071         if ( hist_in -> GetBinContent(Highest_bin_index + search_index) == Highest_bin_content ) {
2072             same_height_bin.push_back( hist_in -> GetBinCenter(Highest_bin_index + search_index) );
2073         }
2074         search_index ++;
2075     }
2076
2077     // cout<<"test, search_index right : "<<search_index<<endl;
2078
2079     search_index = 1;
2080     while( search_range.first < (hist_in -> GetBinCenter(Highest_bin_index - search_index)) ){
2081         if ( hist_in -> GetBinContent(Highest_bin_index - search_index) == Highest_bin_content ) {
2082             same_height_bin.push_back( hist_in -> GetBinCenter(Highest_bin_index - search_index) );
2083         }
2084         search_index ++;
2085     }
2086
2087     // cout<<"test, search_index left : "<<search_index<<endl;
2088     // cout<<"test, same_height_bin.size(): "<<same_height_bin.size()<<endl;
2089
2090     return accumulate( same_height_bin.begin(), same_height_bin.end(), 0.0 ) / double(same_height_bin.size());
2091
2092 }
2093
2094 //  note : N group, group size, group tag, group width ?  // note : {group size, group entry, group tag, group widthL, group widthR}
2095 // note : {N_group, ratio (if two), peak widthL, peak widthR}
2096 vector<double> INTTZvtx::find_Ngroup(TH1F * hist_in)
2097 {
2098     double Highest_bin_Content  = hist_in -> GetBinContent(hist_in -> GetMaximumBin());
2099     double Highest_bin_Center   = hist_in -> GetBinCenter(hist_in -> GetMaximumBin());
2100
2101     int group_Nbin = 0;
2102     int peak_group_ID;
2103     double group_entry = 0;
2104     double peak_group_ratio;
2105     vector<int> group_Nbin_vec; group_Nbin_vec.clear();
2106     vector<double> group_entry_vec; group_entry_vec.clear();
2107     vector<double> group_widthL_vec; group_widthL_vec.clear();
2108     vector<double> group_widthR_vec; group_widthR_vec.clear();
2109
2110     for (int i = 0; i < hist_in -> GetNbinsX(); i++){
2111         // todo : the background rejection is here : Highest_bin_Content/2. for the time being
2112         double bin_content = ( hist_in -> GetBinContent(i+1) <= Highest_bin_Content/2.) ? 0. : ( hist_in -> GetBinContent(i+1) - Highest_bin_Content/2. );
2113
2114         if (bin_content != 0){
2115
2116             if (group_Nbin == 0) {
2117                 group_widthL_vec.push_back(hist_in -> GetBinCenter(i+1) - (hist_in -> GetBinWidth(i+1)/2.));
2118             }
2119
2120             group_Nbin += 1;
2121             group_entry += bin_content;
2122         }
2123         else if (bin_content == 0 && group_Nbin != 0){
2124             group_widthR_vec.push_back(hist_in -> GetBinCenter(i+1) - (hist_in -> GetBinWidth(i+1)/2.));
2125             group_Nbin_vec.push_back(group_Nbin);
2126             group_entry_vec.push_back(group_entry);
2127             group_Nbin = 0;
2128             group_entry = 0;
2129         }
2130     }
2131     if (group_Nbin != 0) {
2132         group_Nbin_vec.push_back(group_Nbin);
2133         group_entry_vec.push_back(group_entry);
2134         group_widthR_vec.push_back(hist_in -> GetXaxis()->GetXmax());
2135     } // note : the last group at the edge
2136
2137     // note : find the peak group
2138     for (int i = 0; i < group_Nbin_vec.size(); i++){
2139         if (group_widthL_vec[i] < Highest_bin_Center && Highest_bin_Center < group_widthR_vec[i]){
2140             peak_group_ID = i;
2141             break;
2142         }
2143     }
2144
2145
2146     peak_group_ratio = group_entry_vec[peak_group_ID] / (accumulate( group_entry_vec.begin(), group_entry_vec.end(), 0.0 ));
2147
2148     // for (int i = 0; i < group_Nbin_vec.size(); i++)
2149     // {
2150     //     cout<<" "<<endl;
2151     //     cout<<"group size : "<<group_Nbin_vec[i]<<endl;
2152     //     cout<<"group entry : "<<group_entry_vec[i]<<endl;
2153     //     cout<<group_widthL_vec[i]<<" "<<group_widthR_vec[i]<<endl;
2154     // }
2155
2156     // cout<<" "<<endl;
2157     // cout<<"N group : "<<group_Nbin_vec.size()<<endl;
2158     // cout<<"Peak group ID : "<<peak_group_ID<<endl;
2159     // cout<<"peak group width : "<<group_widthL_vec[peak_group_ID]<<" "<<group_widthR_vec[peak_group_ID]<<endl;
2160     // cout<<"ratio : "<<peak_group_ratio<<endl;
2161
2162     // note : {N_group, ratio (if two), peak widthL, peak widthR}
2163     return {double(group_Nbin_vec.size()), peak_group_ratio, group_widthL_vec[peak_group_ID], group_widthR_vec[peak_group_ID]};
2164 }
2165
2166 vector<double> INTTZvtx::get_half_hist_vec(TH1F * hist_in)
2167 {
2168     double Highest_bin_Content  = hist_in -> GetBinContent(hist_in -> GetMaximumBin());
2169
2170     vector<double> output_vec; output_vec.clear();
2171
2172     for (int i = 0; i < hist_in -> GetNbinsX(); i++){
2173         // todo : the background rejection is here : Highest_bin_Content/2. for the time being
2174         double bin_content = ( hist_in -> GetBinContent(i+1) <= Highest_bin_Content/2.) ? 0. : ( hist_in -> GetBinContent(i+1) - Highest_bin_Content/2. );
2175         output_vec.push_back(bin_content);
2176     }
2177
2178     return output_vec;
2179 }
2180
2181 double INTTZvtx::get_delta_phi(double angle_1, double angle_2)
2182 {
2183     vector<double> vec_abs = {fabs(angle_1 - angle_2), fabs(angle_1 - angle_2 + 360), fabs(angle_1 - angle_2 - 360)};
2184     vector<double> vec = {(angle_1 - angle_2), (angle_1 - angle_2 + 360), (angle_1 - angle_2 - 360)};
2185     return vec[std::distance(vec_abs.begin(), std::min_element(vec_abs.begin(),vec_abs.end()))];
2186 }
2187
2188 double INTTZvtx::get_track_phi(double inner_clu_phi_in, double delta_phi_in)
2189 {
2190     double track_phi = inner_clu_phi_in - (delta_phi_in/2.);
2191     if (track_phi < 0) {track_phi += 360;}
2192     else if (track_phi > 360) {track_phi -= 360;}
2193     else if (track_phi == 360) {track_phi = 0;}
2194     else {track_phi = track_phi;}
2195     return track_phi;
2196 }
2197
2198 void INTTZvtx::FakeClone1D(TH1F * hist_in, TH1F * hist_out)
2199 {
2200     if (hist_in -> GetNbinsX() != hist_out -> GetNbinsX()) {cout<<"In INTTZvtx::FakeClone1D, the input and output hist have different number of bins!!!"<<endl; return;}
2201
2202     for (int i = 0; i < hist_in -> GetNbinsX(); i++)
2203     {
2204         hist_out -> SetBinContent(i+1, hist_in -> GetBinContent(i+1));
2205     }
2206
2207     return;
2208 }
2209
2210 #endif