analysis_INTT_CW/DST_MC/INTTXYvtxEvt.h

0001 #ifndef INTTXYvtxEvt_h
0002 #define INTTXYvtxEvt_h
0003
0004 #include "INTTXYvtx.h"
0005
0006 class INTTXYvtxEvt : public INTTXYvtx{
0007
0008     public:
0009         INTTXYvtxEvt(string run_type, string out_folder_directory, pair<double,double> beam_origin, int geo_mode_id, double phi_diff_cut = 0.11, pair<double, double> DCA_cut = {-1,1}, int N_clu_cutl = 20, int N_clu_cut = 10000, bool draw_event_display = true, double peek = 3.32405, double angle_diff_new_l = 0.0, double angle_diff_new_r = 3, bool print_message_opt = true, double window_size = 2.5, int quadrant_trial = 9):
0010         INTTXYvtx(run_type, out_folder_directory, beam_origin, geo_mode_id, phi_diff_cut, DCA_cut, N_clu_cutl, N_clu_cut, draw_event_display, peek, angle_diff_new_l, angle_diff_new_r, print_message_opt), window_size(window_size), quadrant_trial(quadrant_trial)
0011         {
0012             good_tracklet_xy_vec.clear();
0013             good_tracklet_rz_vec.clear();
0014             cluster_pair_vec.clear();
0015             evt_vtx_info.clear();
0016
0017             out_reco_vtx_x.clear();
0018             out_reco_vtx_y.clear();
0019             out_reco_vtx_x_stddev.clear();
0020             out_reco_vtx_y_stddev.clear();
0021             out_binwidth_x.clear();
0022             out_binwidth_y.clear();
0023
0024             legend = new TLegend(0.2,0.87,0.64,0.92);
0025             // legend -> SetMargin(0);
0026             legend->SetTextSize(0.03);
0027             legend->SetNColumns(2);
0028
0029             inner_clu_phi_map.clear();
0030             outer_clu_phi_map.clear();
0031             inner_clu_phi_map = vector<vector<pair<bool,clu_info>>>(360);
0032             outer_clu_phi_map = vector<vector<pair<bool,clu_info>>>(360);
0033             InitTreeOut();
0034             InitHist_Evt();
0035             InitGraphEvt();
0036             if (draw_event_display == true) {c2 -> Print(Form("%s/temp_event_display.pdf(",out_folder_directory.c_str()));}
0037             if (draw_event_display == true) {c1 -> Print(Form("%s/evt_LineFill2D.pdf(",out_folder_directory.c_str()));}
0038         };
0039
0040         INTTXYvtxEvt(string run_type, string out_folder_directory):INTTXYvtx(run_type, out_folder_directory)
0041         {
0042             return;
0043         };
0044
0045         // note : the function that new created for this class
0046         virtual void ProcessEvt(
0047             int event_i,
0048             vector<clu_info> temp_sPH_inner_nocolumn_vec,
0049             vector<clu_info> temp_sPH_outer_nocolumn_vec,
0050             vector<vector<double>> temp_sPH_nocolumn_vec,
0051             vector<vector<double>> temp_sPH_nocolumn_rz_vec,
0052             int NvtxMC,
0053             vector<double> TrigvtxMC,
0054             bool PhiCheckTag,
0055             Long64_t bco_full,
0056             pair<double,double> evt_z
0057         );
0058         pair<double,double> GetVtxXYEvt();
0059         void PrintPlots_Evt();
0060         void InitTreeOut() override;
0061         void InitHist_Evt();
0062         vector<pair<double,double>> GetEvtVtxInfo();
0063         int GetReturnTag();
0064
0065         // note : the function that inherited from the INTTXYvtx class
0066         void ClearEvt() override;
0067         void EndRun() override;
0068
0069     protected:
0070         // note : these two functions are from the INTTZvtx class, and I copied them here, not so smart sigh...
0071         pair<double,double> Get_possible_zvtx(double rvtx, vector<double> p0, vector<double> p1);
0072         double Get_extrapolation(double given_y, double p0x, double p0y, double p1x, double p1y);
0073
0074         // note : functions inherited from INTTXYvtx class and would like to modify
0075         vector<pair<double,double>> MacroVTXSquare(double length, int N_trial, bool draw_plot_opt = true, vector<double> TrigvtxMC = {});
0076         void subMacroPlotWorking(bool phi_correction, double cos_fit_rangel, double cos_fit_ranger, double guas_fit_range) override;
0077
0078         // note : for the event display
0079         void temp_bkg(TPad * c1, double peek, pair<double,double> beam_origin, InttConversion * ch_pos_DB);
0080         double get_z_vertex(clu_info inner_clu, clu_info outer_clu, double target_x, double target_y);
0081         void InitGraphEvt();
0082         bool isPointInsideSquare(const std::pair<double, double> point, const std::pair<double, double> square_center, double length);
0083         void TH2LineFill(TH2F * hist_in, std::pair<double,double> point_1, std::pair<double,double> point_2);
0084         TCanvas * c2;
0085         TPad * pad_xy_hist_1;
0086         TPad * pad_xy_hist_1_bkgrm;
0087         TPad * pad_xy_hist_2;
0088         TPad * pad_xy_hist_2_bkgrm;
0089         TPad * pad_xy_hist_3;
0090         TPad * pad_xy_hist_3_bkgrm;
0091         TPad * pad_xy;
0092         TPad * pad_rz;
0093         TPad * pad_dca_inner;
0094         TPad * pad_dca_outer;
0095         TPad * pad_delta_phi_inner;
0096         TPad * pad_delta_phi_outer;
0097         TPad * pad_inner_outer_phi;
0098         TPad * pad_track_phi_1D;
0099         TPad * pad_delta_phi_1D;
0100
0101         TH2F * evt_display_xy_hist_1;
0102         TH2F * evt_display_xy_hist_1_bkgrm; // note : fine binning
0103         TH2F * evt_display_xy_hist_1_cm;
0104         TH2F * evt_display_xy_hist_1_bkgrm_cm; // note : fine binning
0105
0106         TH2F * evt_display_xy_hist_2;
0107         TH2F * evt_display_xy_hist_2_bkgrm; // note: second fine binning
0108         TH2F * evt_display_xy_hist_3;
0109         TH2F * evt_display_xy_hist_3_bkgrm; // note : coarse binning
0110         TH2F * evt_dca_inner_2D;
0111         TH2F * evt_dca_outer_2D;
0112         TH2F * evt_delta_phi_inner_2D;
0113         TH2F * evt_delta_phi_outer_2D;
0114         TH2F * evt_inner_outer_phi_2D;
0115         TH1F * evt_track_phi_1D;
0116         TH1F * evt_delta_phi_1D;
0117         TGraph * evt_display_xy_gr;
0118         TGraph * evt_display_rz_gr;
0119         TGraph * true_vertex_gr;
0120         TGraph * reco_vertex_gr;
0121
0122         TGraphErrors * vtx_evt_grE;
0123
0124         InttConversion * ch_pos_DB;
0125         TGraph * bkg;
0126         int return_tag;
0127         double reco_vtx_x;
0128         double reco_vtx_y;
0129         double reco_vtx_xE;
0130         double reco_vtx_yE;
0131         double reco_vtx_xStdDev;
0132         double reco_vtx_yStdDev;
0133
0134         int print_rate = 50;
0135
0136         vector<vector<pair<bool,clu_info>>> inner_clu_phi_map; // note: phi
0137         vector<vector<pair<bool,clu_info>>> outer_clu_phi_map; // note: phi
0138
0139         TLegend * legend;
0140
0141         TFile * file_out;
0142         TTree * tree_out;
0143         int out_eID;
0144         int out_NClus;
0145         double out_true_vtx_x;
0146         double out_true_vtx_y;
0147         double out_true_vtx_z;
0148         double out_reco_vtx_z;
0149         double out_reco_vtx_z_width;
0150         Long64_t out_bco_full;
0151         vector<double> out_reco_vtx_x;
0152         vector<double> out_reco_vtx_y;
0153         vector<double> out_reco_vtx_x_stddev;
0154         vector<double> out_reco_vtx_y_stddev;
0155         vector<double> out_binwidth_x;
0156         vector<double> out_binwidth_y;
0157
0158         // note : new created functions only for this class
0159         // void Init_Hist_Evt();
0160
0161         // note : new variables declared in this class
0162         double window_size;
0163         int quadrant_trial;
0164         vector<pair<double,double>> evt_vtx_info;
0165         vector< pair< pair<double,double>, pair<double,double> > > good_tracklet_xy_vec;
0166         vector< pair< pair<double,double>, pair<double,double> > > good_tracklet_rz_vec;
0167
0168         // note : new generated plots only for this class
0169         TH2F * VTXx_correlation;
0170         TH2F * VTXy_correlation;
0171
0172         TH1F * VTXx_diff_1D;
0173         TH1F * VTXy_diff_1D;
0174
0175         TH2F * VTXx_diff_Nclus;
0176         TH2F * VTXy_diff_Nclus;
0177
0178         TH2F * Reco_VTXxy_2D;
0179         TH1F * VTXx_1D;
0180         TH1F * VTXy_1D;
0181 };
0182
0183 void INTTXYvtxEvt::InitTreeOut()
0184 {
0185     file_out = new TFile(Form("%s/evt_XY_tree.root",out_folder_directory.c_str()),"RECREATE");
0186     file_out -> cd();
0187
0188     tree_out = new TTree("tree", "tree evt VTX xy");
0189     tree_out -> Branch("eID",&out_eID);
0190     tree_out -> Branch("NClus", &out_NClus);
0191     tree_out -> Branch("bco_full", &out_bco_full);
0192     tree_out -> Branch("true_vtx_x", &out_true_vtx_x);
0193     tree_out -> Branch("true_vtx_y", &out_true_vtx_y);
0194     tree_out -> Branch("true_vtx_z", &out_true_vtx_z);
0195     tree_out -> Branch("reco_vtx_x", &out_reco_vtx_x); // note : vector
0196     tree_out -> Branch("reco_vtx_y", &out_reco_vtx_y); // note : vector
0197     tree_out -> Branch("reco_vtx_z", &out_reco_vtx_z); // note : double
0198     tree_out -> Branch("reco_vtx_z_width", &out_reco_vtx_z_width); // note : double
0199     tree_out -> Branch("reco_vtx_x_stddev", &out_reco_vtx_x_stddev); // note : vector
0200     tree_out -> Branch("reco_vtx_y_stddev", &out_reco_vtx_y_stddev); // note : vector
0201     tree_out -> Branch("binwidth_x", &out_binwidth_x); // note : vector
0202     tree_out -> Branch("binwidth_y", &out_binwidth_y); // note : vector
0203 }
0204
0205 void INTTXYvtxEvt::InitGraphEvt()
0206 {
0207     c2 = new TCanvas("","",4000,2400);
0208     c2 -> cd();
0209
0210     pad_xy = new TPad(Form("pad_xy"), "", 0.0, 0.66, 0.2, 1.0);
0211     Characterize_Pad(pad_xy, 0.15, 0.1, 0.1, 0.2, 0, 0);
0212     pad_xy -> Draw();
0213
0214     pad_rz = new TPad(Form("pad_rz"), "", 0.2, 0.66, 0.40, 1.0);
0215     Characterize_Pad(pad_rz, 0.15, 0.1, 0.1, 0.2, 0, 0);
0216     pad_rz -> Draw();
0217
0218     pad_inner_outer_phi = new TPad(Form("pad_inner_outer_phi"), "", 0.40, 0.66, 0.6, 1.0);
0219     Characterize_Pad(pad_inner_outer_phi, 0.15, 0.1, 0.1, 0.2, 0, 0);
0220     pad_inner_outer_phi -> Draw();
0221
0222     pad_track_phi_1D = new TPad(Form("pad_track_phi_1D"), "", 0.6, 0.66, 0.8, 1.0);
0223     Characterize_Pad(pad_track_phi_1D, 0.15, 0.1, 0.1, 0.2, 0, 0);
0224     pad_track_phi_1D -> Draw();
0225
0226     pad_delta_phi_1D = new TPad(Form("pad_delta_phi_1D"), "", 0.8, 0.66, 1, 1.0);
0227     Characterize_Pad(pad_delta_phi_1D, 0.15, 0.1, 0.1, 0.2, 0, 0);
0228     pad_delta_phi_1D -> Draw();
0229
0230     pad_dca_inner = new TPad(Form("pad_dca_inner"), "", 0.0, 0.33, 0.2, 0.66);
0231     Characterize_Pad(pad_dca_inner, 0.15, 0.1, 0.1, 0.2, 0, 0);
0232     pad_dca_inner -> Draw();
0233
0234     pad_dca_outer = new TPad(Form("pad_dca_outer"), "", 0.2, 0.33, 0.40, 0.66);
0235     Characterize_Pad(pad_dca_outer, 0.15, 0.1, 0.1, 0.2, 0, 0);
0236     pad_dca_outer -> Draw();
0237
0238     pad_delta_phi_inner = new TPad(Form("pad_delta_phi_inner"), "", 0.4, 0.33, 0.60, 0.66);
0239     Characterize_Pad(pad_delta_phi_inner, 0.15, 0.1, 0.1, 0.2, 0, 0);
0240     // pad_delta_phi_inner -> SetLogz(1);
0241     pad_delta_phi_inner -> Draw();
0242
0243     pad_delta_phi_outer = new TPad(Form("pad_delta_phi_outer"), "", 0.6, 0.33, 0.80, 0.66);
0244     Characterize_Pad(pad_delta_phi_outer, 0.15, 0.1, 0.1, 0.2, 0, 0);
0245     // pad_delta_phi_outer -> SetLogz(1);
0246     pad_delta_phi_outer -> Draw();
0247
0248     pad_xy_hist_1 = new TPad(Form("pad_xy_hist_1"), "", 0.8, 0.33, 1.0, 0.66);
0249     Characterize_Pad(pad_xy_hist_1, 0.15, 0.1, 0.1, 0.2, 0, 0);
0250     pad_xy_hist_1 -> Draw();
0251
0252     pad_xy_hist_1_bkgrm = new TPad(Form("pad_xy_hist_1_bkgrm"), "", 0.0, 0.0, 0.2, 0.33);
0253     Characterize_Pad(pad_xy_hist_1_bkgrm, 0.15, 0.1, 0.1, 0.2, 0, 0);
0254     pad_xy_hist_1_bkgrm -> Draw();
0255
0256     pad_xy_hist_2 = new TPad(Form("pad_xy_hist_2"), "", 0.2, 0.0, 0.4, 0.33);
0257     Characterize_Pad(pad_xy_hist_2, 0.15, 0.1, 0.1, 0.2, 0, 0);
0258     pad_xy_hist_2 -> Draw();
0259
0260     pad_xy_hist_2_bkgrm = new TPad(Form("pad_xy_hist_2_bkgrm"), "", 0.4, 0.0, 0.6, 0.33);
0261     Characterize_Pad(pad_xy_hist_2_bkgrm, 0.15, 0.1, 0.1, 0.2, 0, 0);
0262     pad_xy_hist_2_bkgrm -> Draw();
0263
0264     pad_xy_hist_3 = new TPad(Form("pad_xy_hist_3"), "", 0.6, 0.0, 0.8, 0.33);
0265     Characterize_Pad(pad_xy_hist_3, 0.15, 0.1, 0.1, 0.2, 0, 0);
0266     pad_xy_hist_3 -> Draw();
0267
0268     pad_xy_hist_3_bkgrm = new TPad(Form("pad_xy_hist_3_bkgrm"), "", 0.8, 0.0, 1.0, 0.33);
0269     Characterize_Pad(pad_xy_hist_3_bkgrm, 0.15, 0.1, 0.1, 0.2, 0, 0);
0270     pad_xy_hist_3_bkgrm -> Draw();
0271
0272     evt_display_xy_gr = new TGraph();
0273     evt_display_rz_gr = new TGraph();
0274
0275     true_vertex_gr    = new TGraph();
0276     true_vertex_gr -> SetMarkerStyle(50);
0277     true_vertex_gr -> SetMarkerColor(2);
0278     true_vertex_gr -> SetMarkerSize(0.5);
0279
0280     reco_vertex_gr    = new TGraph();
0281     reco_vertex_gr -> SetMarkerStyle(50);
0282     reco_vertex_gr -> SetMarkerColor(4);
0283     reco_vertex_gr -> SetMarkerSize(0.5);
0284
0285     vtx_evt_grE = new TGraphErrors();
0286     vtx_evt_grE -> SetMarkerStyle(20);
0287     vtx_evt_grE -> SetMarkerSize(0.5);
0288
0289
0290     bkg = new TGraph();
0291     bkg -> SetPoint(0,0,0);
0292
0293     ch_pos_DB = new InttConversion("ideal", peek); // todo : I make it simplified here, for the setting of the geometry modes.
0294 }
0295
0296 void INTTXYvtxEvt::InitHist_Evt()
0297 {
0298     VTXx_correlation = new TH2F("","VTXx_correlation",100,-1,0.5,100,-1,0.5);
0299     VTXx_correlation -> SetStats(0);
0300     VTXx_correlation -> GetXaxis() -> SetTitle("True VTX_{x} [mm]");
0301     VTXx_correlation -> GetYaxis() -> SetTitle("Reco VTX_{x} [mm]");
0302     VTXx_correlation -> GetXaxis() -> SetNdivisions(505);
0303
0304     VTXy_correlation = new TH2F("","VTXy_correlation",100,1.7,3,100,1.7,3);
0305     VTXy_correlation -> SetStats(0);
0306     VTXy_correlation -> GetXaxis() -> SetTitle("True VTX_{y} [mm]");
0307     VTXy_correlation -> GetYaxis() -> SetTitle("Reco VTX_{y} [mm]");
0308     VTXy_correlation -> GetXaxis() -> SetNdivisions(505);
0309
0310     VTXx_diff_1D = new TH1F("","VTXx_diff_1D",100,-0.5,0.5);
0311     VTXx_diff_1D -> SetStats(0);
0312     VTXx_diff_1D -> GetXaxis() -> SetTitle("#DeltaX (Reco - True) [mm]");
0313     VTXx_diff_1D -> GetYaxis() -> SetTitle("Entry");
0314     VTXx_diff_1D -> GetXaxis() -> SetNdivisions(505);
0315
0316     VTXy_diff_1D = new TH1F("","VTXy_diff_1D",100,-0.5,0.5);
0317     VTXy_diff_1D -> SetStats(0);
0318     VTXy_diff_1D -> GetXaxis() -> SetTitle("#DeltaY (Reco - True) [mm]");
0319     VTXy_diff_1D -> GetYaxis() -> SetTitle("Entry");
0320     VTXy_diff_1D -> GetXaxis() -> SetNdivisions(505);
0321
0322     VTXx_diff_Nclus = new TH2F("","VTXx_diff_Nclus",100,0,8000,100,-1,1);
0323     VTXx_diff_Nclus -> SetStats(0);
0324     VTXx_diff_Nclus -> GetXaxis() -> SetTitle("# of clusters");
0325     VTXx_diff_Nclus -> GetYaxis() -> SetTitle("#DeltaX (Reco - True) [mm]");
0326     VTXx_diff_Nclus -> GetXaxis() -> SetNdivisions(505);
0327
0328     VTXy_diff_Nclus = new TH2F("","VTXy_diff_Nclus",100,0,8000,100,-1,1);
0329     VTXy_diff_Nclus -> SetStats(0);
0330     VTXy_diff_Nclus -> GetXaxis() -> SetTitle("# of clusters");
0331     VTXy_diff_Nclus -> GetYaxis() -> SetTitle("#DeltaY (Reco - True) [mm]");
0332     VTXy_diff_Nclus -> GetXaxis() -> SetNdivisions(505);
0333
0334
0335     // note : the followings are for the event display
0336     // note : the event display window is now moved to the nominoal beam position
0337     evt_display_xy_hist_1 = new TH2F("","evt_display_xy_hist_1", 100, -2.5 + beam_origin.first, 2.5 + beam_origin.first, 100, -2.5 + beam_origin.second, 2.5 + beam_origin.second);
0338     evt_display_xy_hist_1 -> SetStats(0);
0339     evt_display_xy_hist_1 -> GetXaxis() -> SetTitle("X axis [mm]");
0340     evt_display_xy_hist_1 -> GetYaxis() -> SetTitle("Y axis [mm]");
0341     evt_display_xy_hist_1 -> GetXaxis() -> SetNdivisions(505);
0342
0343     evt_display_xy_hist_1_bkgrm = new TH2F("","evt_display_xy_hist_1_bkgrm", 100, -2.5 + beam_origin.first, 2.5 + beam_origin.first, 100, -2.5 + beam_origin.second, 2.5 + beam_origin.second);
0344     evt_display_xy_hist_1_bkgrm -> SetStats(0);
0345     evt_display_xy_hist_1_bkgrm -> GetXaxis() -> SetTitle("X axis [mm]");
0346     evt_display_xy_hist_1_bkgrm -> GetYaxis() -> SetTitle("Y axis [mm]");
0347     evt_display_xy_hist_1_bkgrm -> GetXaxis() -> SetNdivisions(505);
0348
0349
0350
0351     evt_display_xy_hist_1_cm = new TH2F(
0352         "",
0353         "evt_display_xy_hist_1_cm;X axis [cm];Y axis [cm]",
0354         evt_display_xy_hist_1 -> GetNbinsX(),
0355         evt_display_xy_hist_1 -> GetXaxis() -> GetXmin() * 0.1,
0356         evt_display_xy_hist_1 -> GetXaxis() -> GetXmax() * 0.1,
0357         evt_display_xy_hist_1 -> GetNbinsY(),
0358         evt_display_xy_hist_1 -> GetYaxis() -> GetXmin() * 0.1,
0359         evt_display_xy_hist_1 -> GetYaxis() -> GetXmax() * 0.1
0360     );
0361     evt_display_xy_hist_1_cm -> SetStats(0);
0362     evt_display_xy_hist_1_cm -> GetXaxis() -> SetNdivisions(505);
0363
0364     evt_display_xy_hist_1_bkgrm_cm = new TH2F(
0365         "",
0366         "evt_display_xy_hist_1_bkgrm_cm;X axis [cm];Y axis [cm]",
0367         evt_display_xy_hist_1_bkgrm -> GetNbinsX(),
0368         evt_display_xy_hist_1_bkgrm -> GetXaxis() -> GetXmin() * 0.1,
0369         evt_display_xy_hist_1_bkgrm -> GetXaxis() -> GetXmax() * 0.1,
0370         evt_display_xy_hist_1_bkgrm -> GetNbinsY(),
0371         evt_display_xy_hist_1_bkgrm -> GetYaxis() -> GetXmin() * 0.1,
0372         evt_display_xy_hist_1_bkgrm -> GetYaxis() -> GetXmax() * 0.1
0373     );
0374     evt_display_xy_hist_1_bkgrm_cm -> SetStats(0);
0375     evt_display_xy_hist_1_bkgrm_cm -> GetXaxis() -> SetNdivisions(505);
0376
0377
0378
0379     evt_display_xy_hist_2 = new TH2F("","evt_display_xy_hist_2", 75, -2.5 + beam_origin.first, 2.5 + beam_origin.first, 75, -2.5 + beam_origin.second, 2.5 + beam_origin.second);
0380     evt_display_xy_hist_2 -> SetStats(0);
0381     evt_display_xy_hist_2 -> GetXaxis() -> SetTitle("X axis [mm]");
0382     evt_display_xy_hist_2 -> GetYaxis() -> SetTitle("Y axis [mm]");
0383     evt_display_xy_hist_2 -> GetXaxis() -> SetNdivisions(505);
0384
0385     evt_display_xy_hist_2_bkgrm = new TH2F("","evt_display_xy_hist_2_bkgrm", 75, -2.5 + beam_origin.first, 2.5 + beam_origin.first, 75, -2.5 + beam_origin.second, 2.5 + beam_origin.second);
0386     evt_display_xy_hist_2_bkgrm -> SetStats(0);
0387     evt_display_xy_hist_2_bkgrm -> GetXaxis() -> SetTitle("X axis [mm]");
0388     evt_display_xy_hist_2_bkgrm -> GetYaxis() -> SetTitle("Y axis [mm]");
0389     evt_display_xy_hist_2_bkgrm -> GetXaxis() -> SetNdivisions(505);
0390
0391     evt_display_xy_hist_3 = new TH2F("","evt_display_xy_hist_3", 50, -2.5 + beam_origin.first, 2.5 + beam_origin.first, 50, -2.5 + beam_origin.second, 2.5 + beam_origin.second);
0392     evt_display_xy_hist_3 -> SetStats(0);
0393     evt_display_xy_hist_3 -> GetXaxis() -> SetTitle("X axis [mm]");
0394     evt_display_xy_hist_3 -> GetYaxis() -> SetTitle("Y axis [mm]");
0395     evt_display_xy_hist_3 -> GetXaxis() -> SetNdivisions(505);
0396
0397     evt_display_xy_hist_3_bkgrm = new TH2F("","evt_display_xy_hist_3_bkgrm", 50, -2.5 + beam_origin.first, 2.5 + beam_origin.first, 50, -2.5 + beam_origin.second, 2.5 + beam_origin.second);
0398     evt_display_xy_hist_3_bkgrm -> SetStats(0);
0399     evt_display_xy_hist_3_bkgrm -> GetXaxis() -> SetTitle("X axis [mm]");
0400     evt_display_xy_hist_3_bkgrm -> GetYaxis() -> SetTitle("Y axis [mm]");
0401     evt_display_xy_hist_3_bkgrm -> GetXaxis() -> SetNdivisions(505);
0402
0403     VTXx_1D = new TH1F("","VTXx_1D", 100, -2 + beam_origin.first, 2 + beam_origin.first);
0404     VTXx_1D -> SetStats(0);
0405     VTXx_1D -> GetXaxis() -> SetTitle("X axis [mm]");
0406     VTXx_1D -> GetYaxis() -> SetTitle("Entry");
0407     VTXx_1D -> GetXaxis() -> SetNdivisions(505);
0408
0409     VTXy_1D = new TH1F("","VTXy_1D", 100, -2 + beam_origin.second, 2 + beam_origin.second);
0410     VTXy_1D -> SetStats(0);
0411     VTXy_1D -> GetXaxis() -> SetTitle("Y axis [mm]");
0412     VTXy_1D -> GetYaxis() -> SetTitle("Entry");
0413     VTXy_1D -> GetXaxis() -> SetNdivisions(505);
0414
0415     Reco_VTXxy_2D = new TH2F("","Reco_VTXxy_2D", 100, -2 + beam_origin.first, 2 + beam_origin.first, 100, -2 + beam_origin.second, 2 + beam_origin.second);
0416     // Reco_VTXxy_2D -> SetStats(0);
0417     Reco_VTXxy_2D -> GetXaxis() -> SetTitle("X axis [mm]");
0418     Reco_VTXxy_2D -> GetYaxis() -> SetTitle("Y axis [mm]");
0419     // Reco_VTXxy_2D -> GetXaxis() -> SetNdivisions(505);
0420
0421     evt_dca_inner_2D = new TH2F("","evt_dca_inner_2D",72,0,360,20,-10,10);
0422     evt_dca_inner_2D -> SetStats(0);
0423     evt_dca_inner_2D -> GetXaxis() -> SetTitle("Inner cluster #phi [degree]");
0424     evt_dca_inner_2D -> GetYaxis() -> SetTitle("DCA [mm]");
0425     evt_dca_inner_2D -> GetXaxis() -> SetNdivisions(505);
0426
0427     evt_dca_outer_2D = new TH2F("","evt_dca_outer_2D",72,0,360,20,-10,10);
0428     evt_dca_outer_2D -> SetStats(0);
0429     evt_dca_outer_2D -> GetXaxis() -> SetTitle("Outer cluster #phi [degree]");
0430     evt_dca_outer_2D -> GetYaxis() -> SetTitle("DCA [mm]");
0431     evt_dca_outer_2D -> GetXaxis() -> SetNdivisions(505);
0432
0433     evt_delta_phi_inner_2D = new TH2F("","evt_delta_phi_inner_2D",72,0,360,50,-1.5,1.5);
0434     evt_delta_phi_inner_2D -> SetStats(0);
0435     evt_delta_phi_inner_2D -> GetXaxis() -> SetTitle("Inner cluster #phi [degree]");
0436     evt_delta_phi_inner_2D -> GetYaxis() -> SetTitle("#Delta#phi (inner - outer) [degree]");
0437     evt_delta_phi_inner_2D -> GetXaxis() -> SetNdivisions(505);
0438
0439     evt_delta_phi_outer_2D = new TH2F("","evt_delta_phi_outer_2D",72,0,360,50,-1.5,1.5);
0440     evt_delta_phi_outer_2D -> SetStats(0);
0441     evt_delta_phi_outer_2D -> GetXaxis() -> SetTitle("Outer cluster #phi [degree]");
0442     evt_delta_phi_outer_2D -> GetYaxis() -> SetTitle("#Delta#phi (inner - outer) [degree]");
0443     evt_delta_phi_outer_2D -> GetXaxis() -> SetNdivisions(505);
0444
0445     evt_inner_outer_phi_2D = new TH2F("","evt_inner_outer_phi_2D",72,0,360,72,0,360);
0446     evt_inner_outer_phi_2D -> SetStats(0);
0447     evt_inner_outer_phi_2D -> GetXaxis() -> SetTitle("Inner cluster #phi [degree]");
0448     evt_inner_outer_phi_2D -> GetYaxis() -> SetTitle("Outer cluster #phi [degree]");
0449     evt_inner_outer_phi_2D -> GetXaxis() -> SetNdivisions(505);
0450
0451     evt_track_phi_1D = new TH1F("","evt_track_phi_1D",180,0,360);
0452     evt_track_phi_1D -> SetStats(0);
0453     evt_track_phi_1D -> GetXaxis() -> SetTitle("Tracklet #phi [degree]");
0454     evt_track_phi_1D -> GetYaxis() -> SetTitle("Entry / (2 degree)");
0455     evt_track_phi_1D -> GetXaxis() -> SetNdivisions(505);
0456
0457     evt_delta_phi_1D = new TH1F("","evt_delta_phi_1D",100,-10,10);
0458     evt_delta_phi_1D -> SetStats(0);
0459     evt_delta_phi_1D -> GetXaxis() -> SetTitle("#Delta#phi (inner - outer) [degree]");
0460     evt_delta_phi_1D -> GetYaxis() -> SetTitle("Entry / (0.2 mm)");
0461     evt_delta_phi_1D -> GetXaxis() -> SetNdivisions(505);
0462 }
0463
0464 void INTTXYvtxEvt::PrintPlots_Evt()
0465 {
0466     c1 -> cd();
0467
0468     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
0469     VTXx_correlation -> Draw("colz0");
0470     VTXx_correlation -> Fit(correlation_fit_MC, "NQ");
0471     correlation_fit_MC -> Draw("l same");
0472     draw_text -> DrawLatex(0.21, 0.75, Form("INTT NClus > %d", N_clu_cutl));
0473     draw_text -> DrawLatex(0.21, 0.71, Form("y = %.4f * X +  %.4f",correlation_fit_MC -> GetParameter(1), correlation_fit_MC -> GetParameter(0)));
0474     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
0475     c1 -> Print(Form("%s/VTXx_correlation.pdf",out_folder_directory.c_str()));
0476     c1 -> Clear();
0477
0478     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
0479     VTXy_correlation -> Draw("colz0");
0480     VTXy_correlation -> Fit(correlation_fit_MC, "NQ");
0481     correlation_fit_MC -> Draw("l same");
0482     draw_text -> DrawLatex(0.21, 0.75, Form("INTT NClus > %d", N_clu_cutl));
0483     draw_text -> DrawLatex(0.21, 0.71, Form("y = %.4f * X +  %.4f",correlation_fit_MC -> GetParameter(1), correlation_fit_MC -> GetParameter(0)));
0484     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
0485     c1 -> Print(Form("%s/VTXy_correlation.pdf",out_folder_directory.c_str()));
0486     c1 -> Clear();
0487
0488     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
0489     VTXx_diff_1D -> Draw("hist");
0490     gaus_fit_MC -> SetParameters(VTXx_diff_1D -> GetBinContent( VTXx_diff_1D -> GetMaximumBin() ), VTXx_diff_1D -> GetBinCenter( VTXx_diff_1D -> GetMaximumBin() ), 0.05, 0);
0491     gaus_fit_MC -> SetParLimits(0,0,100000);  // note : size
0492     gaus_fit_MC -> SetParLimits(2,0,10000);   // note : Width
0493     gaus_fit_MC -> SetParLimits(3,0,10000);   // note : offset
0494     VTXx_diff_1D -> Fit(gaus_fit_MC, "NQ", "", VTXx_diff_1D -> GetBinCenter( VTXx_diff_1D -> GetMaximumBin() ) - (2 * VTXx_diff_1D -> GetStdDev() ), VTXx_diff_1D -> GetBinCenter( VTXx_diff_1D -> GetMaximumBin() ) + (2 * VTXx_diff_1D -> GetStdDev() ) );
0495     gaus_fit_MC -> SetRange( gaus_fit_MC->GetParameter(1) - gaus_fit_MC->GetParameter(2) * 2.5, gaus_fit_MC->GetParameter(1) + gaus_fit_MC->GetParameter(2) * 2.5 );
0496     gaus_fit_MC -> Draw("lsame");
0497     draw_text -> DrawLatex(0.21, 0.75, Form("INTT NClus > %d", N_clu_cutl));
0498     draw_text -> DrawLatex(0.21, 0.71, Form("Gaus mean  : %.4f mm",gaus_fit_MC -> GetParameter(1)));
0499     draw_text -> DrawLatex(0.21, 0.67, Form("Gaus width : %.4f mm",fabs(gaus_fit_MC -> GetParameter(2))));
0500     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
0501     c1 -> Print(Form("%s/VTXx_diff_1D.pdf",out_folder_directory.c_str()));
0502     c1 -> Clear();
0503
0504     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
0505     VTXy_diff_1D -> Draw("hist");
0506     gaus_fit_MC -> SetParameters(VTXy_diff_1D -> GetBinContent( VTXy_diff_1D -> GetMaximumBin() ), VTXy_diff_1D -> GetBinCenter( VTXy_diff_1D -> GetMaximumBin() ), 0.05, 0);
0507     gaus_fit_MC -> SetParLimits(0,0,100000);  // note : size
0508     gaus_fit_MC -> SetParLimits(2,0,10000);   // note : Width
0509     gaus_fit_MC -> SetParLimits(3,0,10000);   // note : offset
0510     VTXy_diff_1D -> Fit(gaus_fit_MC, "NQ", "", VTXy_diff_1D -> GetBinCenter( VTXy_diff_1D -> GetMaximumBin() ) - (2 * VTXy_diff_1D -> GetStdDev() ), VTXy_diff_1D -> GetBinCenter( VTXy_diff_1D -> GetMaximumBin() ) + (2 * VTXy_diff_1D -> GetStdDev() ) );
0511     gaus_fit_MC -> SetRange( gaus_fit_MC->GetParameter(1) - gaus_fit_MC->GetParameter(2) * 2.5, gaus_fit_MC->GetParameter(1) + gaus_fit_MC->GetParameter(2) * 2.5 );
0512     gaus_fit_MC -> Draw("lsame");
0513     draw_text -> DrawLatex(0.21, 0.75, Form("INTT NClus > %d", N_clu_cutl));
0514     draw_text -> DrawLatex(0.21, 0.71, Form("Gaus mean  : %.4f mm",gaus_fit_MC -> GetParameter(1)));
0515     draw_text -> DrawLatex(0.21, 0.67, Form("Gaus width : %.4f mm",fabs(gaus_fit_MC -> GetParameter(2))));
0516     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
0517     c1 -> Print(Form("%s/VTXy_diff_1D.pdf",out_folder_directory.c_str()));
0518     c1 -> Clear();
0519
0520     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
0521     VTXx_diff_Nclus -> Draw("colz0");
0522     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
0523     c1 -> Print(Form("%s/VTXx_diff_Nclus.pdf",out_folder_directory.c_str()));
0524     c1 -> Clear();
0525
0526     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
0527     VTXy_diff_Nclus -> Draw("colz0");
0528     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
0529     c1 -> Print(Form("%s/VTXy_diff_Nclus.pdf",out_folder_directory.c_str()));
0530     c1 -> Clear();
0531
0532     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
0533     Reco_VTXxy_2D -> Draw("colz0");
0534     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
0535     draw_text -> DrawLatex(0.21, 0.88, Form("Entries: %.0f", Reco_VTXxy_2D -> GetEntries()));
0536     c1 -> Print(Form("%s/Reco_VTXxy_2D.pdf",out_folder_directory.c_str()));
0537     c1 -> Clear();
0538
0539     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
0540     c1 -> cd();
0541     VTXx_1D -> Draw("hist");
0542     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
0543
0544     gaus_fit_MC -> SetParameters(VTXx_1D -> GetBinContent( VTXx_1D -> GetMaximumBin() ), VTXx_1D -> GetBinCenter( VTXx_1D -> GetMaximumBin() ), 0.05, 0);
0545     gaus_fit_MC -> SetParLimits(0,0,100000);  // note : size
0546     gaus_fit_MC -> SetParLimits(2,0,10000);   // note : Width
0547     gaus_fit_MC -> SetParLimits(3,0,10000);   // note : offset
0548     VTXx_1D -> Fit(gaus_fit_MC, "N", "", VTXx_1D -> GetBinCenter( VTXx_1D -> GetMaximumBin() ) - (2 * VTXx_1D -> GetStdDev() ), VTXx_1D -> GetBinCenter( VTXx_1D -> GetMaximumBin() ) + (2 * VTXx_1D -> GetStdDev() ) );
0549     gaus_fit_MC -> SetRange( gaus_fit_MC->GetParameter(1) - double(gaus_fit_MC->GetParameter(2)) * 2.5, gaus_fit_MC->GetParameter(1) + double(gaus_fit_MC->GetParameter(2)) * 2.5 );
0550     gaus_fit_MC -> Draw("lsame");
0551
0552     draw_text -> DrawLatex(0.21, 0.71, Form("Gaus mean  : %.4f mm",gaus_fit_MC -> GetParameter(1)));
0553     draw_text -> DrawLatex(0.21, 0.67, Form("Gaus width : %.4f mm",fabs(gaus_fit_MC -> GetParameter(2))));
0554     c1 -> Print(Form("%s/VTXx_1D.pdf",out_folder_directory.c_str()));
0555     c1 -> Clear();
0556
0557     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
0558     c1 -> cd();
0559     VTXy_1D -> Draw("hist");
0560     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
0561
0562     gaus_fit_MC -> SetParameters(VTXy_1D -> GetBinContent( VTXy_1D -> GetMaximumBin() ), VTXy_1D -> GetBinCenter( VTXy_1D -> GetMaximumBin() ), 0.05, 0);
0563     gaus_fit_MC -> SetParLimits(0,0,100000);  // note : size
0564     gaus_fit_MC -> SetParLimits(2,0,10000);   // note : Width
0565     gaus_fit_MC -> SetParLimits(3,0,10000);   // note : offset
0566     VTXy_1D -> Fit(gaus_fit_MC, "N", "", VTXy_1D -> GetBinCenter( VTXy_1D -> GetMaximumBin() ) - (2 * VTXy_1D -> GetStdDev() ), VTXy_1D -> GetBinCenter( VTXy_1D -> GetMaximumBin() ) + (2 * VTXy_1D -> GetStdDev() ) );
0567     gaus_fit_MC -> SetRange( gaus_fit_MC->GetParameter(1) - double(gaus_fit_MC->GetParameter(2)) * 2.5, gaus_fit_MC->GetParameter(1) + double(gaus_fit_MC->GetParameter(2)) * 2.5 );
0568     gaus_fit_MC -> Draw("lsame");
0569
0570     draw_text -> DrawLatex(0.21, 0.71, Form("Gaus mean  : %.4f mm",gaus_fit_MC -> GetParameter(1)));
0571     draw_text -> DrawLatex(0.21, 0.67, Form("Gaus width : %.4f mm",fabs(gaus_fit_MC -> GetParameter(2))));
0572     c1 -> Print(Form("%s/VTXy_1D.pdf",out_folder_directory.c_str()));
0573     c1 -> Clear();
0574
0575 }
0576
0577 // note : this function is for the event by event vertex calculation
0578 void INTTXYvtxEvt::ProcessEvt(int event_i, vector<clu_info> temp_sPH_inner_nocolumn_vec, vector<clu_info> temp_sPH_outer_nocolumn_vec, vector<vector<double>> temp_sPH_nocolumn_vec, vector<vector<double>> temp_sPH_nocolumn_rz_vec, int NvtxMC, vector<double> TrigvtxMC, bool PhiCheckTag, Long64_t bco_full, pair<double,double> evt_z) // note : evt_z : {z, width}
0579 {
0580     return_tag = 0;
0581
0582     if (event_i%1 == 0) {cout<<"In INTTXYvtxEvt class, running event : "<<event_i<<endl;}
0583
0584     total_NClus = temp_sPH_inner_nocolumn_vec.size() + temp_sPH_outer_nocolumn_vec.size();
0585
0586
0587     if (total_NClus < zvtx_cal_require) {cout<<"low total_Nclus, "<<total_NClus<<endl; return; }
0588
0589     if (run_type == "MC" && NvtxMC != 1) { cout<<"In INTTXYvtxEvt class, event : "<<event_i<<" Nvtx : "<<NvtxMC<<" Nvtx more than one "<<endl; return; }
0590     if (PhiCheckTag == false)            { cout<<"In INTTXYvtxEvt class, event : "<<event_i<<" Nvtx : "<<NvtxMC<<" Not full phi has hits "<<endl; return; }
0591
0592     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)
0593     {
0594         printf("In INTTXYvtxEvt class, event : %i, low clu continue, NClus : %lu \n", event_i, total_NClus);
0595         return;
0596     }
0597
0598     // cout<<"test_1"<<endl;
0599     // note : put the cluster into the phi map, the first bool is for the cluster usage.
0600     // note : false means the cluster is not used
0601     for ( int inner_i = 0; inner_i < temp_sPH_inner_nocolumn_vec.size(); inner_i++ ){
0602         evt_display_xy_gr -> SetPoint(evt_display_xy_gr -> GetN(), temp_sPH_inner_nocolumn_vec[inner_i].x, temp_sPH_inner_nocolumn_vec[inner_i].y);
0603         evt_display_rz_gr -> SetPoint(evt_display_rz_gr -> GetN(), temp_sPH_inner_nocolumn_vec[inner_i].z, (temp_sPH_inner_nocolumn_vec[inner_i].phi > 180) ? get_radius(temp_sPH_inner_nocolumn_vec[inner_i].x,temp_sPH_inner_nocolumn_vec[inner_i].y) * -1 : get_radius(temp_sPH_inner_nocolumn_vec[inner_i].x,temp_sPH_inner_nocolumn_vec[inner_i].y));
0604
0605         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());
0606         inner_clu_phi_map[ int(Clus_InnerPhi_Offset) ].push_back({false,temp_sPH_inner_nocolumn_vec[inner_i]});
0607     }
0608     // cout<<"test_2"<<endl;
0609     for ( int outer_i = 0; outer_i < temp_sPH_outer_nocolumn_vec.size(); outer_i++ ){
0610         evt_display_xy_gr -> SetPoint(evt_display_xy_gr -> GetN(), temp_sPH_outer_nocolumn_vec[outer_i].x, temp_sPH_outer_nocolumn_vec[outer_i].y);
0611         evt_display_rz_gr -> SetPoint(evt_display_rz_gr -> GetN(), temp_sPH_outer_nocolumn_vec[outer_i].z, (temp_sPH_outer_nocolumn_vec[outer_i].phi > 180) ? get_radius(temp_sPH_outer_nocolumn_vec[outer_i].x,temp_sPH_outer_nocolumn_vec[outer_i].y) * -1 : get_radius(temp_sPH_outer_nocolumn_vec[outer_i].x,temp_sPH_outer_nocolumn_vec[outer_i].y));
0612
0613         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());
0614         outer_clu_phi_map[ int(Clus_OuterPhi_Offset) ].push_back({false,temp_sPH_outer_nocolumn_vec[outer_i]});
0615     }
0616
0617     for (int inner_phi_i = 0; inner_phi_i < 360; inner_phi_i++) // note : each phi cell (1 degree)
0618     {
0619         // note : N cluster in this phi cell
0620         for (int inner_phi_clu_i = 0; inner_phi_clu_i < inner_clu_phi_map[inner_phi_i].size(); inner_phi_clu_i++)
0621         {
0622             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());
0623
0624             // todo: change the outer phi scan range
0625             // note : the outer phi index, -2, -1, 0, 1, 2
0626             for (int scan_i = -2; scan_i < 3; scan_i++)
0627             {
0628                 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;
0629
0630                 // note : N clusters in that outer phi cell
0631                 for (int outer_phi_clu_i = 0; outer_phi_clu_i < outer_clu_phi_map[true_scan_i].size(); outer_phi_clu_i++)
0632                 {
0633
0634                     // cout<<event_i<<" test_5 "<<inner_i<<" "<<outer_i<<endl;
0635                     // note : calculate the cluster phi after the vertex correction which can enhence the purity of the tracklet selection
0636                     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());
0637
0638                     // note : before calculating the possible z vertex range of one tracklet, the vertex correction was applied
0639                     pair<double,double> z_range_info = Get_possible_zvtx(
0640                         0., // get_radius(beam_origin.first,beam_origin.second),
0641                         {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
0642                         {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
0643                     );
0644                     // cout<<event_i<<" test_6 "<<inner_i<<" "<<outer_i<<endl;
0645                     // note : this is a cut to constraint on the z vertex, only if the tracklets with the range that covers the z vertex can pass this cut
0646                     if (z_range_info.first - z_range_info.second > evt_z.first + evt_z.second || z_range_info.first + z_range_info.second < evt_z.first - evt_z.second) {continue;}
0647                     // cout<<"range test: "<<z_range_info.first - z_range_info.second<<" "<<z_range_info.first + z_range_info.second<<" vertex: "<<evt_z.first<<" "<<evt_z.second<<endl;
0648
0649                     // note : the reason to use the DCA calculation with sign is because that the distribution of 1D DCA will be single peak only
0650                     double DCA_sign = calculateAngleBetweenVectors(
0651                         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,
0652                         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,
0653                         beam_origin.first, beam_origin.second
0654                     );
0655
0656                     vector<double> DCA_info_vec = calculateDistanceAndClosestPoint(
0657                         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,
0658                         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,
0659                         beam_origin.first, beam_origin.second
0660                     );
0661
0662                     if (DCA_info_vec[0] != fabs(DCA_sign) && fabs( DCA_info_vec[0] - fabs(DCA_sign) ) > 0.1 ){
0663                         cout<<"different DCA : "<<DCA_info_vec[0]<<" "<<DCA_sign<<" diff : "<<DCA_info_vec[0] - fabs(DCA_sign)<<endl;
0664                     }
0665                     // cout<<event_i<<" test_7 "<<inner_i<<" "<<outer_i<<endl;
0666                     evt_dca_inner_2D -> Fill(inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.phi, DCA_sign);
0667                     evt_dca_outer_2D -> Fill(outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.phi, DCA_sign);
0668                     evt_delta_phi_inner_2D -> Fill(inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.phi, 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);
0669                     evt_delta_phi_outer_2D -> Fill(outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.phi, 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);
0670                     evt_inner_outer_phi_2D -> Fill(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);
0671                     evt_delta_phi_1D -> Fill( 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 );
0672
0673                     // if (fabs(Clus_InnerPhi_Offset - Clus_OuterPhi_Offset) < phi_diff_cut)
0674                     if (true) // note : use the cut based on the cluster map, which means, -2 to 2 degrees
0675                     {
0676                         // if (DCA_cut.first < DCA_sign && DCA_sign < DCA_cut.second)
0677                         if (true)
0678                         {
0679
0680                             // // note : before calculating the possible z vertex range of one tracklet, the vertex correction was applied
0681                             // pair<double,double> z_range_info = Get_possible_zvtx(
0682                             //     0., // get_radius(beam_origin.first,beam_origin.second),
0683                             //     {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
0684                             //     {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
0685                             // );
0686
0687                             // // note : this is a cut to constraint on the z vertex, only if the tracklets with the range that covers the z vertex can pass this cut
0688                             // if (z_range_info.first - z_range_info.second > evt_z.first + evt_z.second || z_range_info.first + z_range_info.second < evt_z.first - evt_z.second) {continue;}
0689                             // // cout<<"range test: "<<z_range_info.first - z_range_info.second<<" "<<z_range_info.first + z_range_info.second<<" vertex: "<<evt_z.first<<" "<<evt_z.second<<endl;
0690
0691                             // cout<<event_i<<" test_8 "<<inner_i<<" "<<outer_i<<endl;
0692
0693                             // note : for the tracklets that pass all the cut listed above can then be used for the vertex XY calculation
0694                             cluster_pair_vec.push_back({{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},{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}});
0695                             // cout<<"we have something here "<<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<<endl;
0696
0697                             good_tracklet_xy_vec.push_back({
0698                                 {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},
0699                                 {DCA_info_vec[1], DCA_info_vec[2]}
0700                             });
0701
0702                             good_tracklet_rz_vec.push_back({
0703                                 {outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.z, (outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.phi > 180) ? get_radius(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) * -1 : get_radius(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)},
0704                                 {get_z_vertex(inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second,outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second,DCA_info_vec[1],DCA_info_vec[2]), (outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.phi > 180) ? get_radius(DCA_info_vec[1],DCA_info_vec[2]) * -1 : get_radius(DCA_info_vec[1],DCA_info_vec[2])}
0705                             });
0706
0707                             // cout<<event_i<<" test_9 "<<inner_i<<" "<<outer_i<<endl;
0708
0709                             TH2FSampleLineFill(evt_display_xy_hist_1, 0.001, {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}, {DCA_info_vec[1], DCA_info_vec[2]});
0710                             // TH2FSampleLineFill(evt_display_xy_hist_2, 0.001, {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}, {DCA_info_vec[1], DCA_info_vec[2]});
0711                             // TH2FSampleLineFill(evt_display_xy_hist_3, 0.001, {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}, {DCA_info_vec[1], DCA_info_vec[2]});
0712
0713                             // TH2LineFill(evt_display_xy_hist_1, {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}, {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});
0714                             evt_track_phi_1D -> Fill((Clus_InnerPhi_Offset + Clus_OuterPhi_Offset)/2.);
0715
0716                             // cout<<event_i<<" test_10 "<<inner_i<<" "<<outer_i<<endl;
0717                         }
0718
0719                     }
0720
0721
0722
0723
0724
0725
0726
0727
0728                 }
0729             }
0730         }
0731     }
0732
0733     // cout<<"test_3"<<endl;
0734     // for ( int inner_i = 0; inner_i < temp_sPH_inner_nocolumn_vec.size(); inner_i++ )
0735     // {
0736     //     // cout<<"test_4"<<endl;
0737     //     for ( int outer_i = 0; outer_i < temp_sPH_outer_nocolumn_vec.size(); outer_i++ )
0738     //     {
0739     //         // cout<<event_i<<" test_5 "<<inner_i<<" "<<outer_i<<endl;
0740     //         // note : calculate the cluster phi after the vertex correction which can enhence the purity of the tracklet selection
0741     //         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());
0742     //         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());
0743
0744     //         // note : before calculating the possible z vertex range of one tracklet, the vertex correction was applied
0745     //         pair<double,double> z_range_info = Get_possible_zvtx(
0746     //             0., // get_radius(beam_origin.first,beam_origin.second),
0747     //             {get_radius(temp_sPH_inner_nocolumn_vec[inner_i].x - beam_origin.first, temp_sPH_inner_nocolumn_vec[inner_i].y - beam_origin.second), temp_sPH_inner_nocolumn_vec[inner_i].z}, // note : unsign radius
0748     //             {get_radius(temp_sPH_outer_nocolumn_vec[outer_i].x - beam_origin.first, temp_sPH_outer_nocolumn_vec[outer_i].y - beam_origin.second), temp_sPH_outer_nocolumn_vec[outer_i].z}  // note : unsign radius
0749     //         );
0750     //         // cout<<event_i<<" test_6 "<<inner_i<<" "<<outer_i<<endl;
0751     //         // note : this is a cut to constraint on the z vertex, only if the tracklets with the range that covers the z vertex can pass this cut
0752     //         if (z_range_info.first - z_range_info.second > evt_z.first + evt_z.second || z_range_info.first + z_range_info.second < evt_z.first - evt_z.second) {continue;}
0753     //         // cout<<"range test: "<<z_range_info.first - z_range_info.second<<" "<<z_range_info.first + z_range_info.second<<" vertex: "<<evt_z.first<<" "<<evt_z.second<<endl;
0754
0755     //         // note : the reason to use the DCA calculation with sign is because that the distribution of 1D DCA will be single peak only
0756     //         double DCA_sign = calculateAngleBetweenVectors(
0757     //             temp_sPH_outer_nocolumn_vec[outer_i].x, temp_sPH_outer_nocolumn_vec[outer_i].y,
0758     //             temp_sPH_inner_nocolumn_vec[inner_i].x, temp_sPH_inner_nocolumn_vec[inner_i].y,
0759     //             beam_origin.first, beam_origin.second
0760     //         );
0761
0762     //         vector<double> DCA_info_vec = calculateDistanceAndClosestPoint(
0763     //             temp_sPH_inner_nocolumn_vec[inner_i].x, temp_sPH_inner_nocolumn_vec[inner_i].y,
0764     //             temp_sPH_outer_nocolumn_vec[outer_i].x, temp_sPH_outer_nocolumn_vec[outer_i].y,
0765     //             beam_origin.first, beam_origin.second
0766     //         );
0767
0768     //         if (DCA_info_vec[0] != fabs(DCA_sign) && fabs( DCA_info_vec[0] - fabs(DCA_sign) ) > 0.1 ){
0769     //             cout<<"different DCA : "<<DCA_info_vec[0]<<" "<<DCA_sign<<" diff : "<<DCA_info_vec[0] - fabs(DCA_sign)<<endl;
0770     //         }
0771     //         // cout<<event_i<<" test_7 "<<inner_i<<" "<<outer_i<<endl;
0772     //         evt_dca_inner_2D -> Fill(temp_sPH_inner_nocolumn_vec[inner_i].phi, DCA_sign);
0773     //         evt_dca_outer_2D -> Fill(temp_sPH_outer_nocolumn_vec[outer_i].phi, DCA_sign);
0774     //         evt_delta_phi_inner_2D -> Fill(temp_sPH_inner_nocolumn_vec[inner_i].phi, temp_sPH_inner_nocolumn_vec[inner_i].phi - temp_sPH_outer_nocolumn_vec[outer_i].phi);
0775     //         evt_delta_phi_outer_2D -> Fill(temp_sPH_outer_nocolumn_vec[outer_i].phi, temp_sPH_inner_nocolumn_vec[inner_i].phi - temp_sPH_outer_nocolumn_vec[outer_i].phi);
0776     //         evt_inner_outer_phi_2D -> Fill(temp_sPH_inner_nocolumn_vec[inner_i].phi, temp_sPH_outer_nocolumn_vec[outer_i].phi);
0777     //         evt_delta_phi_1D -> Fill( temp_sPH_inner_nocolumn_vec[inner_i].phi - temp_sPH_outer_nocolumn_vec[outer_i].phi );
0778
0779     //         if (fabs(Clus_InnerPhi_Offset - Clus_OuterPhi_Offset) < phi_diff_cut)
0780     //         {
0781     //             if (DCA_cut.first < DCA_sign && DCA_sign < DCA_cut.second){
0782
0783     //                 // // note : before calculating the possible z vertex range of one tracklet, the vertex correction was applied
0784     //                 // pair<double,double> z_range_info = Get_possible_zvtx(
0785     //                 //     0., // get_radius(beam_origin.first,beam_origin.second),
0786     //                 //     {get_radius(temp_sPH_inner_nocolumn_vec[inner_i].x - beam_origin.first, temp_sPH_inner_nocolumn_vec[inner_i].y - beam_origin.second), temp_sPH_inner_nocolumn_vec[inner_i].z}, // note : unsign radius
0787     //                 //     {get_radius(temp_sPH_outer_nocolumn_vec[outer_i].x - beam_origin.first, temp_sPH_outer_nocolumn_vec[outer_i].y - beam_origin.second), temp_sPH_outer_nocolumn_vec[outer_i].z}  // note : unsign radius
0788     //                 // );
0789
0790     //                 // // note : this is a cut to constraint on the z vertex, only if the tracklets with the range that covers the z vertex can pass this cut
0791     //                 // if (z_range_info.first - z_range_info.second > evt_z.first + evt_z.second || z_range_info.first + z_range_info.second < evt_z.first - evt_z.second) {continue;}
0792     //                 // // cout<<"range test: "<<z_range_info.first - z_range_info.second<<" "<<z_range_info.first + z_range_info.second<<" vertex: "<<evt_z.first<<" "<<evt_z.second<<endl;
0793
0794     //                 // cout<<event_i<<" test_8 "<<inner_i<<" "<<outer_i<<endl;
0795
0796     //                 // note : for the tracklets that pass all the cut listed above can then be used for the vertex XY calculation
0797     //                 cluster_pair_vec.push_back({{temp_sPH_inner_nocolumn_vec[inner_i].x,temp_sPH_inner_nocolumn_vec[inner_i].y},{temp_sPH_outer_nocolumn_vec[outer_i].x,temp_sPH_outer_nocolumn_vec[outer_i].y}});
0798     //                 // cout<<"we have something here "<<temp_sPH_inner_nocolumn_vec[inner_i].x<<" "<<temp_sPH_inner_nocolumn_vec[inner_i].y<<endl;
0799
0800     //                 good_tracklet_xy_vec.push_back({
0801     //                     {temp_sPH_outer_nocolumn_vec[outer_i].x, temp_sPH_outer_nocolumn_vec[outer_i].y},
0802     //                     {DCA_info_vec[1], DCA_info_vec[2]}
0803     //                 });
0804
0805     //                 good_tracklet_rz_vec.push_back({
0806     //                     {temp_sPH_outer_nocolumn_vec[outer_i].z, (temp_sPH_outer_nocolumn_vec[outer_i].phi > 180) ? get_radius(temp_sPH_outer_nocolumn_vec[outer_i].x,temp_sPH_outer_nocolumn_vec[outer_i].y) * -1 : get_radius(temp_sPH_outer_nocolumn_vec[outer_i].x,temp_sPH_outer_nocolumn_vec[outer_i].y)},
0807     //                     {get_z_vertex(temp_sPH_inner_nocolumn_vec[inner_i],temp_sPH_outer_nocolumn_vec[outer_i],DCA_info_vec[1],DCA_info_vec[2]), (temp_sPH_outer_nocolumn_vec[outer_i].phi > 180) ? get_radius(DCA_info_vec[1],DCA_info_vec[2]) * -1 : get_radius(DCA_info_vec[1],DCA_info_vec[2])}
0808     //                 });
0809
0810     //                 // cout<<event_i<<" test_9 "<<inner_i<<" "<<outer_i<<endl;
0811
0812     //                 TH2FSampleLineFill(evt_display_xy_hist_1, 0.001, {temp_sPH_inner_nocolumn_vec[inner_i].x, temp_sPH_inner_nocolumn_vec[inner_i].y}, {DCA_info_vec[1], DCA_info_vec[2]});
0813     //                 // TH2LineFill(evt_display_xy_hist_1, {temp_sPH_inner_nocolumn_vec[inner_i].x, temp_sPH_inner_nocolumn_vec[inner_i].y}, {temp_sPH_outer_nocolumn_vec[outer_i].x, temp_sPH_outer_nocolumn_vec[outer_i].y});
0814     //                 evt_track_phi_1D -> Fill((temp_sPH_inner_nocolumn_vec[inner_i].phi + temp_sPH_outer_nocolumn_vec[outer_i].phi)/2.);
0815
0816     //                 // cout<<event_i<<" test_10 "<<inner_i<<" "<<outer_i<<endl;
0817     //             }
0818
0819     //         }
0820
0821     //     }
0822     // }
0823
0824     // note : try to remove the background
0825     TH2F_FakeClone(evt_display_xy_hist_1,evt_display_xy_hist_1_bkgrm);
0826     TH2F_threshold_advanced_2(evt_display_xy_hist_1_bkgrm, 0.7);
0827
0828     // note : for the presentation, mm -> cm
0829     TH2F_FakeClone(evt_display_xy_hist_1, evt_display_xy_hist_1_cm);
0830     TH2F_FakeClone(evt_display_xy_hist_1_bkgrm, evt_display_xy_hist_1_bkgrm_cm);
0831
0832     // TH2F_FakeClone(evt_display_xy_hist_2,evt_display_xy_hist_2_bkgrm);
0833     // TH2F_threshold_advanced_2(evt_display_xy_hist_2_bkgrm, 0.7);
0834
0835     // TH2F_FakeClone(evt_display_xy_hist_3,evt_display_xy_hist_3_bkgrm);
0836     // TH2F_threshold_advanced_2(evt_display_xy_hist_3_bkgrm, 0.7);
0837
0838     // cout<<"--- In process event: "<<event_i<<" "<<total_NClus<<" "<<cluster_pair_vec.size()<<endl;
0839
0840     // evt_vtx_info = MacroVTXSquare(window_size, quadrant_trial, false, TrigvtxMC);
0841     // evt_vtx_info = {{1,2},{1,2}};
0842     // double reco_vtx_x = evt_vtx_info[0].first;
0843     // double reco_vtx_y = evt_vtx_info[0].second;
0844
0845     // delete DCA_distance_inner_phi_peak_final;
0846     // delete angle_diff_inner_phi_peak_final;
0847
0848
0849     reco_vtx_x = evt_display_xy_hist_1_bkgrm->GetMean(1); // + (evt_display_xy_hist_1_bkgrm -> GetXaxis() -> GetBinWidth(1) / 2.); // note : the TH2F calculate the GetMean based on the bin center, no need to apply additional offset
0850     reco_vtx_y = evt_display_xy_hist_1_bkgrm->GetMean(2); // + (evt_display_xy_hist_1_bkgrm -> GetYaxis() -> GetBinWidth(1) / 2.); // note : the TH2F calculate the GetMean based on the bin center, no need to apply additional offset
0851     reco_vtx_xE = evt_display_xy_hist_1_bkgrm->GetMeanError(1);
0852     reco_vtx_yE = evt_display_xy_hist_1_bkgrm->GetMeanError(2);
0853     reco_vtx_xStdDev = evt_display_xy_hist_1_bkgrm->GetStdDev(1);
0854     reco_vtx_yStdDev = evt_display_xy_hist_1_bkgrm->GetStdDev(2);
0855
0856
0857
0858     // note : prepare for the tree
0859     out_eID = event_i;
0860     out_NClus = total_NClus;
0861     out_bco_full = bco_full;
0862     out_true_vtx_x = TrigvtxMC[0]*10.;
0863     out_true_vtx_y = TrigvtxMC[1]*10.;
0864     out_true_vtx_z = TrigvtxMC[2]*10.;
0865     out_reco_vtx_z = evt_z.first;
0866     out_reco_vtx_z_width = evt_z.second;
0867     out_reco_vtx_x = {evt_display_xy_hist_1_bkgrm->GetMean(1)}; // , evt_display_xy_hist_2_bkgrm->GetMean(1), evt_display_xy_hist_3_bkgrm->GetMean(1)};
0868     out_reco_vtx_y = {evt_display_xy_hist_1_bkgrm->GetMean(2)}; // , evt_display_xy_hist_2_bkgrm->GetMean(2), evt_display_xy_hist_3_bkgrm->GetMean(2)};
0869     out_reco_vtx_x_stddev = {evt_display_xy_hist_1_bkgrm->GetStdDev(1)}; // , evt_display_xy_hist_2_bkgrm->GetStdDev(1), evt_display_xy_hist_3_bkgrm->GetStdDev(1)};
0870     out_reco_vtx_y_stddev = {evt_display_xy_hist_1_bkgrm->GetStdDev(2)}; // , evt_display_xy_hist_2_bkgrm->GetStdDev(2), evt_display_xy_hist_3_bkgrm->GetStdDev(2)};
0871     out_binwidth_x = {evt_display_xy_hist_1_bkgrm->GetXaxis()->GetBinWidth(1)}; // , evt_display_xy_hist_2_bkgrm->GetXaxis()->GetBinWidth(1), evt_display_xy_hist_3_bkgrm->GetXaxis()->GetBinWidth(1)};
0872     out_binwidth_y = {evt_display_xy_hist_1_bkgrm->GetYaxis()->GetBinWidth(1)}; // , evt_display_xy_hist_2_bkgrm->GetYaxis()->GetBinWidth(1), evt_display_xy_hist_3_bkgrm->GetYaxis()->GetBinWidth(1)};
0873
0874     tree_out -> Fill();
0875
0876     // todo : The high multiplicity cut for the performance cut is here
0877     if (total_NClus > 2500)
0878     {
0879         VTXx_correlation -> Fill(TrigvtxMC[0]*10., reco_vtx_x);
0880         VTXy_correlation -> Fill(TrigvtxMC[1]*10., reco_vtx_y);
0881         VTXx_diff_1D -> Fill(reco_vtx_x - TrigvtxMC[0]*10.);
0882         VTXy_diff_1D -> Fill(reco_vtx_y - TrigvtxMC[1]*10.);
0883
0884         Reco_VTXxy_2D -> Fill(reco_vtx_x, reco_vtx_y);
0885         VTXx_1D -> Fill(reco_vtx_x);
0886         VTXy_1D -> Fill(reco_vtx_y);
0887     }
0888
0889     VTXx_diff_Nclus -> Fill(total_NClus, reco_vtx_x - TrigvtxMC[0]*10.);
0890     VTXy_diff_Nclus -> Fill(total_NClus, reco_vtx_y - TrigvtxMC[1]*10.);
0891
0892     cout<<"note, event: "<<event_i<<"----------------------- ----------------------- ----------------------- ----------------------- -----------------------"<<endl;
0893     cout<<"reco vtxXY by reco zvtx: "<<reco_vtx_x<<" "<<reco_vtx_y<<endl;
0894     if (run_type == "MC")
0895     {
0896         cout<<"true vtxXY : "<<TrigvtxMC[0]*10.<<" "<<TrigvtxMC[1]*10.<<endl;
0897         cout<<"Deviation "<<reco_vtx_x - TrigvtxMC[0]*10.<<" "<<reco_vtx_y - TrigvtxMC[1]*10.<<endl;
0898         cout<<"Distance "<<sqrt(pow(reco_vtx_x - TrigvtxMC[0]*10.,2)+pow(reco_vtx_y - TrigvtxMC[1]*10.,2))<<endl;
0899     }
0900
0901
0902     return_tag = 1;
0903
0904     if (draw_event_display == true && (event_i % print_rate) == 0)
0905     {
0906         c2 -> cd();
0907         // evt_display_xy_gr = new TGraph(temp_sPH_nocolumn_vec[0].size(),&temp_sPH_nocolumn_vec[0][0],&temp_sPH_nocolumn_vec[1][0]);
0908         evt_display_xy_gr -> SetTitle("INTT event display X-Y plane");
0909         evt_display_xy_gr -> GetXaxis() -> SetLimits(-150,150);
0910         evt_display_xy_gr -> GetYaxis() -> SetRangeUser(-150,150);
0911         evt_display_xy_gr -> GetXaxis() -> SetTitle("X [mm]");
0912         evt_display_xy_gr -> GetYaxis() -> SetTitle("Y [mm]");
0913         evt_display_xy_gr -> SetMarkerStyle(20);
0914         evt_display_xy_gr -> SetMarkerColor(2);
0915         evt_display_xy_gr -> SetMarkerSize(1);
0916
0917         // evt_display_rz_gr = new TGraph(temp_sPH_nocolumn_rz_vec[0].size(),&temp_sPH_nocolumn_rz_vec[0][0],&temp_sPH_nocolumn_rz_vec[1][0]);
0918         evt_display_rz_gr -> SetTitle("INTT event display r-Z plane");
0919         evt_display_rz_gr -> GetXaxis() -> SetLimits(-500,500);
0920         evt_display_rz_gr -> GetYaxis() -> SetRangeUser(-150,150);
0921         evt_display_rz_gr -> GetXaxis() -> SetTitle("Z [mm]");
0922         evt_display_rz_gr -> GetYaxis() -> SetTitle("Radius [mm]");
0923         evt_display_rz_gr -> SetMarkerStyle(20);
0924         evt_display_rz_gr -> SetMarkerColor(2);
0925         evt_display_rz_gr -> SetMarkerSize(1);
0926
0927         // note : ----------------------- ----------------------- ----------------------- ----------------------- -----------------------
0928         pad_xy -> cd();
0929         temp_bkg(pad_xy, peek, beam_origin, ch_pos_DB);
0930         evt_display_xy_gr -> Draw("p same");
0931         for (int track_i = 0; track_i < good_tracklet_xy_vec.size(); track_i++){
0932             track_line -> DrawLine(
0933                 good_tracklet_xy_vec[track_i].first.first, good_tracklet_xy_vec[track_i].first.second,
0934                 good_tracklet_xy_vec[track_i].second.first, good_tracklet_xy_vec[track_i].second.second
0935             );
0936         }
0937         // track_line -> Draw("l same");
0938         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()));
0939
0940         // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
0941         pad_rz -> cd();
0942         evt_display_rz_gr -> Draw("ap");
0943         // eff_sig_range_line -> DrawLine(temp_event_zvtx_info[0],-150,temp_event_zvtx_info[0],150);
0944         coord_line -> DrawLine(0,-150,0,150);
0945         coord_line -> DrawLine(-500,0,500,0);
0946         for (int track_i = 0; track_i < good_tracklet_rz_vec.size(); track_i++){
0947             track_line -> DrawLine(
0948                 good_tracklet_rz_vec[track_i].first.first, good_tracklet_rz_vec[track_i].first.second,
0949                 good_tracklet_rz_vec[track_i].second.first, good_tracklet_rz_vec[track_i].second.second
0950             );
0951         }
0952         draw_text -> DrawLatex(0.2, 0.85, Form("Negative radius : Clu_{outer} > 180^{0}"));
0953
0954         // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
0955         pad_dca_inner -> cd();
0956         evt_dca_inner_2D -> Draw("colz0");
0957
0958         // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
0959         pad_dca_outer -> cd();
0960         evt_dca_outer_2D -> Draw("colz0");
0961
0962         // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
0963         pad_delta_phi_inner -> cd();
0964         evt_delta_phi_inner_2D -> Draw("colz0");
0965
0966         // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
0967         pad_delta_phi_outer -> cd();
0968         evt_delta_phi_outer_2D -> Draw("colz0");
0969
0970         // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
0971         pad_inner_outer_phi -> cd();
0972         evt_inner_outer_phi_2D -> Draw("colz0");
0973
0974         // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
0975         pad_track_phi_1D -> cd();
0976         evt_track_phi_1D -> Draw("hist");
0977
0978         // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
0979         pad_delta_phi_1D -> cd();
0980         evt_delta_phi_1D -> Draw("hist");
0981
0982         // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
0983         pad_xy_hist_1 -> cd();
0984         evt_display_xy_hist_1 -> Draw("colz0");
0985         true_vertex_gr -> SetPoint(true_vertex_gr->GetN(), TrigvtxMC[0] * 10, TrigvtxMC[1] * 10);
0986         true_vertex_gr -> Draw("psame");
0987         // cout<<event_i<<" test_1 "<<endl;
0988         // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
0989         pad_xy_hist_1_bkgrm -> cd();
0990         evt_display_xy_hist_1_bkgrm -> Draw("colz0");
0991         reco_vertex_gr -> SetPoint(reco_vertex_gr->GetN(), reco_vtx_x , reco_vtx_y);
0992         true_vertex_gr -> Draw("psame");
0993         reco_vertex_gr -> Draw("psame");
0994
0995         // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
0996         pad_xy_hist_2 -> cd();
0997         evt_display_xy_hist_2 -> Draw("colz0");
0998         true_vertex_gr -> Draw("psame");
0999         // cout<<event_i<<" test_1 "<<endl;
1000         // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1001         pad_xy_hist_2_bkgrm -> cd();
1002         evt_display_xy_hist_2_bkgrm -> Draw("colz0");
1003         true_vertex_gr -> Draw("psame");
1004
1005         // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1006         pad_xy_hist_3 -> cd();
1007         evt_display_xy_hist_1 -> Draw("colz0");
1008         true_vertex_gr -> Draw("psame");
1009         // cout<<event_i<<" test_1 "<<endl;
1010         // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1011         pad_xy_hist_3_bkgrm -> cd();
1012         evt_display_xy_hist_1_bkgrm -> Draw("colz0");
1013         true_vertex_gr -> Draw("psame");
1014
1015         // cout<<event_i<<" test_2 "<<endl;
1016
1017         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()));}
1018
1019         pad_xy -> Clear();
1020         pad_rz -> Clear();
1021         pad_xy_hist_1 -> Clear();
1022         pad_xy_hist_1_bkgrm -> Clear();
1023         pad_xy_hist_2 -> Clear();
1024         pad_xy_hist_2_bkgrm -> Clear();
1025         pad_xy_hist_3 -> Clear();
1026         pad_xy_hist_3_bkgrm -> Clear();
1027         pad_dca_inner -> Clear();
1028         pad_dca_outer -> Clear();
1029         pad_delta_phi_inner -> Clear();
1030         pad_delta_phi_outer -> Clear();
1031         pad_inner_outer_phi -> Clear();
1032         pad_track_phi_1D -> Clear();
1033         pad_delta_phi_1D -> Clear();
1034
1035         // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1036
1037         if (true_vertex_gr -> GetN() != 0) {true_vertex_gr -> Set(0);}
1038         if (reco_vertex_gr -> GetN() != 0) {reco_vertex_gr -> Set(0);}
1039
1040         true_vertex_gr -> SetPoint(true_vertex_gr -> GetN(), TrigvtxMC[0], TrigvtxMC[1]);
1041         reco_vertex_gr -> SetPoint(reco_vertex_gr -> GetN(), reco_vtx_x * 0.1 , reco_vtx_y * 0.1);
1042
1043         true_vertex_gr -> SetMarkerStyle(50);
1044         true_vertex_gr -> SetMarkerColor(2);
1045         true_vertex_gr -> SetMarkerSize(1.);
1046         reco_vertex_gr -> SetMarkerStyle(50);
1047         reco_vertex_gr -> SetMarkerColor(4);
1048         reco_vertex_gr -> SetMarkerSize(1.);
1049
1050         c1 -> cd();
1051         evt_display_xy_hist_1_cm -> Draw("colz0");
1052         ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1053         draw_text -> DrawLatex(0.21, 0.82, Form("event ID: %i", event_i));
1054         if (draw_event_display && (event_i % print_rate) == 0) {c1 -> Print(Form("%s/evt_LineFill2D.pdf",out_folder_directory.c_str()));}
1055         c1 -> Clear();
1056
1057
1058         legend -> AddEntry(reco_vertex_gr, "Reco. vertex XY", "p");
1059
1060         evt_display_xy_hist_1_bkgrm_cm -> Draw("colz0");
1061         ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1062         draw_text -> DrawLatex(0.21, 0.82, Form("event ID: %i", event_i));
1063         draw_text -> DrawLatex(0.21, 0.78, Form("Reco. vertex XY: %.4f cm, %.4f cm", reco_vertex_gr->GetPointX(0), reco_vertex_gr->GetPointY(0)));
1064         reco_vertex_gr -> Draw("psame");
1065         if (run_type == "MC"){
1066             legend -> AddEntry(true_vertex_gr, "True vertex XY", "p");
1067             true_vertex_gr -> Draw("psame");
1068             draw_text -> DrawLatex(0.21, 0.74, Form("True vertex XY: %.4f cm, %.4f cm", true_vertex_gr->GetPointX(0), true_vertex_gr->GetPointY(0)));
1069         }
1070         legend -> Draw("same");
1071         if (draw_event_display && (event_i % print_rate) == 0) {c1 -> Print(Form("%s/evt_LineFill2D.pdf",out_folder_directory.c_str()));}
1072         c1 -> Clear();
1073         legend -> Clear();
1074
1075         if (true_vertex_gr -> GetN() != 0) {true_vertex_gr -> Set(0);}
1076         if (reco_vertex_gr -> GetN() != 0) {reco_vertex_gr -> Set(0);}
1077     }
1078 }
1079
1080 pair<double,double> INTTXYvtxEvt::GetVtxXYEvt() { return evt_vtx_info[0]; }
1081 vector<pair<double,double>> INTTXYvtxEvt::GetEvtVtxInfo() { return {{reco_vtx_x, reco_vtx_y}, {reco_vtx_xE, reco_vtx_yE}, {reco_vtx_xStdDev, reco_vtx_yStdDev}}; }
1082 int INTTXYvtxEvt::GetReturnTag() { return return_tag; }
1083
1084 void INTTXYvtxEvt::ClearEvt()
1085 {
1086     good_tracklet_xy_vec.clear();
1087     good_tracklet_rz_vec.clear();
1088     cluster_pair_vec.clear();
1089     if (evt_display_xy_gr -> GetN() != 0) {evt_display_xy_gr -> Set(0);}
1090     if (evt_display_rz_gr -> GetN() != 0) {evt_display_rz_gr -> Set(0);}
1091
1092     evt_display_xy_hist_1 -> Reset("ICESM");
1093     evt_display_xy_hist_2 -> Reset("ICESM");
1094     evt_display_xy_hist_3 -> Reset("ICESM");
1095     evt_dca_inner_2D -> Reset("ICESM");
1096     evt_dca_outer_2D -> Reset("ICESM");
1097     evt_delta_phi_inner_2D -> Reset("ICESM");
1098     evt_delta_phi_outer_2D -> Reset("ICESM");
1099     evt_inner_outer_phi_2D -> Reset("ICESM");
1100     evt_track_phi_1D -> Reset("ICESM");
1101     evt_delta_phi_1D -> Reset("ICESM");
1102     evt_display_xy_hist_1_bkgrm -> Reset("ICESM");
1103     evt_display_xy_hist_2_bkgrm -> Reset("ICESM");
1104     evt_display_xy_hist_3_bkgrm -> Reset("ICESM");
1105
1106     out_reco_vtx_x.clear();
1107     out_reco_vtx_y.clear();
1108     out_reco_vtx_x_stddev.clear();
1109     out_reco_vtx_y_stddev.clear();
1110     out_binwidth_x.clear();
1111     out_binwidth_y.clear();
1112
1113     inner_clu_phi_map.clear();
1114     outer_clu_phi_map.clear();
1115     inner_clu_phi_map = vector<vector<pair<bool,clu_info>>>(360);
1116     outer_clu_phi_map = vector<vector<pair<bool,clu_info>>>(360);
1117
1118     if (true_vertex_gr -> GetN() != 0) {true_vertex_gr -> Set(0);}
1119     if (reco_vertex_gr -> GetN() != 0) {reco_vertex_gr -> Set(0);}
1120
1121     return;
1122 }
1123
1124 vector<pair<double,double>> INTTXYvtxEvt::MacroVTXSquare(double length, int N_trial, bool draw_plot_opt, vector<double> TrigvtxMC)
1125 {
1126     double original_length = length;
1127     pair<double,double> origin = {0,0};
1128     // pair<double,double> origin = {beam_origin.first,beam_origin.second};
1129
1130     vector<pair<double,double>> vtx_vec = Get4vtx(origin,length); // vtx_vec.push_back(origin);
1131     int small_index;
1132     vector<double> small_info_vec;
1133
1134     vector<pair<double,double>> vtx_vec_axis = Get4vtxAxis(origin,length);
1135     int small_index_axis;
1136     vector<double> small_info_vec_axis[vtx_vec_axis.size()];
1137
1138     vector<double> grr_x; grr_x.clear();
1139     vector<double> grr_E; grr_E.clear();
1140     vector<double> grr_y; grr_y.clear();
1141
1142     if (print_message_opt == true) {cout<<N_trial<<" runs, smart. which gives you the resolution down to "<<length/pow(2,N_trial)<<" mm"<<endl;}
1143
1144     for (int i = 0; i < N_trial; i++)
1145     {
1146         if (print_message_opt == true) {cout<<"~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~~~~"<<" step "<<i<<" ~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~~~~"<<endl;}
1147         for (int i1 = 0; i1 < vtx_vec.size(); i1++)
1148         {
1149             // if (print_message_opt == true) {cout<<"testd vertex : "<<vtx_vec[i1].first<<" "<<vtx_vec[i1].second<<endl;}
1150             current_vtxX = vtx_vec[i1].first;
1151             current_vtxY = vtx_vec[i1].second;
1152             vector<double> info_vec = subMacroVTXxyCorrection(i,i1, draw_plot_opt);
1153             if (print_message_opt == true) {cout<<"DCA fit E: "<<info_vec[3]<<" delta phi fit E: "<<info_vec[5]<<endl;}
1154             if (i1 == 0)
1155             {
1156                 small_info_vec = info_vec;
1157                 small_index = i1;
1158
1159                 // DCA_distance_inner_phi_peak_final = (TH2F*)DCA_distance_inner_phi_peak -> Clone();
1160                 // angle_diff_inner_phi_peak_final = (TH2F*)angle_diff_inner_phi_peak -> Clone();
1161             }
1162             else
1163             {
1164                 // if (info_vec[3] < small_info_vec[3] && info_vec[5] < small_info_vec[5]) // note : the chi2 of the pol0 fit
1165                 if (info_vec[3] + info_vec[5] * 10 < small_info_vec[3] + small_info_vec[5] * 10)
1166                 {
1167                     small_info_vec = info_vec;
1168                     small_index = i1;
1169
1170                     // DCA_distance_inner_phi_peak_final -> Delete();
1171                     // angle_diff_inner_phi_peak_final -> Delete();
1172                     // delete DCA_distance_inner_phi_peak_final;
1173                     // delete angle_diff_inner_phi_peak_final;
1174                     // DCA_distance_inner_phi_peak_final = (TH2F*)DCA_distance_inner_phi_peak -> Clone();
1175                     // angle_diff_inner_phi_peak_final = (TH2F*)angle_diff_inner_phi_peak -> Clone();
1176                 }
1177             }
1178             if (print_message_opt == true){cout<<" "<<endl;}
1179
1180             ClearHist(1);
1181         }
1182
1183         // if (print_message_opt == true) {cout<<"the Quadrant "<<small_index<<" won the competition with the vertex: "<<vtx_vec[small_index].first<<" "<<vtx_vec[small_index].second<<endl;}
1184
1185         grr_x.push_back(i);
1186         grr_y.push_back(small_index);
1187         grr_E.push_back(0);
1188
1189         // note : the for loop for the case of the verteix on the axis
1190         for (int i1 = 0; i1 < vtx_vec_axis.size(); i1++)
1191         {
1192             current_vtxX = vtx_vec_axis[i1].first;
1193             current_vtxY = vtx_vec_axis[i1].second;
1194             small_info_vec_axis[i1] = subMacroVTXxyCorrection(i,i1, draw_plot_opt);
1195
1196             ClearHist(1);
1197         }
1198
1199         int winner_quadrant_axis = axis_quadrant_map[{
1200             (small_info_vec_axis[0][3] + small_info_vec_axis[0][5] * 10 < small_info_vec_axis[1][3] + small_info_vec_axis[1][5] * 10) ? 1 : -1,
1201             (small_info_vec_axis[2][3] + small_info_vec_axis[2][5] * 10 < small_info_vec_axis[3][3] + small_info_vec_axis[3][5] * 10) ? 1 : -1
1202         }];
1203
1204         if (true == true)
1205         {
1206             // cout<<" "<<endl;
1207             cout<<"------------------------------------------------------------------------------------"<<endl;
1208             cout<<"for axis method: "<<endl;
1209             cout<<"inner DCA, r-chi2: "<<small_info_vec_axis[0][2]<<" fit error: "<<small_info_vec_axis[0][3]<<endl;
1210             cout<<"inner delta Phi, r-chi2: "<<small_info_vec_axis[0][4]<<" fit error: "<<small_info_vec_axis[0][5]<<endl;
1211             cout<<"outer DCA, r-chi2: "<<small_info_vec_axis[0][6]<<" fit error: "<<small_info_vec_axis[0][7]<<endl;
1212             cout<<"outer delta Phi, r-chi2: "<<small_info_vec_axis[0][8]<<" fit error: "<<small_info_vec_axis[0][9]<<endl;
1213             cout<<"~~~~      ~~~~      ~~~~      ~~~~      ~~~~      ~~~~      ~~~~      ~~~~"<<endl;
1214             cout<<"inner DCA, r-chi2: "<<small_info_vec_axis[1][2]<<" fit error: "<<small_info_vec_axis[1][3]<<endl;
1215             cout<<"inner delta Phi, r-chi2: "<<small_info_vec_axis[1][4]<<" fit error: "<<small_info_vec_axis[1][5]<<endl;
1216             cout<<"outer DCA, r-chi2: "<<small_info_vec_axis[1][6]<<" fit error: "<<small_info_vec_axis[1][7]<<endl;
1217             cout<<"outer delta Phi, r-chi2: "<<small_info_vec_axis[1][8]<<" fit error: "<<small_info_vec_axis[1][9]<<endl;
1218             cout<<"~~~~      ~~~~      ~~~~      ~~~~      ~~~~      ~~~~      ~~~~      ~~~~"<<endl;
1219             cout<<"inner DCA, r-chi2: "<<small_info_vec_axis[2][2]<<" fit error: "<<small_info_vec_axis[2][3]<<endl;
1220             cout<<"inner delta Phi, r-chi2: "<<small_info_vec_axis[2][4]<<" fit error: "<<small_info_vec_axis[2][5]<<endl;
1221             cout<<"outer DCA, r-chi2: "<<small_info_vec_axis[2][6]<<" fit error: "<<small_info_vec_axis[2][7]<<endl;
1222             cout<<"outer delta Phi, r-chi2: "<<small_info_vec_axis[2][8]<<" fit error: "<<small_info_vec_axis[2][9]<<endl;
1223             cout<<"~~~~      ~~~~      ~~~~      ~~~~      ~~~~      ~~~~      ~~~~      ~~~~"<<endl;
1224             cout<<"inner DCA, r-chi2: "<<small_info_vec_axis[3][2]<<" fit error: "<<small_info_vec_axis[3][3]<<endl;
1225             cout<<"inner delta Phi, r-chi2: "<<small_info_vec_axis[3][4]<<" fit error: "<<small_info_vec_axis[3][5]<<endl;
1226             cout<<"outer DCA, r-chi2: "<<small_info_vec_axis[3][6]<<" fit error: "<<small_info_vec_axis[3][7]<<endl;
1227             cout<<"outer delta Phi, r-chi2: "<<small_info_vec_axis[3][8]<<" fit error: "<<small_info_vec_axis[3][9]<<endl;
1228             cout<<" "<<endl;
1229             cout<<"the trial origin: "<<origin.first<<" "<<origin.second<<" length: "<<length<<endl;
1230             cout<<"the true vertex: "<<TrigvtxMC[0]*10<<" "<<TrigvtxMC[1]*10<<endl;
1231             cout<<"the Quadrant "<<small_index<<" won the competition with the vertex: "<<vtx_vec[small_index].first<<" "<<vtx_vec[small_index].second<<endl;
1232             cout<<"the quadrant, true: "<<find_quadrant({origin},{TrigvtxMC[0]*10, TrigvtxMC[1]*10})<<" the axis method: "<<winner_quadrant_axis<<" the quadrant method: "<<small_index<<endl;
1233             cout<<"------------------------------------------------------------------------------------"<<endl;
1234         }
1235
1236         // if ( winner_quadrant_axis != small_index) {
1237         //     cout<<" In N trail: "<<i<<", the two methods give different results, "<<"the Quadrant, Axis_method:  "<<winner_quadrant_axis<<" cornor_method: "<<small_index<<endl;
1238         //     cout<<" The true answer should be : "<<find_quadrant({origin},{TrigvtxMC[0]*10, TrigvtxMC[1]*10})<<endl;
1239         // }
1240
1241         // note : generating the new 5 vertex for the next comparison
1242         // note : start to shrink the square
1243         if (i != N_trial - 1)
1244         {
1245             origin = {(vtx_vec[small_index].first + origin.first)/2., (vtx_vec[small_index].second + origin.second)/2.};
1246             // cout<<"test : "<<origin.first<<" "<<origin.second<<" length: "<<length<<endl;
1247             length /= 2.;
1248             vtx_vec = Get4vtx(origin,length); // vtx_vec.push_back(origin);
1249             vtx_vec_axis = Get4vtxAxis(origin,length);
1250         }
1251     }
1252
1253     if (draw_plot_opt == true) {DrawTGraphErrors(grr_x, grr_y, grr_E, grr_E, out_folder_directory, {Form("Square_scan_history_%.1fmm_%iTrials", original_length, N_trial), "nth scan", "Winner index", "APL"});}
1254
1255     // note : final vtx winner, the origin in that test, fit error info.,        fit parameter info.
1256     return {vtx_vec[small_index],origin, {small_info_vec[3],small_info_vec[5]}, {small_info_vec[10],small_info_vec[11]}};
1257     // return {{1,2},{1,2}};
1258 }
1259
1260 void INTTXYvtxEvt::subMacroPlotWorking(bool phi_correction, double cos_fit_rangel, double cos_fit_ranger, double guas_fit_range)
1261 {
1262     for (int i = 0; i < cluster_pair_vec.size(); i++)
1263     {
1264         vector<double> DCA_info_vec = calculateDistanceAndClosestPoint(
1265             cluster_pair_vec[i].first.x, cluster_pair_vec[i].first.y,
1266             cluster_pair_vec[i].second.x, cluster_pair_vec[i].second.y,
1267             current_vtxX, current_vtxY
1268         );
1269
1270         double DCA_sign = calculateAngleBetweenVectors(
1271             cluster_pair_vec[i].second.x, cluster_pair_vec[i].second.y,
1272             cluster_pair_vec[i].first.x, cluster_pair_vec[i].first.y,
1273             current_vtxX, current_vtxY
1274         );
1275
1276         if (phi_correction == true)
1277         {
1278             // cout<<"option selected "<<endl;
1279             Clus_InnerPhi_Offset = (cluster_pair_vec[i].first.y - current_vtxY < 0) ? atan2(cluster_pair_vec[i].first.y - current_vtxY, cluster_pair_vec[i].first.x - current_vtxX) * (180./TMath::Pi()) + 360 : atan2(cluster_pair_vec[i].first.y - current_vtxY, cluster_pair_vec[i].first.x - current_vtxX) * (180./TMath::Pi());
1280             Clus_OuterPhi_Offset = (cluster_pair_vec[i].second.y - current_vtxY < 0) ? atan2(cluster_pair_vec[i].second.y - current_vtxY, cluster_pair_vec[i].second.x - current_vtxX) * (180./TMath::Pi()) + 360 : atan2(cluster_pair_vec[i].second.y - current_vtxY, cluster_pair_vec[i].second.x - current_vtxX) * (180./TMath::Pi());
1281         }
1282         else // note : phi_correction == false
1283         {
1284             Clus_InnerPhi_Offset = (cluster_pair_vec[i].first.y < 0) ? atan2(cluster_pair_vec[i].first.y, cluster_pair_vec[i].first.x) * (180./TMath::Pi()) + 360 : atan2(cluster_pair_vec[i].first.y, cluster_pair_vec[i].first.x) * (180./TMath::Pi());
1285             Clus_OuterPhi_Offset = (cluster_pair_vec[i].second.y < 0) ? atan2(cluster_pair_vec[i].second.y, cluster_pair_vec[i].second.x) * (180./TMath::Pi()) + 360 : atan2(cluster_pair_vec[i].second.y, cluster_pair_vec[i].second.x) * (180./TMath::Pi());
1286         }
1287
1288         // double phi_test = (cluster_pair_vec[i].first.y < 0) ? atan2(cluster_pair_vec[i].first.y, cluster_pair_vec[i].first.x) * (180./TMath::Pi()) + 360 : atan2(cluster_pair_vec[i].first.y, cluster_pair_vec[i].first.x) * (180./TMath::Pi());
1289         // double phi_test_offset = (cluster_pair_vec[i].first.y - current_vtxY < 0) ? atan2(cluster_pair_vec[i].first.y - current_vtxY, cluster_pair_vec[i].first.x - current_vtxX) * (180./TMath::Pi()) + 360 : atan2(cluster_pair_vec[i].first.y - current_vtxY, cluster_pair_vec[i].first.x - current_vtxX) * (180./TMath::Pi());
1290         // cout<<"angle offset test : "<<phi_test<<", "<<phi_test_offset<<endl;
1291
1292         angle_correlation -> Fill(Clus_InnerPhi_Offset, Clus_OuterPhi_Offset);
1293         angle_diff_DCA_dist -> Fill(Clus_InnerPhi_Offset - Clus_OuterPhi_Offset, DCA_sign);
1294         angle_diff -> Fill(abs(Clus_InnerPhi_Offset - Clus_OuterPhi_Offset));
1295         angle_diff_inner_phi -> Fill(Clus_InnerPhi_Offset, Clus_InnerPhi_Offset - Clus_OuterPhi_Offset);
1296         angle_diff_outer_phi -> Fill(Clus_OuterPhi_Offset, Clus_InnerPhi_Offset - Clus_OuterPhi_Offset);
1297         DCA_point -> Fill(DCA_info_vec[1], DCA_info_vec[2]);
1298         DCA_distance_inner_phi -> Fill(Clus_InnerPhi_Offset, DCA_sign);
1299         DCA_distance_outer_phi -> Fill(Clus_OuterPhi_Offset, DCA_sign);
1300         DCA_distance -> Fill(DCA_sign);
1301         DCA_distance_inner_X -> Fill(cluster_pair_vec[i].first.x, DCA_sign);
1302         DCA_distance_inner_Y -> Fill(cluster_pair_vec[i].first.y, DCA_sign);
1303         DCA_distance_outer_X -> Fill(cluster_pair_vec[i].second.x, DCA_sign);
1304         DCA_distance_outer_Y -> Fill(cluster_pair_vec[i].second.y, DCA_sign);
1305
1306         angle_diff_new -> Fill(abs(Clus_InnerPhi_Offset - Clus_OuterPhi_Offset));
1307     } // note : end of the loop for the cluster pair
1308
1309     // note : ----------- ----------- ----------- ---------- ----------- ----------- ---------- ----------- ----------- ----------- -----------
1310     DCA_distance_inner_phi_peak = (TH2F*)DCA_distance_inner_phi -> Clone();
1311     TH2F_threshold(DCA_distance_inner_phi_peak, 0.5); // todo : the background cut can be modified, the ratio 0.5
1312     DCA_distance_inner_phi_peak_profile =  DCA_distance_inner_phi_peak->ProfileX("DCA_distance_inner_phi_peak_profile");
1313     // TGraph * DCA_distance_inner_phi_peak_profile_graph = new TGraph();
1314     double point_index = 0;
1315     vector<double> hist_column_content = SumTH2FColumnContent(DCA_distance_inner_phi_peak);
1316     for (int i = 0; i < DCA_distance_inner_phi_peak_profile->GetNbinsX(); i++){
1317         // if (hist_column_content[i] < 5){continue;} // note : in order to remove some remaining background
1318
1319         DCA_distance_inner_phi_peak_profile_graph -> SetPoint(point_index, DCA_distance_inner_phi_peak_profile->GetBinCenter(i+1), DCA_distance_inner_phi_peak_profile->GetBinContent(i+1));
1320         // cout<<"("<<DCA_distance_inner_phi_peak_profile->GetBinCenter(i+1)<<", "<< DCA_distance_inner_phi_peak_profile->GetBinContent(i+1)<<")"<<endl;
1321         point_index += 1;
1322     }
1323
1324     // cos_fit -> SetParameters(4, 1.49143e-02,  185, 0.3); // todo : we may have to apply more constaints on the fitting
1325     // cos_fit -> SetParLimits(0,0,1000); // note : the amplitude has to be positive
1326     // cos_fit -> SetParLimits(2,0,360); // note : the peak location has to be positive
1327
1328     // note : here is the test with a gaus fitting to find the peak
1329     // gaus_fit -> SetParameters(-4.5, cos_fit->GetParameter(2), 50, 0);
1330     // gaus_fit -> SetParLimits(0,-100,0); // note : the gaus distribution points down
1331     // DCA_distance_inner_phi_peak_profile -> Fit(gaus_fit, "N","",100, 260);
1332     // cout<<"test, gaus fit range : "<<gaus_fit->GetParameter(1) - 25<<" "<<gaus_fit->GetParameter(1) + 25<<endl;
1333
1334     horizontal_fit_inner -> SetParameter(0,0);
1335
1336     // todo : the fit range of the gaussian fit can be modified here
1337     DCA_distance_inner_phi_peak_profile_graph -> Fit(horizontal_fit_inner,"NQ","",0,360);
1338     DCA_distance_inner_phi_peak_profile_graph -> Fit(gaus_fit, "NQ","",cos_fit->GetParameter(2) - guas_fit_range, cos_fit->GetParameter(2) + guas_fit_range); // note : what we want and need is the peak position, so we fit the peak again
1339     DCA_distance_inner_phi_peak_profile_graph -> Fit(cos_fit, "NQ","",cos_fit_rangel, cos_fit_ranger);
1340
1341     // note : ----------- ----------- ----------- ---------- ----------- ----------- ---------- ----------- ----------- ----------- -----------
1342     angle_diff_new_bkg_remove = (TH1F*)angle_diff_new -> Clone();
1343     angle_diff_new_bkg_remove -> SetLineColor(2);
1344     Hist_1D_bkg_remove(angle_diff_new_bkg_remove, 1.5);
1345
1346     // note : ----------- ----------- ----------- ---------- ----------- ----------- ---------- ----------- ----------- ----------- -----------
1347     angle_diff_inner_phi_peak = (TH2F*)angle_diff_inner_phi -> Clone();
1348     // TH2F_threshold_advanced(angle_diff_inner_phi_peak, 0.5);
1349     TH2F_threshold_advanced_2(angle_diff_inner_phi_peak, 0.5); // todo : threshold ratio can be modified here
1350     hist_column_content = SumTH2FColumnContent(angle_diff_inner_phi_peak);
1351     angle_diff_inner_phi_peak_profile =  angle_diff_inner_phi_peak->ProfileX("angle_diff_inner_phi_peak_profile");
1352     // TGraph * angle_diff_inner_phi_peak_profile_graph = new TGraph();
1353     point_index = 0;
1354     for (int i = 0; i < angle_diff_inner_phi_peak_profile->GetNbinsX(); i++){
1355         // if (hist_column_content[i] < 5){continue;} // note : in order to remove some remaining background
1356
1357         angle_diff_inner_phi_peak_profile_graph -> SetPoint(point_index, angle_diff_inner_phi_peak_profile->GetBinCenter(i+1), angle_diff_inner_phi_peak_profile->GetBinContent(i+1));
1358         // cout<<"("<<angle_diff_inner_phi_peak_profile->GetBinCenter(i+1)<<", "<< angle_diff_inner_phi_peak_profile->GetBinContent(i+1)<<")"<<endl;
1359         point_index += 1;
1360     }
1361
1362     horizontal_fit_angle_diff_inner -> SetParameter(0,0);
1363     angle_diff_inner_phi_peak_profile_graph -> Fit(horizontal_fit_angle_diff_inner,"NQ","",0,360);
1364
1365     // note : ----------- ----------- ----------- ---------- ----------- ----------- ---------- ----------- ----------- ----------- -----------
1366     DCA_distance_outer_phi_peak = (TH2F*)DCA_distance_outer_phi -> Clone();
1367     TH2F_threshold(DCA_distance_outer_phi_peak, 0.5); // todo : the background cut can be modified, the ratio 0.5
1368     DCA_distance_outer_phi_peak_profile =  DCA_distance_outer_phi_peak->ProfileX("DCA_distance_outer_phi_peak_profile");
1369     // TGraph * DCA_distance_outer_phi_peak_profile_graph = new TGraph();
1370     point_index = 0;
1371     hist_column_content = SumTH2FColumnContent(DCA_distance_outer_phi_peak);
1372     for (int i = 0; i < DCA_distance_outer_phi_peak_profile->GetNbinsX(); i++){
1373         // if (hist_column_content[i] < 5){continue;} // note : in order to remove some remaining background
1374
1375         DCA_distance_outer_phi_peak_profile_graph -> SetPoint(point_index, DCA_distance_outer_phi_peak_profile->GetBinCenter(i+1), DCA_distance_outer_phi_peak_profile->GetBinContent(i+1));
1376         // cout<<"("<<DCA_distance_outer_phi_peak_profile->GetBinCenter(i+1)<<", "<< DCA_distance_outer_phi_peak_profile->GetBinContent(i+1)<<")"<<endl;
1377         point_index += 1;
1378     }
1379
1380     horizontal_fit_outer -> SetParameter(0,0);
1381     // todo : the fit range of the gaussian fit can be modified here
1382     DCA_distance_outer_phi_peak_profile_graph -> Fit(horizontal_fit_outer,"NQ","",0,360);
1383
1384     // note : ----------- ----------- ----------- ---------- ----------- ----------- ---------- ----------- ----------- ----------- -----------
1385     angle_diff_outer_phi_peak = (TH2F*)angle_diff_outer_phi -> Clone();
1386     // TH2F_threshold_advanced(angle_diff_outer_phi_peak, 0.5);
1387     TH2F_threshold_advanced_2(angle_diff_outer_phi_peak, 0.5); // todo : threshold ratio can be modified here
1388     hist_column_content = SumTH2FColumnContent(angle_diff_outer_phi_peak);
1389     angle_diff_outer_phi_peak_profile =  angle_diff_outer_phi_peak->ProfileX("angle_diff_outer_phi_peak_profile");
1390     // TGraph * angle_diff_outer_phi_peak_profile_graph = new TGraph();
1391     point_index = 0;
1392     for (int i = 0; i < angle_diff_outer_phi_peak_profile->GetNbinsX(); i++){
1393         // if (hist_column_content[i] < 5){continue;} // note : in order to remove some remaining background
1394
1395         angle_diff_outer_phi_peak_profile_graph -> SetPoint(point_index, angle_diff_outer_phi_peak_profile->GetBinCenter(i+1), angle_diff_outer_phi_peak_profile->GetBinContent(i+1));
1396         // cout<<"("<<angle_diff_outer_phi_peak_profile->GetBinCenter(i+1)<<", "<< angle_diff_outer_phi_peak_profile->GetBinContent(i+1)<<")"<<endl;
1397         point_index += 1;
1398     }
1399
1400     horizontal_fit_angle_diff_outer -> SetParameter(0,0);
1401     angle_diff_outer_phi_peak_profile_graph -> Fit(horizontal_fit_angle_diff_outer,"NQ","",0,360);
1402     // note : ----------- ----------- ----------- ---------- ----------- ----------- ---------- ----------- ----------- ----------- -----------
1403
1404
1405
1406     angle_diff_inner_phi_peak_profile_graph -> Set(0);
1407     angle_diff_outer_phi_peak_profile_graph -> Set(0);
1408     DCA_distance_inner_phi_peak_profile_graph -> Set(0);
1409     DCA_distance_outer_phi_peak_profile_graph -> Set(0);
1410
1411     if (print_message_opt == true) {cout<<"circle radius : "<<abs(gaus_fit->GetParameter(0) + gaus_fit->GetParameter(3))<<" possible correction angle : "<<gaus_fit->GetParameter(1)<<endl;}
1412     // return {abs(gaus_fit->GetParameter(0) + gaus_fit->GetParameter(3)), gaus_fit->GetParameter(1)}; // note : gaussian based
1413
1414     // return {
1415     //     abs(cos_fit->GetParameter(0) - cos_fit->GetParameter(3)), cos_fit->GetParameter(2), // note : cosine wave based
1416     //     angle_diff_new_bkg_remove -> GetStdDev(), angle_diff_new_bkg_remove -> GetStdDevError() // note : the std dev of the angle diff 1D distribution
1417     // };
1418 }
1419
1420 pair<double,double> INTTXYvtxEvt::Get_possible_zvtx(double rvtx, vector<double> p0, vector<double> p1) // note : inner p0, outer p1, vector {r,z}, -> {y,x}
1421 {
1422     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}
1423     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}
1424
1425     double edge_first  = Get_extrapolation(rvtx,p0_z_edge[0],p0[0],p1_z_edge[1],p1[0]);
1426     double edge_second = Get_extrapolation(rvtx,p0_z_edge[1],p0[0],p1_z_edge[0],p1[0]);
1427
1428     double mid_point = (edge_first + edge_second) / 2.;
1429     double possible_width = fabs(edge_first - edge_second) / 2.;
1430
1431     return {mid_point, possible_width}; // note : first : mid point, second : width
1432
1433 }
1434
1435 double INTTXYvtxEvt::Get_extrapolation(double given_y, double p0x, double p0y, double p1x, double p1y) // note : x : z, y : r
1436 {
1437     if ( fabs(p0x - p1x) < 0.00001 ){ // note : the line is vertical (if z is along the x axis)
1438         return p0x;
1439     }
1440     else {
1441         double slope = (p1y - p0y) / (p1x - p0x);
1442         double yIntercept = p0y - slope * p0x;
1443         double xCoordinate = (given_y - yIntercept) / slope;
1444         return xCoordinate;
1445     }
1446 }
1447
1448 void INTTXYvtxEvt::EndRun()
1449 {
1450     inner_pos_xy -> Reset("ICESM");
1451     outer_pos_xy -> Reset("ICESM");
1452     inner_outer_pos_xy -> Reset("ICESM");
1453     N_cluster_correlation -> Reset("ICESM");
1454     N_cluster_correlation_close -> Reset("ICESM");
1455     if (draw_event_display == true) {c2 -> Print(Form("%s/temp_event_display.pdf)",out_folder_directory.c_str()));};
1456     if (draw_event_display == true) {c1 -> Print(Form("%s/evt_LineFill2D.pdf)",out_folder_directory.c_str()));}
1457
1458     file_out -> cd();
1459     tree_out -> SetDirectory(file_out);
1460     tree_out -> Write();
1461     file_out -> Close();
1462
1463     return;
1464 }
1465
1466
1467 // note : this function shouldn't be used for the official z vertex calculation, it should only be used for the event display
1468 double INTTXYvtxEvt::get_z_vertex(clu_info inner_clu, clu_info outer_clu, double target_x, double target_y)
1469 {
1470     // note : x = z,
1471     // note : y = radius
1472     double inner_clu_r = sqrt(pow(inner_clu.x,2)+ pow(inner_clu.y,2));
1473     double outer_clu_r = sqrt(pow(outer_clu.x,2)+ pow(outer_clu.y,2));
1474     double target_r    = sqrt(pow(target_x,2)   + pow(target_y,2));
1475
1476     // Use the slope equation (y = ax + b) to calculate the x-coordinate for the target y
1477     if ( fabs(outer_clu.z - inner_clu.z) < 0.00001 ){
1478         return outer_clu.z;
1479     }
1480     else {
1481         double slope = (outer_clu_r - inner_clu_r) / (outer_clu.z - inner_clu.z);
1482         double yIntercept = inner_clu_r - slope * inner_clu.z;
1483         double xCoordinate = (target_r - yIntercept) / slope;
1484         return xCoordinate;
1485     }
1486
1487 }
1488
1489 void INTTXYvtxEvt::temp_bkg(TPad * c1, double peek, pair<double,double> beam_origin, InttConversion * ch_pos_DB)
1490 {
1491     c1 -> cd();
1492
1493     int N_ladder[4] = {12, 12, 16, 16};
1494     string ladder_index_string[16] = {"00","01","02","03","04","05","06","07","08","09","10","11","12","13","14","15"};
1495
1496     vector<double> x_vec; x_vec.clear();
1497     vector<double> y_vec; y_vec.clear();
1498
1499     vector<double> x_vec_2; x_vec_2.clear();
1500     vector<double> y_vec_2; y_vec_2.clear();
1501
1502     bkg -> SetTitle("INTT event display X-Y plane");
1503     bkg -> SetMarkerStyle(20);
1504     bkg -> SetMarkerSize(0.1);
1505     // bkg -> SetPoint(1,beam_origin.first,beam_origin.second);
1506     bkg -> GetXaxis() -> SetLimits(-150,150);
1507     bkg -> GetYaxis() -> SetRangeUser(-150,150);
1508     bkg -> GetXaxis() -> SetTitle("X [mm]");
1509     bkg -> GetYaxis() -> SetTitle("Y [mm]");
1510
1511     bkg -> Draw("ap");
1512
1513
1514
1515     for (int server_i = 0; server_i < 4; server_i++)
1516     {
1517         for (int FC_i = 0; FC_i < 14; FC_i++)
1518         {
1519             ladder_line -> DrawLine(
1520                 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,
1521                 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
1522             );
1523         }
1524     }
1525
1526     ladder_line -> Draw("l same");
1527 }
1528
1529 bool INTTXYvtxEvt::isPointInsideSquare(const std::pair<double, double> point, const std::pair<double, double> square_center, double length) {
1530     // Calculate half-length of the square
1531     double halfLength = length / 2.0;
1532
1533     // Calculate boundaries of the square
1534     double x_min = square_center.first - halfLength;
1535     double x_max = square_center.first + halfLength;
1536     double y_min = square_center.second - halfLength;
1537     double y_max = square_center.second + halfLength;
1538
1539     // Check if the point lies inside the square
1540     return (point.first > x_min && point.first < x_max && point.second > y_min && point.second < y_max);
1541 }
1542
1543 void INTTXYvtxEvt::TH2LineFill(TH2F * hist_in, std::pair<double,double> point_1, std::pair<double,double> point_2)
1544 {
1545     double cell_width_x = hist_in -> GetXaxis() -> GetBinWidth(1);
1546     double cell_width_y = hist_in -> GetYaxis() -> GetBinWidth(1);
1547     if (cell_width_x != cell_width_y) {cout<<"the size of the cell is not square!"<<endl; exit(1);}
1548
1549     for (int xi = 1; xi < hist_in -> GetNbinsX()+1; xi++)
1550     {
1551         for (int yi = 1; yi < hist_in -> GetNbinsY()+1; yi++)
1552         {
1553             vector<double> DCA_info = calculateDistanceAndClosestPoint(
1554                 point_1.first, point_1.second,
1555                 point_2.first, point_2.second,
1556                 hist_in -> GetXaxis() -> GetBinCenter(xi), hist_in -> GetYaxis() -> GetBinCenter(yi)
1557             );
1558
1559             // cout<<" "<<endl;
1560             // cout<<xi<<" "<<yi<<" bin center : "<<hist_in -> GetXaxis() -> GetBinCenter(xi)<<" "<< hist_in -> GetYaxis() -> GetBinCenter(yi)<<" DCA: "<<DCA_info[0]<<" DCA point : "<<DCA_info[1]<<" "<<DCA_info[2]<<endl;
1561
1562             if (isPointInsideSquare({DCA_info[1],DCA_info[2]},{hist_in -> GetXaxis() -> GetBinCenter(xi), hist_in -> GetYaxis() -> GetBinCenter(yi)},cell_width_x))
1563             {
1564                 // cout<<"------- before check bin content: "<<xi<<" "<<yi<<" "<<hist_in -> GetBinContent(xi,xi)<<endl;
1565                 // cout<<"------- So the cell is filled "<<xi<<" "<<yi<<endl;
1566                 hist_in -> SetBinContent(xi,yi,hist_in -> GetBinContent(xi,yi) + 1);
1567                 // cout<<"------- after check bin content: "<<xi<<" "<<yi<<" "<<hist_in -> GetBinContent(xi,xi)<<endl;
1568             }
1569         }
1570     }
1571 }
1572
1573
1574
1575
1576 #endif